Näide 1. Kinnitava faktoranalüüsi mudeli rakendamine

Lavaani mooduliga tutvumiseks teeme alustuseks suhteliselt lihtsa kinnitava faktoranalüüsi mudeli, milles üheksa vaadeldud muutuja abil defineerime kolm latentset muutujat (faktorit). Sageli on mõistlik teha enne andmete analüüsimist skeem muutujate vaheliste seoste kohta. Praegusel juhul näeks mudeli diagramm välja järgmine:

Joonis 1.

Jooniselt on näha, et vaadeldud muutujad x1-x3 mõõdavad latentset muutujat nimega visual, x4-x6 latentset muutujat nimega textual ja x7-x9 latentset muutujat nimega speed. Vastavalt tavale tähistavad SEM-i diagrammidel nelinurgad vaadeldud muutujaid ja ringid või ellipsid latentseid muutujaid. Latentsete muutujate vahelised kaared tähistavad nendevahelisi korrelatsioone ja väikesed noolekesed muutujate juures tähistavad dispersioone. Vaadeldud muutujate puhul on tegemist jääkdispersiooniga ehk selle osaga muutuja variatiivsusest, mis ei ole määratud latentse faktori poolt. Latentsetele muutujatele on vaja anda ka skaala, selleks fikseeritakse iga faktori puhul ühe sellele laaduva vaadeldud muutuja laadung. Seda tähistavad number 1-d mõne noole juures.

Mudeli diagramm aitab ka hinnata, kas analüüsi kaasatavas andmestikus on piisavalt infot, et arvutada välja kõik need mudeli parameetrid, mida meil vaja läheb. Jooniselt saame kokku lugeda, et meil on vaja hinnata 6 faktorlaadungit (9 vaadeldavat muutujat miinus need 3, mille laadungid on fikseeritud), 3 korrelatsiooni ja 12 dispersiooni. See teeb kokku 21 parameetrit. Maksimaalne parameetrite arv sõltub mudelisse kaasatavate vaadeldud muutujate arvust ja seda saab arvutada valemi \(p(p+1)/2\) abil, kus “p” on mudelisse kaasatavate vaadeldud muutujate arv. Praegusel juhul anname mudelile ette 9 muutujat ja sellest tulenevalt maksimaalne vabade parameetrite arv mudelis oleks \(9(9+1)/2=45\). Järelikult pole antud juhul meil probleemi sellega, et tahaksime mudelil lasta hinnata liiga suurt arvu vabu parameetreid.

SEM eeldab ka, et analüüsi kaasatavad muutujad oleksid normaaljaotusega ja et nad ei oleks liiga tugevalt seotud. Lisaks sellele võib SEM-i puhul probleeme tekitada see, kui analüüsi kaasatavad muutujad on väga erinevatel numbrilistel skaaladel (nt mõned muutujad ühekohalised komakohtadega arvud ja teised kolmekohalised arvud). Erinevalt uurivast faktoranalüüsist, mis teostab arvutusi standardiseeritud korrelatsioonikordajate peal, kasutab SEM kovariatsioonikordajaid, mis sisaldavad infot ka muutujate skaala kohta. Väga erinevate suurusjärkude korral tasub korrutada või jagada mõned muutujad mingi koefitsiendiga, et kõik muutujad oleksid enam-vähem sarnastes suurusjärkudes. Ühe- ja kolmekohalistest arvudest koosnevate muutujate puhul saab ühekohalised muutujad korrutada läbi 100-ga.

Mudeli kirjeldamine

Järgnevalt peaksime mudeli viima R-ile ja lavaanile arusaadavale kujule. Mudeli kirjeldamisel saame märkida erinevaid seoseid:

1. regressioonseose tähistus:
\(sõltuv muutuja \sim sõltumatu.muutuja\)

2. latentse muutuja defineerimine (kus x1, x2 ja x3 on latentsele tunnusele laaduvad vaadeldud muutujad):

\(latentne.muutuja =\sim x1 + x2 + x3\)

3. muutujate-vahelise korrelatsiooni saame mudelisse kahekordse tilde abil:

\(muutuja1 \sim \sim muutuja2\)

Neid kolme tüüpi seoseid läheb mudelite kirjeldamisel kõige sagedamini vaja. Antud juhul näeks meie 3-faktoriline mudel välja järgmine:

mudel1 <- "
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9"

Mudeli kirjeldus peab paiknema jutumärkide (või ülakomade) vahel ja iga mudeli valem peab olema eraldi real. Pikemate mudelite korral on vahel selguse huvides kasulik mudelit liigendada, jättes valemite vahele tühje ridu või kommentaare. Nagu mujalgi R-i koodis saab sisestada kommentaare sümboli “#”" abil.

Mudeli väljund

Mudeli saab sobitada oma andmetele kasutades funktsiooni cfa (confirmatory factor analysis), millele anname argumentideks mudeli kirjelduse ja andmestiku nime.

fit1 <- cfa(mudel1, data=HolzingerSwineford1939)

Seejärel saame funktsiooni summary abil uurida mudeli väljundit. Lavaani mudelite puhul võib lisade veel kaks argumenti: fit.measures=TRUE (annab lisaks hii-ruudule täiendavaid mudeli sobitusastme näitajaid), standardized=TRUE (annab lisaks standardiseerimata mudeli koefitsientidele ka standardiseeritud koefitsiendid).

summary(fit1, fit.measures=TRUE, standardized=TRUE)

