R/best.rsa.R

In RSAtools: Advanced Response Surface Analysis

#' @title Compare a list of polynomial models against the data
#' @description Compares any number of predefined or user-specific polynomial models and extracts their fit indices, thereby establishing best-fitting solutions.
#' @details This function compares models based on information-theoretic criteria and statistical tests. The cubic saturated polynomial provides a benchmark reference for fit, against which predefined polynomial families (37 to date) or user-specific variants of these families are compared for absolute fit (likelihood ratio test), parsimony (wAIC), explained variance (adjusted R2), and ordinary SEM criteria (e.g., CFI, TLI, RMSEA, SRMR).
#' @param RSA_object x an object of class "RSA_object" generated by RSAmodel(). If x is already generated by best.rsa(), then best.rsa will simply apply reordering according to "order".
#' @param order Single or vector of fit indices used to determine best-fitting polynomial families. The output matrix is ordered based on this fit index
#' @param robust Should robust fit indices should be extracted? (default= TRUE)
#' @return A table containing fit indices for each model 
#' 
#' @examples
#' #####ESTIMATE RSA OBJECT
#' RSA_step1 <-  RSAmodel(engagement ~ needs*supplies,
#' data= sim_NSfit, model= c("CUBIC","FM8_INCONG","FM9_INCONG","FM20_ASYMCONG",
#' "FM21_ASYMCONG","FM26_PARALLELASYMWEAK"))
#' ##### COMPARE POLYNOMIAL FAMILIES FROM THE RSA OBJECT
#' RSA_step1_fit <- best.rsa(RSA_step1,order=c("wAIC"))
#' names(RSA_step1$models)
#' #Inspect best-fitting family model
#' summary(RSA_step1$models$FM26_PARALLELASYMWEAK)
#' @export
#' @importFrom lavaan lavTestLRT fitmeasures inspect nobs

best.rsa <- function(RSA_object, order=c("wAIC","R2adj"),robust=F){

if(inherits(RSA_object, "data.frame")){
best_all <- RSA_object		
best <- best_all[-1,]

if(length(order)==1){
rsa_fit_out <- best[order(best[,order],decreasing=T,na.last=F),]
}

if(length(order)==2){
rsa_fit_out <- best[order(best[,order[[1]] ],best[,order[[2]] ],decreasing=T,na.last=F),]
}

if(length(order)==3){
rsa_fit_out <- best[order(best[,order[[1]] ],best[,order[[2]] ], best[,order[[3]] ],decreasing=T,na.last=F),]
}
rsa_fit_out <- data.frame(rbind(best_all[1,], rsa_fit_out))		
}



if(!inherits(RSA_object, "data.frame")) {
####Lists to compare
list_models_all <- list()
list_models_fitted <- list()

for(i in 1:length(RSA_object[[1]]) ){
	if(RSA_object[[1]][i]!="NULL" )
list_models_all[[names(RSA_object[[1]][i])]]	<- RSA_object[[1]][i]
	else(NA)
}

for(j in 1:length(list_models_all) ){
		if( is.na(lavaan::lavTestLRT(list_models_all[[j]][[1]])[2,"Chisq"])  )
		NA
		
		else if( !is.na(lavaan::lavTestLRT(list_models_all[[j]][[1]])[2,"Chisq"])  )
list_models_fitted[[names(list_models_all[[j]])]]	<- list_models_all[[j]]
	else(NA)
}


#######Complementary fit indices

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

# Validate estimator compatibility with robust indicators
if (robust == TRUE) {
  estimators <- sapply(list_models_fitted, function(x) x[[1]]@Options$estimator)
  invalid_estimators <- estimators[!estimators %in% c("MLR", "WLSMV", "ULSMV", "WLSM")]

  if (length(invalid_estimators) > 0) {
    stop("Error: Robust fit indices requested, but at least one model does not use a robust estimator (e.g., 'MLR', 'WLSMV'). Please set 'robust=FALSE' or use robust estimators.")
  }

  fit_names <- fit_names_robust
} else {
  fit_names <- fit_names_plain
}

# Prepare storage
names_models <- names(list_models_fitted)
matrix_fitind <- matrix(nrow = length(names_models), ncol = length(fit_names), dimnames = list(names_models, fit_names))
list_sic <- list()

# Extract fit indices
for (m in seq_along(list_models_fitted)) {
  fitRSA <- list_models_fitted[[m]][[1]]
  matrix_fitind[m, -c(1:4)] <- lavaan::fitmeasures(fitRSA, fit_names[-c(1:4)])

  # R2 and adjusted R2
  free.max <- .getFreeParameters(fitRSA)
  R2 <- lavaan::inspect(fitRSA, "r2")
  matrix_fitind[m, "R2"] <- R2
  N <- lavaan::nobs(fitRSA)
  k <- free.max - matrix_fitind[m, "df"]
  matrix_fitind[m, "R2.adj"] <- 1 - ((1 - R2)) * ((N - 1) / (N - k - 1))
}

# Standardized column names (optional fallback)
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(matrix_fitind[,c("chisq","df","pvalue")], matrix_fitind[,c("aic_w","bic_w","aic","bic")], R2= matrix_fitind[,c("R2")],R2adj= matrix_fitind[,c("R2.adj")], matrix_fitind[,c("cfi","tli","rmsea","rmsea.pvalue","srmr")])
merge_fit <- round(data.frame(merge_fit),3)
merge_fit[,c("chisq","aic","bic")] <- round(merge_fit[,c("chisq","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","RMSEA","RMSEA_pvalue","SRMR")

# Ordered by "order"
merge_fit_order <- merge_fit[-1,]
order_all <- merge_fit_order[, order[1]]
if(length(order)>1){
order_all <- list()	
for(o in 1:length(order)){
order_all[[ order[o] ]]  <- merge_fit_order[, order[o]]
}
}

if(length(order)==1){
rsa_fit_order <- merge_fit_order[order(order_all,decreasing=T,na.last=F),]
}

if(length(order)==2){
rsa_fit_order <- merge_fit_order[order(order_all[[1]],order_all[[2]],decreasing=T,na.last=F),]
}

if(length(order)==3){
rsa_fit_order <- merge_fit_order[order(order_all[[1]],order_all[[2]],order_all[[3]],decreasing=T,na.last=F),]
}
rsa_fit_out <- data.frame(rbind(merge_fit[1,], rsa_fit_order))
head(rsa_fit_out)

}

##### OUT all
compare_out <- rsa_fit_out[,c("LRT_chi2","LRT_df","LRT_pvalue","AIC","wAIC","R2adj","CFI","TLI","RMSEA","SRMR","BIC","wBIC","R2")]  
compare_out

}