momente şi schiţe de informatică şi matematică
To attain knowledge, write. To attain wisdom, rewrite.

Studiul datelor evaluării naţionale din 2016 (folosind R)

Evaluare naţională | limbajul R
2016 oct

Am lansat LibreOffice pe adresa furnizată de data.gov.ro pentru descărcarea setului de date conţinând "rezultatele anonimizate ale Evaluării Naționale 2016":

vb@Home:~/16_evn$ url_res="http://data.gov.ro/dataset/\
> 5eac43d1-b364-468f-bdff-18e668abfb20/resource/8d31a801-a430-47e3-96de-04bc4cc3ce55/\
> download/evnat20162.ods"
vb@Home:~/16_evn$ libreoffice $url_res

Interpretorul de comenzi bash permite separarea pe mai multe rânduri a unui text (încadrat de ghilimele), folosind '\' (caracterul "backslash") la capătul fiecărui rând; valoarea variabilei se referă apoi folosind '$'.

Am descărcat fişierul în format "ODS" (care - salvat pe disc şi apoi dezarhivat - are aproape 400 MB) nu pentru a-l folosi în vreun fel, ci doar pentru a obţine de sub LibreOffice fişierul în format "CSV" standard (cu acesta se va putea lucra).

Este drept că puteam descărca direct "evnat20162.csv" (în loc să folosesc intermediar, fişierul ".ods") - dar acesta are un format nestandard şi necesita transformări (v. [1] §1); de exemplu, foloseşte ',' şi nu '.' ca separator de zecimale (ceea ce poate că este convenabil pentru afişarea datelor, dar nu convine pentru a opera cu numere în programe); în plus, foloseşte codul specific Windows (doi octeţi, 0x0D0A) pentru separarea liniilor - ori de sub LibreOffice am obţinut direct codul specific Linux (un singur octet, 0x0A).

Inspectăm întâi fişierul CSV obţinut astfel (folosind vizualizatorul încorporat în GNOME Commander):

Trebuind să ne referim în programe la câmpurile respective, vom înlocui scrierea ALL CAPS şi vom simplifica denumirile (evitând în orice caz să folosim caracterul "spaţiu", pentru a nu fi obligaţi ulterior să încadrăm între ghilimele); vom renunţa la unele câmpuri (de exemplu, nu are de ce să ne intereseze codurile "SIIIR" - încât le vom înlocui cu denumirile judeţelor).

Lansăm R şi constituim un obiect de clasă 'data.frame' pentru setul de date din fişierul nostru:

vb@Home:~/16_evn$ R -q  # --quiet  (omite mesajul introductiv)
> evn <- read.csv("evnat20162.csv")
> str(evn)  # inspectăm structura de date creată
'data.frame':	153673 obs. of  21 variables:  # am adăugat aici şi indecşii câmpurilor
[1] $ COD.UNIC.CANDIDAT          : int  158870 136940 106543 152070 121592 106178 ...
[2] $ SEX                        : Factor w/ 2 levels "F","M": 2 2 1 2 1 1 1 1 2 ...
[3] $ MEDIU                      : Factor w/ 2 levels "RURAL","URBAN": 1 2 1 1 1 ...
[4] $ COD.SIIIR                  : num  2.56e+09 3.26e+09 2.56e+09 2.26e+09 ...
[5] $ STATUS.ROMANA              : Factor w/ 3 levels "ABSENT","ELIMINAT",..: 3 3 ...
[6] $ STATUS.LIMBA.MATERNA       : Factor w/ 3 levels "-","ABSENT","PREZENT": 1 1 ...
[7] $ STATUS.MATEMATICA          : Factor w/ 3 levels "ABSENT","ELIMINAT",..: 3 3 ...
[8] $ NOTA.ROMANA                : num  3.2 5.1 6.05 6.4 7.4 9.5 4.1 9.9 8.9 8.65 ...
[9] $ NOTA.LIMBA.MATERNA         : num  NA NA NA NA NA NA NA NA NA NA ...
[10] $ NOTA.MATEMATICA            : num  2.2 6.5 4.15 6.8 3.85 7.95 4.1 8 9.85 7.6 ...
[11] $ CONTESTATIE.ROMANA         : Factor w/ 2 levels "DA","NU": 2 2 2 2 2 2 2 2 1 ...
[12] $ NOTA.CONTESTATIE.ROMANA    : num  NA NA NA NA NA NA NA NA 8.85 NA ...
[13] $ CONTESTATIE.LIMBA.MATERNA  : Factor w/ 2 levels "DA","NU": 2 2 2 2 2 2 2 2 2 ...
[14] $ NOTA.CONTESTATIE.LB.MATERNA: num  NA NA NA NA NA NA NA NA NA NA ...
[15] $ CONTESTATIE.MATEMATICA     : Factor w/ 2 levels "DA","NU": 2 2 2 2 2 2 2 2 2 ...
[16] $ NOTA.CONTESTATIE.MATEMATICA: num  NA NA NA NA NA NA NA NA NA NA ...
[17] $ NOTA.FINALA.ROMANA         : num  3.2 5.1 6.05 6.4 7.4 9.5 4.1 9.9 8.9 8.65 ...
[18] $ NOTA.FINALA.LB.MATERNA     : num  NA NA NA NA NA NA NA NA NA NA ...
[19] $ NOTA.FINALA.MATEMATICA     : num  2.2 6.5 4.15 6.8 3.85 7.95 4.1 8 9.85 7.6 ...
[20] $ MEDIA                      : num  2.7 5.8 5.1 6.6 5.62 8.72 4.1 8.95 9.37 ...
[21] $ MEDIA.V.VIII               : num  6.94 8.69 9.21 8 8.84 9.58 7.9 9.76 8.66 ...

