Replicating Online Tutorial
ram-online-tutorial.Rmd
library(minvariance)
library(tidyverse)
#> ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
#> ✓ ggplot2 3.3.3 ✓ purrr 0.3.4
#> ✓ tibble 3.1.0 ✓ dplyr 1.0.5
#> ✓ tidyr 1.1.3 ✓ stringr 1.4.0
#> ✓ readr 1.4.0 ✓ forcats 0.5.1
#> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
#> x dplyr::filter() masks stats::filter()
#> x dplyr::lag() masks stats::lag()
library(lavaan)
#> This is lavaan 0.6-8
#> lavaan is FREE software! Please report any bugs.
library(semPlot)Load data
filepath <- "https://quantdev.ssri.psu.edu/sites/qdev/files/WISC_MIexample.csv"
df <- read_csv(file = url(filepath))
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> id = col_double(),
#> info1 = col_double(),
#> comp1 = col_double(),
#> simi1 = col_double(),
#> voca1 = col_double(),
#> info6 = col_double(),
#> comp6 = col_double(),
#> simi6 = col_double(),
#> voca6 = col_double()
#> )Remove correlated uniqueness
# Correlated uniqueness can be removed using the "remove" argument
configural <- long_minvariance_syntax(var_list = list(t1 = c("info1", "comp1", "simi1", "voca1"),
t6 = c("info6", "comp6", "simi6", "voca6")),
model = "configural",
remove = list(unique_covar = TRUE))
#> Configural Invariance Model (Pattern Invariance)
fit_configural <- cfa(configural,
data = df)
summary(fit_configural, fit.measures = TRUE)
#> lavaan 0.6-8 ended normally after 75 iterations
#>
#> Estimator ML
#> Optimization method NLMINB
#> Number of model parameters 27
#> Number of equality constraints 2
#>
#> Number of observations 204
#>
#> Model Test User Model:
#>
#> Test statistic 25.968
#> Degrees of freedom 19
#> P-value (Chi-square) 0.131
#>
#> Model Test Baseline Model:
#>
#> Test statistic 847.740
#> Degrees of freedom 28
#> P-value 0.000
#>
#> User Model versus Baseline Model:
#>
#> Comparative Fit Index (CFI) 0.991
#> Tucker-Lewis Index (TLI) 0.987
#>
#> Loglikelihood and Information Criteria:
#>
#> Loglikelihood user model (H0) -5601.115
#> Loglikelihood unrestricted model (H1) -5588.131
#>
#> Akaike (AIC) 11252.229
#> Bayesian (BIC) 11335.182
#> Sample-size adjusted Bayesian (BIC) 11255.975
#>
#> Root Mean Square Error of Approximation:
#>
#> RMSEA 0.042
#> 90 Percent confidence interval - lower 0.000
#> 90 Percent confidence interval - upper 0.079
#> P-value RMSEA <= 0.05 0.588
#>
#> Standardized Root Mean Square Residual:
#>
#> SRMR 0.031
#>
#> Parameter Estimates:
#>
#> Standard errors Standard
#> Information Expected
#> Information saturated (h1) model Structured
#>
#> Latent Variables:
#> Estimate Std.Err z-value P(>|z|)
#> eta1 =~
#> info1 (lmb1) 4.451 0.396 11.244 0.000
#> comp1 6.850 0.637 10.750 0.000
#> simi1 4.590 0.515 8.918 0.000
#> voca1 5.039 0.393 12.809 0.000
#> eta2 =~
#> info6 (lmb1) 4.451 0.396 11.244 0.000
#> comp6 4.006 0.485 8.261 0.000
#> simi6 4.551 0.546 8.342 0.000
#> voca6 4.102 0.448 9.150 0.000
#>
#> Covariances:
#> Estimate Std.Err z-value P(>|z|)
#> eta1 ~~
#> eta2 1.837 0.211 8.703 0.000
#>
#> Intercepts:
#> Estimate Std.Err z-value P(>|z|)
#> eta1 0.000
#> eta2 6.455 0.601 10.743 0.000
#> .info1 (tau1) 19.776 0.427 46.273 0.000
#> .comp1 21.797 0.680 32.036 0.000
#> .simi1 14.903 0.528 28.223 0.000
#> .voca1 20.396 0.439 46.416 0.000
#> .info6 (tau1) 19.776 0.427 46.273 0.000
#> .comp6 19.317 2.281 8.467 0.000
#> .simi6 11.922 2.538 4.697 0.000
#> .voca6 17.970 1.815 9.903 0.000
#>
#> Variances:
#> Estimate Std.Err z-value P(>|z|)
#> eta1 1.000
#> eta2 5.834 1.159 5.035 0.000
#> .info1 17.448 2.186 7.981 0.000
#> .comp1 47.511 5.748 8.266 0.000
#> .simi1 35.810 3.969 9.022 0.000
#> .voca1 13.999 2.085 6.712 0.000
#> .info6 47.096 6.305 7.470 0.000
#> .comp6 73.850 8.395 8.797 0.000
#> .simi6 88.920 10.222 8.699 0.000
#> .voca6 23.267 4.076 5.709 0.000
semPaths(fit_configural,
what="path",
whatLabels = "name",
sizeInt = 4, sizeMan = 5, sizeLat = 6)
Includeing correlated uniqueness
configural_unique <- long_minvariance_syntax(var_list = list(t1 = c("info1", "comp1", "simi1", "voca1"),
t6 = c("info6", "comp6", "simi6", "voca6")),
model = "configural")
#> Configural Invariance Model (Pattern Invariance)
fit_configural_unique <- cfa(configural_unique,
data = df)
summary(fit_configural_unique,
fit.measures = TRUE)
#> lavaan 0.6-8 ended normally after 111 iterations
#>
#> Estimator ML
#> Optimization method NLMINB
#> Number of model parameters 31
#> Number of equality constraints 2
#>
#> Number of observations 204
#>
#> Model Test User Model:
#>
#> Test statistic 24.882
#> Degrees of freedom 15
#> P-value (Chi-square) 0.052
#>
#> Model Test Baseline Model:
#>
#> Test statistic 847.740
#> Degrees of freedom 28
#> P-value 0.000
#>
#> User Model versus Baseline Model:
#>
#> Comparative Fit Index (CFI) 0.988
#> Tucker-Lewis Index (TLI) 0.977
#>
#> Loglikelihood and Information Criteria:
#>
#> Loglikelihood user model (H0) -5600.572
#> Loglikelihood unrestricted model (H1) -5588.131
#>
#> Akaike (AIC) 11259.143
#> Bayesian (BIC) 11355.369
#> Sample-size adjusted Bayesian (BIC) 11263.488
#>
#> Root Mean Square Error of Approximation:
#>
#> RMSEA 0.057
#> 90 Percent confidence interval - lower 0.000
#> 90 Percent confidence interval - upper 0.095
#> P-value RMSEA <= 0.05 0.350
#>
#> Standardized Root Mean Square Residual:
#>
#> SRMR 0.030
#>
#> Parameter Estimates:
#>
#> Standard errors Standard
#> Information Expected
#> Information saturated (h1) model Structured
#>
#> Latent Variables:
#> Estimate Std.Err z-value P(>|z|)
#> eta1 =~
#> info1 (lmb1) 4.470 0.397 11.253 0.000
#> comp1 6.868 0.639 10.747 0.000
#> simi1 4.636 0.516 8.990 0.000
#> voca1 5.000 0.396 12.618 0.000
#> eta2 =~
#> info6 (lmb1) 4.470 0.397 11.253 0.000
#> comp6 4.026 0.483 8.334 0.000
#> simi6 4.564 0.543 8.406 0.000
#> voca6 4.090 0.442 9.244 0.000
#>
#> Covariances:
#> Estimate Std.Err z-value P(>|z|)
#> eta1 ~~
#> eta2 1.817 0.208 8.715 0.000
#> .info1 ~~
#> .info6 0.979 2.568 0.381 0.703
#> .comp1 ~~
#> .comp6 -0.222 4.821 -0.046 0.963
#> .simi1 ~~
#> .simi6 -1.445 4.436 -0.326 0.745
#> .voca1 ~~
#> .voca6 1.741 2.025 0.860 0.390
#>
#> Intercepts:
#> Estimate Std.Err z-value P(>|z|)
#> eta1 0.000
#> eta2 6.429 0.598 10.757 0.000
#> .info1 (tau1) 19.776 0.427 46.267 0.000
#> .comp1 21.797 0.680 32.039 0.000
#> .simi1 14.903 0.528 28.207 0.000
#> .voca1 20.396 0.439 46.428 0.000
#> .info6 (tau1) 19.776 0.427 46.267 0.000
#> .comp6 19.294 2.282 8.453 0.000
#> .simi6 11.957 2.540 4.708 0.000
#> .voca6 18.151 1.821 9.970 0.000
#>
#> Variances:
#> Estimate Std.Err z-value P(>|z|)
#> eta1 1.000
#> eta2 5.797 1.139 5.089 0.000
#> .info1 17.294 2.205 7.845 0.000
#> .comp1 47.245 5.782 8.171 0.000
#> .simi1 35.454 3.964 8.944 0.000
#> .voca1 14.374 2.142 6.710 0.000
#> .info6 46.958 6.398 7.339 0.000
#> .comp6 73.488 8.419 8.729 0.000
#> .simi6 88.773 10.271 8.643 0.000
#> .voca6 23.859 4.191 5.692 0.000
semPaths(fit_configural_unique,
what="path",
whatLabels = "name",
sizeInt = 4, sizeMan = 5, sizeLat = 6)
Fit statistics
minvariance(data = df,
var_list = list(t1 = c("info1", "comp1", "simi1", "voca1"),
t6 = c("info6", "comp6", "simi6", "voca6")),
remove = list(unique_covar = T),
return = "fit_statistics") %>%
knitr::kable(digits = 3)| model | chisq | npar | aic | bic | cfi | rmsea | srmr | tli | converged |
|---|---|---|---|---|---|---|---|---|---|
| configural | 25.968 | 25 | 11252.23 | 11335.18 | 0.991 | 0.042 | 0.031 | 0.987 | TRUE |
| weak | 41.897 | 22 | 11262.16 | 11335.16 | 0.976 | 0.067 | 0.076 | 0.969 | TRUE |
| strong | 53.723 | 19 | 11267.98 | 11331.03 | 0.965 | 0.075 | 0.087 | 0.961 | TRUE |
| strict | 134.559 | 15 | 11340.82 | 11390.59 | 0.871 | 0.134 | 0.169 | 0.876 | TRUE |
Model tests
minvariance(data = df,
var_list = list(t1 = c("info1", "comp1", "simi1", "voca1"),
t6 = c("info6", "comp6", "simi6", "voca6")),
remove = list(unique_covar = T),
return = "model_tests") %>%
knitr::kable(digits = 3)| Df | AIC | BIC | Chisq | Chisq diff | Df diff | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|
| configural | 19 | 11252.23 | 11335.18 | 25.968 | NA | NA | NA |
| weak | 22 | 11262.16 | 11335.16 | 41.897 | 15.929 | 3 | 0.001 |
| strong | 25 | 11267.98 | 11331.03 | 53.723 | 11.826 | 3 | 0.008 |
| strict | 29 | 11340.82 | 11390.59 | 134.559 | 80.836 | 4 | 0.000 |
lavaan syntax
lminvar_syntax <- minvariance(data = df,
var_list = list(t1 = c("info1", "comp1", "simi1", "voca1"),
t6 = c("info6", "comp6", "simi6", "voca6")),
remove = list(unique_covar = T),
return = "lavaan_syntax") Configural
lminvar_syntax$configural %>%
cat()
#> #### CONFIGURAL INVARIANCE MODEL ####
#> # Specify Latent Factors ----
#> eta1 =~ NA * info1 + lambda1 * info1 + comp1 + simi1 + voca1
#> eta2 =~ NA * info6 + lambda1 * info6 + comp6 + simi6 + voca6
#> # Specify Latent Variable Means ----
#> eta1 ~ 0 * 1
#> eta2 ~ 1
#> # Specify Latent Variable Variances ----
#> eta1 ~~ 1 * eta1
#> eta2 ~~ eta2
#> # Specify Latent Variable Covariances ----
#> eta1 ~~ eta2
#> # Specify Observed Variable Intercepts ----
#> info1 ~ tau1 * 1
#> comp1 ~ 1
#> simi1 ~ 1
#> voca1 ~ 1
#> info6 ~ tau1 * 1
#> comp6 ~ 1
#> simi6 ~ 1
#> voca6 ~ 1
#> # Specify Unique Variances ----
#> info1 ~~ info1
#> comp1 ~~ comp1
#> simi1 ~~ simi1
#> voca1 ~~ voca1
#> info6 ~~ info6
#> comp6 ~~ comp6
#> simi6 ~~ simi6
#> voca6 ~~ voca6
#> # Specify Unique Covariances ----
#> # REMOVEDWeak
lminvar_syntax$weak %>%
cat()
#> #### WEAK INVARIANCE MODEL ####
#> # Specify Latent Factors ----
#> eta1 =~ NA * info1 + lambda1 * info1 + lambda2 * comp1 + lambda3 * simi1 + lambda4 * voca1
#> eta2 =~ NA * info6 + lambda1 * info6 + lambda2 * comp6 + lambda3 * simi6 + lambda4 * voca6
#> # Specify Latent Variable Means ----
#> eta1 ~ 0 * 1
#> eta2 ~ 1
#> # Specify Latent Variable Variances ----
#> eta1 ~~ 1 * eta1
#> eta2 ~~ eta2
#> # Specify Latent Variable Covariances ----
#> eta1 ~~ eta2
#> # Specify Observed Variable Intercepts ----
#> info1 ~ tau1 * 1
#> comp1 ~ 1
#> simi1 ~ 1
#> voca1 ~ 1
#> info6 ~ tau1 * 1
#> comp6 ~ 1
#> simi6 ~ 1
#> voca6 ~ 1
#> # Specify Unique Variances ----
#> info1 ~~ info1
#> comp1 ~~ comp1
#> simi1 ~~ simi1
#> voca1 ~~ voca1
#> info6 ~~ info6
#> comp6 ~~ comp6
#> simi6 ~~ simi6
#> voca6 ~~ voca6
#> # Specify Unique Covariances ----
#> # REMOVEDStrong
lminvar_syntax$strong %>%
cat()
#> #### STRONG INVARIANCE MODEL ####
#> # Specify Latent Factors ----
#> eta1 =~ NA * info1 + lambda1 * info1 + lambda2 * comp1 + lambda3 * simi1 + lambda4 * voca1
#> eta2 =~ NA * info6 + lambda1 * info6 + lambda2 * comp6 + lambda3 * simi6 + lambda4 * voca6
#> # Specify Latent Variable Means ----
#> eta1 ~ 0 * 1
#> eta2 ~ 1
#> # Specify Latent Variable Variances ----
#> eta1 ~~ 1 * eta1
#> eta2 ~~ eta2
#> # Specify Latent Variable Covariances ----
#> eta1 ~~ eta2
#> # Specify Observed Variable Intercepts ----
#> info1 ~ tau1 * 1
#> comp1 ~ tau2 * 1
#> simi1 ~ tau3 * 1
#> voca1 ~ tau4 * 1
#> info6 ~ tau1 * 1
#> comp6 ~ tau2 * 1
#> simi6 ~ tau3 * 1
#> voca6 ~ tau4 * 1
#> # Specify Unique Variances ----
#> info1 ~~ info1
#> comp1 ~~ comp1
#> simi1 ~~ simi1
#> voca1 ~~ voca1
#> info6 ~~ info6
#> comp6 ~~ comp6
#> simi6 ~~ simi6
#> voca6 ~~ voca6
#> # Specify Unique Covariances ----
#> # REMOVEDStrict
lminvar_syntax$strict %>%
cat()
#> #### STRICT INVARIANCE MODEL ####
#> # Specify Latent Factors ----
#> eta1 =~ NA * info1 + lambda1 * info1 + lambda2 * comp1 + lambda3 * simi1 + lambda4 * voca1
#> eta2 =~ NA * info6 + lambda1 * info6 + lambda2 * comp6 + lambda3 * simi6 + lambda4 * voca6
#> # Specify Latent Variable Means ----
#> eta1 ~ 0 * 1
#> eta2 ~ 1
#> # Specify Latent Variable Variances ----
#> eta1 ~~ 1 * eta1
#> eta2 ~~ eta2
#> # Specify Latent Variable Covariances ----
#> eta1 ~~ eta2
#> # Specify Observed Variable Intercepts ----
#> info1 ~ tau1 * 1
#> comp1 ~ tau2 * 1
#> simi1 ~ tau3 * 1
#> voca1 ~ tau4 * 1
#> info6 ~ tau1 * 1
#> comp6 ~ tau2 * 1
#> simi6 ~ tau3 * 1
#> voca6 ~ tau4 * 1
#> # Specify Unique Variances ----
#> info1 ~~ theta1 * info1
#> comp1 ~~ theta2 * comp1
#> simi1 ~~ theta3 * simi1
#> voca1 ~~ theta4 * voca1
#> info6 ~~ theta1 * info6
#> comp6 ~~ theta2 * comp6
#> simi6 ~~ theta3 * simi6
#> voca6 ~~ theta4 * voca6
#> # Specify Unique Covariances ----
#> # REMOVED