# To install the latest lavaan version not yet released on CRAN. install.packages("lavaan", repos="http://www.da.ugent.be", type = "source") #To load the installed lavaan. library(lavaan) # Download the R file "ESS5Police.RData" where the data are saved. Assuming that # you have saved in the folder "C:\PL", load it on the R workspace. load("C:\\PL\\ESS5Police.RData") #To print the first rows of the data. head(PoliceDataRc) ####### Section 5.4 ########## # To fit the model. # Extract the GB data from all-country data and save it under the name # "PoliceDataRcGB" by giving the command below. PoliceDataRcGB <- PoliceDataRc[PoliceDataRc$cntry=="GB",] # Note that lavaan requires the data to be saved as data frame which is the case # for "PoliceDataRcGB". To confirm use the command below. is.data.frame(PoliceDataRcGB) # If a data set you wish to analyze with lavaan is not data frame, use the # command as.data.frame(). # Specify that all variables are ordinal except of course for the variables # "idno" and "cntry" in the first two columns of the data. Save the new format # of the data under the name "PoliceDataRcGBOrd". PoliceDataRcGBOrd <- PoliceDataRcGB PoliceDataRcGBOrd[,2:17] <- lapply(PoliceDataRcGBOrd[,2:17], ordered) # Specify the model to be fitted and save it as "Ex1Model". Ex1Model <- ' #Measurement part of the model TrEf =~ plcpvcr + plccbrg + plcarcr TrFa =~ plcrspc + plcfrdc + plcexdc ObOb =~ bplcdc + doplcsy + dpcstrb MoAl =~ plcrgwr + plcipvl + gsupplc WiCo =~ caplcst + widprsn + wevdct #Structural part of the model ObOb ~ TrEf + TrFa MoAl ~ TrEf + TrFa WiCo ~ TrEf + ObOb + MoAl TrEf ~~ TrFa #Cov(TrEf, TrFa) to be estimated ObOb ~~ MoAl #Cov(zeta3, zeta4) to be estimated ' # Fit the model using the function sem. Specify PL as the estimation method by # setting the input argument "estimator" equal to "PML". The argument # "std.lv = TRUE" fixes the variances of TrEf and TrFa, and the variances of the # residuals of ObOb, MoAl, and WiCo to 1 to define the unit of the factor scales. # The argument "verbose = TRUE" prints the progress of the computations on the # R session. If you do not wish to do so, you can omit it. Ex1FittedModel <- sem(model = Ex1Model, data = PoliceDataRcGBOrd, estimator = "PML", std.lv = TRUE, verbose = TRUE) # To print the output of the fitted model in the R session give the command # below. summary(Ex1FittedModel) # To save the parameter estimates, standard errors, p-values, and their # 95% confidence interval as a data frame, use the command below. This enables # you to select certain values for further process. Ex1FittedModel_ParEst <- parameterEstimates(Ex1FittedModel) # To obtain the PLRT for overall fit (where the thresholds are nuisance # parameters). Ex1FittedModel # Alternative ways to obtain the same results: summary(Ex1FittedModel) # and lavaan:::ctr_pml_plrt(Ex1FittedModel) ############################################################################## ####### Section 5.5 ########## # Specify the model. Ex2_MeasEquivModel <- ' #Measurement part of the model TrEf =~ plcpvcr + plccbrg + plcarcr TrFa =~ plcrspc + plcfrdc + plcexdc ObOb =~ bplcdc + doplcsy + dpcstrb MoAl =~ plcrgwr + plcipvl + gsupplc WiCo =~ caplcst + widprsn + wevdct #Structural part of the model ObOb ~ TrEf + TrFa MoAl ~ TrEf + TrFa WiCo ~ TrEf + ObOb + MoAl TrEf ~~ TrFa ObOb ~~ MoAl TrEf ~~ c(1,NA)*TrEf TrFa ~~ c(1,NA)*TrFa ObOb ~~ c(1,NA)*ObOb MoAl ~~ c(1,NA)*MoAl WiCo ~~ c(1,NA)*WiCo ' # The vector c(1,NA) specifies that the variance of each factor is fixed to 1 in # the first group, here GB, and is free to be estimated in the other group, # here IE. In general, the vector c() has as many elements as the number of # groups which are ordered in alphabetical order. We have added this # specification in the model syntax because, when we fit the model below # using the function sem, we specify 'std.lv = TRUE' which fixes the # variances of the factors to 1 in both groups. So, we override this # specification by adding the vector c(1,NA) above. In turn, the argument # 'std.lv = TRUE' needs to be added in the sem function below to override # the default setting which is to fix the loading of one indicator of each # factor to 1 (used to determine the factor scale unit). No need to add such # a vector for the factor means and the variances of the underlying # variables because these parameters are fixed to 0 and 1, respectively, # only in the first group by default. # Extract the GB and IE data from the full data set and save them under the name # PoliceDataRcGBIE, for example. PoliceDataRcGBIE <- PoliceDataRc[PoliceDataRc$cntry=="GB"| PoliceDataRc$cntry=="IE",] # Recall that the data is already saved as data frame. To confirm: is.data.frame(PoliceDataRcGBIE) # Specify that all the variables, except for "idno" and "cntry", are ordinal. # Save the new format of the data under the name PoliceDataRcGBIEOrd. PoliceDataRcGBIEOrd <- PoliceDataRcGBIE PoliceDataRcGBIEOrd[,3:17] <- lapply(PoliceDataRcGBIEOrd[,3:17], ordered) # Fit the model and save the output as Ex2FittedMeasEquivModel. Note that # here we use two more input arguments, 'group' and 'group.equal'. The first # one is to specify the grouping variable, here 'cntry', and this way, state # that a multi-group SEM should be fitted. With the second argument, we # determine the parameters which cross-country equality constraints will be # imposed on; here, loadings and thresholds. Ex2FittedMeasEquivModel <- sem(model = Ex2_MeasEquivModel, data = PoliceDataRcGBIEOrd, std.lv=TRUE, estimator = "PML", group = "cntry", group.equal = c("loadings","thresholds"), verbose=TRUE) # To print the output give the command below. The results for the GB data are # presented first followed by those for the IE data. Parameters constrained to # be equal are marked by the same label printed in parentheses next to the # parameters. summary(Ex2FittedMeasEquivModel) # Wald test. # First save the estimates of the factor means / interecepts in a vector, # e.g. named 'est_IE_means'. Recall that all parameter estimates can be # retrieved using the function parameterEstimates() with input the name of # the fitted model. Giving the command # parameterEstimates(Ex2FittedMeasEquivModel) we see that the parameters of # interest are reported in lines 334 to 338 of the returned table. Out of # all columns of the table, we only need the one that gives the estimates # labelled 'est'. est_IE_means <- parameterEstimates(Ex2FittedMeasEquivModel)[334:338,"est"] # Next obtain the estimated variance-covariance matrix of all parameter # estimates using the function vcov() and input the name of the fitted model. VCOV_Ex2FittedMeasEquivModel <- vcov(Ex2FittedMeasEquivModel) # From the above matrix, we only need the rows and columns which provide the # variances and covariances of the estimated factor means / intercepts. An # easy way to identify the index of these rows and columns is to call the # function rownames(VCOV_Ex2FittedMeasEquivModel) which gives the labels of # the rows. This way, we see that the row indices are 249 to 253. We save # the submatrix of variances-covariances under the name 'VCOV_IE_means'. VCOV_IE_means <- VCOV_Ex2FittedMeasEquivModel[249:253, 249:253] # We use the vector of the estimates and the matrix of the estimated # variances-covariances to compute the Wald test. Recall that the matrix of # the estimated variances-covariances needs to be inverted and for this, we # use the command solve(). Wald_test <- matrix(est_IE_means, nrow=1) %*% solve(VCOV_IE_means) %*% matrix(est_IE_means, ncol=1) # The value of the Wald test should be compared to a chi-squared distribution # with 5 degrees of freedom since we test 5 parameters simultaneously. # The p-value of the test can be computed as follows. pvalue_Wald_test <- 1-pchisq(Wald_test, df=5) # Print the values. Wald_test pvalue_Wald_test # To obtain the PLRT for overall fit. In this example parametric strucure is # imposed on both the polychoric correlations and the thresholds. lavaan:::ctr_pml_plrt2(Ex2FittedMeasEquivModel) ############################################################################## ####### Section 5.6 ########## # Define the models to be compared. # The model with the minimum set of constraints. Ex2_MinConModel <- ' ##Measurement part of the model TrEf =~ c(c1,c1)*plcpvcr + plccbrg + plcarcr TrFa =~ c(c4,c4)*plcrspc + plcfrdc + plcexdc ObOb =~ c(c7,c7)*bplcdc + doplcsy + dpcstrb MoAl =~ c(c10,c10)*plcrgwr + plcipvl + gsupplc WiCo =~ c(c13,c13)*caplcst + widprsn + wevdct plcpvcr | c(p1,p1)*t1 + c(p2,p2)*t2 + t3 + t4 + t5 + t6 + t7 + t8 + t9 + t10 plccbrg | c(p3,p3)*t1 + t2 + t3 + t4 + t5 + t6 + t7 + t8 + t9 + t10 plcarcr | c(p4,p4)*t1 + t2 + t3 + t4 + t5 + t6 + t7 + t8 + t9 + t10 plcrspc | c(p5,p5)*t1 + c(p6,p6)*t2 + t3 plcfrdc | c(p7,p7)*t1 + t2 + t3 plcexdc | c(p8,p8)*t1 + t2 + t3 bplcdc | c(p9,p9)*t1 + c(p10,p10)*t2 + t3 + t4 + t5 + t6 + t7 + t8 + t9 + t10 doplcsy | c(p11,p11)*t1 + t2 + t3 + t4 + t5 + t6 + t7 + t8 + t9 + t10 dpcstrb | c(p12,p12)*t1 + t2 + t3 + t4 + t5 + t6 + t7 + t8 + t9 + t10 plcrgwr | c(p13,p13)*t1 + c(p14,p14)*t2 + t3 + t4 plcipvl | c(p15,p15)*t1 + t2 + t3 + t4 gsupplc | c(p16,p16)*t1 + t2 + t3 + t4 caplcst | c(p17,p17)*t1 + c(p18,p18)*t2 + t3 widprsn | c(p19,p19)*t1 + t2 + t3 wevdct | c(p20,p20)*t1 + t2 + t3 #Define the unit of the scales of the underlying variables. plcpvcr ~*~ c(1,NA)*plcpvcr plccbrg ~*~ c(1,NA)*plccbrg plcarcr ~*~ c(1,NA)*plcarcr plcrspc ~*~ c(1,NA)*plcrspc plcfrdc ~*~ c(1,NA)*plcfrdc plcexdc ~*~ c(1,NA)*plcexdc bplcdc ~*~ c(1,NA)*bplcdc doplcsy ~*~ c(1,NA)*doplcsy dpcstrb ~*~ c(1,NA)*dpcstrb plcrgwr ~*~ c(1,NA)*plcrgwr plcipvl ~*~ c(1,NA)*plcipvl gsupplc ~*~ c(1,NA)*gsupplc caplcst ~*~ c(1,NA)*caplcst widprsn ~*~ c(1,NA)*widprsn wevdct ~*~ c(1,NA)*wevdct ##Structural part of the model ObOb ~ TrEf + TrFa MoAl ~ TrEf + TrFa WiCo ~ TrEf + ObOb + MoAl TrEf ~~ TrFa ObOb ~~ MoAl #Define the origin of the factor scales. TrEf ~ c(0,NA)*1 TrFa ~ c(0,NA)*1 ObOb ~ c(0,NA)*1 MoAl ~ c(0,NA)*1 WiCo ~ c(0,NA)*1 #Define the unit of the factor scales. TrEf ~~ c(1,NA)*TrEf TrFa ~~ c(1,NA)*TrFa ObOb ~~ c(1,NA)*ObOb MoAl ~~ c(1,NA)*MoAl WiCo ~~ c(1,NA)*WiCo ' # Note that the elements of the vector c() can be either a label, a number, or NA. # Same labels imply equality constraint on the corresponding parameters. A number # implies that the corresponding parameter is fixed to that number, and NA impies # that the corresponding parameter is free to be estimated. Whenever there is no # vector, it is implied that distinct parameters are to be estimated for the groups. # The above model where cross-country equality constraints on all of the # loadings have been added. Ex2_MetricInvModel <- ' ##Measurement part of the model #Only the following 5 lines differ from the previous model syntax. TrEf =~ c(c1,c1)*plcpvcr + c(c2,c2)*plccbrg + c(c3,c3)*plcarcr TrFa =~ c(c4,c4)*plcrspc + c(c5,c5)*plcfrdc + c(c6,c6)*plcexdc ObOb =~ c(c7,c7)*bplcdc + c(c8,c8)*doplcsy + c(c9,c9)*dpcstrb MoAl =~ c(c10,c10)*plcrgwr + c(c11,c11)*plcipvl + c(c12,c12)*gsupplc WiCo =~ c(c13,c13)*caplcst + c(c14,c14)*widprsn + c(c15,c15)*wevdct plcpvcr | c(p1,p1)*t1 + c(p2,p2)*t2 + t3 + t4 + t5 + t6 + t7 + t8 + t9 + t10 plccbrg | c(p3,p3)*t1 + t2 + t3 + t4 + t5 + t6 + t7 + t8 + t9 + t10 plcarcr | c(p4,p4)*t1 + t2 + t3 + t4 + t5 + t6 + t7 + t8 + t9 + t10 plcrspc | c(p5,p5)*t1 + c(p6,p6)*t2 + t3 plcfrdc | c(p7,p7)*t1 + t2 + t3 plcexdc | c(p8,p8)*t1 + t2 + t3 bplcdc | c(p9,p9)*t1 + c(p10,p10)*t2 + t3 + t4 + t5 + t6 + t7 + t8 + t9 + t10 doplcsy | c(p11,p11)*t1 + t2 + t3 + t4 + t5 + t6 + t7 + t8 + t9 + t10 dpcstrb | c(p12,p12)*t1 + t2 + t3 + t4 + t5 + t6 + t7 + t8 + t9 + t10 plcrgwr | c(p13,p13)*t1 + c(p14,p14)*t2 + t3 + t4 plcipvl | c(p15,p15)*t1 + t2 + t3 + t4 gsupplc | c(p16,p16)*t1 + t2 + t3 + t4 caplcst | c(p17,p17)*t1 + c(p18,p18)*t2 + t3 widprsn | c(p19,p19)*t1 + t2 + t3 wevdct | c(p20,p20)*t1 + t2 + t3 #Define the unit of the scales of the underlying variables. plcpvcr ~*~ c(1,NA)*plcpvcr plccbrg ~*~ c(1,NA)*plccbrg plcarcr ~*~ c(1,NA)*plcarcr plcrspc ~*~ c(1,NA)*plcrspc plcfrdc ~*~ c(1,NA)*plcfrdc plcexdc ~*~ c(1,NA)*plcexdc bplcdc ~*~ c(1,NA)*bplcdc doplcsy ~*~ c(1,NA)*doplcsy dpcstrb ~*~ c(1,NA)*dpcstrb plcrgwr ~*~ c(1,NA)*plcrgwr plcipvl ~*~ c(1,NA)*plcipvl gsupplc ~*~ c(1,NA)*gsupplc caplcst ~*~ c(1,NA)*caplcst widprsn ~*~ c(1,NA)*widprsn wevdct ~*~ c(1,NA)*wevdct ##Structural part of the model ObOb ~ TrEf + TrFa MoAl ~ TrEf + TrFa WiCo ~ TrEf + ObOb + MoAl TrEf ~~ TrFa ObOb ~~ MoAl #Define the origin of the factor scales. TrEf ~ c(0,NA)*1 TrFa ~ c(0,NA)*1 ObOb ~ c(0,NA)*1 MoAl ~ c(0,NA)*1 WiCo ~ c(0,NA)*1 #Define the unit of the factor scales. TrEf ~~ c(1,NA)*TrEf TrFa ~~ c(1,NA)*TrFa ObOb ~~ c(1,NA)*ObOb MoAl ~~ c(1,NA)*MoAl WiCo ~~ c(1,NA)*WiCo ' # To fit the models, we use the data saved in "PoliceDataRcGBIEOrd" # defined by the R commands given in the previous section. # To fit the model with the minimum set of constraints. Ex2FittedMinConModel <- sem(model = Ex2_MinConModel, data = PoliceDataRcGBIEOrd, std.lv=TRUE, estimator = "PML", group = "cntry", verbose=TRUE) # To fit the model with the additional cross-country equality constraints on # the loadings. Ex2_FittedMetricInvModel <- sem(model = Ex2_MetricInvModel, data = PoliceDataRcGBIEOrd, std.lv = TRUE, estimator = "PML", group = "cntry", verbose=TRUE) # To print the outputs of the fitted models. summary(Ex2FittedMinConModel) summary(Ex2_FittedMetricInvModel) # To compare the models. lavTestLRT(Ex2_FittedMetricInvModel, Ex2FittedMinConModel) # To obtain the PL-AIC and PL-BIC values of the model with the minimum set of # constraints. lavaan:::ctr_pml_aic_bic(Ex2FittedMinConModel) ############################################################################## ####### Section 5.7 ########## # To fit the model of Section 5.4 using CP or AC. # Fitting the model using the CP. Note that the input argument missing="pairwise" # determines that the CP should be used. The argument test="none" has been added # to save computational time from calculating the default PLRT for overall # fit. Extending the PLRT to the case of CP is still under research. Ex3FittedModel_CP <- sem(model = Ex1Model, data = PoliceDataRcGBOrd, estimator = "PML", missing = "pairwise", std.lv = TRUE, test="none", verbose=TRUE) # Fitting the model using the AC. Here we specify missing = "available.cases". # The argument test="none" has been added for the same reason as above. Ex3FittedModel_AC <- sem(model = Ex1Model, data = PoliceDataRcGBOrd, estimator = "PML", missing = "available.cases", std.lv = TRUE, test="none", verbose=TRUE) # To print the number of observations used in each case. Ex3FittedModel_CP Ex3FittedModel_AC # To fit the multi-group model of Section 5.5 using CP or AC. # Using the CP. Ex4FittedMeasEquivModel_CP <- sem(model = Ex2_MeasEquivModel, data = PoliceDataRcGBIEOrd, std.lv = TRUE, estimator = "PML", missing = "pairwise", group = "cntry", group.equal = c("loadings","thresholds"), test = "none", verbose=TRUE) # Using the AC. Ex4FittedMeasEquivModel_AC <- sem(model = Ex2_MeasEquivModel, data = PoliceDataRcGBIEOrd, std.lv = TRUE, estimator = "PML", missing = "available.cases", group = "cntry", group.equal = c("loadings","thresholds"), test = "none", verbose=TRUE) ############################################################################## ####### Appendix ########## # To obtain the full outputs of the above models as presented in the Appendix. summary(Ex1FittedModel) summary(Ex2FittedMeasEquivModel) summary(Ex2FittedMinConModel) summary(Ex2_FittedMetricInvModel) summary(Ex3FittedModel_CP) summary(Ex3FittedModel_AC)