scripts/Simlation study/results-simulations-sx2.R

In joakimwallmark/PolyOptimalIRT:

library(xtable)


# Natmath --------------------------------------
load(file="data/NatMath_fittedmodel_mirt.RData")
badfits <- which(itemfit(natMathGpcm)[5]<0.05)
mirtSX2 <- itemfit(natMathGpcm)[2]
items <- 32

cyc <- c(5, 10)
basis <- c(15) # 12 becomes default for natMath data
bins <- c(25) # 22 becomes default for natMath data
basis <- c(7, 9, 12) # 12 becomes default for natMath data
bins <- c(15, 22) # 22 becomes default for natMath data

sx2res <- matrix(0, nrow = length(cyc)*length(basis)*length(bins), ncol = items)
sx2res <- data.frame(scenario=rep("scenario", 12), sx2res)
sx2kernres <- matrix(0, nrow = length(cyc)*length(basis)*length(bins), ncol = items)
sx2kernres <- data.frame(scenario=rep("scenario", 12), sx2kernres)
i=1
for (cy in cyc) {
    for (bas in basis) {
        for (bin in bins) {
            load(file = paste("data/sx2/natmath",
                              " cyc", cy,
                              " basis", bas,
                              " bins", bin, ".rda", sep = ""))
            sx2res[i, ] <- c(paste("cyc:", cy," basis:", bas," bins:", bin), sx2)
            load(file = paste("data/sx2/natmath ", "kernel",
                              " cyc", cy,
                              " basis", bas,
                              " bins", bin, ".rda", sep = ""))
            sx2kernres[i, ] <- c(paste("cyc:", cy," basis:", bas," bins:", bin), sx2)
            i=i+1
        }
    }
}
names(sx2res) <- names(sx2kernres) <- c("scenario", 1:items)
sx2res[, 2:(items+1)] <- sapply(sx2res[, 2:(items+1)], as.numeric)
sx2kernres[, 2:(items+1)] <- sapply(sx2kernres[, 2:(items+1)], as.numeric)

summary(sx2res)
summary(sx2kernres)

# small s-x2 diff between kernel or TG theta (3.5 and overall avg is ~83-87)
colMeans(sx2res[, 2:(items+1)]) - colMeans(sx2kernres[, 2:(items+1)])
mean(colMeans(sx2res[, 2:(items+1)]) - colMeans(sx2kernres[, 2:(items+1)]))
mean(mirtSX2[, 1]-colMeans(sx2res[, 2:(items+1)])) # avg of testgardener is higher than mirt
mean(mirtSX2[, 1]-colMeans(sx2kernres[, 2:(items+1)]))
mean(mirtSX2[badfits, 1]-colMeans(sx2res[, badfits])) # avg of testgardener is higher for worst mirt fits
mean(mirtSX2[badfits, 1]-colMeans(sx2kernres[, badfits]))

table(apply(sx2res[, 2:(items+1)], 2, which.min)) # 3, 5, 6, 11 and 12 does best
table(apply(sx2kernres[, 2:(items+1)], 2, which.min)) # 6, 11 and 12
# all have 12 basis functions... bins matters less
sx2res[3, 1]
sx2res[5, 1]
sx2res[6, 1]
sx2res[11, 1]
sx2res[12, 1]
sx2kernres[6, 1]
sx2kernres[11, 1]
sx2kernres[12, 1]

# 3 is best for non-kernel
mean(as.vector(mode = "numeric", sx2res[5, 2:(items+1)]-sx2res[3, 2:(items+1)]))
mean(as.vector(mode = "numeric", sx2res[11, 2:(items+1)]-sx2res[3, 2:(items+1)]))
mean(as.vector(mode = "numeric", sx2res[6, 2:(items+1)]-sx2res[3, 2:(items+1)]))
mean(as.vector(mode = "numeric", sx2res[12, 2:(items+1)]-sx2res[3, 2:(items+1)]))
mean(as.vector(mode = "numeric", mirtSX2[, 1]-sx2res[3, 2:(items+1)])) # best is barely lower
mean(as.vector(mode = "numeric", sx2kernres[11, 2:(items+1)]-sx2kernres[6, 2:(items+1)]))
mean(as.vector(mode = "numeric", sx2kernres[12, 2:(items+1)]-sx2kernres[6, 2:(items+1)]))
mean(as.vector(mode = "numeric", mirtSX2[, 1]-sx2kernres[6, 2:(items+1)])) # best is barely lower

