Analiza unui model de regresie liniară folosind R cu lavaan și sem (2/4)

După ce ieri am aflat cum putem folosi Amos pentru a analiza o regresie liniară, astăzi, în al doilea episod, vom folosi limbajul R, împreună cu pachetele „lavaan” și „sem”. Deja pe trei sferturi știm ce avem de făcut. Vom crea un nou fișier de cod, apoi vom încărca structura de date și, pentru a debuta cu pachetul „lavaan”, încărcăm și librăriile cunoscute:

> load("D:/Baze de date/002 - Regresii liniare/Regresii.RData")
> library(lavaan)
> library(semPlot)

Urmează specificarea modelului. Am văzut că „lavaan” respectă standardele R iar corelațiile și covarianțele au fost specificate printr-o tildă dublă, însă în cazul nostru avem și efecte directe, fiind necesară introducerea unei noi reguli:

  • Efectele directe se specifică folosind o singură tildă (~). Pare destul de evident ca săgeata simplă să fie reprezentată cu o singură tildă iar săgeata dublă cu două tilde, așa că n-ar trebui să avem probleme cu reținerea acestor simboluri.

Dacă am aflat cum se reprezintă efectul direct, este deja o joacă să specificăm modelul:

> modelul <- ' # Specificarea ecuatiei de regresie
               vanzari ~ inteligenta + memorie + cunostinte
               # Specificarea covariantelor
               inteligenta ~~ memorie
               inteligenta ~~ cunostinte
               cunostinte ~~ memorie
             '

Desigur, aceasta nu este o specificare completă. Ar trebui să includem și covarianțele celor patru variabile observate cu ele însele pentru a îndeplini criteriul de identificare și ar trebui să stabilim metrica erorii de estimare. Din fericire, de metrică se ocupă automat „lavaan” iar varianțele unice le putem include automat direct în analiză:

> analiza <- lavaan(model = modelul, data = regresii, auto.var=TRUE)
 Warning message:
 In lav_data_full(data = data, group = group, cluster = cluster, :
 lavaan WARNING: some observed variances are (at least) a factor 1000 times larger than others; use varTable(fit) to investigate

Iată că „lavaan” își dă seama că ceva nu este în regulă cu această analiză, lucru de care ne-am dat seama și noi încă de la prelucrarea folosind Amos. Era mult prea mare varianța unică a criteriului. Într-adevăr, „lavaan” observă că varianța unei variabile este de peste 1000 de ori mai mare față de varianțele celorlalte variabile, prin urmare ceva se întâmplă. Mai mult, ne invită să și investigăm, prin utilizarea funcției „varTable()”:

> varTable(analiza)
    name     idx nobs type   exo user      mean          var  nlev lnam
 1 vanzari     1 102 numeric   0    0  6797.902  18027855.555    0
 2 inteligenta 2 102 numeric   0    0    10.738       7.444      0
 3 memorie     3 102 numeric   0    0    71.021    1006.471      0
 4 cunostinte  4 102 numeric   0    0    46.833     295.994      0

N-am putut refuza oferta și am constatat, într-adevăr, o varianță extrem de mare a criteriului. Acest lucru nu poate însemna decât că asumpția normalității distribuției criteriului nu a fost îndeplinită. Fie există agenți de vânzări care încasează extrem de mult în comparație cu ceilalți, fie unii dintre ei iau salariul degeaba pentru că nu reușesc să vândă nici măcar de un leu. O analiză univariată a variabilelor respective, iată, este obligatorie înainte de a trece la lucruri mai complicate. Oricum, această situație nu ne împiedică să afișăm rezultatele și diagrama de cale:

> summary(analiza, fit.measures=TRUE)
lavaan (0.5-23.1097) converged normally after 123 iterations

Number of observations 102

Estimator ML
 Minimum Function Test Statistic 0.000
 Degrees of freedom 0

Model test baseline model:

Minimum Function Test Statistic 247.785
 Degrees of freedom 6
 P-value 0.000

User model versus baseline model:

Comparative Fit Index (CFI) 1.000
 Tucker-Lewis Index (TLI) 1.000

Loglikelihood and Information Criteria:

Loglikelihood user model (H0) -2050.316
 Loglikelihood unrestricted model (H1) -2050.316

Number of free parameters 10
 Akaike (AIC) 4120.631
 Bayesian (BIC) 4146.881
 Sample-size adjusted Bayesian (BIC) 4115.295

Root Mean Square Error of Approximation:

RMSEA 0.000
 90 Percent Confidence Interval 0.000 0.000
 P-value RMSEA <= 0.05 NA

Standardized Root Mean Square Residual:

SRMR 0.000

Parameter Estimates:

Information Expected
 Standard Errors Standard

Regressions:
             Estimate  Std.Err   z-value   P(>|z|)
 vanzari ~
 inteligenta 177.199   183.916     0.963    0.335
 memorie      50.896     8.387     6.069    0.000
 cunostinte  141.435    29.317     4.824    0.000

Covariances:
               Estimate   Std.Err   z-value   P(>|z|)
 inteligenta ~~
 memorie         -5.301     8.503   -0.623      0.533
 cunostinte      39.517     6.041    6.542      0.000
 memorie ~~
 cunostinte      63.955    53.887    1.187      0.235

Variances:
                 Estimate       Std.Err     z-value    P(>|z|)
 .vanzari     6369159.463    891860.715       7.141     0.000
 inteligenta        7.371         1.032       7.141     0.000
 memorie          996.604       139.552       7.141     0.000
 cunostinte       293.092        41.041       7.141     0.000

Avem o mulțime de informații, pe multe dintre ele încă nu știm să le interpretăm, însă vom recunoaște cu ușurință coeficienții de regresie, covarianțele și varianțele unice. Lipsește însă un indicator foarte important al regresiei, anume pătratul coeficientului multiplu de corelație cu ajutorul căruia aflăm puterea predictivă a modelului. Din fericire, știm câteva trucuri în „lavaan” și vom folosi funcția „inspect()”:

> inspect(analiza, 'r2')
 vanzari
 0.643

Am indicat să se inspecteze obiectul rezultat în urma analizei și să se caute indicatorul dorit. Iată că R l-a găsit și ni l-a prezentat. Nu rămâne decât să afișăm diagrama de cale, inclusiv coeficienții nestandardizați, folosind deja cunoscuta funcție „semPaths”:

> semPaths(analiza, whatLabels = 'est',
          nodeLabels = c("Performanta", "Inteligenta",
                         "Memorie", "Cunostinte"),
          sizeMan = 15, sizeMan2 = 5,
          curve = 2, curvature = 0.7,
          color = list(man='green', lat='red'),
          edge.color = "brown",
          normalize = FALSE)

Am obținut, desigur, aceleași rezultate ca și în cazul utilizării Amos, orice comentarii fiind acum inutile, însă nu am terminat lucrul în R. Vom studia aceleași proceduri, folosind pachetul „sem”.

Deoarece pachetele „sem” și „lavaan” implementează funcții cu același nume, ar fi util ca înaintea încărcării pachetului „sem” să eliminăm „lavaan” din memorie:

> detach(package:semPlot, unload=TRUE)
> detach(package:lavaan, unload=TRUE)

Am detașat întâi pachetul „semPlot”, apoi „lavaan” deoarece primul folosește funcții din „lavaan” și ar fi blocat scoaterea sa din memorie. După aceea, desigur, încărcăm pachetul „sem”:

> library(sem)

Ne amintim că specificarea modelului se realizează cu ajutorul funcției „specifyModel()” sau „specifyEquation()”, folosindu-se standardul RAM, după regulile deja cunoscute:

> model.sem <- specifyModel(text='
                            # Specificarea coeficientilor de regresie
                            inteligenta -> vanzari, int_vanz, NA
                            memorie -> vanzari, mem_vanz, NA
                            cunostinte -> vanzari, cun_vanz, NA
                            # Specificarea covariantelor
                            inteligenta <-> memorie, int_mem, NA
                            inteligenta <-> cunostinte, int_cun, NA
                            memorie <-> cunostinte, mem_cun, NA
                            # Specificarea variantelor unice
                            inteligenta <-> inteligenta, int_int, NA
                            memorie <-> memorie, mem_mem, NA
                            cunostinte <-> cunostinte, cun_cun, NA
                            vanzari <-> vanzari, van_van, NA
                            ')

Deoarece nu mai avem opțiunea de identificare și estimare automată a varianțelor, trebuie specificați toți parametrii, apoi se va introduce modelul în memorie și va putea fi analizat.

> analiza <- sem(model = model.sem, data = regresii)
> summary(analiza)

Model Chisquare = -3.588241e-13 Df = 0 Pr(>Chisq) = NA
 AIC = 20
 BIC = -3.588241e-13

Normalized Residuals
    Min.      1st Qu.  Median     Mean    3rd Qu.     Max. 
-1.159e-15 0.000e+00 2.483e-16 7.153e-16 2.050e-15 2.937e-15

R-square for Endogenous Variables
vanzari 
 0.6432 
 Parameter Estimates
              Estimate       Std Error       z value       Pr(>|z|) 
int_vanz    1.771990e+02    1.848246e+02    0.9587414   3.376891e-01 vanzari <--- inteligenta 
mem_vanz    5.089570e+01    8.428155e+00    6.0387713   1.552921e-09 vanzari <--- memorie 
cun_vanz    1.414354e+02    2.946206e+01    4.8005927   1.581967e-06 vanzari <--- cunostinte 
int_mem    -5.353965e+00    8.629473e+00   -0.6204279   5.349761e-01 memorie <--> inteligenta 
int_cun     3.990856e+01    6.130750e+00    6.5095723   7.536506e-11 cunostinte <--> inteligenta 
mem_cun     6.458812e+01    5.468917e+01    1.1810038   2.376012e-01 cunostinte <--> memorie 
int_int     7.444408e+00    1.047573e+00    7.1063352   1.191648e-12 inteligenta <--> inteligenta
mem_mem     1.006471e+03    1.416301e+02    7.1063352   1.191648e-12 memorie <--> memorie 
cun_cun     2.959943e+02    4.165218e+01    7.1063352   1.191648e-12 cunostinte <--> cunostinte 
van_van     6.432220e+06    9.051389e+05    7.1063352   1.191648e-12 vanzari <--> vanzari

Iterations = 0

Nu intrăm în detalii, am obținut aceleași rezultate, exprimate în notație științifică, inclusiv pătratul coeficientului de corelație multiplă.

În cazul în care considerăm că specificarea modelului în notație RAM este prea incomodă, putem folosi, desigur, cea de-a doua variantă, obținându-se aceleași rezultate, care ar putea fi apoi transpuse sub formă de diagramă de cale, cu aceeași funcție, „semPaths”:

> model.sem <- specifyEquations(text='
                                # Specificarea ecuatiei de regresie
                                vanzari = b1*inteligenta + b2*memorie + b3*cunostinte
                                # Specificarea covariantelor
                               C(inteligenta, memorie)=int_mem
                               C(inteligenta, cunostinte)=int_cun
                               C(memorie, cunostinte)=mem_cun
                               # Specificarea variantelor unice
                               V(inteligenta)=int
                               V(memorie)=mem
                               V(cunostinte)=cun
                               V(vanzari)=van
                               ')
 > analiza <- sem(model = model.sem, data = regresii)
 > summary(analiza)

 Dacă această sintaxă vi se pare nefamiliară și dificilă, vă asigur că nu este așa. Deocamdată vă lipsesc unele noțiuni introductive. Răbdare, toate se vor lumina când va apărea cursul eLearnig sau cartea. Data viitoare vom studia cum se poate face același lucru, însă folosind Mplus.

Lasă un răspuns

Adresa ta de email nu va fi publicată. Câmpurile obligatorii sunt marcate cu *