## lavaan (0.5-23.1097) converged normally after  35 iterations
## 
##   Number of observations                           301
## 
##   Estimator                                         ML
##   Minimum Function Test Statistic               85.306
##   Degrees of freedom                                24
##   P-value (Chi-square)                           0.000
## 
## Model test baseline model:
## 
##   Minimum Function Test Statistic              918.852
##   Degrees of freedom                                36
##   P-value                                        0.000
## 
## User model versus baseline model:
## 
##   Comparative Fit Index (CFI)                    0.931
##   Tucker-Lewis Index (TLI)                       0.896
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -3737.745
##   Loglikelihood unrestricted model (H1)      -3695.092
## 
##   Number of free parameters                         21
##   Akaike (AIC)                                7517.490
##   Bayesian (BIC)                              7595.339
##   Sample-size adjusted Bayesian (BIC)         7528.739
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.092
##   90 Percent Confidence Interval          0.071  0.114
##   P-value RMSEA <= 0.05                          0.001
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.065
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Standard Errors                             Standard
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   visual =~                                                             
##     x1                1.000                               0.900    0.772
##     x2                0.554    0.100    5.554    0.000    0.498    0.424
##     x3                0.729    0.109    6.685    0.000    0.656    0.581
##   textual =~                                                            
##     x4                1.000                               0.990    0.852
##     x5                1.113    0.065   17.014    0.000    1.102    0.855
##     x6                0.926    0.055   16.703    0.000    0.917    0.838
##   speed =~                                                              
##     x7                1.000                               0.619    0.570
##     x8                1.180    0.165    7.152    0.000    0.731    0.723
##     x9                1.082    0.151    7.155    0.000    0.670    0.665
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   visual ~~                                                             
##     textual           0.408    0.074    5.552    0.000    0.459    0.459
##     speed             0.262    0.056    4.660    0.000    0.471    0.471
##   textual ~~                                                            
##     speed             0.173    0.049    3.518    0.000    0.283    0.283
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .x1                0.549    0.114    4.833    0.000    0.549    0.404
##    .x2                1.134    0.102   11.146    0.000    1.134    0.821
##    .x3                0.844    0.091    9.317    0.000    0.844    0.662
##    .x4                0.371    0.048    7.779    0.000    0.371    0.275
##    .x5                0.446    0.058    7.642    0.000    0.446    0.269
##    .x6                0.356    0.043    8.277    0.000    0.356    0.298
##    .x7                0.799    0.081    9.823    0.000    0.799    0.676
##    .x8                0.488    0.074    6.573    0.000    0.488    0.477
##    .x9                0.566    0.071    8.003    0.000    0.566    0.558
##     visual            0.809    0.145    5.564    0.000    1.000    1.000
##     textual           0.979    0.112    8.737    0.000    1.000    1.000
##     speed             0.384    0.086    4.451    0.000    1.000    1.000

Mudeli väljundi ülaosas on näha valimi suurus ning hii-ruut-statistik (Minimum Function Test Statistic) koos vabadusastmete arvu ja p-väärtusega. Hii-ruut testib antud juhul nn täpse sobitumise hüpoteesi, mille kohaselt ei erine mudeli põhjal arvutatud kovaritasioonimaatriks andmestiku põhjal arvutatud kovariatsioonimaatriksist. Me tahame, et hii-ruut oleks võimalikult väike ja p-väärtus võimalikult suur. Kui me tüüpiliselt tahame, et p-väärtus oleks alla 0.05-e, siis praegu tahame, et see oleks 0.05-st suurem. Praegusel juhul näitab hii-ruut, et mudel andmetele kuigi hästi ei sobi.

Järgnevalt on ära toodud täiendavad sobitusastme näitajad, mida me funktsioonilt summary tellisime argumendi fit.measures abil. Sobitusastme näitajad võib laias laastus jagada kahte rühma: lisanduva sobitusastme indeksid (incrementalt indices) ja absoluutse sobitusastme indeksid (absolute fit indices). Esimesed mõõdavad sobitusastme paranemist võrreldes sellise mudeliga, milles vaadeldud muutujate vahelised korrelatsioonid puuduvad. Teised mõõdavad, kui hästi suudab meie poolt koostatud mudel taastada tegelike andmete põhjal arvutatud kovariatsioonimaatriksi. Lisanduva sobitusastme indeksistest on väljundis toodud TLI (Tucker-Lewis index) ja CFI (comparative fit index). Nende puhul näitavad mudelis head sobitusastet väärtused 0.95 ja üles. Absoluutse sobitusastme indeksitest on ära toodud RMSEA (root mean square error of approximation) ja SRMR (standardized root mean square residual), mille puhul head sobitusastet näitavad väärtused 0.05 ja alla. RMSEA puhul on ära toodud ka 90% usalduspiirid. Hea mudeli korral oleks alumine piir alla 0.05 ja ülemine mitte üle 0.10. Ka kõigi nende näitajate põhjal võime öelda, et sellisel kujul antud mudel andmetega ei sobi.

Sobitusastme indeksite järel tuleb tabel mudeli parameetrite hinnangutega. Tabelis on kõigepealt toodud fakorlaadungid (alajaotus Latent variables), seejärel faktorite vahelised seosed (alajaotus Covariance) ja dispersioonid (Variances). Standardiseerimata parameetrid on tulbas nimega Estimates ja standardiseeritud parameetrid tulbas nimega Std.all. Lisaks sellele on tabelis toodud ka standardvead (Std.err), z-väärtused (mis on standardiseerimata parameeter jagatuna st.veaga) ning p-väärtused. (Tulbas Std.lv on parameetrite väärtused, siis kui standardiseerida ainult latentsed muutujad).

