R/CacSSMIRT.R
In liuhuan90123/Reliability: Reliability and CSEM in CTT, IRT, and GT
Defines functions
FxTheta
#' @title TestRelSSMIRT
#'
#' @description
#' A function to calculate test reliability of SS MIRT with 2PL model
#'
#' @param itemPara a data frame or matrix with parameters of sequence b, a1, a2,...,ai on the 1.702 metric
#' @param strat a vector containing number of items for each strat
#' @param cormat a correlation matrix for factors
#' @return test reliability of SS MIRT with 2PL model
#'
#' @author {Huan Liu, University of Iowa, \email{huan-liu-1@@uiowa.edu}}
#' @export
# library(LaplacesDemon)
# library(profvis)
# source("R/NormalQuadraPoints.R")
source("R/LordWingersky.R")
library(mvtnorm)
# TestRelMIRT(itemPara_SS_A, strat, cormat_A, convTable_A, numOfQuad)
# Form A
itemPara_SS_A <- read.table("TestData/SpanishLit_prm_A_SS.txt")[,c(7:10)]
cormat_A <- matrix(c(1, 0.9067069, 0.6994119,
0.9067069, 1, 0.4891160,
0.6994119,0.4891160,1), nrow = 3)
strat <- c(13, 12, 6)
# item parameter transformation
names(itemPara_SS_A) <- c("b", "a1","a2","a3")
itemPara_SS_A$a <- c(itemPara_SS_A$a1[1:13], itemPara_SS_A$a2[14:25], itemPara_SS_A$a3[26:31])
itemPara_SS_A[,"b"] <- -itemPara_SS_A[,"b"]/itemPara_SS_A[,"a"]
itemPara_SS_A[,"a"] <- itemPara_SS_A[,"a"]/1.702
itemPara_SS_A$a1[1:13] <- itemPara_SS_A$a[1:13]
itemPara_SS_A$a2[14:25] <- itemPara_SS_A$a[14:25]
itemPara_SS_A$a3[26:31] <- itemPara_SS_A$a[26:31]
# form B
itemPara_SS_B <- read.table("TestData/SpanishLit_prm_B_SS.txt")[,c(7:10)]
cormat_B <- matrix(c(1, 0.9722234, 0.5602197,
0.9722234, 1, 0.4795721,
0.5602197,0.4795721,1), nrow = 3)
strat <- c(13, 12, 6)
# item parameter transformation
names(itemPara_SS_B) <- c("b", "a1","a2","a3")
itemPara_SS_B$a <- c(itemPara_SS_B$a1[1:13], itemPara_SS_B$a2[14:25], itemPara_SS_B$a3[26:31])
itemPara_SS_B[,"b"] <- -itemPara_SS_B[,"b"]/itemPara_SS_B[,"a"]
itemPara_SS_B[,"a"] <- itemPara_SS_B[,"a"]/1.702
itemPara_SS_B$a1[1:13] <- itemPara_SS_B$a[1:13]
itemPara_SS_B$a2[14:25] <- itemPara_SS_B$a[14:25]
itemPara_SS_B$a3[26:31] <- itemPara_SS_B$a[26:31]
# read conversion tables
convTable_A <- read.csv("TestData/conversion_table_Form A.csv")
convTable_A <- convTable_A[1:32, c("RawScore", "roundedSS")]
names(convTable_A) <- c("y","roundedSS")
convTable_B <- read.csv("TestData/conversion_table_Form B.csv")
convTable_B <- convTable_B[1:32, c("RawScore", "roundedSS")]
names(convTable_B) <- c("y","roundedSS")
# test
itemPara <- itemPara_SS_A
cormat <- cormat_A
convTable <- convTable_A
# TestRelMIRT <- function(itemPara, strat, cormat,convTable, numOfQuad = 11){ #modelType,
# num of quadratures
numOfQuad <- 11
# number of factors
numOfFactors <- length(strat)
# set nodes and weights
nodes <- seq(-4, 4, length.out = numOfQuad)
# if(modelType = "bf"){
# numOfFactors <- numOfFactors + 1
# }
nodesM <- as.matrix(switch(numOfFactors,
{expand.grid(nodes)},
{expand.grid(nodes, nodes)},
{expand.grid(nodes, nodes, nodes)},
{expand.grid(nodes, nodes, nodes, nodes)},
{expand.grid(nodes, nodes, nodes, nodes, nodes)}))
weightsUnwtd <- dmvnorm(nodesM, c(rep(0,numOfFactors)), cormat, log=FALSE) # mvtnorm
nodesM <- as.data.frame(nodesM)
nodesM$weightsWtd <- weightsUnwtd / sum(weightsUnwtd)
# item parameters and quadrature points to P/Q
# num of items
numOfItem <- nrow(itemPara)
# fxtheta distribution function
FxTheta <- function(itemPara){
# itemPara <- itemPara[1:strat[1],c("b", "a")]#test
NormalQuadraPoints <- function(n){
# set nodes ranging from -4 to 4
nodes <- seq(-4, 4, length.out = n)
# unnormalized weights
weightsUnwtd <- sapply(nodes, FUN = function(x) dnorm(x))
# normalized weightes
weightsWtd <- weightsUnwtd / sum(weightsUnwtd)
# return nodes and normalized weights
return(list("nodes" = nodes, "weights" = weightsWtd))
}
# names item parameter
names(itemPara) <- c("b", "a")
# num of items
numOfItem <- nrow(itemPara)
# weights and nodes
quadPoints <- NormalQuadraPoints(numOfQuad)
# replicate item parameter and theta
itemParaRep <-itemPara[rep(seq_len(numOfItem), each = numOfQuad),]
itemParaRep$theta <- rep(quadPoints$nodes, each = 1, length.out = numOfQuad*numOfItem)
# calculate information by theta
itemParaRep <- within(itemParaRep, {
P = 0 + (1 - 0) / (1 + exp(-1.702 * a * (theta - b)))
Q = 1 - P
PQ = P * Q
info = 1.702**2 * a**2 * P * Q
})
# order by theta
itemParaRep <- itemParaRep[order(itemParaRep$theta),]
# define matrix of marginal distribution of theta
fxTheta <- matrix(NA, nrow = numOfQuad, ncol = numOfItem + 1) # 41 num of quadratures, 41 num of raw sxores
# for loop to calculate fxTheta
for (i in 1:numOfQuad){
probs <- matrix(c(itemParaRep[(1 + numOfItem * (i - 1)):(numOfItem * i),]$P),
nrow = numOfItem, ncol = 1, byrow = FALSE)
fxTheta[i, ] <- LordWingersky(probs)$probability
}
# transform to data frame
fxTheta <- as.data.frame(fxTheta)
fxTheta
}
# fxtheta distribution
fxTheta1 <- FxTheta(itemPara[1:strat[1],c("b", "a")])
fxTheta2 <- FxTheta(itemPara[(strat[1]+1):(strat[1]+strat[2]),c("b", "a")])
fxTheta3 <- FxTheta(itemPara[(strat[1]+strat[2]+1):(strat[1]+strat[2]+strat[3]),c("b", "a")])
names(fxTheta1) <- c(0:strat[1])
names(fxTheta2) <- c(0:strat[2])
names(fxTheta3) <- c(0:strat[3])
fxTheta1$Var1 <- seq(-4, 4, length.out = numOfQuad)
nodesM <- merge(x = nodesM, y = fxTheta1, by = "Var1", all.x = TRUE)
fxTheta2$Var2 <- seq(-4, 4, length.out = numOfQuad)
nodesM <- merge(x = nodesM, y = fxTheta2, by = "Var2", all.x = TRUE)
fxTheta3$Var3 <- seq(-4, 4, length.out = numOfQuad)
nodesM <- merge(x = nodesM, y = fxTheta3, by = "Var3", all.x = TRUE)
### f(y|theta) ---------------
for(i in 1:nrow(nodesM)){
# i <- 1 # test
# fx distribution
fx1 <- t(nodesM[i, 5:(5+strat[1])])
fx2 <- t(nodesM[i, (5+strat[1]+1):(5+strat[1]+strat[2]+1)])
fx3 <- t(nodesM[i, (5+strat[1]+strat[2]+1+1):(5+strat[1]+strat[2]+1+strat[3]+1)])
rownames(fx1) <- 0:strat[1]
rownames(fx2) <- 0:strat[2]
rownames(fx3) <- 0:strat[3]
xSum <- expand.grid(rownames(fx1), rownames(fx2), rownames(fx3))
names(xSum) <- c("x1", "x2", "x3")
fxSum <- expand.grid(fx1, fx2, fx3)
names(fxSum) <- c("fx1", "fx2", "fx3")
fxThetaSum <- cbind(fxSum, xSum)
fxThetaSum$x1 <- as.numeric(as.character(fxThetaSum$x1))
fxThetaSum$x2 <- as.numeric(as.character(fxThetaSum$x2))
fxThetaSum$x3 <- as.numeric(as.character(fxThetaSum$x3))
# fy distribution
fxThetaSum$y <- fxThetaSum$x1 + fxThetaSum$x2 + fxThetaSum$x3
fxThetaSum$wty <- fxThetaSum$fx1 * fxThetaSum$fx2 * fxThetaSum$fx3
fy <- fxThetaSum[,c("y", "wty")]
fyDist <- aggregate(fy$wty, by=list(Category=fy$y), FUN=sum)
names(fyDist) <- c("y", "wts")
# weighted mean of Obs Y (true y) and variance of Obs Y
weightedMean <- sum(fyDist$y * fyDist$wts)/sum(fyDist$wts)
# varianceY <- sum(fyDist$wts * (fyDist$y - weightedMean)^2)
# save results
nodesM[i,"weightedMean"] <- weightedMean
# nodesM[i,"varianceY"] <- varianceY
# SS
# fySSDist <- merge(fyDist, convTable, by = "y")
#
# # weighted mean of Obs Y (true y) and variance of Obs Y
# weightedMeanSS <- sum(fySSDist$roundedSS * fySSDist$wts)/sum(fySSDist$wts)
# varianceYSS <- sum(fySSDist$wts * (fySSDist$roundedSS - weightedMeanSS)^2)
#
# # store results
# nodesM[i,"weightedMeanSS"] <- weightedMeanSS
# nodesM[i,"varianceYSS"] <- varianceYSS
nodesM[i,c(40:71)] <- t(fyDist$wts)
}
### cutscore
cutscore <- 21
nn <- 1331
sc <- 32
nc <- length(cutscore)
exp.TS <- nodesM$weightedMean
rec.mat <- nodesM[ ,c(40:71)]
esacc <- escon <-matrix(NA,nc,nn, dimnames = list(paste("cut at",round(cutscore,3)))) ## create matrix to save data
for(j in 1:nc){ # j is index of cutscore
# j <- 1 # test
cuts<-c(0, cutscore[j], sc) # cutscores
categ<-cut(exp.TS,cuts,labels=FALSE, right = FALSE) # compare true score to cutscores using function "cut", 1s<29 and 2s>=29, !! 29 is level 2
# cut(c(25, 29, 30), cuts,labels=FALSE, right = FALSE)
bang<-ceiling(cuts) # integer
rec.s <- list(NA) # create list to save data
for(i in 1:2){
# i <- 1 # test
rec.s[[i]] <- as.matrix(rec.mat[ , (bang[i]+1):bang[i+1]]) # columns 1:29, that is score 0:28, separate the rec.s to two lists
}
for(i in 1:nn){
# i =1 # test
esacc[j,i]<- sum(rec.s[[categ[i]]][i,]) # categ[i] = 1, x < 29; categ[i] = 2, x >= 29
}
escon[j,] <- rowSums(rec.s[[1]])^2 + rowSums(rec.s[[2]])^2 ## row sumSums first for each category, then ^2, then sum over category
}
### weighted by weights, because use D method
### Accuracy
sum(nodesM$weightsWtd * esacc)
### Conssitency
sum(nodesM$weightsWtd * escon)
ans<- (list("Marginal" = cbind("Accuracy" = sum(nodesM$weightsWtd * esacc), "Consistency" = sum(nodesM$weightsWtd * escon)), "Conditional" = list("Accuracy" =t(esacc), "Consistency" = t(escon))))
ans
### Guo's method ----------------- P method, no weights needed, use likelihood heights
cutscore <- 21/31 * (15- (-15)) + -15
5.32
# read raw data
rawData_A <- read.table("TestData/FormA_31_3000.txt")
rawData <- rawData_A
# read theta and se
## MLE ---------------------------------------
# number of factors
numOfFactors <- 3
# Form A
scoSS_MLE <- read.table("TestData/SpanishLit_sco_A_SS_MLE.txt")[,c(4:9, 10, 12, 15)]
cor <- c(0.9067069, # 1&2
0.6994119, # 1&3
0.4891160) # 2&3
# Form B
# scoSS_MLE <- read.table("TestData/SpanishLit_sco_B_SS_MLE.txt")[,c(4:9, 10, 12, 15)]
# cor <- c(0.9722234, # 1&2
# 0.5602197, # 1&3
# 0.4795721) # 2&3
# range(scoSS_MLE$theta1)
# Form Test
# scoSS_MLE <- read.table("TestData/SpanishLit-sco-test.txt")[,c(4:9, 10, 12, 15)]
# cor <- c(0.91, # 1&2
# 0.71, # 1&3
# 0.51) # 2&3
# change variable name
names(scoSS_MLE) <- c("theta1", "theta2", "theta3", "se1", "se2", "se3", "var1", "var2", "var3")
nn <- nrow(scoSS_MLE)
scoSS_MLE$thetaMLE <- scoSS_MLE$theta1 + scoSS_MLE$theta2 + scoSS_MLE$theta3
# item parameter
# Form A
itemPara_SS_A <- read.table("TestData/SpanishLit_prm_A_SS.txt")[,c(7:10)]
cormat_A <- matrix(c(1, 0.9067069, 0.6994119,
0.9067069, 1, 0.4891160,
0.6994119,0.4891160,1), nrow = 3)
strat <- c(13, 12, 6)
# item parameter transformation
names(itemPara_SS_A) <- c("b", "a1","a2","a3")
itemPara_SS_A$a <- c(itemPara_SS_A$a1[1:13], itemPara_SS_A$a2[14:25], itemPara_SS_A$a3[26:31])
itemPara_SS_A[,"b"] <- -itemPara_SS_A[,"b"]/itemPara_SS_A[,"a"]
itemPara_SS_A[,"a"] <- itemPara_SS_A[,"a"]/1.702
itemPara_SS_A$a1[1:13] <- itemPara_SS_A$a[1:13]
itemPara_SS_A$a2[14:25] <- itemPara_SS_A$a[14:25]
itemPara_SS_A$a3[26:31] <- itemPara_SS_A$a[26:31]
# test
itemPara <- itemPara_SS_A
cormat <- cormat_A
convTable <- convTable_A
## multivariate normal distribution, nodes and weights
# num of quadratures
numOfQuad <- 11
# number of factors
numOfFactors <- length(strat)
# set nodes and weights
nodes <- seq(-5, 5, length.out = numOfQuad)
# if(modelType = "bf"){
# numOfFactors <- numOfFactors + 1
# }
nodesM <- as.matrix(switch(numOfFactors,
{expand.grid(nodes)},
{expand.grid(nodes, nodes)},
{expand.grid(nodes, nodes, nodes)},
{expand.grid(nodes, nodes, nodes, nodes)},
{expand.grid(nodes, nodes, nodes, nodes, nodes)}))
weightsUnwtd <- dmvnorm(nodesM, c(rep(0,numOfFactors)), cormat, log=FALSE) # mvtnorm
nodesM <- as.data.frame(nodesM)
nodesM$weightsWtd <- weightsUnwtd / sum(weightsUnwtd)
numOfQuadM <- nrow(nodesM)
itemParaRep <-itemPara[rep(seq_len(numOfItem), each = numOfQuadM),]
itemParaRep$theta1 <- rep(nodesM$Var1, each = 1, length.out = numOfQuadM*numOfItem)
itemParaRep$theta2 <- rep(nodesM$Var2, each = 1, length.out = numOfQuadM*numOfItem)
itemParaRep$theta3 <- rep(nodesM$Var3, each = 1, length.out = numOfQuadM*numOfItem)
itemParaRep <- within(itemParaRep, {
P1 = 1 / (1 + exp(-1.702 * a1 * (theta1 - b)))
Q1 = 1 - P1
P2 = 1 / (1 + exp(-1.702 * a2 * (theta2 - b)))
Q2 = 1 - P2
P3 = 1 / (1 + exp(-1.702 * a3 * (theta3 - b)))
Q3 = 1 - P3
})
# 13*1331, 12*1331, 6*1331
itemParaRep$P <- c(itemParaRep$P1[1:(13*numOfQuadM)],
itemParaRep$P2[(13*numOfQuadM+1):(25*numOfQuadM)],
itemParaRep$P3[(25*numOfQuadM+1):(31*numOfQuadM)])
itemParaRep$Q <- 1 - itemParaRep$P
# create table to save Lri matrix
rec.mat <- matrix(NA, nn, 31) ## create matrix to save data # nn should be number of examinees # numofQuadM is not
## right here, should be changed to the number of sum theta score, 31 here, or 1331 for approach 2
#### approach 1 ------------------
# use responses
for (i in 1: nn){
# i <- 1 # test
itemParaRep$resp <- as.numeric(rep(rawData[i,], each = numOfQuadM, length.out = numOfQuadM*numOfItem))
itemParaRep$respi <- 1 - itemParaRep$resp
itemParaRep$lh <- itemParaRep$resp * itemParaRep$P + itemParaRep$respi * itemParaRep$Q
Lri <- aggregate(itemParaRep$lh, by=list(itemParaRep$theta1, itemParaRep$theta2, itemParaRep$theta3), FUN=prod)
# normalize
Lri$x <- Lri$x/sum(Lri$x)
# summarize
Lri$thetaS <- Lri$Group.1 + Lri$Group.2 + Lri$Group.3
# aggregate
LriAgg <- aggregate(Lri$x, by=list(Lri$thetaS), FUN=sum) ## should not aggregate here, otherwise consistency will be wrong
# rank
# sum(LriAgg$x)
rec.mat[i,] <- t(LriAgg$x)
}
rec.mat_cp <- rec.mat
### weighted by weights
# nodesM$thetaS <- nodesM$Var1 + nodesM$Var2 + nodesM$Var3
#
# weightsSum <- aggregate(nodesM$weightsWtd, by=list(nodesM$thetaS), FUN=sum)
#
#
# sum(rec.mat[2,]) * weightsSum$x)
range(scoSS_MLE$thetaMLE)
esacc <- escon <-matrix(NA,nc,nn) ## create matrix to save data
for(j in 1:nc){ # j is index of cutscore
# j <- 1 # test
cuts<-c(-16, cutscore[j], 16) # cutscores : * change to -16, and +16
categ<-cut(scoSS_MLE$thetaMLE, cuts,labels=FALSE, right = FALSE) # compare true score to cutscores using function "cut", 1s<29 and 2s>=29, !! 29 is level 2
# cut(c(25, 29, 30), cuts,labels=FALSE, right = FALSE)
# bang<-ceiling(cuts) # integer
bang <- c(0, 21, 31) ###
rec.s <- list(NA) # create list to save data
for(i in 1:2){
# i <- 1 # test
rec.s[[i]] <- as.matrix(rec.mat[ , (bang[i]+1):bang[i+1]]) # columns 1:29, that is score 0:28, separate the rec.s to two lists
}
for(i in 1:nn){
# i =1 # test
esacc[j,i]<- sum(rec.s[[categ[i]]][i,]) # categ[i] = 1, x < 29; categ[i] = 2, x >= 29
}
escon[j,] <- rowSums(rec.s[[1]])^2 + rowSums(rec.s[[2]])^2 ## row sumSums first for each category, then ^2, then sum over category
}
rowMeans(esacc)
rowMeans(escon)
# > rowMeans(esacc)
# [1] 0.9051184
# >
# > rowMeans(escon)
# [1] 0.875969
ans<- (list("Marginal" = cbind("Accuracy" = rowMeans(esacc), "Consistency" = rowMeans(escon)), "Conditional" = list("Accuracy" =t(esacc), "Consistency" = t(escon))))
ans
#### approach 2 ------------------
# use responses
for (i in 1: nn){
# i <- 2 # test
itemParaRep$resp <- as.numeric(rep(rawData[i,], each = 1331, length.out = numOfQuadM*numOfItem))
itemParaRep$respi <- 1 - itemParaRep$resp
itemParaRep$lh <- itemParaRep$resp * itemParaRep$P + itemParaRep$respi * itemParaRep$Q
Lri <- aggregate(itemParaRep$lh, by=list(itemParaRep$theta1, itemParaRep$theta2, itemParaRep$theta3), FUN=prod)
# normalize
Lri$x <- Lri$x/sum(Lri$x)
# summarize
Lri$thetaS <- Lri$Group.1 + Lri$Group.2 + Lri$Group.3
# aggregate
# LriAgg <- aggregate(Lri$x, by=list(Lri$thetaS), FUN=sum) ## should not aggregate here, otherwise consistency will be wrong
# order by thetaS
Lri <- Lri[order(Lri$thetaS),]
# sum(LriAgg$x)
rec.mat[i,] <- t(Lri$x)
}
# sum(Lri$thetaS <= 5) 1111, 220
range(scoSS_MLE$thetaMLE)
esacc <- escon <-matrix(NA,nc,nn) ## create matrix to save data
for(j in 1:nc){ # j is index of cutscore
# j <- 1 # test
cuts<-c(-16, cutscore[j], 16) # cutscores : * change to -6, and +6
categ<-cut(scoSS_MLE$thetaMLE, cuts,labels=FALSE, right = FALSE) # compare true score to cutscores using function "cut", 1s<29 and 2s>=29, !! 29 is level 2
# cut(c(25, 29, 30), cuts,labels=FALSE, right = FALSE)
# bang<-ceiling(cuts) # integer
bang <- c(0, 1111, 1331) ###
rec.s <- list(NA) # create list to save data
for(i in 1:2){
# i <- 1 # test
rec.s[[i]] <- as.matrix(rec.mat[ , (bang[i]+1):bang[i+1]]) # columns 1:29, that is score 0:28, separate the rec.s to two lists
}
for(i in 1:nn){
# i =2 # test
esacc[j,i]<- sum(rec.s[[categ[i]]][i,]) # categ[i] = 1, x < 29; categ[i] = 2, x >= 29
}
escon[j,] <- rowSums(rec.s[[1]])^2 + rowSums(rec.s[[2]])^2 ## row sumSums first for each category, then ^2, then sum over category
}
rowMeans(esacc)
rowMeans(escon)
ans<- (list("Marginal" = cbind("Accuracy" = rowMeans(esacc), "Consistency" = rowMeans(escon)), "Conditional" = list("Accuracy" =t(esacc), "Consistency" = t(escon))))
ans
# sum(rec.mat[2, 1:1111])^2 + sum(rec.mat[2, 1112:1331])^2
### Rudner's method ------------------
# library(mvtnorm)
# dmvnorm(x=c(0,0))
# dmvnorm(x=c(0,0), mean=c(1,1))
# sigma <- matrix(c(4,2,2,3), ncol=2)
# x <- rmvnorm(n=500, mean=c(1,2), sigma=sigma)
# colMeans(x)
# var(x)
# x <- rmvnorm(n=500, mean=c(1,2), sigma=sigma, method="chol")
# colMeans(x)
# var(x)
# plot(x)
#
#
# library(tmvtnorm)
# library(matrixcalc)
# library(corpcor)
#
# is.positive.definite(sigma)
#
# make.positive.definite(sigma, tol=1e-3)
# 0.649660 0.604941 1.088144
# 0.422058 0.365954 1.184058
sigma <- matrix(c(0.649660, 0.9067069, 0.6994119,
0.9067069, 0.604941, 0.4891160,
0.6994119,0.4891160,1.088144), nrow = 3)
sigma <- make.positive.definite(sigma, tol=1e-3)
rmvnorm(n=1000, mean = c(-1.880281, -1.148598, 0.42946), sigma=sigma)
# # variance of error
# varianceError <- sum(nodesM$weightsWtd * nodesM$varianceY)
# varianceErrorSS <- sum(nodesM$weightsWtd * nodesM$varianceYSS)
#
#
# # weighted mean of True Y
# weightedMean <- sum(nodesM$weightedMean * nodesM$weightsWtd)/sum(nodesM$weightsWtd)
# weightedMeanSS <- sum(nodesM$weightedMeanSS * nodesM$weightsWtd)/sum(nodesM$weightsWtd)
#
#
# # variance of True Y
# varianceTrueY <- sum(nodesM$weightsWtd * (nodesM$weightedMean - weightedMean)^2)
# varianceTrueYSS <- sum(nodesM$weightsWtd * (nodesM$weightedMeanSS - weightedMeanSS)^2)
#
# # variance of Obs Y
# varianceObsY <- varianceError + varianceTrueY
# varianceObsYSS <- varianceErrorSS + varianceTrueYSS
#
#
# # variance of Obs Y Approach 2
#
# # sum of observed score variance
# fyThetaWeighted <- apply(nodesM[,43:(43 + numOfItem)], 2, function(x) x * nodesM[,"weightsWtd"])
#
# # sum weighted distribution
# fyObsDist <- as.data.frame(matrix(colSums(fyThetaWeighted[,1:(1 + numOfItem)]), nrow = (1 + numOfItem), ncol = 1))
# fyObsDist$y <- c(0:numOfItem) # test
# fyObsDist$roundedSS <- convTable$roundedSS
# names(fyObsDist) <- c("wts", "y","roundedSS")
#
# # weighted mean of Obs Y
# weightedMeanY <- sum(fyObsDist$y * fyObsDist$wts)/sum(fyObsDist$wts)
# weightedMeanSSY <- sum(fyObsDist$roundedSS * fyObsDist$wts)/sum(fyObsDist$wts)
#
# # variance of Obs Y
# varianceObsY2 <- sum(fyObsDist$wts * (fyObsDist$y - weightedMeanY)^2)
# varianceObsYSS2 <- sum(fyObsDist$wts * (fyObsDist$roundedSS - weightedMeanSSY)^2)
#
# fyObsDist <- fyObsDist[,c(2,3,1)]
# names(fyObsDist) <- c("Raw", "Scale", "Fitted_Freq")
#
#
# # MIRT test reliability
# TestRelSSMIRT <- 1 - varianceError/(varianceError + varianceTrueY)
# TestRelSSMIRTSS <- 1 - varianceErrorSS/(varianceErrorSS + varianceTrueYSS)
#
#
# VarandRel <- as.data.frame(matrix(c(varianceError, varianceTrueY, varianceObsY, TestRelSSMIRT,
# varianceErrorSS, varianceTrueYSS, varianceObsYSS, TestRelSSMIRTSS
# )))
# rownames(VarandRel) <- c("Overall error variance for raw scores",
# "True score variance for raw scores",
# "Observed score variance for raw scores",
# "Reliability for raw scores",
# "Overall error variance for scale scores",
# "True score variance for scale scores",
# "Observed score variance for scale scores",
# "Reliability for scale scores")
# colnames(VarandRel) <- "coefficient"
#
#
# conditionalSEMs <- nodesM[,c("Var1", "Var2", "Var3", "weightsWtd","weightedMean", "varianceY", "weightedMeanSS", "varianceYSS")]
# names(conditionalSEMs) <- c("Theta1", "Theta2", "Theta3", "weights","Ex_Raw", "Raw_Variance", "Ex_Scale", "Scale_Variance")
# conditionalSEMs$Raw_CSEM <- sqrt(conditionalSEMs$Raw_Variance)
# conditionalSEMs$Scale_CSEM <- sqrt(conditionalSEMs$Scale_Variance)
# conditionalSEMs <- conditionalSEMs[with(conditionalSEMs, order(Theta1, Theta2, Theta3)), ]
# rownames(conditionalSEMs) <- 1:nrow(conditionalSEMs)
#
# return(list("Fitted Frequencies" = fyObsDist,
# "Variance and Reliability" = VarandRel,
# "Conditional SEMs" = conditionalSEMs[,c("Theta1", "Theta2", "Theta3", "weights",
# "Ex_Raw", "Ex_Scale", "Raw_CSEM", "Scale_CSEM")]))
# }
#
# TestRelMIRT(itemPara_SS_A, strat, cormat_A, convTable_A, numOfQuad)
# TestRelMIRT(itemPara_SS_B, strat, cormat_B, convTable_B, numOfQuad)
x <- matrix(c(1,2, 3, 4), 2, 2)
rowSums(x)^2
# profvis({
# TestRelMIRT(itemPara_SS_A, strat, cormat_A, convTable_A, numOfQuad)
# })
#
# profvis({
# TestRelMIRT(itemPara_SS_A, strat, cormat_A, convTable_A, numOfQuad=41)
# })
# library(xlsx)
# write.xlsx(conditionalSEMs, "conditionalSEMs_SS_B.xlsx")
liuhuan90123/Reliability documentation built on Aug. 28, 2021, 1:49 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 FxTheta
#' @title TestRelSSMIRT
#'
#' @description
#' A function to calculate test reliability of SS MIRT with 2PL model
#'
#' @param itemPara a data frame or matrix with parameters of sequence b, a1, a2,...,ai on the 1.702 metric
#' @param strat a vector containing number of items for each strat
#' @param cormat a correlation matrix for factors
#' @return test reliability of SS MIRT with 2PL model
#'
#' @author {Huan Liu, University of Iowa, \email{huan-liu-1@@uiowa.edu}}
#' @export
# library(LaplacesDemon)
# library(profvis)
# source("R/NormalQuadraPoints.R")
source("R/LordWingersky.R")
library(mvtnorm)
# TestRelMIRT(itemPara_SS_A, strat, cormat_A, convTable_A, numOfQuad)
# Form A
itemPara_SS_A <- read.table("TestData/SpanishLit_prm_A_SS.txt")[,c(7:10)]
cormat_A <- matrix(c(1, 0.9067069, 0.6994119,
0.9067069, 1, 0.4891160,
0.6994119,0.4891160,1), nrow = 3)
strat <- c(13, 12, 6)
# item parameter transformation
names(itemPara_SS_A) <- c("b", "a1","a2","a3")
itemPara_SS_A$a <- c(itemPara_SS_A$a1[1:13], itemPara_SS_A$a2[14:25], itemPara_SS_A$a3[26:31])
itemPara_SS_A[,"b"] <- -itemPara_SS_A[,"b"]/itemPara_SS_A[,"a"]
itemPara_SS_A[,"a"] <- itemPara_SS_A[,"a"]/1.702
itemPara_SS_A$a1[1:13] <- itemPara_SS_A$a[1:13]
itemPara_SS_A$a2[14:25] <- itemPara_SS_A$a[14:25]
itemPara_SS_A$a3[26:31] <- itemPara_SS_A$a[26:31]
# form B
itemPara_SS_B <- read.table("TestData/SpanishLit_prm_B_SS.txt")[,c(7:10)]
cormat_B <- matrix(c(1, 0.9722234, 0.5602197,
0.9722234, 1, 0.4795721,
0.5602197,0.4795721,1), nrow = 3)
strat <- c(13, 12, 6)
# item parameter transformation
names(itemPara_SS_B) <- c("b", "a1","a2","a3")
itemPara_SS_B$a <- c(itemPara_SS_B$a1[1:13], itemPara_SS_B$a2[14:25], itemPara_SS_B$a3[26:31])
itemPara_SS_B[,"b"] <- -itemPara_SS_B[,"b"]/itemPara_SS_B[,"a"]
itemPara_SS_B[,"a"] <- itemPara_SS_B[,"a"]/1.702
itemPara_SS_B$a1[1:13] <- itemPara_SS_B$a[1:13]
itemPara_SS_B$a2[14:25] <- itemPara_SS_B$a[14:25]
itemPara_SS_B$a3[26:31] <- itemPara_SS_B$a[26:31]
# read conversion tables
convTable_A <- read.csv("TestData/conversion_table_Form A.csv")
convTable_A <- convTable_A[1:32, c("RawScore", "roundedSS")]
names(convTable_A) <- c("y","roundedSS")
convTable_B <- read.csv("TestData/conversion_table_Form B.csv")
convTable_B <- convTable_B[1:32, c("RawScore", "roundedSS")]
names(convTable_B) <- c("y","roundedSS")
# test
itemPara <- itemPara_SS_A
cormat <- cormat_A
convTable <- convTable_A
# TestRelMIRT <- function(itemPara, strat, cormat,convTable, numOfQuad = 11){ #modelType,
# num of quadratures
numOfQuad <- 11
# number of factors
numOfFactors <- length(strat)
# set nodes and weights
nodes <- seq(-4, 4, length.out = numOfQuad)
# if(modelType = "bf"){
# numOfFactors <- numOfFactors + 1
# }
nodesM <- as.matrix(switch(numOfFactors,
{expand.grid(nodes)},
{expand.grid(nodes, nodes)},
{expand.grid(nodes, nodes, nodes)},
{expand.grid(nodes, nodes, nodes, nodes)},
{expand.grid(nodes, nodes, nodes, nodes, nodes)}))
weightsUnwtd <- dmvnorm(nodesM, c(rep(0,numOfFactors)), cormat, log=FALSE) # mvtnorm
nodesM <- as.data.frame(nodesM)
nodesM$weightsWtd <- weightsUnwtd / sum(weightsUnwtd)
# item parameters and quadrature points to P/Q
# num of items
numOfItem <- nrow(itemPara)
# fxtheta distribution function
FxTheta <- function(itemPara){
# itemPara <- itemPara[1:strat[1],c("b", "a")]#test
NormalQuadraPoints <- function(n){
# set nodes ranging from -4 to 4
nodes <- seq(-4, 4, length.out = n)
# unnormalized weights
weightsUnwtd <- sapply(nodes, FUN = function(x) dnorm(x))
# normalized weightes
weightsWtd <- weightsUnwtd / sum(weightsUnwtd)
# return nodes and normalized weights
return(list("nodes" = nodes, "weights" = weightsWtd))
}
# names item parameter
names(itemPara) <- c("b", "a")
# num of items
numOfItem <- nrow(itemPara)
# weights and nodes
quadPoints <- NormalQuadraPoints(numOfQuad)
# replicate item parameter and theta
itemParaRep <-itemPara[rep(seq_len(numOfItem), each = numOfQuad),]
itemParaRep$theta <- rep(quadPoints$nodes, each = 1, length.out = numOfQuad*numOfItem)
# calculate information by theta
itemParaRep <- within(itemParaRep, {
P = 0 + (1 - 0) / (1 + exp(-1.702 * a * (theta - b)))
Q = 1 - P
PQ = P * Q
info = 1.702**2 * a**2 * P * Q
})
# order by theta
itemParaRep <- itemParaRep[order(itemParaRep$theta),]
# define matrix of marginal distribution of theta
fxTheta <- matrix(NA, nrow = numOfQuad, ncol = numOfItem + 1) # 41 num of quadratures, 41 num of raw sxores
# for loop to calculate fxTheta
for (i in 1:numOfQuad){
probs <- matrix(c(itemParaRep[(1 + numOfItem * (i - 1)):(numOfItem * i),]$P),
nrow = numOfItem, ncol = 1, byrow = FALSE)
fxTheta[i, ] <- LordWingersky(probs)$probability
}
# transform to data frame
fxTheta <- as.data.frame(fxTheta)
fxTheta
}
# fxtheta distribution
fxTheta1 <- FxTheta(itemPara[1:strat[1],c("b", "a")])
fxTheta2 <- FxTheta(itemPara[(strat[1]+1):(strat[1]+strat[2]),c("b", "a")])
fxTheta3 <- FxTheta(itemPara[(strat[1]+strat[2]+1):(strat[1]+strat[2]+strat[3]),c("b", "a")])
names(fxTheta1) <- c(0:strat[1])
names(fxTheta2) <- c(0:strat[2])
names(fxTheta3) <- c(0:strat[3])
fxTheta1$Var1 <- seq(-4, 4, length.out = numOfQuad)
nodesM <- merge(x = nodesM, y = fxTheta1, by = "Var1", all.x = TRUE)
fxTheta2$Var2 <- seq(-4, 4, length.out = numOfQuad)
nodesM <- merge(x = nodesM, y = fxTheta2, by = "Var2", all.x = TRUE)
fxTheta3$Var3 <- seq(-4, 4, length.out = numOfQuad)
nodesM <- merge(x = nodesM, y = fxTheta3, by = "Var3", all.x = TRUE)
### f(y|theta) ---------------
for(i in 1:nrow(nodesM)){
# i <- 1 # test
# fx distribution
fx1 <- t(nodesM[i, 5:(5+strat[1])])
fx2 <- t(nodesM[i, (5+strat[1]+1):(5+strat[1]+strat[2]+1)])
fx3 <- t(nodesM[i, (5+strat[1]+strat[2]+1+1):(5+strat[1]+strat[2]+1+strat[3]+1)])
rownames(fx1) <- 0:strat[1]
rownames(fx2) <- 0:strat[2]
rownames(fx3) <- 0:strat[3]
xSum <- expand.grid(rownames(fx1), rownames(fx2), rownames(fx3))
names(xSum) <- c("x1", "x2", "x3")
fxSum <- expand.grid(fx1, fx2, fx3)
names(fxSum) <- c("fx1", "fx2", "fx3")
fxThetaSum <- cbind(fxSum, xSum)
fxThetaSum$x1 <- as.numeric(as.character(fxThetaSum$x1))
fxThetaSum$x2 <- as.numeric(as.character(fxThetaSum$x2))
fxThetaSum$x3 <- as.numeric(as.character(fxThetaSum$x3))
# fy distribution
fxThetaSum$y <- fxThetaSum$x1 + fxThetaSum$x2 + fxThetaSum$x3
fxThetaSum$wty <- fxThetaSum$fx1 * fxThetaSum$fx2 * fxThetaSum$fx3
fy <- fxThetaSum[,c("y", "wty")]
fyDist <- aggregate(fy$wty, by=list(Category=fy$y), FUN=sum)
names(fyDist) <- c("y", "wts")
# weighted mean of Obs Y (true y) and variance of Obs Y
weightedMean <- sum(fyDist$y * fyDist$wts)/sum(fyDist$wts)
# varianceY <- sum(fyDist$wts * (fyDist$y - weightedMean)^2)
# save results
nodesM[i,"weightedMean"] <- weightedMean
# nodesM[i,"varianceY"] <- varianceY
# SS
# fySSDist <- merge(fyDist, convTable, by = "y")
#
# # weighted mean of Obs Y (true y) and variance of Obs Y
# weightedMeanSS <- sum(fySSDist$roundedSS * fySSDist$wts)/sum(fySSDist$wts)
# varianceYSS <- sum(fySSDist$wts * (fySSDist$roundedSS - weightedMeanSS)^2)
#
# # store results
# nodesM[i,"weightedMeanSS"] <- weightedMeanSS
# nodesM[i,"varianceYSS"] <- varianceYSS
nodesM[i,c(40:71)] <- t(fyDist$wts)
}
### cutscore
cutscore <- 21
nn <- 1331
sc <- 32
nc <- length(cutscore)
exp.TS <- nodesM$weightedMean
rec.mat <- nodesM[ ,c(40:71)]
esacc <- escon <-matrix(NA,nc,nn, dimnames = list(paste("cut at",round(cutscore,3)))) ## create matrix to save data
for(j in 1:nc){ # j is index of cutscore
# j <- 1 # test
cuts<-c(0, cutscore[j], sc) # cutscores
categ<-cut(exp.TS,cuts,labels=FALSE, right = FALSE) # compare true score to cutscores using function "cut", 1s<29 and 2s>=29, !! 29 is level 2
# cut(c(25, 29, 30), cuts,labels=FALSE, right = FALSE)
bang<-ceiling(cuts) # integer
rec.s <- list(NA) # create list to save data
for(i in 1:2){
# i <- 1 # test
rec.s[[i]] <- as.matrix(rec.mat[ , (bang[i]+1):bang[i+1]]) # columns 1:29, that is score 0:28, separate the rec.s to two lists
}
for(i in 1:nn){
# i =1 # test
esacc[j,i]<- sum(rec.s[[categ[i]]][i,]) # categ[i] = 1, x < 29; categ[i] = 2, x >= 29
}
escon[j,] <- rowSums(rec.s[[1]])^2 + rowSums(rec.s[[2]])^2 ## row sumSums first for each category, then ^2, then sum over category
}
### weighted by weights, because use D method
### Accuracy
sum(nodesM$weightsWtd * esacc)
### Conssitency
sum(nodesM$weightsWtd * escon)
ans<- (list("Marginal" = cbind("Accuracy" = sum(nodesM$weightsWtd * esacc), "Consistency" = sum(nodesM$weightsWtd * escon)), "Conditional" = list("Accuracy" =t(esacc), "Consistency" = t(escon))))
ans
### Guo's method ----------------- P method, no weights needed, use likelihood heights
cutscore <- 21/31 * (15- (-15)) + -15
5.32
# read raw data
rawData_A <- read.table("TestData/FormA_31_3000.txt")
rawData <- rawData_A
# read theta and se
## MLE ---------------------------------------
# number of factors
numOfFactors <- 3
# Form A
scoSS_MLE <- read.table("TestData/SpanishLit_sco_A_SS_MLE.txt")[,c(4:9, 10, 12, 15)]
cor <- c(0.9067069, # 1&2
0.6994119, # 1&3
0.4891160) # 2&3
# Form B
# scoSS_MLE <- read.table("TestData/SpanishLit_sco_B_SS_MLE.txt")[,c(4:9, 10, 12, 15)]
# cor <- c(0.9722234, # 1&2
# 0.5602197, # 1&3
# 0.4795721) # 2&3
# range(scoSS_MLE$theta1)
# Form Test
# scoSS_MLE <- read.table("TestData/SpanishLit-sco-test.txt")[,c(4:9, 10, 12, 15)]
# cor <- c(0.91, # 1&2
# 0.71, # 1&3
# 0.51) # 2&3
# change variable name
names(scoSS_MLE) <- c("theta1", "theta2", "theta3", "se1", "se2", "se3", "var1", "var2", "var3")
nn <- nrow(scoSS_MLE)
scoSS_MLE$thetaMLE <- scoSS_MLE$theta1 + scoSS_MLE$theta2 + scoSS_MLE$theta3
# item parameter
# Form A
itemPara_SS_A <- read.table("TestData/SpanishLit_prm_A_SS.txt")[,c(7:10)]
cormat_A <- matrix(c(1, 0.9067069, 0.6994119,
0.9067069, 1, 0.4891160,
0.6994119,0.4891160,1), nrow = 3)
strat <- c(13, 12, 6)
# item parameter transformation
names(itemPara_SS_A) <- c("b", "a1","a2","a3")
itemPara_SS_A$a <- c(itemPara_SS_A$a1[1:13], itemPara_SS_A$a2[14:25], itemPara_SS_A$a3[26:31])
itemPara_SS_A[,"b"] <- -itemPara_SS_A[,"b"]/itemPara_SS_A[,"a"]
itemPara_SS_A[,"a"] <- itemPara_SS_A[,"a"]/1.702
itemPara_SS_A$a1[1:13] <- itemPara_SS_A$a[1:13]
itemPara_SS_A$a2[14:25] <- itemPara_SS_A$a[14:25]
itemPara_SS_A$a3[26:31] <- itemPara_SS_A$a[26:31]
# test
itemPara <- itemPara_SS_A
cormat <- cormat_A
convTable <- convTable_A
## multivariate normal distribution, nodes and weights
# num of quadratures
numOfQuad <- 11
# number of factors
numOfFactors <- length(strat)
# set nodes and weights
nodes <- seq(-5, 5, length.out = numOfQuad)
# if(modelType = "bf"){
# numOfFactors <- numOfFactors + 1
# }
nodesM <- as.matrix(switch(numOfFactors,
{expand.grid(nodes)},
{expand.grid(nodes, nodes)},
{expand.grid(nodes, nodes, nodes)},
{expand.grid(nodes, nodes, nodes, nodes)},
{expand.grid(nodes, nodes, nodes, nodes, nodes)}))
weightsUnwtd <- dmvnorm(nodesM, c(rep(0,numOfFactors)), cormat, log=FALSE) # mvtnorm
nodesM <- as.data.frame(nodesM)
nodesM$weightsWtd <- weightsUnwtd / sum(weightsUnwtd)
numOfQuadM <- nrow(nodesM)
itemParaRep <-itemPara[rep(seq_len(numOfItem), each = numOfQuadM),]
itemParaRep$theta1 <- rep(nodesM$Var1, each = 1, length.out = numOfQuadM*numOfItem)
itemParaRep$theta2 <- rep(nodesM$Var2, each = 1, length.out = numOfQuadM*numOfItem)
itemParaRep$theta3 <- rep(nodesM$Var3, each = 1, length.out = numOfQuadM*numOfItem)
itemParaRep <- within(itemParaRep, {
P1 = 1 / (1 + exp(-1.702 * a1 * (theta1 - b)))
Q1 = 1 - P1
P2 = 1 / (1 + exp(-1.702 * a2 * (theta2 - b)))
Q2 = 1 - P2
P3 = 1 / (1 + exp(-1.702 * a3 * (theta3 - b)))
Q3 = 1 - P3
})
# 13*1331, 12*1331, 6*1331
itemParaRep$P <- c(itemParaRep$P1[1:(13*numOfQuadM)],
itemParaRep$P2[(13*numOfQuadM+1):(25*numOfQuadM)],
itemParaRep$P3[(25*numOfQuadM+1):(31*numOfQuadM)])
itemParaRep$Q <- 1 - itemParaRep$P
# create table to save Lri matrix
rec.mat <- matrix(NA, nn, 31) ## create matrix to save data # nn should be number of examinees # numofQuadM is not
## right here, should be changed to the number of sum theta score, 31 here, or 1331 for approach 2
#### approach 1 ------------------
# use responses
for (i in 1: nn){
# i <- 1 # test
itemParaRep$resp <- as.numeric(rep(rawData[i,], each = numOfQuadM, length.out = numOfQuadM*numOfItem))
itemParaRep$respi <- 1 - itemParaRep$resp
itemParaRep$lh <- itemParaRep$resp * itemParaRep$P + itemParaRep$respi * itemParaRep$Q
Lri <- aggregate(itemParaRep$lh, by=list(itemParaRep$theta1, itemParaRep$theta2, itemParaRep$theta3), FUN=prod)
# normalize
Lri$x <- Lri$x/sum(Lri$x)
# summarize
Lri$thetaS <- Lri$Group.1 + Lri$Group.2 + Lri$Group.3
# aggregate
LriAgg <- aggregate(Lri$x, by=list(Lri$thetaS), FUN=sum) ## should not aggregate here, otherwise consistency will be wrong
# rank
# sum(LriAgg$x)
rec.mat[i,] <- t(LriAgg$x)
}
rec.mat_cp <- rec.mat
### weighted by weights
# nodesM$thetaS <- nodesM$Var1 + nodesM$Var2 + nodesM$Var3
#
# weightsSum <- aggregate(nodesM$weightsWtd, by=list(nodesM$thetaS), FUN=sum)
#
#
# sum(rec.mat[2,]) * weightsSum$x)
range(scoSS_MLE$thetaMLE)
esacc <- escon <-matrix(NA,nc,nn) ## create matrix to save data
for(j in 1:nc){ # j is index of cutscore
# j <- 1 # test
cuts<-c(-16, cutscore[j], 16) # cutscores : * change to -16, and +16
categ<-cut(scoSS_MLE$thetaMLE, cuts,labels=FALSE, right = FALSE) # compare true score to cutscores using function "cut", 1s<29 and 2s>=29, !! 29 is level 2
# cut(c(25, 29, 30), cuts,labels=FALSE, right = FALSE)
# bang<-ceiling(cuts) # integer
bang <- c(0, 21, 31) ###
rec.s <- list(NA) # create list to save data
for(i in 1:2){
# i <- 1 # test
rec.s[[i]] <- as.matrix(rec.mat[ , (bang[i]+1):bang[i+1]]) # columns 1:29, that is score 0:28, separate the rec.s to two lists
}
for(i in 1:nn){
# i =1 # test
esacc[j,i]<- sum(rec.s[[categ[i]]][i,]) # categ[i] = 1, x < 29; categ[i] = 2, x >= 29
}
escon[j,] <- rowSums(rec.s[[1]])^2 + rowSums(rec.s[[2]])^2 ## row sumSums first for each category, then ^2, then sum over category
}
rowMeans(esacc)
rowMeans(escon)
# > rowMeans(esacc)
# [1] 0.9051184
# >
# > rowMeans(escon)
# [1] 0.875969
ans<- (list("Marginal" = cbind("Accuracy" = rowMeans(esacc), "Consistency" = rowMeans(escon)), "Conditional" = list("Accuracy" =t(esacc), "Consistency" = t(escon))))
ans
#### approach 2 ------------------
# use responses
for (i in 1: nn){
# i <- 2 # test
itemParaRep$resp <- as.numeric(rep(rawData[i,], each = 1331, length.out = numOfQuadM*numOfItem))
itemParaRep$respi <- 1 - itemParaRep$resp
itemParaRep$lh <- itemParaRep$resp * itemParaRep$P + itemParaRep$respi * itemParaRep$Q
Lri <- aggregate(itemParaRep$lh, by=list(itemParaRep$theta1, itemParaRep$theta2, itemParaRep$theta3), FUN=prod)
# normalize
Lri$x <- Lri$x/sum(Lri$x)
# summarize
Lri$thetaS <- Lri$Group.1 + Lri$Group.2 + Lri$Group.3
# aggregate
# LriAgg <- aggregate(Lri$x, by=list(Lri$thetaS), FUN=sum) ## should not aggregate here, otherwise consistency will be wrong
# order by thetaS
Lri <- Lri[order(Lri$thetaS),]
# sum(LriAgg$x)
rec.mat[i,] <- t(Lri$x)
}
# sum(Lri$thetaS <= 5) 1111, 220
range(scoSS_MLE$thetaMLE)
esacc <- escon <-matrix(NA,nc,nn) ## create matrix to save data
for(j in 1:nc){ # j is index of cutscore
# j <- 1 # test
cuts<-c(-16, cutscore[j], 16) # cutscores : * change to -6, and +6
categ<-cut(scoSS_MLE$thetaMLE, cuts,labels=FALSE, right = FALSE) # compare true score to cutscores using function "cut", 1s<29 and 2s>=29, !! 29 is level 2
# cut(c(25, 29, 30), cuts,labels=FALSE, right = FALSE)
# bang<-ceiling(cuts) # integer
bang <- c(0, 1111, 1331) ###
rec.s <- list(NA) # create list to save data
for(i in 1:2){
# i <- 1 # test
rec.s[[i]] <- as.matrix(rec.mat[ , (bang[i]+1):bang[i+1]]) # columns 1:29, that is score 0:28, separate the rec.s to two lists
}
for(i in 1:nn){
# i =2 # test
esacc[j,i]<- sum(rec.s[[categ[i]]][i,]) # categ[i] = 1, x < 29; categ[i] = 2, x >= 29
}
escon[j,] <- rowSums(rec.s[[1]])^2 + rowSums(rec.s[[2]])^2 ## row sumSums first for each category, then ^2, then sum over category
}
rowMeans(esacc)
rowMeans(escon)
ans<- (list("Marginal" = cbind("Accuracy" = rowMeans(esacc), "Consistency" = rowMeans(escon)), "Conditional" = list("Accuracy" =t(esacc), "Consistency" = t(escon))))
ans
# sum(rec.mat[2, 1:1111])^2 + sum(rec.mat[2, 1112:1331])^2
### Rudner's method ------------------
# library(mvtnorm)
# dmvnorm(x=c(0,0))
# dmvnorm(x=c(0,0), mean=c(1,1))
# sigma <- matrix(c(4,2,2,3), ncol=2)
# x <- rmvnorm(n=500, mean=c(1,2), sigma=sigma)
# colMeans(x)
# var(x)
# x <- rmvnorm(n=500, mean=c(1,2), sigma=sigma, method="chol")
# colMeans(x)
# var(x)
# plot(x)
#
#
# library(tmvtnorm)
# library(matrixcalc)
# library(corpcor)
#
# is.positive.definite(sigma)
#
# make.positive.definite(sigma, tol=1e-3)
# 0.649660 0.604941 1.088144
# 0.422058 0.365954 1.184058
sigma <- matrix(c(0.649660, 0.9067069, 0.6994119,
0.9067069, 0.604941, 0.4891160,
0.6994119,0.4891160,1.088144), nrow = 3)
sigma <- make.positive.definite(sigma, tol=1e-3)
rmvnorm(n=1000, mean = c(-1.880281, -1.148598, 0.42946), sigma=sigma)
# # variance of error
# varianceError <- sum(nodesM$weightsWtd * nodesM$varianceY)
# varianceErrorSS <- sum(nodesM$weightsWtd * nodesM$varianceYSS)
#
#
# # weighted mean of True Y
# weightedMean <- sum(nodesM$weightedMean * nodesM$weightsWtd)/sum(nodesM$weightsWtd)
# weightedMeanSS <- sum(nodesM$weightedMeanSS * nodesM$weightsWtd)/sum(nodesM$weightsWtd)
#
#
# # variance of True Y
# varianceTrueY <- sum(nodesM$weightsWtd * (nodesM$weightedMean - weightedMean)^2)
# varianceTrueYSS <- sum(nodesM$weightsWtd * (nodesM$weightedMeanSS - weightedMeanSS)^2)
#
# # variance of Obs Y
# varianceObsY <- varianceError + varianceTrueY
# varianceObsYSS <- varianceErrorSS + varianceTrueYSS
#
#
# # variance of Obs Y Approach 2
#
# # sum of observed score variance
# fyThetaWeighted <- apply(nodesM[,43:(43 + numOfItem)], 2, function(x) x * nodesM[,"weightsWtd"])
#
# # sum weighted distribution
# fyObsDist <- as.data.frame(matrix(colSums(fyThetaWeighted[,1:(1 + numOfItem)]), nrow = (1 + numOfItem), ncol = 1))
# fyObsDist$y <- c(0:numOfItem) # test
# fyObsDist$roundedSS <- convTable$roundedSS
# names(fyObsDist) <- c("wts", "y","roundedSS")
#
# # weighted mean of Obs Y
# weightedMeanY <- sum(fyObsDist$y * fyObsDist$wts)/sum(fyObsDist$wts)
# weightedMeanSSY <- sum(fyObsDist$roundedSS * fyObsDist$wts)/sum(fyObsDist$wts)
#
# # variance of Obs Y
# varianceObsY2 <- sum(fyObsDist$wts * (fyObsDist$y - weightedMeanY)^2)
# varianceObsYSS2 <- sum(fyObsDist$wts * (fyObsDist$roundedSS - weightedMeanSSY)^2)
#
# fyObsDist <- fyObsDist[,c(2,3,1)]
# names(fyObsDist) <- c("Raw", "Scale", "Fitted_Freq")
#
#
# # MIRT test reliability
# TestRelSSMIRT <- 1 - varianceError/(varianceError + varianceTrueY)
# TestRelSSMIRTSS <- 1 - varianceErrorSS/(varianceErrorSS + varianceTrueYSS)
#
#
# VarandRel <- as.data.frame(matrix(c(varianceError, varianceTrueY, varianceObsY, TestRelSSMIRT,
# varianceErrorSS, varianceTrueYSS, varianceObsYSS, TestRelSSMIRTSS
# )))
# rownames(VarandRel) <- c("Overall error variance for raw scores",
# "True score variance for raw scores",
# "Observed score variance for raw scores",
# "Reliability for raw scores",
# "Overall error variance for scale scores",
# "True score variance for scale scores",
# "Observed score variance for scale scores",
# "Reliability for scale scores")
# colnames(VarandRel) <- "coefficient"
#
#
# conditionalSEMs <- nodesM[,c("Var1", "Var2", "Var3", "weightsWtd","weightedMean", "varianceY", "weightedMeanSS", "varianceYSS")]
# names(conditionalSEMs) <- c("Theta1", "Theta2", "Theta3", "weights","Ex_Raw", "Raw_Variance", "Ex_Scale", "Scale_Variance")
# conditionalSEMs$Raw_CSEM <- sqrt(conditionalSEMs$Raw_Variance)
# conditionalSEMs$Scale_CSEM <- sqrt(conditionalSEMs$Scale_Variance)
# conditionalSEMs <- conditionalSEMs[with(conditionalSEMs, order(Theta1, Theta2, Theta3)), ]
# rownames(conditionalSEMs) <- 1:nrow(conditionalSEMs)
#
# return(list("Fitted Frequencies" = fyObsDist,
# "Variance and Reliability" = VarandRel,
# "Conditional SEMs" = conditionalSEMs[,c("Theta1", "Theta2", "Theta3", "weights",
# "Ex_Raw", "Ex_Scale", "Raw_CSEM", "Scale_CSEM")]))
# }
#
# TestRelMIRT(itemPara_SS_A, strat, cormat_A, convTable_A, numOfQuad)
# TestRelMIRT(itemPara_SS_B, strat, cormat_B, convTable_B, numOfQuad)
x <- matrix(c(1,2, 3, 4), 2, 2)
rowSums(x)^2
# profvis({
# TestRelMIRT(itemPara_SS_A, strat, cormat_A, convTable_A, numOfQuad)
# })
#
# profvis({
# TestRelMIRT(itemPara_SS_A, strat, cormat_A, convTable_A, numOfQuad=41)
# })
# library(xlsx)
# write.xlsx(conditionalSEMs, "conditionalSEMs_SS_B.xlsx")
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.