# Packages ---- library(lavaan) library(semTools) library(dplyr) library(tidyr) library(future) library(future.apply) library(mice) library(lavaan.mi) library(ggplot2) # Helper functions ---- # Compare model fit with AFI ---- compare_fit <- function(fit_1, fit_2, indices = c('cfi')){ fitmeasures(fit_2, fit.measures = indices)[] - fitmeasures(fit_1, fit.measures = indices)[] } # Calculate RMSEAd ---- # based on Savalei et al. 2024 # Input can be a list of models or individual lavaan models # e.g., # compare_nested_models_RMSEA(fit0, fit1) # or # model_list <- list(fit0, fit1, fit2) # compare_nested_models_RMSEA(model_list) compare_nested_models_RMSEA <- function(..., model_names = NULL, CI = .90) { require(lavaan) # Capture all input arguments models_input <- list(...) # Check if the first argument is a list of models if (length(models_input) == 1 && is.list(models_input[[1]]) && all(sapply(models_input[[1]], function(x) inherits(x, "lavaan")))) { models <- models_input[[1]] } else { # Assume that individual models are passed directly models <- models_input } # Check if at least two models are provided if (length(models) < 2) { stop("Please provide at least two nested lavaan models for comparison.") } # Assign default model names if not provided if (is.null(model_names)) { model_names <- paste0("Model_", seq_along(models)) } else { # Check if the number of names matches the number of models if (length(model_names) != length(models)) { stop("Length of model_names must match the number of models provided.") } } # Ensure that all elements in the models list are lavaan objects if (!all(sapply(models, function(x) inherits(x, c("lavaan","lavaan.mi"))))) { stop("All elements in the models list must be lavaan model objects.") } # Extract number of groups from each model n_groups <- sapply(models, function(x) lavaan::inspect(x, "ngroups")) # Check if all models have the same number of groups if(length(unique(n_groups)) != 1){ stop("All models must have the same number of groups.") } G <- n_groups[1] # Number of groups # Internal RMSEA.CI function RMSEA.CI <- function(T, df, N, G, CI){ lower.lambda <- function(lambda) { (pchisq(T, df=df, ncp=lambda) - (1 - (1-CI)/2)) } upper.lambda <- function(lambda) { (pchisq(T, df=df, ncp=lambda) - ((1-CI)/2)) } # RMSEA CI Lower Bound lambda.l <- try(uniroot(f=lower.lambda, lower=0, upper=T)$root, silent=TRUE) if(inherits(lambda.l, "try-error")) { lambda.l <- NA RMSEA.CI.l <- NA } else { if(lambda.l < 0){ RMSEA.CI.l <- 0 } else { RMSEA.CI.l <- sqrt(lambda.l * G / ((N - G) * df)) # Appendix of Savalei paper suggests 1/(N-1) should be replaced with 1/(N-G), which also seems consistent with Table 1, but OSF syntax uses 1/(N-1) } } # RMSEA CI Upper Bound N.RMSEA <- max(N, T * 4) lambda.u <- try(uniroot(f=upper.lambda, lower=0, upper=N.RMSEA)$root, silent=TRUE) if(inherits(lambda.u, "try-error")) { lambda.u <- NA RMSEA.CI.u <- NA } else { if(lambda.u < 0){ RMSEA.CI.u <- 0 } else { RMSEA.CI.u <- sqrt(lambda.u * G / ((N - G) * df)) # See comment on RMSEA.CI.l } } RMSEA.CI <- c(RMSEA.CI.l, RMSEA.CI.u) names(RMSEA.CI) <- c("Lower_CI", "Upper_CI") return(RMSEA.CI) } # Internal pvals function pvals <- function(T, df, N, G){ RMSEA0 <- c(0, 0.01, 0.05, 0.06, 0.08, 0.1) eps0 <- df * RMSEA0^2 / G nonc <- eps0 * (N - G) pvals <- pchisq(T, df=df, ncp=nonc) names(pvals) <- c("RMSEA>0", "RMSEA>.01", "RMSEA>.05", "RMSEA>.06", "RMSEA>.08", "RMSEA>.10") return(pvals) } # Extract fit indices, sample size, and number of parameters from each model fit_indices <- data.frame( Model = model_names, ChiSq = sapply(models, function(x) fitMeasures(x, "chisq")), df = sapply(models, function(x) fitMeasures(x, "df")), RMSEA = sapply(models, function(x) fitMeasures(x, "rmsea")), stringsAsFactors = FALSE ) # Extract sample size (assuming all models have the same sample size) N_values <- lapply(models, function(x) inspect(x, "nobs")) # Check that all models have the same number of groups n_groups <- sapply(N_values, length) if (length(unique(n_groups)) != 1) { stop("All models must have the same number of groups.") } # Check that sample sizes per group are the same across all models for (g in 1:n_groups[1]) { group_sizes <- sapply(N_values, function(n) n[g]) if (length(unique(group_sizes)) != 1) { stop(paste("Sample size for group", g, "differs across models.")) } } # If all checks pass, extract sample sizes for each group N <- nobs(models[[1]]) # should return the same as above # Initialize vectors to store comparison results comparison <- paste(model_names[-length(model_names)], "vs", model_names[-1]) Delta_ChiSq <- vector("numeric", length = length(models) - 1) Delta_df <- vector("numeric", length = length(models) - 1) Delta_RMSEA <- vector("numeric", length = length(models) - 1) RMSEA_d <- vector("numeric", length = length(models) - 1) RMSEA_d_Lower_CI <- vector("numeric", length = length(models) - 1) RMSEA_d_Upper_CI <- vector("numeric", length = length(models) - 1) pvals_list <- vector("list", length = length(models) - 1) # Loop through consecutive model pairs for(i in 2:length(models)) { # Chi-square and df differences Delta_ChiSq[i - 1] <- fit_indices$ChiSq[i] - fit_indices$ChiSq[i - 1] Delta_df[i - 1] <- fit_indices$df[i] - fit_indices$df[i - 1] # G is number of groups # Already extracted before the loop # Compute RMSEAd numerator <- G*(Delta_ChiSq[i - 1] - Delta_df[i - 1]) #Savalei has 1/G in denom which = G* in numerator denominator <- (N - G) * Delta_df[i - 1] # N - G reduces to N - 1 in single-group models if(denominator <= 0) { warning(paste("Non-positive denominator for RMSEAd computation between", model_names[i - 1], "and", model_names[i], ". Setting RMSEAd to NA.")) rmsea_d_value <- NA } else { if(numerator <= 0){ rmsea_d_value <- 0 } else { rmsea_d_value <- sqrt(numerator / denominator) } } RMSEA_d[i - 1] <- rmsea_d_value # Compute RMSEA CI for RMSEAd if(!is.na(rmsea_d_value)) { ci <- RMSEA.CI(T = Delta_ChiSq[i - 1], df = Delta_df[i - 1], N = N, G = G, CI = CI) RMSEA_d_Lower_CI[i - 1] <- ci["Lower_CI"] RMSEA_d_Upper_CI[i - 1] <- ci["Upper_CI"] } else { RMSEA_d_Lower_CI[i - 1] <- NA RMSEA_d_Upper_CI[i - 1] <- NA } # Equivalence tests based on RMSEAd if(!is.na(rmsea_d_value)) { pvals_comp <- pvals(T = Delta_ChiSq[i - 1], df = Delta_df[i - 1], N = N, G = G) pvals_list[[i - 1]] <- round(pvals_comp, 3) } else { pvals_list[[i - 1]] <- rep(NA, 6) names(pvals_list[[i - 1]]) <- c("RMSEA>0", "RMSEA>.01", "RMSEA>.05", "RMSEA>.06", "RMSEA>.08", "RMSEA>.10") } # Compute Delta RMSEA Delta_RMSEA[i - 1] <- fit_indices$RMSEA[i] - fit_indices$RMSEA[i - 1] } # Compile p-values into separate columns pvals_df <- do.call(rbind, pvals_list) rownames(pvals_df) <- NULL # Create the final comparison data frame comparison_df <- data.frame( Comparison = comparison, Delta_ChiSq = Delta_ChiSq, Delta_df = Delta_df, Delta_RMSEA = Delta_RMSEA, RMSEA_d = RMSEA_d, RMSEA_d_Lower_CI = RMSEA_d_Lower_CI, RMSEA_d_Upper_CI = RMSEA_d_Upper_CI, #p_RMSEA_gt0 = pvals_df[, "RMSEA>0"], #p_RMSEA_gt0.01 = pvals_df[, "RMSEA>.01"], p_RMSEA_gt0.05 = pvals_df[, "RMSEA>.05"], p_RMSEA_gt0.06 = pvals_df[, "RMSEA>.06"], p_RMSEA_gt0.08 = pvals_df[, "RMSEA>.08"], p_RMSEA_gt0.10 = pvals_df[, "RMSEA>.10"], stringsAsFactors = FALSE, row.names = NULL ) return(comparison_df) } # Search Function ---- # Wrapper for Regex-based search for variables in dataset; returns string vector #example:varlist(data,pattern="hsns",exclude="5|oj") #above example identifies all variable names that include 'hsns' but #excludes those that have a 5 or "oj". Pipe key "|" is the "or" operator #other operators (i.e,. "&") can also be used. varlist <- function (df=NULL, pattern=NULL, exclude=NULL, type=c("numeric","factor","character", "double", "logical", "integer"), ignore.case=TRUE) { vars <- character(0) if (any(type %in% c("numeric", "double"))) { vars <- c(vars,names(df)[sapply(df,is.double)]) #had is.numeric and separate is.double, but that was resulting in duplicates } if (any(type %in% "factor")) { vars <- c(vars,names(df)[sapply(df,is.factor)]) } if (any(type %in% "character")) { vars <- c(vars,names(df)[sapply(df,is.character)]) } if (any(type %in% "logical")) { vars <- c(vars,names(df)[sapply(df,is.logical)]) } if (any(type %in% "integer")) { vars <- c(vars,names(df)[sapply(df,is.integer)]) } if(!is.null(exclude)){ if(ignore.case==TRUE){ list<-vars[(!vars %in% vars[grepl(vars,pattern=exclude,ignore.case=TRUE)]) & grepl(vars,pattern=pattern, ignore.case=TRUE)] } else{ list<-vars[(!vars %in% vars[grepl(vars,pattern=exclude,ignore.case=FALSE)]) & grepl(vars,pattern=pattern, ignore.case=FALSE)] } } if(is.null(exclude)){ if(ignore.case==TRUE){ list<-vars[(!vars %in% exclude) & grepl(vars,pattern=pattern, ignore.case=TRUE)] } else{ list<-vars[(!vars %in% exclude) & grepl(vars,pattern=pattern, ignore.case=FALSE)] } } list } #More examples of how varlist function works, df = dataframe object: #can also look up data on the grepl function, all of that should work ## All variable starting with cred: # varlist(df,pattern="^cred") ### All variables labeled neo followed by any number (e.g., 'neo1','neo2', etc.): # varlist(df,pattern="neo[1-9]") ## All numeric variables: # varlist(df,type="numeric") ## All factor variable except variable gb and variables starting with c: # varlist(df,type="factor",exclude=c("gb|^c")) ## Can use in conjunction with sapply: # sapply(df[,varlist(df,type="numeric",pattern="credit")], summary) # Wrapper for identifying number of NAs per row ---- numNAs <- function(x) { sum(is.na(x)) } # Data cleaning ---- # Replace the NAs with 0 where where item 2 was denied df <- mutate(df, VACSS_3 = ifelse(VACSS_2 == 0, 0, VACSS_3), VACSS_4 = ifelse(VACSS_2 == 0, 0, VACSS_4), VACSS_5 = ifelse(VACSS_2 == 0, 0, VACSS_5), VACSS_5a = ifelse(VACSS_5 == 0, 0, VACSS_5a), VACSS_6a = ifelse(VACSS_6 == 0, 0, VACSS_6a)) follow_df <- mutate(follow_df, VACSS_3 = ifelse(VACSS_2 == 0, 0, VACSS_3), VACSS_4 = ifelse(VACSS_2 == 0, 0, VACSS_4), VACSS_5 = ifelse(VACSS_2 == 0, 0, VACSS_5), VACSS_5a = ifelse(VACSS_5 == 0, 0, VACSS_5a), VACSS_6a = ifelse(VACSS_6 == 0, 0, VACSS_6a)) # Convert follow-up data from long to wide form # Identifying variables that are constant within each subject ID constant_vars <- follow_df %>% group_by(SubjID) %>% summarise_all(~ length(unique(na.omit(.))) == 1) %>% summarise_all(all) constant_vars <- names(constant_vars)[constant_vars == T] follow_df <- follow_df %>% pivot_wider( id_cols = SubjID, names_from = studyweek, values_from = !c(constant_vars, studyweek), unused_fn = first ) # merge the dataframes df <- full_join(df, follow_df) ## Recode the C-SSRS items to remove local dependence ---- # make polytomous C-SSRS suicidal ideation item at each time point df <- df %>% mutate(VACSS_2r = ifelse(VACSS_5a == 1, 5, ifelse(VACSS_5 == 1, 4, ifelse(VACSS_4 == 1, 3, ifelse(VACSS_3 == 1, 2, ifelse(VACSS_2 == 1, 1, 0))))), VACSS_2r_4 = ifelse(VACSS_5a_4 == 1, 5, ifelse(VACSS_5_4 == 1, 4, ifelse(VACSS_4_4 == 1, 3, ifelse(VACSS_3_4 == 1, 2, ifelse(VACSS_2_4 == 1, 1, 0))))), VACSS_2r_8 = ifelse(VACSS_5a_8 == 1, 5, ifelse(VACSS_5_8 == 1, 4, ifelse(VACSS_4_8 == 1, 3, ifelse(VACSS_3_8 == 1, 2, ifelse(VACSS_2_8 == 1, 1, 0))))), VACSS_2r_12 = ifelse(VACSS_5a_12 == 1, 5, ifelse(VACSS_5_12 == 1, 4, ifelse(VACSS_4_12 == 1, 3, ifelse(VACSS_3_12 == 1, 2, ifelse(VACSS_2_12 == 1, 1, 0))))), VACSS_2r_18 = ifelse(VACSS_5a_18 == 1, 5, ifelse(VACSS_5_18 == 1, 4, ifelse(VACSS_4_18 == 1, 3, ifelse(VACSS_3_18 == 1, 2, ifelse(VACSS_2_18 == 1, 1, 0))))), VACSS_2r_24 = ifelse(VACSS_5a_24 == 1, 5, ifelse(VACSS_5_24 == 1, 4, ifelse(VACSS_4_24 == 1, 3, ifelse(VACSS_3_24 == 1, 2, ifelse(VACSS_2_24 == 1, 1, 0))))), VACSS_2r_36 = ifelse(VACSS_5a_36 == 1, 5, ifelse(VACSS_5_36 == 1, 4, ifelse(VACSS_4_36 == 1, 3, ifelse(VACSS_3_36 == 1, 2, ifelse(VACSS_2_36 == 1, 1, 0)))))) # Polytomous C-SSRS suicidal behavior item at each time point df <- df %>% mutate(VACSS_6r = ifelse(VACSS_6a == 1, 2, ifelse(VACSS_6 == 1, 1, 0)), VACSS_6r_4 = ifelse(VACSS_6a_4 == 1, 2, ifelse(VACSS_6_4 == 1, 1, 0)), VACSS_6r_8 = ifelse(VACSS_6a_8 == 1, 2, ifelse(VACSS_6_8 == 1, 1, 0)), VACSS_6r_12 = ifelse(VACSS_6a_12 == 1, 2, ifelse(VACSS_6_12 == 1, 1, 0)), VACSS_6r_18 = ifelse(VACSS_6a_18 == 1, 2, ifelse(VACSS_6_18 == 1, 1, 0)), VACSS_6r_24 = ifelse(VACSS_6a_24 == 1, 2, ifelse(VACSS_6_24 == 1, 1, 0)), VACSS_6r_36 = ifelse(VACSS_6a_36 == 1, 2, ifelse(VACSS_6_36 == 1, 1, 0))) # Remove participants with all C-SSRS items missing at all time points from baseline to wk 24 ---- vars <- select(df, matches("vacss"), -matches('r_|r$|36$')) %>% names df <- df[-which(apply(select(df, all_of(vars)), 1, numNAs) == length(vars)), ] # *********** ----- # Data analyses ---- # *********** ----- # Response Distribution ---- response_dist <- data.frame(matrix(ncol = 6, nrow = 8)) names(response_dist) <- c('wk1', 'wk4', 'wk8', 'wk12', 'wk18', 'wk24') row.names(response_dist) <- c('VACSS_1', 'VACSS_2', 'VACSS_3', 'VACSS_4', 'VACSS_5', 'VACSS_5a', 'VACSS_6', 'VACSS_6a') response_dist[,1] <- select(df, matches('VACSS_[1-6]$|VACSS_[1-6]a$')) %>% lapply((function(x){ n <- table(x)[[2]] p <- table(x) %>% prop.table() %>% .[[2]] paste0(n, " (", sprintf("%1.2f%%", 100*p), ")") })) %>% do.call(rbind, .) response_dist[,2] <- select(df, matches('VACSS_[1-6]_4$|VACSS_[1-6]a_4$')) %>% lapply((function(x){ n <- table(x)[[2]] p <- table(x) %>% prop.table() %>% .[[2]] paste0(n, " (", sprintf("%1.2f%%", 100*p), ")") })) %>% do.call(rbind, .) response_dist[,3] <- select(df, matches('VACSS_[1-6]_8$|VACSS_[1-6]a_8$')) %>% lapply((function(x){ n <- table(x)[[2]] p <- table(x) %>% prop.table() %>% .[[2]] paste0(n, " (", sprintf("%1.2f%%", 100*p), ")") })) %>% do.call(rbind, .) response_dist[,4] <- select(df, matches('VACSS_[1-6]_12$|VACSS_[1-6]a_12$')) %>% lapply((function(x){ n <- table(x)[[2]] p <- table(x) %>% prop.table() %>% .[[2]] paste0(n, " (", sprintf("%1.2f%%", 100*p), ")") })) %>% do.call(rbind, .) response_dist[,5] <- select(df, matches('VACSS_[1-6]_18$|VACSS_[1-6]a_18$')) %>% lapply((function(x){ n <- table(x)[[2]] p <- table(x) %>% prop.table() %>% .[[2]] paste0(n, " (", sprintf("%1.2f%%", 100*p), ")") })) %>% do.call(rbind, .) response_dist[,6] <- select(df, matches('VACSS_[1-6]_24$|VACSS_[1-6]a_24$')) %>% lapply((function(x){ n <- table(x)[[2]] p <- table(x) %>% prop.table() %>% .[[2]] paste0(n, " (", sprintf("%1.2f%%", 100*p), ")") })) %>% do.call(rbind, .) # df: wide format dataframe containing all timepoints and covariates # Full sample models ---- ## Baseline ---- long_cfa <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + VACSS_2r + VACSS_6r risk_4 =~ NA*VACSS_1_4 + VACSS_2r_4 + VACSS_6r_4 risk_8 =~ NA*VACSS_1_8 + VACSS_2r_8 + VACSS_6r_8 risk_12 =~ NA*VACSS_1_12 + VACSS_2r_12 + VACSS_6r_12 risk_18 =~ NA*VACSS_1_18 + VACSS_2r_18 + VACSS_6r_18 risk_24 =~ NA*VACSS_1_24 + VACSS_2r_24 + VACSS_6r_24 ## INTERCEPTS VACSS_1 ~ 0*1 VACSS_2r ~ 0*1 VACSS_6r ~ 0*1 VACSS_1_4 ~ 0*1 VACSS_2r_4 ~ 0*1 VACSS_6r_4 ~ 0*1 VACSS_1_8 ~ 0*1 VACSS_2r_8 ~ 0*1 VACSS_6r_8 ~ 0*1 VACSS_1_12 ~ 0*1 VACSS_2r_12 ~ 0*1 VACSS_6r_12 ~ 0*1 VACSS_1_18 ~ 0*1 VACSS_2r_18 ~ 0*1 VACSS_6r_18 ~ 0*1 VACSS_1_24 ~ 0*1 VACSS_2r_24 ~ 0*1 VACSS_6r_24 ~ 0*1 ## SCALING PARAMETER VACSS_1 ~*~ 1*VACSS_1 VACSS_2r ~*~ 1*VACSS_2r VACSS_6r ~*~ 1*VACSS_6r VACSS_1_4 ~*~ 1*VACSS_1_4 VACSS_1_8 ~*~ 1*VACSS_1_8 VACSS_1_12 ~*~ 1*VACSS_1_12 VACSS_1_18 ~*~ 1*VACSS_1_18 VACSS_1_24 ~*~ 1*VACSS_1_24 VACSS_2r_4 ~*~ 1*VACSS_2r_4 VACSS_6r_4 ~*~ 1*VACSS_6r_4 VACSS_2r_8 ~*~ 1*VACSS_2r_8 VACSS_6r_8 ~*~ 1*VACSS_6r_8 VACSS_2r_12 ~*~ 1*VACSS_2r_12 VACSS_6r_12 ~*~ 1*VACSS_6r_12 VACSS_2r_18 ~*~ 1*VACSS_2r_18 VACSS_6r_18 ~*~ 1*VACSS_6r_18 VACSS_2r_24 ~*~ 1*VACSS_2r_24 VACSS_6r_24 ~*~ 1*VACSS_6r_24 ## LV VARIANCE risk_1 ~~ 1*risk_1 risk_4 ~~ 1*risk_4 risk_8 ~~ 1*risk_8 risk_12 ~~ 1*risk_12 risk_18 ~~ 1*risk_18 risk_24 ~~ 1*risk_24 ## LATENT MEANS risk_1 ~ 0*1 risk_4 ~ 0*1 risk_8 ~ 0*1 risk_12 ~ 0*1 risk_18 ~ 0*1 risk_24 ~ 0*1 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_8 ~~ VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_18 + VACSS_1_24 VACSS_1_18 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_8 ~~ VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_18 + VACSS_2r_24 VACSS_2r_18 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_8 ~~ VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_18 + VACSS_6r_24 VACSS_6r_18 ~~ VACSS_6r_24 " fit0 <- cfa(long_cfa, data = df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') summary(fit0, fit.measures = T, standardized = T) ### Check convergent validity ---- factor_scores <- as.data.frame(lavPredict(fit0)) criteria <- select(df, CSSRS_VA_MOD_Suicide_Risk, PHQ9_9, PHQ_SCORE, ANX_SCORE, PCL5_SCORE) converge_cors <- cbind(factor_scores$risk_1, criteria) %>% cor(use = 'pairwise') ## Threshold Invariance ---- thresh_model <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + VACSS_2r + VACSS_6r risk_4 =~ NA*VACSS_1_4 + VACSS_2r_4 + VACSS_6r_4 risk_8 =~ NA*VACSS_1_8 + VACSS_2r_8 + VACSS_6r_8 risk_12 =~ NA*VACSS_1_12 + VACSS_2r_12 + VACSS_6r_12 risk_18 =~ NA*VACSS_1_18 + VACSS_2r_18 + VACSS_6r_18 risk_24 =~ NA*VACSS_1_24 + VACSS_2r_24 + VACSS_6r_24 ## THRESHOLDS: VACSS_1 | v1_thr*t1 VACSS_1_4 | v1_thr*t1 VACSS_1_8 | v1_thr*t1 VACSS_1_12 | v1_thr*t1 VACSS_1_18 | v1_thr*t1 VACSS_1_24 | v1_thr*t1 VACSS_2r | v2_thr1*t1 VACSS_2r | v2_thr2*t2 VACSS_2r | v2_thr3*t3 VACSS_2r | v2_thr4*t4 VACSS_2r | v2_thr5*t5 VACSS_2r_4 | v2_thr1*t1 VACSS_2r_4 | v2_thr2*t2 VACSS_2r_4 | v2_thr3*t3 VACSS_2r_4 | v2_thr4*t4 VACSS_2r_4 | v2_thr5*t5 VACSS_2r_8 | v2_thr1*t1 VACSS_2r_8 | v2_thr2*t2 VACSS_2r_8 | v2_thr3*t3 VACSS_2r_8 | v2_thr4*t4 VACSS_2r_8 | v2_thr5*t5 VACSS_2r_12 | v2_thr1*t1 VACSS_2r_12 | v2_thr2*t2 VACSS_2r_12 | v2_thr3*t3 VACSS_2r_12 | v2_thr4*t4 VACSS_2r_12 | v2_thr5*t5 VACSS_2r_18 | v2_thr1*t1 VACSS_2r_18 | v2_thr2*t2 VACSS_2r_18 | v2_thr3*t3 VACSS_2r_18 | v2_thr4*t4 VACSS_2r_18 | v2_thr5*t5 VACSS_2r_24 | v2_thr1*t1 VACSS_2r_24 | v2_thr2*t2 VACSS_2r_24 | v2_thr3*t3 VACSS_2r_24 | v2_thr4*t4 VACSS_2r_24 | v2_thr5*t5 VACSS_6r | v6_thr1*t1 VACSS_6r | v6_thr2*t2 VACSS_6r_4 | v6_thr1*t1 VACSS_6r_4 | v6_thr2*t2 VACSS_6r_8 | v6_thr1*t1 VACSS_6r_8 | v6_thr2*t2 VACSS_6r_12 | v6_thr1*t1 VACSS_6r_12 | v6_thr2*t2 VACSS_6r_18 | v6_thr1*t1 VACSS_6r_18 | v6_thr2*t2 VACSS_6r_24 | v6_thr1*t1 VACSS_6r_24 | v6_thr2*t2 ## INTERCEPTS VACSS_1 ~ 0*1 VACSS_2r ~ 0*1 VACSS_6r ~ 0*1 VACSS_1_4 ~ NA*1 VACSS_2r_4 ~ NA*1 VACSS_6r_4 ~ NA*1 VACSS_1_8 ~ NA*1 VACSS_2r_8 ~ NA*1 VACSS_6r_8 ~ NA*1 VACSS_1_12 ~ NA*1 VACSS_2r_12 ~ NA*1 VACSS_6r_12 ~ NA*1 VACSS_1_18 ~ NA*1 VACSS_2r_18 ~ NA*1 VACSS_6r_18 ~ NA*1 VACSS_1_24 ~ NA*1 VACSS_2r_24 ~ NA*1 VACSS_6r_24 ~ NA*1 ## SCALING PARAMETER VACSS_1 ~*~ 1*VACSS_1 VACSS_2r ~*~ 1*VACSS_2r VACSS_6r ~*~ 1*VACSS_6r # scaling factor of 1 item at each time point to 1 VACSS_1_4 ~*~ 1*VACSS_1_4 VACSS_1_8 ~*~ 1*VACSS_1_8 VACSS_1_12 ~*~ 1*VACSS_1_12 VACSS_1_18 ~*~ 1*VACSS_1_18 VACSS_1_24 ~*~ 1*VACSS_1_24 # free remaining scaling factors VACSS_2r_4 ~*~ NA*VACSS_2r_4 VACSS_6r_4 ~*~ NA*VACSS_6r_4 VACSS_2r_8 ~*~ NA*VACSS_2r_8 VACSS_6r_8 ~*~ NA*VACSS_6r_8 VACSS_2r_12 ~*~ NA*VACSS_2r_12 VACSS_6r_12 ~*~ NA*VACSS_6r_12 VACSS_2r_18 ~*~ NA*VACSS_2r_18 VACSS_6r_18 ~*~ NA*VACSS_6r_18 VACSS_2r_24 ~*~ NA*VACSS_2r_24 VACSS_6r_24 ~*~ NA*VACSS_6r_24 ## LATENT MEANS risk_1 ~ 0*1 risk_4 ~ 0*1 risk_8 ~ 0*1 risk_12 ~ 0*1 risk_18 ~ 0*1 risk_24 ~ 0*1 ## LV VARIANCE risk_1 ~~ 1*risk_1 risk_4 ~~ 1*risk_4 risk_8 ~~ 1*risk_8 risk_12 ~~ 1*risk_12 risk_18 ~~ 1*risk_18 risk_24 ~~ 1*risk_24 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_8 ~~ VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_18 + VACSS_1_24 VACSS_1_18 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_8 ~~ VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_18 + VACSS_2r_24 VACSS_2r_18 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_8 ~~ VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_18 + VACSS_6r_24 VACSS_6r_18 ~~ VACSS_6r_24 " fit1 <- cfa(thresh_model, data = df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') summary(fit1, fit.measures = T, standardized = T) ### Test Model net(fit0, fit1) #confirm the models are nested lavTestLRT(fit0, fit1) # chi-square test compare_fit(fit0, fit1) # CFI difference compare_nested_models_RMSEA(fit0, fit1) # RMSEAd ## Metric Invariance ---- # The following does not match the syntax you would get from measEq.syntax() # Proposition 7 in Wu & Estabrook (2016, p. 1029)) defines the following identification conditions for theshold and loading invariance: ## scaling parameter = 1 @ reference time point ## intercepts = 0 @ reference time point ## latent variable variances = 1 @ reference time point ## latent variable means = 0 @ all time points # This meets those standards. Note that it is not placing any constraints on the dichotomous item relative to the baseline model. # Wu & Estabrook point out that there is no way to test threshold only or threshold and loading invariance for dichotomous items. metric_model <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + l1*VACSS_1 + l2*VACSS_2r + l3*VACSS_6r risk_4 =~ NA*VACSS_1_4 + l1*VACSS_1_4 + l2*VACSS_2r_4 + l3*VACSS_6r_4 risk_8 =~ NA*VACSS_1_8 + l1*VACSS_1_8 + l2*VACSS_2r_8 + l3*VACSS_6r_8 risk_12 =~ NA*VACSS_1_12 + l1*VACSS_1_12 + l2*VACSS_2r_12 + l3*VACSS_6r_12 risk_18 =~ NA*VACSS_1_18 + l1*VACSS_1_18 + l2*VACSS_2r_18 + l3*VACSS_6r_18 risk_24 =~ NA*VACSS_1_24 + l1*VACSS_1_24 + l2*VACSS_2r_24 + l3*VACSS_6r_24 ## THRESHOLDS: VACSS_1 | v1_thr*t1 VACSS_1_4 | v1_thr*t1 VACSS_1_8 | v1_thr*t1 VACSS_1_12 | v1_thr*t1 VACSS_1_18 | v1_thr*t1 VACSS_1_24 | v1_thr*t1 VACSS_2r | v2_thr1*t1 VACSS_2r | v2_thr2*t2 VACSS_2r | v2_thr3*t3 VACSS_2r | v2_thr4*t4 VACSS_2r | v2_thr5*t5 VACSS_2r_4 | v2_thr1*t1 VACSS_2r_4 | v2_thr2*t2 VACSS_2r_4 | v2_thr3*t3 VACSS_2r_4 | v2_thr4*t4 VACSS_2r_4 | v2_thr5*t5 VACSS_2r_8 | v2_thr1*t1 VACSS_2r_8 | v2_thr2*t2 VACSS_2r_8 | v2_thr3*t3 VACSS_2r_8 | v2_thr4*t4 VACSS_2r_8 | v2_thr5*t5 VACSS_2r_12 | v2_thr1*t1 VACSS_2r_12 | v2_thr2*t2 VACSS_2r_12 | v2_thr3*t3 VACSS_2r_12 | v2_thr4*t4 VACSS_2r_12 | v2_thr5*t5 VACSS_2r_18 | v2_thr1*t1 VACSS_2r_18 | v2_thr2*t2 VACSS_2r_18 | v2_thr3*t3 VACSS_2r_18 | v2_thr4*t4 VACSS_2r_18 | v2_thr5*t5 VACSS_2r_24 | v2_thr1*t1 VACSS_2r_24 | v2_thr2*t2 VACSS_2r_24 | v2_thr3*t3 VACSS_2r_24 | v2_thr4*t4 VACSS_2r_24 | v2_thr5*t5 VACSS_6r | v6_thr1*t1 VACSS_6r | v6_thr2*t2 VACSS_6r_4 | v6_thr1*t1 VACSS_6r_4 | v6_thr2*t2 VACSS_6r_8 | v6_thr1*t1 VACSS_6r_8 | v6_thr2*t2 VACSS_6r_12 | v6_thr1*t1 VACSS_6r_12 | v6_thr2*t2 VACSS_6r_18 | v6_thr1*t1 VACSS_6r_18 | v6_thr2*t2 VACSS_6r_24 | v6_thr1*t1 VACSS_6r_24 | v6_thr2*t2 ## INTERCEPTS VACSS_1 ~ 0*1 VACSS_2r ~ 0*1 VACSS_6r ~ 0*1 VACSS_1_4 ~ NA*1 VACSS_2r_4 ~ NA*1 VACSS_6r_4 ~ NA*1 VACSS_1_8 ~ NA*1 VACSS_2r_8 ~ NA*1 VACSS_6r_8 ~ NA*1 VACSS_1_12 ~ NA*1 VACSS_2r_12 ~ NA*1 VACSS_6r_12 ~ NA*1 VACSS_1_18 ~ NA*1 VACSS_2r_18 ~ NA*1 VACSS_6r_18 ~ NA*1 VACSS_1_24 ~ NA*1 VACSS_2r_24 ~ NA*1 VACSS_6r_24 ~ NA*1 ## SCALING PARAMETER VACSS_1 ~*~ 1*VACSS_1 VACSS_2r ~*~ 1*VACSS_2r VACSS_6r ~*~ 1*VACSS_6r VACSS_1_4 ~*~ NA*VACSS_1_4 VACSS_1_8 ~*~ NA*VACSS_1_8 VACSS_1_12 ~*~ NA*VACSS_1_12 VACSS_1_18 ~*~ NA*VACSS_1_18 VACSS_1_24 ~*~ NA*VACSS_1_24 VACSS_2r_4 ~*~ NA*VACSS_2r_4 VACSS_6r_4 ~*~ NA*VACSS_6r_4 VACSS_2r_8 ~*~ NA*VACSS_2r_8 VACSS_6r_8 ~*~ NA*VACSS_6r_8 VACSS_2r_12 ~*~ NA*VACSS_2r_12 VACSS_6r_12 ~*~ NA*VACSS_6r_12 VACSS_2r_18 ~*~ NA*VACSS_2r_18 VACSS_6r_18 ~*~ NA*VACSS_6r_18 VACSS_2r_24 ~*~ NA*VACSS_2r_24 VACSS_6r_24 ~*~ NA*VACSS_6r_24 ## LATENT MEANS risk_1 ~ 0*1 risk_4 ~ 0*1 risk_8 ~ 0*1 risk_12 ~ 0*1 risk_18 ~ 0*1 risk_24 ~ 0*1 ## LV VARIANCE risk_1 ~~ 1*risk_1 risk_4 ~~ NA*risk_4 risk_8 ~~ NA*risk_8 risk_12 ~~ NA*risk_12 risk_18 ~~ NA*risk_18 risk_24 ~~ NA*risk_24 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_8 ~~ VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_18 + VACSS_1_24 VACSS_1_18 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_8 ~~ VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_18 + VACSS_2r_24 VACSS_2r_18 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_8 ~~ VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_18 + VACSS_6r_24 VACSS_6r_18 ~~ VACSS_6r_24 " fit2 <- cfa(metric_model, data = df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') summary(fit2, fit.measures = T, standardized = T) ## Test Model net(fit1, fit2) lavTestLRT(fit1, fit2) compare_fit(fit1, fit2) compare_nested_models_RMSEA(fit1, fit2) ## Scalar Invariance ---- scalar_model <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + l1*VACSS_1 + l2*VACSS_2r + l3*VACSS_6r risk_4 =~ NA*VACSS_1_4 + l1*VACSS_1_4 + l2*VACSS_2r_4 + l3*VACSS_6r_4 risk_8 =~ NA*VACSS_1_8 + l1*VACSS_1_8 + l2*VACSS_2r_8 + l3*VACSS_6r_8 risk_12 =~ NA*VACSS_1_12 + l1*VACSS_1_12 + l2*VACSS_2r_12 + l3*VACSS_6r_12 risk_18 =~ NA*VACSS_1_18 + l1*VACSS_1_18 + l2*VACSS_2r_18 + l3*VACSS_6r_18 risk_24 =~ NA*VACSS_1_24 + l1*VACSS_1_24 + l2*VACSS_2r_24 + l3*VACSS_6r_24 ## THRESHOLDS: VACSS_1 | v1_thr*t1 VACSS_1_4 | v1_thr*t1 VACSS_1_8 | v1_thr*t1 VACSS_1_12 | v1_thr*t1 VACSS_1_18 | v1_thr*t1 VACSS_1_24 | v1_thr*t1 VACSS_2r | v2_thr1*t1 VACSS_2r | v2_thr2*t2 VACSS_2r | v2_thr3*t3 VACSS_2r | v2_thr4*t4 VACSS_2r | v2_thr5*t5 VACSS_2r_4 | v2_thr1*t1 VACSS_2r_4 | v2_thr2*t2 VACSS_2r_4 | v2_thr3*t3 VACSS_2r_4 | v2_thr4*t4 VACSS_2r_4 | v2_thr5*t5 VACSS_2r_8 | v2_thr1*t1 VACSS_2r_8 | v2_thr2*t2 VACSS_2r_8 | v2_thr3*t3 VACSS_2r_8 | v2_thr4*t4 VACSS_2r_8 | v2_thr5*t5 VACSS_2r_12 | v2_thr1*t1 VACSS_2r_12 | v2_thr2*t2 VACSS_2r_12 | v2_thr3*t3 VACSS_2r_12 | v2_thr4*t4 VACSS_2r_12 | v2_thr5*t5 VACSS_2r_18 | v2_thr1*t1 VACSS_2r_18 | v2_thr2*t2 VACSS_2r_18 | v2_thr3*t3 VACSS_2r_18 | v2_thr4*t4 VACSS_2r_18 | v2_thr5*t5 VACSS_2r_24 | v2_thr1*t1 VACSS_2r_24 | v2_thr2*t2 VACSS_2r_24 | v2_thr3*t3 VACSS_2r_24 | v2_thr4*t4 VACSS_2r_24 | v2_thr5*t5 VACSS_6r | v6_thr1*t1 VACSS_6r | v6_thr2*t2 VACSS_6r_4 | v6_thr1*t1 VACSS_6r_4 | v6_thr2*t2 VACSS_6r_8 | v6_thr1*t1 VACSS_6r_8 | v6_thr2*t2 VACSS_6r_12 | v6_thr1*t1 VACSS_6r_12 | v6_thr2*t2 VACSS_6r_18 | v6_thr1*t1 VACSS_6r_18 | v6_thr2*t2 VACSS_6r_24 | v6_thr1*t1 VACSS_6r_24 | v6_thr2*t2 ## INTERCEPTS VACSS_1 ~ 0*1 VACSS_2r ~ 0*1 VACSS_6r ~ 0*1 VACSS_1_4 ~ 0*1 VACSS_2r_4 ~ 0*1 VACSS_6r_4 ~ 0*1 VACSS_1_8 ~ 0*1 VACSS_2r_8 ~ 0*1 VACSS_6r_8 ~ 0*1 VACSS_1_12 ~ 0*1 VACSS_2r_12 ~ 0*1 VACSS_6r_12 ~ 0*1 VACSS_1_18 ~ 0*1 VACSS_2r_18 ~ 0*1 VACSS_6r_18 ~ 0*1 VACSS_1_24 ~ 0*1 VACSS_2r_24 ~ 0*1 VACSS_6r_24 ~ 0*1 ## SCALING PARAMETER # Variances at T1 to 1 VACSS_1 ~*~ 1*VACSS_1 VACSS_2r ~*~ 1*VACSS_2r VACSS_6r ~*~ 1*VACSS_6r VACSS_1_4 ~*~ NA*VACSS_1_4 VACSS_1_8 ~*~ NA*VACSS_1_8 VACSS_1_12 ~*~ NA*VACSS_1_12 VACSS_1_18 ~*~ NA*VACSS_1_18 VACSS_1_24 ~*~ NA*VACSS_1_24 VACSS_2r_4 ~*~ NA*VACSS_2r_4 VACSS_6r_4 ~*~ NA*VACSS_6r_4 VACSS_2r_8 ~*~ NA*VACSS_2r_8 VACSS_6r_8 ~*~ NA*VACSS_6r_8 VACSS_2r_12 ~*~ NA*VACSS_2r_12 VACSS_6r_12 ~*~ NA*VACSS_6r_12 VACSS_2r_18 ~*~ NA*VACSS_2r_18 VACSS_6r_18 ~*~ NA*VACSS_6r_18 VACSS_2r_24 ~*~ NA*VACSS_2r_24 VACSS_6r_24 ~*~ NA*VACSS_6r_24 ## LATENT MEANS risk_1 ~ 0*1 risk_4 ~ 1 risk_8 ~ 1 risk_12 ~ 1 risk_18 ~ 1 risk_24 ~ 1 ## LV VARIANCE risk_1 ~~ 1*risk_1 risk_4 ~~ NA*risk_4 risk_8 ~~ NA*risk_8 risk_12 ~~ NA*risk_12 risk_18 ~~ NA*risk_18 risk_24 ~~ NA*risk_24 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_8 ~~ VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_18 + VACSS_1_24 VACSS_1_18 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_8 ~~ VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_18 + VACSS_2r_24 VACSS_2r_18 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_8 ~~ VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_18 + VACSS_6r_24 VACSS_6r_18 ~~ VACSS_6r_24 " fit3 <- cfa(scalar_model, data = df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') summary(fit3, fit.measures = T, standardized = T) ## Test Model net(fit2, fit3) lavTestLRT(fit2, fit3) compare_fit(fit2, fit3) compare_nested_models_RMSEA(fit2, fit3) modindices(fit3, op = '~1') %>% arrange(desc(`mi`)) # largest MI is VACSS_6r # VACSS_6r is the largest MI, so we will free that parameter ## Partial Scalar Invariance ---- scalar_model_v2 <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + l1*VACSS_1 + l2*VACSS_2r + l3*VACSS_6r risk_4 =~ NA*VACSS_1_4 + l1*VACSS_1_4 + l2*VACSS_2r_4 + l3*VACSS_6r_4 risk_8 =~ NA*VACSS_1_8 + l1*VACSS_1_8 + l2*VACSS_2r_8 + l3*VACSS_6r_8 risk_12 =~ NA*VACSS_1_12 + l1*VACSS_1_12 + l2*VACSS_2r_12 + l3*VACSS_6r_12 risk_18 =~ NA*VACSS_1_18 + l1*VACSS_1_18 + l2*VACSS_2r_18 + l3*VACSS_6r_18 risk_24 =~ NA*VACSS_1_24 + l1*VACSS_1_24 + l2*VACSS_2r_24 + l3*VACSS_6r_24 ## THRESHOLDS: VACSS_1 | v1_thr*t1 VACSS_1_4 | v1_thr*t1 VACSS_1_8 | v1_thr*t1 VACSS_1_12 | v1_thr*t1 VACSS_1_18 | v1_thr*t1 VACSS_1_24 | v1_thr*t1 VACSS_2r | v2_thr1*t1 VACSS_2r | v2_thr2*t2 VACSS_2r | v2_thr3*t3 VACSS_2r | v2_thr4*t4 VACSS_2r | v2_thr5*t5 VACSS_2r_4 | v2_thr1*t1 VACSS_2r_4 | v2_thr2*t2 VACSS_2r_4 | v2_thr3*t3 VACSS_2r_4 | v2_thr4*t4 VACSS_2r_4 | v2_thr5*t5 VACSS_2r_8 | v2_thr1*t1 VACSS_2r_8 | v2_thr2*t2 VACSS_2r_8 | v2_thr3*t3 VACSS_2r_8 | v2_thr4*t4 VACSS_2r_8 | v2_thr5*t5 VACSS_2r_12 | v2_thr1*t1 VACSS_2r_12 | v2_thr2*t2 VACSS_2r_12 | v2_thr3*t3 VACSS_2r_12 | v2_thr4*t4 VACSS_2r_12 | v2_thr5*t5 VACSS_2r_18 | v2_thr1*t1 VACSS_2r_18 | v2_thr2*t2 VACSS_2r_18 | v2_thr3*t3 VACSS_2r_18 | v2_thr4*t4 VACSS_2r_18 | v2_thr5*t5 VACSS_2r_24 | v2_thr1*t1 VACSS_2r_24 | v2_thr2*t2 VACSS_2r_24 | v2_thr3*t3 VACSS_2r_24 | v2_thr4*t4 VACSS_2r_24 | v2_thr5*t5 VACSS_6r | v6_thr1*t1 VACSS_6r | v6_thr2*t2 VACSS_6r_4 | v6_thr1*t1 VACSS_6r_4 | v6_thr2*t2 VACSS_6r_8 | v6_thr1*t1 VACSS_6r_8 | v6_thr2*t2 VACSS_6r_12 | v6_thr1*t1 VACSS_6r_12 | v6_thr2*t2 VACSS_6r_18 | v6_thr1*t1 VACSS_6r_18 | v6_thr2*t2 VACSS_6r_24 | v6_thr1*t1 VACSS_6r_24 | v6_thr2*t2 ## INTERCEPTS VACSS_1 ~ 0*1 VACSS_2r ~ 0*1 VACSS_6r ~ NA*1 # Freed this parameter VACSS_1_4 ~ 0*1 VACSS_2r_4 ~ 0*1 VACSS_6r_4 ~ 0*1 VACSS_1_8 ~ 0*1 VACSS_2r_8 ~ 0*1 VACSS_6r_8 ~ 0*1 VACSS_1_12 ~ 0*1 VACSS_2r_12 ~ 0*1 VACSS_6r_12 ~ 0*1 VACSS_1_18 ~ 0*1 VACSS_2r_18 ~ 0*1 VACSS_6r_18 ~ 0*1 VACSS_1_24 ~ 0*1 VACSS_2r_24 ~ 0*1 VACSS_6r_24 ~ 0*1 ## SCALING PARAMETER VACSS_1 ~*~ 1*VACSS_1 VACSS_2r ~*~ 1*VACSS_2r VACSS_6r ~*~ 1*VACSS_6r VACSS_1_4 ~*~ NA*VACSS_1_4 VACSS_1_8 ~*~ NA*VACSS_1_8 VACSS_1_12 ~*~ NA*VACSS_1_12 VACSS_1_18 ~*~ NA*VACSS_1_18 VACSS_1_24 ~*~ NA*VACSS_1_24 VACSS_2r_4 ~*~ NA*VACSS_2r_4 VACSS_6r_4 ~*~ NA*VACSS_6r_4 VACSS_2r_8 ~*~ NA*VACSS_2r_8 VACSS_6r_8 ~*~ NA*VACSS_6r_8 VACSS_2r_12 ~*~ NA*VACSS_2r_12 VACSS_6r_12 ~*~ NA*VACSS_6r_12 VACSS_2r_18 ~*~ NA*VACSS_2r_18 VACSS_6r_18 ~*~ NA*VACSS_6r_18 VACSS_2r_24 ~*~ NA*VACSS_2r_24 VACSS_6r_24 ~*~ NA*VACSS_6r_24 ## LATENT MEANS risk_1 ~ 0*1 risk_4 ~ 1 risk_8 ~ 1 risk_12 ~ 1 risk_18 ~ 1 risk_24 ~ 1 ## LV VARIANCE risk_1 ~~ 1*risk_1 risk_4 ~~ NA*risk_4 risk_8 ~~ NA*risk_8 risk_12 ~~ NA*risk_12 risk_18 ~~ NA*risk_18 risk_24 ~~ NA*risk_24 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_8 ~~ VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_18 + VACSS_1_24 VACSS_1_18 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_8 ~~ VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_18 + VACSS_2r_24 VACSS_2r_18 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_8 ~~ VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_18 + VACSS_6r_24 VACSS_6r_18 ~~ VACSS_6r_24 " fit3.1 <- cfa(scalar_model_v2, data = df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') summary(fit3.1, fit.measures = T, standardized = T) ## Test Model net(fit2, fit3.1) lavTestLRT(fit2, fit3.1) compare_fit(fit2, fit3.1) compare_nested_models_RMSEA(fit2, fit3.1) # Quantifying non-invariance ---- ## Pull out the relevant parameters ---- est <- parameterEstimates(fit3.1) # Loading for item 6 (constrained equal across time) lambda_6 <- est %>% filter(op == "=~", rhs == "VACSS_6r") %>% pull(est) # Thresholds for item 6 (constrained equal across time) tau_6 <- est %>% filter(op == "|", lhs == "VACSS_6r") %>% pull(est) # tau_6[1] = threshold between cat 0 and 1 (no behavior vs. lifetime) # tau_6[2] = threshold between cat 1 and 2 (lifetime vs. past 3 months) # Intercepts: freed at T1, fixed to 0 at follow-ups nu_6_T1 <- est %>% filter(op == "~1", lhs == "VACSS_6r") %>% pull(est) nu_6_FU <- 0 # fixed at follow-ups # Residual variance theta_T1 <- est %>% filter(op == "~~", lhs == "VACSS_6r", rhs == "VACSS_6r") %>% pull(est) theta_FU <- est %>% filter(op == "~~", lhs == rhs, lhs %in% c("VACSS_6r_4", "VACSS_6r_8", "VACSS_6r_12", "VACSS_6r_18", "VACSS_6r_24")) %>% pull(est) %>% mean # avg across follow-ups ## Define a function to compute P(Y >= c | eta) ---- prob_at_or_above_cond <- function(eta, lambda, nu, tau_c, theta) { 1 - pnorm((tau_c - lambda * eta - nu) / sqrt(theta)) } ## Choose meaningful percentiles of latent risk ---- # Use the empirical distribution of factor scores fscores <- lavPredict(fit3.1) # Use T1 factor scores as the reference distribution eta_T1 <- fscores[, "risk_1"] percentiles <- c(0.50, 0.75, 0.90, 0.95, 0.99) eta_vals <- quantile(eta_T1, percentiles, na.rm = TRUE) ## ---- Compute implied response probabilities ---- # P(any lifetime behavior) = P(Y >= 1) P_any_T1 <- sapply(eta_vals, function(e) prob_at_or_above_cond(e, lambda_6, nu_6_T1, tau_6[1], theta_T1)) P_any_FU <- sapply(eta_vals, function(e) prob_at_or_above_cond(e, lambda_6, nu_6_FU, tau_6[1], theta_FU)) P_any_FU - P_any_T1 # Subgroup Analyses ---- ## Sex ---- ### Males ---- male_df <- filter(df, female == 0) #### Baseline ---- fit0 <- cfa(long_cfa, data = male_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') summary(fit0, fit.measures = T, standardized = T) #### Threshold Invariance ---- fit1 <- cfa(thresh_model, data = male_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') ##### Test Model net(fit0, fit1) lavTestLRT(fit0, fit1) compare_nested_models_RMSEA(fit0, fit1) #### Metric Invariance ---- fit2 <- cfa(metric_model, data = male_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') ##### Test Model lavTestLRT(fit1, fit2) compare_nested_models_RMSEA(fit1, fit2) #### Scalar Invariance ---- fit3 <- cfa(scalar_model, data = male_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') summary(fit3, standardized = T) ##### Test Model lavTestLRT(fit2, fit3) compare_nested_models_RMSEA(fit2, fit3) modindices(fit3, op = '~1') %>% arrange(desc(`mi`)) # largest MI is VACSS_6r #### Partial Scalar Invariance ---- fit3.1 <- cfa(scalar_model_v2, data = male_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') summary(fit3.1, fit.measures = T, standardized = T) ##### Test Model lavTestLRT(fit2, fit3.1) compare_nested_models_RMSEA(fit2, fit3.1) ### Females ---- female_df <- filter(df, female == 1) #### Baseline ---- # excluding weeks 8 and 18 due to missing response options long_cfa_f <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + VACSS_2r + VACSS_6r risk_4 =~ NA*VACSS_1_4 + VACSS_2r_4 + VACSS_6r_4 risk_12 =~ NA*VACSS_1_12 + VACSS_2r_12 + VACSS_6r_12 risk_24 =~ NA*VACSS_1_24 + VACSS_2r_24 + VACSS_6r_24 ## INTERCEPTS VACSS_1 ~ 0*1 VACSS_2r ~ 0*1 VACSS_6r ~ 0*1 VACSS_1_4 ~ 0*1 VACSS_2r_4 ~ 0*1 VACSS_6r_4 ~ 0*1 VACSS_1_12 ~ 0*1 VACSS_2r_12 ~ 0*1 VACSS_6r_12 ~ 0*1 VACSS_1_24 ~ 0*1 VACSS_2r_24 ~ 0*1 VACSS_6r_24 ~ 0*1 ## SCALING PARAMETER VACSS_1 ~*~ 1*VACSS_1 VACSS_2r ~*~ 1*VACSS_2r VACSS_6r ~*~ 1*VACSS_6r VACSS_1_4 ~*~ 1*VACSS_1_4 VACSS_1_12 ~*~ 1*VACSS_1_12 VACSS_1_24 ~*~ 1*VACSS_1_24 VACSS_2r_4 ~*~ 1*VACSS_2r_4 VACSS_6r_4 ~*~ 1*VACSS_6r_4 VACSS_2r_12 ~*~ 1*VACSS_2r_12 VACSS_6r_12 ~*~ 1*VACSS_6r_12 VACSS_2r_24 ~*~ 1*VACSS_2r_24 VACSS_6r_24 ~*~ 1*VACSS_6r_24 ## LV VARIANCE risk_1 ~~ 1*risk_1 risk_4 ~~ 1*risk_4 risk_12 ~~ 1*risk_12 risk_24 ~~ 1*risk_24 ## LATENT MEANS risk_1 ~ 0*1 risk_4 ~ 0*1 risk_12 ~ 0*1 risk_24 ~ 0*1 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_12 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_12 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_12 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_12 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_12 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_12 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_24 " fit0_f <- cfa(long_cfa_f, data = female_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') summary(fit0_f, fit.measures = T, standardized = T) #### Threshold Invariance ---- thresh_model_f <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + VACSS_2r + VACSS_6r risk_4 =~ NA*VACSS_1_4 + VACSS_2r_4 + VACSS_6r_4 risk_12 =~ NA*VACSS_1_12 + VACSS_2r_12 + VACSS_6r_12 risk_24 =~ NA*VACSS_1_24 + VACSS_2r_24 + VACSS_6r_24 ## THRESHOLDS: VACSS_1 | v1_thr*t1 VACSS_1_4 | v1_thr*t1 VACSS_1_12 | v1_thr*t1 VACSS_1_24 | v1_thr*t1 VACSS_2r | v2_thr1*t1 VACSS_2r | v2_thr2*t2 VACSS_2r | v2_thr3*t3 VACSS_2r | v2_thr4*t4 VACSS_2r | v2_thr5*t5 VACSS_2r_4 | v2_thr1*t1 VACSS_2r_4 | v2_thr2*t2 VACSS_2r_4 | v2_thr3*t3 VACSS_2r_4 | v2_thr4*t4 VACSS_2r_4 | v2_thr5*t5 VACSS_2r_12 | v2_thr1*t1 VACSS_2r_12 | v2_thr2*t2 VACSS_2r_12 | v2_thr3*t3 VACSS_2r_12 | v2_thr4*t4 VACSS_2r_12 | v2_thr5*t5 VACSS_2r_24 | v2_thr1*t1 VACSS_2r_24 | v2_thr2*t2 VACSS_2r_24 | v2_thr3*t3 VACSS_2r_24 | v2_thr4*t4 VACSS_2r_24 | v2_thr5*t5 VACSS_6r | v6_thr1*t1 VACSS_6r | v6_thr2*t2 VACSS_6r_4 | v6_thr1*t1 VACSS_6r_4 | v6_thr2*t2 VACSS_6r_12 | v6_thr1*t1 VACSS_6r_12 | v6_thr2*t2 VACSS_6r_24 | v6_thr1*t1 VACSS_6r_24 | v6_thr2*t2 ## INTERCEPTS VACSS_1 ~ 0*1 VACSS_2r ~ 0*1 VACSS_6r ~ 0*1 VACSS_1_4 ~ NA*1 VACSS_2r_4 ~ NA*1 VACSS_6r_4 ~ NA*1 VACSS_1_12 ~ NA*1 VACSS_2r_12 ~ NA*1 VACSS_6r_12 ~ NA*1 VACSS_1_24 ~ NA*1 VACSS_2r_24 ~ NA*1 VACSS_6r_24 ~ NA*1 ## SCALING PARAMETER VACSS_1 ~*~ 1*VACSS_1 VACSS_2r ~*~ 1*VACSS_2r VACSS_6r ~*~ 1*VACSS_6r # scaling factor of 1 item at each time point to 1 VACSS_1_4 ~*~ 1*VACSS_1_4 VACSS_1_12 ~*~ 1*VACSS_1_12 VACSS_1_24 ~*~ 1*VACSS_1_24 # free remaining scaling factors VACSS_2r_4 ~*~ NA*VACSS_2r_4 VACSS_6r_4 ~*~ NA*VACSS_6r_4 VACSS_2r_12 ~*~ NA*VACSS_2r_12 VACSS_6r_12 ~*~ NA*VACSS_6r_12 VACSS_2r_24 ~*~ NA*VACSS_2r_24 VACSS_6r_24 ~*~ NA*VACSS_6r_24 ## LATENT MEANS risk_1 ~ 0*1 risk_4 ~ 0*1 risk_12 ~ 0*1 risk_24 ~ 0*1 ## LV VARIANCE risk_1 ~~ 1*risk_1 risk_4 ~~ 1*risk_4 risk_12 ~~ 1*risk_12 risk_24 ~~ 1*risk_24 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_12 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_12 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_12 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_12 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_12 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_12 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_24 " fit1_f <- cfa(thresh_model_f, data = female_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') summary(fit1_f, fit.measures = T, standardized = T) ##### Test Model net(fit0_f, fit1_f) lavTestLRT(fit0_f, fit1_f) compare_fit(fit0_f, fit1_f) compare_nested_models_RMSEA(fit0_f, fit1_f) #### Metric Invariance ---- metric_model_f <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + l1*VACSS_1 + l2*VACSS_2r + l3*VACSS_6r risk_4 =~ NA*VACSS_1_4 + l1*VACSS_1_4 + l2*VACSS_2r_4 + l3*VACSS_6r_4 risk_12 =~ NA*VACSS_1_12 + l1*VACSS_1_12 + l2*VACSS_2r_12 + l3*VACSS_6r_12 risk_24 =~ NA*VACSS_1_24 + l1*VACSS_1_24 + l2*VACSS_2r_24 + l3*VACSS_6r_24 ## THRESHOLDS: VACSS_1 | v1_thr*t1 VACSS_1_4 | v1_thr*t1 VACSS_1_12 | v1_thr*t1 VACSS_1_24 | v1_thr*t1 VACSS_2r | v2_thr1*t1 VACSS_2r | v2_thr2*t2 VACSS_2r | v2_thr3*t3 VACSS_2r | v2_thr4*t4 VACSS_2r | v2_thr5*t5 VACSS_2r_4 | v2_thr1*t1 VACSS_2r_4 | v2_thr2*t2 VACSS_2r_4 | v2_thr3*t3 VACSS_2r_4 | v2_thr4*t4 VACSS_2r_4 | v2_thr5*t5 VACSS_2r_12 | v2_thr1*t1 VACSS_2r_12 | v2_thr2*t2 VACSS_2r_12 | v2_thr3*t3 VACSS_2r_12 | v2_thr4*t4 VACSS_2r_12 | v2_thr5*t5 VACSS_2r_24 | v2_thr1*t1 VACSS_2r_24 | v2_thr2*t2 VACSS_2r_24 | v2_thr3*t3 VACSS_2r_24 | v2_thr4*t4 VACSS_2r_24 | v2_thr5*t5 VACSS_6r | v6_thr1*t1 VACSS_6r | v6_thr2*t2 VACSS_6r_4 | v6_thr1*t1 VACSS_6r_4 | v6_thr2*t2 VACSS_6r_12 | v6_thr1*t1 VACSS_6r_12 | v6_thr2*t2 VACSS_6r_24 | v6_thr1*t1 VACSS_6r_24 | v6_thr2*t2 ## INTERCEPTS VACSS_1 ~ 0*1 VACSS_2r ~ 0*1 VACSS_6r ~ 0*1 VACSS_1_4 ~ NA*1 VACSS_2r_4 ~ NA*1 VACSS_6r_4 ~ NA*1 VACSS_1_12 ~ NA*1 VACSS_2r_12 ~ NA*1 VACSS_6r_12 ~ NA*1 VACSS_1_24 ~ NA*1 VACSS_2r_24 ~ NA*1 VACSS_6r_24 ~ NA*1 ## SCALING PARAMETER VACSS_1 ~*~ 1*VACSS_1 VACSS_2r ~*~ 1*VACSS_2r VACSS_6r ~*~ 1*VACSS_6r VACSS_1_4 ~*~ NA*VACSS_1_4 VACSS_1_12 ~*~ NA*VACSS_1_12 VACSS_1_24 ~*~ NA*VACSS_1_24 VACSS_2r_4 ~*~ NA*VACSS_2r_4 VACSS_6r_4 ~*~ NA*VACSS_6r_4 VACSS_2r_12 ~*~ NA*VACSS_2r_12 VACSS_6r_12 ~*~ NA*VACSS_6r_12 VACSS_2r_24 ~*~ NA*VACSS_2r_24 VACSS_6r_24 ~*~ NA*VACSS_6r_24 ## LATENT MEANS risk_1 ~ 0*1 risk_4 ~ 0*1 risk_12 ~ 0*1 risk_24 ~ 0*1 ## LV VARIANCE risk_1 ~~ 1*risk_1 risk_4 ~~ NA*risk_4 risk_12 ~~ NA*risk_12 risk_24 ~~ NA*risk_24 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_12 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_12 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_12 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_12 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_12 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_12 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_24 " fit2_f <- cfa(metric_model_f, data = female_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') summary(fit2_f, fit.measures = T, standardized = T) ##### Test Model net(fit1_f, fit2_f) lavTestLRT(fit1_f, fit2_f) compare_fit(fit1_f, fit2_f) compare_nested_models_RMSEA(fit1_f, fit2_f) #### Scalar Invariance ---- scalar_model_f <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + l1*VACSS_1 + l2*VACSS_2r + l3*VACSS_6r risk_4 =~ NA*VACSS_1_4 + l1*VACSS_1_4 + l2*VACSS_2r_4 + l3*VACSS_6r_4 risk_12 =~ NA*VACSS_1_12 + l1*VACSS_1_12 + l2*VACSS_2r_12 + l3*VACSS_6r_12 risk_24 =~ NA*VACSS_1_24 + l1*VACSS_1_24 + l2*VACSS_2r_24 + l3*VACSS_6r_24 ## THRESHOLDS: VACSS_1 | v1_thr*t1 VACSS_1_4 | v1_thr*t1 VACSS_1_12 | v1_thr*t1 VACSS_1_24 | v1_thr*t1 VACSS_2r | v2_thr1*t1 VACSS_2r | v2_thr2*t2 VACSS_2r | v2_thr3*t3 VACSS_2r | v2_thr4*t4 VACSS_2r | v2_thr5*t5 VACSS_2r_4 | v2_thr1*t1 VACSS_2r_4 | v2_thr2*t2 VACSS_2r_4 | v2_thr3*t3 VACSS_2r_4 | v2_thr4*t4 VACSS_2r_4 | v2_thr5*t5 VACSS_2r_12 | v2_thr1*t1 VACSS_2r_12 | v2_thr2*t2 VACSS_2r_12 | v2_thr3*t3 VACSS_2r_12 | v2_thr4*t4 VACSS_2r_12 | v2_thr5*t5 VACSS_2r_24 | v2_thr1*t1 VACSS_2r_24 | v2_thr2*t2 VACSS_2r_24 | v2_thr3*t3 VACSS_2r_24 | v2_thr4*t4 VACSS_2r_24 | v2_thr5*t5 VACSS_6r | v6_thr1*t1 VACSS_6r | v6_thr2*t2 VACSS_6r_4 | v6_thr1*t1 VACSS_6r_4 | v6_thr2*t2 VACSS_6r_12 | v6_thr1*t1 VACSS_6r_12 | v6_thr2*t2 VACSS_6r_24 | v6_thr1*t1 VACSS_6r_24 | v6_thr2*t2 ## INTERCEPTS VACSS_1 ~ 0*1 VACSS_2r ~ 0*1 VACSS_6r ~ 0*1 VACSS_1_4 ~ 0*1 VACSS_2r_4 ~ 0*1 VACSS_6r_4 ~ 0*1 VACSS_1_12 ~ 0*1 VACSS_2r_12 ~ 0*1 VACSS_6r_12 ~ 0*1 VACSS_1_24 ~ 0*1 VACSS_2r_24 ~ 0*1 VACSS_6r_24 ~ 0*1 ## SCALING PARAMETER VACSS_1 ~*~ 1*VACSS_1 VACSS_2r ~*~ 1*VACSS_2r VACSS_6r ~*~ 1*VACSS_6r VACSS_1_4 ~*~ NA*VACSS_1_4 VACSS_1_12 ~*~ NA*VACSS_1_12 VACSS_1_24 ~*~ NA*VACSS_1_24 VACSS_2r_4 ~*~ NA*VACSS_2r_4 VACSS_6r_4 ~*~ NA*VACSS_6r_4 VACSS_2r_12 ~*~ NA*VACSS_2r_12 VACSS_6r_12 ~*~ NA*VACSS_6r_12 VACSS_2r_24 ~*~ NA*VACSS_2r_24 VACSS_6r_24 ~*~ NA*VACSS_6r_24 ## LATENT MEANS risk_1 ~ 0*1 risk_4 ~ 1 risk_12 ~ 1 risk_24 ~ 1 ## LV VARIANCE risk_1 ~~ 1*risk_1 risk_4 ~~ NA*risk_4 risk_12 ~~ NA*risk_12 risk_24 ~~ NA*risk_24 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_12 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_12 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_12 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_12 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_12 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_12 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_24 " fit3_f <- cfa(scalar_model_f, data = female_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') summary(fit3_f, fit.measures = T, standardized = T) ##### Test Model net(fit2_f, fit3_f) lavTestLRT(fit2_f, fit3_f) compare_fit(fit2_f, fit3_f) compare_nested_models_RMSEA(fit2_f, fit3_f) modindices(fit3_f, op = '~1') %>% arrange(desc(`mi`)) # largest MI is VACSS_6r #### Partial Scalar Invariance ---- scalar_model_v2_f <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + l1*VACSS_1 + l2*VACSS_2r + l3*VACSS_6r risk_4 =~ NA*VACSS_1_4 + l1*VACSS_1_4 + l2*VACSS_2r_4 + l3*VACSS_6r_4 risk_12 =~ NA*VACSS_1_12 + l1*VACSS_1_12 + l2*VACSS_2r_12 + l3*VACSS_6r_12 risk_24 =~ NA*VACSS_1_24 + l1*VACSS_1_24 + l2*VACSS_2r_24 + l3*VACSS_6r_24 ## THRESHOLDS: VACSS_1 | v1_thr*t1 VACSS_1_4 | v1_thr*t1 VACSS_1_12 | v1_thr*t1 VACSS_1_24 | v1_thr*t1 VACSS_2r | v2_thr1*t1 VACSS_2r | v2_thr2*t2 VACSS_2r | v2_thr3*t3 VACSS_2r | v2_thr4*t4 VACSS_2r | v2_thr5*t5 VACSS_2r_4 | v2_thr1*t1 VACSS_2r_4 | v2_thr2*t2 VACSS_2r_4 | v2_thr3*t3 VACSS_2r_4 | v2_thr4*t4 VACSS_2r_4 | v2_thr5*t5 VACSS_2r_12 | v2_thr1*t1 VACSS_2r_12 | v2_thr2*t2 VACSS_2r_12 | v2_thr3*t3 VACSS_2r_12 | v2_thr4*t4 VACSS_2r_12 | v2_thr5*t5 VACSS_2r_24 | v2_thr1*t1 VACSS_2r_24 | v2_thr2*t2 VACSS_2r_24 | v2_thr3*t3 VACSS_2r_24 | v2_thr4*t4 VACSS_2r_24 | v2_thr5*t5 VACSS_6r | v6_thr1*t1 VACSS_6r | v6_thr2*t2 VACSS_6r_4 | v6_thr1*t1 VACSS_6r_4 | v6_thr2*t2 VACSS_6r_12 | v6_thr1*t1 VACSS_6r_12 | v6_thr2*t2 VACSS_6r_24 | v6_thr1*t1 VACSS_6r_24 | v6_thr2*t2 ## INTERCEPTS VACSS_1 ~ 0*1 VACSS_2r ~ 0*1 VACSS_6r ~ NA*1 # freed intercept VACSS_1_4 ~ 0*1 VACSS_2r_4 ~ 0*1 VACSS_6r_4 ~ 0*1 VACSS_1_12 ~ 0*1 VACSS_2r_12 ~ 0*1 VACSS_6r_12 ~ 0*1 VACSS_1_24 ~ 0*1 VACSS_2r_24 ~ 0*1 VACSS_6r_24 ~ 0*1 ## SCALING PARAMETER VACSS_1 ~*~ 1*VACSS_1 VACSS_2r ~*~ 1*VACSS_2r VACSS_6r ~*~ 1*VACSS_6r VACSS_1_4 ~*~ NA*VACSS_1_4 VACSS_1_12 ~*~ NA*VACSS_1_12 VACSS_1_24 ~*~ NA*VACSS_1_24 VACSS_2r_4 ~*~ NA*VACSS_2r_4 VACSS_6r_4 ~*~ NA*VACSS_6r_4 VACSS_2r_12 ~*~ NA*VACSS_2r_12 VACSS_6r_12 ~*~ NA*VACSS_6r_12 VACSS_2r_24 ~*~ NA*VACSS_2r_24 VACSS_6r_24 ~*~ NA*VACSS_6r_24 ## LATENT MEANS risk_1 ~ 0*1 risk_4 ~ 1 risk_12 ~ 1 risk_24 ~ 1 ## LV VARIANCE risk_1 ~~ 1*risk_1 risk_4 ~~ NA*risk_4 risk_12 ~~ NA*risk_12 risk_24 ~~ NA*risk_24 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_12 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_12 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_12 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_12 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_12 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_12 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_24 " fit3.1_f <- cfa(scalar_model_v2_f, data = female_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') summary(fit3.1_f, fit.measures = T, standardized = T) ##### Test Model lavTestLRT(fit2_f, fit3.1_f) compare_fit(fit2_f, fit3.1_f) compare_nested_models_RMSEA(fit2_f, fit3.1_f) ## Race ---- ### White ---- white_df <- filter(df, Race == 1) #### Baseline ---- fit0 <- cfa(long_cfa, data = white_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') summary(fit0, fit.measures = T, standardized = T) #### Threshold Invariance ---- fit1 <- cfa(thresh_model, data = white_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') ##### Test Model net(fit0, fit1) lavTestLRT(fit0, fit1) compare_nested_models_RMSEA(fit0, fit1) #### Metric Invariance ---- fit2 <- cfa(metric_model, data = white_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') ##### Test Model lavTestLRT(fit1, fit2) compare_nested_models_RMSEA(fit1, fit2) #### Scalar Invariance ---- fit3 <- cfa(scalar_model, data = white_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') ##### Test Model lavTestLRT(fit2, fit3) compare_nested_models_RMSEA(fit2, fit3) modindices(fit3, op = '~1') %>% arrange(desc(`mi`)) # largest MI is VACSS_6r #### Partial Scalar Invariance ---- fit3.1 <- cfa(scalar_model_v2, data = white_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') summary(fit3.1, fit.measures = T, standardized = T) ##### Test Model lavTestLRT(fit2, fit3.1) compare_nested_models_RMSEA(fit2, fit3.1) ### BIPOC ---- # Note that `bipoc_df` is shorthand for all individuals self-identifying their race as non-White bipoc_df <- filter(df, Race != 1 & Race != 6) #### Baseline ---- fit0 <- cfa(long_cfa, data = bipoc_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') summary(fit0, fit.measures = T, standardized = T) #### Threshold Invariance ---- fit1 <- cfa(thresh_model, data = bipoc_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') ##### Test Model net(fit0, fit1) lavTestLRT(fit0, fit1) compare_nested_models_RMSEA(fit0, fit1) #### Metric Invariance ---- fit2 <- cfa(metric_model, data = bipoc_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') ##### Test Model lavTestLRT(fit1, fit2) compare_nested_models_RMSEA(fit1, fit2) #### Scalar Invariance ---- fit3 <- cfa(scalar_model, data = bipoc_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') ##### Test Model lavTestLRT(fit2, fit3) compare_nested_models_RMSEA(fit2, fit3) modindices(fit3, op = '~1') %>% arrange(desc(`mi`)) # largest MI is VACSS_6r #### Partial Scalar Invariance ---- fit3.1 <- cfa(scalar_model_v2, data = bipoc_df, ordered = T, estimator = "WLSMV", missing = 'pairwise', parameterization = 'delta') summary(fit3.1, fit.measures = T, standardized = T) ##### Test Model lavTestLRT(fit2, fit3.1) compare_nested_models_RMSEA(fit2, fit3.1) # *********** ----- # Supplemental Analyses ---- # *********** ----- # Distribution of Recoded Scores ---- # Define the recoded variables grouped by item and wave recoded_vars <- list( VACSS_1 = c("VACSS_1", "VACSS_1_4", "VACSS_1_8", "VACSS_1_12", "VACSS_1_18", "VACSS_1_24"), VACSS_2r = c("VACSS_2r", "VACSS_2r_4", "VACSS_2r_8", "VACSS_2r_12", "VACSS_2r_18", "VACSS_2r_24"), VACSS_6r = c("VACSS_6r", "VACSS_6r_4", "VACSS_6r_8", "VACSS_6r_12", "VACSS_6r_18", "VACSS_6r_24") ) wave_labels <- c("Baseline", "Week 4", "Week 8", "Week 12", "Week 18", "Week 24") # Response category labels for each recoded item response_labels <- list( VACSS_1 = c( "0" = "0 (No wish to be dead)", "1" = "1 (Wish to be dead)" ), VACSS_2r = c( "0" = "0 (No active SI)", "1" = "1 (Non-specific active SI)", "2" = "2 (Active SI with method, no plan)", "3" = "3 (Active SI with some intent, no plan)", "4" = "4 (Active SI with specific plan)", "5" = "5 (Active SI with specific plan and intent)" ), VACSS_6r = c( "0" = "0 (No history of suicidal behavior)", "1" = "1 (Lifetime history of suicidal behavior)", "2" = "2 (Suicidal behavior in past 3 months)" ) ) # ---- Helper: build a per-item frequency table get_freq_table <- function(item_vars, wave_labels, data, response_lookup) { all_levels <- names(response_lookup) result_list <- lapply(seq_along(item_vars), function(i) { v <- item_vars[i] x <- as.character(data[[v]]) x <- x[!is.na(x)] n_total <- length(x) counts <- sapply(all_levels, function(lvl) sum(x == lvl)) pcts <- round(counts / n_total * 100, 2) formatted <- sprintf("%d (%.2f%%)", counts, pcts) data.frame( Response = response_lookup[all_levels], Wave = wave_labels[i], formatted = formatted, stringsAsFactors = FALSE ) }) combined <- do.call(rbind, result_list) out <- combined %>% pivot_wider(names_from = Wave, values_from = formatted) # Preserve wave order out <- out[, c("Response", wave_labels)] out } # Generate tables table_VACSS_1 <- get_freq_table(recoded_vars$VACSS_1, wave_labels, df, response_labels$VACSS_1) table_VACSS_2r <- get_freq_table(recoded_vars$VACSS_2r, wave_labels, df, response_labels$VACSS_2r) table_VACSS_6r <- get_freq_table(recoded_vars$VACSS_6r, wave_labels, df, response_labels$VACSS_6r) combined_table <- bind_rows( table_VACSS_1, table_VACSS_2r, table_VACSS_6r ) # Group x Time models ---- ## BIPOC x Time model ---- ### Baseline ---- long_cfa_g <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + VACSS_2r + VACSS_6r risk_4 =~ NA*VACSS_1_4 + VACSS_2r_4 + VACSS_6r_4 risk_8 =~ NA*VACSS_1_8 + VACSS_2r_8 + VACSS_6r_8 risk_12 =~ NA*VACSS_1_12 + VACSS_2r_12 + VACSS_6r_12 risk_18 =~ NA*VACSS_1_18 + VACSS_2r_18 + VACSS_6r_18 risk_24 =~ NA*VACSS_1_24 + VACSS_2r_24 + VACSS_6r_24 ## INTERCEPTS # all set to zero VACSS_1 ~ 0*1 VACSS_2r ~ 0*1 VACSS_6r ~ 0*1 VACSS_1_4 ~ 0*1 VACSS_2r_4 ~ 0*1 VACSS_6r_4 ~ 0*1 VACSS_1_8 ~ 0*1 VACSS_2r_8 ~ 0*1 VACSS_6r_8 ~ 0*1 VACSS_1_12 ~ 0*1 VACSS_2r_12 ~ 0*1 VACSS_6r_12 ~ 0*1 VACSS_1_18 ~ 0*1 VACSS_2r_18 ~ 0*1 VACSS_6r_18 ~ 0*1 VACSS_1_24 ~ 0*1 VACSS_2r_24 ~ 0*1 VACSS_6r_24 ~ 0*1 ## SCALING PARAMETER # Scaling at T1 to 1 VACSS_1 ~*~ 1*VACSS_1 VACSS_2r ~*~ 1*VACSS_2r VACSS_6r ~*~ 1*VACSS_6r VACSS_1_4 ~*~ 1*VACSS_1_4 VACSS_1_8 ~*~ 1*VACSS_1_8 VACSS_1_12 ~*~ 1*VACSS_1_12 VACSS_1_18 ~*~ 1*VACSS_1_18 VACSS_1_24 ~*~ 1*VACSS_1_24 VACSS_2r_4 ~*~ 1*VACSS_2r_4 VACSS_6r_4 ~*~ 1*VACSS_6r_4 VACSS_2r_8 ~*~ 1*VACSS_2r_8 VACSS_6r_8 ~*~ 1*VACSS_6r_8 VACSS_2r_12 ~*~ 1*VACSS_2r_12 VACSS_6r_12 ~*~ 1*VACSS_6r_12 VACSS_2r_18 ~*~ 1*VACSS_2r_18 VACSS_6r_18 ~*~ 1*VACSS_6r_18 VACSS_2r_24 ~*~ 1*VACSS_2r_24 VACSS_6r_24 ~*~ 1*VACSS_6r_24 ## LV VARIANCE risk_1 ~~ 1*risk_1 risk_4 ~~ 1*risk_4 risk_8 ~~ 1*risk_8 risk_12 ~~ 1*risk_12 risk_18 ~~ 1*risk_18 risk_24 ~~ 1*risk_24 ## LATENT MEANS ## At the reference measurement occassion, the common factor mean is constrained to zero risk_1 ~ 0*1 risk_4 ~ 0*1 risk_8 ~ 0*1 risk_12 ~ 0*1 risk_18 ~ 0*1 risk_24 ~ 0*1 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_8 ~~ VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_18 + VACSS_1_24 VACSS_1_18 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_8 ~~ VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_18 + VACSS_2r_24 VACSS_2r_18 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_8 ~~ VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_18 + VACSS_6r_24 VACSS_6r_18 ~~ VACSS_6r_24 " fit0_bipoc <- cfa(long_cfa_g, data = df, ordered = T, estimator = "WLSMV", group = "bipoc", missing = 'pairwise', parameterization = 'delta') sum_fit0_bipoc <- summary(fit0_bipoc, fit.measures = T, standardized = T) ### Theshold Invariance ---- thresh_model_g <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + VACSS_2r + VACSS_6r risk_4 =~ NA*VACSS_1_4 + VACSS_2r_4 + VACSS_6r_4 risk_8 =~ NA*VACSS_1_8 + VACSS_2r_8 + VACSS_6r_8 risk_12 =~ NA*VACSS_1_12 + VACSS_2r_12 + VACSS_6r_12 risk_18 =~ NA*VACSS_1_18 + VACSS_2r_18 + VACSS_6r_18 risk_24 =~ NA*VACSS_1_24 + VACSS_2r_24 + VACSS_6r_24 ## THRESHOLDS: VACSS_1 | c(v1_thr, v1_thr)*t1 VACSS_1_4 | c(v1_thr, v1_thr)*t1 VACSS_1_8 | c(v1_thr, v1_thr)*t1 VACSS_1_12 | c(v1_thr, v1_thr)*t1 VACSS_1_18 | c(v1_thr, v1_thr)*t1 VACSS_1_24 | c(v1_thr, v1_thr)*t1 VACSS_2r | c(v2_thr1, v2_thr1)*t1 VACSS_2r | c(v2_thr2, v2_thr2)*t2 VACSS_2r | c(v2_thr3, v2_thr3)*t3 VACSS_2r | c(v2_thr4, v2_thr4)*t4 VACSS_2r | c(v2_thr5, v2_thr5)*t5 VACSS_2r_4 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_4 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_4 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_4 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_4 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_8 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_8 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_8 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_8 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_8 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_12 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_12 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_12 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_12 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_12 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_18 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_18 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_18 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_18 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_18 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_24 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_24 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_24 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_24 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_24 | c(v2_thr5, v2_thr5)*t5 VACSS_6r | c(v6_thr1, v6_thr1)*t1 VACSS_6r | c(v6_thr2, v6_thr2)*t2 VACSS_6r_4 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_4 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_8 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_8 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_12 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_12 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_18 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_18 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_24 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_24 | c(v6_thr2, v6_thr2)*t2 ## INTERCEPTS VACSS_1 ~ c(0, NA)*1 VACSS_2r ~ c(0, NA)*1 VACSS_6r ~ c(0, NA)*1 VACSS_1_4 ~ c(NA, NA)*1 VACSS_2r_4 ~ c(NA, NA)*1 VACSS_6r_4 ~ c(NA, NA)*1 VACSS_1_8 ~ c(NA, NA)*1 VACSS_2r_8 ~ c(NA, NA)*1 VACSS_6r_8 ~ c(NA, NA)*1 VACSS_1_12 ~ c(NA, NA)*1 VACSS_2r_12 ~ c(NA, NA)*1 VACSS_6r_12 ~ c(NA, NA)*1 VACSS_1_18 ~ c(NA, NA)*1 VACSS_2r_18 ~ c(NA, NA)*1 VACSS_6r_18 ~ c(NA, NA)*1 VACSS_1_24 ~ c(NA, NA)*1 VACSS_2r_24 ~ c(NA, NA)*1 VACSS_6r_24 ~ c(NA, NA)*1 ## SCALING PARAMETER VACSS_1 ~*~ c(1, 1)*VACSS_1 VACSS_2r ~*~ c(1, NA)*VACSS_2r VACSS_6r ~*~ c(1, NA)*VACSS_6r # scaling factor of 1 item at each time point to 1 VACSS_1_4 ~*~ c(1, 1)*VACSS_1_4 VACSS_1_8 ~*~ c(1, 1)*VACSS_1_8 VACSS_1_12 ~*~ c(1, 1)*VACSS_1_12 VACSS_1_18 ~*~ c(1, 1)*VACSS_1_18 VACSS_1_24 ~*~ c(1, 1)*VACSS_1_24 # free remaining scaling factors VACSS_2r_4 ~*~ c(NA, NA)*VACSS_2r_4 VACSS_6r_4 ~*~ c(NA, NA)*VACSS_6r_4 VACSS_2r_8 ~*~ c(NA, NA)*VACSS_2r_8 VACSS_6r_8 ~*~ c(NA, NA)*VACSS_6r_8 VACSS_2r_12 ~*~ c(NA, NA)*VACSS_2r_12 VACSS_6r_12 ~*~ c(NA, NA)*VACSS_6r_12 VACSS_2r_18 ~*~ c(NA, NA)*VACSS_2r_18 VACSS_6r_18 ~*~ c(NA, NA)*VACSS_6r_18 VACSS_2r_24 ~*~ c(NA, NA)*VACSS_2r_24 VACSS_6r_24 ~*~ c(NA, NA)*VACSS_6r_24 ## LATENT MEANS risk_1 ~ 0*1 risk_4 ~ 0*1 risk_8 ~ 0*1 risk_12 ~ 0*1 risk_18 ~ 0*1 risk_24 ~ 0*1 ## LV VARIANCE risk_1 ~~ 1*risk_1 risk_4 ~~ 1*risk_4 risk_8 ~~ 1*risk_8 risk_12 ~~ 1*risk_12 risk_18 ~~ 1*risk_18 risk_24 ~~ 1*risk_24 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_8 ~~ VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_18 + VACSS_1_24 VACSS_1_18 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_8 ~~ VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_18 + VACSS_2r_24 VACSS_2r_18 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_8 ~~ VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_18 + VACSS_6r_24 VACSS_6r_18 ~~ VACSS_6r_24 " fit1_bipoc <- cfa(thresh_model_g, data = df, ordered = T, estimator = "WLSMV", group = "bipoc", missing = 'pairwise', parameterization = 'delta') sum_fit1_bipoc <- summary(fit1_bipoc, fit.measures = T, standardized = T) ##### Test Model net(fit0_bipoc, fit1_bipoc) fit0_fit1_lrt_bipoc <- lavTestLRT(fit0_bipoc, fit1_bipoc) fit0_fit1_cfi_bipoc <- compare_fit(fit0_bipoc, fit1_bipoc) fit0_fit1_rmsea_bipoc <- compare_nested_models_RMSEA(fit0_bipoc, fit1_bipoc) ### Metric Invariance ---- metric_model_g <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + l1*VACSS_1 + l2*VACSS_2r + l3*VACSS_6r risk_4 =~ NA*VACSS_1_4 + l1*VACSS_1_4 + l2*VACSS_2r_4 + l3*VACSS_6r_4 risk_8 =~ NA*VACSS_1_8 + l1*VACSS_1_8 + l2*VACSS_2r_8 + l3*VACSS_6r_8 risk_12 =~ NA*VACSS_1_12 + l1*VACSS_1_12 + l2*VACSS_2r_12 + l3*VACSS_6r_12 risk_18 =~ NA*VACSS_1_18 + l1*VACSS_1_18 + l2*VACSS_2r_18 + l3*VACSS_6r_18 risk_24 =~ NA*VACSS_1_24 + l1*VACSS_1_24 + l2*VACSS_2r_24 + l3*VACSS_6r_24 ## THRESHOLDS: VACSS_1 | c(v1_thr, v1_thr)*t1 VACSS_1_4 | c(v1_thr, v1_thr)*t1 VACSS_1_8 | c(v1_thr, v1_thr)*t1 VACSS_1_12 | c(v1_thr, v1_thr)*t1 VACSS_1_18 | c(v1_thr, v1_thr)*t1 VACSS_1_24 | c(v1_thr, v1_thr)*t1 VACSS_2r | c(v2_thr1, v2_thr1)*t1 VACSS_2r | c(v2_thr2, v2_thr2)*t2 VACSS_2r | c(v2_thr3, v2_thr3)*t3 VACSS_2r | c(v2_thr4, v2_thr4)*t4 VACSS_2r | c(v2_thr5, v2_thr5)*t5 VACSS_2r_4 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_4 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_4 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_4 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_4 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_8 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_8 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_8 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_8 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_8 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_12 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_12 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_12 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_12 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_12 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_18 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_18 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_18 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_18 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_18 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_24 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_24 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_24 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_24 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_24 | c(v2_thr5, v2_thr5)*t5 VACSS_6r | c(v6_thr1, v6_thr1)*t1 VACSS_6r | c(v6_thr2, v6_thr2)*t2 VACSS_6r_4 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_4 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_8 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_8 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_12 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_12 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_18 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_18 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_24 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_24 | c(v6_thr2, v6_thr2)*t2 ## INTERCEPTS VACSS_1 ~ c(0, NA)*1 VACSS_2r ~ c(0, NA)*1 VACSS_6r ~ c(0, NA)*1 VACSS_1_4 ~ c(NA, NA)*1 VACSS_2r_4 ~ c(NA, NA)*1 VACSS_6r_4 ~ c(NA, NA)*1 VACSS_1_8 ~ c(NA, NA)*1 VACSS_2r_8 ~ c(NA, NA)*1 VACSS_6r_8 ~ c(NA, NA)*1 VACSS_1_12 ~ c(NA, NA)*1 VACSS_2r_12 ~ c(NA, NA)*1 VACSS_6r_12 ~ c(NA, NA)*1 VACSS_1_18 ~ c(NA, NA)*1 VACSS_2r_18 ~ c(NA, NA)*1 VACSS_6r_18 ~ c(NA, NA)*1 VACSS_1_24 ~ c(NA, NA)*1 VACSS_2r_24 ~ c(NA, NA)*1 VACSS_6r_24 ~ c(NA, NA)*1 ## SCALING PARAMETER VACSS_1 ~*~ c(1, NA)*VACSS_1 VACSS_2r ~*~ c(1, NA)*VACSS_2r VACSS_6r ~*~ c(1, NA)*VACSS_6r VACSS_1_4 ~*~ c(NA, NA)*VACSS_1_4 VACSS_1_8 ~*~ c(NA, NA)*VACSS_1_8 VACSS_1_12 ~*~ c(NA, NA)*VACSS_1_12 VACSS_1_18 ~*~ c(NA, NA)*VACSS_1_18 VACSS_1_24 ~*~ c(NA, NA)*VACSS_1_24 VACSS_2r_4 ~*~ c(NA, NA)*VACSS_2r_4 VACSS_6r_4 ~*~ c(NA, NA)*VACSS_6r_4 VACSS_2r_8 ~*~ c(NA, NA)*VACSS_2r_8 VACSS_6r_8 ~*~ c(NA, NA)*VACSS_6r_8 VACSS_2r_12 ~*~ c(NA, NA)*VACSS_2r_12 VACSS_6r_12 ~*~ c(NA, NA)*VACSS_6r_12 VACSS_2r_18 ~*~ c(NA, NA)*VACSS_2r_18 VACSS_6r_18 ~*~ c(NA, NA)*VACSS_6r_18 VACSS_2r_24 ~*~ c(NA, NA)*VACSS_2r_24 VACSS_6r_24 ~*~ c(NA, NA)*VACSS_6r_24 ## LATENT MEANS risk_1 ~ 0*1 risk_4 ~ 0*1 risk_8 ~ 0*1 risk_12 ~ 0*1 risk_18 ~ 0*1 risk_24 ~ 0*1 ## LV VARIANCE risk_1 ~~ c(1, NA)*risk_1 risk_4 ~~ c(NA, NA)*risk_4 risk_8 ~~ c(NA, NA)*risk_8 risk_12 ~~ c(NA, NA)*risk_12 risk_18 ~~ c(NA, NA)*risk_18 risk_24 ~~ c(NA, NA)*risk_24 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_8 ~~ VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_18 + VACSS_1_24 VACSS_1_18 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_8 ~~ VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_18 + VACSS_2r_24 VACSS_2r_18 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_8 ~~ VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_18 + VACSS_6r_24 VACSS_6r_18 ~~ VACSS_6r_24 " fit2_bipoc <- cfa(metric_model_g, data = df, ordered = T, estimator = "WLSMV", group = "bipoc", missing = 'pairwise', parameterization = 'delta') sum_fit2_bipoc <- summary(fit2_bipoc, fit.measures = T, standardized = T) ##### Test Model net(fit1_bipoc, fit2_bipoc) fit1_fit2_lrt_bipoc <- lavTestLRT(fit1_bipoc, fit2_bipoc) fit1_fit2_cfi_bipoc <- compare_fit(fit1_bipoc, fit2_bipoc) fit1_fit2_rmsea_bipoc <- compare_nested_models_RMSEA(fit1_bipoc, fit2_bipoc) ### Scalar Invariance ---- scalar_model_g <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + l1*VACSS_1 + l2*VACSS_2r + l3*VACSS_6r risk_4 =~ NA*VACSS_1_4 + l1*VACSS_1_4 + l2*VACSS_2r_4 + l3*VACSS_6r_4 risk_8 =~ NA*VACSS_1_8 + l1*VACSS_1_8 + l2*VACSS_2r_8 + l3*VACSS_6r_8 risk_12 =~ NA*VACSS_1_12 + l1*VACSS_1_12 + l2*VACSS_2r_12 + l3*VACSS_6r_12 risk_18 =~ NA*VACSS_1_18 + l1*VACSS_1_18 + l2*VACSS_2r_18 + l3*VACSS_6r_18 risk_24 =~ NA*VACSS_1_24 + l1*VACSS_1_24 + l2*VACSS_2r_24 + l3*VACSS_6r_24 ## THRESHOLDS: VACSS_1 | c(v1_thr, v1_thr)*t1 VACSS_1_4 | c(v1_thr, v1_thr)*t1 VACSS_1_8 | c(v1_thr, v1_thr)*t1 VACSS_1_12 | c(v1_thr, v1_thr)*t1 VACSS_1_18 | c(v1_thr, v1_thr)*t1 VACSS_1_24 | c(v1_thr, v1_thr)*t1 VACSS_2r | c(v2_thr1, v2_thr1)*t1 VACSS_2r | c(v2_thr2, v2_thr2)*t2 VACSS_2r | c(v2_thr3, v2_thr3)*t3 VACSS_2r | c(v2_thr4, v2_thr4)*t4 VACSS_2r | c(v2_thr5, v2_thr5)*t5 VACSS_2r_4 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_4 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_4 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_4 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_4 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_8 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_8 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_8 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_8 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_8 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_12 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_12 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_12 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_12 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_12 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_18 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_18 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_18 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_18 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_18 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_24 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_24 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_24 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_24 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_24 | c(v2_thr5, v2_thr5)*t5 VACSS_6r | c(v6_thr1, v6_thr1)*t1 VACSS_6r | c(v6_thr2, v6_thr2)*t2 VACSS_6r_4 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_4 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_8 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_8 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_12 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_12 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_18 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_18 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_24 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_24 | c(v6_thr2, v6_thr2)*t2 ## INTERCEPTS VACSS_1 ~ 0*1 VACSS_2r ~ 0*1 VACSS_6r ~ 0*1 VACSS_1_4 ~ 0*1 VACSS_2r_4 ~ 0*1 VACSS_6r_4 ~ 0*1 VACSS_1_8 ~ 0*1 VACSS_2r_8 ~ 0*1 VACSS_6r_8 ~ 0*1 VACSS_1_12 ~ 0*1 VACSS_2r_12 ~ 0*1 VACSS_6r_12 ~ 0*1 VACSS_1_18 ~ 0*1 VACSS_2r_18 ~ 0*1 VACSS_6r_18 ~ 0*1 VACSS_1_24 ~ 0*1 VACSS_2r_24 ~ 0*1 VACSS_6r_24 ~ 0*1 ## SCALING PARAMETER VACSS_1 ~*~ c(1, NA)*VACSS_1 VACSS_2r ~*~ c(1, NA)*VACSS_2r VACSS_6r ~*~ c(1, NA)*VACSS_6r VACSS_1_4 ~*~ c(NA, NA)*VACSS_1_4 VACSS_1_8 ~*~ c(NA, NA)*VACSS_1_8 VACSS_1_12 ~*~ c(NA, NA)*VACSS_1_12 VACSS_1_18 ~*~ c(NA, NA)*VACSS_1_18 VACSS_1_24 ~*~ c(NA, NA)*VACSS_1_24 VACSS_2r_4 ~*~ c(NA, NA)*VACSS_2r_4 VACSS_6r_4 ~*~ c(NA, NA)*VACSS_6r_4 VACSS_2r_8 ~*~ c(NA, NA)*VACSS_2r_8 VACSS_6r_8 ~*~ c(NA, NA)*VACSS_6r_8 VACSS_2r_12 ~*~ c(NA, NA)*VACSS_2r_12 VACSS_6r_12 ~*~ c(NA, NA)*VACSS_6r_12 VACSS_2r_18 ~*~ c(NA, NA)*VACSS_2r_18 VACSS_6r_18 ~*~ c(NA, NA)*VACSS_6r_18 VACSS_2r_24 ~*~ c(NA, NA)*VACSS_2r_24 VACSS_6r_24 ~*~ c(NA, NA)*VACSS_6r_24 ## LATENT MEANS risk_1 ~ c(0, NA)*1 risk_4 ~ c(NA, NA)*1 risk_8 ~ c(NA, NA)*1 risk_12 ~ c(NA, NA)*1 risk_18 ~ c(NA, NA)*1 risk_24 ~ c(NA, NA)*1 ## LV VARIANCE risk_1 ~~ c(1, NA)*risk_1 risk_4 ~~ c(NA, NA)*risk_4 risk_8 ~~ c(NA, NA)*risk_8 risk_12 ~~ c(NA, NA)*risk_12 risk_18 ~~ c(NA, NA)*risk_18 risk_24 ~~ c(NA, NA)*risk_24 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_8 ~~ VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_18 + VACSS_1_24 VACSS_1_18 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_8 ~~ VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_18 + VACSS_2r_24 VACSS_2r_18 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_8 ~~ VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_18 + VACSS_6r_24 VACSS_6r_18 ~~ VACSS_6r_24 " fit3_bipoc <- cfa(scalar_model_g, data = df, ordered = T, estimator = "WLSMV", group = "bipoc", missing = 'pairwise', parameterization = 'delta') sum_fit3_bipoc <- summary(fit3_bipoc, fit.measures = T, standardized = T) ##### Test Model net(fit2_bipoc, fit3_bipoc) fit2_fit3_lrt_bipoc <- lavTestLRT(fit2_bipoc, fit3_bipoc) fit2_fit3_cfi_bipoc <- compare_fit(fit2_bipoc, fit3_bipoc) fit2_fit3_rmsea_bipoc <- compare_nested_models_RMSEA(fit2_bipoc, fit3_bipoc) modindices(fit3_bipoc, op = '~1') %>% arrange(desc(`mi`)) # largest MI is VACSS_6r ### Partial Scalar Invariance ---- scalar_model_v2_g <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + l1*VACSS_1 + l2*VACSS_2r + l3*VACSS_6r risk_4 =~ NA*VACSS_1_4 + l1*VACSS_1_4 + l2*VACSS_2r_4 + l3*VACSS_6r_4 risk_8 =~ NA*VACSS_1_8 + l1*VACSS_1_8 + l2*VACSS_2r_8 + l3*VACSS_6r_8 risk_12 =~ NA*VACSS_1_12 + l1*VACSS_1_12 + l2*VACSS_2r_12 + l3*VACSS_6r_12 risk_18 =~ NA*VACSS_1_18 + l1*VACSS_1_18 + l2*VACSS_2r_18 + l3*VACSS_6r_18 risk_24 =~ NA*VACSS_1_24 + l1*VACSS_1_24 + l2*VACSS_2r_24 + l3*VACSS_6r_24 ## THRESHOLDS: VACSS_1 | c(v1_thr, v1_thr)*t1 VACSS_1_4 | c(v1_thr, v1_thr)*t1 VACSS_1_8 | c(v1_thr, v1_thr)*t1 VACSS_1_12 | c(v1_thr, v1_thr)*t1 VACSS_1_18 | c(v1_thr, v1_thr)*t1 VACSS_1_24 | c(v1_thr, v1_thr)*t1 VACSS_2r | c(v2_thr1, v2_thr1)*t1 VACSS_2r | c(v2_thr2, v2_thr2)*t2 VACSS_2r | c(v2_thr3, v2_thr3)*t3 VACSS_2r | c(v2_thr4, v2_thr4)*t4 VACSS_2r | c(v2_thr5, v2_thr5)*t5 VACSS_2r_4 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_4 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_4 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_4 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_4 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_8 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_8 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_8 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_8 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_8 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_12 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_12 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_12 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_12 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_12 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_18 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_18 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_18 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_18 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_18 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_24 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_24 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_24 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_24 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_24 | c(v2_thr5, v2_thr5)*t5 VACSS_6r | c(v6_thr1, v6_thr1)*t1 VACSS_6r | c(v6_thr2, v6_thr2)*t2 VACSS_6r_4 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_4 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_8 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_8 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_12 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_12 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_18 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_18 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_24 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_24 | c(v6_thr2, v6_thr2)*t2 ## INTERCEPTS VACSS_1 ~ 0*1 VACSS_2r ~ 0*1 VACSS_6r ~ c(mu6, mu6)*1 # freeing intercept, but stil making it equal across groups VACSS_1_4 ~ 0*1 VACSS_2r_4 ~ 0*1 VACSS_6r_4 ~ 0*1 VACSS_1_8 ~ 0*1 VACSS_2r_8 ~ 0*1 VACSS_6r_8 ~ 0*1 VACSS_1_12 ~ 0*1 VACSS_2r_12 ~ 0*1 VACSS_6r_12 ~ 0*1 VACSS_1_18 ~ 0*1 VACSS_2r_18 ~ 0*1 VACSS_6r_18 ~ 0*1 VACSS_1_24 ~ 0*1 VACSS_2r_24 ~ 0*1 VACSS_6r_24 ~ 0*1 ## SCALING PARAMETER # Variances at T1 to 1 VACSS_1 ~*~ c(1, NA)*VACSS_1 VACSS_2r ~*~ c(1, NA)*VACSS_2r VACSS_6r ~*~ c(1, NA)*VACSS_6r VACSS_1_4 ~*~ c(NA, NA)*VACSS_1_4 VACSS_1_8 ~*~ c(NA, NA)*VACSS_1_8 VACSS_1_12 ~*~ c(NA, NA)*VACSS_1_12 VACSS_1_18 ~*~ c(NA, NA)*VACSS_1_18 VACSS_1_24 ~*~ c(NA, NA)*VACSS_1_24 VACSS_2r_4 ~*~ c(NA, NA)*VACSS_2r_4 VACSS_6r_4 ~*~ c(NA, NA)*VACSS_6r_4 VACSS_2r_8 ~*~ c(NA, NA)*VACSS_2r_8 VACSS_6r_8 ~*~ c(NA, NA)*VACSS_6r_8 VACSS_2r_12 ~*~ c(NA, NA)*VACSS_2r_12 VACSS_6r_12 ~*~ c(NA, NA)*VACSS_6r_12 VACSS_2r_18 ~*~ c(NA, NA)*VACSS_2r_18 VACSS_6r_18 ~*~ c(NA, NA)*VACSS_6r_18 VACSS_2r_24 ~*~ c(NA, NA)*VACSS_2r_24 VACSS_6r_24 ~*~ c(NA, NA)*VACSS_6r_24 ## LATENT MEANS risk_1 ~ c(0, NA)*1 risk_4 ~ c(NA, NA)*1 risk_8 ~ c(NA, NA)*1 risk_12 ~ c(NA, NA)*1 risk_18 ~ c(NA, NA)*1 risk_24 ~ c(NA, NA)*1 ## LV VARIANCE risk_1 ~~ c(1, NA)*risk_1 risk_4 ~~ c(NA, NA)*risk_4 risk_8 ~~ c(NA, NA)*risk_8 risk_12 ~~ c(NA, NA)*risk_12 risk_18 ~~ c(NA, NA)*risk_18 risk_24 ~~ c(NA, NA)*risk_24 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_8 + VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_8 ~~ VACSS_1_12 + VACSS_1_18 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_18 + VACSS_1_24 VACSS_1_18 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_8 + VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_8 ~~ VACSS_2r_12 + VACSS_2r_18 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_18 + VACSS_2r_24 VACSS_2r_18 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_8 + VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_8 ~~ VACSS_6r_12 + VACSS_6r_18 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_18 + VACSS_6r_24 VACSS_6r_18 ~~ VACSS_6r_24 " fit3.1_bipoc <- cfa(scalar_model_v2_g, data = df, ordered = T, estimator = "WLSMV", group = "bipoc", missing = 'pairwise', parameterization = 'delta') sum_fit3.1_bipoc <- summary(fit3.1_bipoc, fit.measures = T, standardized = T) ##### Test Model net(fit2_bipoc, fit3.1_bipoc) fit2_fit3.1_lrt_bipoc <- lavTestLRT(fit2_bipoc, fit3.1_bipoc) fit2_fit3.1_cfi_bipoc <- compare_fit(fit2_bipoc, fit3.1_bipoc) fit2_fit3.1_rmsea_bipoc <- compare_nested_models_RMSEA(fit2_bipoc, fit3.1_bipoc) ## Randomization Group Invarinace ---- ### Baseline ---- fit0_rand <- cfa(long_cfa_g, data = df, ordered = T, estimator = "WLSMV", group = "randomization_group", missing = 'pairwise', parameterization = 'delta') sum_fit0_rand <- summary(fit0_rand, fit.measures = T, standardized = T) ### Threshold Invariance ---- fit1_rand <- cfa(thresh_model_g, data = df, ordered = T, estimator = "WLSMV", group = "randomization_group", missing = 'pairwise', parameterization = 'delta') sum_fit1_rand <- summary(fit1_rand, fit.measures = T, standardized = T) ##### Test Model net(fit0_rand, fit1_rand) fit0_fit1_lrt_rand <- lavTestLRT(fit0_rand, fit1_rand) fit0_fit1_cfi_rand <- compare_fit(fit0_rand, fit1_rand) fit0_fit1_rmsea_rand <- compare_nested_models_RMSEA(fit0_rand, fit1_rand) ### Metric Invariance ---- fit2_rand <- cfa(metric_model_g, data = df, ordered = T, estimator = "WLSMV", group = "randomization_group", missing = 'pairwise', parameterization = 'delta') sum_fit2_rand <- summary(fit2_rand, fit.measures = T, standardized = T) ##### Test Model net(fit1_rand, fit2_rand) fit1_fit2_lrt_rand <- lavTestLRT(fit1_rand, fit2_rand) fit1_fit2_cfi_rand <- compare_fit(fit1_rand, fit2_rand) fit1_fit2_rmsea_rand <- compare_nested_models_RMSEA(fit1_rand, fit2_rand) ### Scalar Invariance ---- fit3_rand <- cfa(scalar_model_g, data = df, ordered = T, estimator = "WLSMV", group = "randomization_group", missing = 'pairwise', parameterization = 'delta') sum_fit3_rand <- summary(fit3_rand, fit.measures = T, standardized = T) ##### Test Model net(fit2_rand, fit3_rand) fit2_fit3_lrt_rand <- lavTestLRT(fit2_rand, fit3_rand) fit2_fit3_cfi_rand <- compare_fit(fit2_rand, fit3_rand) fit2_fit3_rmsea_rand <- compare_nested_models_RMSEA(fit2_rand, fit3_rand) modindices(fit3_rand, op = '~1') %>% arrange(desc(`mi`)) # largest MI is VACSS_6r ### Partial Scalar Invariance ---- fit3.1_rand <- cfa(scalar_model_v2_g, data = df, ordered = T, estimator = "WLSMV", group = "randomization_group", missing = 'pairwise', parameterization = 'delta') sum_fit3.1_rand <- summary(fit3.1_rand, fit.measures = T, standardized = T) ##### Test Model net(fit2_rand, fit3.1_rand) fit2_fit3.1_lrt_rand <- lavTestLRT(fit2_rand, fit3.1_rand) fit2_fit3.1_cfi_rand <- compare_fit(fit2_rand, fit3.1_rand) fit2_fit3.1_rmsea_rand <- compare_nested_models_RMSEA(fit2_rand, fit3.1_rand) ## Sex x Time model ---- ### Baseline ---- long_cfa <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + VACSS_2r + VACSS_6r risk_4 =~ NA*VACSS_1_4 + VACSS_2r_4 + VACSS_6r_4 risk_12 =~ NA*VACSS_1_12 + VACSS_2r_12 + VACSS_6r_12 risk_24 =~ NA*VACSS_1_24 + VACSS_2r_24 + VACSS_6r_24 ## INTERCEPTS VACSS_1 ~ 0*1 VACSS_2r ~ 0*1 VACSS_6r ~ 0*1 VACSS_1_4 ~ 0*1 VACSS_2r_4 ~ 0*1 VACSS_6r_4 ~ 0*1 VACSS_1_12 ~ 0*1 VACSS_2r_12 ~ 0*1 VACSS_6r_12 ~ 0*1 VACSS_1_24 ~ 0*1 VACSS_2r_24 ~ 0*1 VACSS_6r_24 ~ 0*1 ## SCALING PARAMETER VACSS_1 ~*~ 1*VACSS_1 VACSS_2r ~*~ 1*VACSS_2r VACSS_6r ~*~ 1*VACSS_6r VACSS_1_4 ~*~ 1*VACSS_1_4 VACSS_1_12 ~*~ 1*VACSS_1_12 VACSS_1_24 ~*~ 1*VACSS_1_24 VACSS_2r_4 ~*~ 1*VACSS_2r_4 VACSS_6r_4 ~*~ 1*VACSS_6r_4 VACSS_2r_12 ~*~ 1*VACSS_2r_12 VACSS_6r_12 ~*~ 1*VACSS_6r_12 VACSS_2r_24 ~*~ 1*VACSS_2r_24 VACSS_6r_24 ~*~ 1*VACSS_6r_24 ## LV VARIANCE risk_1 ~~ 1*risk_1 risk_4 ~~ 1*risk_4 risk_12 ~~ 1*risk_12 risk_24 ~~ 1*risk_24 ## LATENT MEANS risk_1 ~ 0*1 risk_4 ~ 0*1 risk_12 ~ 0*1 risk_24 ~ 0*1 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_12 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_12 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_12 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_12 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_12 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_12 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_24 " fit0_female <- cfa(long_cfa, data = df, ordered = T, estimator = "WLSMV", group = "female", missing = 'pairwise', parameterization = 'delta') sum_fit0_female <- summary(fit0_female, fit.measures = T, standardized = T) ### Threshold Invariance ---- thresh_model <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + VACSS_2r + VACSS_6r risk_4 =~ NA*VACSS_1_4 + VACSS_2r_4 + VACSS_6r_4 risk_12 =~ NA*VACSS_1_12 + VACSS_2r_12 + VACSS_6r_12 risk_24 =~ NA*VACSS_1_24 + VACSS_2r_24 + VACSS_6r_24 ## THRESHOLDS: VACSS_1 | c(v1_thr, v1_thr)*t1 VACSS_1_4 | c(v1_thr, v1_thr)*t1 VACSS_1_12 | c(v1_thr, v1_thr)*t1 VACSS_1_24 | c(v1_thr, v1_thr)*t1 VACSS_2r | c(v2_thr1, v2_thr1)*t1 VACSS_2r | c(v2_thr2, v2_thr2)*t2 VACSS_2r | c(v2_thr3, v2_thr3)*t3 VACSS_2r | c(v2_thr4, v2_thr4)*t4 VACSS_2r | c(v2_thr5, v2_thr5)*t5 VACSS_2r_4 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_4 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_4 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_4 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_4 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_12 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_12 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_12 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_12 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_12 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_24 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_24 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_24 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_24 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_24 | c(v2_thr5, v2_thr5)*t5 VACSS_6r | c(v6_thr1, v6_thr1)*t1 VACSS_6r | c(v6_thr2, v6_thr2)*t2 VACSS_6r_4 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_4 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_12 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_12 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_24 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_24 | c(v6_thr2, v6_thr2)*t2 ## INTERCEPTS VACSS_1 ~ c(0, NA)*1 VACSS_2r ~ c(0, NA)*1 VACSS_6r ~ c(0, NA)*1 VACSS_1_4 ~ c(NA, NA)*1 VACSS_2r_4 ~ c(NA, NA)*1 VACSS_6r_4 ~ c(NA, NA)*1 VACSS_1_12 ~ c(NA, NA)*1 VACSS_2r_12 ~ c(NA, NA)*1 VACSS_6r_12 ~ c(NA, NA)*1 VACSS_1_24 ~ c(NA, NA)*1 VACSS_2r_24 ~ c(NA, NA)*1 VACSS_6r_24 ~ c(NA, NA)*1 ## SCALING PARAMETER VACSS_1 ~*~ c(1, 1)*VACSS_1 VACSS_2r ~*~ c(1, NA)*VACSS_2r VACSS_6r ~*~ c(1, NA)*VACSS_6r # scaling factor of 1 item at each time point to 1 VACSS_1_4 ~*~ c(1, 1)*VACSS_1_4 VACSS_1_12 ~*~ c(1, 1)*VACSS_1_12 VACSS_1_24 ~*~ c(1, 1)*VACSS_1_24 # free remaining scaling factors VACSS_2r_4 ~*~ c(NA, NA)*VACSS_2r_4 VACSS_6r_4 ~*~ c(NA, NA)*VACSS_6r_4 VACSS_2r_12 ~*~ c(NA, NA)*VACSS_2r_12 VACSS_6r_12 ~*~ c(NA, NA)*VACSS_6r_12 VACSS_2r_24 ~*~ c(NA, NA)*VACSS_2r_24 VACSS_6r_24 ~*~ c(NA, NA)*VACSS_6r_24 ## LATENT MEANS risk_1 ~ 0*1 risk_4 ~ 0*1 risk_12 ~ 0*1 risk_24 ~ 0*1 ## LV VARIANCE risk_1 ~~ 1*risk_1 risk_4 ~~ 1*risk_4 risk_12 ~~ 1*risk_12 risk_24 ~~ 1*risk_24 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_12 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_12 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_12 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_12 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_12 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_12 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_24 " fit1_female <- cfa(thresh_model, data = df, ordered = T, estimator = "WLSMV", group = "female", missing = 'pairwise', parameterization = 'delta') sum_fit1_female <- summary(fit1_female, fit.measures = T, standardized = T) ##### Test Model net(fit0_female, fit1_female) fit0_fit1_lrt_female <- lavTestLRT(fit0_female, fit1_female) fit0_fit1_cfi_female <- compare_fit(fit0_female, fit1_female) fit0_fit1_rmsea_female <- compare_nested_models_RMSEA(fit0_female, fit1_female) ### Metric Invariance ---- metric_model <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + l1*VACSS_1 + l2*VACSS_2r + l3*VACSS_6r risk_4 =~ NA*VACSS_1_4 + l1*VACSS_1_4 + l2*VACSS_2r_4 + l3*VACSS_6r_4 risk_12 =~ NA*VACSS_1_12 + l1*VACSS_1_12 + l2*VACSS_2r_12 + l3*VACSS_6r_12 risk_24 =~ NA*VACSS_1_24 + l1*VACSS_1_24 + l2*VACSS_2r_24 + l3*VACSS_6r_24 ## THRESHOLDS: VACSS_1 | c(v1_thr, v1_thr)*t1 VACSS_1_4 | c(v1_thr, v1_thr)*t1 VACSS_1_12 | c(v1_thr, v1_thr)*t1 VACSS_1_24 | c(v1_thr, v1_thr)*t1 VACSS_2r | c(v2_thr1, v2_thr1)*t1 VACSS_2r | c(v2_thr2, v2_thr2)*t2 VACSS_2r | c(v2_thr3, v2_thr3)*t3 VACSS_2r | c(v2_thr4, v2_thr4)*t4 VACSS_2r | c(v2_thr5, v2_thr5)*t5 VACSS_2r_4 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_4 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_4 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_4 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_4 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_12 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_12 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_12 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_12 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_12 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_24 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_24 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_24 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_24 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_24 | c(v2_thr5, v2_thr5)*t5 VACSS_6r | c(v6_thr1, v6_thr1)*t1 VACSS_6r | c(v6_thr2, v6_thr2)*t2 VACSS_6r_4 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_4 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_12 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_12 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_24 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_24 | c(v6_thr2, v6_thr2)*t2 ## INTERCEPTS VACSS_1 ~ c(0, NA)*1 VACSS_2r ~ c(0, NA)*1 VACSS_6r ~ c(0, NA)*1 VACSS_1_4 ~ c(NA, NA)*1 VACSS_2r_4 ~ c(NA, NA)*1 VACSS_6r_4 ~ c(NA, NA)*1 VACSS_1_12 ~ c(NA, NA)*1 VACSS_2r_12 ~ c(NA, NA)*1 VACSS_6r_12 ~ c(NA, NA)*1 VACSS_1_24 ~ c(NA, NA)*1 VACSS_2r_24 ~ c(NA, NA)*1 VACSS_6r_24 ~ c(NA, NA)*1 ## SCALING PARAMETER VACSS_1 ~*~ c(1, NA)*VACSS_1 VACSS_2r ~*~ c(1, NA)*VACSS_2r VACSS_6r ~*~ c(1, NA)*VACSS_6r VACSS_1_4 ~*~ c(NA, NA)*VACSS_1_4 VACSS_1_12 ~*~ c(NA, NA)*VACSS_1_12 VACSS_1_24 ~*~ c(NA, NA)*VACSS_1_24 VACSS_2r_4 ~*~ c(NA, NA)*VACSS_2r_4 VACSS_6r_4 ~*~ c(NA, NA)*VACSS_6r_4 VACSS_2r_12 ~*~ c(NA, NA)*VACSS_2r_12 VACSS_6r_12 ~*~ c(NA, NA)*VACSS_6r_12 VACSS_2r_24 ~*~ c(NA, NA)*VACSS_2r_24 VACSS_6r_24 ~*~ c(NA, NA)*VACSS_6r_24 ## LATENT MEANS risk_1 ~ 0*1 risk_4 ~ 0*1 risk_12 ~ 0*1 risk_24 ~ 0*1 ## LV VARIANCE risk_1 ~~ c(1, NA)*risk_1 risk_4 ~~ c(NA, NA)*risk_4 risk_12 ~~ c(NA, NA)*risk_12 risk_24 ~~ c(NA, NA)*risk_24 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_12 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_12 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_12 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_12 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_12 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_12 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_24 " fit2_female <- cfa(metric_model, data = df, ordered = T, estimator = "WLSMV", group = "female", missing = 'pairwise', parameterization = 'delta') sum_fit2_female <- summary(fit2_female, fit.measures = T, standardized = T) ##### Test Model net(fit1_female, fit2_female) fit1_fit2_lrt_female <- lavTestLRT(fit1_female, fit2_female) fit1_fit2_cfi_female <- compare_fit(fit1_female, fit2_female) fit1_fit2_rmsea_female <- compare_nested_models_RMSEA(fit1_female, fit2_female) ### Scalar Invariance ---- scalar_model <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + l1*VACSS_1 + l2*VACSS_2r + l3*VACSS_6r risk_4 =~ NA*VACSS_1_4 + l1*VACSS_1_4 + l2*VACSS_2r_4 + l3*VACSS_6r_4 risk_12 =~ NA*VACSS_1_12 + l1*VACSS_1_12 + l2*VACSS_2r_12 + l3*VACSS_6r_12 risk_24 =~ NA*VACSS_1_24 + l1*VACSS_1_24 + l2*VACSS_2r_24 + l3*VACSS_6r_24 ## THRESHOLDS: VACSS_1 | c(v1_thr, v1_thr)*t1 VACSS_1_4 | c(v1_thr, v1_thr)*t1 VACSS_1_12 | c(v1_thr, v1_thr)*t1 VACSS_1_24 | c(v1_thr, v1_thr)*t1 VACSS_2r | c(v2_thr1, v2_thr1)*t1 VACSS_2r | c(v2_thr2, v2_thr2)*t2 VACSS_2r | c(v2_thr3, v2_thr3)*t3 VACSS_2r | c(v2_thr4, v2_thr4)*t4 VACSS_2r | c(v2_thr5, v2_thr5)*t5 VACSS_2r_4 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_4 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_4 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_4 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_4 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_12 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_12 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_12 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_12 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_12 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_24 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_24 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_24 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_24 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_24 | c(v2_thr5, v2_thr5)*t5 VACSS_6r | c(v6_thr1, v6_thr1)*t1 VACSS_6r | c(v6_thr2, v6_thr2)*t2 VACSS_6r_4 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_4 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_12 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_12 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_24 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_24 | c(v6_thr2, v6_thr2)*t2 ## INTERCEPTS VACSS_1 ~ 0*1 VACSS_2r ~ 0*1 VACSS_6r ~ 0*1 VACSS_1_4 ~ 0*1 VACSS_2r_4 ~ 0*1 VACSS_6r_4 ~ 0*1 VACSS_1_12 ~ 0*1 VACSS_2r_12 ~ 0*1 VACSS_6r_12 ~ 0*1 VACSS_1_24 ~ 0*1 VACSS_2r_24 ~ 0*1 VACSS_6r_24 ~ 0*1 ## SCALING PARAMETER VACSS_1 ~*~ c(1, NA)*VACSS_1 VACSS_2r ~*~ c(1, NA)*VACSS_2r VACSS_6r ~*~ c(1, NA)*VACSS_6r VACSS_1_4 ~*~ c(NA, NA)*VACSS_1_4 VACSS_1_12 ~*~ c(NA, NA)*VACSS_1_12 VACSS_1_24 ~*~ c(NA, NA)*VACSS_1_24 VACSS_2r_4 ~*~ c(NA, NA)*VACSS_2r_4 VACSS_6r_4 ~*~ c(NA, NA)*VACSS_6r_4 VACSS_2r_12 ~*~ c(NA, NA)*VACSS_2r_12 VACSS_6r_12 ~*~ c(NA, NA)*VACSS_6r_12 VACSS_2r_24 ~*~ c(NA, NA)*VACSS_2r_24 VACSS_6r_24 ~*~ c(NA, NA)*VACSS_6r_24 ## LATENT MEANS risk_1 ~ c(0, NA)*1 risk_4 ~ c(NA, NA)*1 risk_12 ~ c(NA, NA)*1 risk_24 ~ c(NA, NA)*1 ## LV VARIANCE risk_1 ~~ c(1, NA)*risk_1 risk_4 ~~ c(NA, NA)*risk_4 risk_12 ~~ c(NA, NA)*risk_12 risk_24 ~~ c(NA, NA)*risk_24 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_12 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_12 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_12 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_12 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_12 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_12 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_24 " fit3_female <- cfa(scalar_model, data = df, ordered = T, estimator = "WLSMV", group = "female", missing = 'pairwise', parameterization = 'delta') sum_fit3_female <- summary(fit3_female, fit.measures = T, standardized = T) ##### Test Model net(fit2_female, fit3_female) fit2_fit3_lrt_female <- lavTestLRT(fit2_female, fit3_female) fit2_fit3_cfi_female <- compare_fit(fit2_female, fit3_female) fit2_fit3_rmsea_female <- compare_nested_models_RMSEA(fit2_female, fit3_female) modindices(fit3_female, op = '~1') %>% arrange(desc(`mi`)) # largest MI is VACSS_6r ### Partial Scalar Invariance ---- scalar_model_v2 <- " ## LOADINGS risk_1 =~ NA*VACSS_1 + l1*VACSS_1 + l2*VACSS_2r + l3*VACSS_6r risk_4 =~ NA*VACSS_1_4 + l1*VACSS_1_4 + l2*VACSS_2r_4 + l3*VACSS_6r_4 risk_12 =~ NA*VACSS_1_12 + l1*VACSS_1_12 + l2*VACSS_2r_12 + l3*VACSS_6r_12 risk_24 =~ NA*VACSS_1_24 + l1*VACSS_1_24 + l2*VACSS_2r_24 + l3*VACSS_6r_24 ## THRESHOLDS: VACSS_1 | c(v1_thr, v1_thr)*t1 VACSS_1_4 | c(v1_thr, v1_thr)*t1 VACSS_1_12 | c(v1_thr, v1_thr)*t1 VACSS_1_24 | c(v1_thr, v1_thr)*t1 VACSS_2r | c(v2_thr1, v2_thr1)*t1 VACSS_2r | c(v2_thr2, v2_thr2)*t2 VACSS_2r | c(v2_thr3, v2_thr3)*t3 VACSS_2r | c(v2_thr4, v2_thr4)*t4 VACSS_2r | c(v2_thr5, v2_thr5)*t5 VACSS_2r_4 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_4 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_4 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_4 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_4 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_12 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_12 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_12 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_12 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_12 | c(v2_thr5, v2_thr5)*t5 VACSS_2r_24 | c(v2_thr1, v2_thr1)*t1 VACSS_2r_24 | c(v2_thr2, v2_thr2)*t2 VACSS_2r_24 | c(v2_thr3, v2_thr3)*t3 VACSS_2r_24 | c(v2_thr4, v2_thr4)*t4 VACSS_2r_24 | c(v2_thr5, v2_thr5)*t5 VACSS_6r | c(v6_thr1, v6_thr1)*t1 VACSS_6r | c(v6_thr2, v6_thr2)*t2 VACSS_6r_4 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_4 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_12 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_12 | c(v6_thr2, v6_thr2)*t2 VACSS_6r_24 | c(v6_thr1, v6_thr1)*t1 VACSS_6r_24 | c(v6_thr2, v6_thr2)*t2 ## INTERCEPTS VACSS_1 ~ 0*1 VACSS_2r ~ 0*1 VACSS_6r ~ c(mu6, mu6)*1 VACSS_1_4 ~ 0*1 VACSS_2r_4 ~ 0*1 VACSS_6r_4 ~ 0*1 VACSS_1_12 ~ 0*1 VACSS_2r_12 ~ 0*1 VACSS_6r_12 ~ 0*1 VACSS_1_24 ~ 0*1 VACSS_2r_24 ~ 0*1 VACSS_6r_24 ~ 0*1 ## SCALING PARAMETER VACSS_1 ~*~ c(1, NA)*VACSS_1 VACSS_2r ~*~ c(1, NA)*VACSS_2r VACSS_6r ~*~ c(1, NA)*VACSS_6r VACSS_1_4 ~*~ c(NA, NA)*VACSS_1_4 VACSS_1_12 ~*~ c(NA, NA)*VACSS_1_12 VACSS_1_24 ~*~ c(NA, NA)*VACSS_1_24 VACSS_2r_4 ~*~ c(NA, NA)*VACSS_2r_4 VACSS_6r_4 ~*~ c(NA, NA)*VACSS_6r_4 VACSS_2r_12 ~*~ c(NA, NA)*VACSS_2r_12 VACSS_6r_12 ~*~ c(NA, NA)*VACSS_6r_12 VACSS_2r_24 ~*~ c(NA, NA)*VACSS_2r_24 VACSS_6r_24 ~*~ c(NA, NA)*VACSS_6r_24 ## LATENT MEANS risk_1 ~ c(0, NA)*1 risk_4 ~ c(NA, NA)*1 risk_12 ~ c(NA, NA)*1 risk_24 ~ c(NA, NA)*1 ## LV VARIANCE risk_1 ~~ c(1, NA)*risk_1 risk_4 ~~ c(NA, NA)*risk_4 risk_12 ~~ c(NA, NA)*risk_12 risk_24 ~~ c(NA, NA)*risk_24 ## OV AUTO-CORRELATIONS VACSS_1 ~~ VACSS_1_4 + VACSS_1_12 + VACSS_1_24 VACSS_1_4 ~~ VACSS_1_12 + VACSS_1_24 VACSS_1_12 ~~ VACSS_1_24 VACSS_2r ~~ VACSS_2r_4 + VACSS_2r_12 + VACSS_2r_24 VACSS_2r_4 ~~ VACSS_2r_12 + VACSS_2r_24 VACSS_2r_12 ~~ VACSS_2r_24 VACSS_6r ~~ VACSS_6r_4 + VACSS_6r_12 + VACSS_6r_24 VACSS_6r_4 ~~ VACSS_6r_12 + VACSS_6r_24 VACSS_6r_12 ~~ VACSS_6r_24 " fit3.1_female <- cfa(scalar_model_v2, data = df, ordered = T, estimator = "WLSMV", group = "female", missing = 'pairwise', parameterization = 'delta') sum_fit3.1_female <- summary(fit3.1_female, fit.measures = T, standardized = T) ##### Test Model net(fit2_female, fit3.1_female) fit2_fit3.1_lrt_female <- lavTestLRT(fit2_female, fit3.1_female) fit2_fit3.1_cfi_female <- compare_fit(fit2_female, fit3.1_female) fit2_fit3.1_rmsea_female <- compare_nested_models_RMSEA(fit2_female, fit3.1_female) # *********** ----- # Multiple Imputation Analyses ---- # *********** ----- # See # of missing responses for each item varlist(df, pattern = 'vacss_1|vacss_2r|vacss_6r', exclude = '_36$') %>% {df[,.]} %>% sapply(function(x) mean(is.na(x))) # Proportion missing per variable # Prepare data for imputation ---- # Strategy: Impute raw dichotomous items, then derive composite variables (VACSS_2r, VACSS_6r) # This preserves the logical dependencies between items. analysis_vars <- c( "female", "bipoc", # Item 1 "VACSS_1", "VACSS_1_4", "VACSS_1_8", "VACSS_1_12", "VACSS_1_18", "VACSS_1_24", # Raw items for Item 2 Composite "VACSS_2", "VACSS_3", "VACSS_4", "VACSS_5", "VACSS_5a", "VACSS_2_4", "VACSS_3_4", "VACSS_4_4", "VACSS_5_4", "VACSS_5a_4", "VACSS_2_8", "VACSS_3_8", "VACSS_4_8", "VACSS_5_8", "VACSS_5a_8", "VACSS_2_12", "VACSS_3_12", "VACSS_4_12", "VACSS_5_12", "VACSS_5a_12", "VACSS_2_18", "VACSS_3_18", "VACSS_4_18", "VACSS_5_18", "VACSS_5a_18", "VACSS_2_24", "VACSS_3_24", "VACSS_4_24", "VACSS_5_24", "VACSS_5a_24", # Raw items for Item 6 Composite "VACSS_6", "VACSS_6a", "VACSS_6_4", "VACSS_6a_4", "VACSS_6_8", "VACSS_6a_8", "VACSS_6_12", "VACSS_6a_12", "VACSS_6_18", "VACSS_6a_18", "VACSS_6_24", "VACSS_6a_24" ) # Auxiliary variables aux_vars <- c( "age_consent", #cont "siteid", #nominal "randomization_group", #dich "Ethnicity", #nominal varlist(df, pattern = 'phq9', exclude = '_36$'), #cont varlist(df, pattern = 'gad7', exclude = '_36$'), #cont varlist(df, pattern = 'pcl5_[0-9].*', exclude = '_36$'), #cont varlist(df, pattern = 'vr12', exclude = 'mcs|pcs|_36$'), #cont varlist(df, pattern = 'socialsupport', exclude = '_36$'), #cont varlist(df, pattern = 'aud', exclude = '_'), #dich varlist(df, 'drug_2'), #cont varlist(df, 'tobacco', exclude = 'tobacco_2|tobacco_3|_36$'), #cont varlist(df, "Current_Med_Non_Adherence", exclude = '_36$'), #dich varlist(df, "Current_Med_Side_Effects_Significant", exclude = '_36$'), #dich varlist(df, "sat_", exclude = '_36$'), #cont "prov_degree", #nom "ref_facility", #nom "prac_loc", #nom "prov_type", #nom "Covid", #dich "treatment_confidence", #cont "patient_confidence", #cont "Marital", #nom "DEM_5", #dich "full_time", #dich "DEM_6", #nom "DEM_7", #cont "years_serv" #cont ) # dichotomous vars dich <- c( analysis_vars, # All primary outcomes are dichotomous (0,1) "randomization_group", varlist(df, pattern = 'aud', exclude = '_'), varlist(df, "Current_Med_Non_Adherence", exclude = '_36$'), #dich varlist(df, "Current_Med_Side_Effects_Significant", exclude = '_36$'), #dich "Covid", "DEM_5", "full_time" ) nom <- c( "siteid", #nominal "Ethnicity", #nominal "prov_degree", #nom "ref_facility", #nom "prac_loc", #nom "prov_type", #nom "Marital", #nom "DEM_6" #nom ) cont <- c( "age_consent", #cont varlist(df, pattern = 'phq9', exclude = '_36$'), #cont varlist(df, pattern = 'gad7', exclude = '_36$'), #cont varlist(df, pattern = 'pcl5_[0-9].*', exclude = '_36$'), #cont varlist(df, pattern = 'vr12', exclude = 'mcs|pcs|_36$'), #cont varlist(df, pattern = 'socialsupport', exclude = '_36$'), #cont varlist(df, 'drug_2'), #cont varlist(df, 'tobacco', exclude = 'tobacco_2|tobacco_3|_36$'), #cont varlist(df, "sat_", exclude = '_36$'), #cont "treatment_confidence", #cont "patient_confidence", #cont "DEM_7", #cont "years_serv" #cont ) # Union of VACSS + auxiliaries for imputation input (avoids having same variable present twice) dat <- df %>% select(all_of(union(analysis_vars, aux_vars))) # Vars only measured at T1 (include in all predictions) t1_only <- c( "female", "bipoc", "age_consent", #cont "siteid", #nominal "randomization_group", #dich "Ethnicity", #nominal "prov_degree", #nom "ref_facility", #nom "prac_loc", #nom "prov_type", #nom "Covid", #dich "treatment_confidence", #cont "patient_confidence", #cont "Marital", #nom "DEM_5", #dich "full_time", #dich "DEM_6", #nom "DEM_7", #cont "years_serv", #cont varlist(df, pattern = 'aud', exclude = '_'), #dich varlist(df, 'drug_2') #cont ) # Baseline vars baseline <- c( varlist(dat, pattern = 'vacss', exclude = 'vacss_.{1,2}_[0-9]{1,2}'), varlist(dat, pattern = 'phq9', exclude = 'phq9_.{1,2}_[0-9]{1,2}'), #cont varlist(dat, pattern = 'gad7_[1-8]$'), #cont varlist(dat, pattern = 'pcl5_[0-9]{1,2}$*'), #cont varlist(dat, pattern = 'vr12_.{1,2}$'), #cont varlist(dat, '^tobacco$') #cont ) # wk 4 wk4 <- c( varlist(dat, pattern = 'vacss_.{1,2}_4$'), varlist(dat, pattern = 'phq9_[0-9]{1,2}_4$'), #cont varlist(dat, pattern = 'gad7_[0-9]_4$'), #cont varlist(dat, pattern = 'pcl5_[0-9]{1,2}_4$'), #cont varlist(dat, pattern = 'vr12_.{1,2}_4$'), #cont varlist(dat, pattern = 'socialsupport_[0-9]_4$'), #cont varlist(dat, 'tobacco_1_4'), #cont varlist(dat, "Current_Med_Non_Adherence_4"), #dich varlist(dat, "Current_Med_Side_Effects_Significant_4"), #dich varlist(dat, "sat_[1-2]_4$") #cont ) # wk 8 wk8 <- c( varlist(dat, pattern = 'vacss_.{1,2}_8$'), varlist(dat, pattern = 'phq9_[0-9]{1,2}_8$'), #cont varlist(dat, pattern = 'gad7_[0-9]_8$'), #cont varlist(dat, pattern = 'pcl5_[0-9]{1,2}_8$'), #cont varlist(dat, pattern = 'vr12_.{1,2}_8$'), #cont varlist(dat, pattern = 'socialsupport_[0-9]_8$'), #cont varlist(dat, 'tobacco_1_8'), #cont varlist(dat, "Current_Med_Non_Adherence_8"), #dich varlist(dat, "Current_Med_Side_Effects_Significant_8"), #dich varlist(dat, "sat_[1-2]_8$") #cont ) # wk 12 wk12 <- c( varlist(dat, pattern = 'vacss_.{1,2}_12$'), varlist(dat, pattern = 'phq9_[0-9]{1,2}_12$'), #cont varlist(dat, pattern = 'gad7_[0-9]_12$'), #cont varlist(dat, pattern = 'pcl5_[0-9]{1,2}_12$'), #cont varlist(dat, pattern = 'vr12_.{1,2}_12$'), #cont varlist(dat, pattern = 'socialsupport_[0-9]_12$'), #cont varlist(dat, 'tobacco_1_12'), #cont varlist(dat, "Current_Med_Non_Adherence_12"), #dich varlist(dat, "Current_Med_Side_Effects_Significant_12"), #dich varlist(dat, "sat_[1-2]_12$") #cont ) # wk 18 wk18 <- c( varlist(dat, pattern = 'vacss_.{1,2}_18$'), varlist(dat, pattern = 'phq9_[0-9]{1,2}_18$'), #cont varlist(dat, pattern = 'gad7_[0-9]_18$'), #cont varlist(dat, pattern = 'pcl5_[0-9]{1,2}_18$'), #cont varlist(dat, pattern = 'vr12_.{1,2}_18$'), #cont varlist(dat, pattern = 'socialsupport_[0-9]_18$'), #cont varlist(dat, 'tobacco_1_18'), #cont varlist(dat, "Current_Med_Non_Adherence_18"), #dich varlist(dat, "Current_Med_Side_Effects_Significant_18"), #dich varlist(dat, "sat_[1-2]_18$") #cont ) # wk 24 wk24 <- c( varlist(dat, pattern = 'vacss_.{1,2}_24$'), varlist(dat, pattern = 'phq9_[0-9]{1,2}_24$'), #cont varlist(dat, pattern = 'gad7_[0-9]_24$'), #cont varlist(dat, pattern = 'pcl5_[0-9]{1,2}_24$'), #cont varlist(dat, pattern = 'vr12_.{1,2}_24$'), #cont varlist(dat, pattern = 'socialsupport_[0-9]_24$'), #cont varlist(dat, 'tobacco_1_24'), #cont varlist(dat, "Current_Med_Non_Adherence_24"), #dich varlist(dat, "Current_Med_Side_Effects_Significant_24"), #dich varlist(dat, "sat_[1-2]_24$") #cont ) # Data Cleaning & Type Conversion ---- # mice uses the variable class to determine the default imputation method. # Converting dichotomous variables to factors ensures 'logreg' is used. dat <- dat %>% mutate(across(all_of(dich), as.factor), across(all_of(nom), as.factor), across(all_of(cont), as.numeric)) # Setup MICE Configuration ---- # Initialize with 0 iterations to extract the default setup ini <- mice(dat, maxit = 0) meth <- ini$method pred <- ini$predictorMatrix pred[,] <- 0 # start with prediction matrix that is all zeros pred[c(t1_only, baseline), c(t1_only, baseline)] <- 1 # allow t1 and baseline items to predict themselves pred[c(t1_only, wk4), c(t1_only, wk4)] <- 1 # allow t1 and wk4 items to predict themselves pred[c(t1_only, wk8), c(t1_only, wk8)] <- 1 pred[c(t1_only, wk12), c(t1_only, wk12)] <- 1 pred[c(t1_only, wk18), c(t1_only, wk18)] <- 1 pred[c(t1_only, wk24), c(t1_only, wk24)] <- 1 diag(pred) <- 0 # Identify all variables that have missingness vars_with_missing <- names(which(colSums(is.na(dat)) > 0)) # Assign methods only to analysis variables that have missing data meth[intersect(vars_with_missing, dich)] <- "logreg" meth[intersect(vars_with_missing, analysis_vars)] <- "pmm" # using pmm for VACSS data, provided better match to observed distribution meth[intersect(vars_with_missing, nom)] <- "polyreg" meth[intersect(vars_with_missing, cont)] <- "pmm" # Run Imputation ---- # Using futuremice for parallel processing to speed up the 25 imputations imp <- futuremice( dat, m = 25, parallelseed = 123, n.core = 5, future.plan = "multisession", method = meth, predictorMatrix = pred, maxit = 20 ) # Post-Imputation: Derive Composites ---- # After imputation, we calculate VACSS_2r and VACSS_6r within each imputed dataset imp_long <- complete(imp, action = "long", include = T) # Recode scores to 0 where where item 2 or 6 were denied imp_long <- mutate(imp_long, VACSS_3 = ifelse(VACSS_2 == 0, 0, VACSS_3), VACSS_4 = ifelse(VACSS_2 == 0, 0, VACSS_4), VACSS_5 = ifelse(VACSS_2 == 0, 0, VACSS_5), VACSS_5a = ifelse(VACSS_5 == 0, 0, VACSS_5a), VACSS_6a = ifelse(VACSS_6 == 0, 0, VACSS_6a), #wk 4 VACSS_3_4 = ifelse(VACSS_2_4 == 0, 0, VACSS_3_4), VACSS_4_4 = ifelse(VACSS_2_4 == 0, 0, VACSS_4_4), VACSS_5_4 = ifelse(VACSS_2_4 == 0, 0, VACSS_5_4), VACSS_5a_4 = ifelse(VACSS_5_4 == 0, 0, VACSS_5a_4), VACSS_6a_4 = ifelse(VACSS_6_4 == 0, 0, VACSS_6a_4), #wk 8 VACSS_3_8 = ifelse(VACSS_2_8 == 0, 0, VACSS_3_8), VACSS_4_8 = ifelse(VACSS_2_8 == 0, 0, VACSS_4_8), VACSS_5_8 = ifelse(VACSS_2_8 == 0, 0, VACSS_5_8), VACSS_5a_8 = ifelse(VACSS_5_8 == 0, 0, VACSS_5a_8), VACSS_6a_8 = ifelse(VACSS_6_8 == 0, 0, VACSS_6a_8), #wk 12 VACSS_3_12 = ifelse(VACSS_2_12 == 0, 0, VACSS_3_12), VACSS_4_12 = ifelse(VACSS_2_12 == 0, 0, VACSS_4_12), VACSS_5_12 = ifelse(VACSS_2_12 == 0, 0, VACSS_5_12), VACSS_5a_12 = ifelse(VACSS_5_12 == 0, 0, VACSS_5a_12), VACSS_6a_12 = ifelse(VACSS_6_12 == 0, 0, VACSS_6a_12), #wk 18 VACSS_3_18 = ifelse(VACSS_2_18 == 0, 0, VACSS_3_18), VACSS_4_18 = ifelse(VACSS_2_18 == 0, 0, VACSS_4_18), VACSS_5_18 = ifelse(VACSS_2_18 == 0, 0, VACSS_5_18), VACSS_5a_18 = ifelse(VACSS_5_18 == 0, 0, VACSS_5a_18), VACSS_6a_18 = ifelse(VACSS_6_18 == 0, 0, VACSS_6a_18), #wk 24 VACSS_3_24 = ifelse(VACSS_2_24 == 0, 0, VACSS_3_24), VACSS_4_24 = ifelse(VACSS_2_24 == 0, 0, VACSS_4_24), VACSS_5_24 = ifelse(VACSS_2_24 == 0, 0, VACSS_5_24), VACSS_5a_24 = ifelse(VACSS_5_24 == 0, 0, VACSS_5a_24), VACSS_6a_24 = ifelse(VACSS_6_24 == 0, 0, VACSS_6a_24)) # Recoded C-SSRS #2 imp_long <- imp_long %>% mutate(VACSS_2r = ifelse(VACSS_5a == 1, 5, ifelse(VACSS_5 == 1, 4, ifelse(VACSS_4 == 1, 3, ifelse(VACSS_3 == 1, 2, ifelse(VACSS_2 == 1, 1, 0))))), VACSS_2r_4 = ifelse(VACSS_5a_4 == 1, 5, ifelse(VACSS_5_4 == 1, 4, ifelse(VACSS_4_4 == 1, 3, ifelse(VACSS_3_4 == 1, 2, ifelse(VACSS_2_4 == 1, 1, 0))))), VACSS_2r_8 = ifelse(VACSS_5a_8 == 1, 5, ifelse(VACSS_5_8 == 1, 4, ifelse(VACSS_4_8 == 1, 3, ifelse(VACSS_3_8 == 1, 2, ifelse(VACSS_2_8 == 1, 1, 0))))), VACSS_2r_12 = ifelse(VACSS_5a_12 == 1, 5, ifelse(VACSS_5_12 == 1, 4, ifelse(VACSS_4_12 == 1, 3, ifelse(VACSS_3_12 == 1, 2, ifelse(VACSS_2_12 == 1, 1, 0))))), VACSS_2r_18 = ifelse(VACSS_5a_18 == 1, 5, ifelse(VACSS_5_18 == 1, 4, ifelse(VACSS_4_18 == 1, 3, ifelse(VACSS_3_18 == 1, 2, ifelse(VACSS_2_18 == 1, 1, 0))))), VACSS_2r_24 = ifelse(VACSS_5a_24 == 1, 5, ifelse(VACSS_5_24 == 1, 4, ifelse(VACSS_4_24 == 1, 3, ifelse(VACSS_3_24 == 1, 2, ifelse(VACSS_2_24 == 1, 1, 0)))))) # Recoded C-SSRS #6 imp_long <- imp_long %>% mutate(VACSS_6r = ifelse(VACSS_6a == 1, 2, ifelse(VACSS_6 == 1, 1, 0)), VACSS_6r_4 = ifelse(VACSS_6a_4 == 1, 2, ifelse(VACSS_6_4 == 1, 1, 0)), VACSS_6r_8 = ifelse(VACSS_6a_8 == 1, 2, ifelse(VACSS_6_8 == 1, 1, 0)), VACSS_6r_12 = ifelse(VACSS_6a_12 == 1, 2, ifelse(VACSS_6_12 == 1, 1, 0)), VACSS_6r_18 = ifelse(VACSS_6a_18 == 1, 2, ifelse(VACSS_6_18 == 1, 1, 0)), VACSS_6r_24 = ifelse(VACSS_6a_24 == 1, 2, ifelse(VACSS_6_24 == 1, 1, 0))) imp_long <- imp_long %>% mutate(VACSS_1 = factor(VACSS_1, levels = 0:1, ordered = T), VACSS_1_4 = factor(VACSS_1_4, levels = 0:1, ordered = T), VACSS_1_8 = factor(VACSS_1_8, levels = 0:1, ordered = T), VACSS_1_12 = factor(VACSS_1_12, levels = 0:1, ordered = T), VACSS_1_18 = factor(VACSS_1_18, levels = 0:1, ordered = T), VACSS_1_24 = factor(VACSS_1_24, levels = 0:1, ordered = T), VACSS_2r = factor(VACSS_2r, levels = 0:5, ordered = T), VACSS_2r_4 = factor(VACSS_2r_4, levels = 0:5, ordered = T), VACSS_2r_8 = factor(VACSS_2r_8, levels = 0:5, ordered = T), VACSS_2r_12 = factor(VACSS_2r_12, levels = 0:5, ordered = T), VACSS_2r_18 = factor(VACSS_2r_18, levels = 0:5, ordered = T), VACSS_2r_24 = factor(VACSS_2r_24, levels = 0:5, ordered = T), VACSS_6r = factor(VACSS_6r, levels = 0:2, ordered = T), VACSS_6r_4 = factor(VACSS_6r_4, levels = 0:2, ordered = T), VACSS_6r_8 = factor(VACSS_6r_8, levels = 0:2, ordered = T), VACSS_6r_12 = factor(VACSS_6r_12, levels = 0:2, ordered = T), VACSS_6r_18 = factor(VACSS_6r_18, levels = 0:2, ordered = T), VACSS_6r_24 = factor(VACSS_6r_24, levels = 0:2, ordered = T) ) imp_final <- as.mids(imp_long) # Diagnostic plot ---- # ---- Helper: extract long data for a set of variables prep_imp_data <- function(mids_obj, vars, wave_labels) { long <- complete(mids_obj, action = "long", include = TRUE) %>% dplyr::select(.imp, all_of(vars)) %>% tidyr::pivot_longer(cols = all_of(vars), names_to = "variable", values_to = "value") %>% dplyr::mutate( source = ifelse(.imp == 0, "Observed", "Imputed"), variable = factor(variable, levels = vars, labels = wave_labels), value = as.numeric(as.character(value)) ) %>% dplyr::filter(!is.na(value)) long } plot_theme <- theme_bw(base_family = "Times New Roman", base_size = 11) + theme( strip.background = element_rect(fill = "grey92", color = NA), strip.text = element_text(face = "bold", size = 10), panel.grid.minor = element_blank(), panel.grid.major = element_line(color = "grey90", linewidth = 0.3), legend.position = "bottom", legend.title = element_blank(), axis.title = element_text(size = 11), plot.title = element_text(face = "bold", hjust = 0.5) ) # ---- Variable sets and wave labels wave_labels <- c("Baseline", "Week 4", "Week 8", "Week 12", "Week 18", "Week 24") vars_item1 <- c("VACSS_1", "VACSS_1_4", "VACSS_1_8", "VACSS_1_12", "VACSS_1_18", "VACSS_1_24") vars_item2 <- c("VACSS_2r", "VACSS_2r_4", "VACSS_2r_8", "VACSS_2r_12", "VACSS_2r_18", "VACSS_2r_24") vars_item6 <- c("VACSS_6r", "VACSS_6r_4", "VACSS_6r_8", "VACSS_6r_12", "VACSS_6r_18", "VACSS_6r_24") make_bar_plot <- function(mids_obj, vars, wave_labels, title) { dat <- prep_imp_data(mids_obj, vars, wave_labels) dat_props <- dat %>% dplyr::group_by(source, .imp, variable, value) %>% dplyr::summarize(n = dplyr::n(), .groups = "drop") %>% dplyr::group_by(source, .imp, variable) %>% dplyr::mutate(prop = n / sum(n)) %>% dplyr::ungroup() %>% dplyr::group_by(source, variable, value) %>% dplyr::summarize(prop = mean(prop), .groups = "drop") %>% dplyr::mutate(value = factor(value)) ggplot(dat_props, aes(x = value, y = prop, fill = source)) + geom_col(position = position_dodge(width = 0.8), width = 0.7) + scale_fill_manual(values = c("Observed" = "#0072B2", "Imputed" = "#D55E00")) + scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + facet_wrap(~ variable, ncol = 3) + labs(title = title, x = "Response category", y = "Proportion") + plot_theme } p1_bar <- make_bar_plot(imp_final, vars_item1, wave_labels, "C-SSRS-SV Item 1 (Wish to be Dead)") p2_bar <- make_bar_plot(imp_final, vars_item2, wave_labels, "C-SSRS-SV Recoded Active Suicidal Ideation") p6_bar <- make_bar_plot(imp_final, vars_item6, wave_labels, "C-SSRS-SV Recoded Suicidal Behavior")