Nagu näha on muutujate x1, x4 ja x7 standardiseerimata laadung 1.000 ja mõned näitajad on neist ridadest puudu. See tuleneb sellest, et struktuurivõrrandite mudelite puhul tuleb latentsele muutujale anda skaala. Kõige sagedamini antakse mõne mõõdetud muutuja skaala fikseerides selle standardiseerimata laadung ühega. Lavaan, fikseerib automaatselt iga latentse muutuja puhul esimese indikaatori standardiseerimata laadungi.

Kuidas parameetreid tõlgendada? Standardiseerimata laadungid on tõlgendatavad regressioonikoefitsientidena. Kui muutuja x2 standardiseerimata laadung faktorile visual on 0.554 siis võib faktori väärtuse 1-punktilise suurenemise korral oodata, et x2 suureneb 0.554 punkti võrra.

Ainult ühele faktorile laaduvate vaadeldud muutujate puhul võib standardiseeritud laadungeid tõlgendada faktori ja muutuja vahelise korrelatsioonikordajana. Nagu korrelatsioonikordaja puhul ikka saame seda ruutu võttes kindlaks teha kui suure osa tunnuse variatiivsusest seos ära määrab. Tunnuse x2 puhul on korrelatsioon faktoriga visual 0.424 ja faktor määrab ära 0.42422=0.179 ehk umbes 18 protsenti tunnuse x2 variatiivsusest. Muutujate puhul, mis laaduvad rohkem kui ühele faktorile (praeguses mudelis meil selliseid pole) ei saa standardiseeritud laadungeid tõlgendada korrelatsioonikordajatena vaid standardiseeritud regressioonikordajatena, mille puhul hoitakse korreleeritud faktoreid kontrolli all. Neid ruutu võttes ei saa kindlaks teha seletatava variatiivsuse osakaalu.

Alajaotuses Covariances olevad standardiseerimata parameetrid kujutavad endas faktorite vahelisi kovariatsioonikordajaid ja standardiseeritud parameetrid korrelatsioonikordajaid. Antud mudeli puhul võib näiteks näha et faktori visual korrelatsioon faktoriga textual on 0.459 ja faktoriga speed 0.471. Vaadeldud muutujate puhul näitavad dispersioonid (Variances) standardiseeritud kujul seda osa tunnuse variatiivsusest, mida faktor ei seleta. Nt tunnuse x2 puhul ei suuda faktor visual seletada tervelt 82.1 protsenti variatiivsusest.

Parameetrite tabeli puhul tasub ka kontrollida, et ei esineks negatiivseid jääkdispersioone ja standardiseeritud parameetreid, mille absoluutväärtus on suurem kui üks.. Selliseid väärtusi ei ole võimalik tõlgendada ja nende esinemine annab märku probleemidest mudeli või andmetega (nt liiga väike valim või liiga tugevalt korreleeritud muutujad). Mudeli väljundit saab tellida ka osade kaupa. Seda võib vaja minna näiteks siis kui soovime teada saada mõnda sellist sobitusastme näitajat, mida funktsioon summary välja ei trükkinud. Nt üks suhteliselt sageli raporteeritav sobitusastme näitaja on GFI (goodness of fit index, hea väärtus 0.95 ja üles). Selle ja suure hulga muid sobitusindekseid saame kätte funktsiooni fitMeasures abil:

fitMeasures(fit1)

Aga eraldi saab tellida ka standardiseerimata parameetrite hinnanguid…

parameterEstimates(fit1)

… ja standardiseeritud hinnanguid.

standardizedSolution(fit1)

Kui mudel andmetega ei sobi, siis üks võimalus probleeme diagnoosida, on vaadata mudeli jääkide maatriksit. Jäägid kujutavad endast erinevusi tegelike andmete põhjal arvutatud ja mudeli põhjal taastatud kovariatsioonikordajate vahel. Jääke, mida struktuurivõrrandite mudelite puhul uurida saab, on erinevat tüüpi. Vaatame selliseid, mida nimetatakse korrelatsioonijääkideks (ehk siis kovariatsioonikordajad arvutatakse ümber korrelatsioonideks, mida on lihtsam tõlgendada). Korrelatsioonijäägid saame funktsiooni residuals abil, määrates argumendi type väärtuseks “cor”. Tabelist peaksime üles otsima jääkkorrelatsioonid absoluutväärtusega > 0.10. Need nn suured jäägid näitavad meile, milliste muutujate vaheliste seoste seletamisega mudel väga hästi hakkama ei saa.

residuals(fit1, type="cor")

## $type
## [1] "cor.bollen"
## 
## $cor
##    x1     x2     x3     x4     x5     x6     x7     x8     x9    
## x1  0.000                                                        
## x2 -0.030  0.000                                                 
## x3 -0.008  0.094  0.000                                          
## x4  0.071 -0.012 -0.068  0.000                                   
## x5 -0.009 -0.027 -0.151  0.005  0.000                            
## x6  0.060  0.030 -0.026 -0.009  0.003  0.000                     
## x7 -0.140 -0.189 -0.084  0.037 -0.036 -0.014  0.000              
## x8 -0.039 -0.052 -0.012 -0.067 -0.036 -0.022  0.075  0.000       
## x9  0.149  0.073  0.147  0.048  0.067  0.056 -0.038 -0.032  0.000
## 
## $mean
## x1 x2 x3 x4 x5 x6 x7 x8 x9 
##  0  0  0  0  0  0  0  0  0