Pecizăm că mai toate comenzile pe care le folosim aici au fost deja descrise în [1] şi [2]; în redarea secvenţelor de comenzi vom omite de obicei, promptul consolei R (caracterele iniţiale "> "); de altfel, secvenţele mai lungi le scriem într-un fişier "prog.R" separat (deci, fără prompt) şi folosim mecanismul asigurat de comanda source("prog.R").

Faţă de 2015 - v. [3] - numărul candidaţilor s-a micşorat cu aproape 10000. Ca şi în 2015, avem 21 de câmpuri - dar acum lipseşte câmpul "DATA NASTERE" şi în schimb, avem câmpul (cu siguranţă mai interesant) "MEDIA V-VIII".

Ne întrebăm dacă n-am putea renunţa la câmpurile "STATUS" (de indecşi 5-7):

summary(evn[5:7])
  STATUS.ROMANA    STATUS.LIMBA.MATERNA STATUS.MATEMATICA
 ABSENT  :  4543   -      :144528       ABSENT  :  5020  
 ELIMINAT:     2   ABSENT :   238       ELIMINAT:     5  
 PREZENT :149128   PREZENT:  8907       PREZENT :148648  

Au fost eliminaţi numai 7 elevi, iar proporţia totală a absenţilor este 9563 / 153673 ≈ 6.22%; nu zicem că ar fi mică, dar situaţia de absent (sau eliminat) este reflectată pe câmpul "MEDIA": în aceste cazuri, media lipseşte - fiind înlocuită de către R cu valoarea specifică 'NA'. Prin urmare, putem colecta indicii acestor câmpuri, pentru eliminare:

col_to_remove <- c(5:7)  # colectează indicii coloanelor pe care le vom elimina în final

De asemenea, putem elimina cele trei câmpuri "CONTESTATIE" (cu valori "DA/NU"): pentru "DA" avem înscrisă nota corespunzătoare contestaţiei în câmpul asociat "NOTA CONTESTATIE" (iar altfel, acest câmp conţine valoarea 'NA'):

append(col_to_remove, c(11, 13, 15)) -> col_to_remove

Câmpul "COD SIIIR" conţine pentru fiecare candidat, codul SIIIR corespunzător acestuia; extragem codul judeţului (prima sau primele două cifre, după cum SIIIR are 9 sau 10 cifre) - folosind funcţia jud_from_siiir() din [1] - şi înfiinţăm coloana 'jud' (de clasă factor) conţinând denumirile judeţelor - folosind pentru aceasta datele din [1]:

evn$jud <- jud_from_siiir(evn[[4]])  # codul "NN" al judeţului candidatului (din COD.SIIIR)
evn$jud <- as.factor(evn$jud)
load("../15_bac/bac5.RData")  # recuperează 'bac5', pentru nivelurile factorului 'jud'
levels(evn$jud) <- levels(bac5$jud)
rm(bac5)  # elimină 'bac5' din sesiunea de lucru curentă (rm() = "remove")

Noul câmp 'jud' a fost adăugat la sfârşit, încât nu afectează indicii redaţi la început ai câmpurilor existente deja. Ştergem acum, coloanele colectate în acest scop (incluzând şi 'COD.SIIIR'):

