Nothing
#' Log-normal response time modelling
#'
#' @param RT
#' a Person-x-Item matrix of log-response times (time spent on solving an item).
#' @param data
#' either a list or a simLNIRT object containing the response time matrix.
#' If a simLNIRT object is provided, in the summary the simulated time parameters are shown alongside of the estimates.
#' If the RT variable cannot be found in the list, or if no data object is given, then the RT variable is taken
#' from the environment from which LNRT is called.
#' @param XG
#' the number of MCMC iterations to perform (default: 1000).
#' @param burnin
#' the percentage of MCMC iterations to discard as burn-in period (default: 10).
#' @param XGresid
#' the number of MCMC iterations to perform before residuals are computed (default: 1000).
#' @param residual
#' compute residuals, >1000 iterations are recommended (default: false).
#' @param td
#' estimate the time-discrimination parameter (default: true).
#' @param WL
#' define the time-discrimination parameter as measurement error variance parameter (default: false).
#' @param XPT
#' an optional matrix of predictors for the person speed parameters.
#' @param XIT
#' an optional matrix of predictors for the item time intensity parameters.
#'
#' @return
#' an object of class LNRT.
#'
#' @examples
#' \dontrun{
#' # Log-normal response time modelling
#' data <- simLNIRT(N = 500, K = 20, rho = 0.8, WL = FALSE)
#' out <- LNRT(RT = RT, data = data, XG = 1500, residual = TRUE, td = TRUE, WL = FALSE)
#' summary(out) # Print results
#' out$Post.Means$Time.Intensity # Extract posterior mean estimates
#'
#' library(coda)
#' mcmc.object <- as.mcmc(out$MCMC.Samples$Time.Intensity) # Extract MCMC samples for coda
#' summary(mcmc.object)
#' plot(mcmc.object)
#' }
#' @export
LNRT <-
function(RT,
data,
XG = 1000,
burnin = 10,
XGresid = 1000,
residual = FALSE,
td = TRUE,
WL = FALSE,
XPT = NULL,
XIT = NULL) {
## ident = 1: Identification : fix mean item difficulty(intensity) and product item (time) discrimination responses and response times
## ident = 2: Identification : fix mean ability and speed and product item discrimination responses and response times
# ident <- 1
ident <- 2 # (to investigate person fit using latent scores)
if (XG <= 0) {
stop("XG must be > 0.")
}
if ((burnin <= 0) || (burnin >= 100)) {
stop("burn-in period must be between 0% and 100%.")
}
if (residual && (XGresid >= XG || XGresid <= 0)) {
warning("XGresid must be < XG and > 0. Residuals will not be computed.")
residual <- FALSE
}
if (!missing(data) && !is.null(data)) {
# Try to find RT in the data set first
tryCatch(
RT <- eval(substitute(RT), data),
error = function(e)
NULL
)
# Try to find predictors in the data set first
tryCatch(
XPT <- eval(substitute(XPT), data),
error = function(e)
NULL
)
tryCatch(
XIT <- eval(substitute(XIT), data),
error = function(e)
NULL
)
} else {
data <- NULL
}
RT <- as.matrix(RT)
N <- nrow(RT)
K <- ncol(RT)
if (!is.null(XPT)) {
XPT <- as.matrix(XPT)
if (nrow(XPT) != N) {
stop("nrow(XPT) must be equal to the number of persons.")
}
}
if (!is.null(XIT)) {
XIT <- as.matrix(XIT)
if (nrow(XIT) != K) {
stop("nrow(XIT) must be equal to the number of items.")
}
}
if (WL) {
td <- TRUE #WL <- 1 #time discrimination = 1/sqrt(error variance)
}
cat (" \n")
cat (" LNIRT v", packageDescription("LNIRT")$Version, "\n", sep = "")
cat (" ", rep('-', 20), "\n\n", sep = "")
cat (
" * MCMC sampler initialized (XG:",
XG,
", Burnin:",
paste(burnin, "%", sep = ""),
")\n",
sep = ""
)
cat (" * Response time matrix loaded (", N, "x", K, ") \n\n", sep = "")
# Initialize progress bar
cat (" MCMC progress: \n")
pb <-
txtProgressBar(
min = 1,
max = XG,
initial = 1,
style = 3,
width = 45,
char = "="
)
## Predictors person
if (is.null(XPT)) {
nopredictorp <- TRUE
XPT <-
matrix(1, ncol = 1, nrow = N) #default intercept for speed
}
else {
nopredictorp <- FALSE
if (is.null(XPT)) {
XPT <- matrix(1, ncol = 1, nrow = N) #default intercept for speed
}
}
MmuP <- matrix(0, nrow = XG, ncol = ncol(XPT))
## Predictors item
if (is.null(XIT)) {
nopredictori <- TRUE
MmuI <-
matrix(0, nrow = XG, ncol = 2) # time discrimination + time intensity
XIT <- matrix(NA, nrow = K, ncol = 0)
}
else {
nopredictori <- FALSE
if (!is.null(XIT)) {
if (ident == 2) {
XIT <- cbind(rep(1, K), XIT)
}
}
else {
XIT <- matrix(1, ncol = 1, nrow = K) #default intercept for speed
}
MmuI <-
matrix(0, nrow = XG, ncol = 1 + ncol(XIT)) # time discrimination + XIT
}
## population theta (ability - speed)
theta <- matrix(rnorm(N))
muP <- matrix(0, N, 1)
SigmaP <- diag(1)
muP0 <- matrix(0, 1, 1)
SigmaP0 <- diag(1) / 100
ingroup <- rep(1, N)
Mingroup <- matrix(0, ncol = 2, nrow = N)
if (XG > XGresid) {
flagged <- matrix(0, ncol = 1, nrow = N)
ss <- 1
}
## population item (ability - speed)
ab <- matrix(rnorm(K * 2), ncol = 2)
ab[, 1] <- 1
muI <- t(matrix(rep(c(1, 0), K), ncol = K))
muI0 <- muI[1,]
SigmaI <- diag(2)
SigmaI0 <- diag(2) * 10
for (ii in 1:2) {
SigmaI0[ii,] <- SigmaI0[ii,] * c(0.5, 3)
}
## storage
MT <- MT2 <- array(0, dim = c(N))
MAB <- array(0, dim = c(XG, K, 2))
#MmuP <- array(0, dim = c(XG, 1))
#MmuI <- array(0, dim = c(XG, 2))
MSP <- array(0, dim = c(XG, 1, 1))
MSI <- array(0, dim = c(XG, 2, 2))
sigma2 <- rep(1, K)
Msigma2 <- matrix(0, ncol = K, nrow = XG)
lZP <- 0
lZPT <- 0
lZI <- 0
EAPresid <- matrix(0, ncol = K, nrow = N)
EAPKS <- matrix(0, ncol = 1, nrow = K)
iis <- 1
EAPphi <- matrix(0, ncol = 1, nrow = K)
EAPlambda <- matrix(0, ncol = 1, nrow = K)
EAPtheta <- matrix(0, ncol = 1, nrow = N)
EAPsigma2 <- matrix(0, ncol = 1, nrow = K)
DT <- matrix(1, ncol = 1, nrow = N * K)
DT[which(is.na(RT))] <- 0
DT <- matrix(DT, nrow = N, ncol = K)
EAPCP <- matrix(0, ncol = 1, nrow = N)
# Output
MCMC.Samples <- list()
MCMC.Samples$Person.Speed <- matrix(NA, nrow = XG, ncol = N)
## Start MCMC algorithm
for (ii in 1:XG) {
if (sum(DT == 0) > 0) {
if (WL) {
RT <-
SimulateRT(
RT = RT,
zeta = theta,
lambda = ab[, 2],
phi = rep(1, K),
sigma2 = sigma2,
DT = DT
)
} else {
RT <-
SimulateRT(
RT = RT,
zeta = theta,
lambda = ab[, 2],
phi = ab[, 1],
sigma2 = sigma2,
DT = DT
)
}
}
theta <-
DrawZeta(RT, ab[, 1], ab[, 2], sigma2, muP, SigmaP[1, 1])
theta[1:N] <- theta[1:N] - mean(theta)
MCMC.Samples$Person.Speed[ii,] <- theta
MT[1:N] <- MT[1:N] + theta[1:N]
MT2[1:N] <- MT2[1:N] + theta[1:N] ^ 2
if ((WL)) {
# no time discrimination, 1/(sqrt(error variance)) = discrimination on MVN prior
ab[, 1] <- 1
ab1 <- cbind(1 / sqrt(sigma2), ab[, 2])
dum <-
Conditional(
kk = 2,
Mu = muI,
Sigma = SigmaI,
Z = ab1
)
ab[, 2] <-
DrawLambda_LNRT(
RT = RT,
zeta = theta,
sigma2 = sigma2,
mu = dum$CMU,
sigma = dum$CVAR[1, 1]
)
} else {
if (td) {
dum <- DrawLambdaPhi_LNRT(RT, theta, sigma2, muI, SigmaI, ingroup)
ab[, 1] <- dum$phi
ab[, 1] <- ab[, 1] / (prod(ab[, 1]) ^ (1 / K))
ab[, 2] <- dum$lambda
} else {
ab[, 1] <- rep(1, K)
ab[, 2] <-
DrawLambda_LNRT(
RT = RT,
zeta = theta,
sigma2 = sigma2,
mu = muI[1, 2],
sigma = SigmaI[2, 2]
)
}
}
MAB[ii, 1:K, 1:2] <- ab
sigma2 <- SampleS_LNRT(RT, theta, ab[, 2], ab[, 1], ingroup)
if (WL) {
Msigma2[ii, 1:K] <- 1 / sqrt(sigma2)
} else {
Msigma2[ii, 1:K] <- sigma2
}
# X <- matrix(1, N, 1)
# muP <- SampleB(theta, X, SigmaP, muP0, SigmaP0)
# MmuP[ii, ] <- muP$B
# muP <- muP$pred
## Predictors persons
if (nopredictorp) {
# Population mean estimate for person ability and speed
X <- matrix(1, N, 1)
muP <-
SampleB(
Y = theta,
X = X,
Sigma = SigmaP,
B0 = muP0,
V0 = SigmaP0
)
MmuP[ii,] <- muP$B
muP <- muP$pred
} else{
muPP <- SampleBX_LNRT(Y = theta, XPT = XPT)
MmuP[ii,] <- muPP$B
muP <- muPP$pred
}
# Covariance matrix person parameters
SS <- crossprod(theta - muP) + SigmaP0
SigmaP <- rwishart(1 + N, chol2inv(chol(SS)))$IW
MSP[ii, ,] <- SigmaP
X <- matrix(1, K, 1)
if (WL) {
ab1 <- cbind(1 / sqrt(sigma2), ab[, 2])
} else {
ab1 <- ab
}
# muI2 <- SampleB(ab1, X, SigmaI, muI0, SigmaI0)
# MmuI[ii, 1] <- muI2$B[1]
# MmuI[ii, 2] <- muI2$B[2]
# muI[, 1] <- muI2$pred[, 1]
# muI[, 2] <- muI2$pred[, 2]
#
# SS <- crossprod(ab1 - muI) + SigmaI0
# SigmaI <- rwishart(2 + K, chol2inv(chol(SS)))$IW
# MSI[ii, , ] <- SigmaI
## Predictors items
if (nopredictori) {
# Population mean estimates for item parameters
## Adjust sampling mean item parameters
if ((!WL) && (!td)) {
meanmuI2 <-
(sum(ab[, 2]) / SigmaI[2, 2] + muI0[2] / SigmaI0[2, 2]) / (1 / SigmaI0[2, 2] +
1 / (SigmaI[2, 2]))
sdmuI2 <- sqrt(1 / (1 / SigmaI0[2, 2] + K / (SigmaI[2, 2])))
muI2 <- rnorm(1, mean = meanmuI2, sd = sdmuI2)
MmuI[ii, 1] <- 1
MmuI[ii, 2] <- muI2
muI[, 1] <- 1
muI[, 2] <- muI2
} else{
muI2 <- SampleB(ab1, X, SigmaI, muI0, SigmaI0)
MmuI[ii, 1] <- muI2$B[1]
MmuI[ii, 2] <- muI2$B[2]
muI[, 1] <- muI2$pred[, 1]
muI[, 2] <- muI2$pred[, 2]
}
} else {
#### XIT should include intercept when ident=2
set2 <- 2 # time intensity
ab2 <- matrix(ab1[, set2], ncol = 1)
dum <- SampleBX_LNRT(Y = ab2, XPT = XIT)
MmuI[ii, 2:(ncol(XIT) + 1)] <- dum$B
muI[, set2] <- dum$pred
#mean discrimination and time discrimination
if ((!WL) && (!td)) {
set2 <- 1 # time discrimination
MmuI[ii, 1] <- 1
muI[, 1] <- 1
}
else {
set2 <- 1 # time discrimination
SigmaI1 <- SigmaI[set2, set2]
muI01 <- muI0[set2]
SigmaI01 <- SigmaI0[set2, set2]
ab2 <- matrix(ab1[, set2], ncol = 1)
muI2 <-
SampleB(
Y = ab2,
X = X,
Sigma = SigmaI1,
B0 = muI01,
V0 = SigmaI01
)
MmuI[ii, 1] <- muI2$B
muI[, set2] <- muI2$pred
}
}
## Adjust sampling covariance matrix item parameters
if (!td) {
SS <- sum((ab1[, 2] - muI[, 2]) ** 2) + SigmaI0[2, 2]
SigmaI[2, 2] <- SS / rgamma(1, (K + 1) / 2, 1 / 2) #change sampling
SigmaI[1, 1] <- 1
MSI[ii, , ] <- SigmaI
} else{
SS <- crossprod(ab1 - muI) + SigmaI0
SigmaI <- rwishart(2 + K, chol2inv(chol(SS)))$IW
MSI[ii, , ] <- SigmaI
}
if (ii > XGresid && residual) {
EAPphi <- (ab[, 1] + (iis - 1) * EAPphi) / iis
EAPlambda <- (ab[, 2] + (iis - 1) * EAPlambda) / iis
EAPtheta <- (theta + (iis - 1) * EAPtheta) / iis
EAPsigma2 <- (sigma2 + (iis - 1) * EAPsigma2) / iis
dum <-
personfitLN(
RT = RT,
theta = theta,
phi = ab[, 1],
lambda = ab[, 2],
sigma2 = sigma2
)
lZP <- lZP + dum$lZP
lZPT <- lZPT + dum$lZPT
CF <-
ifelse(dum$lZP < 0.05, 1, 0) #significance level = .05
EAPCP <- (CF + (iis - 1) * EAPCP) / iis
dum <-
itemfitLN(
RT = RT,
theta = theta,
phi = ab[, 1],
lambda = ab[, 2],
sigma2 = sigma2
)
lZI <- lZI + dum$lZI
dum <-
residualLN(
RT = RT,
theta = theta,
phi = ab[, 1],
lambda = ab[, 2],
sigma2 = sigma2,
EAPtheta = EAPtheta,
EAPlambda = EAPlambda,
EAPphi = EAPphi,
EAPsigma2 = EAPsigma2
)
EAPresid <- EAPresid + dum$presid
EAPKS <- (dum$KS[1:K, 1] + (iis - 1) * EAPKS) / iis
iis <- iis + 1
}
# Update progress bar
setTxtProgressBar(pb, ii)
# if (ii%%100 == 0)
# cat("Iteration ", ii, " ", "\n")
# flush.console()
}
MT <- MT / XG
MT2 <- sqrt(MT2 / XG - MT ^ 2)
if (ii > XGresid && residual) {
lZP <- lZP / (XG - XGresid)
lZPT <- lZPT / (XG - XGresid)
lZI <- lZI / (XG - XGresid)
EAPresid <- EAPresid / (XG - XGresid)
}
kit <- 0
if (ncol(XIT) > 0)
kit <- ncol(XIT) - 1
# Output
#MCMC.Samples <- list()
#MCMC.Samples$Person.Speed <- matrix(NA, nrow = XG, ncol = N)
MCMC.Samples$Mu.Person.Speed <- MmuP
MCMC.Samples$Var.Person.Speed <- MSP
MCMC.Samples$Time.Discrimination <- MAB[, , 1]
MCMC.Samples$Time.Intensity <- MAB[, , 2]
MCMC.Samples$Mu.Time.Discrimination <- MmuI[, 1]
MCMC.Samples$Mu.Time.Intensity <- MmuI[, 2:(ncol(MmuI))]
MCMC.Samples$Sigma2 <- Msigma2
MCMC.Samples$CovMat.Item <- MSI
XGburnin <- round(XG * burnin / 100, 0)
Post.Means <- list()
Post.Means$Person.Speed <- MT
if (ncol(XPT) == 1)
Post.Means$Mu.Person.Speed <- mean(MmuP[XGburnin:XG, ])
else
Post.Means$Mu.Person.Speed <- colMeans(MmuP[XGburnin:XG, ])
Post.Means$Var.Person.Speed <- mean(MSP[XGburnin:XG, , ])
Post.Means$Time.Discrimination <- colMeans(MAB[XGburnin:XG, , 1])
Post.Means$Time.Intensity <- colMeans(MAB[XGburnin:XG, , 2])
Post.Means$Mu.Time.Discrimination <- mean(MmuI[XGburnin:XG, 1])
if (kit == 0)
Post.Means$Mu.Time.Intensity <-
mean(MmuI[XGburnin:XG, 2:(ncol(MmuI))])
else
Post.Means$Mu.Time.Intensity <-
colMeans(MmuI[XGburnin:XG, 2:(ncol(MmuI))])
Post.Means$Sigma2 <- colMeans(Msigma2[XGburnin:XG,])
Post.Means$CovMat.Item <-
c(round(apply(MSI[XGburnin:XG, , 1], 2, mean), 3), round(apply(MSI[XGburnin:XG, , 2], 2, mean), 3))
if (!(is(data, "simLNIRT"))) {
data <- NULL # only attach sim data for summary function
}
if (XG > XGresid && residual) {
out <-
list(
Post.Means = Post.Means,
MCMC.Samples = MCMC.Samples,
Mtheta = MT,
MTSD = MT2,
MAB = MAB,
MmuP = MmuP,
MSP = MSP,
MmuI = MmuI,
MSI = MSI,
lZP = lZP,
lZPT = lZPT,
Msigma2 = Msigma2,
theta = theta,
sigma2 = sigma2,
lZI = lZI,
EAPresid = EAPresid,
EAPKS = EAPKS,
RT = RT,
EAPCP = EAPCP,
td = td,
WL = WL,
data = data,
XPT = XPT,
XIT = XIT,
XG = XG,
burnin = burnin,
ident = ident,
residual = residual,
XGresid = XGresid
)
} else {
out <-
list(
Post.Means = Post.Means,
MCMC.Samples = MCMC.Samples,
Mtheta = MT,
MTSD = MT2,
MAB = MAB,
MmuP = MmuP,
MSP = MSP,
MmuI = MmuI,
MSI = MSI,
Msigma2 = Msigma2,
theta = theta,
sigma2 = sigma2,
RT = RT,
td = td,
WL = WL,
data = data,
XPT = XPT,
XIT = XIT,
XG = XG,
burnin = burnin,
ident = ident,
residual = residual,
XGresid = XGresid
)
}
cat("\n\n")
class(out) <- c("LNRT", "list")
return(out)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.