Antud juhul näeme, et muutuja x3 omab suurt jääkkorrelatsiooni muutujaga x5, muutuja x7 muutujatega x1 ja x2 ning muutuja x9 muutujatega x1 ja x3. Kui muutuja omab suuri jääkkorrelatsioone mõnele teisele faktorile laaduvate muutujatega, on võimalik, et muutuja mõõdab rohkem kui ühte konstrukti ja tal tuleks lubada laaduda rohkem kui ühele faktorile. Teine võimalus on, et muutujate vahelisi seoseid ei põhjusta mitte ainult mudeli faktor(id), vaid ka mingi muu tegur. Sellise võimaluse arvesse võtmiseks saab mudelis määrata jääkdispersioonide vahelised korrelatsioonid. Teeme proovi võttes uue mudeli, milles määrame ka muutujate vahelised korrelatsioonid, mida faktorid ei seleta.

mudel2 <- "
# latentsed muutujad
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9
# jääkdispersioonide vahelised korrelatsioonid
x7 ~~ x1 + x2
x9 ~~ x1 + x3
x3 ~~ x5"

Vaadake selle mudeli väljundit. Kas sobitusastme indeksid läksid paremaks või mitte?

fit2 <- cfa(mudel2, data=HolzingerSwineford1939)
summary(fit2, fit.measures=TRUE, standardized=TRUE)

## lavaan (0.5-23.1097) converged normally after  36 iterations
## 
##   Number of observations                           301
## 
##   Estimator                                         ML
##   Minimum Function Test Statistic               40.161
##   Degrees of freedom                                19
##   P-value (Chi-square)                           0.003
## 
## Model test baseline model:
## 
##   Minimum Function Test Statistic              918.852
##   Degrees of freedom                                36
##   P-value                                        0.000
## 
## User model versus baseline model:
## 
##   Comparative Fit Index (CFI)                    0.976
##   Tucker-Lewis Index (TLI)                       0.955
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -3715.173
##   Loglikelihood unrestricted model (H1)      -3695.092
## 
##   Number of free parameters                         26
##   Akaike (AIC)                                7482.346
##   Bayesian (BIC)                              7578.730
##   Sample-size adjusted Bayesian (BIC)         7496.273
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.061
##   90 Percent Confidence Interval          0.034  0.087
##   P-value RMSEA <= 0.05                          0.227
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.043
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Standard Errors                             Standard
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   visual =~                                                             
##     x1                1.000                               0.859    0.747
##     x2                0.562    0.108    5.217    0.000    0.483    0.412
##     x3                0.710    0.118    6.039    0.000    0.610    0.546
##   textual =~                                                            
##     x4                1.000                               0.986    0.849
##     x5                1.105    0.065   16.916    0.000    1.090    0.852
##     x6                0.931    0.056   16.763    0.000    0.918    0.840
##   speed =~                                                              
##     x7                1.000                               0.680    0.623
##     x8                1.094    0.150    7.295    0.000    0.744    0.736
##     x9                0.913    0.126    7.265    0.000    0.621    0.616
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .x1 ~~                                                                 
##    .x7               -0.148    0.058   -2.546    0.011   -0.148   -0.226
##  .x2 ~~                                                                 
##    .x7               -0.241    0.063   -3.800    0.000   -0.241   -0.265
##  .x1 ~~                                                                 
##    .x9                0.173    0.059    2.941    0.003    0.173    0.285
##  .x3 ~~                                                                 
##    .x9                0.189    0.055    3.455    0.001    0.189    0.254
##    .x5               -0.135    0.047   -2.911    0.004   -0.135   -0.216
##   visual ~~                                                             
##     textual           0.398    0.071    5.638    0.000    0.470    0.470
##     speed             0.257    0.059    4.376    0.000    0.440    0.440
##   textual ~~                                                            
##     speed             0.187    0.053    3.526    0.000    0.278    0.278
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .x1                0.584    0.117    4.980    0.000    0.584    0.442
##    .x2                1.139    0.102   11.115    0.000    1.139    0.830
##    .x3                0.874    0.091    9.622    0.000    0.874    0.702
##    .x4                0.378    0.048    7.938    0.000    0.378    0.280
##    .x5                0.449    0.058    7.715    0.000    0.449    0.274
##    .x6                0.353    0.043    8.273    0.000    0.353    0.295
##    .x7                0.728    0.083    8.740    0.000    0.728    0.612
##    .x8                0.469    0.075    6.219    0.000    0.469    0.459
##    .x9                0.631    0.069    9.162    0.000    0.631    0.621
##     visual            0.738    0.144    5.125    0.000    1.000    1.000
##     textual           0.973    0.112    8.707    0.000    1.000    1.000
##     speed             0.462    0.096    4.830    0.000    1.000    1.000

Näide 2: regressiooniseostega mudeli rakendamine

Proovime teha sellise mudeli, mis sisaldab lisaks latentsetele muutujatele ka regressioonseoseid. Selleks vaatame teist lavaaniga kaasa tulevat näidisandmestikku nimega “PoliticalDemocracy”. See andemestik sisaldab 75 riigi kohta 11 poliitilist ja majanduslikku arengut puudutavat näitajat.

Muutujad y1 kuni y4 on vastavalt ajakirjandusvabaduse tase, vaba politilise opositsiooni olemasolu, vabade valimiste toimumine ja valitud parlamendi efektiivsus mõõdetuna aastal 1960. Muutujad y5 kuni y8 on samad näitajad mõõdetuna aastal 1965. Muutujad x1 kuni x3 on vastavalt riigi sisemajanduse koguprodukt, energiatarbimine elaniku kohta ja tööstussektoris töötava tööjõu osakaal mõõdetuna aastal 1960. Andmestiku tutvustus avaneb alljärgneva käsuga:

help(PoliticalDemocracy)