evn[, c(col_to_remove, 4)] <- NULL
names(evn)  # numele coloanelor rămase
 [1] "COD.UNIC.CANDIDAT"           "SEX"                         "MEDIU"                      
 [4] "NOTA.ROMANA"                 "NOTA.LIMBA.MATERNA"          "NOTA.MATEMATICA"            
 [7] "NOTA.CONTESTATIE.ROMANA"     "NOTA.CONTESTATIE.LB.MATERNA" "NOTA.CONTESTATIE.MATEMATICA"
[10] "NOTA.FINALA.ROMANA"          "NOTA.FINALA.LB.MATERNA"      "NOTA.FINALA.MATEMATICA"     
[13] "MEDIA"                       "MEDIA.V.VIII"                "jud"                        

Ne bazăm de-acum pe noii indici de coloane, dar modificăm denumirile astfel:

names(evn) <- c("id", "Sex", "Mediu", 
                "Română1", "Maternă1", "Matematică1",  # notele iniţiale
                "Română2", "Maternă2", "Matematică2",  # note contestaţii
                "Română", "Maternă", "Matematică",     # note finale
                "Media", "Media8", "Jud")

Mai departe - rămâne să vedem ce spun datele colectate în setul de date "evn"; dacă ar fi să ne luăm numai după ce ştim din [1]-[3] - fără alte investigaţii - am zice de exemplu, că rezultatele examenului sunt mai bune pentru fete decât pentru băieţi şi la fel, pentru mediul "urban" faţă de cel "rural" (şi la fel pentru judeţul Cluj faţă de judeţul Ilfov, ş.a.m.d.).

1. Repartiţia candidaţilor după sex, mediu, sau judeţ

Faţă de 2015 - v. [3] - proporţia "F" = 51%, "M" = 49% se menţine, dar a scăzut cu câte 1% la "RURAL" şi a crescut cu câte 1% la "URBAN":

mediu_sex <- addmargins(table(evn[c('Mediu', 'Sex')]))
nrcd <- mediu_sex[3,3]  # total: 153673 candidaţi
mediu_sex <- 100*round(mediu_sex / nrcd, 2)

mediu_sex # 2016
       Sex  # valori procentuale (aproximative)
Mediu     F   M Sum
  RURAL  23  22  45  # cu câte 1% mai puţin ca în 2015
  URBAN  28  27  55  # cu câte 1% mai mult ca în 2015
  Sum    51  49 100

Graficul redat mai sus l-am obţinut (v. şi [2] §5) prin:

require(lattice)
barchart(mediu_sex, aspect=1, stack=FALSE, xlab="Procente",
    auto.key = list(columns=3, cex=0.8, size=1.5, points=FALSE, rectangles=TRUE),
    panel = function(...) { panel.barchart(...)
                            panel.abline(v=seq(0,100,10), lwd=0.6, col="lightgrey") }
)

Să vedem şi repartiţia pe judeţe, a candidaţilor:

nrcd <- nrow(evn)
rep_jud <- as.data.frame(sort(table(evn$Jud) / nrcd))
require(ggplot2)
require(scales)  # scales::percent() (formatare procentuală)
ggplot(rep_jud, aes(x = Var1, y = Freq)) + 
    geom_bar(stat='identity', col="white", fill="lightgreen", width=0.85) + 
    labs(x="", y="", title = paste0("Repartiţia pe judeţe (", 
                                    nrcd, " candidaţi, 2016)")) + 
    scale_y_continuous(labels = percent) + coord_flip()

Candidaţii din "M.Bucureşti" sunt cu aproape 4% până la 7% mai mulţi, decât în celelalte judeţe. Următoarea secvenţă produce diferenţele faţă de 2015 (dar de data aceasta nu ne convine calculul procentual):

