R/hier2.R
In melissagwolf/dynamic: DFI Cutoffs for Latent Variable Models
Defines functions
print.hier2
hier2
Documented in
hier2
print.hier2
#' @title Dynamic fit index (DFI) cutoffs for structural misspecification in hierarchical factor models
#'
#' @description This function generates DFI cutoffs for structural misspecification in hierarchical factor models using, by default, ML estimation.
#' The default argument is a singular argument: a \code{\link{lavaan}} object from the \code{\link{cfa}} function.
#' The function can also accommodate manual entry of the model statement and sample size.
#'
#' The app-based version of this function can be found at \href{https://dynamicfit.app/}{dynamicfit.app}.
#'
#' @param model This can either be a \code{\link{lavaan}} object from the \code{\link{cfa}} function,
#' OR a model statement written in \code{\link{lavaan}} \code{\link{model.syntax}} with standardized loadings.
#' @param n If you entered a \code{\link{lavaan}} object for model, leave this blank.
#' Otherwise, enter your sample size (numeric).
#' @param plot Displays distributions of fit indices for each level of misspecification.
#' @param manual If you entered a \code{\link{lavaan}} object, keep this set to FALSE.
#' If you manually entered standardized loadings and sample size, set this to TRUE.
#' @param reps The number of replications used in your simulation. This is set to 500 by default in both the
#' R package and the corresponding Shiny App.
#' @param estimator Which estimator to use within the simulations (enter in quotes). The default is maximum likelihood.
#'
#' @import dplyr lavaan simstandard ggplot2 stringr
#' @importFrom purrr map map_dfr map2
#' @importFrom tidyr unite extract
#' @importFrom patchwork plot_layout plot_annotation wrap_plots
#'
#' @author Daniel McNeish, Melissa G Wolf, & Patrick D Manapat
#'
#' Maintainer: Daniel McNeish <dmcneish@asu.edu>
#'
#' @rdname hier2
#'
#' @return Dynamic fit index (DFI) cutoffs for SRMR, RMSEA, and CFI.
#' @export
#'
#' @examples
#' #Manual entry example for a sample size of 2200 (manual=TRUE), from Reynolds & Keith (2017)
#'
#' manmod <- "G =~ .51*F1 + .83*F2 + .97*F3 + .83*F4 + .87*F5 + .55*Y7
#' F1 =~ .41*Y1 + .81*Y2 + .71*Y3
#' F2 =~ .79*Y4 + .64*Y5 + .81*Y6 + .22*Y7
#' F3 =~ .53*Y8 + .68*Y9 + .66*Y10
#' F4 =~ .79*Y11 + .76*Y12
#' F5 =~ .82*Y13 + .71*Y14 + .85*Y15 + .81*Y16
#' F3 ~~ .77*F4"
#' \donttest{hier2(model=manmod,n=2200,manual=TRUE)}
#'
hier2 <- function(model,n=NULL,estimator="ML",plot=FALSE,manual=FALSE,reps=500){
#If manual, expect manual (a la Shiny app)
if(manual){
tryCatch(cleanmodel(model),
error=function(err5){
if (grepl("no method for coercing this S4 class to a vector", err5)){
stop("dynamic Error: Did you accidentally include 'manual=TRUE' with a non-manually entered lavaan object?")
}
})
n <- n
model9 <- model
tryCatch(defre(model9,n),
error=function(err4){
if (grepl("non-numeric matrix extent", err4)){
stop("dynamic Error: Did you forget to include a sample size with your manually entered model?")
}
})
}else{
#Use this to rewrite error message for when someone forgot to use manual=TRUE
#But entered in model statement and sample size
#This is hacky but works, although traceback might confuse people
#https://community.rstudio.com/t/create-custom-error-messages/39058/4
tryCatch(cfa_n(model),
error=function(err){
if (grepl("trying to get slot", err)) {
stop("dynamic Error: Did you forget to use manual=TRUE?")
}
})
#Error for when someone enters an object that doesn't exist, or a non-lavaan object
tryCatch(cfa_n(model),
error=function(err2){
if (grepl("Error in base::unlist", err2)){
stop("dynamic Error: Did you enter a lavaan object? Confirm that it is a lavaan object using class(). If you do not have a lavaan object, enter the arguments manually and select manual=TRUE.")
}
})
#Use these functions to convert to manual (input is a lavaan object)
#Probably what we should expect for people using R
#need 'n' first because otherwise model will overwrite
n <- cfa_n(model)
model9 <- cfa_lavmod(model)
}
if (unstandardized(model9)>0){
stop("dynamic Error: One of your loadings or correlations has an absolute value of 1 or above (an impossible value). Please use standardized loadings. If all of your loadings are under 1, try looking for a missing decimal somewhere in your model statement.")
}
if (number_factor_first(model9)==0){
stop("dynamic Error: The model you entered is not hierarchical - please use a different function. For example, cfaOne for one-factor models or cfaHB for multi-factor models.")
}
if (defre(model9,n)==0){
stop("dynamic Error: It is impossible to add misspecifications to a just identified model.")
}
if (number_factor_first(model9)<4){
stop("dynamic Error: The first-order factor covariance matrix is saturated and hypothetical misspecifications cannot be tested.")
}
#Create list to store outputs (table and plot)
res <- list()
#Output fit indices if someone used manual=F
#Will ignore in print statement if manual=T
#Exclamation point is how we indicate if manual = T (because default is F)
if(!manual){
if (model@Options$test=="satorra.bentler" |model@Options$test=="yuan.bentler.mplus" | model@Options$test=="yuan.bentler.mplus"){
fitted <- round(lavaan::fitmeasures(model,c("chisq.scaled","df","pvalue.scaled","srmr","rmsea.robust","cfi.robust")),3)
} else if (model@Options$test=="scaled.shifted" | model@Options$test=="mean.var.adusted"){
fitted <- round(lavaan::fitmeasures(model,c("chisq","df","pvalue","srmr","rmsea.scaled","cfi.scaled")),3)
} else if(model@Options$test=="standard" ){
fitted <- round(lavaan::fitmeasures(model,c("chisq","df","pvalue","srmr","rmsea","cfi")),3)}
fitted_m <- as.matrix(fitted)
fitted_t <- t(fitted_m)
fitted_t <- as.data.frame(fitted_t)
colnames(fitted_t) <- c("Chi-Square"," df","p-value"," SRMR"," RMSEA"," CFI")
rownames(fitted_t) <- c("")
res$fit <- fitted_t
}
#Run simulation
results <- one_df_hier2(model9,n,estimator,reps)
#Save the data and make it exportable
res$data <- fit_data(results)
#For each list element (misspecification) compute the cutoffs
misspec_sum <- purrr::map(results,~dplyr::reframe(.,SRMR_M=stats::quantile(SRMR_M, c(.05,.1)),
RMSEA_M=stats::quantile(RMSEA_M, c(.05,.1)),
CFI_M=stats::quantile(CFI_M, c(.95,.9))))
#For the true model, compute the cutoffs (these will all be the same - just need in list form)
true_sum <- purrr::map(results,~dplyr::reframe(.,SRMR_T=stats::quantile(SRMR_T, c(.95,.9)),
RMSEA_T=stats::quantile(RMSEA_T, c(.95,.9)),
CFI_T=stats::quantile(CFI_T, c(.05,.1))))
#Bind each of the misspecified cutoffs to the true cutoffs, listwise
Table <- purrr::map(misspec_sum,~base::cbind(.,true_sum[[1]]) %>%
dplyr::mutate(SRMR_R=base::round(SRMR_M,3),
RMSEA_R=base::round(RMSEA_M,3),
CFI_R=base::round(CFI_M,3),
SRMR=base::ifelse(SRMR_T<SRMR_M,SRMR_R,"NONE"),
RMSEA=base::ifelse(RMSEA_T<RMSEA_M,RMSEA_R,"NONE"),
CFI=base::ifelse(CFI_T>CFI_M,CFI_R,"NONE")) %>%
dplyr::select(SRMR,RMSEA,CFI))
#This is to clean up the table for presentation
#list is a function within mutate to apply function lead across each element
Row2 <- purrr::map_dfr(Table,~dplyr::mutate(.,SRMR_1=SRMR,
RMSEA_1=RMSEA,
CFI_1=CFI) %>%
dplyr::mutate_at(c("SRMR_1","RMSEA_1","CFI_1"),base::list(dplyr::lead)) %>%
dplyr::slice(1) %>%
dplyr::mutate(SRMR=base::ifelse(base::is.character(SRMR),SRMR_1,"--"),
RMSEA=base::ifelse(base::is.character(RMSEA),RMSEA_1,"--"),
CFI=base::ifelse(base::is.character(CFI),CFI_1,"--"),
SRMR=stringr::str_replace_all(base::as.character(SRMR),"0\\.","."),
RMSEA=stringr::str_replace_all(base::as.character(RMSEA),"0\\.","."),
CFI=stringr::str_replace_all(base::as.character(CFI),"0\\.",".")) %>%
dplyr::select(SRMR,RMSEA,CFI))
#Still cleaning
#Unlist Table
Table_C <- purrr::map_dfr(Table,~dplyr::mutate(.,SRMR=stringr::str_replace_all(base::as.character(SRMR),"0\\.","."),
RMSEA=stringr::str_replace_all(base::as.character(RMSEA),"0\\.","."),
CFI=stringr::str_replace_all(base::as.character(CFI),"0\\.",".")))
#Cleaning
Table_C[base::seq(2,nrow(Table_C),by=2),] <- Row2
#Create row names for level
Table_C$levelnum <- base::paste("Level", base::rep(1:(base::nrow(Table_C)/2),each=2))
#Create row names for proportions
Table_C$cut <- base::rep(c("95/5","90/10"))
#Add rownames to final table
Final_Table <- Table_C %>%
tidyr::unite(Cut,levelnum,cut,sep=": ") %>%
tibble::column_to_rownames(var='Cut')
#Put into list
res$cutoffs <- Final_Table
#If user selects plot = T
if(plot){
#For each list element (misspecification) compute the cutoffs
misspec_sum <- purrr::map(results,~dplyr::reframe(.,SRMR_M=stats::quantile(SRMR_M, c(.05,.1)),
RMSEA_M=stats::quantile(RMSEA_M, c(.05,.1)),
CFI_M=stats::quantile(CFI_M, c(.95,.9))))
#For the true model, compute the cutoffs (these will all be the same - just need in list form)
true_sum <- purrr::map(results,~dplyr::reframe(.,SRMR_T=stats::quantile(SRMR_T, c(.95,.9)),
RMSEA_T=stats::quantile(RMSEA_T, c(.95,.9)),
CFI_T=stats::quantile(CFI_T, c(.05,.1))))
#Select just those variables and rename columns to be the same
Misspec_dat <- purrr::map(results,~dplyr::select(.,SRMR_M:Type_M) %>%
`colnames<-`(c("SRMR","RMSEA","CFI","Model")))
#Select just those variables and rename columns to be the same
True_dat <- purrr::map(results,~dplyr::select(.,SRMR_T:Type_T) %>%
`colnames<-`(c("SRMR","RMSEA","CFI","Model")))
#For each element in the list, bind the misspecified cutoffs to the true cutoffs
#rbind doesn't work well with lists (needs do.call statement)
plot <- base::lapply(base::seq(base::length(Misspec_dat)),function(x) dplyr::bind_rows(Misspec_dat[x],True_dat[x]))
#Plot SRMR. Need map2 and data=.x (can't remember why).
SRMR_plot <- purrr::map2(plot,misspec_sum,~ggplot(data=.x,aes(x=SRMR,fill=Model))+
geom_histogram(position="identity",
alpha=.5, bins=30)+
scale_fill_manual(values=c("#E9798C","#66C2F5"))+
geom_vline(aes(xintercept=.y$SRMR_M[1],
linetype="misspec_sum$SRMR_M[1]",color="misspec_sum$SRMR_M[1]"),
linewidth=.6)+
geom_vline(aes(xintercept=.08,
linetype=".08",color=".08"),
linewidth=.75)+
scale_color_manual(name="Cutoff Values",
labels=c("Hu & Bentler Cutoff","Dynamic Cutoff"),
values=c("misspec_sum$SRMR_M[1]"="black",
".08"="black"))+
scale_linetype_manual(name="Cutoff Values",
labels=c("Hu & Bentler Cutoff","Dynamic Cutoff"),
values=c("misspec_sum$SRMR_M[1]"="longdash",
".08"="dotted"))+
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.background = element_blank(),
axis.line = element_line(color="black"),
legend.position = "none",
legend.title = element_blank(),
legend.box = "vertical"))
#Plot RMSEA. Need map2 and data=.x (can't remember why).
RMSEA_plot <- purrr::map2(plot,misspec_sum,~ggplot(data=.x,aes(x=RMSEA,fill=Model))+
geom_histogram(position="identity",
alpha=.5, bins=30)+
scale_fill_manual(values=c("#E9798C","#66C2F5"))+
geom_vline(aes(xintercept=.y$RMSEA_M[1],
linetype="misspec_sum$RMSEA_M[1]",color="misspec_sum$RMSEA_M[1]"),
linewidth=.6)+
geom_vline(aes(xintercept=.06,
linetype=".06",color=".06"),
linewidth=.75)+
scale_color_manual(name="Cutoff Values",
labels=c("Hu & Bentler Cutoff","Dynamic Cutoff"),
values=c("misspec_sum$RMSEA_M[1]"="black",
".06"="black"))+
scale_linetype_manual(name="Cutoff Values",
labels=c("Hu & Bentler Cutoff","Dynamic Cutoff"),
values=c("misspec_sum$RMSEA_M[1]"="longdash",
".06"="dotted"))+
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.background = element_blank(),
axis.line = element_line(color="black"),
legend.position = "none",
legend.title = element_blank(),
legend.box = "vertical"))
#Plot CFI. Need map2 and data=.x (can't remember why).
CFI_plot <- purrr::map2(plot,misspec_sum,~ggplot(data=.x,aes(x=CFI,fill=Model))+
geom_histogram(position="identity",
alpha=.5, bins=30)+
scale_fill_manual(values=c("#E9798C","#66C2F5"))+
geom_vline(aes(xintercept=.y$CFI_M[1],
linetype="misspec_sum$CFI_M[1]",color="misspec_sum$CFI_M[1]"),
linewidth=.6)+
geom_vline(aes(xintercept=.95,
linetype=".95",color=".95"),
linewidth=.75)+
scale_color_manual(name="Cutoff Values",
labels=c("Hu & Bentler Cutoff","Dynamic Cutoff"),
values=c("misspec_sum$CFI_M[1]"="black",
".95"="black"))+
scale_linetype_manual(name="Cutoff Values",
labels=c("Hu & Bentler Cutoff","Dynamic Cutoff"),
values=c("misspec_sum$CFI_M[1]"="longdash",
".95"="dotted"))+
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.background = element_blank(),
axis.line = element_line(color="black"),
legend.position = "none",
legend.title = element_blank(),
legend.box = "vertical"))
#Create a list with the plots combined for each severity level
plots_combo <- base::lapply(base::seq(base::length(plot)),function(x) c(SRMR_plot[x],RMSEA_plot[x],CFI_plot[x]))
#Add a collective legend and title with the level indicator
plots <- base::lapply(base::seq(base::length(plots_combo)), function(x) patchwork::wrap_plots(plots_combo[[x]])+
plot_layout(guides = "collect")+
plot_annotation(title=paste("Level", x))
& theme(legend.position = 'bottom'))
#Put into list
res$plots <- plots
}
#Create object (necessary for subsequent print statement)
class(res) <- 'hier2'
return(res)
}
#' @method print hier2
#' @param x hier2 object
#' @param ... other print parameters
#' @rdname hier2
#' @export
#Print suppression/organization statement for list
#Needs same name as class, not function name
#Need to add ... param or will get error message in CMD check
print.hier2 <- function(x,...){
#Automatically return fit cutoffs
base::cat("Your DFI cutoffs: \n")
base::print(x$cutoffs)
#Only print fit indices from lavaan object if someone submits a lavaan object
if(!is.null(x$fit)){
base::cat("\n")
base::cat("Empirical fit indices: \n")
base::print(x$fit)
}
if(!is.null(x$plots)){
base::cat("\n The distributions for each level are in the Plots tab \n")
base::print(x$plots)
}
#Hides this function
base::invisible()
}
melissagwolf/dynamic documentation built on June 29, 2024, 6:24 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions print.hier2 hier2
Documented in hier2 print.hier2
#' @title Dynamic fit index (DFI) cutoffs for structural misspecification in hierarchical factor models
#'
#' @description This function generates DFI cutoffs for structural misspecification in hierarchical factor models using, by default, ML estimation.
#' The default argument is a singular argument: a \code{\link{lavaan}} object from the \code{\link{cfa}} function.
#' The function can also accommodate manual entry of the model statement and sample size.
#'
#' The app-based version of this function can be found at \href{https://dynamicfit.app/}{dynamicfit.app}.
#'
#' @param model This can either be a \code{\link{lavaan}} object from the \code{\link{cfa}} function,
#' OR a model statement written in \code{\link{lavaan}} \code{\link{model.syntax}} with standardized loadings.
#' @param n If you entered a \code{\link{lavaan}} object for model, leave this blank.
#' Otherwise, enter your sample size (numeric).
#' @param plot Displays distributions of fit indices for each level of misspecification.
#' @param manual If you entered a \code{\link{lavaan}} object, keep this set to FALSE.
#' If you manually entered standardized loadings and sample size, set this to TRUE.
#' @param reps The number of replications used in your simulation. This is set to 500 by default in both the
#' R package and the corresponding Shiny App.
#' @param estimator Which estimator to use within the simulations (enter in quotes). The default is maximum likelihood.
#'
#' @import dplyr lavaan simstandard ggplot2 stringr
#' @importFrom purrr map map_dfr map2
#' @importFrom tidyr unite extract
#' @importFrom patchwork plot_layout plot_annotation wrap_plots
#'
#' @author Daniel McNeish, Melissa G Wolf, & Patrick D Manapat
#'
#' Maintainer: Daniel McNeish <dmcneish@asu.edu>
#'
#' @rdname hier2
#'
#' @return Dynamic fit index (DFI) cutoffs for SRMR, RMSEA, and CFI.
#' @export
#'
#' @examples
#' #Manual entry example for a sample size of 2200 (manual=TRUE), from Reynolds & Keith (2017)
#'
#' manmod <- "G =~ .51*F1 + .83*F2 + .97*F3 + .83*F4 + .87*F5 + .55*Y7
#' F1 =~ .41*Y1 + .81*Y2 + .71*Y3
#' F2 =~ .79*Y4 + .64*Y5 + .81*Y6 + .22*Y7
#' F3 =~ .53*Y8 + .68*Y9 + .66*Y10
#' F4 =~ .79*Y11 + .76*Y12
#' F5 =~ .82*Y13 + .71*Y14 + .85*Y15 + .81*Y16
#' F3 ~~ .77*F4"
#' \donttest{hier2(model=manmod,n=2200,manual=TRUE)}
#'
hier2 <- function(model,n=NULL,estimator="ML",plot=FALSE,manual=FALSE,reps=500){
#If manual, expect manual (a la Shiny app)
if(manual){
tryCatch(cleanmodel(model),
error=function(err5){
if (grepl("no method for coercing this S4 class to a vector", err5)){
stop("dynamic Error: Did you accidentally include 'manual=TRUE' with a non-manually entered lavaan object?")
}
})
n <- n
model9 <- model
tryCatch(defre(model9,n),
error=function(err4){
if (grepl("non-numeric matrix extent", err4)){
stop("dynamic Error: Did you forget to include a sample size with your manually entered model?")
}
})
}else{
#Use this to rewrite error message for when someone forgot to use manual=TRUE
#But entered in model statement and sample size
#This is hacky but works, although traceback might confuse people
#https://community.rstudio.com/t/create-custom-error-messages/39058/4
tryCatch(cfa_n(model),
error=function(err){
if (grepl("trying to get slot", err)) {
stop("dynamic Error: Did you forget to use manual=TRUE?")
}
})
#Error for when someone enters an object that doesn't exist, or a non-lavaan object
tryCatch(cfa_n(model),
error=function(err2){
if (grepl("Error in base::unlist", err2)){
stop("dynamic Error: Did you enter a lavaan object? Confirm that it is a lavaan object using class(). If you do not have a lavaan object, enter the arguments manually and select manual=TRUE.")
}
})
#Use these functions to convert to manual (input is a lavaan object)
#Probably what we should expect for people using R
#need 'n' first because otherwise model will overwrite
n <- cfa_n(model)
model9 <- cfa_lavmod(model)
}
if (unstandardized(model9)>0){
stop("dynamic Error: One of your loadings or correlations has an absolute value of 1 or above (an impossible value). Please use standardized loadings. If all of your loadings are under 1, try looking for a missing decimal somewhere in your model statement.")
}
if (number_factor_first(model9)==0){
stop("dynamic Error: The model you entered is not hierarchical - please use a different function. For example, cfaOne for one-factor models or cfaHB for multi-factor models.")
}
if (defre(model9,n)==0){
stop("dynamic Error: It is impossible to add misspecifications to a just identified model.")
}
if (number_factor_first(model9)<4){
stop("dynamic Error: The first-order factor covariance matrix is saturated and hypothetical misspecifications cannot be tested.")
}
#Create list to store outputs (table and plot)
res <- list()
#Output fit indices if someone used manual=F
#Will ignore in print statement if manual=T
#Exclamation point is how we indicate if manual = T (because default is F)
if(!manual){
if (model@Options$test=="satorra.bentler" |model@Options$test=="yuan.bentler.mplus" | model@Options$test=="yuan.bentler.mplus"){
fitted <- round(lavaan::fitmeasures(model,c("chisq.scaled","df","pvalue.scaled","srmr","rmsea.robust","cfi.robust")),3)
} else if (model@Options$test=="scaled.shifted" | model@Options$test=="mean.var.adusted"){
fitted <- round(lavaan::fitmeasures(model,c("chisq","df","pvalue","srmr","rmsea.scaled","cfi.scaled")),3)
} else if(model@Options$test=="standard" ){
fitted <- round(lavaan::fitmeasures(model,c("chisq","df","pvalue","srmr","rmsea","cfi")),3)}
fitted_m <- as.matrix(fitted)
fitted_t <- t(fitted_m)
fitted_t <- as.data.frame(fitted_t)
colnames(fitted_t) <- c("Chi-Square"," df","p-value"," SRMR"," RMSEA"," CFI")
rownames(fitted_t) <- c("")
res$fit <- fitted_t
}
#Run simulation
results <- one_df_hier2(model9,n,estimator,reps)
#Save the data and make it exportable
res$data <- fit_data(results)
#For each list element (misspecification) compute the cutoffs
misspec_sum <- purrr::map(results,~dplyr::reframe(.,SRMR_M=stats::quantile(SRMR_M, c(.05,.1)),
RMSEA_M=stats::quantile(RMSEA_M, c(.05,.1)),
CFI_M=stats::quantile(CFI_M, c(.95,.9))))
#For the true model, compute the cutoffs (these will all be the same - just need in list form)
true_sum <- purrr::map(results,~dplyr::reframe(.,SRMR_T=stats::quantile(SRMR_T, c(.95,.9)),
RMSEA_T=stats::quantile(RMSEA_T, c(.95,.9)),
CFI_T=stats::quantile(CFI_T, c(.05,.1))))
#Bind each of the misspecified cutoffs to the true cutoffs, listwise
Table <- purrr::map(misspec_sum,~base::cbind(.,true_sum[[1]]) %>%
dplyr::mutate(SRMR_R=base::round(SRMR_M,3),
RMSEA_R=base::round(RMSEA_M,3),
CFI_R=base::round(CFI_M,3),
SRMR=base::ifelse(SRMR_T<SRMR_M,SRMR_R,"NONE"),
RMSEA=base::ifelse(RMSEA_T<RMSEA_M,RMSEA_R,"NONE"),
CFI=base::ifelse(CFI_T>CFI_M,CFI_R,"NONE")) %>%
dplyr::select(SRMR,RMSEA,CFI))
#This is to clean up the table for presentation
#list is a function within mutate to apply function lead across each element
Row2 <- purrr::map_dfr(Table,~dplyr::mutate(.,SRMR_1=SRMR,
RMSEA_1=RMSEA,
CFI_1=CFI) %>%
dplyr::mutate_at(c("SRMR_1","RMSEA_1","CFI_1"),base::list(dplyr::lead)) %>%
dplyr::slice(1) %>%
dplyr::mutate(SRMR=base::ifelse(base::is.character(SRMR),SRMR_1,"--"),
RMSEA=base::ifelse(base::is.character(RMSEA),RMSEA_1,"--"),
CFI=base::ifelse(base::is.character(CFI),CFI_1,"--"),
SRMR=stringr::str_replace_all(base::as.character(SRMR),"0\\.","."),
RMSEA=stringr::str_replace_all(base::as.character(RMSEA),"0\\.","."),
CFI=stringr::str_replace_all(base::as.character(CFI),"0\\.",".")) %>%
dplyr::select(SRMR,RMSEA,CFI))
#Still cleaning
#Unlist Table
Table_C <- purrr::map_dfr(Table,~dplyr::mutate(.,SRMR=stringr::str_replace_all(base::as.character(SRMR),"0\\.","."),
RMSEA=stringr::str_replace_all(base::as.character(RMSEA),"0\\.","."),
CFI=stringr::str_replace_all(base::as.character(CFI),"0\\.",".")))
#Cleaning
Table_C[base::seq(2,nrow(Table_C),by=2),] <- Row2
#Create row names for level
Table_C$levelnum <- base::paste("Level", base::rep(1:(base::nrow(Table_C)/2),each=2))
#Create row names for proportions
Table_C$cut <- base::rep(c("95/5","90/10"))
#Add rownames to final table
Final_Table <- Table_C %>%
tidyr::unite(Cut,levelnum,cut,sep=": ") %>%
tibble::column_to_rownames(var='Cut')
#Put into list
res$cutoffs <- Final_Table
#If user selects plot = T
if(plot){
#For each list element (misspecification) compute the cutoffs
misspec_sum <- purrr::map(results,~dplyr::reframe(.,SRMR_M=stats::quantile(SRMR_M, c(.05,.1)),
RMSEA_M=stats::quantile(RMSEA_M, c(.05,.1)),
CFI_M=stats::quantile(CFI_M, c(.95,.9))))
#For the true model, compute the cutoffs (these will all be the same - just need in list form)
true_sum <- purrr::map(results,~dplyr::reframe(.,SRMR_T=stats::quantile(SRMR_T, c(.95,.9)),
RMSEA_T=stats::quantile(RMSEA_T, c(.95,.9)),
CFI_T=stats::quantile(CFI_T, c(.05,.1))))
#Select just those variables and rename columns to be the same
Misspec_dat <- purrr::map(results,~dplyr::select(.,SRMR_M:Type_M) %>%
`colnames<-`(c("SRMR","RMSEA","CFI","Model")))
#Select just those variables and rename columns to be the same
True_dat <- purrr::map(results,~dplyr::select(.,SRMR_T:Type_T) %>%
`colnames<-`(c("SRMR","RMSEA","CFI","Model")))
#For each element in the list, bind the misspecified cutoffs to the true cutoffs
#rbind doesn't work well with lists (needs do.call statement)
plot <- base::lapply(base::seq(base::length(Misspec_dat)),function(x) dplyr::bind_rows(Misspec_dat[x],True_dat[x]))
#Plot SRMR. Need map2 and data=.x (can't remember why).
SRMR_plot <- purrr::map2(plot,misspec_sum,~ggplot(data=.x,aes(x=SRMR,fill=Model))+
geom_histogram(position="identity",
alpha=.5, bins=30)+
scale_fill_manual(values=c("#E9798C","#66C2F5"))+
geom_vline(aes(xintercept=.y$SRMR_M[1],
linetype="misspec_sum$SRMR_M[1]",color="misspec_sum$SRMR_M[1]"),
linewidth=.6)+
geom_vline(aes(xintercept=.08,
linetype=".08",color=".08"),
linewidth=.75)+
scale_color_manual(name="Cutoff Values",
labels=c("Hu & Bentler Cutoff","Dynamic Cutoff"),
values=c("misspec_sum$SRMR_M[1]"="black",
".08"="black"))+
scale_linetype_manual(name="Cutoff Values",
labels=c("Hu & Bentler Cutoff","Dynamic Cutoff"),
values=c("misspec_sum$SRMR_M[1]"="longdash",
".08"="dotted"))+
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.background = element_blank(),
axis.line = element_line(color="black"),
legend.position = "none",
legend.title = element_blank(),
legend.box = "vertical"))
#Plot RMSEA. Need map2 and data=.x (can't remember why).
RMSEA_plot <- purrr::map2(plot,misspec_sum,~ggplot(data=.x,aes(x=RMSEA,fill=Model))+
geom_histogram(position="identity",
alpha=.5, bins=30)+
scale_fill_manual(values=c("#E9798C","#66C2F5"))+
geom_vline(aes(xintercept=.y$RMSEA_M[1],
linetype="misspec_sum$RMSEA_M[1]",color="misspec_sum$RMSEA_M[1]"),
linewidth=.6)+
geom_vline(aes(xintercept=.06,
linetype=".06",color=".06"),
linewidth=.75)+
scale_color_manual(name="Cutoff Values",
labels=c("Hu & Bentler Cutoff","Dynamic Cutoff"),
values=c("misspec_sum$RMSEA_M[1]"="black",
".06"="black"))+
scale_linetype_manual(name="Cutoff Values",
labels=c("Hu & Bentler Cutoff","Dynamic Cutoff"),
values=c("misspec_sum$RMSEA_M[1]"="longdash",
".06"="dotted"))+
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.background = element_blank(),
axis.line = element_line(color="black"),
legend.position = "none",
legend.title = element_blank(),
legend.box = "vertical"))
#Plot CFI. Need map2 and data=.x (can't remember why).
CFI_plot <- purrr::map2(plot,misspec_sum,~ggplot(data=.x,aes(x=CFI,fill=Model))+
geom_histogram(position="identity",
alpha=.5, bins=30)+
scale_fill_manual(values=c("#E9798C","#66C2F5"))+
geom_vline(aes(xintercept=.y$CFI_M[1],
linetype="misspec_sum$CFI_M[1]",color="misspec_sum$CFI_M[1]"),
linewidth=.6)+
geom_vline(aes(xintercept=.95,
linetype=".95",color=".95"),
linewidth=.75)+
scale_color_manual(name="Cutoff Values",
labels=c("Hu & Bentler Cutoff","Dynamic Cutoff"),
values=c("misspec_sum$CFI_M[1]"="black",
".95"="black"))+
scale_linetype_manual(name="Cutoff Values",
labels=c("Hu & Bentler Cutoff","Dynamic Cutoff"),
values=c("misspec_sum$CFI_M[1]"="longdash",
".95"="dotted"))+
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.background = element_blank(),
axis.line = element_line(color="black"),
legend.position = "none",
legend.title = element_blank(),
legend.box = "vertical"))
#Create a list with the plots combined for each severity level
plots_combo <- base::lapply(base::seq(base::length(plot)),function(x) c(SRMR_plot[x],RMSEA_plot[x],CFI_plot[x]))
#Add a collective legend and title with the level indicator
plots <- base::lapply(base::seq(base::length(plots_combo)), function(x) patchwork::wrap_plots(plots_combo[[x]])+
plot_layout(guides = "collect")+
plot_annotation(title=paste("Level", x))
& theme(legend.position = 'bottom'))
#Put into list
res$plots <- plots
}
#Create object (necessary for subsequent print statement)
class(res) <- 'hier2'
return(res)
}
#' @method print hier2
#' @param x hier2 object
#' @param ... other print parameters
#' @rdname hier2
#' @export
#Print suppression/organization statement for list
#Needs same name as class, not function name
#Need to add ... param or will get error message in CMD check
print.hier2 <- function(x,...){
#Automatically return fit cutoffs
base::cat("Your DFI cutoffs: \n")
base::print(x$cutoffs)
#Only print fit indices from lavaan object if someone submits a lavaan object
if(!is.null(x$fit)){
base::cat("\n")
base::cat("Empirical fit indices: \n")
base::print(x$fit)
}
if(!is.null(x$plots)){
base::cat("\n The distributions for each level are in the Plots tab \n")
base::print(x$plots)
}
#Hides this function
base::invisible()
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.