momente şi schiţe de informatică şi matematică
anti point—and—click

Populaţie şi şcoli (de la INSse la Python şi R)

Python | R
2020 apr

Câtă populaţie şi câte şcoli are fiecare judeţ ? Cât de adevărat ar fi, că numărul de locuitori influenţează „liniar” numărul de şcoli ?

Sunt oare importante (sau… corect puse), aceste întrebări ⁈ Mai degrabă nu; dacă ar fi, probabil că am găsi uşor datele brute necesare investigării—dar constatăm că nu prea este cazul şi pentru a obţine totuşi numerele necesare, avem de conceput anumite programe „intermediare” (iar aceasta, împreună cu ce-i prin jur, face ca şi lucrurile neimportante să devină măcar interesante).

Căutarea firească a datelor brute

Pe //data.gov.ro găsisem în 2016 – în setul de date "Demografie" – fişierul Excel "populatia-romaniei-pe-judete.xls" (descris cu tot dichisul şi exploatat în [1]); un asemenea fişier (dar actualizat pentru 2019), ar fi necesar acum. Astăzi însă (prin martie 2020) – n-am reuşit să mai dau de vreun set de date demografice.

În schimb, am remarcat desigur noua „faţă” a site-ului; se foloseşte ca şi mai înainte, CKAN ("a tool for making open data websites", bazat intern pe Python, Apache, PostgreSQL şi pe JDK) şi se folosesc bineînţeles Bootstrap, jQuery şi Jinja2. Lucrurile arată nemaipomenit, totuşi decent, dar nicidecum simplu şi sobru ca pentru sarcina de a furniza celor interesaţi seturi de date; pe de altă parte, probabil că "open data websites" implică şi asigurarea către terţi (instituţii) a posibilităţii de "upload"-are directă a seturilor de date (în diverse formate) – încât o asemenea complexitate a instrumentării lucrurilor poate să pară justificată (cu atât mai mult desigur, cu cât se permit şi comentarii din parte utilizatorilor).

Pe //insse.ro/cms/ (Institutul Naţional de Statistică) am avut mai mult noroc: deja pe pagina de bază găsim afişat mare "Populaţia rezidentă 19,4 mil. persoane" (la 1 ianuarie 2019); însă nu am reuşit imediat, să găsesc şi vreo statistică pe judeţe.

În schimb, am remarcat ştirea despre conferinţa internaţională Use of R in Official Statistics, în care se recunoaşte că „R a devenit «lingua franca» pentru statisticieni, metodologi și cercetători de date din întreaga lume”. Această ştire ar trebui să încurajeze şi la noi, preocupările faţă de R; dar ştiind aplecarea obstinentă (şi în învăţământ şi în instituţii) faţă de produsele Microsoft şi tehnologia "point-and-click" – estimez că la noi, conferinţa respectivă nu va avea cine ştie ce urmări…

Totuşi, pe de o parte, am identificat pagina Site-urile Direcțiilor Regionale și Județene de Statistică, pe care se precizează că site-ul asociat unui judeţ se poate accesa prin www.[nume_judet].insse.ro. De exemplu, //vaslui.insse.ro/ ne conduce la site-ul "Direcţia Judeţeană de Statistică VASLUI" – iar aici vedem imediat (afişat mare) ceea ce ne interesează: "Populaţia la 01.01.2019" este 495779.

Pe de altă parte, am identificat link-ul TEMPO-ONLINE („baze de date statistice”) şi acesta iată că furnizează datele dorite, anume cam aşa:

Tabelul furnizat poate că este frumos colorat (şi nobil liniat), dar altfel este oribil (îi lipseşte doar celebra coloană iniţială "Nr.Crt.", ca să fie cotat ca funcţionăresc); un tabel de vreo 60 de rânduri merită prezentat mai simplu şi mai convenabil. Dar nu „cum arată” ne interesează, ci cum putem extrage datele din acest tabel (spre a le folosi eventual în programele proprii; preferam desigur, să ni se dea un fişier CSV simplu).

Dar să nu ne grăbim să extragem de aici datele… Pe linia judeţului Vaslui găsim în coloana "Numar persoane" valoarea 373863 – ori aceasta diferă semnificativ faţă de aceea menţionată la "Populaţie" pe site-ul //vaslui.insse.ro (valoare pe care am redat-o mai sus). Diferenţele între valorile din tabel şi cele furnizate de site-urile judeţene s-ar explica prin faptul că "Numar persoane" exclude persoanele care au emigrat, în timp ce "Populaţie" vizează toate persoanele cu domiciliul în judeţ.

