Skip to content

Commit

Permalink
MLE calibration for IRTest_Cont
Browse files Browse the repository at this point in the history
  • Loading branch information
@SeewooLi
SeewooLi committed May 31, 2024
1 parent 4035b92 commit 96639b7
Show file tree
Hide file tree
Showing 6 changed files with 224 additions and 60 deletions.
2 changes: 1 addition & 1 deletion NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,12 +55,12 @@ importFrom(stats,dnorm)
importFrom(stats,integrate)
importFrom(stats,median)
importFrom(stats,nlminb)
importFrom(stats,optim)
importFrom(stats,pchisq)
importFrom(stats,qbeta)
importFrom(stats,rbeta)
importFrom(stats,rbinom)
importFrom(stats,rchisq)
importFrom(stats,rgamma)
importFrom(stats,rnorm)
importFrom(stats,runif)
importFrom(stats,xtabs)
Expand Down
3 changes: 2 additions & 1 deletion NEWS.md
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@
# IRTest 2.1.0

* Weighted likelihood estimates (WLE) are available for ability parameter estimation.
* `adaptive_test` function is added for ability parameter estimation with fixed item parameters
* `adaptive_test` function is added for ability parameter estimation with fixed item parameters.
* `IRTest_Cont` utilizes MLE in item parameter estimation for a faster estimation.

# IRTest 2.0.0

Expand Down
4 changes: 2 additions & 2 deletions R/DataGeneration.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
#' @description This function generates an artificial item response dataset allowing various options.
#'
#' @importFrom betafunctions rBeta.4P
#' @importFrom stats dnorm rbinom rchisq rnorm runif rbeta
#' @importFrom stats dnorm rbinom rchisq rnorm runif rbeta rgamma
#'
#' @param seed A numeric value that is used for random sampling.
#' Seed number can guarantee a replicability of the result.
Expand Down Expand Up @@ -300,7 +300,7 @@ DataGeneration <- function(seed=1, N=2000,
)
initialitem_C[,1] <- (a_l+a_u)/2

item_C[,3] <- 10
item_C[,3] <- round(rgamma(nitem_C, shape = 10, scale = 1), 2)
initialitem_C[,3] <- 10

