#Tutorial: http://lavaan.ugent.be/tutorial/ #Manuál: http://users.ugent.be/~yrosseel/lavaan/lavaanIntroduction.pdf #Lavaan cheat sheet: http://jeromyanglim.tumblr.com/post/33556941601/lavaan-cheat-sheet library(lavaan) library(semPlot) library(semTools) library(psych) setwd("C:/Users/JEZEK/Dropbox/!Výuka/PSY_232_532/Lekce/12.CFA") raw_data<-read.csv2("NfCCsousyp.csv") # Preface ----------------------------------------------------------------- # CFA model s různými prvky na ukázku. pretest2015<-raw_data[raw_data$zdroj=="pretest2015",] pretest2016<-raw_data[raw_data$zdroj=="pretest2016",] pretesty<-rbind(pretest2015, pretest2016) data<-pretesty[,2:16] M <-'nfo =~ nfc01+nfc02+nfc03 #need for order prda =~ nfc04+nfc05+nfc06+nfc10+nfc11+nfc12 #predictability, discomfort with ambiguity deci =~ nfc07+nfc08+nfc09 +nfc10 #decisiveness s crossloadující 10 clos =~ nfc13+nfc14+nfc15 +nfc12 #closed-mindedness nfc11~~nfc12 ' fit.M<-cfa(M, data=data, estimator="MLR", missing="FIML") windows(height=8, width=12) semPaths(fit.M, what="path", #std rozliší tloušťkou, eq rozliší barevně constraints whatLabels = "std", #par=nestandardizované rotation=4, intercepts=FALSE, #průměry style="lisrel",residScale = 10, #rozptyly maluj jako latentní, "ram" dává šipkroužky, residScale=délka šipky layout="tree2", #"circle2" někdy zajímavé residuals=TRUE, curve=2, edge.label.cex = 0.6, #velikost písma parametrů sizeLat=4, #velikost latentních kroužků sizeMan=3, #velikost manifestních čtverečků edge.color="black" #ouvrrajd, když nechci barvičky ) #### Konec preface. #Znění položek v1 - kde zdroj je autonomie. items<-c("1. žít dobře uspořádaný život mi sedne", "2. pevný režim pomáhá užívat si života", "3. líbí se mi jasný a uspořádaný způsob života", "4. nerada se vystavuju sictuacím, kde nevím, co mohu očekávat", "5. nerada trávím čas ve společnosti lidí, kteří jsou schopni jednat nepředvídatelně", "6. nemám rád nepředvídatelné situace", "7. obyčejně se mi uleví, jakmile se pro něco rozhodnu", "8. když stojím před problémem, obyvykle se snažím ho co nejrychleji vyřešit", "9. pokud nemůžu najít řešení problému okamžitě, jsem netrpělivý a podrážděný", "10.nemám rád nejisté situace", "11.je mi nepříjemné, když nechápu důvod nějaké události, která se mi přihodila", "12.nemám rád, když nějaký výrok může znamenat více věcí", "13.obvykle mě dráždí, když jeden člověk nesouhlasí s něčím, co si myslí všichni ostatní", "14.obvykle se na věci dívám z mnoha různých úhlů pohledu, než si na ně udělám svůj názor", "15.nemám rád otázky, na které lze odpovědět mnoha různými způsoby") #Znění položek v2 -kde zdroj je pretest. items<-c("žít dobře uspořádaný život mi sedne", "pevný režim pomáhá užívat si života", "líbí se mi jasný a uspořádaný způsob života", "nerada se vystavuju sictuacím, kde nevím, co mohu očekávat", "nerada trávím čas ve společnosti lidí, kteří jsou schopni jednat nepředvídatelně", "nemám rád nepředvídatelné situace", "obyčejně se mi uleví, jakmile se pro něco rozhodnu", "když stojím před problémem, obyvykle se snažím ho co nejrychleji vyřešit", "pokud nemůžu najít řešení problému okamžitě, jsem netrpělivý a podrážděný", "nemám rád nejisté situace", "je mi nepříjemné, když nechápu důvod nějaké události, která se mi přihodila", "nemám rád, když nějaký výrok může znamenat více věcí", "obvykle mě dráždí, když jeden člověk nesouhlasí s něčím, co si myslí všichni ostatní", "nepotřebuji se zabývat všemi možnými úhly pohledu na to, abych si udělal(a) vlastní názor", "nemám rád otázky, na které lze odpovědět mnoha různými způsoby") #data pretest2015<-raw_data[raw_data$zdroj=="pretest2015",] pretest2016<-raw_data[raw_data$zdroj=="pretest2016",] pretesty<-rbind(pretest2015, pretest2016) autP1<-raw_data[raw_data$zdroj=="autonomie_P1",] autP9<-raw_data[raw_data$zdroj=="autonomie_P9",] autonomie<-rbind(autP1, autP9) # Výběr dat --------------------------------------------------------------- data<-pretesty[,2:16] # Popisné statistiky ------------------------------------------------------ #Popisné statistiky položek describe(data) summary(as.data.frame(lapply(data, ordered))) #za běhu prohlásit proměnné za ordinální #Vstupní korelační a kovarianční matice C<-cor(data, use="complete.obs") V<-cov(data, use="complete.obs") #Symbolické zobrazení korelační matice symnum(C, cutpoints = c(-1, 0, 0.2, 0.4, 0.6, 0.8, 1)) # Definice modelů --------------------------------------------------------- #Model 1 - jednoduchý, jednodimenzionální M1<-'nfcc =~ nfc01+nfc02+nfc03+nfc04+nfc05+nfc06+nfc07+nfc08+nfc09+nfc10 +nfc11+nfc12+nfc13+nfc14+nfc15' #Model 2 - jednoduchý, pětifaktorový M2<-'nfo =~ nfc01+nfc02+nfc03 #need for order pred =~ nfc04+nfc05+nfc06 #predictability deci =~ nfc07+nfc08+nfc09 #decisiveness disa =~ nfc10+nfc11+nfc12 #discomfort with ambiguity clos =~ nfc13+nfc14+nfc15 #closed-mindedness ' #Model 3 - pětifaktorový s korelací reziduí M3<-'nfo =~ nfc01+nfc02+nfc03 pred =~ nfc04+nfc05+nfc06 deci =~ nfc07+nfc08+nfc09 disa =~ nfc10+nfc11+nfc12 clos =~ nfc13+nfc14+nfc15+nfc12 #crossloading 12 nfc01 ~~ nfc02 #korelace mezi reziduálními rozptyly ' #Model 4 - jednoduchý, čtyřfaktorový M4<-'nfo =~ nfc01+nfc02+nfc03 #need for order prda =~ nfc04+nfc05+nfc06+nfc10+nfc11+nfc12 #predictability, discomfort with ambiguity deci =~ nfc07+nfc08+nfc09 #decisiveness clos =~ nfc13+nfc14+nfc15 #closed-mindedness ' #Model 4a - čtyřfaktorový s komplikacemi M4a<-'nfo =~ nfc01+nfc02+nfc03 #need for order prda =~ nfc04+nfc05+nfc06+nfc10+nfc11+nfc12 #predictability, discomfort with ambiguity deci =~ nfc07+nfc08+nfc09 +nfc10 #decisiveness s crossloadující 10 clos =~ nfc13+nfc14+nfc15 +nfc12 #closed-mindedness nfc11~~nfc12 ' #Model 5 - bifaktor, pětifaktorový - ortogonalita faktorů je zajištěna voláním cfa #První položka už nenastavuje škálu,místo toho se fixuje rozptyl latentních. M5<-'nfo =~ NA*nfc01+nfc02+nfc03 pred =~ NA*nfc04+nfc05+nfc06 deci =~ NA*nfc07+nfc08+nfc09 disa =~ NA*nfc10+nfc11+nfc12 clos =~ NA*nfc13+nfc14+nfc15 nfcc =~ NA*nfc01+nfc02+nfc03+nfc04+nfc05+nfc06+nfc07+nfc08+nfc09+nfc10 +nfc11+nfc12+nfc13+nfc14+nfc15 nfo ~~ 1*nfo pred ~~ 1*pred deci ~~ 1*deci disa ~~ 1*disa clos ~~ 1*clos nfcc ~~ 1*nfcc ' # Odhad parametrů --------------------------------------------------------- fit.M1<-cfa(M1, data=data, estimator="MLR", #(ML: MLR, MLM) (DWLS: WLS, WLSM, WLSMV) missing="FIML" # alternativou je "listwise" ) fit.M2<-cfa(M2, data=data, estimator="MLR", missing="FIML") fit.M3<-cfa(M3, data=data, estimator="MLR", missing="FIML") fit.M4<-cfa(M4, data=data, estimator="MLR", missing="FIML") fit.M4a<-cfa(M4a, data=data, estimator="MLR", missing="FIML") fit.M5<-cfa(M5, data=data, estimator="MLR", missing="FIML", orthogonal=TRUE) #Pokud chci prohlásit položky za ordinální, tak takto. fit.M1o<-cfa(M1, data=data, ordered = names(data)[2:16], estimator="WLSMV") fit.M2o<-cfa(M2, data=data, ordered = names(data)[2:16], estimator="WLSMV") fit.M5o<-cfa(M3, data=data, ordered = names(data)[2:16], estimator="WLSMV", orthogonal=TRUE) ##Warnings about empty cells-see http://www.tandfonline.com/doi/abs/10.1080/10705511.2011.557339 # Výsledky ---------------------------------------------------------------- #Stručný souhrn výsledků. summary(fit.M1, standardized = TRUE, fit.measures = TRUE, ci=TRUE) summary(fit.M2, standardized = TRUE, fit.measures = TRUE, ci=TRUE) summary(fit.M3, standardized = TRUE, fit.measures = TRUE, ci=TRUE) summary(fit.M4, standardized = TRUE, fit.measures = TRUE, ci=TRUE) summary(fit.M4a, standardized = TRUE, fit.measures = TRUE, ci=TRUE) summary(fit.M5, standardized = TRUE, fit.measures = TRUE, ci=TRUE) #Přehled odhadnutých parametrů - nepřehledný, technický. parTable(fit.M3o) #Všechny indikátory shody počítané lavaanem. fitmeasures(fit.M4a) #A ještě nějaké navíc. moreFitIndices(fit.M1) #V případě nízkého nullRMSEA vyplicne i odkaz na Kennyho nullRMSEA(fit.M1) #Spojení fitů do matice - bude jich hodně. fity<-cbind(fitmeasures(fit.M1)) fity<-cbind(fity,fitmeasures(fit.M2)) fity<-cbind(fity,fitmeasures(fit.M3)) fity<-cbind(fity,fitmeasures(fit.M4)) fity<-cbind(fity,fitmeasures(fit.M4a)) fityO<-cbind(fitmeasures(fit.M1o)) #ordinální fity mají jinou strukturu fityO<-cbind(fityO,fitmeasures(fit.M2o)) # Diagram ----------------------------------------------------------------- windows(height=8, width=12) semPaths(fit.M4a, what="path", #std rozliší tloušťkou, eq rozliší barevně constraints whatLabels = "std", #par=nestandardizované rotation=4, intercepts=FALSE, #průměry style="lisrel",residScale = 10, #rozptyly maluj jako latentní, "ram" dává šipkroužky, residScale=délka šipky layout="tree2", #"circle2" někdy zajímavé residuals=TRUE, curve=2, edge.label.cex = 0.6, #velikost písma parametrů sizeLat=4, #velikost latentních kroužků sizeMan=3, #velikost manifestních čtverečků edge.color="black" #ouvrrajd, když nechci barvičky ) semPaths(fit.M2, "std", rotation=1, style="lisrel", layout="tree2") semPaths(fit.M2, "std", rotation=1, style="lisrel", layout="circle2") semPaths(fit.M3o, "std", rotation=1, style="lisrel", layout="tree2", intercepts=FALSE, residuals = FALSE, bifactor="nfcc") # Detaily řešení extrahované inspectem ------------------------------------ #Modelem implikovaná kovarianční matice. fitted(fit.M1) #Vstupní kovarianční matice a vektor průměrů inspect(fit.M1, 'sampstat')$cov inspect(fit.M1, 'sampstat')$mean #Koeficienty inspect(fit.M4a, "coefficients") #vše inspect(fit.M1, "coefficients")$lambda #náboje inspect(fit.M1, "coefficients")$theta #residua inspect(fit.M2, "coefficients")$psi #kovariance faktorů cov2cor(inspect(fit.M2, "coefficients")$psi) #korelace faktorů #Residuální matice. residuals(fit.M2) residuals(fit.M2, type = "standardized") #Nahlédnutí do objektu lavaan. Defaultně dá přehled parametrů. inspect(fit.M2) #Fity (ty samé, co už tu byly) inspect(fit.M, "fit") #selmTools::clipboard vykopíruje, cokoli, co inspect, do schránky clipboard(fit.M4, "summary") # Modifikační indexy. ----------------------------------------------------- #modindices produkuje dataframe, aby usnadnil filtrování. Zde ala Mplus. mi.fit.M4<-modindices(fit.M4) mi.fit.M4[mi.fit.M4$mi>10,] #Vyfiltruj MI>10 # Porovnávání modelů ------------------------------------------------------ #Jednoduchý LR chisquare test pro vnořené modely. anova(fit.M4, fit.M4a) #semTools::compareFit https://rdrr.io/cran/semTools/ fitdif<-compareFit(fit.M4, fit.M4a, nested=TRUE) summary(fitdif) #Neparameticky pomocí lavaan::bootstrapLRT bootstrapLRT(h0=fit.M4, h1=fit.M4a, R=1000, parallel="multicore", ncpus=4, type="ordinary", verbose=FALSE, return.LRT=TRUE) # Reliability ------------------------------------------------------------- #semTools počítají reliabilitu - Raykov rhó, McDonald omega reliability(fit.M4a) # Multivariate normality -------------------------------------------------- # Normalita jen zřídka naplněna, multivariační zvlášť. mardia(pretesty[2:16]) # Multigroup -------------------------------------------------------------- #Spustit řádky 5-11 a 85-88 #17. sloupec je "student_psy" data<-pretesty[,2:17] #Defaultně se multigroup model odhaduje tak, že všechny parametry se mohou mezi skupinami lišit. fit.M4a.g1<-cfa(M4a, data=data, estimator="MLR", missing="FIML", group="student_psy") summary(fit.M4a.g1, standardized = TRUE, fit.measures = TRUE, ci=TRUE) windows(height=8, width=12) #nutno zapnout "recording" v okně semPaths(fit.M4a.g1, what="path", #std rozliší tloušťkou, eq rozliší barevně constraints whatLabels = "std", #par=nestandardizované rotation=4, intercepts=FALSE, #průměry style="lisrel",residScale = 10, #rozptyly maluj jako latentní, "ram" dává šipkroužky, residScale=délka šipky layout="tree2", #"circle2" někdy zajímavé residuals=TRUE, curve=2, edge.label.cex = 0.6, #velikost písma parametrů sizeLat=4, #velikost latentních kroužků sizeMan=3, #velikost manifestních čtverečků edge.color="black" #ouvrrajd, když nechci barvičky ) mi.fit.M4a.g1<-modindices(fit.M4a.g1) mi.fit.M4a.g1[mi.fit.M4a.g1$mi>5,] #Vyfiltruj MI>10 #Na základě rozdílů mezi skupinami můžeme na zkoušku "vyladit" specifikaci modelu v každé skupině zvlášť. #Model 4b - čtyřfaktorový s komplikacemi, upravený pro každou skupinu zvlášť M4b<-'nfo =~ nfc01 +nfc02 +nfc03 prda =~ nfc04 +nfc05 +nfc06+nfc10+nfc11+c(NA,0)*nfc07 deci =~ nfc07 +nfc08 +nfc09 +c(NA,0)*nfc10+c(0,NA)*nfc11 #loading nfc10 je v 1. skupině free, v druhé skupině fixed clos =~ nfc13 +nfc14 +nfc15 +nfc12 nfc11~~nfc12 nfc06 ~~ c(NA,0)*nfc14 nfc09 ~~ c(0,NA)*nfc12 ' fit.M4b.g1<-cfa(M4b, data=data, estimator="MLR", missing="FIML", group="student_psy") summary(fit.M4b.g1, standardized = TRUE, fit.measures = TRUE, ci=TRUE) mi.fit.M4b.g1<-modindices(fit.M4b.g1) mi.fit.M4b.g1[mi.fit.M4b.g1$mi>5,] #Vyfiltruj MI>10 windows(height=8, width=12) #nutno zapnout "recording" v okně semPaths(fit.M4b.g1, what="path", #std rozliší tloušťkou, eq rozliší barevně constraints whatLabels = "std", #par=nestandardizované rotation=4, intercepts=FALSE, #průměry style="lisrel",residScale = 10, #rozptyly maluj jako latentní, "ram" dává šipkroužky, residScale=délka šipky layout="tree2", #"circle2" někdy zajímavé residuals=TRUE, curve=2, edge.label.cex = 0.6, #velikost písma parametrů sizeLat=4, #velikost latentních kroužků sizeMan=3, #velikost manifestních čtverečků edge.color="black" #ouvrrajd, když nechci barvičky ) #Výsledkem je, že u stucdentů psychologie jsou faktory definovány jinak než u ostatních :-( #Jenže takhle se to nedělá (nebo dělat nemá). M4b je výsledkem rybaření na malých skupinách. #Hypotéza, kterou máme touhu podpořit, je jestli je struktura v obou skupinách stejná - INVARIANCE/EKVIVALENCE/DIF. #Co to ale znamená? No v prvé řadě to, jestli lze říci, že faktorové náboje jsou stejné (a ne lovit mezi nimi rozdíly). #Jednotlivě jde shodu loadingů nastavit v definici modelu (př. nfo =~ c(a1,a1)*nfc01 +c(a2,a2)*nfc02 +c(a3,a3)*nfc03 ) M4c<-'nfo =~ nfc01 +c(L2,L2)*nfc02 +c(L3,L3)*nfc03 prda =~ nfc04 +nfc05 +nfc06+nfc10+nfc11+c(NA,0)*nfc07 deci =~ nfc07 +nfc08 +nfc09 +c(NA,0)*nfc10+c(0,NA)*nfc11 #loading nfc10 je v 1. skupině free, v druhé skupině fixed clos =~ nfc13 +nfc14 +nfc15 +nfc12 nfc11~~nfc12 nfc06 ~~ c(NA,0)*nfc14 nfc09 ~~ c(0,NA)*nfc12 ' fit.M4c.g1<-cfa(M4c, data=data, estimator="MLR", missing="FIML", group="student_psy") summary(fit.M4c.g1, standardized = TRUE, fit.measures = TRUE, ci=TRUE) #Ale spíš to děláme po skupinách, typech parametrů, a to při volání odhadu parametrů. #Zde nastavíme náboje náboje položek na faktorech jako stejné napříč skupinami - WEAK invariance. fit.M4a.g2<-cfa(M4a, data=data, estimator="MLR", missing="FIML", group="student_psy", group.equal="loadings") #A porovnáme fit modelu, kde se parametry mohou lišit, s modelem, kde se lišit nemohou. summary(compareFit(fit.M4a.g1, fit.M4a.g2, nested = TRUE)) #2 důvody, proč se modely meliší - malé skupiny (i zde power!) nebo oba fitují blbě #Jinak je ale shoda mezi modely s free loadingy a equal loadingy dobrou zprávou. #Podívejme se, jak vypadá model s equal loadingy. #Vidíme, že podmínka(constraint) platí pro nestandardizované koeficienty. summary(fit.M4a.g2, standardized = TRUE, fit.measures = TRUE, ci=TRUE) #STRONG invariance znamená shodu nejen loadingů, ale i průsečíků ("obtížnosti"). fit.M4a.g3<-cfa(M4a, data=data, estimator="MLR", missing="FIML", group="student_psy", group.equal=c("loadings","intercepts")) summary(compareFit(fit.M4a.g1, fit.M4a.g2, fit.M4a.g3, nested = TRUE)) summary(fit.M4a.g3, standardized = TRUE, fit.measures = TRUE, ci=TRUE) #Model s shodnými i průsečíky je sice signifikantně horší, ale prizmatem některých ukazatelů, # které penalizují složitost modelu je vlastně nejlepší. #STRICT invariance znamená zafixovat navíc i reziduální rozptyly položek. fit.M4a.g4<-cfa(M4a, data=data, estimator="MLR", missing="FIML", group="student_psy", group.equal=c("loadings","intercepts","residuals")) summary(compareFit(fit.M4a.g1, fit.M4a.g2, fit.M4a.g3, fit.M4a.g4, nested = TRUE)) summary(fit.M4a.g4, standardized = TRUE, fit.measures = TRUE, ci=TRUE) #Dále lze hromadně fixovat ještě: # intercepts: the intercepts of the observed variables # means: the intercepts/means of the latent variables # residuals: the residual variances of the observed variables # residual.covariances: the residual covariances of the observed variables # lv.variances: the (residual) variances of the latent variables # lv.covariances: the (residual) covariances of the latent varibles # regressions: all regression coefficients in the model #Hledání příčin non-invariance #Řádky s existujícími loadingy ukazují na možný efekt uvolnění equality constraint mi.fit.M4a.g3<-modindices(fit.M4a.g3) #Pak lze testovat, zda uvolnění podmínky u některé položky model nezlepší. fit.M4a.g5<-cfa(M4a, data=data, estimator="MLR", missing="FIML", group="student_psy", group.equal=c("loadings","intercepts","residuals"), group.partial = c("nfo =~ nfc02","nfc02~1", "nfc02~~nfc02")) #uvolnění loadingu, interceptu a rezidua summary(fit.M4a.g5, standardized = TRUE, fit.measures = TRUE, ci=TRUE) # Measurement invariance -------------------------------------------------- measurementInvariance(M4a, data=data, group="student_psy")