Model Fit Index

Load file

load(file = "C:/Users/DELL/Documents/integration-simulations/sampleSize.RData")
library(pracma)
library(ggplot2)

BIC

Compute BIC

df$BIC = 2 * df$NLL + 2 * log(nTrials)
df$BIC_alt = 2 * df$NLL_alt + 3 * log(nTrials)

AIC

Compute AIC

df$AIC = 2 * df$NLL + 2 * 2
df$AIC_alt = 2 * df$NLL_alt + 3 * 2

Now determine

plottingInfoAll = data.frame(Outcome = rep(c("Incorrect & N.S.", 
                                             "Correct & N.S.", 
                                             "Incorrect & p < 0.05",
                                             "Correct & p < 0.05"), times = 2),
                             Count = 0,
                             MPI = rep(c("AIC", "BIC"), each = 4))

index = "AIC"
sampleSize = 60
nTrials = 60
test = vector('logical', 100000)
better = vector('logical', 100000)
for (i in 1:length(test)){
  subsample = df[sample(1:nrow(df), sampleSize, replace = T), ]
  better[i] = sum(subsample$AIC) < sum(subsample$AIC_alt)
  test[i] = t.test(subsample$AIC, subsample$AIC_alt, paired = T)$p.value < 0.05
}
plottingInfoAll$Count[which(plottingInfoAll$MPI == index & 
                              plottingInfoAll$Outcome == "Incorrect & N.S.")] = 
    sum(better==F & test==F)
plottingInfoAll$Count[which(plottingInfoAll$MPI == index & 
                            plottingInfoAll$Outcome == "Correct & N.S.")] = 
  sum(better==T & test==F)
plottingInfoAll$Count[which(plottingInfoAll$MPI == index & 
                            plottingInfoAll$Outcome == "Incorrect & p < 0.05")] = 
  sum(better==F & test==T)
plottingInfoAll$Count[which(plottingInfoAll$MPI == index & 
                            plottingInfoAll$Outcome == "Correct & p < 0.05")] = 
  sum(better==T & test==T)

index = "BIC"
test = vector('logical', 100000)
better = vector('logical', 100000)
for (i in 1:length(test)){
  subsample = df[sample(1:nrow(df), sampleSize, replace = T), ]
  better[i] = sum(subsample$BIC) < sum(subsample$BIC_alt)
  test[i] = t.test(subsample$BIC, subsample$BIC_alt, paired = T)$p.value < 0.05
}
plottingInfoAll$Count[which(plottingInfoAll$MPI == index & 
                              plottingInfoAll$Outcome == "Incorrect & N.S.")] = 
    sum(better==F & test==F)
plottingInfoAll$Count[which(plottingInfoAll$MPI == index & 
                            plottingInfoAll$Outcome == "Correct & N.S.")] = 
  sum(better==T & test==F)
plottingInfoAll$Count[which(plottingInfoAll$MPI == index & 
                            plottingInfoAll$Outcome == "Incorrect & p < 0.05")] = 
  sum(better==F & test==T)
plottingInfoAll$Count[which(plottingInfoAll$MPI == index & 
                            plottingInfoAll$Outcome == "Correct & p < 0.05")] = 
  sum(better==T & test==T)

plottingInfoAll$MPI = as.factor(plottingInfoAll$MPI)
plottingInfoAll$Outcome = factor(plottingInfoAll$Outcome,
                                 levels = c("Correct & p < 0.05",
                                            "Correct & N.S.", 
                                            "Incorrect & N.S.",
                                            "Incorrect & p < 0.05"))
plottingInfoAll$Proportion = (plottingInfoAll$Count/length(test))

these = which(plottingInfoAll$Outcome == "Correct & p < 0.05")
plottingInfoAll$Position[these] = 1.175

those = which(plottingInfoAll$Outcome == "Correct & N.S.")
plottingInfoAll$Position[those] = 1 - plottingInfoAll$Proportion[these] -
  plottingInfoAll$Proportion[those]/2

plottingInfoAll$Position[those[which(plottingInfoAll$Proportion[these] == 0)]] = 1.175

these = which(plottingInfoAll$Outcome == "Incorrect & p < 0.05")
plottingInfoAll$Position[these] = plottingInfoAll$Proportion[these]/2
those = which(plottingInfoAll$Outcome == "Incorrect & N.S.")
plottingInfoAll$Position[those] = plottingInfoAll$Proportion[these] +
  plottingInfoAll$Proportion[those]/2

if (sum(plottingInfoAll$Proportion == 0) > 0){
  final = plottingInfoAll[-which(plottingInfoAll$Proportion == 0),]
} else {final = plottingInfoAll}

ggplot(data = final, aes(x = Proportion, y = MPI, group = Outcome,
                                   fill = Outcome)) + 
  geom_col() + theme_minimal() + 
  scale_fill_manual(values = c("green", "forestgreen", "firebrick4", "red")) +
  labs(x = "Proportion", y = "Model Fit Index", 
       fill = "Result",
       label = NULL) + 
  geom_label(aes(label = Proportion, x = Position)) + lims(x = c(-0.025, 1.25))

Now save workspace for later reference

save.image(file = "C:/Users/DELL/Documents/integration-simulations/modelFitIndex.RData")