load("evna.RData")  # recuperăm datele din 2015 (163418 candidaţi) - v. [3]
rep_jud5 <- as.data.frame(table(evna$Jud))
rep_jud6 <- as.data.frame(table(evn$Jud))
rep_jud6$difer <- rep_jud6$Freq - rep_jud5$Freq
rm(list = c("evna", "rep_jud5"))  # elimină obiectele de care nu mai avem nevoie
print(data.frame(rep_jud6[1:14, ], rep_jud6[15:28, ], rep_jud6[29:42, ]))
            Var1 Freq difer    Var1.1 Freq.1 difer.1      Var1.2 Freq.2 difer.2
            Alba 2588  -256 Dâmboviţa   4290    -234     Prahova   5279    -284
            Arad 3172  -191      Dolj   4571    -765   Satu-Mare   2565    -194
           Argeş 4930  -149    Galaţi   3917    -595       Sălaj   1997    -136
           Bacău 5226  -326      Gorj   3169    -214       Sibiu   2808    -142
           Bihor 4704  -234  Harghita   2774       7     Suceava   6678    -466
 Bistriţa-Năsăud 2557  -244 Hunedoara   3097    -155   Teleorman   2530    -327
        Botoşani 3986  -121  Ialomiţa   2049    -193       Timiş   4477    -170
          Braşov 3622   -88      Iaşi   6332    -414      Tulcea   1594     -26
          Brăila 2023   -53     Ilfov   2257     -59      Vaslui   3593    -407
           Buzău 3380  -302 Maramureş   4053    -231      Vâlcea   2970    -209
   Caraş-Severin 2040  -242 Mehedinţi   2030     -80     Vrancea   2817    -144
            Cluj 4102  -291     Mureş   4298    -190 M.Bucureşti  11983     -69
       Constanţa 5555  -161     Neamţ   4372    -486    Călăraşi   2085    -288
         Covasna 1532   -54       Olt   3280    -524     Giurgiu   2391     -38

Exceptând numai Harghita (unde a crescut cu 7), în toate judeţele numărul candidaţilor a scăzut faţă de 2015 (cel mai mult la Dolj (-765), Galaţi, Neamţ, Olt, Iaşi, Vaslui (-407)) - scăderea totală (obţinută prin sum(rep_jud6$difer)) fiind -9745.

Dacă am fi lucrat procentual - raportând numărul candidaţilor din fiecare judeţ la numărul total de candidaţi din 2016 şi respectiv, la cel din 2015 - atunci diferenţele procentelor respective ar avea suma 0 (reducându-se la diferenţa dintre suma procentelor judeţene din 2016 şi aceea din 2015, deci la 100% - 100%) şi ar fi fost greu de imaginat cum să le interpretăm corect.

2. Distribuţia mediilor de examen

Să verificăm întâi, notele şi mediile (amintindu-ne surpriza din [1] §9.5, când în coloana de note am găsit şi valori ca "-1" şi "-2" - însemnând calitatea de "absent", sau "eliminat"); aplicăm funcţia summary() fiecăreia dintre coloanele respective şi apoi transpunem matricea rezultată:

sapply(evn[, 4:14], FUN = summary) -> brev
t(brev) -> brev
cbind(brev, brev[, 3] - brev[, 4])  # adaugă coloana diferenţelor Median-Mean
            Min. 1st Qu. Median  Mean 3rd Qu. Max.   NA's  Median - Mean     
Română1     1.00    5.40   7.25 6.870    8.60   10   4545      0.380  # note iniţiale
Maternă1    1.00    6.35   7.80 7.415    8.80   10 144766      0.385
Matematică1 1.00    4.30   6.10 6.056    7.95   10   5025      0.044
Română2     1.20    7.05   8.15 7.803    8.85   10 141739      0.347  # contestaţii
Maternă2    3.45    7.10   8.10 7.821    8.70   10 153281      0.279
Matematică2 1.00    5.90   7.30 6.971    8.20   10 144774      0.329
Română      1.00    5.40   7.25 6.884    8.65   10   4545      0.366  # note finale
Maternă     1.00    6.40   7.80 7.432    8.80   10 144766      0.368
Matematică  1.00    4.30   6.10 6.058    7.95   10   5025      0.042
Media       1.00    5.00   6.70 6.509    8.22   10   5025      0.191  # medii finale
Media8      5.13    7.89   8.79 8.628    9.48   10      5      0.162  # medii V-VIII

De data aceasta, toate notele sunt ceea ce trebuie să fie - valori din gama [1, 10]. Observăm că pentru toate cele 11 coloane, valoarea medie 'Mean' este mai mică decât valoarea mediană Median; aceasta înseamnă de regulă, că valorile respective au o distribuţie asimetrică, având o plajă spre stânga medianei ceva mai întinsă faţă de coada (mai abruptă) spre mediile mai mari ca mediana (vom vedea acest aspect pe graficele distribuţiilor, mai jos).

Am adăugat (prin cbind()) şi o coloană pe care am calculat diferenţele dintre valorile mediane şi cele medii - mult mai mari la română decât la matematică (0.38 faţă de 0.04), încât probabil că asimetria intuită mai sus este mai pronunţată la română, faţă de matematică.

