Nothing
context("Tests for make_standata")
test_that(paste("make_standata returns correct data names ",
"for fixed and random effects"), {
expect_equal(sort(names(make_standata(rating ~ treat + period + carry
+ (1|subject), data = inhaler))),
sort(c("N", "Y", "K", "Kc", "X", "Z_1_1",
"J_1", "N_1", "M_1", "NC_1", "prior_only")))
expect_equal(sort(names(make_standata(rating ~ treat + period + carry
+ (1+treat|id|subject), data = inhaler,
family = "categorical"))),
sort(c("N", "Y", "ncat", "K_mu2", "Kc_mu2", "X_mu2", "Z_1_mu2_1",
"Z_1_mu2_2", "K_mu3", "Kc_mu3", "X_mu3", "Z_1_mu3_3", "Z_1_mu3_4",
"K_mu4", "Kc_mu4", "X_mu4", "Z_1_mu4_5", "Z_1_mu4_6",
"J_1", "N_1", "M_1", "NC_1", "prior_only")))
expect_equal(sort(names(make_standata(rating ~ treat + period + carry
+ (1+treat|subject), data = inhaler))),
sort(c("N", "Y", "K", "Kc", "X", "Z_1_1", "Z_1_2", "J_1", "N_1", "M_1",
"NC_1", "prior_only")))
dat <- data.frame(y = 1:10, g = 1:10, h = 11:10, x = rep(0,10))
expect_equal(sort(names(make_standata(y ~ 0 + Intercept + x + (1|g) + (1|h),
dat, "poisson"))),
sort(c("N", "Y", "K", "X", "Z_1_1", "Z_2_1",
"J_1", "J_2", "N_1", "M_1", "NC_1", "N_2", "M_2", "NC_2",
"prior_only")))
expect_true(all(c("Z_1_1", "Z_1_2", "Z_2_1", "Z_2_2") %in%
names(make_standata(y ~ x + (1+x|g/h), dat))))
expect_equal(make_standata(y ~ x + (1+x|g+h), dat),
make_standata(y ~ x + (1+x|g) + (1+x|h), dat))
})
test_that(paste("make_standata handles variables used as fixed effects",
"and grouping factors at the same time"), {
data <- data.frame(y = 1:9, x = factor(rep(c("a","b","c"), 3)))
standata <- make_standata(y ~ x + (1|x), data = data)
expect_equal(colnames(standata$X), c("Intercept", "xb", "xc"))
expect_equal(standata$J_1, as.array(rep(1:3, 3)))
standata2 <- make_standata(y ~ x + (1|x), data = data,
control = list(not4stan = TRUE))
expect_equal(colnames(standata2$X), c("Intercept", "xb", "xc"))
})
test_that("make_standata returns correct data names for addition terms", {
dat <- data.frame(y = 1:10, w = 1:10, t = 1:10, x = rep(0,10),
c = sample(-1:1,10,TRUE))
expect_equal(names(make_standata(y | se(w) ~ x, dat, gaussian())),
c("N", "Y", "se", "K", "Kc", "X", "sigma", "prior_only"))
expect_equal(names(make_standata(y | weights(w) ~ x, dat, "gaussian")),
c("N", "Y", "weights", "K", "Kc", "X", "prior_only"))
expect_equal(names(make_standata(y | cens(c) ~ x, dat, "student")),
c("N", "Y", "cens", "K", "Kc", "X", "prior_only"))
expect_equal(names(make_standata(y | trials(t) ~ x, dat, "binomial")),
c("N", "Y", "trials", "K", "Kc", "X", "prior_only"))
expect_equal(names(make_standata(y | trials(10) ~ x, dat, "binomial")),
c("N", "Y", "trials", "K", "Kc", "X", "prior_only"))
expect_equal(names(make_standata(y | thres(11) ~ x, dat, "acat")),
c("N", "Y", "nthres", "K", "Kc", "X", "disc", "prior_only"))
expect_equal(names(make_standata(y | thres(10) ~ x, dat, cumulative())),
c("N", "Y", "nthres", "K", "Kc", "X", "disc", "prior_only"))
sdata <- make_standata(y | trunc(0,20) ~ x, dat, "gaussian")
expect_true(all(sdata$lb == 0) && all(sdata$ub == 20))
sdata <- make_standata(y | trunc(ub = 21:30) ~ x, dat)
expect_true(all(all(sdata$ub == 21:30)))
})
test_that(paste("make_standata accepts correct response variables",
"depending on the family"), {
expect_equal(make_standata(y ~ 1, data = data.frame(y = seq(-9.9,0,0.1)),
family = "student")$Y, as.array(seq(-9.9,0,0.1)))
expect_equal(make_standata(y | trials(10) ~ 1, data = data.frame(y = 1:10),
family = "binomial")$Y, as.array(1:10))
expect_equal(make_standata(y ~ 1, data = data.frame(y = 10:20),
family = "poisson")$Y, as.array(10:20))
expect_equal(make_standata(y ~ 1, data = data.frame(y = rep(-c(1:2),5)),
family = "bernoulli")$Y, as.array(rep(1:0,5)))
expect_equal(make_standata(y ~ 1, data = data.frame(y = rep(c(TRUE, FALSE),5)),
family = "bernoulli")$Y, as.array(rep(1:0,5)))
expect_equal(make_standata(y ~ 1, data = data.frame(y = rep(1,5)),
family = "bernoulli")$Y, as.array(rep(1, 5)))
expect_equal(make_standata(y ~ 1, data = data.frame(y = rep(0,5)),
family = "bernoulli")$Y, as.array(rep(0, 5)))
expect_equal(make_standata(y ~ 1, data = data.frame(y = rep(1:10,5)),
family = "categorical")$Y, as.array(rep(1:10,5)))
expect_equal(make_standata(y ~ 1, data = data.frame(y = rep(11:20,5)),
family = "categorical")$Y, as.array(rep(1:10,5)))
expect_equal(make_standata(y ~ 1, data = data.frame(y = factor(rep(11:20,5))),
family = "categorical")$Y, as.array(rep(1:10,5)))
expect_equal(make_standata(y ~ 1, data = data.frame(y = rep(1:10,5)),
family = "cumulative")$Y, as.array(rep(1:10,5)))
dat <- data.frame(y = factor(rep(-4:5,5), order = TRUE))
expect_equal(make_standata(y ~ 1, data = dat, family = "acat")$Y,
as.array(rep(1:10,5)))
expect_equal(make_standata(y ~ 1, data = data.frame(y = seq(1,10,0.1)),
family = "exponential")$Y, as.array(seq(1,10,0.1)))
dat <- data.frame(y1 = 1:10, y2 = 11:20, x = rep(0,10))
form <- bf(mvbind(y1, y2) ~ x) + set_rescor(TRUE)
sdata <- make_standata(form, data = dat)
expect_equal(sdata$Y_y1, as.array(1:10))
expect_equal(sdata$Y_y2, as.array(11:20))
})
test_that(paste("make_standata rejects incorrect response variables",
"depending on the family"), {
expect_error(make_standata(y ~ 1, data = data.frame(y = factor(1:10)),
family = "student"),
"Family 'student' requires numeric responses")
expect_error(make_standata(y ~ 1, data = data.frame(y = -5:5),
family = "geometric"),
"Family 'geometric' requires response greater than or equal to 0")
expect_error(make_standata(y ~ 1, data = data.frame(y = -1:1),
family = "bernoulli"),
"contain only two different values")
expect_error(make_standata(y ~ 1, data = data.frame(y = factor(-1:1)),
family = "cratio"),
"Family 'cratio' requires either positive integers or ordered factors")
expect_error(make_standata(y ~ 1, data = data.frame(y = rep(0.5:7.5), 2),
family = "sratio"),
"Family 'sratio' requires either positive integers or ordered factors")
expect_error(make_standata(y ~ 1, data = data.frame(y = rep(-7.5:7.5), 2),
family = "gamma"),
"Family 'gamma' requires response greater than 0")
expect_error(make_standata(y ~ 1, data = data.frame(y = c(0.1, 0.5, 1)),
family = Beta()),
"Family 'beta' requires response smaller than 1")
expect_error(make_standata(y ~ 1, data = data.frame(y = c(0, 0.5, 4)),
family = von_mises()),
"Family 'von_mises' requires response smaller than or equal to 3.14")
expect_error(make_standata(y ~ 1, data = data.frame(y = c(-1, 2, 5)),
family = hurdle_gamma()),
"Family 'hurdle_gamma' requires response greater than or equal to 0")
})
test_that("make_standata suggests using family bernoulli if appropriate", {
expect_message(make_standata(y | trials(1) ~ 1, data = list(y = rep(0:1,5)),
family = "binomial"),
"family 'bernoulli' might be a more efficient choice.")
expect_message(make_standata(y ~ 1, data = data.frame(y = rep(1:2, 5)),
family = "acat"),
"family 'bernoulli' might be a more efficient choice.")
expect_message(make_standata(y ~ 1, data = data.frame(y = rep(0:1,5)),
family = "categorical"),
"family 'bernoulli' might be a more efficient choice.")
})
test_that("make_standata returns correct values for addition terms", {
dat <- data.frame(y = rnorm(9), s = 1:9, w = 1:9, c1 = rep(-1:1, 3),
c2 = rep(c("left","none","right"), 3),
c3 = c(rep(c(TRUE, FALSE), 4), FALSE),
c4 = c(sample(-1:1, 5, TRUE), rep(2, 4)),
t = 11:19)
expect_equivalent(make_standata(y | se(s) ~ 1, data = dat)$se,
as.array(1:9))
expect_equal(make_standata(y | weights(w) ~ 1, data = dat)$weights,
as.array(1:9))
expect_equal(make_standata(y | cens(c1) ~ 1, data = dat)$cens,
as.array(rep(-1:1, 3)))
expect_equal(make_standata(y | cens(c2) ~ 1, data = dat)$cens,
as.array(rep(-1:1, 3)))
expect_equal(make_standata(y | cens(c3) ~ 1, data = dat)$cens,
as.array(c(rep(1:0, 4), 0)))
expect_equal(make_standata(y | cens(c4, y + 2) ~ 1, data = dat)$rcens,
as.array(c(rep(0, 5), dat$y[6:9] + 2)))
expect_equal(make_standata(s | trials(10) ~ 1, dat,
family = "binomial")$trials,
as.array(rep(10, 9)))
expect_equal(make_standata(s | trials(t) ~ 1, data = dat,
family = "binomial")$trials,
as.array(11:19))
expect_equal(SW(make_standata(s | cat(19) ~ 1, data = dat,
family = "cumulative"))$nthres,
18)
})
test_that("make_standata rejects incorrect addition terms", {
dat <- data.frame(y = rnorm(9), s = -(1:9), w = -(1:9),
c = rep(-2:0, 3), t = 9:1, z = 1:9)
expect_error(make_standata(y | se(s) ~ 1, data = dat),
"Standard errors must be non-negative")
expect_error(make_standata(y | weights(w) ~ 1, data = dat),
"Weights must be non-negative")
expect_error(make_standata(y | cens(c) ~ 1, data = dat))
expect_error(make_standata(z | trials(t) ~ 1, data = dat,
family = "binomial"),
"Number of trials is smaller than the number of events")
})
test_that("make_standata handles multivariate models", {
dat <- data.frame(
y1 = 1:10, y2 = 11:20,
x = rep(0, 10), g = rep(1:2, 5),
censi = sample(0:1, 10, TRUE),
tim = 10:1, w = 1:10
)
form <- bf(mvbind(y1, y2) | weights(w) ~ x) + set_rescor(TRUE)
sdata <- make_standata(form, data = dat)
expect_equal(sdata$Y_y1, as.array(dat$y1))
expect_equal(sdata$Y_y2, as.array(dat$y2))
expect_equal(sdata$weights_y1, as.array(1:10))
expect_error(make_standata(bf(mvbind(y1, y2, y2) ~ x) + set_resor(FALSE),
data = dat),
"Cannot use the same response variable twice")
form <- bf(mvbind(y1 / y2, y2, y1 * 3) ~ x) + set_rescor(FALSE)
sdata <- make_standata(form, data = dat)
expect_equal(sdata$Y_y1y2, as.array(dat$y1 / dat$y2))
sdata <- suppressWarnings(
make_standata(mvbind(y1, y2) ~ x, dat, autocor = cor_ar(~ tim | g))
)
target1 <- c(seq(9, 1, -2), seq(10, 2, -2))
expect_equal(sdata$Y_y1, as.array(target1))
target2 <- c(seq(19, 11, -2), seq(20, 12, -2))
expect_equal(sdata$Y_y2, as.array(target2))
# models without residual correlations
expect_warning(
bform <- bf(y1 | cens(censi) ~ x + y2 + (x|2|g)) +
gaussian() + cor_ar() +
(bf(x ~ 1) + mixture(poisson, nmix = 2)) +
(bf(y2 ~ s(y2) + (1|2|g)) + skew_normal()),
"Using 'cor_brms' objects for 'autocor' is deprecated"
)
bprior <- prior(normal(0, 5), resp = y1) +
prior(normal(0, 10), resp = y2) +
prior(dirichlet(2, 1), theta, resp = x)
sdata <- make_standata(bform, dat, prior = bprior)
sdata_names <- c(
"N", "J_1_y1", "cens_y1", "Kma_y1", "Z_1_y2_3",
"Zs_y2_1_1", "Y_y2", "con_theta_x", "X_mu2_x"
)
expect_true(all(sdata_names %in% names(sdata)))
expect_equal(sdata$con_theta_x, as.array(c(2, 1)))
})
test_that("make_standata removes NAs correctly", {
dat <- data.frame(y = c(rnorm(9), NA))
sdata <- suppressWarnings(make_standata(y ~ 1, dat))
expect_equal(as.numeric(sdata$Y), dat$y[1:9])
})
test_that("make_standata handles the 'subset' addition argument correctly", {
dat1 <- data.frame(
y1 = rnorm(15), y2 = NA,
x1 = rnorm(15), x2 = NA, x3 = rnorm(15),
sub1 = 1, sub2 = 0
)
dat2 <- data.frame(
y1 = NA, y2 = rnorm(10),
x1 = NA, x2 = rnorm(10), x3 = NA,
sub1 = 0, sub2 = 1
)
dat <- rbind(dat1, dat2)
bform <-
bf(y1 | subset(sub1) ~ x1*x3 + sin(x1), family = gaussian()) +
bf(y2 | subset(sub2) ~ x2, family = gaussian()) +
set_rescor(FALSE)
sdata <- make_standata(bform, dat)
nsub1 <- sum(dat$sub1)
nsub2 <- sum(dat$sub2)
expect_equal(sdata$N_y1, nsub1)
expect_equal(sdata$N_y2, nsub2)
expect_equal(length(sdata$Y_y1), nsub1)
expect_equal(nrow(sdata$X_y2), nsub2)
})
test_that("make_standata returns correct data for ARMA terms", {
dat <- data.frame(y = 1:10, x = rep(0, 10), tim = 10:1, g = rep(3:4, 5))
sdata <- make_standata(y ~ x + ma(tim, g), data = dat)
expect_equal(sdata$J_lag, as.array(c(1, 1, 1, 1, 0, 1, 1, 1, 1, 0)))
sdata <- make_standata(y ~ x + ar(tim, g, p = 2), data = dat)
expect_equal(sdata$J_lag, as.array(c(1, 2, 2, 2, 0, 1, 2, 2, 2, 0)))
sdata <- make_standata(y ~ x + ar(tim, g, cov = TRUE), data = dat)
expect_equal(sdata$begin_tg, as.array(c(1, 6)))
expect_equal(sdata$nobs_tg, as.array(c(5, 5)))
sdata <- make_standata(y ~ x + ar(tim), data = dat, family = poisson(),
prior = prior(horseshoe(), class = sderr))
expect_equal(sdata$Kscales, 1)
bform <- bf(y ~ exp(b * x), b ~ 1, nl = TRUE, autocor = ~arma())
sdata <- make_standata(bform, dat)
})
test_that("make_standata returns correct data for UNSTR covariance terms", {
dat <- data.frame(y = 1:12, x = rnorm(12), tim = c(5:1, 1:5, c(0, 4)),
g = c(rep(3:4, 5), rep(2, 2)))
sdata <- make_standata(y ~ x + unstr(tim, g), data = dat)
expect_equal(sdata$n_unique_t, 6)
expect_equal(sdata$n_unique_cortime, 15)
Jtime <- rbind(c(1, 5, 0, 0, 0), 2:6, 2:6)
expect_equal(sdata$Jtime_tg, Jtime)
})
test_that("make_standata allows to retrieve the initial data order", {
dat <- data.frame(y1 = rnorm(100), y2 = rnorm(100),
id = sample(1:10, 100, TRUE),
time = sample(1:100, 100))
# univariate model
sdata1 <- make_standata(y1 ~ ar(time, id), data = dat, internal = TRUE)
expect_equal(dat$y1, as.numeric(sdata1$Y[attr(sdata1, "old_order")]))
# multivariate model
form <- bf(mvbind(y1, y2) ~ ma(time, id)) + set_rescor(FALSE)
sdata2 <- make_standata(form, data = dat, internal = TRUE)
expect_equal(sdata2$Y_y1[attr(sdata2, "old_order")], as.array(dat$y1))
expect_equal(sdata2$Y_y2[attr(sdata2, "old_order")], as.array(dat$y2))
})
test_that("make_standata handles covariance matrices correctly", {
A <- structure(diag(1, 4), dimnames = list(1:4, NULL))
sdata <- make_standata(count ~ Trt + (1|gr(visit, cov = A)),
data = epilepsy, data2 = list(A = A))
expect_equivalent(sdata$Lcov_1, t(chol(A)))
B <- structure(diag(1:5), dimnames = list(c(1,5,2,4,3), NULL))
sdata <- make_standata(count ~ Trt + (1|gr(visit, cov = B)),
data = epilepsy, data2 = list(B = B))
expect_equivalent(sdata$Lcov_1, t(chol(B[c(1,3,5,4), c(1,3,5,4)])))
B <- diag(1, 4)
expect_error(make_standata(count ~ Trt + (1|gr(visit, cov = B)),
data = epilepsy, data2 = list(B = B)),
"Row or column names are required")
B <- structure(diag(1, 4), dimnames = list(2:5, NULL))
expect_error(make_standata(count ~ Trt + (1|gr(visit, cov = B)),
data = epilepsy, data2 = list(B = B)),
"Levels of .* do not match")
B <- A
B[1,2] <- 0.5
expect_error(make_standata(count ~ Trt + (1|gr(visit, cov = B)),
data = epilepsy, data2 = list(B = B)),
"must be symmetric")
expect_warning(
sdata <- make_standata(count ~ Trt + (1|visit), data = epilepsy,
cov_ranef = list(visit = A)),
"Argument 'cov_ranef' is deprecated"
)
expect_equivalent(sdata$Lcov_1, t(chol(A)))
})
test_that("make_standata correctly prepares data for non-linear models", {
flist <- list(a ~ x + (1|1|g), b ~ mo(z) + (1|1|g))
dat <- data.frame(
y = rnorm(9), x = rnorm(9), z = sample(1:9, 9), g = rep(1:3, 3)
)
bform <- bf(y ~ a - b^z, flist = flist, nl = TRUE)
sdata <- make_standata(bform, data = dat)
expect_equal(names(sdata),
c("N", "Y", "C_1", "K_a", "X_a", "Z_1_a_1",
"K_b", "X_b", "Ksp_b", "Imo_b", "Xmo_b_1", "Jmo_b",
"con_simo_b_1", "Z_1_b_2", "J_1", "N_1",
"M_1", "NC_1", "prior_only")
)
expect_equal(colnames(sdata$X_a), c("Intercept", "x"))
expect_equal(sdata$J_1, as.array(dat$g))
bform <- bf(y ~ x) +
nlf(sigma ~ a1 * exp(-x/(a2 + z))) +
lf(a1 ~ 1, a2 ~ z + (x|g)) +
lf(alpha ~ x)
sdata <- make_standata(bform, dat, family = skew_normal())
sdata_names <- c("C_sigma_1", "X_a2", "Z_1_a2_1")
expect_true(all(sdata_names %in% names(sdata)))
})
test_that("make_standata correctly prepares data for monotonic effects", {
data <- data.frame(
y = rpois(120, 10), x1 = rep(1:4, 30), z = rnorm(10),
x2 = factor(rep(c("a", "b", "c"), 40), ordered = TRUE)
)
sdata <- make_standata(y ~ mo(x1)*mo(x2)*y, data = data)
sdata_names <- c("Xmo_1", "Imo", "Jmo", "con_simo_8", "con_simo_5")
expect_true(all(sdata_names %in% names(sdata)))
expect_equivalent(sdata$Xmo_1, as.array(data$x1 - 1))
expect_equivalent(sdata$Xmo_2, as.array(as.numeric(data$x2) - 1))
expect_equal(
as.vector(unname(sdata$Jmo)),
rep(c(max(data$x1) - 1, length(unique(data$x2)) - 1), 4)
)
expect_equal(sdata$con_simo_1, as.array(rep(1, 3)))
prior <- set_prior("dirichlet(1:3)", coef = "mox11",
class = "simo", dpar = "sigma")
sdata <- make_standata(bf(y ~ 1, sigma ~ mo(x1)),
data = data, prior = prior)
expect_equal(sdata$con_simo_sigma_1, as.array(1:3))
prior <- c(
set_prior("normal(0,1)", class = "b", coef = "mox1"),
set_prior("dirichlet(c(1, 0.5, 2))", class = "simo", coef = "mox11"),
prior_(~dirichlet(c(1, 0.5, 2)), class = "simo", coef = "mox1:mox21")
)
sdata <- make_standata(y ~ mo(x1)*mo(x2), data = data, prior = prior)
expect_equal(sdata$con_simo_1, as.array(c(1, 0.5, 2)))
expect_equal(sdata$con_simo_3, as.array(c(1, 0.5, 2)))
# test issue #924 (conditional monotonicity)
prior <- c(prior(dirichlet(c(1, 0.5, 2)), simo, coef = "v"),
prior(dirichlet(c(1,3)), simo, coef = "w"))
sdata <- make_standata(y ~ y*mo(x1, id = "v")*mo(x2, id = "w"),
data, prior = prior)
expect_equal(sdata$con_simo_1, as.array(c(1, 0.5, 2)))
expect_equal(sdata$con_simo_2, as.array(c(1, 3)))
expect_true(!"sdata$con_simo_3" %in% names(sdata))
expect_error(
make_standata(y ~ mo(z), data = data),
"Monotonic predictors must be integers or ordered factors"
)
prior <- c(set_prior("dirichlet(c(1,0.5,2))", class = "simo", coef = "mox21"))
expect_error(
make_standata(y ~ mo(x2), data = data, prior = prior),
"Invalid Dirichlet prior"
)
})
test_that("make_standata returns FCOR covariance matrices", {
data <- data.frame(y = 1:5)
data2 <- list(V = diag(5))
expect_equal(make_standata(y ~ fcor(V), data, data2 = data2)$Mfcor,
data2$V, check.attributes = FALSE)
expect_warning(
expect_error(
make_standata(y~1, data, autocor = cor_fixed(diag(2))),
"Dimensions of 'M' for FCOR terms must be equal"
),
"Using 'cor_brms' objects for 'autocor' is deprecated"
)
})
test_that("make_standata returns data for GAMMs", {
dat <- data.frame(y = rnorm(10), x1 = rnorm(10), x2 = rnorm(10),
x3 = rnorm(10), z = rnorm(10), g = factor(rep(1:2, 5)))
sdata <- make_standata(y ~ s(x1) + z + s(x2, by = x3), data = dat)
expect_equal(sdata$nb_1, 1)
expect_equal(as.vector(sdata$knots_2), 8)
expect_equal(dim(sdata$Zs_1_1), c(10, 8))
expect_equal(dim(sdata$Zs_2_1), c(10, 8))
bform <- bf(y ~ lp, lp ~ s(x1) + z + s(x2, by = x3), nl = TRUE)
sdata <- make_standata(bform, dat)
expect_equal(sdata$nb_lp_1, 1)
expect_equal(as.vector(sdata$knots_lp_2), 8)
expect_equal(dim(sdata$Zs_lp_1_1), c(10, 8))
expect_equal(dim(sdata$Zs_lp_2_1), c(10, 8))
sdata <- make_standata(y ~ g + s(x2, by = g), data = dat)
expect_true(all(c("knots_1", "knots_2") %in% names(sdata)))
# test issue #562
dat$g <- as.character(dat$g)
sdata <- make_standata(y ~ g + s(x2, by = g), data = dat)
expect_true(all(c("knots_1", "knots_2") %in% names(sdata)))
sdata <- make_standata(y ~ t2(x1, x2), data = dat)
expect_equal(sdata$nb_1, 3)
expect_equal(as.vector(sdata$knots_1), c(9, 6, 6))
expect_equal(dim(sdata$Zs_1_1), c(10, 9))
expect_equal(dim(sdata$Zs_1_3), c(10, 6))
expect_error(make_standata(y ~ te(x1, x2), data = dat),
"smooths 'te' and 'ti' are not yet implemented")
})
test_that("make_standata returns correct group ID data", {
form <- bf(count ~ Trt + (1+Trt|3|visit) + (1|patient),
shape ~ (1|3|visit) + (Trt||patient))
sdata <- make_standata(form, data = epilepsy, family = negbinomial())
expect_true(all(c("Z_1_1", "Z_2_2", "Z_3_shape_1", "Z_2_shape_3") %in%
names(sdata)))
form <- bf(count ~ a, sigma ~ (1|3|visit) + (Trt||patient),
a ~ Trt + (1+Trt|3|visit) + (1|patient), nl = TRUE)
sdata <- make_standata(form, data = epilepsy, family = student())
expect_true(all(c("Z_1_sigma_1", "Z_2_a_3", "Z_2_sigma_1",
"Z_3_a_1") %in% names(sdata)))
})
test_that("make_standata handles population-level intercepts", {
dat <- data.frame(y = 10:1, x = 1:10)
sdata <- make_standata(y ~ 0 + x, data = dat)
expect_equal(unname(sdata$X[, 1]), dat$x)
sdata <- make_standata(y ~ x, dat, cumulative(),
control = list(not4stan = TRUE))
expect_equal(unname(sdata$X[, 1]), dat$x)
sdata <- make_standata(y ~ 0 + Intercept + x, data = dat)
expect_equal(unname(sdata$X), cbind(1, dat$x))
})
test_that("make_standata handles category specific effects", {
sdata <- make_standata(rating ~ period + carry + cse(treat),
data = inhaler, family = sratio())
expect_equivalent(sdata$Xcs, matrix(inhaler$treat))
sdata <- make_standata(rating ~ period + carry + cs(treat) + (cs(1)|subject),
data = inhaler, family = acat())
expect_equivalent(sdata$Z_1_3, as.array(rep(1, nrow(inhaler))))
sdata <- make_standata(rating ~ period + carry + (cs(treat)|subject),
data = inhaler, family = cratio())
expect_equivalent(sdata$Z_1_4, as.array(inhaler$treat))
expect_warning(
make_standata(rating ~ 1 + cs(treat), data = inhaler,
family = "cumulative"),
"Category specific effects for this family should be considered experimental"
)
expect_error(make_standata(rating ~ 1 + (treat + cs(1)|subject),
data = inhaler, family = "cratio"),
"category specific effects in separate group-level terms")
})
test_that("make_standata handles wiener diffusion models", {
dat <- data.frame(q = 1:10, resp = sample(0:1, 10, TRUE), x = rnorm(10))
dat$dec <- ifelse(dat$resp == 0, "lower", "upper")
dat$test <- "a"
sdata <- make_standata(q | dec(resp) ~ x, data = dat, family = wiener())
expect_equal(sdata$dec, as.array(dat$resp))
sdata <- make_standata(q | dec(dec) ~ x, data = dat, family = wiener())
expect_equal(sdata$dec, as.array(dat$resp))
expect_error(make_standata(q | dec(test) ~ x, data = dat, family = wiener()),
"Decisions should be 'lower' or 'upper'")
})
test_that("make_standata handles noise-free terms", {
N <- 30
dat <- data.frame(
y = rnorm(N), x = rnorm(N), z = rnorm(N),
xsd = abs(rnorm(N, 1)), zsd = abs(rnorm(N, 1)),
ID = rep(1:5, each = N / 5)
)
sdata <- make_standata(
bf(y ~ me(x, xsd)*me(z, zsd)*x, sigma ~ me(x, xsd)),
data = dat
)
expect_equal(sdata$Xn_1, as.array(dat$x))
expect_equal(sdata$noise_2, as.array(dat$zsd))
expect_equal(unname(sdata$Csp_3), as.array(dat$x))
expect_equal(sdata$Ksp, 6)
expect_equal(sdata$NCme_1, 1)
})
test_that("make_standata handles noise-free terms with grouping factors", {
dat <- data.frame(
y = rnorm(10),
x1 = rep(1:5, each = 2),
sdx = rep(1:5, each = 2),
g = rep(c("b", "c", "a", "d", 1), each = 2)
)
sdata <- make_standata(y ~ me(x1, sdx, gr = g), dat)
expect_equal(unname(sdata$Xn_1), as.array(c(5, 3, 1, 2, 4)))
expect_equal(unname(sdata$noise_1), as.array(c(5, 3, 1, 2, 4)))
dat$sdx[2] <- 10
expect_error(
make_standata(y ~ me(x1, sdx, gr = g), dat),
"Measured values and measurement error should be unique"
)
})
test_that("make_standata handles missing value terms", {
dat = data.frame(y = rnorm(10), x = rnorm(10), g = 1:10)
miss <- c(1, 3, 9)
dat$x[miss] <- NA
bform <- bf(y ~ mi(x)*g) + bf(x | mi() ~ g) + set_rescor(FALSE)
sdata <- make_standata(bform, dat)
expect_equal(sdata$Jmi_x, as.array(miss))
expect_true(all(is.infinite(sdata$Y_x[miss])))
# dots in variable names are correctly handled #452
dat$x.2 <- dat$x
bform <- bf(y ~ mi(x.2)*g) + bf(x.2 | mi() ~ g) + set_rescor(FALSE)
sdata <- make_standata(bform, dat)
expect_equal(sdata$Jmi_x, as.array(miss))
dat$z <- rbeta(10, 1, 1)
dat$z[miss] <- NA
bform <- bf(exp(y) ~ mi(z)*g) + bf(z | mi() ~ g, family = Beta()) +
set_rescor(FALSE)
sdata <- make_standata(bform, dat)
expect_equal(sdata$Jmi_z, as.array(miss))
})
test_that("make_standata handles overimputation", {
dat = data.frame(y = rnorm(10), x = rnorm(10), g = 1:10, sdy = 1)
miss <- c(1, 3, 9)
dat$x[miss] <- dat$sdy[miss] <- NA
bform <- bf(y ~ mi(x)*g) + bf(x | mi(sdy) ~ g) + set_rescor(FALSE)
sdata <- make_standata(bform, dat)
expect_equal(sdata$Jme_x, as.array(setdiff(1:10, miss)))
expect_true(all(is.infinite(sdata$Y_x[miss])))
expect_true(all(is.infinite(sdata$noise_x[miss])))
})
test_that("make_standata handles 'mi' terms with 'subset'", {
dat <- data.frame(
y = rnorm(10), x = c(rnorm(9), NA), z = rnorm(10),
g1 = sample(1:5, 10, TRUE), g2 = 10:1, g3 = 1:10,
s = c(FALSE, rep(TRUE, 9))
)
bform <- bf(y ~ mi(x, idx = g1)) +
bf(x | mi() + index(g2) + subset(s) ~ 1) +
set_rescor(FALSE)
sdata <- make_standata(bform, dat)
expect_true(all(sdata$idxl_y_x_1 %in% 9:5))
# test a bunch of errors
# fails on CRAN for some reason
# bform <- bf(y ~ mi(x, idx = g1)) +
# bf(x | mi() + index(g3) + subset(s) ~ 1) +
# set_rescor(FALSE)
# expect_error(make_standata(bform, dat),
# "Could not match all indices in response 'x'"
# )
bform <- bf(y ~ mi(x, idx = g1)) +
bf(x | mi() + subset(s) ~ 1) +
set_rescor(FALSE)
expect_error(make_standata(bform, dat),
"Response 'x' needs to have an 'index' addition term"
)
bform <- bf(y ~ mi(x)) +
bf(x | mi() + subset(s) + index(g2) ~ 1) +
set_rescor(FALSE)
expect_error(make_standata(bform, dat),
"mi() terms of subsetted variables require the 'idx' argument",
fixed = TRUE
)
bform <- bf(y | mi() ~ mi(x, idx = g1)) +
bf(x | mi() + subset(s) + index(g2) ~ mi(y)) +
set_rescor(FALSE)
expect_error(make_standata(bform, dat),
"mi() terms in subsetted formulas require the 'idx' argument",
fixed = TRUE
)
})
test_that("make_standata handles multi-membership models", {
dat <- data.frame(y = rnorm(10), g1 = c(7:2, rep(10, 4)),
g2 = 1:10, w1 = rep(1, 10),
w2 = rep(abs(rnorm(10))))
sdata <- make_standata(y ~ (1|mm(g1,g2,g1,g2)), data = dat)
expect_true(all(paste0(c("W_1_", "J_1_"), 1:4) %in% names(sdata)))
expect_equal(sdata$W_1_4, as.array(rep(0.25, 10)))
expect_equal(unname(sdata$Z_1_1_1), as.array(rep(1, 10)))
expect_equal(unname(sdata$Z_1_1_2), as.array(rep(1, 10)))
# this checks whether combintation of factor levels works as intended
expect_equal(sdata$J_1_1, as.array(c(6, 5, 4, 3, 2, 1, 7, 7, 7, 7)))
expect_equal(sdata$J_1_2, as.array(c(8, 1, 2, 3, 4, 5, 6, 9, 10, 7)))
sdata <- make_standata(y ~ (1|mm(g1,g2, weights = cbind(w1, w2))), dat)
expect_equal(sdata$W_1_1, as.array(dat$w1 / (dat$w1 + dat$w2)))
# tests mmc terms
sdata <- make_standata(y ~ (1+mmc(w1, w2)|mm(g1,g2)), data = dat)
expect_equal(unname(sdata$Z_1_2_1), as.array(dat$w1))
expect_equal(unname(sdata$Z_1_2_2), as.array(dat$w2))
expect_error(
make_standata(y ~ (mmc(w1, w2, y)|mm(g1,g2)), data = dat),
"Invalid term 'mmc(w1, w2, y)':", fixed = TRUE
)
expect_error(
make_standata(y ~ (mmc(w1, w2)*y|mm(g1,g2)), data = dat),
"The term 'mmc(w1, w2):y' is invalid", fixed = TRUE
)
# tests if ":" works in multi-membership models
sdata <- make_standata(y ~ (1|mm(w1:g1,w1:g2)), dat)
expect_true(all(c("J_1_1", "J_1_2") %in% names(sdata)))
})
test_that("by variables in grouping terms are handled correctly", {
gvar <- c("1A", "1B", "2A", "2B", "3A", "3B", "10", "100", "2", "3")
gvar <- rep(gvar, each = 10)
g_order <- order(gvar)
byvar <- c(0, 4.5, 3, 2, "x 1")
byvar <- factor(rep(byvar, each = 20))
dat <- data.frame(
y = rnorm(100), x = rnorm(100),
g = gvar,
g2 = gvar[g_order],
z = byvar,
z2 = byvar[g_order],
z3 = factor(1:2)
)
sdata <- make_standata(y ~ x + (x | gr(g, by = z)), dat)
expect_equal(sdata$Nby_1, 5)
expect_equal(sdata$Jby_1, as.array(c(2, 2, 1, 1, 5, 4, 4, 5, 3, 3)))
sdata <- make_standata(y ~ x + (x | mm(g, g2, by = cbind(z, z2))), dat)
expect_equal(sdata$Nby_1, 5)
expect_equal(sdata$Jby_1, as.array(c(2, 2, 1, 1, 5, 4, 4, 5, 3, 3)))
expect_error(make_standata(y ~ x + (1|gr(g, by = z3)), dat),
"Some levels of 'g' correspond to multiple levels of 'z3'")
})
test_that("make_standata handles calls to the 'poly' function", {
dat <- data.frame(y = rnorm(10), x = rnorm(10))
expect_equal(colnames(make_standata(y ~ 1 + poly(x, 3), dat)$X),
c("Intercept", "polyx31", "polyx32", "polyx33"))
})
test_that("make_standata allows fixed distributional parameters", {
dat <- list(y = 1:10)
expect_equal(make_standata(bf(y ~ 1, nu = 3), dat, student())$nu, 3)
expect_equal(make_standata(y ~ 1, dat, acat())$disc, 1)
expect_error(make_standata(bf(y ~ 1, bias = 0.5), dat),
"Invalid fixed parameters: 'bias'")
})
test_that("Cell-mean coding can be disabled", {
df <- data.frame(y = 1:10, g = rep(c("a", "b"), 5))
bform <- bf(y ~ g) +
lf(disc ~ 0 + g + (0 + g | y), cmc = FALSE) +
cumulative()
sdata <- make_standata(bform, df)
target <- matrix(rep(0:1, 5), dimnames = list(1:10, "gb"))
expect_equal(sdata$X_disc, target)
expect_equal(unname(sdata$Z_1_disc_1), as.array(rep(0:1, 5)))
expect_true(!"Z_1_disc_2" %in% names(sdata))
bform <- bf(y ~ 0 + g + (1 | y), cmc = FALSE)
sdata <- make_standata(bform, df)
expect_equal(sdata$X, target)
expect_equal(unname(sdata$Z_1_1), as.array(rep(1, 10)))
})
test_that("make_standata correctly includes offsets", {
data <- data.frame(y = rnorm(10), x = rnorm(10), c = 1)
sdata <- make_standata(bf(y ~ x + offset(c), sigma ~ offset(c + 1)), data)
expect_equal(sdata$offsets, as.array(data$c))
expect_equal(sdata$offsets_sigma, as.array(data$c + 1))
sdata <- make_standata(y ~ x + offset(c) + offset(x), data)
expect_equal(sdata$offsets, as.array(data$c + data$x))
})
test_that("make_standata includes data for mixture models", {
data <- data.frame(y = rnorm(10), x = rnorm(10), c = 1)
form <- bf(y ~ x, mu1 ~ 1, family = mixture(gaussian, gaussian))
sdata <- make_standata(form, data)
expect_equal(sdata$con_theta, as.array(c(1, 1)))
expect_equal(dim(sdata$X_mu1), c(10, 1))
expect_equal(dim(sdata$X_mu2), c(10, 2))
form <- bf(y ~ x, family = mixture(gaussian, gaussian))
sdata <- make_standata(form, data, prior = prior(dirichlet(10, 2), theta))
expect_equal(sdata$con_theta, as.array(c(10, 2)))
form <- bf(y ~ x, theta1 = 1, theta2 = 3, family = mixture(gaussian, gaussian))
sdata <- make_standata(form, data)
expect_equal(sdata$theta1, 1/4)
expect_equal(sdata$theta2, 3/4)
})
test_that("make_standata includes data for Gaussian processes", {
dat <- data.frame(y = rnorm(10), x1 = rnorm(10),
z = factor(c(2, 2, 2, 3, 4, rep(5, 5))))
sdata <- make_standata(y ~ gp(x1), dat)
expect_equal(max(sdata$Xgp_1) - min(sdata$Xgp_1), 1)
sdata <- make_standata(y ~ gp(x1, scale = FALSE), dat)
expect_equal(max(sdata$Xgp_1) - min(sdata$Xgp_1), max(dat$x1) - min(dat$x1))
sdata <- SW(make_standata(y ~ gp(x1, by = z, gr = TRUE, scale = FALSE), dat))
expect_equal(sdata$Igp_1_2, as.array(4))
expect_equal(sdata$Jgp_1_4, as.array(1:5))
expect_equal(sdata$Igp_1_4, as.array(6:10))
sdata <- SW(make_standata(y ~ gp(x1, by = y, gr = TRUE), dat))
expect_equal(sdata$Cgp_1, as.array(dat$y))
})
test_that("make_standata includes data for approximate Gaussian processes", {
dat <- data.frame(y = rnorm(10), x1 = sample(1:10, 10),
z = factor(c(2, 2, 2, 3, 4, rep(5, 5))))
sdata <- make_standata(y ~ gp(x1, k = 5, c = 5/4), dat)
expect_equal(sdata$NBgp_1, 5)
expect_equal(dim(sdata$Xgp_1), c(10, 5))
expect_equal(dim(sdata$slambda_1), c(5, 1))
sdata <- SW(make_standata(y ~ gp(x1, by = z, k = 5, c = 5/4, scale = FALSE), dat))
expect_equal(sdata$Igp_1_2, as.array(4))
expect_equal(sdata$Cgp_1_2, as.array(1))
expect_equal(sdata$Igp_1_4, as.array(6:10))
})
test_that("make_standata includes data for SAR models", {
dat <- data.frame(y = rnorm(10), x = rnorm(10))
W <- matrix(0, nrow = 10, ncol = 10)
dat2 <- list(W = W)
sdata <- make_standata(y ~ x + sar(W), data = dat, data2 = dat2)
expect_equal(dim(sdata$M), rep(nrow(W), 2))
dat2 <- list(W = matrix(0, 2, 2))
expect_error(
make_standata(y ~ x + sar(W), data = dat, data2 = dat2),
"Dimensions of 'M' for SAR terms must be equal"
)
})
test_that("make_standata includes data for CAR models", {
dat = data.frame(y = rnorm(10), x = rnorm(10), obs = 1:10)
edges <- cbind(1:10, 10:1)
W <- matrix(0, nrow = 10, ncol = 10)
for (i in seq_len(nrow(edges))) {
W[edges[i, 1], edges[i, 2]] <- 1
}
rownames(W) <- 1:nrow(W)
dat2 <- list(W = W)
sdata <- make_standata(y ~ x + car(W, gr = obs), dat, data2 = dat2)
expect_equal(sdata$Nloc, 10)
expect_equal(unname(sdata$Nneigh), rep(1, 10))
expect_equal(unname(sdata$edges1), as.array(10:6))
expect_equal(unname(sdata$edges2), as.array(1:5))
sdata_old <- SW(make_standata(y ~ x, dat, autocor = cor_car(W)))
expect_equal(sdata, sdata_old)
rownames(dat2$W) <- c("a", 2:9, "b")
dat$group <- rep(c("a", "b"), each = 5)
sdata <- make_standata(y ~ x + car(W, gr = group), dat, data2 = dat2,
prior = prior(horseshoe(), class = sdcar))
expect_equal(sdata$Nloc, 2)
expect_equal(sdata$edges1, as.array(2))
expect_equal(sdata$edges2, as.array(1))
expect_equal(sdata$Kscales, 1)
sdata <- make_standata(y ~ x + car(W, group, type = "bym2"),
data = dat, data2 = dat2)
expect_equal(length(sdata$car_scale), 1L)
dat2$W[1, 10] <- 4
dat2$W[10, 1] <- 4
expect_message(make_standata(y ~ car(W, gr = group), dat, data2 = dat2),
"Converting all non-zero values in 'M' to 1")
# test error messages
rownames(dat2$W) <- c(1:9, "a")
expect_error(make_standata(y ~ car(W, gr = group), dat, data2 = dat2),
"Row names of 'M' for CAR terms do not match")
rownames(dat2$W) <- NULL
expect_error(make_standata(y ~ car(W, gr = group), dat, data2 = dat2),
"Row names are required for 'M'")
dat2$W[1, 10] <- 0
expect_error(make_standata(y ~ car(W), dat, data2 = dat2),
"'M' for CAR terms must be symmetric")
dat2$W[10, 1] <- 0
expect_error(SW(make_standata(y ~ x + car(W), dat, data2 = dat2)),
"all locations should have at least one neighbor")
})
test_that("make_standata includes data of special priors", {
dat <- data.frame(y = 1:10, x1 = rnorm(10), x2 = rnorm(10),
g = rep(1:2, each = 5), x3 = sample(1:5, 10, TRUE))
# horseshoe prior
hs <- horseshoe(7, scale_global = 2, df_global = 3,
df_slab = 6, scale_slab = 3)
sdata <- make_standata(y ~ x1*x2, data = dat,
prior = set_prior(hs))
expect_equal(sdata$hs_df, 7)
expect_equal(sdata$hs_df_global, 3)
expect_equal(sdata$hs_df_slab, 6)
expect_equal(sdata$hs_scale_global, 2)
expect_equal(sdata$hs_scale_slab, 3)
hs <- horseshoe(par_ratio = 0.1)
sdata <- make_standata(y ~ x1*x2, data = dat, prior = set_prior(hs))
expect_equal(sdata$hs_scale_global, 0.1 / sqrt(nrow(dat)))
# R2D2 prior
sdata <- make_standata(y ~ x1*x2, data = dat,
prior = prior(R2D2(0.5, 10)))
expect_equal(sdata$R2D2_mean_R2, 0.5)
expect_equal(sdata$R2D2_prec_R2, 10)
expect_equal(sdata$R2D2_cons_D2, as.array(rep(0.5, 3)))
# horseshoe and R2D2 prior applied in a non-linear model
hs_a1 <- horseshoe(7, scale_global = 2, df_global = 3)
R2D2_a2 <- R2D2(0.5, 10)
sdata <- make_standata(
bf(y ~ a1 + a2, a1 ~ x1, a2 ~ 0 + x2, nl = TRUE),
data = dat, sample_prior = TRUE,
prior = c(set_prior(hs_a1, nlpar = "a1"),
set_prior(R2D2_a2, nlpar = "a2"))
)
expect_equal(sdata$hs_df_a1, 7)
expect_equal(sdata$R2D2_mean_R2_a2, 0.5)
bform <- bf(y ~ x1*mo(x3) + (1|g) + gp(x3) + s(x2) +
arma(p = 2, q = 2, gr = g))
bprior <- prior(R2D2(cons_D2 = 11:1, main = TRUE), class = b) +
prior(R2D2(), class = sd) +
prior(R2D2(), class = sds) +
prior(R2D2(), class = sdgp) +
prior(R2D2(), class = ar) +
prior(R2D2(), class = ma)
sdata <- make_standata(bform, data = dat, prior = bprior)
expect_equal(sdata$Kscales, 11)
expect_equal(sdata$R2D2_cons_D2, as.array(11:1))
})
test_that("dots in formula are correctly expanded", {
dat <- data.frame(y = 1:10, x1 = 1:10, x2 = 1:10)
sdata <- make_standata(y ~ ., dat)
expect_equal(colnames(sdata$X), c("Intercept", "x1", "x2"))
})
test_that("argument 'stanvars' is handled correctly", {
bprior <- prior(normal(mean_intercept, 10), class = "Intercept")
mean_intercept <- 5
stanvars <- stanvar(mean_intercept)
sdata <- make_standata(count ~ Trt, data = epilepsy,
prior = bprior, stanvars = stanvars)
expect_equal(sdata$mean_intercept, 5)
# define a multi_normal prior with known covariance matrix
bprior <- prior(multi_normal(M, V), class = "b")
stanvars <- stanvar(rep(0, 2), "M", scode = " vector[K] M;") +
stanvar(diag(2), "V", scode = " matrix[K, K] V;")
sdata <- make_standata(count ~ Trt + zBase, epilepsy,
prior = bprior, stanvars = stanvars)
expect_equal(sdata$M, rep(0, 2))
expect_equal(sdata$V, diag(2))
})
test_that("addition arguments 'vint' and 'vreal' work correctly", {
dat <- data.frame(size = 10, y = sample(0:10, 20, TRUE), x = rnorm(20))
beta_binomial2 <- custom_family(
"beta_binomial2",
dpars = c("mu", "tau"),
links = c("logit", "log"),
lb = c(NA, 0),
type = "int",
vars = c("vint1[n]", "vreal1[n]")
)
sdata <- make_standata(
y | vint(size) + vreal(x, size) ~ 1,
data = dat, family = beta_binomial2,
)
expect_equal(sdata$vint1, as.array(rep(10, 20)))
expect_equal(sdata$vreal1, as.array(dat$x))
expect_equal(sdata$vreal2, as.array(rep(10, 20)))
})
test_that("reserved variables 'Intercept' is handled correctly", {
dat <- data.frame(y = 1:10)
expect_warning(
sdata <- make_standata(y ~ 0 + intercept, dat),
"Reserved variable name 'intercept' is deprecated."
)
expect_true(all(sdata$X[, "intercept"] == 1))
sdata <- make_standata(y ~ 0 + Intercept, dat)
expect_true(all(sdata$X[, "Intercept"] == 1))
})
test_that("data for multinomial and dirichlet models is correct", {
N <- 15
dat <- as.data.frame(rdirichlet(N, c(3, 2, 1)))
names(dat) <- c("y1", "y2", "y3")
dat$t1 <- round(dat$y1 * rpois(N, 10))
dat$t2 <- round(dat$y2 * rpois(N, 10))
dat$t3 <- round(dat$y3 * rpois(N, 10))
dat$x <- rnorm(N)
dat$y <- with(dat, cbind(y1, y2, y3))
dat$t <- with(dat, cbind(t1, t2, t3))
dat$size <- rowSums(dat$t)
sdata <- make_standata(t | trials(size) ~ x, dat, multinomial())
expect_equal(sdata$trials, as.array(dat$size))
expect_equal(sdata$ncat, 3)
expect_equal(sdata$Y, unname(dat$t))
sdata <- make_standata(y ~ x, data = dat, family = dirichlet())
expect_equal(sdata$ncat, 3)
expect_equal(sdata$Y, unname(dat$y))
expect_error(
make_standata(t | trials(10) ~ x, data = dat, family = multinomial()),
"Number of trials does not match the number of events"
)
expect_error(make_standata(t ~ x, data = dat, family = dirichlet()),
"Response values in simplex models must sum to 1")
})
test_that("make_standata handles cox models correctly", {
data <- data.frame(y = rexp(100), x = rnorm(100))
bform <- bf(y ~ x)
bprior <- prior(dirichlet(3), sbhaz)
sdata <- make_standata(bform, data, brmsfamily("cox"), prior = bprior)
expect_equal(dim(sdata$Zbhaz), c(100, 5))
expect_equal(dim(sdata$Zcbhaz), c(100, 5))
expect_equal(sdata$con_sbhaz, as.array(rep(3, 5)))
sdata <- make_standata(bform, data, brmsfamily("cox", bhaz = list(df = 6)))
expect_equal(dim(sdata$Zbhaz), c(100, 6))
expect_equal(dim(sdata$Zcbhaz), c(100, 6))
})
test_that("make_standata handles addition term 'rate' is correctly", {
data <- data.frame(y = rpois(10, 1), x = rnorm(10), time = 1:10)
sdata <- make_standata(y | rate(time) ~ x, data, poisson())
expect_equal(sdata$denom, as.array(data$time))
})
test_that("make_standata handles grouped ordinal thresholds correctly", {
dat <- data.frame(
y = c(1:5, 1:4, 4),
gr = rep(c("a", "b"), each = 5),
th = rep(5:6, each = 5),
x = rnorm(10)
)
# thresholds without a grouping factor
sdata <- make_standata(y ~ x, dat, cumulative())
expect_equal(sdata$nthres, 4)
sdata <- make_standata(y | thres(5) ~ x, dat, cumulative())
expect_equal(sdata$nthres, 5)
expect_error(
make_standata(y | thres(th) ~ x, dat, cumulative()),
"Number of thresholds needs to be a single value"
)
# thresholds with a grouping factor
sdata <- make_standata(y | thres(th, gr) ~ x, dat, cumulative())
expect_equal(sdata$nthres, as.array(c(5, 6)))
expect_equal(sdata$ngrthres, 2)
expect_equal(unname(sdata$Jthres[1, ]), c(1, 5))
expect_equal(unname(sdata$Jthres[10, ]), c(6, 11))
sdata <- make_standata(y | thres(gr = gr) ~ x, dat, cumulative())
expect_equal(sdata$nthres, as.array(c(4, 3)))
expect_equal(sdata$ngrthres, 2)
sdata <- make_standata(y | thres(6, gr = gr) ~ x, dat, cumulative())
expect_equal(sdata$nthres, as.array(c(6, 6)))
expect_equal(sdata$ngrthres, 2)
})
test_that("information for threading is handled correctly", {
dat <- data.frame(y = 1:10)
sdata <- make_standata(y ~ 1, dat, threads = threading(2, grainsize = 3))
expect_equal(sdata$grainsize, 3)
})
test_that("variables in data2 can be used in population-level effects", {
dat <- data.frame(y = 1:10, x1 = rnorm(10), x2 = rnorm(10), x3 = rnorm(10))
foo <- function(..., idx = NULL) {
out <- cbind(...)
if (!is.null(idx)) {
out <- out[, idx, drop = FALSE]
}
out
}
sdata <- make_standata(y ~ foo(x1, x2, x3, idx = id), data = dat,
data2 = list(id = c(3, 1)))
target <- c("Intercept", "foox1x2x3idxEQidx3", "foox1x2x3idxEQidx1")
expect_equal(colnames(sdata$X), target)
expect_equivalent(sdata$X[, 2], dat$x3)
expect_equivalent(sdata$X[, 3], dat$x1)
})
test_that("NAs are allowed in unused interval censoring variables", {
dat <- data.frame(y = rnorm(10), ce = c(1, rep(2, 9)))
dat$y2 <- dat$y + 2
dat$y2[1] <- NA
sdata <- make_standata(y | cens(ce, y2 = y2) ~ 1, data = dat)
expect_equal(sdata$N, 10L)
expect_equal(sdata$rcens[1], 0)
dat$ce[1] <- 2
expect_error(
make_standata(y | cens(ce, y2 = y2) ~ 1, data = dat),
"'y2' should not be NA for interval censored observations"
)
})
test_that("drop_unused_factor levels works correctly", {
dat <- data.frame(y = rnorm(10), x = factor(c("a", "b"), levels = c("a", "b", "c")))
# should drop level "c"
sdata <- make_standata(y ~ x, data = dat)
expect_equal(colnames(sdata$X), c("Intercept", "xb"))
# should not drop level "c"
sdata <- make_standata(y ~ x, data = dat, drop_unused_levels = FALSE)
expect_equal(colnames(sdata$X), c("Intercept", "xb", "xc"))
})
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.