Üritame andmetele sobitada alljärgneval joonisel kujutatud mudelit. Defineerime kolm latentset muutujat. Tunnused x1 kuni x3 mõõdavad latentset muutujat maj60 (st majanduslik arengutase aastal 1960), tunnused y1 kuni y4 mõõdavad latentset muutujat dem60 (st demokraatia aastal 1960) ning tunnused y5 kuni y8 mõõdavad latentset muutujat dem65 (st demokraatia aastal 1965). Latentsete muutujate vahel on regressioonseosed. Demokraatiat aastal 1960 mõjutab majanduse arengutase aastal 1960. Demokraatiat aastal 1965 mõjutavad nii majandus kui demokraatia aastal 1960. Latentsed muutujad ei suuda antud juhul seletada kõiki demokraatia indikaatorite vahelisi seoseid. Seetõttu on mudelisse lisatud ka nende jääkdispersioonide vahelised korrelatsioonid. Muutuja y1 korreleerub y5-ga, y2 korreleerub y4-ga ja y6-ga, y3 korreleerub y7-ga, y8 korreleerub y4-ga ja y6-ga. Mudelisse lisatud 11 vaadeldud muutujat võimaldavad hinnata \(11(11+1)/2 = 66\) vaba parameetrit. Meil on vaja hinnata 31 parameetrit: 8 faktorlaadungit, 14 dispersiooni, 3 regressiooni- ja 6 korrelatsioonikordajat. (Nagu esimese näite puhul fikseeritakse ka praegu iga latentse muutuja puhul esimese indikaatori standardiseerimata laadung, selleks et omistada latentsele muutujale skaala. Seetõttu ei olegi meil vaja hinnata 11 vaid ainult 8 faktorlaadungit). Lavaanis näeks selline mudel välja järgnevalt:

Joonis 2.

mudel3 <- "
# defineerime latentsed tunnused
maj60 =~ x1 + x2 + x3
dem60 =~ y1 + y2 + y3 + y4
dem65 =~ y5 + y6 + y7 + y8

# regressioonseosed latentsete tunnuste vahel
dem60 ~ maj60
dem65 ~ maj60 + dem60
# jääkdispersioonide vahelised korrelatsioonid
y1 ~~ y5
y2 ~~ y4 + y6
y3 ~~ y7
y4 ~~ y8
y6 ~~ y8"

Sobitame mudeli andmetele ja uurime väljundit. Kuna mudel sisaldab ka regressioonseoseid, siis kasutame sobitamiseks funktsiooni sem. Kas mudel sobib andmetega?

fit3 <- sem(mudel3, data=PoliticalDemocracy)

## Warning in lav_data_full(data = data, group = group, cluster = cluster, :
## lavaan WARNING: unordered factor(s) with more than 2 levels detected in
## data: x1 x2 x3 y1 y2 y3 y4 y5 y6 y7 y8

summary(fit3, fit.measures=TRUE, standardized=TRUE)

## lavaan (0.5-23.1097) converged normally after 223 iterations
## 
##   Number of observations                            75
## 
##   Estimator                                         ML
##   Minimum Function Test Statistic               34.580
##   Degrees of freedom                                35
##   P-value (Chi-square)                           0.488
## 
## Model test baseline model:
## 
##   Minimum Function Test Statistic              542.972
##   Degrees of freedom                                55
##   P-value                                        0.000
## 
## User model versus baseline model:
## 
##   Comparative Fit Index (CFI)                    1.000
##   Tucker-Lewis Index (TLI)                       1.001
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)                     NA
##   Loglikelihood unrestricted model (H1)             NA
## 
##   Number of free parameters                         31
##   Akaike (AIC)                                      NA
##   Bayesian (BIC)                                    NA
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.000
##   90 Percent Confidence Interval          0.000  0.082
##   P-value RMSEA <= 0.05                          0.751
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.049
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Standard Errors                             Standard
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   maj60 =~                                                              
##     x1                1.000                              18.271    0.912
##     x2                1.079    0.068   15.983    0.000   19.720    0.987
##     x3                0.801    0.067   11.923    0.000   14.632    0.873
##   dem60 =~                                                              
##     y1                1.000                               5.168    0.635
##     y2                0.770    0.164    4.689    0.000    3.978    0.700
##     y3                0.772    0.167    4.617    0.000    3.988    0.720
##     y4                0.613    0.151    4.066    0.000    3.170    0.587
##   dem65 =~                                                              
##     y5                1.000                               3.042    0.725
##     y6                1.115    0.207    5.392    0.000    3.391    0.675
##     y7                1.750    0.277    6.313    0.000    5.323    0.840
##     y8                1.133    0.249    4.548    0.000    3.446    0.571
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   dem60 ~                                                               
##     maj60             0.090    0.039    2.316    0.021    0.320    0.320
##   dem65 ~                                                               
##     maj60             0.039    0.015    2.597    0.009    0.231    0.231
##     dem60             0.490    0.092    5.305    0.000    0.832    0.832
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .y1 ~~                                                                 
##    .y5                6.348    2.762    2.298    0.022    6.348    0.349
##  .y2 ~~                                                                 
##    .y4                2.687    2.217    1.212    0.226    2.687    0.151
##    .y6                7.515    2.297    3.272    0.001    7.515    0.499
##  .y3 ~~                                                                 
##    .y7                1.664    2.589    0.643    0.520    1.664    0.126
##  .y4 ~~                                                                 
##    .y8                6.355    2.799    2.270    0.023    6.355    0.293
##  .y6 ~~                                                                 
##    .y8                1.882    2.011    0.936    0.349    1.882    0.102
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .x1               67.509   15.010    4.498    0.000   67.509    0.168
##    .x2               10.704   11.920    0.898    0.369   10.704    0.027
##    .x3               66.654   12.724    5.238    0.000   66.654    0.237
##    .y1               39.621    7.758    5.107    0.000   39.621    0.597
##    .y2               16.489    3.517    4.689    0.000   16.489    0.510
##    .y3               14.787    3.484    4.244    0.000   14.787    0.482
##    .y4               19.146    3.631    5.272    0.000   19.146    0.656
##    .y5                8.330    1.708    4.876    0.000    8.330    0.474
##    .y6               13.736    2.631    5.221    0.000   13.736    0.544
##    .y7               11.801    3.499    3.373    0.001   11.801    0.294
##    .y8               24.606    4.399    5.594    0.000   24.606    0.674
##     maj60           333.828   65.403    5.104    0.000    1.000    1.000
##    .dem60            23.977    8.735    2.745    0.006    0.898    0.898
##    .dem65             1.216    0.795    1.530    0.126    0.131    0.131

