Updating demo scripts

demos/2_mi/mi_demo.R

@@ -1,7 +1,7 @@
 ### Title:    Missing Data in R: MI Demonstration
 ### Author:   Kyle M. Lang
 ### Created:  2015-10-04
-### Modified: 2023-01-31
+### Modified: 2024-01-30
 
 rm(list = ls(all = TRUE)) # Clear workspace
 
@@ -39,7 +39,7 @@ hist(cover)
 cover[cover < 1] %>% range()
 
 ## How many unique response patterns?
-md.pattern(missData) %>% nrow() - 1
+md.pattern(missData, plot = FALSE) %>% nrow() - 1
 
 
 ###-Missing Data Imputation--------------------------------------------------###
@@ -74,7 +74,7 @@ miceOut <- mice(missData, m = 10, maxit = 20, method = meth, seed = 235711)
 
 ###--------------------------------------------------------------------------###
 
-## Use a dummy run of mice() to initialize the meta data:
+## Use a dummy run of mice() to initialize the metadata:
 init <- mice(missData, maxit = 0)
 
 ls(init)
@@ -130,6 +130,7 @@ Rhat.mice(miceOut)
 plot(miceOut)
 
 ## Construct density plots of imputed vs. observed:
+p <- list()
 targets <- names(pm)[pm > 0]
 for(v in targets) 
   p[[v]] <- ggmice(miceOut, aes(.data[[v]], group = .imp)) + geom_density()

demos/3_fiml/fiml_demo.R

@@ -1,14 +1,15 @@
 ### Title:    Missing Data in R: FIML Demonstration
 ### Author:   Kyle M. Lang
 ### Created:  2015-10-26
-### Modified: 2023-01-31
+### Modified: 2024-01-30
 
 rm(list = ls(all = TRUE))
 
 dataDir <- "data/"
 
 library(lavaan)
 library(semTools)
+library(dplyr)
 
 ## Load the data:
 compData <- readRDS(paste0(dataDir, "adams_klps_data-synthetic.rds"))
@@ -20,7 +21,7 @@ missData <- readRDS(paste0(dataDir, "adams_klps_data-incomplete.rds"))
 ## Simple CFA:
 cfaMod1 <- "
 indRac =~ riae1 + riae4 + riae5 + riae6 + riae10
-policy =~ policy1 + policy2 + policy3 + policy4 + policy5 + policy6
+policy =~ policy1 + policy3 + policy4 + policy5 + policy6
 "
 
 ## Estimate the model:
@@ -30,15 +31,15 @@ cfaFit1 <- cfa(model   = cfaMod1,
                missing = "fiml")
 
 summary(cfaFit1)
-fitMeasures(cfaFit1)
+(fm1 <- fitMeasures(cfaFit1))
 
 
 ##----------------------------------------------------------------------------##
 
 ## CFA w/ auxiliary variables:
 cfaMod2 <- "
 indRac =~ riae1 + riae4 + riae5 + riae6 + riae10
-policy =~ policy1 + policy2 + policy3 + policy4 + policy5 + policy6
+policy =~ policy1 + policy3 + policy4 + policy5 + policy6
 
 # Correlate the auxiliary variables with each other and all residuals:
 nori1 ~~
@@ -49,7 +50,6 @@ riae5 +
 riae6 +
 riae10 +
 policy1 +
-policy2 +
 policy3 +
 policy4 +
 policy5 +
@@ -62,7 +62,6 @@ riae5 +
 riae6 +
 riae10 +
 policy1 +
-policy2 +
 policy3 +
 policy4 +
 policy5 +
@@ -78,25 +77,77 @@ cfaFit2 <- cfa(model   = cfaMod2,
 summary(cfaFit2)
 fitMeasures(cfaFit2)
 
+### NOTE: The incremental fit measures above are computed with the wrong
+###       baseline model.
+
+
+##----------------------------------------------------------------------------##
+
+## Define the model syntax for the correct baseline model:
+baseMod <- "
+nori1 ~~
+nori4 +
+riae1 +
+riae4 +
+riae5 +
+riae6 +
+riae10 +
+policy1 +
+policy3 +
+policy4 +
+policy5 +
+policy6
+
+nori4 ~~
+riae1 +
+riae4 +
+riae5 +
+riae6 +
+riae10 +
+policy1 +
+policy3 +
+policy4 +
+policy5 +
+policy6
+"
+
+## Estimate the baseline model:
+baseFit <- cfa(baseMod, data = missData, missing = "FIML")
+summary(baseFit)
+
+## Compute CFA fit measures using the correct baseline model:
+(fm2 <- fitMeasures(cfaFit2, baseline.model = baseFit))
+
 
 ##----------------------------------------------------------------------------##
 
 ## Compare fit of the auxiliary and the no-auxiliary models:
-fitMeasures(cfaFit2) - fitMeasures(cfaFit1)
+(fm2 - fm1)[1:8] %>% round(4)
 
 ## Try the same comparison without any missing data:
-compFit1 <- cfa(model = cfaMod1, data = compData, std.lv = TRUE)
-compFit2 <- cfa(model = cfaMod2, data = compData, std.lv = TRUE)
+compFit1 <- cfa(model         = cfaMod1,
+                data          = compData, 
+                std.lv        = TRUE,
+                meanstructure = TRUE)
+compFit2 <- cfa(model         = cfaMod2,
+                data          = compData, 
+                std.lv        = TRUE,
+                meanstructure = TRUE)
+
+compFm2 <- 
+  fitMeasures(compFit2,
+              baseline.model = cfa(baseMod, data = compData)
+             )
 
-fitMeasures(compFit2) - fitMeasures(compFit1)
+(compFm2 - fitMeasures(compFit1))[1:8] %>% round(4)
 
 
 ###-Structural Equation Model------------------------------------------------###
 
 ## Simple SEM:
 semMod1 <- "
 indRac =~ riae1 + riae4 + riae5 + riae6 + riae10
-policy =~ policy1 + policy2 + policy3 + policy4 + policy5 + policy6
+policy =~ policy1 + policy3 + policy4 + policy5 + policy6
 
 policy + indRac ~ polv
 "
@@ -115,7 +166,7 @@ summary(semFit1, fit.measures = TRUE, fmi = TRUE)
 ## SEM w/ auxiliary variables:
 semMod2 <- "
 indRac =~ riae1 + riae4 + riae5 + riae6 + riae10
-policy =~ policy1 + policy2 + policy3 + policy4 + policy5 + policy6
+policy =~ policy1 + policy3 + policy4 + policy5 + policy6
 
 policy + indRac ~ polv
 
@@ -128,7 +179,6 @@ riae5 +
 riae6 +
 riae10 +
 policy1 +
-policy2 +
 policy3 +
 policy4 +
 policy5 +
@@ -142,7 +192,6 @@ riae5 +
 riae6 +
 riae10 +
 policy1 +
-policy2 +
 policy3 +
 policy4 +
 policy5 +
@@ -234,6 +283,7 @@ summary(paFit3, fit.measures = TRUE, fmi = TRUE)
 inspect(cfaFit3, "coef")$lambda - inspect(cfaFit2, "coef")$lambda
 inspect(cfaFit3, "theta") - inspect(cfaFit2, "theta")
 inspect(cfaFit3, "coef")$psi - inspect(cfaFit2, "coef")$psi
+fitMeasures(cfaFit3) - fm2
 
 
 ###-END----------------------------------------------------------------------###