context("d2LL")
load("cat_objects.Rdata")
d2LL_test <- function(cat, theta, usePrior) {
answered_questions <- which(!is.na(cat@answers))
prior_shift <- 1 / cat@priorParams[2]^2
if(length(answered_questions) == 0) {
L_theta <- prior_shift
}
sum_this <- rep(0, length(answered_questions))
if(cat@model == "ltm" | cat@model == "tpm"){
for(i in 1:length(answered_questions)){
item <- answered_questions[i]
P <- probability(cat, theta, item)
Q <- 1 - P
sum_this[i] <- cat@discrimination[i]^2 * ((P-cat@guessing[i]) /
(1-cat@guessing[i]))^2 * (Q/P)
}
L_theta <- (-1 * sum(sum_this))
}
if(cat@model == "grm"){
for(i in 1:length(answered_questions)){
item <- answered_questions[i]
answer_k <- cat@answers[item]
probs <- probability(cat, theta, item)
Pstar1 <- probs[answer_k+1]
Qstar1 <- 1-Pstar1
Pstar2 <- probs[answer_k]
Qstar2 <- 1 - Pstar2
P <- Pstar1 - Pstar2
W2 <- Pstar2 * Qstar2
W1 <- Pstar1 * Qstar1
sum_this[i] <- cat@discrimination[i]^2 * ((-W2*(Qstar2-Pstar2)
+ W1*(Qstar1-Pstar1)) /
P - ((W1 - W2)^2/P^2))
}
L_theta <- sum(sum_this)
}
if(cat@model == "gpcm"){
for(i in 1:length(answered_questions)){
## first and second derivatives of probability
item <- answered_questions[i]
discrimination = cat@discrimination[item]
categoryparams = c(0, cat@difficulty[[item]])
f <- f_prime <- f_primeprime <- rep(NA,length(categoryparams))
for(ans in 1:length(categoryparams)){
f[ans] <- exp(sum(discrimination * (theta - categoryparams[1:ans])))
f_prime[ans] <- f[ans] * (discrimination * ans)
f_primeprime[ans] <- f[ans] * (discrimination * ans)^2
}
g <- sum(f)
g_prime <- sum(f * (discrimination * 1:length(categoryparams)))
g_primeprime <- sum(f * (discrimination * 1:length(categoryparams))^2)
p_prime <- (g * f_prime - f * g_prime) / g^2
a <- (g * f_prime - f * g_prime)
a_prime <- f_primeprime * g - g_primeprime * f
b <- g^2
b_prime <- (2 * g) * g_prime
p_primeprime <- (b * a_prime - a * b_prime) / b^2
answer_k <- cat@answers[item]
p <- probability(cat, theta, item)
sum_this[i] <- ((p_prime[answer_k]^2 / p[answer_k]^2)
- (p_primeprime[answer_k] / p[answer_k]))
}
L_theta <- - sum(sum_this)
}
if(usePrior == TRUE){
L_theta <- L_theta - prior_shift
}
return(L_theta)
}
test_that("ltm d2LL calculates correctly", {
ltm_cat@answers[1:5] <- c(0, 1, 0, 0, 1)
package_d2LL <- d2LL(ltm_cat, 1, FALSE)
test_d2LL <- d2LL_test(ltm_cat, 1, FALSE)
expect_equal(package_d2LL, test_d2LL)
})
test_that("grm d2LL calculates correctly", {
grm_cat@answers[1:5] <- c(4, 5, 2, 4, 4)
package_d2LL <- d2LL(grm_cat, 1, FALSE)
test_d2LL <- d2LL_test(grm_cat, 1, FALSE)
expect_equal(package_d2LL, test_d2LL)
})
test_that("gpcm d2LL calculates correctly", {
gpcm_cat@answers[1:5] <- c(4, 5, 2, 4, 4)
package_d2LL <- d2LL(gpcm_cat, 1, FALSE)
test_d2LL <- d2LL_test(gpcm_cat, 1, FALSE)
expect_equal(package_d2LL, test_d2LL)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.