# item responses
Expand Down
21 changes: 15 additions & 6 deletions R/IRTest_Cont.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -155,7 +155,8 @@ IRTest_Cont <- function(data, range = c(-6,6), q = 121, initialitem=NULL,
iter <- iter +1

E <- Estep_Cont(item = initialitem, data = data, range = range, q = q, Xk=Xk, Ak=Ak)
initialitem <- Mstep_Cont(E, initialitem, data)
M1 <- Mstep_Cont(E, initialitem, data)
initialitem <- M1[[1]]

diff <- max(
max(abs(I[,-3]-initialitem[,-3]), na.rm = TRUE),
Expand All @@ -175,7 +176,8 @@ IRTest_Cont <- function(data, range = c(-6,6), q = 121, initialitem=NULL,
iter <- iter +1

E <- Estep_Cont(item=initialitem, data=data, q=q, prob=0.5, d=0, sd_ratio=1, range=range, Xk=Xk, Ak=Ak)
initialitem <- Mstep_Cont(E, initialitem, data)
M1 <- Mstep_Cont(E, initialitem, data)
initialitem <- M1[[1]]

ld_est <- latent_dist_est(method = latent_dist, Xk = E$Xk, posterior = E$fk, range=range)
Xk <- ld_est$Xk
Expand All @@ -198,7 +200,9 @@ IRTest_Cont <- function(data, range = c(-6,6), q = 121, initialitem=NULL,
iter <- iter +1

E <- Estep_Cont(item=initialitem, data=data, q=q, prob=prob, d=d, sd_ratio=sd_ratio, range = range)
initialitem <- Mstep_Cont(E, initialitem, data)
M1 <- Mstep_Cont(E, initialitem, data)
initialitem <- M1[[1]]

M2 <- M2step(E)
prob = M2$prob; d = M2$d; sd_ratio = M2$sd_ratio

Expand All @@ -222,7 +226,8 @@ IRTest_Cont <- function(data, range = c(-6,6), q = 121, initialitem=NULL,

E <- Estep_Cont(item=initialitem, data=data, q=q, prob=0.5, d=0, sd_ratio=1,
range=range, Xk=Xk, Ak=Ak)
initialitem <- Mstep_Cont(E, initialitem, data)
M1 <- Mstep_Cont(E, initialitem, data)
initialitem <- M1[[1]]

ld_est <- latent_dist_est(method = latent_dist, Xk = E$Xk, posterior = E$fk, range=range, bandwidth=bandwidth, N=N, q=q)
Xk <- ld_est$Xk
Expand Down Expand Up @@ -251,7 +256,8 @@ IRTest_Cont <- function(data, range = c(-6,6), q = 121, initialitem=NULL,
iter <- iter +1

E <- Estep_Cont(item=initialitem, data=data, q=q, range=range, Xk=Xk, Ak=Ak)
initialitem <- Mstep_Cont(E, initialitem, data)
M1 <- Mstep_Cont(E, initialitem, data)
initialitem <- M1[[1]]

ld_est <- latent_dist_est(method = latent_dist, Xk = E$Xk, posterior = E$fk, range=range, par=density_par,N=N)
Xk <- ld_est$Xk
Expand All @@ -277,7 +283,8 @@ IRTest_Cont <- function(data, range = c(-6,6), q = 121, initialitem=NULL,

E <- Estep_Cont(item=initialitem, data=data, q=q, prob=0.5, d=0, sd_ratio=1,
range=range, Xk=Xk, Ak=Ak)
initialitem <- Mstep_Cont(E, initialitem, data)
M1 <- Mstep_Cont(E, initialitem, data)
initialitem <- M1[[1]]

ld_est <- latent_dist_est(method = latent_dist, Xk = E$Xk, posterior = E$fk, range=range, par=density_par, N=N)
Xk <- ld_est$Xk
Expand Down Expand Up @@ -312,6 +319,7 @@ IRTest_Cont <- function(data, range = c(-6,6), q = 121, initialitem=NULL,
}
dn <- list(colnames(data),c("a", "b", "nu"))
dimnames(initialitem) <- dn
dimnames(M1[[2]]) <- list(colnames(data),c("a", "b", "log(nu)"))

# preparation for outputs

Expand All @@ -324,6 +332,7 @@ IRTest_Cont <- function(data, range = c(-6,6), q = 121, initialitem=NULL,
logL <- logL + as.numeric(E$fk%*%log(Ak)) - sum(E$Pk*log(E$Pk))
return(structure(
list(par_est=initialitem,
se=M1[[2]],
fk=E$fk,
iter=iter,
quad=Xk,
Expand Down
138 changes: 101 additions & 37 deletions R/non_exporting_functions.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -668,29 +668,56 @@ PDs <- function(probab, param, pmat, pcummat, a_supp, par, tcum){
}
}

#' @importFrom stats optim
#'
Mstep_Cont <- function(E, item, data){
# Mstep_Cont2 <- function(E, item, data){
# item_estimated <- item
# item[,3] <- log(item[,3])
# for(i in 1:nrow(item)){
# item_estimated[i,] <- optim(item[i,], fn = LL_Cont, gr = grad_Cont, theta=E$Xk, data=data[,i], Pk = E$Pk, method = "BFGS")$par
# }
# item_estimated[,3] <- exp(item_estimated[,3])
# return(list(item_estimated,NULL))
# }

Mstep_Cont <- function(E, item, data, threshold = 1e-7, max_iter = 20){
nitem <- nrow(item)
item_estimated <- item
item[,3] <- log(item[,3])
for(i in 1:nrow(item)){
item_estimated[i,] <- optim(item[i,], fn = LL_Cont, gr = grad_Cont, theta=E$Xk, data=data[,i], Pk = E$Pk, method = "BFGS")$par
se <- matrix(nrow = nrow(item), ncol = ncol(item))
for(i in 1:nitem){
par <- item[i,]
par[3] <- log(par[3])
iter <- 0
div <- 3
repeat{
iter <- iter + 1

l1l2 <- cont_L1L2(item = par, Xk = E$Xk, data = data[,i], Pk = E$Pk, fk = E$fk)
diff <- l1l2[[1]]%*%l1l2[[2]]

if(is.infinite(sum(abs(diff)))|is.na(sum(abs(diff)))){
par <- par
} else{
if( sum(abs(diff)) > div){
if(max(abs(diff[-1]))/abs(diff[1])>1000){
par[1] <- -par[1]
} else{
par <- par-diff/2
}
} else {
par <- par-diff
div <- sum(abs(diff))
}
}
if( div <= threshold | iter > max_iter) break
}
par[3] <- exp(par[3])
item_estimated[i,] <- par
se[i,] <- suppressWarnings(sqrt(-diag(l1l2[[2]])))
}
item_estimated[,3] <- exp(item_estimated[,3])
return(item_estimated)

return(list(item_estimated, se))
}

# LL_Cont <- function(item, theta, data, Pk){
# nu <- exp(item[3])
# LL <- 0
# for(i in 1:length(theta)){
# mu <- P(theta = theta[i], a = item[1], b = item[2])
# alpha <- mu*nu
# beta <- nu*(1-mu)
# LL <- LL+sum(Pk[,i]*log(dbeta(data, alpha, beta)), na.rm = TRUE)
# }
# return(-LL)
# }

#' @importFrom stats dbeta
#'
LL_Cont <- function(item, theta, data, Pk){
Expand All @@ -702,24 +729,7 @@ LL_Cont <- function(item, theta, data, Pk){
return(-LL)
}

# grad_Cont <- function(item, theta, data, Pk){
# item[3] <- exp(item[3])
# La <- 0
# Lb <- 0
# Lnu <- 0
# for(i in 1:length(theta)){
# mu <- P(theta = theta[i], a = item[1], b = item[2])
# alpha <- mu*item[3]
# beta <- item[3]*(1-mu)
# v1 <- log(data)-digamma(alpha)
# v2 <- log(1-data)-digamma(beta)
# La <- La + sum(Pk[,i]*(theta[i]-item[2])*item[3]*mu*(1-mu)*(v1-v2), na.rm = TRUE)
# Lb <- Lb + sum(-Pk[,i]*item[1]*item[3]*mu*(1-mu)*(v1-v2), na.rm = TRUE)
# Lnu <- Lnu + sum(Pk[,i]*(mu*v1+(1-mu)*v2+digamma(item[3])), na.rm = TRUE)
# }
#
# return(-c(La, Lb, Lnu))
# }

grad_Cont <- function(item, theta, data, Pk){
nu <- exp(item[3])
mu <- P(theta = rep(theta, each=nrow(Pk)), a = item[1], b = item[2])
Expand All @@ -733,6 +743,60 @@ grad_Cont <- function(item, theta, data, Pk){
return(-c(La, Lb, Lnu))
}

cont_L1L2 <- function(item, Xk, data, Pk, fk){
nu <- exp(item[3])
mu <- P(theta = Xk, a = item[1], b = item[2])
aph <- mu*nu
bt <- nu*(1-mu)
s1 <- as.vector(crossprod(log(data[!is.na(data)]), Pk[!is.na(data),]))
s2 <- as.vector(crossprod(log(1-data[!is.na(data)]), Pk[!is.na(data),]))
La <- sum(
(Xk-item[2])*aph*bt/nu*(s1-s2-fk*(digamma(aph)-digamma(bt))),
na.rm = TRUE)
Lb <- -item[1]*sum(
aph*bt/nu*(s1-s2-fk*(digamma(aph)-digamma(bt))),
na.rm = TRUE)
Lxi <- nu*sum(
fk*digamma(nu)+mu*(s1-fk*digamma(aph))+(1-mu)*(s2-fk*digamma(bt)),
na.rm = TRUE)

Laa <- -sum(
(((Xk-item[2])*aph*bt/nu)^2)*fk*(trigamma(aph)+trigamma(bt)),
na.rm = TRUE
)
Lab <- item[1]*sum(
(Xk-item[2])*((aph*bt/nu)^2)*fk*(trigamma(aph)+trigamma(bt)),
na.rm = TRUE
)
Lbb <- -(item[1]^2)*sum(
((aph*bt/nu)^2)*fk*(trigamma(aph)+trigamma(bt)),
na.rm = TRUE
)
Laxi <- -sum(
(Xk-item[2])*aph*bt*fk*(mu*trigamma(aph)-(1-mu)*trigamma(bt)),
na.rm = TRUE
)
Lbxi <- item[1]*sum(
aph*bt*fk*(mu*trigamma(aph)-(1-mu)*trigamma(bt)),
na.rm = TRUE
)
Lxixi <- (nu^2)*sum(fk)*trigamma(nu) - sum(
fk*((aph^2)*trigamma(aph)+(bt^2)*trigamma(bt)),
na.rm = TRUE)

LL_matrix <- matrix(c(
Laa, Lab, Laxi,
Lab, Lbb, Lbxi,
Laxi, Lbxi, Lxixi
),
nrow = 3,
byrow = TRUE)

return(list(L1 = c(La, Lb, Lxi),
L2 = solve(LL_matrix))
)
}

#################################################################################################################
# M2 step
#################################################################################################################
Expand Down
Loading

0 comments on commit 96639b7

Please sign in to comment.