Last updated: 2020-06-10

Checks: 6 1

Knit directory: mcfa-para-est/

This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.


The R Markdown file has unstaged changes. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish to commit the R Markdown file and build the HTML.

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

The command set.seed(20190614) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version eecb366. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    data/compiled_para_results.txt
    Ignored:    data/results_bias_est.csv
    Ignored:    data/results_bias_se.csv
    Ignored:    fig/
    Ignored:    manuscript/
    Ignored:    output/fact-cov-converge-largeN.pdf
    Ignored:    output/fact-cov-converge-medN.pdf
    Ignored:    output/fact-cov-converge-smallN.pdf
    Ignored:    output/loading-converge-largeN.pdf
    Ignored:    output/loading-converge-medN.pdf
    Ignored:    output/loading-converge-smallN.pdf
    Ignored:    references/
    Ignored:    sera-presentation/

Untracked files:
    Untracked:  analysis/ml-cfa-parameter-anova-estimates.Rmd
    Untracked:  analysis/ml-cfa-parameter-anova-relative-bias.Rmd
    Untracked:  analysis/ml-cfa-parameter-bias-latent-ICC.Rmd
    Untracked:  analysis/ml-cfa-parameter-bias-observed-ICC.Rmd
    Untracked:  analysis/ml-cfa-parameter-bias-pub-figure.Rmd
    Untracked:  analysis/ml-cfa-parameter-convergence-ARD-L1-factor-covariance.Rmd
    Untracked:  analysis/ml-cfa-parameter-convergence-ARD-L2-factor-covariance.Rmd
    Untracked:  analysis/ml-cfa-parameter-convergence-ARD-L2-factor-variance.Rmd
    Untracked:  analysis/ml-cfa-parameter-convergence-ARD-L2-residual-variance.Rmd
    Untracked:  analysis/ml-cfa-parameter-convergence-ARD-factor-loadings.Rmd
    Untracked:  analysis/ml-cfa-parameter-convergence-ARD-latent-ICC.Rmd
    Untracked:  analysis/ml-cfa-parameter-convergence-ARD-observed-ICC.Rmd
    Untracked:  analysis/ml-cfa-parameter-convergence-correlation-pubfigure.Rmd
    Untracked:  analysis/ml-cfa-parameter-convergence-trace-plots-factor-loadings.Rmd
    Untracked:  analysis/ml-cfa-standard-error-anova-logSE.Rmd
    Untracked:  analysis/ml-cfa-standard-error-anova-relative-bias.Rmd
    Untracked:  analysis/ml-cfa-standard-error-bias-factor-loadings.Rmd
    Untracked:  analysis/ml-cfa-standard-error-bias-level1-factor-covariances.Rmd
    Untracked:  analysis/ml-cfa-standard-error-bias-level2-factor-covariances.Rmd
    Untracked:  analysis/ml-cfa-standard-error-bias-level2-factor-variances.Rmd
    Untracked:  analysis/ml-cfa-standard-error-bias-level2-residual-variances.Rmd
    Untracked:  analysis/ml-cfa-standard-error-bias-overview.Rmd
    Untracked:  code/r_functions.R
    Untracked:  renv.lock
    Untracked:  renv/

Unstaged changes:
    Modified:   .gitignore
    Modified:   analysis/index.Rmd
    Modified:   analysis/ml-cfa-convergence-summary.Rmd
    Modified:   analysis/ml-cfa-parameter-bias-factor-loadings.Rmd
    Modified:   analysis/ml-cfa-parameter-bias-level1-factor-covariances.Rmd
    Modified:   analysis/ml-cfa-parameter-bias-level2-factor-covariances.Rmd
    Modified:   analysis/ml-cfa-parameter-bias-level2-factor-variances.Rmd
    Modified:   analysis/ml-cfa-parameter-bias-level2-residual-variances.Rmd
    Modified:   analysis/ml-cfa-parameter-convergence-correlation-factor-loadings.Rmd
    Modified:   code/get_data.R
    Modified:   code/load_packages.R

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.


These are the previous versions of the repository in which changes were made to the R Markdown (analysis/ml-cfa-parameter-bias-level2-factor-variances.Rmd) and HTML (docs/ml-cfa-parameter-bias-level2-factor-variances.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File Version Author Date Message
html 95dd5a6 noah-padgett 2020-03-31 Build site.
Rmd 9f8001a noah-padgett 2020-03-31 fixed figure size
html 9f8001a noah-padgett 2020-03-31 fixed figure size
html 51afcfa noah-padgett 2020-03-31 Build site.
Rmd d5ffe1d noah-padgett 2020-03-31 updated bias pages
html d5ffe1d noah-padgett 2020-03-31 updated bias pages

rm(list=ls())
source(paste0(getwd(),"/code/load_packages.R"))
#source(paste0(getwd(),"/code/get_data.R"))
sessionInfo()
R version 3.6.3 (2020-02-29)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18362)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252 
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] xtable_1.8-4          kableExtra_1.1.0      cowplot_1.0.0        
 [4] MplusAutomation_0.7-3 data.table_1.12.8     patchwork_1.0.0      
 [7] forcats_0.5.0         stringr_1.4.0         dplyr_0.8.5          
[10] purrr_0.3.4           readr_1.3.1           tidyr_1.1.0          
[13] tibble_3.0.1          ggplot2_3.3.0         tidyverse_1.3.0      
[16] workflowr_1.6.2      

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.4.6      lubridate_1.7.8   lattice_0.20-38   assertthat_0.2.1 
 [5] rprojroot_1.3-2   digest_0.6.25     R6_2.4.1          cellranger_1.1.0 
 [9] plyr_1.8.6        backports_1.1.7   reprex_0.3.0      evaluate_0.14    