Prin urmare, avem de ales între cele două categorii de date - cele furnizate prin tabelul redat parţial mai sus şi respectiv, cele pe care le putem prelua de pe site-urile //nume_judet.insse.ro. Dar nu ne gândim aici, la o alegere „ştiinţifică”, ci doar la aceea care este mai interesantă pentru programare: vom accesa cele 42 de site-uri judeţene, preluând de pe fiecare valoarea înscrisă la "Populaţie".

Extragerea datelor din tabele „invizibile”

Ca să accesăm site-urile //nume_judet.insse.ro avem nevoie de "nume_judet"; nu putem obţine împreună aceste valori, prin procedeul obişnuit (sau vulgar): selectăm cu mouse-ul zona respectivă din pagina "Site-urile Direcțiilor Regionale..." afişată în fereastra browserului şi o copiem prin Ctrl-C – am obţine doar un text precum "Arad, Caras-Severin, Timis, Hunedoara", ori "nume_judet" nu este "Caras-Severin" ci "carasseverin" (şi nu este "Municipiul Bucuresti", ci "bucuresti").
Ca să obţinem link-urile respective trebuie să vizualizăm sursa paginii, folosind Ctrl-U; cel mai probabil, ele sunt conţinute într-un element <table> (nu-i obligatoriu; ar putea fi înscrise şi într-un <div>) – şi căutând prin Ctrl-F descoperim că există un <table> (şi numai unul singur) care într-adevăr, conţine link-urile pe care le vrem.

Tabelul HTML depistat astfel este „stufos”, specificând şi stilări sau atribute care nu ne interesează (precum <td style="background-color:rgb(...)"><a class="link"...); în loc să-l „periem” în gedit cum am face de obicei, preferăm să folosim Python astfel:

from bs4 import BeautifulSoup as BS
import requests, csv
f = csv.writer(open("nume_jud.csv", "w"))
f.writerow(["Jud", "Link"])  # antetul fişierului CSV
url = 'https://insse.ro/cms/ro/content/' + \
      'site-urile-directiilor-regionale-si-judetene-de-statistica'
req = requests.get(url, headers = {'User-Agent': "and I'm a Browser"})
page = BS(req.content, 'html.parser')
tabel = page.find('table')
links = tabel.find_all('a')
for lk in links:
    f.writerow([lk.get_text(), lk.get('href')])

Încărcăm în gedit fişierul CSV rezultat prin execuţia acestui program şi eliminăm toate porţiunile "http://www." şi ".insse.ro" (acestea încadrează "nume_judet"); deasemenea, corectăm denumirile judeţelor unde este cazul – se ştie că de obicei, pe site-urile instituţiilor publice lipsesc „diacriticele” (în plus, spaţiul este folosit defectuos: de exemplu apare "Nord - Est", dar apoi fară spaţiu "Sud-Vest").
Redăm fişierul final "nume_jud.csv", dar formatându-l aici pe coloane (după regiuni, pentru care am păstrat numai denumirea):

Jud,Link
Nord-Est,           Sud-Est,              Sud,                  Sud-Vest,
Botoşani,botosani   Brăila,braila         Argeş,arges           Dolj,dolj
Neamţ,neamt         Buzău,buzau           Călăraşi,calarasi     Gorj,gorj
Iaşi,iasi           Galaţi,galati         Dâmboviţa,dambovita   Mehedinţi,mehedinti
Bacău,bacau         Tulcea,tulcea         Giurgiu,giurgiu       Olt,olt
Suceava,suceava     Constanţa,constanta   Ialomiţa,ialomita     Vâlcea,valcea
Vaslui,vaslui       Vrancea,vrancea       Prahova,prahova
                                          Teleorman,teleorman

Vest,                        Nord-Vest,                 Centru,           Bucureşti,
Arad,arad                    Bihor,bihor                Alba,alba         Ilfov,ilfov
Caraş-Severin,carasseverin   Bistriţa-Năsăud,bistrita   Braşov,brasov     Bucureşti,bucuresti
Timiş,timis                  Cluj,cluj                  Covasna,covasna
Hunedoara,hunedoara          Maramureş,maramures        Harghita,harghita
                             Satu-Mare,satumare         Mureş,mures
                             Sălaj,salaj                Sibiu,sibiu

