Nothing
library(testthat)
mcmc.summ <- function(x) {
c(m = mean(x), sd = sd(x), quantile(x, c(.025, .975)))
}
################################ dichotomous
N <- 5
M <- 3
# Rasch / 1PL predictions (special case of ADISOP / ACM)
rasch <- outer(
seq(-2, 2, length.out = N), -seq(-1, 2, length.out = M),
function(x, y) (1 + exp(-x + y))^-1
)
test_that("NIRT axioms match with Rasch predictions", {
skip_on_cran()
expect_silent(tmp <- nirt_to_Ab(N, M))
# constraints: 0<theta<1
A <- tmp$A
b <- tmp$b
expect_equal(nrow(A), length(b))
expect_true(inside(c(rasch), A, b))
n <- 500
k <- rpbinom(c(rasch), n)
# k <- rbinom(M*N, n, .5)
# CRAN issue with starting values (2022-11-21):
start <- find_inside(A, b)
cpost <- count_binom(k, n, A, b,
M = 5000, progress = FALSE,
steps = seq(2, nrow(A) - 1, 2), start = start
)
cprior <- count_binom(0, 0, A, b,
M = 5000, progress = FALSE,
steps = seq(2, nrow(A) - 1, 2), start = start
)
count_to_bf(cpost, cprior)
# multinomineq:::count_auto_bin(k, n, A, b, M=5000, eps = .05)
expect_silent(tt <- sampling_binom(k, n, A, b,
start = start,
M = 2000, progress = FALSE
))
plot(tt[, 1], ty = "l")
expect_silent(ppp <- ppp_binom(tt, k, n))
expect_gt(ppp[3], .01)
#### prior distribution:
# tt <- sampling_binom(0,0, A, b, M = 2000, prior = rep(1/(N*M), 2))
# matrix(apply(tt, 2, mean), N)
################### strength of DC constraint:
skip("functions not yet fully implemented")
DC <- nirt_to_Ab(N, M, axioms = c("W1", "W2", "DC"))
# count_binom(0,0, DC$exclude[[1]]$A, rep(0,3), M=5e4)
# relative to W1 / W2
ss <- sampling_binom(0, 0, DC$A, DC$b, M = 3e4)
dc <- sapply(DC$exclude, function(ee) inside(ss, ee$A, ee$b)) # violation DC
table(rowSums(dc))
mean(apply(dc, 1, any))
spost <- sampling_binom(k, n, DC$A, DC$b, M = 3e4)
dc.post <- sapply(DC$exclude, function(ee) inside(spost, ee$A, ee$b))
mean(apply(dc.post, 1, any))
spost.dc <- spost[!apply(dc.post, 1, any), ]
matrix(apply(spost.dc, 2, mean), N)
matrix(apply(spost, 2, mean), N)
})
#############################################################
# example in Karabatsos 2001, Figure 4 p. 409:
# Perline, Wright, and Wainer (1979), Table 2, p. 244.
ppw <- matrix(c(
11, 20, 37,
44, 43, 54,
60, 56, 56
), 3, 3, byrow = TRUE)
N <- matrix(c(61, 84, 82), 3, 3)
dimnames(ppw) <- dimnames(N) <- list(
score = paste0("s", 3:5),
item = paste0("i", c(4, 9, 2))
)
p.obs <- round(ppw / N, 2)
test_that("Examples in Karabatsos (2004, Figure 4) give identical results", {
skip("functions not yet fully implemented")
IRT_W1 <- nirt_to_Ab(3, 3, axioms = "W1")
pp <- sampling_binom(c(ppw), c(N), IRT_W1$A, IRT_W1$b,
M = 10000, burnin = 1000
)
summ <- t(apply(pp, 2, mcmc.summ))
irt_w1 <- c(.23, .52, .73, .33, .51, .68, .55, .64, .71)
round(cbind(p.obs = c(p.obs), irt_w1, summ), 2)
# expect_equal(unname(summ[,1]), irt_w1, tol = .01) ### TODO
### accept-reject sampling: same results
# M <- 2e5
# bb <- matrix(rbeta(M*length(N), c(ppw)+1, c(N-ppw)+1), M, length(N), byrow=TRUE)
# sel <- inside(bb, IRT_W1$A, IRT_W1$b)
# t(apply(bb[sel,], 2, mcmc.summ))
IRT_W2 <- nirt_to_Ab(3, 3, axioms = c("W1", "W2"))
pp <- sampling_binom(c(ppw), c(N), IRT_W2$A, IRT_W2$b,
M = 10000, burnin = 1000
)
summ <- t(apply(pp, 2, mcmc.summ))
matrix(summ[, 1], 3)
irt_w2 <- c(.18, .47, .66, .35, .55, .70, .57, .66, .74)
round(cbind(p.obs = c(p.obs), irt_w2, summ), 2)
# expect_equal(unname(summ[,1]), irt_w2, tol = .01) ### TODO
})
############# With double cancellation
# Karabatsos Figure 7, p. 415 (W1, W2, Co)
N2 <- matrix(c(84, 82, 86), 3, 3)
ppw2 <- matrix(c(
15, 20, 10,
11, 27, 25,
11, 24, 55
), 3, 3, byrow = TRUE)
dimnames(ppw2) <- dimnames(N2) <- list(
score = paste0("s", 4:6),
item = paste0("i", c(6, 1, 8))
)
p.obs2 <- round(ppw2 / N2, 2)
test_that("Examples in Karabatsos (2001, Figure 7) give identical results", {
skip("functions not yet fully implemented")
IRT_all <- nirt_to_Ab(3, 3, axioms = c("W1", "W2", "DC"))
pp <- sampling_binom(c(ppw2), c(N2),
IRT_all$A, IRT_all$b, # prior = c(.5,.5),
M = 15000, burnin = 1000
)
dc.post <- sapply(IRT_all$exclude, function(ee) inside(pp, ee$A, ee$b))
dc.violated <- apply(dc.post, 1, any)
mean(dc.post)
pp.dc <- pp[!dc.violated, ]
irt_all <- c(.11, .15, .21, .22, .28, .34, .33, .42, .64)
summ <- t(apply(pp.dc, 2, mcmc.summ))
round(cbind(p.obs2 = c(p.obs2), irt_all, summ), 2)
# expect_equal(unname(summ[,1]), irt_all, tol = .01) ## TODO
})
##### karabatsos 2001, 2-PL data, p. 419, Figure 9
# p <- matrix(c( .00, .28, .00, .06, .17, .50,
# .00, .36, .00, .43, .50, .71,
# .08, .46, .38, .69, .77, .62,
# .38, .50, .81, .44, .94, .94,
# .80, .45, 1.0, .80, .95, 1.0), 5, byrow = TRUE)
# N <- matrix(c(18, 14,13,16, 20), 5,6)
# k <- round(p*N)
#
# test_that("Examples in Karabatsos (2001, Figure 9) give identical results",{
# irt <- nirt_to_Ab(nrow(k), ncol(k), axioms = c("W1", "W2"))
# pp <- sampling_binom(c(k), c(N), irt$A, irt$b, prior=c(1,1),
# M = 10000, burnin = 1000)
# summ <- t(apply(pp, 2, mcmc.summ))
# summ
# matrix(summ[,1], 5)
#
# means <- c(.13, .64, .70, .21, .36, .85, .29, .73, .39, .80)
# subs <- c(2, 5, 10, 11, 12, 15, 16, 19, 21, 28)
# round(cbind(pp = c(p)[subs], means, summ[subs,]), 2)
# expect_equal(unname(summ[subs,1]), means, tol = .01) # TODO
# })
# karabatsos2004, table 3: monotonicity estimates
pmean <- matrix(c(
.05, .09, .17, .07, .45, .44,
.19, .13, .25, .28, .57, .59,
.26, .23, .41, .46, .67, .73,
.36, .43, .50, .66, .72, .81,
.44, .53, .58, .86, .75, .86,
.55, .63, .81, .96, .80, .96
), 6, byrow = TRUE)
ppp.ij <- matrix(c(
.98, .09, .30, .94, .49, .60,
.60, .39, .56, .62, .70, .63,
.72, .62, .62, .62, .51, .63,
.66, .57, .70, .58, .42, .69,
.68, .69, .54, .64, .33, .50,
.49, .39, .74, .98, .08, .98
), 6, byrow = TRUE)
test_that("Karabatsos, 2004 (Table 3) matches (posterior mean)", {
expect_silent(data(karabatsos2004))
IJ <- dim(karabatsos2004$k.M)
monotonicity <- nirt_to_Ab(IJ[1], IJ[2], axioms = "W1")
pp <- sampling_binom(
k = c(karabatsos2004$k.M),
n = c(karabatsos2004$n.M),
A = monotonicity$A, b = monotonicity$b,
prior = c(.5, .5), M = 10000
)
post.mean <- matrix(apply(pp, 2, mean), IJ[1],
dimnames = dimnames(karabatsos2004$k.M)
)
expect_equal(pmean, unname(post.mean), tol = .01)
ppp <- ppp_binom(pp, c(karabatsos2004$k.M), c(karabatsos2004$n.M),
by = 1:prod(IJ)
)[, 3]
expect_equal(unname(ppp), c(ppp.ij), tol = .03)
})
# karabatsos 2004: PPP for IIO
ppp.ij <- matrix(c(
1.0, 1.0, .98, .91, .75, .45,
.94, .46, .72, .20, .61, .51,
.12, .25, .62, .62, .61, .60,
.59, .59, .49, .49, .59, .59,
.60, .36, .26, .14, .39, .58,
.61, .68, .44, .29, .10, .90,
.47, .76, .92, .98, 1.0, 1.0
), 7, byrow = TRUE)
ppp.items <- c(.18, .26, .47, .21, .40, .37)
test_that("Karabatsos, 2004 (Table 6) matches (PPP values)", {
expect_silent(data(karabatsos2004))
IJ <- dim(karabatsos2004$k.IIO)
monotonicity <- nirt_to_Ab(IJ[1], IJ[2], axioms = "W2")
pp <- sampling_binom(
k = c(karabatsos2004$k.IIO),
n = c(karabatsos2004$n.IIO),
A = monotonicity$A, b = monotonicity$b,
prior = c(.5, .5), M = 10000
)
ppp <- ppp_binom(pp, c(karabatsos2004$k.IIO), c(karabatsos2004$n.IIO),
by = 1:prod(IJ)
)
round(matrix(ppp[, 3], 7), 2)
expect_equal(unname(ppp[, 3]), c(ppp.ij), tol = .03)
ppp <- ppp_binom(pp, c(karabatsos2004$k.IIO), c(karabatsos2004$n.IIO),
by = rep(1:IJ[2], each = IJ[1])
)
expect_equal(unname(ppp[, 3]), c(ppp.items), tol = .03)
})
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.