Cei care nu au notă la una sau alta dintre probe (fiind absenţi sau eliminaţi) au fost contorizaţi în coloana 'NA's şi constatăm că numărul acestora la oricare probă este sub 10000. Pe de altă parte, observăm că valorile medianei şi mediei pentru notele iniţiale aproape coincid, cu cele pentru notele finale (în cazul fiecăreia dintre cele trei probe). Prin urmare - proporţia de contestaţii pe fiecare probă fiind relativ mică şi în plus, neavând în urma rezolvării contestaţiilor, modificări sensibile ale valorilor centrale - putem să ignorăm de-acum "contestaţiile" (considerându-le ca nerelevante din punct de vedere statistic).

Valorile medianei şi mediei pentru coloana 'Media' sunt cu 2 puncte mai mici decât cele pentru 'Media8' - dar nu prea putem interpreta acum, această diferenţă: valorile respective încep într-un caz de la 1.00, iar în celălalt de la 5.13 (iar pe de altă parte, valorile din 'Medie8' acoperă toate obiectele şcolare, nu numai pe cele două sau trei din care au rezultat valorile din 'Medie').

Există obiceiul de a categorisi notele pe intervale: "sub 5", [5, 6), [6, 7), ..., [9, 10]; dar dacă vrem să evidenţiem distribuţia lor, cel mai convenabil este să le considerăm ca atare, lăsând gruparea acestora în seama funcţiilor care produc histograme sau curbe de densitate.

Programul următor (scris într-un fişier "Media.R" şi lansat din consola R prin comanda source("Media.R")) trasează o histogramă a mediilor finale (grupându-le pe intervale de câte o zecime), curba densităţii estimate pentru acestea şi curba distribuţiei normale având ca parametri valoarea medie şi dispersia acestora:

opar <- par(mar=c(2,2,0,0)+0.1, cex.axis=0.9, xpd=FALSE)
    # histograma:
hist(evn$Media, breaks = "FD", freq = FALSE, right = FALSE, 
     col="cornsilk", border="lightgray", ylim=c(0, 0.21), main="")
grid(lwd = 1.1)
    # estimarea densităţii:
dst <- density(evn$Media, from = 1, to = 10, adjust = 0.75, na.rm=TRUE)
lines(dst, col="blue", lwd=1.5)
    # densitatea distribuţiei normale N(mev, sev):
mev <- mean(evn$Media, na.rm=TRUE)  # media notelor (= 6.509183)
sev <- sd(evn$Media, na.rm=TRUE)  # dispersia (= 2.044166)
lines(dst$x, dnorm(dst$x, mean = mev, sd = sev),  lwd=1.5, col="magenta")
    # zona mev ± sev cuprinde 68.27% dintre valori (în distribuţia normală):
msev <- c(mev-sev, mev, mev+sev)
rmsev <- round(msev, 2)
abline(v = msev, lwd=1, lty="dotted", col="red")
    # legende şi etichetări:
legend("topleft", 
       legend=c("densitatea estimată", 
                expression("distribuţia normală N("*mu*","~sigma*")"), "",
                as.expression(bquote(mu ~"="~ .(rmsev[2]))),
                bquote(sigma ~"="~ .(rmsev[2]-rmsev[1]))),
       text.col=c("blue", "magenta", "white", "black", "black"), 
       bty="n", cex=0.9, inset=c(0, 0.15))
text(c(mev-sev, mev, mev+sev), 0, 
     labels = c(as.expression(bquote(mu ~ "-" ~ sigma ~ "=" ~ .(rmsev[1]))), 
                bquote(mu == .(rmsev[2])), 
                bquote(mu ~ "+" ~sigma == .(rmsev[3]))),
              cex=0.9, col="red", pos=1, offset=0.22)
legend("top", legend="Densitatea distribuţiei mediilor, Evaluare Naţională 2016", 
       bty="n", text.font=3, cex=1.2, xjust=0.5, yjust=0.5)
par(opar)  # reconstituie parametrii grafici standard

În maniera în care este folosită mai sus, funcţia expression() îşi interpretează argumentul ca "expresie matematică" (cu sensul şi sintaxa din TeX; numele 'mu', 'sigma' produc la afişarea expresiei literele greceşti astfel denumite, iar operatorii '*' şi '~' asigură concatenarea şi spaţierea).
Funcţia bquote() permite (prin construcţia .( )) combinarea unei expresii cu valori ale unor variabile. Dar există o problemă la folosirea ei într-un vector de expresii construit cu funcţia c(), cum avem mai sus în funcţia text() (anume, se returnează un obiect de clasă 'list', în loc de 'expression'); la stackoverflow am găsit ca soluţie convertirea prin as.expression() a fiecăreia dintre aplicările bquote() existente în construcţia c() respectivă (folosind eventual mecanismul sapply()) - dar în secvenţa de mai sus am demonstrat deja că este suficient să facem aceasta doar pentru una singură dintre componente.