# average s-x2 is lower for 5 cycles 21, 23, 26, 27 are exceptions
colMeans(sx2res[grepl("cyc: 5", sx2res$scenario), 2:(items+1)])-colMeans(sx2res[grepl("cyc: 10", sx2res$scenario), 2:(items+1)])
mean(colMeans(sx2res[grepl("cyc: 5", sx2res$scenario), 2:(items+1)])-colMeans(sx2res[grepl("cyc: 10", sx2res$scenario), 2:(items+1)]))
# similar for kernel density
colMeans(sx2kernres[grepl("cyc: 5", sx2kernres$scenario), 2:(items+1)])-colMeans(sx2kernres[grepl("cyc: 10", sx2kernres$scenario), 2:(items+1)])
mean(colMeans(sx2kernres[grepl("cyc: 5", sx2kernres$scenario), 2:(items+1)])-colMeans(sx2kernres[grepl("cyc: 10", sx2kernres$scenario), 2:(items+1)]))




# Hads anxiety ------------------------------------------------------------
load(file="data/hads_anxiety_fittedmodel_mirt.RData")
mirtSX2 <- itemfit(hadsAnxietyGrm)[2]
badfits <- which(itemfit(hadsAnxietyGrm)[5]<0.05)

items <- 7
cyc <- c(5, 10)
basis <- c(10, 12) # 8 becomes default for hads data
bins <- c(18, 22) # 14 becomes default for 90% of hads data
basis <- c(6, 8, 10) # 8 becomes default for hads data
bins <- c(12, 14) # 14 becomes default for 90% of hads data

sx2res <- matrix(0, nrow = length(cyc)*length(basis)*length(bins), ncol = items)
sx2res <- data.frame(scenario=rep("scenario", 12), sx2res)
sx2kernres <- matrix(0, nrow = length(cyc)*length(basis)*length(bins), ncol = items)
sx2kernres <- data.frame(scenario=rep("scenario", 12), sx2kernres)
i=1
for (cy in cyc) {
    for (bas in basis) {
        for (bin in bins) {
            load(file = paste("data/sx2/hads anxiety cyc", cy,
                              " basis", bas,
                              " bins", bin, ".rda", sep = ""))
            sx2res[i, ] <- c(paste("cyc:", cy," basis:", bas," bins:", bin), sx2)
            load(file = paste("data/sx2/hads anxiety kernel cyc", cy,
                              " basis", bas,
                              " bins", bin, ".rda", sep = ""))
            sx2kernres[i, ] <- c(paste("cyc:", cy," basis:", bas," bins:", bin), sx2)
            i=i+1
        }
    }
}
names(sx2res) <- names(sx2kernres) <- c("scenario", 1:7)
sx2res[, 2:8] <- sapply(sx2res[, 2:8], as.numeric)
sx2kernres[, 2:8] <- sapply(sx2kernres[, 2:8], as.numeric)

summary(sx2res)
summary(sx2kernres)

# kernel density gives better fit for all items
colMeans(sx2res[, 2:(items+1)]) - colMeans(sx2kernres[, 2:(items+1)])
mean(colMeans(sx2res[, 2:(items+1)]) - colMeans(sx2kernres[, 2:(items+1)]))
mean(mirtSX2[, 1]-colMeans(sx2res[, 2:(items+1)])) # avg of testgardener is higher than mirt
mean(mirtSX2[, 1]-colMeans(sx2kernres[, 2:(items+1)]))
mean(mirtSX2[badfits, 1]-colMeans(sx2res[, badfits])) # avg of testgardener is lower for worst mirt fits
mean(mirtSX2[badfits, 1]-colMeans(sx2kernres[, badfits]))

# average s-x2 is lower for 5 cycles
colMeans(sx2res[grepl("cyc: 5", sx2res$scenario), 2:(items+1)])-colMeans(sx2res[grepl("cyc: 10", sx2res$scenario), 2:(items+1)])
mean(colMeans(sx2res[grepl("cyc: 5", sx2res$scenario), 2:(items+1)])-colMeans(sx2res[grepl("cyc: 10", sx2res$scenario), 2:(items+1)]))
# similar for kernel density
colMeans(sx2kernres[grepl("cyc: 5", sx2kernres$scenario), 2:(items+1)])-colMeans(sx2kernres[grepl("cyc: 10", sx2kernres$scenario), 2:(items+1)])
mean(colMeans(sx2kernres[grepl("cyc: 5", sx2kernres$scenario), 2:(items+1)])-colMeans(sx2kernres[grepl("cyc: 10", sx2kernres$scenario), 2:(items+1)]))

table(apply(sx2res[, 2:(items+1)], 2, which.min)) # 11 is best
table(apply(sx2kernres[, 2:(items+1)], 2, which.min)) # 11 is best, 6 second
# 11 uses 12 bins 10 basis 10 cycles (highest basis used)
sx2res[11, 1]
sx2kernres[11, 1]
sx2kernres[6, 1] # 14 bins 5 cycles...