[13] coda_0.19-3       httr_1.4.1        pillar_1.4.4      rlang_0.4.6      
[17] readxl_1.3.1      rstudioapi_0.11   whisker_0.4       blob_1.2.1       
[21] texreg_1.36.23    gsubfn_0.7        rmarkdown_2.1     proto_1.0.0      
[25] webshot_0.5.2     pander_0.6.3      munsell_0.5.0     broom_0.5.6      
[29] compiler_3.6.3    httpuv_1.5.2      modelr_0.1.8      xfun_0.14        
[33] pkgconfig_2.0.3   htmltools_0.4.0   tidyselect_1.1.0  viridisLite_0.3.0
[37] fansi_0.4.1       crayon_1.3.4      dbplyr_1.4.4      withr_2.2.0      
[41] later_1.0.0       grid_3.6.3        nlme_3.1-144      jsonlite_1.6.1   
[45] gtable_0.3.0      lifecycle_0.2.0   DBI_1.1.0         git2r_0.27.1     
[49] magrittr_1.5      scales_1.1.1      cli_2.0.2         stringi_1.4.6    
[53] fs_1.4.1          promises_1.1.0    xml2_1.3.2        ellipsis_0.3.1   
[57] generics_0.0.2    vctrs_0.3.0       boot_1.3-24       tools_3.6.3      
[61] glue_1.4.1        hms_0.5.3         parallel_3.6.3    yaml_2.2.1       
[65] colorspace_1.4-1  rvest_0.3.5       knitr_1.28        haven_2.3.0      
# general options
theme_set(theme_bw())
options(digits=3)
# set up vectors of variable names
pvec <- c(paste0('lambda1',1:6), paste0('lambda2',6:10), 'psiW12','psiB1', 'psiB2', 'psiB12', paste0('thetaB',1:10), 'icc_lv1_est', 'icc_lv2_est', paste0('icc_ov',1:10,'_est'))
# stored "true" values of parameters by each condition
ptvec <- c(rep('lambdaT',11), 'psiW12T', 'psiB1T', 'psiB2T', 'psiB12T', rep("thetaBT", 10), rep('icc_lv',2), rep('icc_ov',10))

result <- read_csv(paste0(w.d, "/data/results_bias_est.csv"))
Parsed with column specification:
cols(
  N1 = col_double(),
  N2 = col_double(),
  ICC_LV = col_double(),
  ICC_OV = col_double(),
  Variable = col_character(),
  Estimator = col_character(),
  TrueValue = col_double(),
  RB = col_double(),
  RMSE = col_double(),
  Bias = col_double(),
  SampVar = col_double(),
  muRE = col_double(),
  mwRE = col_double(),
  uwRE = col_double(),
  nRep = col_double(),
  estMean = col_double(),
  estSD = col_double()
)
# Set conditions levels as categorical values
result <- result %>%
  mutate(N1 = factor(N1, c("5", "10", "30")),
         N2 = factor(N2, c("30", "50", "100", "200")),
         ICC_OV = factor(ICC_OV, c("0.1","0.3", "0.5")),
         ICC_LV = factor(ICC_LV, c("0.1", "0.5")),
         wi = nRep/500)

Summarizing Results

First, we will plot estimates (botxplots) to show how these estimates changed across conditions. To summarize the results we will average over the parameters that only differ y indices. Meaning we will describe the “average factor loading bias” by reporting the average bias for factor loadings. Additionally, different conditions resultedin different “sample sizes.” By this we mean the number of uses replications. The different number of cases per condition was accounted for by creating a “weight” variable for each row of the result object. This meant that conditions that had more usable replications counted more towards to averages reported (or count as much as if we averaged over the individual replications).

*Click here for more details

Level-2 Factor Variance

sdat <- filter(result, Variable %in% c("psiB1","psiB2"))
sdat <- sdat %>%
  mutate(TrueValue = factor(TrueValue))

TRUEVALUE <- as.numeric(levels(sdat$TrueValue))