Graficul obţinut mai sus arată că distribuţia mediilor finale nu poate fi considerată ca rezultând din cea normală. Avem o confirmare a acestui fapt şi printr-un calcul mai simplu: într-o distribuţie N(μ, σ) între limitele μ-σ şi μ+σ avem 68.27% dintre valorile respective, pe când în distribuţia mediilor finale avem numai 62.74%:

medii <- subset(evn, Media >= msev[1] & Media <= msev[3])
nrow(medii) / (nrow(evn) - 5025)  # minus 5025 valori 'NA' (v. tabelul 'brev' mai sus)
[1] 0.6273815  # 62.74% medii în intervalul [4.47, 8.55] ("normal" ar fi 68.27%)

Calcule similare celui redat mai sus ne arată că avem 18.45% medii sub 4.47 şi 18.81% medii peste 8.55.

3. 'F' versus 'M' (distribuţia mediilor condiţionată de sex)

Programul următor determină şi plotează densitatea mediilor pentru fiecare valoare a variabilei 'Sex':

dens <- density(evn$Media[evn$Sex == "F"], from = 1, to = 10, na.rm = TRUE)
dens1 <- density(evn$Media[evn$Sex == "M"], from = 1, to = 10, na.rm = TRUE)
plot(dens, col = "firebrick2", lwd = 1.5, cex.axis = 0.9, xlab = "", ylab = "",
     main = "Densitatea mediilor după sex", cex.main = 1, font.main = 3)
lines(dens1, col="blue", lwd=1.5)
legend("topleft", bty = "n", inset = c(0, 0.05), horiz=TRUE,
       legend=c("F", "M"), title="Sex", title.col="black",
                text.col= c("firebrick2", "blue"))
grid()

În privinţa mediilor finale, avem un decalaj important între 'M' şi 'F'; proporţia mediilor sub 6 este sensibil mai mare, iar a acelora peste 6 este sensibil mai mică, la 'M' faţă de 'F'. Nu avem cum explica acest decalaj; observăm doar că distribuţiile redate mai sus seamănă (iarăşi inexplicabil) cu cele pe care le obţinem mai jos pentru notele de la probele de bază.

4. Română versus Matematică

Tabelul "brev" constituit mai sus indică deosebiri importante, între categoriile de note finale (dar precizăm că ignorăm 'Maternă' - proba suplimentară, susţinută de doar 8907 candidaţi, reprezentând abia 5.8% din total). Quartilele întâia şi a treia (în spatele cărora sunt 25% şi respectiv 75% dintre valori) sunt cu 1.1 şi respectiv 0.7 mai mici pentru 'Matematică', decât pentru 'Română' iar diferenţa 'Median - Mean' este mult mai mare la 'Română' (0.366 faţă de 0.042).

În secvenţa următoare preluăm quartilele celor două probe din tabelul "brev" (pentru a le marca ulterior pe graficul densităţilor), folosim mecanismul lapply() pentru a obţine "într-un singur pas" (spre deosebire de programul redat în §3) cele două densităţi (pentru coloanele 2 şi respectiv 10 din "evn", corespunzătoare notelor finale pentru cele două probe) şi apoi plotăm curbele densităţilor:

q13 <- c(brev[9, 2], brev[9, 5], brev[7, 2], brev[7, 5])  # quartilele 1 (25%) şi 3 (75%)
cmr <- c("red", "red", "blue", "blue")  # pentru linii şi etichete, la cele două probe
dst <- lapply(evn[c(10, 12)], FUN = density,  # estimarea densităţilor, pe cele două probe
              from = 1, to = 10, adjust = 1.75, na.rm = TRUE)
plot(dst[[1]], col = "blue", lwd = 1.5, cex.axis = 0.9, xlab = "", ylab = "",
     main = "Română vs. Matematică (Ev. Nat. 2016)", cex.main = 1, font.main = 3)
lines(dst[[2]], col = "red", lwd = 1.5)
abline(v = q13, lty="dotted", col = cmr)  # marchează quartilele fiecărei probe
mtext(c("q1", "q3", "q1", "q3"), side=1, line=0, cex=0.8, at = q13, col = cmr)
legend("topleft", bty = "n", cex = 0.9, inset = c(0, 0.05),
       legend=c("Română", "Matematică", "", "q1 = 25%, q3 = 75%"),
                text.col= c("blue", "red", "white", "black"))
grid()