Modifikatsiooniindeksid

Modifikatsiooniindeksid näitavad, kui palju väheneks mudeli hii-ruut statistik mingi vaba parameetri mudelisse lisamisel. (Mida väiksem hiiruut, seda paremini mudel andmetega sobib.) Hoiatusena tuleks märkida, et indeksite puhul on tegemist umbkaudsete oletustega. Statistikaprogramm ei proovi kõigi nende parameetrite mudelisse lisamist ka tegelikult läbi, vaid kasutab nende hindamiseks mingit suhteliselt jämedat rusikareeglit.

Defineerime sama 3-faktorilise mudeli:

HS.model <- "visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9"

Sobitame mudeli andmetele ja vaatame mudeli väljundit.

fit1 <- cfa(HS.model, data=HolzingerSwineford1939)
summary(fit1, fit.measures=TRUE)

## lavaan (0.5-23.1097) converged normally after  35 iterations
## 
##   Number of observations                           301
## 
##   Estimator                                         ML
##   Minimum Function Test Statistic               85.306
##   Degrees of freedom                                24
##   P-value (Chi-square)                           0.000
## 
## Model test baseline model:
## 
##   Minimum Function Test Statistic              918.852
##   Degrees of freedom                                36
##   P-value                                        0.000
## 
## User model versus baseline model:
## 
##   Comparative Fit Index (CFI)                    0.931
##   Tucker-Lewis Index (TLI)                       0.896
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -3737.745
##   Loglikelihood unrestricted model (H1)      -3695.092
## 
##   Number of free parameters                         21
##   Akaike (AIC)                                7517.490
##   Bayesian (BIC)                              7595.339
##   Sample-size adjusted Bayesian (BIC)         7528.739
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.092
##   90 Percent Confidence Interval          0.071  0.114
##   P-value RMSEA <= 0.05                          0.001
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.065
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Standard Errors                             Standard
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   visual =~                                           
##     x1                1.000                           
##     x2                0.554    0.100    5.554    0.000
##     x3                0.729    0.109    6.685    0.000
##   textual =~                                          
##     x4                1.000                           
##     x5                1.113    0.065   17.014    0.000
##     x6                0.926    0.055   16.703    0.000
##   speed =~                                            
##     x7                1.000                           
##     x8                1.180    0.165    7.152    0.000
##     x9                1.082    0.151    7.155    0.000
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   visual ~~                                           
##     textual           0.408    0.074    5.552    0.000
##     speed             0.262    0.056    4.660    0.000
##   textual ~~                                          
##     speed             0.173    0.049    3.518    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .x1                0.549    0.114    4.833    0.000
##    .x2                1.134    0.102   11.146    0.000
##    .x3                0.844    0.091    9.317    0.000
##    .x4                0.371    0.048    7.779    0.000
##    .x5                0.446    0.058    7.642    0.000
##    .x6                0.356    0.043    8.277    0.000
##    .x7                0.799    0.081    9.823    0.000
##    .x8                0.488    0.074    6.573    0.000
##    .x9                0.566    0.071    8.003    0.000
##     visual            0.809    0.145    5.564    0.000
##     textual           0.979    0.112    8.737    0.000
##     speed             0.384    0.086    4.451    0.000

Nagu näha, ei sobi mudel andmetega väga hästi. Hii-ruudu p-väärtus on alla 0.05-e, TLI ja CFI on alla 0.95-e ning RMSEA ja SRMR on üle 0.05-e. Järgnevalt vaatame mudeli modifikatsiooniindeksite tabelit ja otsime sellest üles kõige suurema indeksiga parameetri.

modindices(fit1)

