Nothing
context("blmer numerical results with cov prior")
source(system.file("common", "lmmData.R", package = "blme"))
lme4Version <- packageVersion("lme4")
control <- lmerControl(optimizer = "Nelder_Mead")
test_that("blmer fits test data with gamma prior(TRUE), matching previous version", {
cov.prior <- "g.1 ~ gamma(rate = 0.5)"
fit <- blmer(y ~ x.1 + x.2 + (1 | g.1), testData, cov.prior = cov.prior, control = control)
if (lme4Version < "1.1-8") {
expect_equal(fit@theta, 0.626025390625)
} else {
expect_equal(fit@theta, 0.626021723159)
}
})
test_that("blmer fits test data with invgamma prior(TRUE), matching previous version", {
cov.prior <- "g.1 ~ invgamma(scale = 2.0)"
fit <- blmer(y ~ x.1 + x.2 + (1 | g.1), testData, cov.prior = cov.prior, control = control)
if (lme4Version < "1.1-4") {
expect_equal(fit@theta, 0.93956054688)
} else if (lme4Version < "1.1-8") {
expect_equal(fit@theta, 0.93955078125)
} else {
expect_equal(fit@theta, 0.93955687941)
}
})
test_that("blmer fits test data with wishart prior(TRUE), matching previous version", {
cov.prior <- "g.1 ~ wishart(scale = 2)"
fit <- blmer(y ~ x.1 + x.2 + (1 + x.1 | g.1), testData, cov.prior = cov.prior, control = control)
expect_equal(fit@theta, c(0.677745102365688, -0.439777135132983, 1.48026251108622))
})
test_that("blmer fits test data with invwishart prior(TRUE), matching previous version", {
cov.prior <- "g.1 ~ invwishart(scale = 2)"
fit <- blmer(y ~ x.1 + x.2 + (1 + x.1 | g.1), testData, cov.prior = cov.prior, control = control)
expect_equal(fit@theta, c(0.627739008945695, -0.137563742254117, 1.05679359569432))
})
test_that("blmer fits test data with gamma prior('var', FALSE), matching previous version", {
cov.prior <- "g.1 ~ gamma(shape = 1.75, rate = 2, posterior.scale = 'var', common.scale = FALSE)"
fit <- blmer(y ~ x.1 + x.2 + (1 | g.1), testData, cov.prior = cov.prior, control = control)
if (lme4Version < "1.1-8") {
expect_equal(fit@theta, 0.435458984375)
} else {
expect_equal(fit@theta, 0.435465082534)
}
})
test_that("blmer fits test data with invgamma prior('var', FALSE), matching previous version", {
cov.prior <- "g.1 ~ invgamma(scale = 0.5, posterior.scale = 'var', common.scale = FALSE)"
fit <- blmer(y ~ x.1 + x.2 + (1 | g.1), testData, cov.prior = cov.prior, control = control)
if (lme4Version < "1.1-4") {
expect_equal(fit@theta, 0.460400390625)
} else if (lme4Version < "1.1-8") {
expect_equal(fit@theta, 0.460410156250)
} else {
expect_equal(fit@theta, 0.460406488784)
}
})
test_that("blmer fits test data with gamma prior('sd', FALSE), matching previous version", {
cov.prior <- "g.1 ~ gamma(rate = 2, posterior.scale = 'sd', common.scale = FALSE)"
fit <- blmer(y ~ x.1 + x.2 + (1 | g.1), testData, cov.prior = cov.prior, control = control)
if (lme4Version < "1.1-8") {
expect_equal(fit@theta, 0.476779702210)
} else {
expect_equal(fit@theta, 0.476774614593)
}
})
test_that("blmer fits test data with invgamma prior('sd', FALSE), matching previous version", {
cov.prior <- "g.1 ~ invgamma(scale = 0.5, posterior.scale = 'sd', common.scale = FALSE)"
fit <- blmer(y ~ x.1 + x.2 + (1 | g.1), testData, cov.prior = cov.prior, control = control)
if (lme4Version < "1.1-4") {
expect_equal(fit@theta, 0.452841796875)
} else if (lme4Version < "1.1-8") {
expect_equal(fit@theta, 0.452832031250)
} else {
expect_equal(fit@theta, 0.452838129409)
}
})
test_that("blmer fits test data with wishart prior(FALSE), matching previous version", {
cov.prior <- "g.1 ~ wishart(scale = 2, common.scale = FALSE)"
fit <- blmer(y ~ x.1 + x.2 + (1 + x.1 | g.1), testData, cov.prior = cov.prior, control = control)
expect_equal(fit@pp$theta, c(0.63996739265564, -0.340538787006457, 1.34228986794088))
})
test_that("blmer fits test data with invwishart prior(FALSE), matching previous version", {
cov.prior <- "g.1 ~ invwishart(scale = 2, common.scale = FALSE)"
fit <- blmer(y ~ x.1 + x.2 + (1 + x.1 | g.1), testData, cov.prior = cov.prior, control = control)
expect_equal(fit@pp$theta, c(0.505864621989816, -0.137623340382083, 0.979903012179649))
})
test_that("blmer fits test data with custom prior, matching builtin wishart", {
dwish <- function(R) {
d <- nrow(R)
nu <- d + 1 + 1.5
R.scale.inv <- diag(1e-2, d)
const <- nu * (d * log(2) - 2 * sum(log(diag(R.scale.inv)))) +
0.5 * d * (d - 1) * log(pi)
for (i in 1:d) const <- const + 2 * lgamma(0.5 * (nu + 1.0 - i))
det <- 2 * sum(log(diag(R)))
const - (nu - d - 1) * det + sum((R %*% R.scale.inv)^2)
}
fit.prof <- blmer(y ~ x.1 + x.2 + (1 + x.1 | g.1), testData, control = control,
cov.prior = wishart(scale = diag(1e4, q.k)))
fit.cust <- blmer(y ~ x.1 + x.2 + (1 + x.1 | g.1), testData, control = control,
cov.prior = custom(dwish, chol = TRUE, scale = "dev"))
expect_equal(fit.prof@theta, fit.cust@theta, tolerance = 1e-6)
fit.prof <- blmer(y ~ x.1 + x.2 + (1 + x.1 | g.1), testData, control = control,
cov.prior = wishart(scale = diag(1e4, q.k), common.scale = FALSE))
fit.cust <- blmer(y ~ x.1 + x.2 + (1 + x.1 | g.1), testData, control = control,
cov.prior = custom(dwish, chol = TRUE, scale = "dev", common.scale = FALSE))
expect_equal(c(fit.prof@pp$theta, fit.prof@devcomp$cmp[["sigmaREML"]]),
c(fit.cust@pp$theta, fit.cust@devcomp$cmp[["sigmaREML"]]), tolerance = 5e-5)
})
test_that("blmer not confused by cov priors supplied in different forms", {
control$optCtrl <- list(maxfun = 1L)
control$checkConv <- NULL
cov.prior <- list("g.1 ~ invgamma(scale = 2)", "g.2 ~ invgamma(scale = 3)")
expect_is(suppressWarnings(blmer(y ~ x.1 + x.2 + (1 | g.1) + (1 | g.2), testData, cov.prior = cov.prior, control = control)),
"blmerMod")
cov.prior <- list(g.1 ~ invgamma(scale = 2), g.2 ~ invgamma(scale = 3))
expect_is(suppressWarnings(blmer(y ~ x.1 + x.2 + (1 | g.1) + (1 | g.2), testData, cov.prior = cov.prior, control = control)),
"blmerMod")
expect_is(suppressWarnings(blmer(y ~ x.1 + x.2 + (1 | g.1), testData, cov.prior = wishart, control = control)), "blmerMod")
expect_is(suppressWarnings(blmer(y ~ x.1 + x.2 + (1 | g.1), testData, cov.prior = wishart(), control = control)), "blmerMod")
})
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.