Asimetria intuită când am constituit tabelul 'brev' este într-adevăr, vizibil mai pronunţată la curba albastră (pentru 'Română') decât la cealaltă. Cele două curbe au un singur punct de intersecţie, pe la abscisa 6.25; observând zonele subântinse de cele două curbe în stânga şi în dreapta acestei abscise, găsim că proporţia notelor sub 6.25 este de aproape două ori mai mare la "Matematică" decât la "Română", iar pentru notele mai mari ca 6.25 raportul tocmai menţionat se inversează.

Dar trebuie menţionat că modificând parametrul "adjust" al funcţiei density() (mai sus am ales 'adjust=1.75'), putem obţine curbe "mai dese"; de exemplu, pentru valoarea 0.3 obţinem curbele:

Graficul rezultat astfel nu contrazice observaţiile anterioare (exceptând faptul că acum, pe lângă intersecţia de abscisă 6.25, apar şi alte puncte de intersecţie). "Zimţii" ca de pieptene, care au apărut pe cele două curbe reflectă faptul că notele respectă un anumit şablon - de exemplu, apar note de 5.50, 5.55, 5.60, 5.65 etc., dar nu şi note "intermediare" (ca 5.51, 5.52, etc.); ne putem edifica asupra acestui aspect folosind funcţia unique(), prin care obţinem reprezentanţii unici ai notelor:

sort(unique(evn$Matematică))  # curba colorată cu roşu, pe graficul de mai sus
  [1]  1.00  1.50  1.65  1.70  1.90  1.95  2.00  2.05  2.10  2.15  2.20  2.25
 [13]  2.30  2.35  2.40  2.45  2.50  2.55  2.60  2.65  2.70  2.75  2.80  2.85
# ... ... ...
[169]  9.95 10.00  # numai 4.57, 6.66, 7.88 sunt în afara şablonului "din 5 în 5 sutimi" 

În cazul redat mai sus, exceptând intervalul [1, 1.50] şi încă trei valori (4.57, 6.66, 7.88) - secvenţa notelor respectă şablonul "din 5 în 5 sutimi" (conducând la forma zimţată a curbei).

5. Densităţile mediilor finale dintr-un interval dat, pe judeţe

Vrem să vizualizăm curbele care estimează densitatea distribuţiei mediilor finale cuprinse într-un interval dat, pe fiecare judeţ, în ordinea valorilor mediei judeţene pe acel interval. Selectăm din "evn" numai datele care ne interesează acum (valorile "Media" diferite de 'NA' şi câmpul "Jud") şi definim o funcţie în care mai întâi, restrângem datele corespunzător intervalului indicat şi apoi procedăm ca în [2] §9: creem un "factor" ordonat după mediile judeţene (pentru datele restrânse), ordonăm crescător vectorul din atributul 'scores' al acestuia şi instituim ca etichete pentru afişări alipirile de nume de judeţe şi respectiv, medii judeţene:

medii <- subset(evn, !is.na(Media), select = c("id", "Media", "Jud"))
dens_gap <- function(from = 1, to = 10) {
    medii <- subset(medii, Media >= from & Media <= to)
    rejud <- reorder(medii$Jud, medii$Media, FUN = mean, order = TRUE)
    medii$Jud <- factor(rejud, levels = levels(rejud), 
                        labels = paste(levels(rejud), 
                                       round(sort(attr(rejud,"scores")), 2)))
    medii <- medii[order(medii$Jud), ]
    require(lattice)
    densityplot(
        ~ Media | Jud,  # densitatea mediilor finale, pe fiecare judeţ în parte
        data = medii, 
        from = from, to = to, n = 128,  # pentru calculul densităţii
        type = "g", aspect = 0.5,  # grid, înălţime/lăţime=0.5
        par.strip.text = list(cex=0.8),  # reduce mărimea de scriere a etichetelor
        layout = c(7, 6),  # 6 linii a câte 7 panouri (42 judeţe)
        as.table = TRUE,  # de la stânga la dreapta, de sus în jos
        main = list(label = bquote("Medii ["*.(from)*","~.(to)*"] în judeţe)"),
                    font = 3, col = "cornsilk4")
    ) -> dst_list  # listă cu toate datele şi setările necesare trasării graficului
    plot(dst_list)
}

Apelul dens_gap() produce densităţile respective pentru intervalul implicit [1, 10]:

Pe măsură ce media judeţeană creşte - adică, vizând panourile de la stânga spre dreapta de sus în jos (sau mai simplu, parcurgând oricare coloană de panouri de sus în jos) - curba densităţii se aplatizează spre stânga (scade proporţia de medii mici, vizibil în special între gradaţiile 1 şi 4) şi se înalţă spre dreapta (culminând la "M.Bucureşti 7.41", unde înălţimea depăşeşte puţin 25%, la o abscisă apropiată de media 9). Observând înălţimea curbelor în punctul de abscisă 10, deducem că proporţia mediilor 10 în judeţele respective este cel mult 5% (în câteva cazuri fiind totuşi puţin peste 5%) şi doar la "M.Bucureşti 7.41" se înalţă până pe la 10%.
Să mai observăm că pe intervalul vizat [1, 10], mediile judeţene variază de-a lungul a două puncte (de la 'Giurgiu 5.39' la 'M.Bucureşti 7.41').

Apelul dens_gap(5) produce densităţile respective pentru intervalul [5, 10]:

Pentru intervalul vizat acum [5, 10], mediile judeţene variază de-a lungul unui punct - de la 'Harghita 6.95', la 'M.Bucureşti 7.93'.

Cu dens_gap(1, 5) obţinem densităţile respective pentru intervalul [1, 5], cu dens_gap(7.5, 9) pe cele corespunzătoare intervalului [7.50, 9], ş.a.m.d.

Într-un jurnal de ştiri am găsit acest aspect interesant: intervalul care cuprinde cele mai multe note este [8, 8.49]; pentru a releva acest aspect putem folosi funcţia cut(), clasificând mediile din jumătate în jumătate de punct:

medii$gap <- cut(medii$Media, breaks = c(seq(1, 9.5, by=0.5), 10.01), right=FALSE)
as.data.frame(sort(table(medii$gap))) -> tst
print(data.frame(tst[1:6, ], tst[7:12, ], tst[13:18, ]), row.names=FALSE)
     Var1 Freq   Var1.1 Freq.1   Var1.2 Freq.2
  [1,1.5)  823  [3.5,4)   6684  [6,6.5)  11083
  [1.5,2) 1276  [4,4.5)   7968  [6.5,7)  11441
  [2,2.5) 2276  [4.5,5)   9278  [7,7.5)  11835
  [2.5,3) 3624  [5,5.5)  10302  [7.5,8)  12541
  [3,3.5) 5121  [9,9.5)  10482  [8.5,9)  13290
 [9.5,10] 6344  [5.5,6)  10758  [8,8.5)  13522

Constatăm că într-adevăr, intervalul [8, 8.49] conţine cel mai mare număr de note (anume, 13522 note), dintre toate intervalele "din 0.5 în 0.5 puncte" (urmat de [8.50, 9) cu 13290 note şi de [7.50, 8) cu 12541 note). Cu dens_gap(8, 8.49) putem vizualiza curbele de densitate corespunzătoare pe judeţe; menţionăm doar că în acest caz, plaja mediilor judeţene (calculate pentru acest interval de note) măsoară numai 5 sutimi - de la 'Covasna 8.21' la 'Sălaj 8.26'.

Putem folosi dens_gap() şi pentru mediile V-VIII: redefinim 'medii', selectând 'Media8' în loc de 'Media' (şi păstrăm filtrarea valorilor 'NA', fiindcă în tabelul 'brev' din §2 vedem că există 5 candidaţi la care lipseşte valoarea mediei V-VIII); renumim 'medii$Media8' cu 'medii$Media' (evitând necesitatea unor modificări în corpul funcţiei dens_gap()) şi apoi putem apela funcţia de mai sus:

medii <- subset(evn, !is.na(Media8), select = c("id", "Media8", "Jud"))
names(medii)[2] <- "Media"  # redenumeşte "Media8", pentru a invoca direct dens_gap()
dens_gap(6)

Am optat pentru intervalul de note [6, 10] fiindcă dacă există medii generale V-VIII mai mici de 6, cu siguranţă ele sunt în număr foarte mic. Redăm doar primul si ultimul panou dintre cele 42 ale graficului (acestea fiind obţinute după includerea în funcţia densityplot() a parametrului subset = as.numeric(Jud) %in% c(1, 42), şi modificarea layout = c(2, 1) (pentru rând grafic cu două panouri)):

Vedem astfel că mediile judeţene V-VIII sunt cuprinse într-un interval de 60 de sutimi, de la 'Tulcea 8.33', la 'M.Bucureşti 8.94'. În general (observând toate cele 42 de panouri grafice), mediile judeţene de examen sunt cam cu 1-1.50 puncte mai mici decât mediile judeţene V-VIII.

vezi Cărţile mele (de programare)

docerpro | Prev | Next