##        lhs op rhs     mi    epc sepc.lv sepc.all sepc.nox
## 25  visual =~  x4  1.211  0.077   0.069    0.059    0.059
## 26  visual =~  x5  7.441 -0.210  -0.189   -0.147   -0.147
## 27  visual =~  x6  2.843  0.111   0.100    0.092    0.092
## 28  visual =~  x7 18.631 -0.422  -0.380   -0.349   -0.349
## 29  visual =~  x8  4.295 -0.210  -0.189   -0.187   -0.187
## 30  visual =~  x9 36.411  0.577   0.519    0.515    0.515
## 31 textual =~  x1  8.903  0.350   0.347    0.297    0.297
## 32 textual =~  x2  0.017 -0.011  -0.011   -0.010   -0.010
## 33 textual =~  x3  9.151 -0.272  -0.269   -0.238   -0.238
## 34 textual =~  x7  0.098 -0.021  -0.021   -0.019   -0.019
## 35 textual =~  x8  3.359 -0.121  -0.120   -0.118   -0.118
## 36 textual =~  x9  4.796  0.138   0.137    0.136    0.136
## 37   speed =~  x1  0.014  0.024   0.015    0.013    0.013
## 38   speed =~  x2  1.580 -0.198  -0.123   -0.105   -0.105
## 39   speed =~  x3  0.716  0.136   0.084    0.075    0.075
## 40   speed =~  x4  0.003 -0.005  -0.003   -0.003   -0.003
## 41   speed =~  x5  0.201 -0.044  -0.027   -0.021   -0.021
## 42   speed =~  x6  0.273  0.044   0.027    0.025    0.025
## 43      x1 ~~  x2  3.606 -0.184  -0.184   -0.134   -0.134
## 44      x1 ~~  x3  0.935 -0.139  -0.139   -0.105   -0.105
## 45      x1 ~~  x4  3.554  0.078   0.078    0.058    0.058
## 46      x1 ~~  x5  0.522 -0.033  -0.033   -0.022   -0.022
## 47      x1 ~~  x6  0.048  0.009   0.009    0.007    0.007
## 48      x1 ~~  x7  5.420 -0.129  -0.129   -0.102   -0.102
## 49      x1 ~~  x8  0.634 -0.041  -0.041   -0.035   -0.035
## 50      x1 ~~  x9  7.335  0.138   0.138    0.117    0.117
## 51      x2 ~~  x3  8.532  0.218   0.218    0.164    0.164
## 52      x2 ~~  x4  0.534 -0.034  -0.034   -0.025   -0.025
## 53      x2 ~~  x5  0.023 -0.008  -0.008   -0.005   -0.005
## 54      x2 ~~  x6  0.785  0.039   0.039    0.031    0.031
## 55      x2 ~~  x7  8.918 -0.183  -0.183   -0.143   -0.143
## 56      x2 ~~  x8  0.054 -0.012  -0.012   -0.010   -0.010
## 57      x2 ~~  x9  1.895  0.075   0.075    0.063    0.063
## 58      x3 ~~  x4  0.142 -0.016  -0.016   -0.012   -0.012
## 59      x3 ~~  x5  7.858 -0.130  -0.130   -0.089   -0.089
## 60      x3 ~~  x6  1.855  0.055   0.055    0.044    0.044
## 61      x3 ~~  x7  0.638 -0.044  -0.044   -0.036   -0.036
## 62      x3 ~~  x8  0.059 -0.012  -0.012   -0.011   -0.011
## 63      x3 ~~  x9  4.126  0.102   0.102    0.089    0.089
## 64      x4 ~~  x5  2.534  0.186   0.186    0.124    0.124
## 65      x4 ~~  x6  6.220 -0.235  -0.235   -0.185   -0.185
## 66      x4 ~~  x7  5.920  0.098   0.098    0.078    0.078
## 67      x4 ~~  x8  3.805 -0.069  -0.069   -0.059   -0.059
## 68      x4 ~~  x9  0.196 -0.016  -0.016   -0.014   -0.014
## 69      x5 ~~  x6  0.916  0.101   0.101    0.072    0.072
## 70      x5 ~~  x7  1.233 -0.049  -0.049   -0.035   -0.035
## 71      x5 ~~  x8  0.347  0.023   0.023    0.018    0.018
## 72      x5 ~~  x9  0.999  0.040   0.040    0.031    0.031
## 73      x6 ~~  x7  0.259 -0.020  -0.020   -0.017   -0.017
## 74      x6 ~~  x8  0.275  0.018   0.018    0.016    0.016
## 75      x6 ~~  x9  0.097 -0.011  -0.011   -0.010   -0.010
## 76      x7 ~~  x8 34.145  0.536   0.536    0.488    0.488
## 77      x7 ~~  x9  5.183 -0.187  -0.187   -0.170   -0.170
## 78      x8 ~~  x9 14.946 -0.423  -0.423   -0.415   -0.415

Modifikatsiooniindeksite tabeli 3 esimest veergu näitavad, millise parameetri kohta indeks käib. Tulbas lhs (sõnadest left hand side) on vasakpoolne muutuja ja tulbas rhs (right hand side) on sama parempoolne muutuja. Tulp op (operator) näitab meile, millist tüüpi seosega on tegu: =\(\sim\) on latentsele muutujale laadumine ning \(\sim\) \(\sim\) puhul on tegemist muutujate vahelise korrelatsiooniga.
Parameetri modifikatsiooniindeksi väärtus on toodud tulbas mi. Antud juhul puudutavad tabeli 24 esimest rida seoseid, mis on mudelis juba vabade parameetritena olemas ja pole seetõttu väga huvitavad. 25. real olev modifikatsiooniindeksi väärtus 1.211 näitab seda, et kui me lubaksime muutujal x4 laaduda ka faktorile visual, võiksime oodata, et mudeli hii-ruut väheneb umbes 1.2 võrra.
Lisaks modifikatsiooniindeksile arvutatakse välja ka parameetri oodatav väärtus (tulp epc, expected parameter change) ja selle standardiseeritud väärtus (meid huvitab eelkõige tulp sepc.all). Mudelisse tasub vabade parameetritena lisade eelkõige neid, mille puhul on suured nii mi kui epc väärtused. Modifikatsiooniindeksite puhul loetakse statistiliselt oluliseks väärtusi üle 3.84. Väiksema indeksi väärtusega parameetrite mudelisse lisamine üldiselt ära ei tasu. Indeksid ei ole üksteisest sõltumatud. See tähendab, et kui tabelis on mitu olulist indeksit, siis ei saa me neid korraga mudelisse lisada ja oodata, et hii-ruut kahaneb nende indeksite summa võrra. Kui me seda teeksime, võivad erinevused oodatavate ja tegelike muutuste vahel olla väga suured. Seega tuleks indeksite poolt soovitatavaid parameetreid mudelisse lisada ühekaupa. Otsime tabelist üles kõige suurema modifikatsiooniindeksiga parameetri ja lisame selle mudelisse. Praegusel juhul on kõige suurema indeksiga parameeter 30. real, mis ütleb, et kui lubaksime muutujal x9 laaduda ka faktorile visual võiks oodata et hii-ruut kahaneb umbes 36 võrra. Teeme uue mudeli, milles see muutus on sisse viidud.