Printr-un program analog celui redat mai sus, vom accesa site-urile judeţene sintetizate în "nume_jud.csv" – urmând să vedem cum „extragem” de pe fiecare, valoarea pentru "Populatie".

Populaţia din fiecare judeţ

Accesând aleatoriu câteva site-uri //judet.insse.ro, conchidem că toate aceste 42 de site-uri au aceeaşi structură – bazată pe un acelaşi fişier-şablon, reflectând pe de o parte un conţinut comun tuturor judeţelor şi pe de altă parte, un conţinut adaptat judeţului. Pe toate, apare o „linie” cu trei coloane precum:

sau, pentru încă o exemplificare:

Ne interesează numai prima coloană ("populaţia după domiciliu") – ultimele două coloane variază după judeţ, în privinţa categoriei statistice reflectate. Majoritatea site-urilor reflectă populaţia de la 01.01.2019 şi probabil (neintenţionând o analiză sociologică riguroasă), o să ignorăm faptul că pe unele au fost făcute actualizări, vizând 01.07.2019 sau chiar "1 ian 2020".

Desigur, datele sunt redate încântător (ne mirăm că purtând mouse-ul pe prima coloană nu şi auzim ceva melodie "PO-pu-LAT-i-AA-aa--AA---LAaa- ..."); ce obiecte să fie, cele trei coloane? Sunt "link"-uri, cum vedem purtând mouse-ul pe fiecare şi observând "bara de stare" a browser-ului; şi anume, toate trei sunt link-uri identice, spre pagina care le cuprinde (https://[nume_judet].insse.ro/#) – cu una sau poate mai multe excepţii, când click-ul trimite la o pagină "OOPS! THAT PAGE CAN’T BE FOUND."

Vizualizând sursa paginii putem identifica <section class="frontpage-statistical"> care conţine o diviziune <div class="container"> care conţine o diviziune <div class="row"> care conţine trei diviziuni <div class="col-sm-12 col-md-12 col-lg-4" style="margin: 0.5em 0; z-index: 999;"> – câte una pentru fiecare coloană; fiecare coloană conţine un link <a class="btn btn-primary text-center" style="width: 100% !important; margin: 1.5em auto 0 auto !important; height: 10em !important;" href="#"> şi apoi, un element ca <h3 class="text-center">Populația la 01.01.2019</h3> urmat de un element ca <h1 class="text-center" style="font-size: 4em !important;">495779</h1>.
Esenţiale sunt cele două elemente (ultimul, de fapt) cu class="text-center", din fiecare coloană. Celelalte aşa de numeroase specificaţii (pe care m-am străduit să le redau mai sus), au un rol decorativ dubios, ţinând mai degrabă de spectacol şi de lipsa bunului simţ, decât de vreo sugestie de „ştiinţă”.

Următorul program Python mapează fişierul "nume_jud.csv" pe un dicţionar având drept chei "Jud" şi "Link" (din antetul CSV); pentru fiecare linie pentru care valoarea din "Link" nu este nulă (deci pentru judeţe, nu şi pentru regiuni) se accesează (folosind requests.get()) site-ul corespunzător acelui link (folosind un nume fictiv de "User-Agent", fără de care accesul ar fi interzis) şi folosind metode din BeautifulSoup, se extrage de pe pagina respectivă conţinutul primului dintre elementele <h1 class="text-center" ...>, scriind rezultatul obţinut într-un fişier CSV:

from bs4 import BeautifulSoup as BFS
import requests, csv
readNJ = csv.DictReader(open("nume_jud.csv"))
f = csv.writer(open("jud_pop.csv", "w"))
f.writerow(["Jud", "Pop"])
for row in readNJ:
    Jud = row['Jud']; Link = row['Link']
    if Link:  # ignorăm "regiunea"
        url = 'https://'+ Link +'.insse.ro/'
        req = requests.get(url, headers = {'User-Agent': "and I'm a Browser"})
        page = BFS(req.content, 'html.parser')
        Pop = page.find('h1', class_='text-center')  # numai primul <h1>
        if Pop:  # la Hunedoara lipseşte <h1 class="text-center">
            f.writerow([Jud, Pop.contents[0]])
        else:
            f.writerow([Jud, ''])

Redăm parţial fişierul "jud_pop.csv" rezultat, constatând anumite gafe şi greşeli:

Jud,Pop
Botoşani,454919
Neamţ,566870
Iaşi,<h2>956.216 persoane</h2>  # H2 conţinut în H1 (??)
...
Brăila,<h1>342663 pers.</h1>
Buzău,462942 pers.
...
Timiş,<h4> noiembrie 2019</h4>
Hunedoara,
Bihor,616264
Bistriţa-Năsăud,COVID-19
...
Sibiu,"<font color=""Coral""><h2>467.568 persoane</h2></font>"
...

Este o dovadă de proastă învăţare, să pui drept conţinut al unui element <h1> un element de aceeaşi natură – cum vedem pe liniile Iaşi, Brăila, Timiş, Sibiu redate mai sus (e drept că browserele nu te reclamă în mod explicit pentru gafele tale din codul HTML).

Bineînţeles că vom elimina de peste tot "persoane" şi "pers.". Pe de altă parte, constatăm că avem de corectat cumva (inspectând eventual chiar site-ul judeţului) în trei locuri: la Timiş nu este indicată populaţia, ci "Efectivul salariaţilor" (pentru "noiembrie 2019"); la Hunedoara nu există elemente <h1> (încât "Pop" a rămas necompletat); iar la Bistriţa-Năsăud primul element <h1> (dintre cele nouă aliniate acolo) are conţinutul "COVID-19" şi populaţia apare pe al treilea rând de elemente <h1>.

După corecturi (căutând eventual prin alte locuri, datele lipsă) – vom folosi fişierul "jud_pop.csv" ceva mai încolo, încercând mica analiză statistică propusă la început (vizând populaţia şi numărul de şcoli din judeţe).

Şcolile din fiecare judeţ

Mai avem nevoie de un fişier "jud_scl.csv", conţinând numărul de şcoli pe fiecare judeţ. Avem noroc (desigur): pe //data.gov.ro găsim un fişier CSV de 4.4MB denumit "datedeschise-specializari-2019-2020.csv" – o să-l redenumim "scoli19.csv" – care listează toate unităţile şcolare „subordonate Ministerului Educației Naționale care funcționează în anul școlar 2019-2020”. N-avem decât să numărăm liniile din acest fişier care corespund unui aceluiaşi judeţ (repetând pentru fiecare judeţ); secvenţa de operaţii necesară poate fi modelată foarte uşor în R, astfel:

vb@Home:~/20mar$ R -q   # sesiune interactivă de lucru (workspace) în R
> scoli19 <- read.csv("scoli19.csv")
> jud_scl <- as.data.frame(table(scoli19$Judet))
> write.csv(jud_scl, file="jud_scl.csv", row.names=FALSE)
> q()
Save workspace image? [y/n/c]: y
vb@Home:~/20mar$

Precizăm că dacă este cazul, se pot obţine uşor lămuriri pentru oricare comandă, tastând de exemplu help(data.frame), sau help(table), sau ?table (sau ?read.csv, etc.).

Redăm primele câteva linii din fişierul "jud_scl.csv" rezultat, pentru a lămuri ce mai avem de făcut:

"Var1","Freq"
"ALBA",803
"ARAD",782
"ARGE?",1121
"BAC?U",1188
"BIHOR",1256
"BISTRI?A-N?S?UD",640

Încă se folosesc majuscule pentru denumiri (obicei general în instituţiile noastre, adoptat pe când se născuse problema semnelor "diacritice"); literele româneşti apar ca "?". Dar nu luăm mai mult seama, la niciunul dintre aceste aspecte: avem deja denumirile corecte ale judeţelor, în fişierul "jud_pop.csv"; n-avem decât să ordonăm la fel (după judeţ) cele două fişiere şi să păstrăm din "jud_scl.csv" numai coloana desemnată prin "Freq" (reprezentând numărul de şcoli din judeţ).

Populaţie, şcoli… şi greşeli

Reconstituim sesiunea de lucru în R anterioară şi ordonăm „la fel” (după coloana "Jud" şi respectiv "Var1") cele două fişiere:

vb@Home:~/20mar$ R -q
> jud_pop <- read.csv("jud_pop.csv")
> write.csv(jud_pop[order(jud_pop$Jud), ], file="jud_pop.csv", row.names=FALSE)
> write.csv(jud_scl[order(jud_scl$Var1), ], file="jud_scl.csv", row.names=FALSE)

Deschizând în gedit cele două fişiere ordonate rezultate şi urmărindu-le „în paralel”, constatăm că în "jud_scl.csv" avem de inversat numai două sau trei linii (de exemplu, linia "MUNICIPIUL BUCURE?TI",1491 trebuie mutată după linia "BRA?OV",798).
Rămâne să „mutăm” coloana "Freq" din "jud_scl.csv" în "jud_pop.csv" şi folosim iarăşi R (nu fiindcă n-ar fi simplu de făcut şi altfel):

vb@Home:~/20mar$ R -q
> rm(list=ls())  # elimină datele păstrate din sesiunea anterioară
> jud_pop <- read.csv("jud_pop.csv")
> jud_scl <- read.csv("jud_scl.csv")
> jud_pop$Scl <- jud_scl$Freq  # adaugă coloana Freq (numind-o Scl) din jud_scl
> write.csv(jud_pop, file="judPopScl.csv", row.names=FALSE)
> q()
Save workspace image? [y/n/c]: n
vb@Home:~/20mar$ 

Am preferat să scriem rezultatul în fişierul "judPopScl.csv"; de-acum nu vom mai avea nevoie de cele două fişiere CSV din care a provenit acesta.

Să listăm acum, datele grupate mai sus (judeţ, populaţie, şcoli):

vb@Home:~/20mar$ R -q
> jPS <- read.csv("judPopScl.csv") 
> jPS
               Jud     Pop  Scl
1             Alba  375586  803
2             Arad  470846  782
...           ...   ...     ...
21            Gorj  362610  798
22        Harghita   71070  829
23       Hunedoara  464109  731
...           ...   ...     ...
42         Vrancea  388495  713
>

De ce ne-am temut (dat fiind modul în care am obţinut datele) nu am scăpat: pe linia 22 vedem că în judeţul Harghita ar fi numai 71070 de locuitori, ceea ce nu poate fi adevărat! Verificând //harghita.insse.ro vedem că numărul respectiv reprezintă doar "Efectiv salariaţi - dec 2019" (nu populaţia "cu domiciliul în judeţ", ca în majoritatea celorlalte cazuri – cum încă sperăm…).
Ce greşală! De fapt, greşala constă în faptul că n-am observat la timp situaţia constatată acum (trebuia observată atunci când am produs fişierul "jud_pop.csv").

Acuma, în loc să mai caut ce populaţie trebuie să înscriu pe linia 22 la Harghita (şi poate, să mai verific şi datele de la celelalte judeţe), am revenit într-o doară pe TEMPO-ONLINE ("baze de date statistice") şi am avut surpriza să văd că acum apare şi "POPULATIA DUPA DOMICILIU" (itemul de meniu "POP107A") şi încă, inclusiv la data de 01.01.2020 (la momentul când începusem să formulez cele de mai sus, apărea numai "populaţia rezidentă" şi numai până în 2019).

Nu mai descriem aici (nu că n-ar fi interesant), cum am „extras” datele furnizate prin "POP107A" (anume, cele pentru anul 2020) şi nici cum le-am înscris în locul celor existente în coloana "Pop" a fişierului "judPopScl.csv".

Şcoli vs. Populaţie (demersuri statistice)

Să listăm acum, datele grupate mai sus (judeţ, populaţie la 01.01.2020, şcoli în anul şcolar 2019-2020):

vb@Home:~/20mar$ R -q
> jPS <- read.csv("judPopScl.csv") 
> jPS
               Jud     Pop  Scl
1             Alba  373710  803
2             Arad  470093  782
3            Argeş  632124 1121
4            Bacău  737629 1188
5            Bihor  616264 1256
6  Bistriţa-Năsăud  327886  640
7         Botoşani  456130  956
8           Brăila  340804  500
9           Braşov  638000  798
10       Bucureşti 2151665 1491
11           Buzău  461039  899
12        Călăraşi  306820  497
13   Caraş-Severin  315473  714
14            Cluj  736945 1053
15       Constanţa  763549  982
16         Covasna  225743  511
17       Dâmboviţa  518645  848
18            Dolj  686350 1129
19          Galaţi  628104  842
20         Giurgiu  269221  458
21            Gorj  355358  798
         Jud    Pop  Scl
22  Harghita 330473  829
23 Hunedoara 453431  731
24  Ialomiţa 285100  424
25      Iaşi 965634 1473
26     Ilfov 451839  444
27 Maramureş 520551  940
28 Mehedinţi 277017  542
29     Mureş 590305 1142
30     Neamţ 566080  880
31       Olt 431248  837
32   Prahova 787529 1256
33     Sălaj 243525  530
34 Satu-Mare 386110  630
35     Sibiu 469285  725
36   Suceava 764123 1280
37 Teleorman 365976  570
38     Timiş 758380 1024
39    Tulcea 234243  406
40    Vâlcea 394725  718
41    Vaslui 504879 1003
42   Vrancea 382688  713

Am listat în ordinea alfabetică a valorilor "Jud" (uşurând astfel, eventuale confruntări cu alte surse); am redat lista pe două coloane, pentru a facilita observarea şi compararea datelor – "coloanele" fiind obţinute prin jPS[1:21, ] şi respectiv jPS[22:42, ].

Însă lucrurile pot deveni interesante ordonând nu alfabetic, ci după "Pop", sau după "Scl". Desigur, putem observa uşor şi direct pe lista redată pe două coloane mai sus, anumite discrepanţe – de exemplu, vedem că la Bihor şi Braşov populaţia diferă puţin (cu vreo 20000), dar numărul de şcoli diferă mult (cu vreo 350); probabil, explicaţia acestei „discrepanţe” ţine cel mai mult de relief (pare clar că pe o zonă muntoasă apar mai puţine construcţii decât pe o zonă întinsă). La fel poate, despre Tulcea şi Sălaj – au cam aceeaşi populaţie, dar diferenţă de vreo sută de şcoli; analog putem exersa pentru alte câteva cazuri.

Dar metoda cea mai bogată în sugestii, dacă e vorba de a confrunta una cu alta şi a sesiza anumite şabloane, constă în a reprezenta grafic mărimile respective.

Dar să vedem întâi cam ce fel de mărimi avem:

> summary(jPS[, c(2,3)])  # [, c(2,3)] selectează din toate liniile, coloanele 2 şi 3
      Pop               Scl        
 Min.   : 225743   Min.   : 406.0  
 1st Qu.: 344442   1st Qu.: 632.5  
 Median : 458584   Median : 816.0  
 Mean   : 527969   Mean   : 842.0  
 3rd Qu.: 631119   3rd Qu.:1018.8  
 Max.   :2151665   Max.   :1491.0

> boxplot(jPS$Pop, horizontal=TRUE, cex.axis=0.8)

Valoarea cea mai mică de pe coloana "Pop" corespunde judeţului Covasna, iar cea mai mare – municipiului Bucureşti; aproximativ un sfert dintre cele 42 de linii au valoarea din "Pop" mai mică decât "prima quartilă" şi un alt sfert – mai mare decât a treia quartilă; celelalte 50% valori sunt reprezentate pe graficul furnizat de boxplot() printr-un dreptunghi care (pentru redarea orizontală) are latura stângă la gradaţia dată de "1st Qu" (să notăm Q1) şi latura din dreapta la abscisa "3rd Qu." (fie Q2); bara figurată în interiorul dreptunghiului are abscisa dată de mediana valorilor.
Liniile punctate care „prelungesc” dreptunghiul spre stânga şi spre dreapta se întind până la cea mai mare dintre valorile "Min." şi Q1-1.5*IRQ, respectiv până la min("Max.", Q3+1.5*IRQ) – unde IRQ este lungimea dreptunghiului (=Q3-Q1). Valorile rămase în afară (cum este cazul valorii "Max.", reprezentate prin cerculeţul din dreapta ultimei gradaţii de pe graficul redat mai sus) sunt considerate „excepţii” faţă de tendinţa principală a distribuţiei datelor.

Vom elimina linia 10, constatând că "Bucureşti" este „excepţie” (sau "outlier"); pe de altă parte, "Iaşi" (pe linia 25) nu este un "outlier", dar are valoarea "Pop" foarte apropiată de limita dreaptă a "boxplot"-ului şi în plus, dacă ne uităm şi la coloana "Scl" – observăm aceeaşi situaţie de apropiere. Am spune atunci, că Bucureşti şi judeţul Iaşi sunt oarecum „excepţionale” în raport cu celelalte, atât ca "Pop" cât şi ca "Scl"; ignorându-le, se poate constata că pentru cele 40 de linii rămase obţinem un "boxplot" mai simetric (sau echilibrat), atât pentru "Pop" cât şi pentru "Scl".

Programul următor (lansat tastând source("work.R"), în sesiunea R curentă) produce scatterplot-ul redat dedesubt:

jPS2 <- jPS[-c(10,25), ]  # exclude Bucureşti ("outliner") şi Iaşi
with(jPS2, 
     { plot(Scl ~ Pop, 
            pch=19, cex=0.5, cex.axis=0.8, col="blue", 
            xlab="Populaţie", ylab="Şcoli")
       text(Scl ~ Pop,
            labels=rownames(jPS2),
            cex=0.7, pos=3, offset=0.2)
     })

Cele 40 de judeţe reţinute în jPS2 sunt reprezentate prin puncte etichetate (prin funcţiile text() şi rownames()) cu rangul liniei corespunzătoare în jPS şi amplasate după valorile aferente lor din coloanele "Pop" şi "Scl" (dar nu cu aceeaşi unitate de măsură).

Uitându-ne pe verticală, sau pe orizontală, putem observa „discrepanţe” (mult mai uşor, decât exersam mai sus direct pe tabel); de exemplu, 19 şi 3 (adică Galaţi şi Argeş) au cam aceeaşi populaţie (punctele sunt cam pe o aceeaşi verticală, ceva mai la dreapta abscisei 6e+05=600000), dar diferă mult ca număr de şcoli (ordonatele punctelor fiind aproximativ 800 şi 1100); la fel 37 şi 1, sau 6 şi 22. Iar 7 (Botoşani) şi 15 (Constanţa) nu diferă mult pe verticală (adică au cam acelaşi număr de şcoli), dar diferă mult pe orizontală (ca populaţie); la fel de exemplu, 1 şi 9, etc. Avem şi perechi de judeţe apropiate şi ca "Pop" şi ca "Scl": 1 (Alba) şi 21 (Gorj), etc.

Dar să încercăm să vedem şi în ansamblu, „micşorând” cumva lucrurile (mai ales pe verticală, unde gradaţiile exprimă abia sute şi nu 105 cum avem pe orizontală) şi „netezind” sau ignorând unele asperităţi (precum punctul 26 (Ilfov) care este clar izolat, faţă de „norul” de deasupra lui):

Linia roşie punctată sintetizează cumva „norul” punctelor (şi exprimă tendinţa principală a acestora), fiind obţinută printr-un proces de înlocuire din aproape în aproape a câte unui grup de puncte vecine prin câte un anumit punct considerat (în baza unui aceluiaşi algoritm de aproximare) ca reprezentativ pentru grup, unind apoi aceşti reprezentanţi; graficul este identic cu cel anterior, dar micşorat pe verticală (trăgând cu mouse-ul), iar linia roşie a fost calculată şi produsă înlocuind în programul redat mai sus funcţia plot() cu funcţia scatter.smooth().

În principiu, linia roşie este astfel calculată încât suma abaterilor pe verticală ale punctelor de deasupra liniei se compensează cu suma abaterilor de la linie ale punctelor aflate dedesubt (şi putem să ne dăm seama „din ochi” că „asperitatea” 26 semnalată mai sus contribuie totuşi la această compensare).

„Netezind” lucrurile şi mai mult, am zice că linia roşie este formată din două „drepte”; prima dreaptă merge până pe la punctul 27 (unde "Pop" este puţin peste 105), iar a doua continuă din acest punct şi are panta ceva mai mică decât prima. Interpretarea ar fi aceasta: până la 500000 de locuitori, "Scl" depinde liniar (crescător) de "Pop"; de la 500000 încolo, "Scl" depinde iarăşi liniar de "Pop", dar cu o rată de creştere ceva mai mică decât în primul caz.

Desigur, putem determina aproximativ pantele respective (raportând diferenţa de şcoli la diferenţa de populaţie, pentru câte două puncte suficient de apropiate de dreptele respective) – pentru a vedea dacă n-am greşi prea mult „unificând” cumva, cele două drepte – scurtând prin "Scl" depinde liniar de "Pop":

Aici, graficului anterior (dar acum, mărit puţin pe verticală) i-am adăugat dreapta colorată cu verde, care „unifică” cele două drepte roşii – ca şi cum am trage puţin de capetele fiecăreia, în direcţii contrare; „se vede” pe figură că abaterile celor două drepte roşii de la dreapta verde se compensează, deci dreapta verde „unifică” într-adevăr, cele două tendinţe.

Desigur, n-am „tras cu mouse-ul” (de capetele dreptelor roşii), ci am folosit funcţia lm(), prin care se constituie un „model liniar” pentru variabilele indicate:

vb@Home:~/20mar$ R -q
> source("work.R")  # produce scatterplot-ul
> LM <- with(jPS2, lm(formula = Scl ~ Pop))
> abline(LM)  # adaugă "dreapta verde" (dreapta de regresie LM)

Să vedem cât de acceptabil este acest model liniar:

> summary(LM)

Residuals:
    Min      1Q  Median      3Q     Max 
-334.65  -90.48   14.40  104.84  267.91 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.031e+02  6.246e+01   3.252  0.00241 ** 
Pop         1.274e-03  1.237e-04  10.296 1.51e-12 ***
---
Signif. codes:  0***0.001**0.01*0.05.0.1 ‘ ’ 1

Residual standard error: 130.7 on 38 degrees of freedom
Multiple R-squared:  0.7361,	Adjusted R-squared:  0.7292 
F-statistic:   106 on 1 and 38 DF,  p-value: 1.507e-12

"Residuals" tabelează quartilele abaterilor verticale ale punctelor reale faţă de punctele de pe verticalele repective care sunt rezidente pe drepta de regresie; reziduurile pot fi plotate (toate) prin plot(jPS2$Pop, residuals(LM), ...); abline(0,0). Cea mai bună situaţie ar fi aceea în care "Min" şi "Max" n-ar diferi „prea mult” în valoare absolută şi "Median" ar fi cât mai mică (spre zero); cazul nostru parcă nu prea este în „cea mai bună” situaţie (dar plotând, vedem că totuşi, reziduurile par distribuite aleatoriu şi nu în jurul vreunei anumite curbe).

Tabelul "Coefficients" ne spune mai multe. De pe coloana "Estimate" avem ecuaţia Scl = 203.1 + 0.001274*Pop (dreapta de regresie); aceasta ne permite să estimăm cam câte şcoli s-ar potrivi pentru o anumită valoare "Pop" (care are sens: în mod firesc, "Pop" este de ordinul 105 sau cel puţin, de ordinul 104). De exemplu, pentru Pop=500000 am găsi Scl = 203+127.4*5 = 840 de şcoli (şi dacă ţinem cont de "St. error" – până la cel mult 840+"St. error"=840+62.46=900 şcoli, aproximativ).

Putem obţine direct estimarea de bază, prin funcţia predict(); de exemplu, pentru valorile "Pop" = k*105 (k de la 1 la 9):

> round(predict(LM, data.frame(Pop = (1:9)*100000)))
   1    2    3    4    5    6    7    8    9  # *100000 (Pop)
 330  458  585  713  840  967 1095 1222 1350  # estimare Scl (prin LM)

Pentru k=5 (adică Pop=500000) regăsim estimarea 840 pe care am calculat-o direct, mai sus.

Celelalte elemente din tabelele furnizate de summary() ("RSError", "R-squared", "p-value", etc.) măsoară acurateţea potrivirii datelor cu modelul "LM"; în principiu, cu cât "R-squared" este mai apropiat de 1 şi "p-value" este foarte mic – cu atât mai bună este estimarea asigurată de "LM" pentru dependenţa unei variabile faţă de cealaltă.

Ceva mai precis, fiindcă "p-value" este sub 0.001, iar "R-squared" este aproximativ 0.74 – putem zice că "Scl" depinde liniar de "Pop", dar numai circa 74% din variaţia variabilei "Scl" se explică liniar prin factorul "Pop"; restul de 26% trebuie „explicat” eventual prin alţi factori (de exemplu, prin forma de relief specifică judeţului).

docerpro | Prev | Next