p1 <- ggplot(sdat, aes(y=estMean, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  geom_hline(yintercept = TRUEVALUE, color="red")+
  labs(y="Average Level-2 Factor Covariance")

p2 <- ggplot(sdat, aes(y=estSD, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="SD of Level-2 Factor Covariance")

p3 <- ggplot(sdat, aes(y=RB, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  geom_hline(yintercept=-10, color="red", linetype="dashed")+
  geom_hline(yintercept=10, color="red", linetype="dashed")+
  labs(y="Relative Bias")

p4 <- ggplot(sdat, aes(y=RMSE, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Root Mean Square Error")

p5 <- ggplot(sdat, aes(y=Bias, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Sqaured Bias")

p6 <- ggplot(sdat, aes(y=SampVar, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Sampling Variance of Estimates")


p <- (p1 + p2 + p3)/(p4 + p5 + p6) + 
  plot_annotation(title="Summarizing bias indices of LEVEL-2 FACTOR VARIANCE")
p

Single Condition Breakdown

Estimation Method

ggplot(sdat, aes(y=estMean, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  geom_hline(yintercept = TRUEVALUE, color="red")+
  labs(y="Average Level-2 Factor Covariance",
       title="LEVEL-2 FACTOR COVARIANCE by Estimation Method",
       subtitle="Parameter Estimates")+
  facet_wrap(.~Estimator)

ggplot(sdat, aes(y=estSD, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="SD of Level-2 Factor Covariances",
       title="LEVEL-2 FACTOR COVARIANCE by Estimation Method",
       subtitle="Standard Deviation of Estimates")+
  facet_wrap(.~Estimator)

ggplot(sdat, aes(y=RB, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  geom_hline(yintercept=-10, color="red", linetype="dashed")+
  geom_hline(yintercept=10, color="red", linetype="dashed")+
  labs(y="Relative Bias",
       title="LEVEL-2 FACTOR COVARIANCE by Estimation Method",
       subtitle="Relative Bias of Estimates")+
  facet_wrap(.~Estimator)

ggplot(sdat, aes(y=RMSE, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Root Mean Square Error",
       title="LEVEL-2 FACTOR COVARIANCE by Estimation Method",
       subtitle="Root Mean Square Error of Estimates")+
  facet_wrap(.~Estimator)

ggplot(sdat, aes(y=Bias, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Sqaured Bias",
       title="LEVEL-2 FACTOR COVARIANCE by Estimation Method",
       subtitle="Squared Bias of Estiamtes")+
  facet_wrap(.~Estimator)

ggplot(sdat, aes(y=SampVar, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Sampling Variance",
       title="LEVEL-2 FACTOR COVARIANCE by Estimation Method",
       subtitle="Sampling Variance of Estimates")+
  facet_wrap(.~Estimator)

c <- sdat %>%
  group_by(Estimator, TrueValue) %>%
  summarise(est = weighted.mean(estMean, wi),
            RB = weighted.mean(RB, wi),
            RMSE = weighted.mean(RMSE, wi),
            Bias = weighted.mean(Bias, wi),
            SampVar =weighted.mean(SampVar, wi))

kable(c, format='html', digits=3,
      caption="Summary Indices of LEVEL-2 FACTOR COVARIANCE by Estimation Method") %>%
  kable_styling(full_width = T)
Summary Indices of LEVEL-2 FACTOR COVARIANCE by Estimation Method
Estimator TrueValue est RB RMSE Bias SampVar
MLR 0.111 0.135 21.952 0.011 0.002 0.009
MLR 1 0.999 -0.124 0.099 0.000 0.099
ULSMV 0.111 0.138 24.341 0.015 0.003 0.012
ULSMV 1 1.035 3.463 0.112 0.003 0.109
WLSMV 0.111 0.119 6.781 0.007 0.001 0.006
WLSMV 1 0.940 -6.012 0.095 0.007 0.089

Level-2 Sample Size

ggplot(sdat, aes(y=estMean, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  geom_hline(yintercept = TRUEVALUE, color="red")+
  labs(y="Average Level-2 Factor Covariance",
       title="LEVEL-2 FACTOR COVARIANCE by Level-2 Sample Size",
       subtitle="Parameter Estimates")+
  facet_wrap(.~N2)

ggplot(sdat, aes(y=estSD, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="SD of Level-2 Factor Covariances",
       title="LEVEL-2 FACTOR COVARIANCE by Level-2 Sample Size",
       subtitle="Standard Deviation of Parameter Estimates")+
  facet_wrap(.~N2)

ggplot(sdat, aes(y=RB, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  geom_hline(yintercept=-10, color="red", linetype="dashed")+
  geom_hline(yintercept=10, color="red", linetype="dashed")+
  labs(y="Relative Bias",
       title="LEVEL-2 FACTOR COVARIANCE by Level-2 Sample Size",
       subtitle="Relative Bias Parameter Estimates")+
  facet_wrap(.~N2)

ggplot(sdat, aes(y=RMSE, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Root Mean Square Error",
       title="LEVEL-2 FACTOR COVARIANCE by Level-2 Sample Size",
       subtitle="Root Mean Square Error")+
  facet_wrap(.~N2)

ggplot(sdat, aes(y=Bias, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Sqaured Bias",
       title="LEVEL-2 FACTOR COVARIANCE by Level-2 Sample Size",
       subtitle="Squared Bias of Parameter Estimates")+
  facet_wrap(.~N2)

ggplot(sdat, aes(y=SampVar, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Sampling Variance of Estimates",
       title="LEVEL-2 FACTOR COVARIANCE by Level-2 Sample Size",
       subtitle="Sampling Variance of Parameter Estimates")+
  facet_wrap(.~N2)

c <- sdat %>%
  group_by(N2, TrueValue) %>%
  summarise(est = weighted.mean(estMean, wi),
            RB = weighted.mean(RB, wi),
            RMSE = weighted.mean(RMSE, wi),
            Bias = weighted.mean(Bias, wi),
            SampVar =weighted.mean(SampVar, wi))

kable(c, format='html', digits=3, 
      caption="Summary Indices of LEVEL-2 FACTOR COVARIANCE by Level-2 Sample Size") %>%
  kable_styling(full_width = T)
Summary Indices of LEVEL-2 FACTOR COVARIANCE by Level-2 Sample Size
N2 TrueValue est RB RMSE Bias SampVar
30 0.111 0.161 44.650 0.030 0.007 0.023
30 1 0.987 -1.316 0.228 0.010 0.218
50 0.111 0.140 26.115 0.013 0.002 0.011
50 1 0.986 -1.435 0.126 0.003 0.123
100 0.111 0.124 11.862 0.006 0.001 0.006
100 1 0.995 -0.480 0.059 0.001 0.058
200 0.111 0.115 3.720 0.004 0.000 0.003
200 1 1.000 -0.031 0.028 0.000 0.028

Level-1 Sample Size

ggplot(sdat, aes(y=estMean, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  geom_hline(yintercept = TRUEVALUE, color="red")+
  labs(y="Average Level-2 Factor Covariance",
       title="LEVEL-2 FACTOR COVARIANCE by Level-1",
       subtitle="Parameter Estimates")+
  facet_wrap(.~N1)

ggplot(sdat, aes(y=estSD, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="SD of Level-2 Factor Covariances",
       title="LEVEL-2 FACTOR COVARIANCE by Level-1 Sample Size",
       subtitle="Standard Deviation of Parameter Estimates")+
  facet_wrap(.~N1)

ggplot(sdat, aes(y=RB, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  geom_hline(yintercept=-10, color="red", linetype="dashed")+
  geom_hline(yintercept=10, color="red", linetype="dashed")+
  labs(y="Relative Bias",
       title="LEVEL-2 FACTOR COVARIANCE by Level-1 Sample Size",
       subtitle="Relative Bias of Parameter Estimates")+
  facet_wrap(.~N1)

ggplot(sdat, aes(y=RMSE, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Root Mean Square Error",
       title="LEVEL-2 FACTOR COVARIANCE by Level-1 Sample Size",
       subtitle="Root Mean Square Error")+
  facet_wrap(.~N1)

ggplot(sdat, aes(y=Bias, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Sqaured Bias",
       title="LEVEL-2 FACTOR COVARIANCE by Level-1 Sample Size",
       subtitle="Squared Bias of Parameter Estimates")+
  facet_wrap(.~N1)

ggplot(sdat, aes(y=SampVar, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Sampling Variance of Estimates",
       title="LEVEL-2 FACTOR COVARIANCE by Level-1 Sample Size",
       subtitle="Sampling Variance of Parameter Estimates")+
  facet_wrap(.~N1)

c <- sdat %>%
  group_by(N1, TrueValue) %>%
  summarise(est = weighted.mean(estMean, wi),
            RB = weighted.mean(RB, wi),
            RMSE = weighted.mean(RMSE, wi),
            Bias = weighted.mean(Bias, wi),
            SampVar =weighted.mean(SampVar, wi))

kable(c, format='html', digits=3,
      caption="Summary Indices of LEVEL-2 FACTOR COVARIANCE  by Level-1 Sample Size") %>%
  kable_styling(full_width = T)
Summary Indices of LEVEL-2 FACTOR COVARIANCE by Level-1 Sample Size
N1 TrueValue est RB RMSE Bias SampVar
5 0.111 0.149 34.048 0.021 0.004 0.017
5 1 1.000 -0.035 0.155 0.006 0.149
10 0.111 0.129 15.895 0.010 0.001 0.008
10 1 0.988 -1.221 0.092 0.003 0.088
30 0.111 0.121 9.228 0.006 0.001 0.005
30 1 0.991 -0.923 0.070 0.001 0.068

ICC Observed Variables

ggplot(sdat, aes(y=estMean, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  geom_hline(yintercept = TRUEVALUE, color="red")+
  labs(y="Average Level-2 Factor Covariance",
       title="LEVEL-2 FACTOR COVARIANCE by ICC of Observed Variables",
       subtitle="Parameter Estimates")+
  facet_wrap(.~ICC_OV)

ggplot(sdat, aes(y=estSD, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="SD of Level-2 Factor Covariances",
       title="LEVEL-2 FACTOR COVARIANCE by ICC of Observed Variables",
       subtitle="Standard Deviation of Parameter Estimates")+
  facet_wrap(.~ICC_OV)

ggplot(sdat, aes(y=RB, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  geom_hline(yintercept=-10, color="red", linetype="dashed")+
  geom_hline(yintercept=10, color="red", linetype="dashed")+
  labs(y="Relative Bias",
       title="LEVEL-2 FACTOR COVARIANCE by ICC of Observed Variables",
       subtitle="Relative Bias of Parameter Estimates")+
  facet_wrap(.~ICC_OV)

ggplot(sdat, aes(y=RMSE, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Root Mean Square Error",
       title="LEVEL-2 FACTOR COVARIANCE by ICC of Observed Variables",
       subtitle="Root Mean Square Error of Parameter Estimates")+
  facet_wrap(.~ICC_OV)

ggplot(sdat, aes(y=Bias, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Sqaured Bias",
       title="LEVEL-2 FACTOR COVARIANCE by ICC of Observed Variables",
       subtitle="Squared Bias of Parameter Estimates")+
  facet_wrap(.~ICC_OV)

ggplot(sdat, aes(y=SampVar, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Sampling Variance of Estimates",
       title="LEVEL-2 FACTOR COVARIANCE by ICC of Observed Variables",
       subtitle="Sampling Variance of Parameter Estimates")+
  facet_wrap(.~ICC_OV)

c <- sdat %>%
  group_by(ICC_OV, TrueValue) %>%
  summarise(est = weighted.mean(estMean, wi),
            RB = weighted.mean(RB, wi),
            RMSE = weighted.mean(RMSE, wi),
            Bias = weighted.mean(Bias, wi),
            SampVar =weighted.mean(SampVar, wi))

kable(c, format='html', digits=3, caption="Summary Indices of LEVEL-2 FACTOR COVARIANCE by ICC of Observed Variables") %>%
  kable_styling(full_width = T)
Summary Indices of LEVEL-2 FACTOR COVARIANCE by ICC of Observed Variables
ICC_OV TrueValue est RB RMSE Bias SampVar
0.1 0.111 0.108 -2.336 0.003 0.000 0.003
0.1 1 0.989 -1.112 0.050 0.001 0.049
0.3 0.111 0.126 13.120 0.008 0.001 0.008
0.3 1 0.985 -1.500 0.092 0.002 0.089
0.5 0.111 0.167 50.375 0.025 0.006 0.019
0.5 1 1.002 0.175 0.144 0.006 0.138

ICC Latent Variables

ggplot(sdat, aes(y=estMean, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  geom_hline(yintercept = TRUEVALUE[1], color="red")+
  labs(y="Average Level-2 Factor Covariance",
       title="LEVEL-2 FACTOR COVARIANCE by ICC of Latent Variables",
       subtitle="Parameter Estimates")+
  facet_wrap(.~ICC_LV)

ggplot(sdat, aes(y=estSD, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="SD of Level-2 Factor Covariances",
       title="LEVEL-2 FACTOR COVARIANCE by ICC of Latent Variables",
       subtitle="Standard Deviation of Parameter Estimates")+
  facet_wrap(.~ICC_LV)

ggplot(sdat, aes(y=RB, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  geom_hline(yintercept=-10, color="red", linetype="dashed")+
  geom_hline(yintercept=10, color="red", linetype="dashed")+
  labs(y="Relative Bias",
       title="LEVEL-2 FACTOR COVARIANCE by ICC of Latent Variables",
       subtitle="Relative Bias of Parameter Estimates")+
  facet_wrap(.~ICC_LV)

ggplot(sdat, aes(y=RMSE, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Root Mean Square Error",
       title="LEVEL-2 FACTOR COVARIANCE by ICC of Latent Variables",
       subtitle="Root Mean Square Error of Parameter Estimates")+
  facet_wrap(.~ICC_LV)

ggplot(sdat, aes(y=Bias, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Sqaured Bias",
       title="LEVEL-2 FACTOR COVARIANCE by ICC of Latent Variables",
       subtitle="Squared Bias of Parameter Estimates")+
  facet_wrap(.~ICC_LV)

ggplot(sdat, aes(y=SampVar, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Sampling Variance of Estimates",
       title="LEVEL-2 FACTOR COVARIANCE by ICC of Latent Variables",
       subtitle="Sampling Variance of Parameter Estimates")+
  facet_wrap(.~ICC_LV)

c <- sdat %>%
  group_by(ICC_LV, TrueValue) %>%
  summarise(est = weighted.mean(estMean, wi),
            RB = weighted.mean(RB, wi),
            RMSE = weighted.mean(RMSE, wi),
            Bias = weighted.mean(Bias, wi),
            SampVar =weighted.mean(SampVar, wi))

kable(c, format='html', digits=3,
      caption="Summary Indices of LEVEL-2 FACTOR COVARIANCE by ICC of Latent Variables") %>%
  kable_styling(full_width = T)
Summary Indices of LEVEL-2 FACTOR COVARIANCE by ICC of Latent Variables
ICC_LV TrueValue est RB RMSE Bias SampVar
0.1 0.111 0.131 17.968 0.011 0.002 0.009
0.5 1 0.992 -0.759 0.102 0.003 0.099

Loadings by Estimation Method and Sample Sizes

Estimation Method & Level-2 Sample Size

ggplot(sdat, aes(y=estMean, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Average Level-2 Factor Covariance")+
  facet_grid(N2~Estimator)

ggplot(sdat, aes(y=RB, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  geom_hline(yintercept=-10, color="red", linetype="dashed")+
  geom_hline(yintercept=10, color="red", linetype="dashed")+
  labs(y="Relative Bias")+
  facet_grid(N2~Estimator)

ggplot(sdat, aes(y=RMSE, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Root Mean Square Error")+
  facet_grid(N2~Estimator)

c <- sdat %>%
  group_by(Estimator, N2, TrueValue) %>%
  summarise(est = weighted.mean(estMean, wi),
            RB = weighted.mean(RB, wi),
            RMSE = weighted.mean(RMSE, wi),
            Bias = weighted.mean(Bias, wi),
            SampVar =weighted.mean(SampVar, wi))

c1 <- cbind(c[ c$Estimator == 'MLR', c( 'N2','TrueValue', 'est', 'RB', 'RMSE')], 
           c[ c$Estimator == 'ULSMV', c('est', 'RB', 'RMSE')], 
           c[ c$Estimator == 'WLSMV', c('est', 'RB', 'RMSE')])
colnames(c1) <- c('N2','True Value', rep(c('est', 'RB', 'RMSE'), 3))

kable(c1, format='html', digits=3, row.names = F) %>%
  kable_styling(full_width = T) %>%
  add_header_above(c(' '=2, 'MLR'=3, 'ULSMV'=3, 'WLSMV'=3))
MLR
ULSMV
WLSMV
N2 True Value est RB RMSE est RB RMSE est RB RMSE
30 0.111 0.165 48.357 0.026 0.182 64.00 0.049 0.130 17.220 0.015
30 1 1.001 0.081 0.210 1.083 8.35 0.268 0.862 -13.789 0.208
50 0.111 0.145 30.967 0.014 0.149 34.03 0.016 0.124 11.740 0.010
50 1 0.995 -0.515 0.119 1.038 3.77 0.135 0.918 -8.243 0.126
100 0.111 0.128 14.945 0.007 0.128 15.21 0.007 0.117 5.134 0.006
100 1 0.999 -0.118 0.057 1.019 1.89 0.060 0.966 -3.415 0.059
200 0.111 0.117 5.246 0.004 0.117 5.46 0.004 0.111 0.418 0.003
200 1 1.001 0.058 0.028 1.011 1.15 0.029 0.987 -1.338 0.028

Estimation Method & Level-1 Sample Size

ggplot(sdat, aes(y=estMean, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  geom_hline(yintercept = TRUEVALUE, color="red")+
  labs(y="Average Level-2 Factor Covariance")+
  facet_grid(N1~Estimator)

ggplot(sdat, aes(y=RB, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  geom_hline(yintercept=-10, color="red", linetype="dashed")+
  geom_hline(yintercept=10, color="red", linetype="dashed")+
  labs(y="Relative Bias")+
  facet_grid(N1~Estimator)

ggplot(sdat, aes(y=RMSE, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Root Mean Square Error")+
  facet_grid(N1~Estimator)

c <- sdat %>%
  group_by(Estimator, N1, TrueValue) %>%
  summarise(est = weighted.mean(estMean, wi),
            RB = weighted.mean(RB, wi),
            RMSE = weighted.mean(RMSE, wi),
            Bias = weighted.mean(Bias, wi),
            SampVar =weighted.mean(SampVar, wi))

c1 <- cbind(c[ c$Estimator == 'MLR', c( 'N1','TrueValue', 'est', 'RB', 'RMSE')], 
           c[ c$Estimator == 'ULSMV', c('est', 'RB', 'RMSE')], 
           c[ c$Estimator == 'WLSMV', c('est', 'RB', 'RMSE')])
colnames(c1) <- c('N1', "True Value", rep(c('est', 'RB', 'RMSE'), 3))

kable(c1, format='html', digits=3, row.names = F) %>%
  kable_styling(full_width = T) %>%
  add_header_above(c(' '=2, 'MLR'=3, 'ULSMV'=3, 'WLSMV'=3))
MLR
ULSMV
WLSMV
N1 True Value est RB RMSE est RB RMSE est RB RMSE
5 0.111 0.153 37.632 0.019 0.161 44.86 0.033 0.130 17.34 0.011
5 1 1.010 0.990 0.148 1.058 5.78 0.181 0.926 -7.39 0.135
10 0.111 0.133 19.947 0.010 0.135 21.98 0.012 0.117 5.03 0.007
10 1 0.995 -0.495 0.088 1.031 3.14 0.098 0.932 -6.76 0.089
30 0.111 0.124 11.812 0.006 0.126 13.59 0.007 0.113 2.13 0.005
30 1 0.993 -0.733 0.067 1.020 1.96 0.072 0.958 -4.24 0.070

Estimation Method, Level-2 Sample Size & Level-1 Sample Size

ggplot(sdat, aes(y=estMean,, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  geom_hline(yintercept = TRUEVALUE, color="red")+
  labs(y="Average Level-12 Factor Covariance")+
  facet_grid(N2+N1~Estimator)

ggplot(sdat, aes(y=RB, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  geom_hline(yintercept=-10, color="red", linetype="dashed")+
  geom_hline(yintercept=10, color="red", linetype="dashed")+
  labs(y="Relative Bias")+
  facet_grid(N2+N1~Estimator)

ggplot(sdat, aes(y=RMSE, x=TrueValue, group=TrueValue))+
  geom_boxplot()+
  labs(y="Root Mean Square Error")+
  facet_grid(N2+N1~Estimator)

c <- sdat %>%
  group_by(Estimator, N2, N1, TrueValue) %>%
  summarise(est = weighted.mean(estMean, wi),
            RB = weighted.mean(RB, wi),
            RMSE = weighted.mean(RMSE, wi),
            Bias = weighted.mean(Bias, wi),
            SampVar =weighted.mean(SampVar, wi))

c1 <- cbind(c[ c$Estimator == 'MLR', c( 'N2','N1','TrueValue', 'est', 'RB', 'RMSE')], 
           c[ c$Estimator == 'ULSMV', c('est', 'RB', 'RMSE')], 
           c[ c$Estimator == 'WLSMV', c('est', 'RB', 'RMSE')])
colnames(c1) <- c('N2','N1','True Value',  rep(c('est', 'RB', 'RMSE'), 3))

kable(c1, format='html', digits=3, row.names = F) %>%
  kable_styling(full_width = T) %>%
  add_header_above(c(' '=3, 'MLR'=3, 'ULSMV'=3, 'WLSMV'=3))
MLR
ULSMV
WLSMV
N2 N1 True Value est RB RMSE est RB RMSE est RB RMSE
30 5 0.111 0.203 82.898 0.049 0.293 164.389 0.178 0.174 56.858 0.032
30 5 1 1.042 4.158 0.357 1.157 15.708 0.486 0.825 -17.523 0.303
30 10 0.111 0.165 48.761 0.025 0.181 62.624 0.035 0.132 18.851 0.016
30 10 1 0.985 -1.541 0.182 1.070 7.024 0.229 0.838 -16.199 0.195
30 30 0.111 0.140 25.817 0.012 0.145 30.538 0.015 0.117 5.703 0.009
30 30 1 0.985 -1.476 0.129 1.044 4.443 0.153 0.911 -8.907 0.151
50 5 0.111 0.173 55.947 0.027 0.193 73.451 0.035 0.155 39.529 0.020
50 5 1 0.998 -0.230 0.175 1.056 5.559 0.215 0.886 -11.355 0.183
50 10 0.111 0.143 28.541 0.012 0.145 30.550 0.014 0.120 7.834 0.009
50 10 1 0.995 -0.505 0.109 1.042 4.158 0.125 0.912 -8.821 0.118
50 30 0.111 0.129 16.498 0.007 0.131 17.883 0.008 0.115 3.995 0.006
50 30 1 0.992 -0.760 0.082 1.022 2.176 0.086 0.946 -5.413 0.089
100 5 0.111 0.142 27.766 0.011 0.143 28.669 0.012 0.126 13.316 0.009
100 5 1 1.002 0.170 0.081 1.032 3.219 0.091 0.961 -3.892 0.085
100 10 0.111 0.123 11.200 0.006 0.124 11.666 0.006 0.113 1.866 0.005
100 10 1 1.000 0.031 0.051 1.017 1.744 0.055 0.963 -3.673 0.055
100 30 0.111 0.120 8.373 0.004 0.121 9.254 0.004 0.114 2.669 0.004
100 30 1 0.995 -0.529 0.041 1.010 0.985 0.041 0.972 -2.807 0.043
200 5 0.111 0.123 10.756 0.005 0.122 9.923 0.005 0.115 3.748 0.005
200 5 1 1.005 0.513 0.039 1.019 1.946 0.041 0.988 -1.236 0.039
200 10 0.111 0.116 4.706 0.003 0.117 5.295 0.003 0.111 0.303 0.003
200 10 1 0.999 -0.147 0.026 1.009 0.911 0.027 0.984 -1.591 0.026
200 30 0.111 0.112 0.783 0.002 0.113 1.675 0.002 0.108 -2.360 0.002
200 30 1 0.998 -0.185 0.020 1.007 0.658 0.020 0.988 -1.186 0.020

Relative Efficiency by Sample Sizes

c <- sdat %>%
  group_by(Estimator, N2, N1) %>%
  summarise(mu = weighted.mean(muRE, wi),
            mw = weighted.mean(mwRE, wi),
            uw = weighted.mean(uwRE, wi))

c1 <- cbind(c[ c$Estimator == 'MLR', c( 'N2','N1', 'mu', 'mw', 'uw')])
colnames(c1) <- c('N2','N1',c('MLR/ULSMV', 'MLR/WLSMV', 'ULSMV/WLSMV'))

kable(c1, format='html', digits=3, row.names = F) %>%
  kable_styling(full_width = T)
N2 N1 MLR/ULSMV MLR/WLSMV ULSMV/WLSMV
30 5 1.958 2.32 1.51
30 10 1.127 1.39 1.25
30 30 1.037 1.22 1.17
50 5 1.337 1.66 1.25
50 10 1.054 1.22 1.15
50 30 0.998 1.08 1.08
100 5 1.108 1.20 1.09
100 10 1.019 1.10 1.08
100 30 0.992 1.02 1.03
200 5 1.023 1.07 1.05
200 10 0.995 1.03 1.04
200 30 0.992 1.00 1.01

Relative Efficiency by All Conditions

c <- sdat %>%
  group_by(Estimator, N2, N1, ICC_OV, ICC_LV) %>%
  summarise(mu = weighted.mean(muRE, wi),
            mw = weighted.mean(mwRE, wi),
            uw = weighted.mean(uwRE, wi))

c1 <- cbind(c[ c$Estimator == 'MLR', c( 'N2','N1','ICC_OV', 'ICC_LV',  'mu', 'mw', 'uw')])
colnames(c1) <- c('N2','N1', 'ICC_OV', 'ICC_LV',c('MLR/ULSMV', 'MLR/WLSMV', 'ULSMV/WLSMV'))

kable(c1, format='html', digits=3, row.names = F) %>%
  kable_styling(full_width = T)
N2 N1 ICC_OV ICC_LV MLR/ULSMV MLR/WLSMV ULSMV/WLSMV
30 5 0.1 0.1 6.458 5.380 0.844
30 5 0.1 0.5 2.940 4.333 1.474
30 5 0.3 0.1 1.106 1.879 1.697
30 5 0.3 0.5 0.930 1.090 1.171
30 5 0.5 0.1 0.568 1.901 3.363
30 5 0.5 0.5 0.891 1.228 1.375
30 10 0.1 0.1 1.576 1.506 0.956
30 10 0.1 0.5 1.748 2.396 1.368
30 10 0.3 0.1 1.004 1.502 1.497
30 10 0.3 0.5 0.985 1.066 1.082
30 10 0.5 0.1 0.845 1.539 1.822
30 10 0.5 0.5 0.883 1.006 1.140
30 30 0.1 0.1 1.076 1.067 0.991
30 30 0.1 0.5 1.300 1.856 1.427
30 30 0.3 0.1 1.013 1.293 1.276
30 30 0.3 0.5 0.979 0.995 1.016
30 30 0.5 0.1 0.838 1.169 1.395
30 30 0.5 0.5 0.921 0.988 1.074
50 5 0.1 0.1 2.573 2.613 1.016
50 5 0.1 0.5 2.270 3.491 1.536
50 5 0.3 0.1 1.013 1.440 1.422
50 5 0.3 0.5 0.962 1.050 1.091
50 5 0.5 0.1 0.936 1.482 1.583
50 5 0.5 0.5 0.915 1.023 1.118
50 10 0.1 0.1 1.150 1.114 0.969
50 10 0.1 0.5 1.322 1.834 1.386
50 10 0.3 0.1 1.028 1.246 1.212
50 10 0.3 0.5 0.982 1.011 1.029
50 10 0.5 0.1 0.895 1.282 1.433
50 10 0.5 0.5 0.935 0.983 1.052
50 30 0.1 0.1 1.011 0.988 0.977
50 30 0.1 0.5 1.078 1.365 1.266
50 30 0.3 0.1 0.992 1.098 1.107
50 30 0.3 0.5 0.995 0.975 0.980
50 30 0.5 0.1 0.907 1.107 1.220
50 30 0.5 0.5 0.966 0.977 1.011
100 5 0.1 0.1 1.367 1.276 0.933
100 5 0.1 0.5 1.477 1.663 1.125
100 5 0.3 0.1 0.991 1.164 1.174
100 5 0.3 0.5 0.981 0.991 1.010
100 5 0.5 0.1 0.951 1.194 1.256
100 5 0.5 0.5 0.954 1.030 1.080
100 10 0.1 0.1 1.032 1.032 1.000
100 10 0.1 0.5 1.121 1.379 1.230
100 10 0.3 0.1 1.011 1.096 1.085
100 10 0.3 0.5 0.988 1.002 1.013
100 10 0.5 0.1 0.963 1.105 1.147
100 10 0.5 0.5 0.979 0.996 1.018
100 30 0.1 0.1 1.004 0.994 0.990
100 30 0.1 0.5 1.013 1.096 1.081
100 30 0.3 0.1 1.002 1.038 1.037
100 30 0.3 0.5 0.985 0.992 1.007
100 30 0.5 0.1 0.945 1.046 1.107
100 30 0.5 0.5 0.991 0.979 0.987
200 5 0.1 0.1 1.048 1.046 0.998
200 5 0.1 0.5 1.123 1.244 1.108
200 5 0.3 0.1 1.019 1.056 1.037
200 5 0.3 0.5 0.984 1.003 1.019
200 5 0.5 0.1 0.981 1.081 1.102
200 5 0.5 0.5 0.980 1.012 1.034
200 10 0.1 0.1 1.008 1.004 0.996
200 10 0.1 0.5 1.007 1.112 1.105
200 10 0.3 0.1 0.997 1.023 1.026
200 10 0.3 0.5 0.991 1.005 1.014
200 10 0.5 0.1 0.975 1.051 1.078
200 10 0.5 0.5 0.986 1.004 1.018
200 30 0.1 0.1 1.002 1.006 1.004
200 30 0.1 0.5 1.002 1.010 1.008
200 30 0.3 0.1 0.989 1.011 1.022
200 30 0.3 0.5 0.990 0.988 0.999
200 30 0.5 0.1 0.984 1.024 1.041
200 30 0.5 0.5 0.984 0.988 1.005

sessionInfo()
R version 3.6.3 (2020-02-29)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18362)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252 
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] xtable_1.8-4          kableExtra_1.1.0      cowplot_1.0.0        
 [4] MplusAutomation_0.7-3 data.table_1.12.8     patchwork_1.0.0      
 [7] forcats_0.5.0         stringr_1.4.0         dplyr_0.8.5          
[10] purrr_0.3.4           readr_1.3.1           tidyr_1.1.0          
[13] tibble_3.0.1          ggplot2_3.3.0         tidyverse_1.3.0      
[16] workflowr_1.6.2      

loaded via a namespace (and not attached):
 [1] httr_1.4.1        jsonlite_1.6.1    viridisLite_0.3.0 gsubfn_0.7       
 [5] modelr_0.1.8      assertthat_0.2.1  highr_0.8         pander_0.6.3     
 [9] blob_1.2.1        cellranger_1.1.0  yaml_2.2.1        pillar_1.4.4     
[13] backports_1.1.7   lattice_0.20-38   glue_1.4.1        digest_0.6.25    
[17] promises_1.1.0    rvest_0.3.5       colorspace_1.4-1  htmltools_0.4.0  
[21] httpuv_1.5.2      plyr_1.8.6        pkgconfig_2.0.3   broom_0.5.6      
[25] haven_2.3.0       scales_1.1.1      webshot_0.5.2     whisker_0.4      
[29] later_1.0.0       git2r_0.27.1      farver_2.0.3      generics_0.0.2   
[33] ellipsis_0.3.1    withr_2.2.0       cli_2.0.2         proto_1.0.0      
[37] magrittr_1.5      crayon_1.3.4      readxl_1.3.1      evaluate_0.14    
[41] fs_1.4.1          fansi_0.4.1       nlme_3.1-144      xml2_1.3.2       
[45] tools_3.6.3       hms_0.5.3         lifecycle_0.2.0   munsell_0.5.0    
[49] reprex_0.3.0      compiler_3.6.3    rlang_0.4.6       grid_3.6.3       
[53] rstudioapi_0.11   texreg_1.36.23    labeling_0.3      rmarkdown_2.1    
[57] boot_1.3-24       gtable_0.3.0      DBI_1.1.0         R6_2.4.1         
[61] lubridate_1.7.8   knitr_1.28        rprojroot_1.3-2   stringi_1.4.6    
[65] parallel_3.6.3    Rcpp_1.0.4.6      vctrs_0.3.0       dbplyr_1.4.4     
[69] tidyselect_1.1.0  xfun_0.14         coda_0.19-3