Nothing
tests/testthat/test-galamm-lmm.R
In galamm: Generalized Additive Latent and Mixed Models
#' @srrstats {G5.0} Datasets from PLmixed package are used for testing, and
#' results from the functions in this package are precomputed for comparison,
#' in cases where PLmixed and galamm support the same models.
#' @srrstats {G5.4} It has been confirmed that PLmixed returns the same results.
#' PLmixed is not run inside the tests, because it is too slow for that.
#' @noRd
NULL
data("IRTsim", package = "PLmixed")
test_that("LMM with simple factor works", {
IRTsub <- IRTsim[IRTsim$item < 4, ] # Select items 1-3
set.seed(12345)
IRTsub <- IRTsub[sample(nrow(IRTsub), 300), ] # Randomly sample 300 responses
IRTsub <- IRTsub[order(IRTsub$item), ] # Order by item
irt.lam <- matrix(c(1, NA, NA), ncol = 1) # Specify the lambda matrix
mod <- galamm(
formula = y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub,
load.var = "item",
factor = "abil.sid",
lambda = irt.lam
)
expect_equal(nobs(mod), 300L)
expect_equal(
as.character(formula(mod)),
c("~", "y", "0 + as.factor(item) + (0 + abil.sid | school/sid)")
)
pdf(NULL)
expect_invisible(plot(mod))
dev.off()
class(IRTsub) <- c("tbl_df", "tbl", "data.frame")
mod1 <- galamm(
y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub, load.var = "item",
factor = "abil.sid", lambda = irt.lam
)
expect_equal(deviance(mod1), deviance(mod))
class(IRTsub) <- c("data.table", "data.frame")
mod2 <- galamm(
y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub, , load.var = "item",
factor = "abil.sid", lambda = irt.lam
)
expect_equal(deviance(mod2), deviance(mod))
IRTsub <- as.data.frame(IRTsub)
expect_equal(mod$hessian, mod1$hessian)
expect_equal(mod$hessian, mod2$hessian)
expect_equal(
sd(ranef(mod)$school$abil.sid),
0.180652430814172
)
expect_output(
lme4::.prt.call(mod$call),
"Formula: y ~ 0 + as.factor(item) + (0 + abil.sid | school/sid)",
fixed = TRUE
)
expect_equal(mod$model$loglik, -193.563337783604)
expect_equal(
summary(mod)$AICtab,
c(
AIC = 403.126675567209, BIC = 432.756935364458,
logLik = -193.563337783604,
deviance = 387.126675567209, df.resid = 292
)
)
expect_equal(
factor_loadings(mod),
structure(c(
1, 1.05448712819873, 1.02126746190855, NA, 0.217881890204074,
0.236814881042725
), dim = 3:2, dimnames = list(
c("lambda1", "lambda2", "lambda3"),
c("abil.sid", "SE")
))
)
expect_equal(
residuals(mod)[c(4, 8, 11)],
c(0.0513522294535425, -0.181269847807669, 0.0759916652950277)
)
expect_snapshot(print(VarCorr(mod), digits = 2))
expect_snapshot(round(confint(mod, parm = "beta"), 2))
expect_snapshot(round(confint(mod, parm = "lambda"), 2))
expect_snapshot(round(confint(mod, parm = "theta"), 2))
mod2 <- galamm(
formula = y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub,
load.var = "item",
factor = "abil.sid",
lambda = irt.lam,
start = list(
theta = mod$parameters$parameter_estimates[mod$parameters$theta_inds],
beta = mod$parameters$parameter_estimates[mod$parameters$beta_inds],
lambda = mod$parameters$parameter_estimates[mod$parameters$lambda_inds]
),
control = galamm_control(reduced_hessian = TRUE)
)
expect_snapshot(print(summary(mod2), digits = 2))
})
test_that("LMM with simple factor works with Nelder-Mead", {
IRTsub <- IRTsim[IRTsim$item < 4, ] # Select items 1-3
set.seed(12345)
IRTsub <- IRTsub[sample(nrow(IRTsub), 300), ] # Randomly sample 300 responses
IRTsub <- IRTsub[order(IRTsub$item), ] # Order by item
irt.lam <- matrix(c(1, NA, NA), ncol = 1) # Specify the lambda matrix
mod <- galamm(
formula = y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub,
load.var = "item",
factor = "abil.sid",
lambda = irt.lam,
control = galamm_control(method = "Nelder-Mead")
)
expect_equal(
sd(ranef(mod)$school$abil.sid),
0.180650156300753
)
expect_output(
lme4::.prt.call(mod$call),
"Formula: y ~ 0 + as.factor(item) + (0 + abil.sid | school/sid)",
fixed = TRUE
)
expect_equal(mod$model$loglik, -193.56333776973)
expect_equal(
summary(mod)$AICtab,
c(
AIC = 403.12667553946, BIC = 432.756935336709, logLik = -193.56333776973,
deviance = 387.12667553946, df.resid = 292
)
)
expect_equal(
factor_loadings(mod),
structure(
c(
1, 1.05449523890597, 1.02128079445268, NA, 0.217885861133189,
0.236820565746212
),
dim = 3:2,
dimnames = list(
c("lambda1", "lambda2", "lambda3"),
c("abil.sid", "SE")
)
)
)
expect_equal(
residuals(mod)[c(4, 8, 11)],
c(0.051356304666937, -0.181273547504248, 0.0759965034510736)
)
expect_snapshot(print(VarCorr(mod), digits = 2))
expect_snapshot(round(confint(mod, parm = "beta"), 2))
expect_snapshot(round(confint(mod, parm = "lambda"), 2))
expect_snapshot(round(confint(mod, parm = "theta"), 2))
mod2 <- galamm(
formula = y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub,
load.var = "item",
factor = "abil.sid",
lambda = irt.lam,
start = list(
theta = mod$parameters$parameter_estimates[mod$parameters$theta_inds],
beta = mod$parameters$parameter_estimates[mod$parameters$beta_inds],
lambda = mod$parameters$parameter_estimates[mod$parameters$lambda_inds]
),
control = galamm_control(reduced_hessian = TRUE)
)
expect_snapshot(print(summary(mod2), digits = 2))
})
data("KYPSsim", package = "PLmixed")
test_that("LMM with two factors works", {
# Making data small for it to run faster
dat <- subset(KYPSsim, hid <= 5 & mid <= 5)
dat$time <- factor(dat$time)
kyps.lam <- rbind(
c(1, 0),
c(NA, 0),
c(NA, 1),
c(NA, NA)
)
form <- esteem ~ time + (0 + hs | hid) + (0 + ms | mid) + (1 | sid)
kyps_model <- galamm(
formula = form,
data = dat,
factor = c("ms", "hs"),
load.var = "time",
lambda = kyps.lam
)
expect_equal(kyps_model$model$loglik, -33.7792632108483)
expect_equal(
kyps_model$parameters$parameter_estimates[
kyps_model$parameters$lambda_inds
],
c(
0.631763067142624, -0.114603879076418, 0.0145214343084242,
1.07489949334559
),
tolerance = 1e-4
)
expect_equal(
llikAIC(kyps_model),
c(
AIC = 91.5585264216966, BIC = 115.202029384322,
logLik = -33.7792632108483,
deviance = 67.5585264216966, df.resid = 41
)
)
expect_equal(
residuals(kyps_model)[c(3, 22)],
c(-0.411647523673306, 1.30437183262668),
tolerance = 1e-4
)
expect_equal(
residuals(kyps_model, type = "deviance")[c(3, 22)],
c(-0.411647523673306, 1.30437183262668),
tolerance = 1e-4
)
expect_warning(vcov(kyps_model), "Rank deficient Hessian matrix")
})
data("JUDGEsim", package = "PLmixed")
test_that("LMM with two raters works", {
dat <- subset(JUDGEsim, item %in% 1:6 & stu < 20)
dat$item <- factor(dat$item)
judge.lam <- rbind(
c(1, 0),
c(NA, 0),
c(NA, 0),
c(0, 1),
c(0, NA),
c(0, NA)
)
judge_galamm <- galamm(
formula = response ~ 0 + item + (1 | class) +
(0 + teacher1 + teacher2 | tch),
data = dat,
lambda = judge.lam,
load.var = "item",
factor = c("teacher1", "teacher2")
)
expect_equal(deviance(judge_galamm), 2158.16155409666)
expect_equal(
logLik(judge_galamm),
structure(-1079.08077704833, nobs = 840L, df = 15L, class = "logLik")
)
expect_equal(
llikAIC(judge_galamm),
c(
AIC = 2188.16155409666, BIC = 2259.16258247423,
logLik = -1079.08077704833,
deviance = 2158.16155409666, df.resid = 825
)
)
expect_equal(factor_loadings(judge_galamm),
structure(c(
1, 1.24788910616366, 0.947028051478771, 0, 0, 0,
NA, 0.264534156451991, 0.228633827429379, NA, NA, NA, 0, 0, 0,
1, 1.00426975887412, 1.10790222766263, NA, NA, NA, NA, 0.262801414362859,
0.279745804150548
), dim = c(6L, 4L), dimnames = list(c(
"lambda1",
"lambda2", "lambda3", "lambda4", "lambda5", "lambda6"
), c(
"teacher1",
"SE", "teacher2", "SE"
))),
tolerance = 1e-4
)
expect_equal(quantile(residuals(judge_galamm)),
c(
`0%` = -2.26827478113357, `25%` = -0.578443286096374,
`50%` = -0.00403503141561323,
`75%` = 0.55979417430233, `100%` = 2.22456795070249
),
tolerance = 1e-4
)
})
test_that("Complex LMM works", {
skip_extended()
data(JUDGEsim, package = "PLmixed")
JUDGEsim$item <- factor(JUDGEsim$item)
judge.lam <- rbind(
c(1, 0, 1, 0, 0, 0),
c(NA, 0, NA, 0, 0, 0),
c(NA, 0, NA, 0, 0, 0),
c(0, 1, 0, 1, 0, 0),
c(0, NA, 0, NA, 0, 0),
c(0, NA, 0, NA, 0, 0),
c(0, 0, 0, 0, 1, 0),
c(0, 0, 0, 0, NA, 0),
c(0, 0, 0, 0, NA, 0),
c(0, 0, 0, 0, 0, 1),
c(0, 0, 0, 0, 0, NA),
c(0, 0, 0, 0, 0, NA)
)
judge_galamm <- galamm(
formula = response ~ 0 + item + (1 | class) +
(0 + trait1.t + trait2.t + trait1.s + trait2.s | stu) +
(0 + teacher1 + teacher2 | tch),
data = JUDGEsim,
lambda = judge.lam,
load.var = "item",
factor = c(
"teacher1", "teacher2", "trait1.t",
"trait2.t", "trait1.s", "trait2.s"
)
)
expect_equal(judge_galamm$model$loglik, -56553.2785661794)
expect_equal(judge_galamm$parameters$parameter_estimates,
c(
0.784430334881896, 0.566133764367859, 0.398568132799786,
0.200370294453492,
0.334085183988745, 0.00647882691731993, 0.235677264194959,
0.773891841821503,
0.337997663719164, 0.869131637864806, 0.498543199885178,
0.239475682251201,
0.482285028560489, 0, 3.39914518474943, 3.36263585640721,
3.36210081563546,
2.81984223040228, 2.93869388713431, 2.8771346729679,
3.4260794096106,
3.55223187960471, 3.59726951301179, 2.33415461561759,
2.90902811819496,
2.47043867252208, 1.12783179637647, 0.998580402300312,
0.972482963140881,
1.2190594859786, 1.09151685852923, 1.06570210817052,
1.05362302739479,
0.958118368867656, 1.32218004851917, 1.14483478564702,
0.873958594688374,
1.09623955388905
),
tolerance = 1e-4
)
tmp <- summary(judge_galamm)
expect_equal(tmp$Lambda,
structure(c(
1, 1.12783179637647, 0.998580402300312, 0, 0, 0,
0, 0, 0, 0, 0, 0, NA, 0.0358385979432342, 0.0334494645747354,
NA, NA, NA, NA, NA, NA, NA, NA, NA, 0, 0, 0, 1, 0.972482963140881,
1.2190594859786, 0, 0, 0, 0, 0, 0, NA, NA, NA, NA, 0.0303605948898014,
0.0343626808674974, NA, NA, NA, NA, NA, NA, 1, 1.09151685852923,
1.06570210817052, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA, 0.0217214373104262,
0.0214398508062957, NA, NA, NA, NA, NA, NA, NA, NA, NA, 0, 0,
0, 1, 1.05362302739479, 0.958118368867656, 0, 0, 0, 0, 0, 0,
NA, NA, NA, NA, 0.0257563811709945, 0.0246027494628801, NA, NA,
NA, NA, NA, NA, 0, 0, 0, 0, 0, 0, 1, 1.32218004851917, 1.14483478564702,
0, 0, 0, NA, NA, NA, NA, NA, NA, NA, 0.0610912384587967, 0.0563316301466295,
NA, NA, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.873958594688374,
1.09623955388905, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 0.0441749396116879,
0.0487381242824818
), dim = c(12L, 12L), dimnames = list(c(
"lambda1",
"lambda2", "lambda3", "lambda4", "lambda5", "lambda6", "lambda7",
"lambda8", "lambda9", "lambda10", "lambda11", "lambda12"
), c(
"teacher1",
"SE", "teacher2", "SE", "trait1.t", "SE", "trait2.t", "SE", "trait1.s",
"SE", "trait2.s", "SE"
))),
tolerance = 1e-4
)
})
test_that("multiple factors in fixed effects works", {
path <-
system.file("testdata", "test_multiple_factors.rds", package = "galamm")
dat <- as.data.frame(readRDS(path))
lmat <- matrix(c(
1, NA, NA, 0, 0, 0,
0, 0, 0, 1, NA, NA
), ncol = 2)
mod <- galamm(
formula = y ~ 0 + x:domain1:lambda1 + x:domain2:lambda2 +
(0 + 1 | id),
data = dat,
load.var = "item",
lambda = lmat,
factor = c("lambda1", "lambda2"),
start = list(
theta = 0.744468091602185,
beta = c(1.03995169865897, 1.87422267819485),
lambda = c(
0.478791387562245, 1.94779433858618,
0.466484983394861, 2.02985361769537
),
weights = numeric(0)
),
control = galamm_control(
optim_control = list(
FtolAbs = 1000,
FtolRel = 1000, XtolRel = 1000,
warnOnly = TRUE, xt = rep(1000, 7)
),
method = "Nelder-Mead"
)
)
expect_equal(deviance(mod), 7891.36597569292, tolerance = .001)
})
#' @srrstats {G5.9b} Algorithms are determinstic, so changing random seeds does
#' not affect the results. This is tested here.
#' @srrstats {G5.9} Noise susceptibility tests here.
#' @srrstats {G5.9a} Adding trivial noise and testing here.
#' @noRd
NULL
test_that(
"Algorithm is robust to trivial noise",
{
IRTsub <- IRTsim[IRTsim$item < 4, ] # Select items 1-3
set.seed(12345)
IRTsub <- IRTsub[sample(nrow(IRTsub), 300), ] # Randomly sample 300 responses
IRTsub <- IRTsub[order(IRTsub$item), ] # Order by item
irt.lam <- matrix(c(1, NA, NA), ncol = 1) # Specify the lambda matrix
set.seed(1)
mod <- galamm(
formula = y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub,
load.var = "item",
factor = "abil.sid",
lambda = irt.lam
)
set.seed(2)
mod0 <- galamm(
formula = y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub,
load.var = "item",
factor = "abil.sid",
lambda = irt.lam
)
expect_equal(deviance(mod), deviance(mod0))
expect_equal(coef(mod), coef(mod0))
expect_equal(vcov(mod), vcov(mod0))
expect_equal(factor_loadings(mod), factor_loadings(mod0))
# Add trivial noise
dat <- IRTsub
dat$y <- dat$y + rnorm(nrow(dat), sd = .Machine$double.eps)
mod0 <- galamm(
formula = y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub,
load.var = "item",
factor = "abil.sid",
lambda = irt.lam
)
expect_equal(deviance(mod), deviance(mod0), tolerance = 1e-8)
expect_equal(coef(mod), coef(mod0), tolerance = 1e-8)
expect_equal(vcov(mod), vcov(mod0), tolerance = 1e-8)
expect_equal(factor_loadings(mod), factor_loadings(mod0), tolerance = 1e-8)
}
)
Try the galamm package in your browser
Any scripts or data that you put into this service are public.
galamm documentation built on June 8, 2025, 12:42 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.
#' @srrstats {G5.0} Datasets from PLmixed package are used for testing, and
#' results from the functions in this package are precomputed for comparison,
#' in cases where PLmixed and galamm support the same models.
#' @srrstats {G5.4} It has been confirmed that PLmixed returns the same results.
#' PLmixed is not run inside the tests, because it is too slow for that.
#' @noRd
NULL
data("IRTsim", package = "PLmixed")
test_that("LMM with simple factor works", {
IRTsub <- IRTsim[IRTsim$item < 4, ] # Select items 1-3
set.seed(12345)
IRTsub <- IRTsub[sample(nrow(IRTsub), 300), ] # Randomly sample 300 responses
IRTsub <- IRTsub[order(IRTsub$item), ] # Order by item
irt.lam <- matrix(c(1, NA, NA), ncol = 1) # Specify the lambda matrix
mod <- galamm(
formula = y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub,
load.var = "item",
factor = "abil.sid",
lambda = irt.lam
)
expect_equal(nobs(mod), 300L)
expect_equal(
as.character(formula(mod)),
c("~", "y", "0 + as.factor(item) + (0 + abil.sid | school/sid)")
)
pdf(NULL)
expect_invisible(plot(mod))
dev.off()
class(IRTsub) <- c("tbl_df", "tbl", "data.frame")
mod1 <- galamm(
y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub, load.var = "item",
factor = "abil.sid", lambda = irt.lam
)
expect_equal(deviance(mod1), deviance(mod))
class(IRTsub) <- c("data.table", "data.frame")
mod2 <- galamm(
y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub, , load.var = "item",
factor = "abil.sid", lambda = irt.lam
)
expect_equal(deviance(mod2), deviance(mod))
IRTsub <- as.data.frame(IRTsub)
expect_equal(mod$hessian, mod1$hessian)
expect_equal(mod$hessian, mod2$hessian)
expect_equal(
sd(ranef(mod)$school$abil.sid),
0.180652430814172
)
expect_output(
lme4::.prt.call(mod$call),
"Formula: y ~ 0 + as.factor(item) + (0 + abil.sid | school/sid)",
fixed = TRUE
)
expect_equal(mod$model$loglik, -193.563337783604)
expect_equal(
summary(mod)$AICtab,
c(
AIC = 403.126675567209, BIC = 432.756935364458,
logLik = -193.563337783604,
deviance = 387.126675567209, df.resid = 292
)
)
expect_equal(
factor_loadings(mod),
structure(c(
1, 1.05448712819873, 1.02126746190855, NA, 0.217881890204074,
0.236814881042725
), dim = 3:2, dimnames = list(
c("lambda1", "lambda2", "lambda3"),
c("abil.sid", "SE")
))
)
expect_equal(
residuals(mod)[c(4, 8, 11)],
c(0.0513522294535425, -0.181269847807669, 0.0759916652950277)
)
expect_snapshot(print(VarCorr(mod), digits = 2))
expect_snapshot(round(confint(mod, parm = "beta"), 2))
expect_snapshot(round(confint(mod, parm = "lambda"), 2))
expect_snapshot(round(confint(mod, parm = "theta"), 2))
mod2 <- galamm(
formula = y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub,
load.var = "item",
factor = "abil.sid",
lambda = irt.lam,
start = list(
theta = mod$parameters$parameter_estimates[mod$parameters$theta_inds],
beta = mod$parameters$parameter_estimates[mod$parameters$beta_inds],
lambda = mod$parameters$parameter_estimates[mod$parameters$lambda_inds]
),
control = galamm_control(reduced_hessian = TRUE)
)
expect_snapshot(print(summary(mod2), digits = 2))
})
test_that("LMM with simple factor works with Nelder-Mead", {
IRTsub <- IRTsim[IRTsim$item < 4, ] # Select items 1-3
set.seed(12345)
IRTsub <- IRTsub[sample(nrow(IRTsub), 300), ] # Randomly sample 300 responses
IRTsub <- IRTsub[order(IRTsub$item), ] # Order by item
irt.lam <- matrix(c(1, NA, NA), ncol = 1) # Specify the lambda matrix
mod <- galamm(
formula = y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub,
load.var = "item",
factor = "abil.sid",
lambda = irt.lam,
control = galamm_control(method = "Nelder-Mead")
)
expect_equal(
sd(ranef(mod)$school$abil.sid),
0.180650156300753
)
expect_output(
lme4::.prt.call(mod$call),
"Formula: y ~ 0 + as.factor(item) + (0 + abil.sid | school/sid)",
fixed = TRUE
)
expect_equal(mod$model$loglik, -193.56333776973)
expect_equal(
summary(mod)$AICtab,
c(
AIC = 403.12667553946, BIC = 432.756935336709, logLik = -193.56333776973,
deviance = 387.12667553946, df.resid = 292
)
)
expect_equal(
factor_loadings(mod),
structure(
c(
1, 1.05449523890597, 1.02128079445268, NA, 0.217885861133189,
0.236820565746212
),
dim = 3:2,
dimnames = list(
c("lambda1", "lambda2", "lambda3"),
c("abil.sid", "SE")
)
)
)
expect_equal(
residuals(mod)[c(4, 8, 11)],
c(0.051356304666937, -0.181273547504248, 0.0759965034510736)
)
expect_snapshot(print(VarCorr(mod), digits = 2))
expect_snapshot(round(confint(mod, parm = "beta"), 2))
expect_snapshot(round(confint(mod, parm = "lambda"), 2))
expect_snapshot(round(confint(mod, parm = "theta"), 2))
mod2 <- galamm(
formula = y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub,
load.var = "item",
factor = "abil.sid",
lambda = irt.lam,
start = list(
theta = mod$parameters$parameter_estimates[mod$parameters$theta_inds],
beta = mod$parameters$parameter_estimates[mod$parameters$beta_inds],
lambda = mod$parameters$parameter_estimates[mod$parameters$lambda_inds]
),
control = galamm_control(reduced_hessian = TRUE)
)
expect_snapshot(print(summary(mod2), digits = 2))
})
data("KYPSsim", package = "PLmixed")
test_that("LMM with two factors works", {
# Making data small for it to run faster
dat <- subset(KYPSsim, hid <= 5 & mid <= 5)
dat$time <- factor(dat$time)
kyps.lam <- rbind(
c(1, 0),
c(NA, 0),
c(NA, 1),
c(NA, NA)
)
form <- esteem ~ time + (0 + hs | hid) + (0 + ms | mid) + (1 | sid)
kyps_model <- galamm(
formula = form,
data = dat,
factor = c("ms", "hs"),
load.var = "time",
lambda = kyps.lam
)
expect_equal(kyps_model$model$loglik, -33.7792632108483)
expect_equal(
kyps_model$parameters$parameter_estimates[
kyps_model$parameters$lambda_inds
],
c(
0.631763067142624, -0.114603879076418, 0.0145214343084242,
1.07489949334559
),
tolerance = 1e-4
)
expect_equal(
llikAIC(kyps_model),
c(
AIC = 91.5585264216966, BIC = 115.202029384322,
logLik = -33.7792632108483,
deviance = 67.5585264216966, df.resid = 41
)
)
expect_equal(
residuals(kyps_model)[c(3, 22)],
c(-0.411647523673306, 1.30437183262668),
tolerance = 1e-4
)
expect_equal(
residuals(kyps_model, type = "deviance")[c(3, 22)],
c(-0.411647523673306, 1.30437183262668),
tolerance = 1e-4
)
expect_warning(vcov(kyps_model), "Rank deficient Hessian matrix")
})
data("JUDGEsim", package = "PLmixed")
test_that("LMM with two raters works", {
dat <- subset(JUDGEsim, item %in% 1:6 & stu < 20)
dat$item <- factor(dat$item)
judge.lam <- rbind(
c(1, 0),
c(NA, 0),
c(NA, 0),
c(0, 1),
c(0, NA),
c(0, NA)
)
judge_galamm <- galamm(
formula = response ~ 0 + item + (1 | class) +
(0 + teacher1 + teacher2 | tch),
data = dat,
lambda = judge.lam,
load.var = "item",
factor = c("teacher1", "teacher2")
)
expect_equal(deviance(judge_galamm), 2158.16155409666)
expect_equal(
logLik(judge_galamm),
structure(-1079.08077704833, nobs = 840L, df = 15L, class = "logLik")
)
expect_equal(
llikAIC(judge_galamm),
c(
AIC = 2188.16155409666, BIC = 2259.16258247423,
logLik = -1079.08077704833,
deviance = 2158.16155409666, df.resid = 825
)
)
expect_equal(factor_loadings(judge_galamm),
structure(c(
1, 1.24788910616366, 0.947028051478771, 0, 0, 0,
NA, 0.264534156451991, 0.228633827429379, NA, NA, NA, 0, 0, 0,
1, 1.00426975887412, 1.10790222766263, NA, NA, NA, NA, 0.262801414362859,
0.279745804150548
), dim = c(6L, 4L), dimnames = list(c(
"lambda1",
"lambda2", "lambda3", "lambda4", "lambda5", "lambda6"
), c(
"teacher1",
"SE", "teacher2", "SE"
))),
tolerance = 1e-4
)
expect_equal(quantile(residuals(judge_galamm)),
c(
`0%` = -2.26827478113357, `25%` = -0.578443286096374,
`50%` = -0.00403503141561323,
`75%` = 0.55979417430233, `100%` = 2.22456795070249
),
tolerance = 1e-4
)
})
test_that("Complex LMM works", {
skip_extended()
data(JUDGEsim, package = "PLmixed")
JUDGEsim$item <- factor(JUDGEsim$item)
judge.lam <- rbind(
c(1, 0, 1, 0, 0, 0),
c(NA, 0, NA, 0, 0, 0),
c(NA, 0, NA, 0, 0, 0),
c(0, 1, 0, 1, 0, 0),
c(0, NA, 0, NA, 0, 0),
c(0, NA, 0, NA, 0, 0),
c(0, 0, 0, 0, 1, 0),
c(0, 0, 0, 0, NA, 0),
c(0, 0, 0, 0, NA, 0),
c(0, 0, 0, 0, 0, 1),
c(0, 0, 0, 0, 0, NA),
c(0, 0, 0, 0, 0, NA)
)
judge_galamm <- galamm(
formula = response ~ 0 + item + (1 | class) +
(0 + trait1.t + trait2.t + trait1.s + trait2.s | stu) +
(0 + teacher1 + teacher2 | tch),
data = JUDGEsim,
lambda = judge.lam,
load.var = "item",
factor = c(
"teacher1", "teacher2", "trait1.t",
"trait2.t", "trait1.s", "trait2.s"
)
)
expect_equal(judge_galamm$model$loglik, -56553.2785661794)
expect_equal(judge_galamm$parameters$parameter_estimates,
c(
0.784430334881896, 0.566133764367859, 0.398568132799786,
0.200370294453492,
0.334085183988745, 0.00647882691731993, 0.235677264194959,
0.773891841821503,
0.337997663719164, 0.869131637864806, 0.498543199885178,
0.239475682251201,
0.482285028560489, 0, 3.39914518474943, 3.36263585640721,
3.36210081563546,
2.81984223040228, 2.93869388713431, 2.8771346729679,
3.4260794096106,
3.55223187960471, 3.59726951301179, 2.33415461561759,
2.90902811819496,
2.47043867252208, 1.12783179637647, 0.998580402300312,
0.972482963140881,
1.2190594859786, 1.09151685852923, 1.06570210817052,
1.05362302739479,
0.958118368867656, 1.32218004851917, 1.14483478564702,
0.873958594688374,
1.09623955388905
),
tolerance = 1e-4
)
tmp <- summary(judge_galamm)
expect_equal(tmp$Lambda,
structure(c(
1, 1.12783179637647, 0.998580402300312, 0, 0, 0,
0, 0, 0, 0, 0, 0, NA, 0.0358385979432342, 0.0334494645747354,
NA, NA, NA, NA, NA, NA, NA, NA, NA, 0, 0, 0, 1, 0.972482963140881,
1.2190594859786, 0, 0, 0, 0, 0, 0, NA, NA, NA, NA, 0.0303605948898014,
0.0343626808674974, NA, NA, NA, NA, NA, NA, 1, 1.09151685852923,
1.06570210817052, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA, 0.0217214373104262,
0.0214398508062957, NA, NA, NA, NA, NA, NA, NA, NA, NA, 0, 0,
0, 1, 1.05362302739479, 0.958118368867656, 0, 0, 0, 0, 0, 0,
NA, NA, NA, NA, 0.0257563811709945, 0.0246027494628801, NA, NA,
NA, NA, NA, NA, 0, 0, 0, 0, 0, 0, 1, 1.32218004851917, 1.14483478564702,
0, 0, 0, NA, NA, NA, NA, NA, NA, NA, 0.0610912384587967, 0.0563316301466295,
NA, NA, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.873958594688374,
1.09623955388905, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 0.0441749396116879,
0.0487381242824818
), dim = c(12L, 12L), dimnames = list(c(
"lambda1",
"lambda2", "lambda3", "lambda4", "lambda5", "lambda6", "lambda7",
"lambda8", "lambda9", "lambda10", "lambda11", "lambda12"
), c(
"teacher1",
"SE", "teacher2", "SE", "trait1.t", "SE", "trait2.t", "SE", "trait1.s",
"SE", "trait2.s", "SE"
))),
tolerance = 1e-4
)
})
test_that("multiple factors in fixed effects works", {
path <-
system.file("testdata", "test_multiple_factors.rds", package = "galamm")
dat <- as.data.frame(readRDS(path))
lmat <- matrix(c(
1, NA, NA, 0, 0, 0,
0, 0, 0, 1, NA, NA
), ncol = 2)
mod <- galamm(
formula = y ~ 0 + x:domain1:lambda1 + x:domain2:lambda2 +
(0 + 1 | id),
data = dat,
load.var = "item",
lambda = lmat,
factor = c("lambda1", "lambda2"),
start = list(
theta = 0.744468091602185,
beta = c(1.03995169865897, 1.87422267819485),
lambda = c(
0.478791387562245, 1.94779433858618,
0.466484983394861, 2.02985361769537
),
weights = numeric(0)
),
control = galamm_control(
optim_control = list(
FtolAbs = 1000,
FtolRel = 1000, XtolRel = 1000,
warnOnly = TRUE, xt = rep(1000, 7)
),
method = "Nelder-Mead"
)
)
expect_equal(deviance(mod), 7891.36597569292, tolerance = .001)
})
#' @srrstats {G5.9b} Algorithms are determinstic, so changing random seeds does
#' not affect the results. This is tested here.
#' @srrstats {G5.9} Noise susceptibility tests here.
#' @srrstats {G5.9a} Adding trivial noise and testing here.
#' @noRd
NULL
test_that(
"Algorithm is robust to trivial noise",
{
IRTsub <- IRTsim[IRTsim$item < 4, ] # Select items 1-3
set.seed(12345)
IRTsub <- IRTsub[sample(nrow(IRTsub), 300), ] # Randomly sample 300 responses
IRTsub <- IRTsub[order(IRTsub$item), ] # Order by item
irt.lam <- matrix(c(1, NA, NA), ncol = 1) # Specify the lambda matrix
set.seed(1)
mod <- galamm(
formula = y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub,
load.var = "item",
factor = "abil.sid",
lambda = irt.lam
)
set.seed(2)
mod0 <- galamm(
formula = y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub,
load.var = "item",
factor = "abil.sid",
lambda = irt.lam
)
expect_equal(deviance(mod), deviance(mod0))
expect_equal(coef(mod), coef(mod0))
expect_equal(vcov(mod), vcov(mod0))
expect_equal(factor_loadings(mod), factor_loadings(mod0))
# Add trivial noise
dat <- IRTsub
dat$y <- dat$y + rnorm(nrow(dat), sd = .Machine$double.eps)
mod0 <- galamm(
formula = y ~ 0 + as.factor(item) + (0 + abil.sid | school / sid),
data = IRTsub,
load.var = "item",
factor = "abil.sid",
lambda = irt.lam
)
expect_equal(deviance(mod), deviance(mod0), tolerance = 1e-8)
expect_equal(coef(mod), coef(mod0), tolerance = 1e-8)
expect_equal(vcov(mod), vcov(mod0), tolerance = 1e-8)
expect_equal(factor_loadings(mod), factor_loadings(mod0), tolerance = 1e-8)
}
)
Try the galamm package in your browser
Any scripts or data that you put into this service are public.
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.