R/best.rsa2.R

In RSAtools: Advanced Response Surface Analysis

#' @title Compare two polynomial models against the data
#' @description Compare two polynomial models, for example to test parametric constraints in STEP2 of the 3-step identification strategy (see \code{\link{RSAmodel}}).
#' @param RSA_object x an object of class "RSA_object" generated by RSAmodel()
#' @param m1 First model to be compared (contained in RSA_object$models)
#' @param m2 Second model to be compared (contained in RSA_object$models)
#' @param order Fit index used to determine best-fitting model.
#' @param robust A boolean stating whether robust fit indices should be extracted (default= TRUE)
#' @return A table containing fit indices for each model 
#' 
#' @examples
#' ##### Test a variant within a family (e.g., FM26_PARALLELASYMWEAK)
#' ##Define variant as constraints
#' list_variant <- list()
#' list_variant[["variant1"]] <- c('
#' ####First-order polynomials: variant-specific			
#' 			b1 == -1/4*b2
#' ####Second-order polynomials: variant-specific			
#' 			b5 == b3/2
#' 			b4 == 0
#' ####Third-order polynomials: FM26_PARALLELASYMWEAK			
#' 			b6 == 0
#' 			b7 == 0
#' 			b9 == b8/-3		
#' ')
#' RSA_NSfit  <- RSAmodel(formula= engagement ~ needs*supplies,
#' data= sim_NSfit, model= c("FM26_PARALLELASYMWEAK","USER"),
#' user_model= list_variant)
#' ##Compare variant to best-fitting family (e.g., LRT_pvalue p > .05)
#' best.rsa2(RSA_NSfit,m1="variant1",m2="FM26_PARALLELASYMWEAK")[2,1:3]
#' @export

best.rsa2 <- function(RSA_object,m1,m2,order=c("wAIC"), robust=TRUE){

names_models <- c(m1,m2)

list_models_fitted <- list()
for(i in names_models ){
	if(RSA_object[[1]][i]!="NULL"  )
list_models_fitted[[names(RSA_object[[1]][i])]]	<- RSA_object[[1]][i]
	else(NA)
}
names(list_models_fitted)

######## Compare models using LRT
compRSA_2  <- round(.compare2(RSA_object,m1,m2,verbose=FALSE)[,-c(14:15)],3)


#######Complementary fit indices

#Robust indicators
fit_names_plain	<- c("aic_w","bic_w","aic","bic","df","chisq","pvalue","cfi","tli","rmsea","rmsea.pvalue","srmr")
fit_names_robust	<- c("aic_w","bic_w","aic","bic","df","chisq.scaled","pvalue.scaled","cfi.robust","tli.robust","rmsea.robust","rmsea.pvalue.robust","srmr")	

if(robust==F){ fit_names	<- fit_names_plain }
if(robust==T){fit_names	<- fit_names_robust	}
names_models <- names(list_models_fitted)
matrix_fitind <- matrix(nrow=length(names_models),ncol=length(fit_names),dimnames=list(names_models,fit_names))

#Extract fit indices
for(i in 1:length(list_models_fitted)){
matrix_fitind[i,-c(1:2)] <-  lavaan::fitmeasures(list_models_fitted[[i]][[1]],fit_names[-c(1:2)])
}
colnames(matrix_fitind) <- fit_names_plain

#Akaike and bayesian weights
delta_i_aic <- matrix_fitind[, "aic"] - min(matrix_fitind[, "aic"],na.rm=T)
matrix_fitind[,"aic_w"] <- exp(-1/2* delta_i_aic)/sum(exp(-1/2* delta_i_aic))
delta_i_bic <- matrix_fitind[, "bic"] - min(matrix_fitind[, "bic"],na.rm=T)
matrix_fitind[,"bic_w"] <- exp(-1/2* delta_i_bic)/sum(exp(-1/2* delta_i_bic))
matrix_fitind <- round(data.frame(matrix_fitind),3)
head(matrix_fitind)

#### Merge fit indices
merge_fit <- cbind(compRSA_2[,c("Chisq diff","Df diff","Pr(>Chisq)")], matrix_fitind[,c("aic_w","bic_w","aic","bic")], R2= compRSA_2[,c("R2")],R2adj= compRSA_2[,c("R2.adj")], matrix_fitind[,c("cfi","tli","rmsea","rmsea.pvalue","srmr")])
merge_fit <- round(data.frame(merge_fit),3)
merge_fit[,c("aic","bic")] <- round(merge_fit[,c("aic","bic")],1)
head(merge_fit)

#### Fit indices out
# Column names
head(merge_fit)
colnames(merge_fit) <- c("LRT_chi2","LRT_df","LRT_pvalue","wAIC","wBIC","AIC","BIC","R2","R2adj","CFI","TLI","RSMEA","RMSEA_pvalue","SRMR")

# Ordered by "order"
merge_fit_order <- merge_fit[-1,]
rsa_fit_order <- merge_fit_order[order(merge_fit_order[, order],decreasing=T,na.last=F),]
rsa_fit_out <- data.frame(rbind(merge_fit[1,], rsa_fit_order))
head(rsa_fit_out)

##### OUT all
compare2_out <- rsa_fit_out 
compare2_out
}