HS.model2 <- "visual =~ x1 + x2 + x3 + x9
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9"
fit2 <- cfa(HS.model2, data=HolzingerSwineford1939)
summary(fit2, fit.measures=TRUE)

Selle mudeli sobitusastme näitajad on läinud natuke paremaks. CFI ja SRMR on isegi ületanud selle väärtuse piiri, mida võib pidada heaks, teised näitajad aga mitte. Nüüd vaatame selle uue mudeli modifikatsiooniindekseid ja otsime sealt üles kõige suurema indeksiga parameetri.

modindices(fit2)

Kõige suurem mi väärtus on real 60, mis ütleb, et lisades mudelisse muutujate x3 ja x5 jääkdispersioonide vahelise korrelatsiooni, võiksime oodata hii-ruudu kahanemist umbes 8.6-e võrra. Lisame jälle uue parameetri mudelisse.

HS.model3 <- "visual =~ x1 + x2 + x3 + x9
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9
x3 ~~ x5"
fit3 <- cfa(HS.model3, data=HolzingerSwineford1939)

summary(fit3, fit.measures=TRUE)

Sobitusastmenäitajad on läinud veel natuke paremaks. Nüüd on ka TLI üle 0.95-e, kuid hii-ruut ja RMSEA veel väga head sobivust ei näita. Uurime jällegi uue mudeli modifiktatsiooniindekseid ja otsime üles suurima.

modindices(fit3)

Suurim väärtus tulbas mi on real 71, mis ütleb, et lisades mudelisse muutujate x4 ja x7 jääkdispersioonide vahelise korrelatsiooni, võiksime oodata hii-ruudu kahanemist umbes 8 võrra. Lisame taas uue parameetri mudelisse.

HS.model4 <- "visual =~ x1 + x2 + x3 + x9
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9
x3 ~~ x5
x4 ~~ x7"
fit4 <- cfa(HS.model4, data=HolzingerSwineford1939)

summary(fit4, fit.measures=TRUE)

Sobitusastme näitajad on jällegi natuke paranenud. Nüüd on ka RMSEA alla 0.05, kuid hii-ruut pole jätkuvalt veel head taset saavutanud. Põhimõtteliselt võiksime parameetrite lisamise protsessi samal viisil jätkata. Selle kohta, kas parameetrite lisamine muutis sobitusastet oluliselt paremaks, saame arvutada ka p-väärtuse. Seda saame teha funktsiooni anova abil, millele anname ette mudelid, mida omavahel võrrelda soovime. Seda funktsiooni saame aga kasutada ainult hierarhiliste mudelite võrdlemiseks. Hierarhilised on mudelid siis, kui üks mudel sisaldub teises ja me saame ühe teisest tuletada vabade parameetrite lisamise teel.

anova(fit3, fit4)

## Chi Square Difference Test
## 
##      Df    AIC    BIC  Chisq Chisq diff Df diff Pr(>Chisq)   
## fit4 21 7473.3 7562.3 35.144                                 
## fit3 22 7479.5 7564.8 43.325     8.1817       1   0.004232 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Antud juhul näeme tabeli kõige parempoolsemast tulbast, et p-väärtus on väike (0.004 ehk <0.01), mis tähendab, et neljas mudel on kolmandast parem. Mudelite indeksite alusel modifitseerimisele tuleks üldiselt läheneda väga ettevaatlikult. Mudeli märkimisväärse muutmisega kaasneb oht, et me ülesobitame mudeli just sellele konkreetsele andmestikule ja leitud mudelit ei ole hiljem võimalik üldistada teistele valimitele. Igal valimi juures on mõningad omapärasused ja kui me oma mudelit mingi valimi andmete põhjal muudame, on oht, et muutused kajastavad just neid omapärasusi ning mitte üldiselt esinevaid seoseid meid huvitavate muutujate vahel.

Ülesanne

1. Andmefailis “suur.viisik” on 30 Suure Viisiku kitsama alaomaduse skoori ehk iga isiksuseomaduse kohta 6 alaomadust. Tehke selle andmestiku kohta 5-faktoriline kinnitava faktoranalüüusi mudel, nii nagu Suure Viisiku teooria seda ette näeb. Hinnake sellise mudeli sobivust antud andmetele. Mudeli koostamiseks on vaja teada tabelis olevate muutujate nimesid. Need saab kätte funktsiooni* names* abil.

names(suur.viisik)

Uurige ka mudeli parameetrite tabelit. Vaadake iga faktori puhul, millised vaadeldud muutujatest laaduvad faktorile kõige tugevamalt ja millised kõige nõrgemalt? Millised faktorid on omavahel kõige tugevamalt ja millised kõige nõrgemalt korreleeritud? Uurige jääkdispersioonide abil, milliseid vaadeldud muutujad seletab antud mudel kõige halvemini?

LISAD

fit2 <- cfa(mudel1, data=HolzingerSwineford1939, estimator = "GLS")
fit3 <- cfa(mudel1, data=HolzingerSwineford1939, estimator = "WLS")