# 3 is best for non-kernel
mean(as.vector(mode = "numeric", sx2res[11, 2:(items+1)]-sx2res[6, 2:(items+1)]))
mean(as.vector(mode = "numeric", mirtSX2[, 1]-sx2res[11, 2:(items+1)])) # mirt fits better
mean(as.vector(mode = "numeric", mirtSX2[, 1]-sx2kernres[11, 2:(items+1)]))



# Hads depression ---------------------------------------------------------
load(file="data/hads_depression_fittedmodel_mirt.RData")
mirtSX2 <- itemfit(hadsDepressionGrm)[2]
badfits <- which(itemfit(hadsDepressionGrm)[5]<0.05) # no bad fits

items <- 7
cyc <- c(10)
basis <- c(10, 12) # 8 becomes default for hads data
bins <- c(18, 22) # 14 becomes default for 90% of hads data
basis <- c(6, 8, 10) # 8 becomes default for hads data
bins <- c(12, 14) # 14 becomes default for 90% of hads data

sx2res <- matrix(0, nrow = length(cyc)*length(basis)*length(bins), ncol = items)
sx2res <- data.frame(scenario=rep("scenario", length(cyc)*length(basis)*length(bins)), sx2res)
sx2kernres <- matrix(0, nrow = length(cyc)*length(basis)*length(bins), ncol = items)
sx2kernres <- data.frame(scenario=rep("scenario", length(cyc)*length(basis)*length(bins)), sx2kernres)
i=1
for (cy in cyc) {
    for (bas in basis) {
        for (bin in bins) {
            load(file = paste("data/sx2/hads depression cyc", cy,
                              " basis", bas,
                              " bins", bin, ".rda", sep = ""))
            sx2res[i, ] <- c(paste("cyc:", cy," basis:", bas," bins:", bin), sx2)
            load(file = paste("data/sx2/hads depression kernel cyc", cy,
                              " basis", bas,
                              " bins", bin, ".rda", sep = ""))
            sx2kernres[i, ] <- c(paste("cyc:", cy," basis:", bas," bins:", bin), sx2)
            i=i+1
        }
    }
}
names(sx2res) <- names(sx2kernres) <- c("scenario", 1:7)
sx2res[, 2:8] <- sapply(sx2res[, 2:8], as.numeric)
sx2kernres[, 2:8] <- sapply(sx2kernres[, 2:8], as.numeric)

summary(sx2res)
summary(sx2kernres)

# similar for both density methods
colMeans(sx2res[, 2:(items+1)]) - colMeans(sx2kernres[, 2:(items+1)])
mean(colMeans(sx2res[, 2:(items+1)]) - colMeans(sx2kernres[, 2:(items+1)]))
mean(mirtSX2[, 1]-colMeans(sx2res[, 2:(items+1)])) # avg of testgardener is higher than mirt
mean(mirtSX2[, 1]-colMeans(sx2kernres[, 2:(items+1)]))
mean(mirtSX2[badfits, 1]-colMeans(sx2res[, badfits])) # no bad fits
mean(mirtSX2[badfits, 1]-colMeans(sx2kernres[, badfits]))

# average s-x2 is lower for 5 cycles
colMeans(sx2res[grepl("cyc: 5", sx2res$scenario), 2:(items+1)])-colMeans(sx2res[grepl("cyc: 10", sx2res$scenario), 2:(items+1)])
mean(colMeans(sx2res[grepl("cyc: 5", sx2res$scenario), 2:(items+1)])-colMeans(sx2res[grepl("cyc: 10", sx2res$scenario), 2:(items+1)]))
# similar for kernel density
colMeans(sx2kernres[grepl("cyc: 5", sx2kernres$scenario), 2:(items+1)])-colMeans(sx2kernres[grepl("cyc: 10", sx2kernres$scenario), 2:(items+1)])
mean(colMeans(sx2kernres[grepl("cyc: 5", sx2kernres$scenario), 2:(items+1)])-colMeans(sx2kernres[grepl("cyc: 10", sx2kernres$scenario), 2:(items+1)]))

table(apply(sx2res[, 2:(items+1)], 2, which.min)) # 5 is best
table(apply(sx2kernres[, 2:(items+1)], 2, which.min)) # 5 is best
# 5 uses 12 bins 10 basis 10 cycles (highest basis used)
sx2res[5, 1]
sx2kernres[5, 1]

# 3 is best for non-kernel
mean(as.vector(mode = "numeric", mirtSX2[, 1]-sx2res[5, 2:(items+1)])) # mirt fits better
mean(as.vector(mode = "numeric", mirtSX2[, 1]-sx2kernres[5, 2:(items+1)]))