tests/testthat/test-Model-generators.R
In StatisticsNZ/demest: Bayesian Demographic Estimation and Forecasting
context("Model-generators")
n.test <- 5
test.identity <- FALSE
test.extended <- FALSE
## addAg ########################################################################
## SpecAgPlaceholder
test_that("addAg works with BinomialVarying and SpecAgPlaceholder", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Binomial(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
ans.obtained <- addAg(model = model,
aggregate = new("SpecAgPlaceholder"))
ans.expected <- model
expect_identical(ans.obtained, ans.expected)
expect_is(ans.obtained, "BinomialVarying")
})
test_that("addAg works with NormalVaryingVarsigmaKnown and SpecAgPlaceholder", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region, sd = 0.5))
set.seed(1)
model <- initialModel(spec, y = y, weights = weights)
ans.obtained <- addAg(model = model,
aggregate = new("SpecAgPlaceholder"))
ans.expected <- model
expect_identical(ans.obtained, ans.expected)
expect_is(ans.obtained, "NormalVaryingVarsigmaKnown")
})
test_that("addAg works with NormalVaryingVarsigmaUnknown and SpecAgPlaceholder", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, weights = weights)
ans.obtained <- addAg(model = model,
aggregate = new("SpecAgPlaceholder"))
ans.expected <- model
expect_identical(ans.obtained, ans.expected)
expect_is(ans.obtained, "NormalVaryingVarsigmaUnknown")
})
test_that("addAg works with PoissonVaryingNotUseExp and SpecAgPlaceholder", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE))
set.seed(1)
model <- initialModel(spec, y = y, exposure = NULL)
ans.obtained <- addAg(model = model,
aggregate = new("SpecAgPlaceholder"))
ans.expected <- model
expect_identical(ans.obtained, ans.expected)
expect_is(ans.obtained, "PoissonVaryingNotUseExp")
})
test_that("addAg works with PoissonVaryingNotUseExp and SpecAgPlaceholder", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
exposure <- y + 1L
spec <- Model(y ~ Poisson(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
ans.obtained <- addAg(model = model,
aggregate = new("SpecAgPlaceholder"))
ans.expected <- model
expect_identical(ans.obtained, ans.expected)
expect_is(ans.obtained, "PoissonVaryingUseExp")
})
## SpecAgCertain
test_that("addAg works with BinomialVarying and SpecAgCertain", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Binomial(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
aggregate <- AgCertain(value = 0.3)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = exposure)
expect_is(ans.obtained, "BinomialVaryingAgCertain")
expect_identical(ans.obtained@slotsToExtract,
new("BinomialVaryingAgCertain")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("BinomialVaryingAgCertain")@iMethodModel)
})
test_that("addAg works with NormalVaryingVarsigmaKnown and SpecAgCertain", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region, sd = 0.5))
set.seed(1)
model <- initialModel(spec, y = y, weights = weights)
aggregate <- AgCertain(value = 0.3)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = weights)
expect_is(ans.obtained, "NormalVaryingVarsigmaKnownAgCertain")
expect_identical(ans.obtained@slotsToExtract,
new("NormalVaryingVarsigmaKnownAgCertain")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("NormalVaryingVarsigmaKnownAgCertain")@iMethodModel)
})
test_that("addAg works with NormalVaryingVarsigmaUnknown and SpecAgCertain", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, weights = weights)
aggregate <- AgCertain(value = 0.3)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = weights)
expect_is(ans.obtained, "NormalVaryingVarsigmaUnknownAgCertain")
expect_identical(ans.obtained@slotsToExtract,
new("NormalVaryingVarsigmaUnknownAgCertain")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("NormalVaryingVarsigmaUnknownAgCertain")@iMethodModel)
})
test_that("addAg works with PoissonVaryingNotUseExp and SpecAgCertain", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
defaultWeights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE))
set.seed(1)
model <- initialModel(spec, y = y, exposure = NULL)
aggregate <- AgCertain(value = 0.3)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = defaultWeights)
expect_is(ans.obtained, "PoissonVaryingNotUseExpAgCertain")
expect_identical(ans.obtained@slotsToExtract,
new("PoissonVaryingNotUseExpAgCertain")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("PoissonVaryingNotUseExpAgCertain")@iMethodModel)
})
test_that("addAg works with PoissonVaryingNotUseExp and SpecAgCertain", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
aggregate <- AgCertain(value = 0.3)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = exposure)
expect_is(ans.obtained, "PoissonVaryingUseExpAgCertain")
expect_identical(ans.obtained@slotsToExtract,
new("PoissonVaryingUseExpAgCertain")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("PoissonVaryingUseExpAgCertain")@iMethodModel)
})
## SpecAgNormal
test_that("addAg works with BinomialVarying and SpecAgNormal", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Binomial(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
aggregate <- AgNormal(value = 0.3, sd = 0.2)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = exposure)
expect_is(ans.obtained, "BinomialVaryingAgNormal")
expect_identical(ans.obtained@slotsToExtract,
new("BinomialVaryingAgNormal")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("BinomialVaryingAgNormal")@iMethodModel)
})
test_that("addAg works with NormalVaryingVarsigmaKnown and SpecAgNormal", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region, sd = 0.5))
set.seed(1)
model <- initialModel(spec, y = y, weights = weights)
aggregate <- AgNormal(value = 0.3, sd = 0.2)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = weights)
expect_is(ans.obtained, "NormalVaryingVarsigmaKnownAgNormal")
expect_identical(ans.obtained@slotsToExtract,
new("NormalVaryingVarsigmaKnownAgNormal")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("NormalVaryingVarsigmaKnownAgNormal")@iMethodModel)
})
test_that("addAg works with NormalVaryingVarsigmaUnknown and SpecAgNormal", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, weights = weights)
aggregate <- AgNormal(value = 0.3, sd = 0.2)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = weights)
expect_is(ans.obtained, "NormalVaryingVarsigmaUnknownAgNormal")
expect_identical(ans.obtained@slotsToExtract,
new("NormalVaryingVarsigmaUnknownAgNormal")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("NormalVaryingVarsigmaUnknownAgNormal")@iMethodModel)
})
test_that("addAg works with PoissonVaryingNotUseExp and SpecAgNormal", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
defaultWeights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE))
set.seed(1)
model <- initialModel(spec, y = y, exposure = NULL)
aggregate <- AgNormal(value = 0.3, sd = 0.2)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = defaultWeights)
expect_is(ans.obtained, "PoissonVaryingNotUseExpAgNormal")
expect_identical(ans.obtained@slotsToExtract,
new("PoissonVaryingNotUseExpAgNormal")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("PoissonVaryingNotUseExpAgNormal")@iMethodModel)
})
test_that("addAg works with PoissonVaryingUseExp and SpecAgNormal", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
aggregate <- AgNormal(value = 0.3, sd = 0.2)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = exposure)
expect_is(ans.obtained, "PoissonVaryingUseExpAgNormal")
expect_identical(ans.obtained@slotsToExtract,
new("PoissonVaryingUseExpAgNormal")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("PoissonVaryingUseExpAgNormal")@iMethodModel)
})
## SpecAgPoisson
test_that("addAg works with BinomialVarying and SpecAgPoisson", {
addAg <- demest:::addAg
expect_error(addAg(model = new("BinomialVarying"),
aggregate = new("SpecAgPoisson")),
"Poisson model for aggregate values can only be used with Poisson likelihood")
})
test_that("addAg works with NormalVaryingVarsigmaKnown and SpecAgPoisson", {
addAg <- demest:::addAg
expect_error(addAg(model = new("NormalVaryingVarsigmaKnown"),
aggregate = new("SpecAgPoisson")),
"Poisson model for aggregate values can only be used with Poisson likelihood")
})
test_that("addAg works with NormalVaryingVarsigmaUnknown and SpecAgPoisson", {
addAg <- demest:::addAg
expect_error(addAg(model = new("NormalVaryingVarsigmaUnknown"),
aggregate = new("SpecAgPoisson")),
"Poisson model for aggregate values can only be used with Poisson likelihood")
})
test_that("addAg works with PoissonVaryingNotUseExp and SpecAgPoisson", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
defaultWeights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE))
set.seed(1)
model <- initialModel(spec, y = y, exposure = NULL)
aggregate <- AgPoisson(value = 0.3)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = defaultWeights)
expect_is(ans.obtained, "PoissonVaryingNotUseExpAgPoisson")
expect_identical(ans.obtained@slotsToExtract,
new("PoissonVaryingNotUseExpAgPoisson")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("PoissonVaryingNotUseExpAgPoisson")@iMethodModel)
})
test_that("addAg works with PoissonVaryingNotUseExp and SpecAgPoisson", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
aggregate <- AgPoisson(value = 0.3)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = exposure)
expect_is(ans.obtained, "PoissonVaryingUseExpAgPoisson")
expect_identical(ans.obtained@slotsToExtract,
new("PoissonVaryingUseExpAgPoisson")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("PoissonVaryingUseExpAgPoisson")@iMethodModel)
})
## SpecAgFun
test_that("addAg works with BinomialVarying and SpecAgFun", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Binomial(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
aggregate <- AgFun(value = 0.3, sd = 0.2, FUN = function(x, weights) 1)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = exposure)
expect_is(ans.obtained, "BinomialVaryingAgFun")
expect_identical(ans.obtained@slotsToExtract,
new("BinomialVaryingAgFun")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("BinomialVaryingAgFun")@iMethodModel)
})
test_that("addAg works with NormalVaryingVarsigmaKnown and SpecAgFun", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region, sd = 0.5))
set.seed(1)
model <- initialModel(spec, y = y, weights = weights)
aggregate <- AgFun(value = 0.3, sd = 0.2, FUN = function(x, weights) 1)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = weights)
expect_is(ans.obtained, "NormalVaryingVarsigmaKnownAgFun")
expect_identical(ans.obtained@slotsToExtract,
new("NormalVaryingVarsigmaKnownAgFun")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("NormalVaryingVarsigmaKnownAgFun")@iMethodModel)
})
test_that("addAg works with NormalVaryingVarsigmaUnknown and SpecAgFun", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, weights = weights)
aggregate <- AgFun(value = 0.3, sd = 0.2, FUN = function(x, weights) 1)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = weights)
expect_is(ans.obtained, "NormalVaryingVarsigmaUnknownAgFun")
expect_identical(ans.obtained@slotsToExtract,
new("NormalVaryingVarsigmaUnknownAgFun")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("NormalVaryingVarsigmaUnknownAgFun")@iMethodModel)
})
test_that("addAg works with PoissonVaryingUseExp and SpecAgFun", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rbinom(n = 20, size = 20, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])),
dimscales = c(age = "Intervals"))
defaultWeights <- Counts(array(rbinom(n = 20, size = 20, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])),
dimscales = c(age = "Intervals"))
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE))
set.seed(1)
model <- initialModel(spec, y = y, exposure = NULL)
aggregate <- AgFun(value = 0.3, sd = 0.2, FUN = function(x, weights) 1)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = defaultWeights)
expect_is(ans.obtained, "PoissonVaryingNotUseExpAgFun")
expect_identical(ans.obtained@slotsToExtract,
new("PoissonVaryingNotUseExpAgFun")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("PoissonVaryingNotUseExpAgFun")@iMethodModel)
})
test_that("addAg works with PoissonVaryingUseExp and SpecAgFun", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
aggregate <- AgFun(value = 0.3, sd = 0.2, FUN = function(x, weights) 1)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = exposure)
expect_is(ans.obtained, "PoissonVaryingUseExpAgFun")
expect_identical(ans.obtained@slotsToExtract,
new("PoissonVaryingUseExpAgFun")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("PoissonVaryingUseExpAgFun")@iMethodModel)
})
## initialModel - Varying ##################################################################
test_that("initialModel creates object of class BinomialVarying from valid inputs", {
initialModel <- demest:::initialModel
makeLinearBetas <- demest:::makeLinearBetas
jitterBetas <- demest:::jitterBetas
makePriors <- demest:::makePriors
## main effect
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Binomial(mean ~ age + region))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
m <- sum(y) / sum(exposure)
nu <- m * (1 - m)
sigma <- new("Scale", runif(1, 0, nu))
theta <- rbeta(n = 20,
shape1 = m * (nu/sigma@.Data^2 - 1) + y,
shape2 = (1-m) * (nu/sigma@.Data^2 - 1) + exposure - y)
logit.theta <- array(log(theta / (1 - theta)), dim = dim(y), dimnames = dimnames(y))
betas <- makeLinearBetas(theta = logit.theta, formula = prob ~ age + region)
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = NULL,
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priors)
expect_identical(x@sigma, sigma)
expect_identical(x@theta, theta)
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@ASigma@.Data, 1)
expect_identical(x@nuSigma@.Data, 7)
expect_identical(x@sigmaMax@.Data, qhalft(p = 0.999, df = 7, scale = 1))
expect_identical(x@lower, -Inf)
expect_identical(x@upper, Inf)
expect_identical(x@tolerance, 1e-5)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@nFailedPropTheta, new("Counter", 0L))
## intercept only
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Binomial(mean ~ 1))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
m <- sum(y) / sum(exposure)
nu <- m * (1 - m)
sigma <- runif(1, 0, nu)
theta <- rbeta(n = 20,
shape1 = m * (nu/sigma^2 -1) + y,
shape2 = (1-m) * (nu/sigma^2 - 1) + exposure - y)
sigma <- new("Scale", sigma)
betas <- list("(Intercept)" = mean(log(theta / (1 - theta))))
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L),
y = y,
sY = NULL,
strucZeroArray = strucZeroArray)
betas <- unname(betas)
betas <- jitterBetas(betas, priors)
expect_identical(x@sigma, sigma)
expect_identical(x@theta, theta)
expect_identical(x@betas, betas)
expect_identical(sapply(x@priorsBetas, class), "ExchFixed")
## lower and upper specified
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Binomial(mean ~ age + region), lower = 0.4, upper = 0.8)
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(all(x@theta >= 0.4))
expect_true(all(x@theta <= 0.8))
expect_equal(x@lower, log(0.4 / 0.6))
expect_equal(x@upper, log(0.8 / 0.2))
## missing values
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y[1:4] <- NA
exposure[1:3] <- NA
spec <- Model(y ~ Binomial(mean ~ age + region))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
m <- sum(y[5:20]) / sum(exposure[5:20])
nu <- m * (1 - m)
sigma <- runif(1, 0, nu)
sigma <- new("Scale", sigma)
y[1:4] <- 0
exposure[1:4] <- 0
theta <- rbeta(n = 20,
shape1 = m * (nu/sigma@.Data^2 - 1) + y,
shape2 = (1-m) * (nu/sigma@.Data^2 - 1) + exposure - y)
expect_identical(x@sigma, sigma)
expect_identical(x@theta, theta)
## saturated model
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Binomial(mean ~ age * region),
age:region ~ Exch(error = Error(scale = HalfT(df = 3, scale = 0.1, max = 0.6))))
x <- initialModel(spec, y = y, exposure = exposure)
expect_identical(x@ASigma, x@priorsBetas[[4]]@ATau)
expect_identical(x@sigmaMax, x@priorsBetas[[4]]@tauMax)
expect_identical(x@nuSigma, x@priorsBetas[[4]]@nuTau)
})
test_that("initialModel methods for Varying model can cope with orig-dest dimtypes", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 25, lambda = 20),
dim = c(5, 5),
dimnames = list(reg_orig = c("a", "b", "c", "d", "e"),
reg_dest = c("a", "b", "c", "d", "e"))))
y <- Counts(array(rbinom(n = 25, size = exposure, prob = 0.5),
dim = c(5, 5),
dimnames = list(reg_orig = c("a", "b", "c", "d", "e"),
reg_dest = c("a", "b", "c", "d", "e"))))
spec <- Model(y ~ Binomial(mean ~ reg_orig + reg_dest))
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
})
test_that("initialModel creates object of class CMPVaryingUseExp from valid inputs", {
initialModel <- demest:::initialModel
makeSigmaInitialPoisson <- demest:::makeSigmaInitialPoisson
loglm <- MASS:::loglm
makePriors <- demest:::makePriors
jitterBetas <- demest:::jitterBetas
## no missing values
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ CMP(mean ~ age + region))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
theta <- rgamma(n = 20,
shape = 0.05 * mean(y) + 0.95 * y,
rate = 0.05 * mean(exposure) + 0.95 * exposure)
sigma <- runif(1, max = 1)
theta <- array(theta, dim = dim(y), dimnames = dimnames(y))
betas <- loglm(~ age + region, theta)$param
theta <- as.numeric(theta)
logNuCMP <- stats::rnorm(n = length(theta))
nuCMP <- exp(logNuCMP)
nuCMP[nuCMP < 0.5] <- 0.5
nuCMP[nuCMP > 2] <- 2
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = NULL,
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priorsBetas = priors)
expect_is(x, "CMPVaryingUseExp")
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@ASigma, new("Scale", 1))
expect_identical(x@nuSigma, new("DegreesFreedom", 7))
expect_identical(x@sigmaMax, new("Scale", qhalft(0.999, df = 7, scale = 1)))
expect_identical(x@tolerance, 1e-5)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@nFailedPropTheta, new("Counter", 0L))
expect_identical(x@nuCMP, new("ParameterVector", nuCMP))
expect_identical(x@meanLogNuCMP, new("Parameter", 0))
expect_identical(x@sdLogNuCMP, new("Scale", 1))
## intercept only
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ CMP(mean ~ 1))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
theta <- rgamma(n = 20,
shape = 0.05 * mean(y) + 0.95 * y,
rate = 0.05 * mean(exposure) + 0.95 * exposure)
sigma <- runif(1, max = 1)
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_identical(sapply(x@priorsBetas, class), "ExchFixed")
## specify upper and lower
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ CMP(mean ~ age + region), lower = 0.01, upper = 2)
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(all(x@theta >= 0.01))
expect_true(all(x@theta <= 2))
expect_equal(x@lower, log(0.01))
expect_equal(x@upper, log(2))
## has missing values
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y[1:5] <- NA
spec <- Model(y ~ CMP(mean ~ age + region))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
mean.y <- mean(y[6:20])
mean.expose <- mean(exposure[6:20])
shape <- 0.05 * mean.y + 0.95 * y
shape[1:5] <- mean.y
rate <- 0.05 * mean.expose + 0.95 * exposure
rate[1:5] <- mean.expose
theta <- rgamma(n = 20, shape = shape, rate = rate)
sigma <- runif(1)
theta <- array(theta, dim = dim(y), dimnames = dimnames(y))
betas <- loglm(~ age + region, theta)$param
theta <- as.numeric(theta)
nuCMP <- rnorm(n = length(theta), mean = 0, sd = 1)
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = NULL,
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priorsBetas = priors)
expect_is(x, "CMPVaryingUseExp")
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@nuSigma, new("DegreesFreedom", 7))
expect_identical(x@ASigma, new("Scale", 1))
expect_identical(x@tolerance, 1e-5)
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@nFailedPropTheta, new("Counter", 0L))
## saturated model
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ CMP(mean ~ age * region),
age:region ~ Exch(error = Error(scale = HalfT(df = 3, scale = 0.1, max = 0.6))))
x <- initialModel(spec, y = y, exposure = exposure)
expect_identical(x@ASigma, x@priorsBetas[[4]]@ATau)
expect_identical(x@sigmaMax, x@priorsBetas[[4]]@tauMax)
expect_identical(x@nuSigma, x@priorsBetas[[4]]@nuTau)
## must have exposure
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ CMP(mean ~ age * region))
expect_error(initialModel(spec, y = y, exposure = NULL),
"model 'y ~ CMP\\(mean ~ age \\* region\\)' uses exposure, but no 'exposure' argument supplied")
})
test_that("initialModel creates object of class NormalVaryingVarsigmaKnown from valid inputs", {
initialModel <- demest:::initialModel
makeLinearBetas <- demest:::makeLinearBetas
makePriors <- demest:::makePriors
jitterBetas <- demest:::jitterBetas
## y has no missing values
varsigma <- 1.3
y <- Counts(array(rnorm(n = 20, mean = 0, sd = varsigma),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region, sd = varsigma))
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
set.seed(1)
sigma <- runif(n = 1, max = sd(y))
theta <- rnorm(n = 20, mean = 0.5 * mean(y) + 0.5 * y, sd = sigma)
theta <- array(theta, dim = dim(y), dimnames = dimnames(y))
betas <- makeLinearBetas(theta = theta, formula = mean ~ age + region)
theta <- as.numeric(theta)
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = sd(y),
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priors)
expect_is(x, "NormalVaryingVarsigmaKnown")
expect_identical(x@theta, theta)
expect_identical(x@varsigma, new("Scale", varsigma))
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@nuSigma, new("DegreesFreedom", 7))
expect_identical(x@ASigma, new("Scale", sd(y)))
expect_identical(x@sigmaMax@.Data, qhalft(p = 0.999, df = 7, scale = sd(y)))
expect_identical(x@lower, -Inf)
expect_identical(x@upper, Inf)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@tolerance, 1e-5)
## intercept only
varsigma <- 2.1
y <- Counts(array(rnorm(n = 20, mean = 0, sd = varsigma),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ 1, sd = varsigma))
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
set.seed(1)
sigma <- runif(1, max = sd(y))
theta <- rnorm(n = 20, mean = 0.5 * mean(y) + 0.5 * y, sd = sigma)
betas <- list(mean(theta))
betas <- jitterBetas(betas)
expect_identical(x@theta, theta)
expect_identical(x@betas, betas)
expect_identical(sapply(x@priorsBetas, class), "ExchFixed")
## lower and upper specified
varsigma <- 1.3
y <- Counts(array(rnorm(n = 20, mean = 0, sd = varsigma),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region, sd = varsigma), lower = 0, upper = 3)
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
expect_true(all(x@theta > 0))
expect_true(all(x@theta < 3))
expect_equal(x@lower, 0)
expect_equal(x@upper, 3)
## y has two values
varsigma <- 2
y <- Counts(array(rnorm(n = 20, mean = 0, sd = varsigma),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y[-(1:2)] <- NA
weights[-(1:2)] <- NA
spec <- Model(y ~ Normal(mean ~ age + region, sd = 2))
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
set.seed(1)
sigma <- runif(1, max = sd(y, na.rm = T))
y.bar <- mean(y, na.rm = T)
mean <- rep(y.bar, 20)
mean[1:2] <- 0.5 * y.bar + 0.5 * y[1:2]
theta <- rnorm(20, mean = mean, sd = sigma)
theta <- array(theta, dim = dim(y), dimnames = dimnames(y))
betas <- makeLinearBetas(theta = theta, formula = mean ~ age + region)
theta <- as.numeric(theta)
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = sd(y, na.rm = T),
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priors)
expect_is(x, "NormalVaryingVarsigmaKnown")
expect_identical(x@theta, theta)
expect_identical(x@varsigma, new("Scale", varsigma))
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@nuSigma, new("DegreesFreedom", 7))
expect_identical(x@ASigma, new("Scale", sd(y, na.rm = T)))
expect_identical(x@lower, -Inf)
expect_identical(x@upper, Inf)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@tolerance, 1e-5)
## maxSigma specified
varsigma <- 1.3
y <- Counts(array(rnorm(n = 20, mean = 0, sd = varsigma),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region, sd = varsigma),
priorSD = HalfT(max = 1000))
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
expect_identical(x@sigmaMax@.Data, 1000)
## saturated model
y <- Values(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age * region, sd = 1),
age:region ~ Exch(error = Error(scale = HalfT(df = 3, scale = 0.1, max = 0.6))))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
x <- initialModel(spec, y = y, weights = weights)
expect_identical(x@ASigma, x@priorsBetas[[4]]@ATau)
expect_identical(x@sigmaMax, x@priorsBetas[[4]]@tauMax)
expect_identical(x@nuSigma, x@priorsBetas[[4]]@nuTau)
## varsigma is 0
y <- Values(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age * region, sd = 0),
age:region ~ Exch(error = Error(scale = HalfT(df = 3, scale = 0.1, max = 0.6))))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
x <- initialModel(spec, y = y, weights = weights)
expect_identical(x@varsigma@.Data, 0)
expect_identical(x@theta, as.numeric(y))
})
test_that("initialModel creates object of class NormalVaryingVarsigmaUnknown from valid inputs", {
initialModel <- demest:::initialModel
rinvchisq1 <- demest:::rinvchisq1
makeLinearBetas <- demest:::makeLinearBetas
makePriors <- demest:::makePriors
jitterBetas <- demest:::jitterBetas
## no missing values
y <- Counts(array(rnorm(n = 20, mean = 0),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region),
priorSD = HalfT(max = 20))
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
set.seed(1)
I <- length(y)
varsigma <- runif(1, max = sd(y))
sigma <- runif(1, max = sd(y))
theta <- rnorm(n = 20, mean = 0.5 * mean(y) + 0.5 * y, sd = sigma)
theta <- array(theta, dim = dim(y), dimnames = dimnames(y))
betas <- makeLinearBetas(theta = theta, formula = mean ~ age + region)
theta <- as.numeric(theta)
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = sd(y),
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priors)
expect_is(x, "NormalVaryingVarsigmaUnknown")
expect_identical(x@theta, theta)
expect_identical(x@varsigma, new("Scale", varsigma))
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@sigmaMax, new("Scale", 20))
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@lower, -Inf)
expect_identical(x@upper, Inf)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@tolerance, 1e-5)
## intercept only
y <- Counts(array(rnorm(n = 20, mean = 0, sd = varsigma),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ 1))
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
set.seed(1)
I <- length(y)
varsigma <- runif(1, max = sd(y))
sigma <- runif(1, max = sd(y))
theta <- rnorm(n = 20, mean = 0.5 * mean(y) + 0.5 * y, sd = sigma)
betas <- list(mean(theta))
betas <- jitterBetas(betas)
expect_identical(x@theta, theta)
expect_identical(x@betas, betas)
expect_identical(sapply(x@priorsBetas, class), "ExchFixed")
## lower and upper specified
y <- Counts(array(rnorm(n = 20, mean = 0, sd = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region), lower = 0, upper = 3)
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
expect_true(all(x@theta > 0))
expect_true(all(x@theta < 3))
expect_equal(x@lower, 0)
expect_equal(x@upper, 3)
## has missing values
y <- Counts(array(rnorm(n = 20, mean = 0),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y[1:10] <- NA
weights[1:5] <- NA
spec <- Model(y ~ Normal(mean ~ age + region))
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
set.seed(1)
varsigma <- runif(n = 1, max = sd(y, na.rm = T))
sigma <- runif(n = 1, max = sd(y, na.rm = T))
mean <- 0.5 * mean(y[11:20]) + 0.5 * y
mean[1:10] <- mean(y[11:20])
theta <- rnorm(n = 20, mean = mean, sd = sigma)
theta <- array(theta, dim = dim(y), dimnames = dimnames(y))
betas <- makeLinearBetas(theta = theta, formula = mean ~ age + region)
theta <- as.numeric(theta)
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = sd(y, na.rm = T),
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priors)
expect_is(x, "NormalVaryingVarsigmaUnknown")
expect_identical(x@theta, theta)
expect_identical(x@varsigma, new("Scale", varsigma))
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@lower, -Inf)
expect_identical(x@upper, Inf)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@tolerance, 1e-5)
## saturated model
y <- Values(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age * region),
age:region ~ Exch(error = Error(scale = HalfT(df = 3, scale = 0.1, max = 0.6))))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
x <- initialModel(spec, y = y, weights = weights)
expect_identical(x@ASigma, x@priorsBetas[[4]]@ATau)
expect_identical(x@sigmaMax, x@priorsBetas[[4]]@tauMax)
expect_identical(x@nuSigma, x@priorsBetas[[4]]@nuTau)
})
test_that("initialModel creates object of class PoissonVaryingUseExp from valid inputs", {
initialModel <- demest:::initialModel
makeSigmaInitialPoisson <- demest:::makeSigmaInitialPoisson
loglm <- MASS:::loglm
makePriors <- demest:::makePriors
jitterBetas <- demest:::jitterBetas
## no missing values
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
theta <- rgamma(n = 20,
y + 1,
rate = exposure + 1)
sigma <- runif(1, max = 1)
theta <- array(theta, dim = dim(y), dimnames = dimnames(y))
betas <- loglm(~ age + region, theta)$param
theta <- as.numeric(theta)
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = NULL,
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priors)
expect_is(x, "PoissonVaryingUseExp")
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@ASigma, new("Scale", 1))
expect_identical(x@nuSigma, new("DegreesFreedom", 7))
expect_identical(x@sigmaMax, new("Scale", qhalft(0.999, df = 7, scale = 1)))
expect_identical(x@tolerance, 1e-5)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@nFailedPropTheta, new("Counter", 0L))
## intercept only
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ 1))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
theta <- rgamma(n = 20,
shape = y + 1,
rate = exposure + 1)
sigma <- runif(1, max = 1)
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_identical(sapply(x@priorsBetas, class), "ExchFixed")
## specify upper and lower
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region), lower = 0.01, upper = 2)
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(all(x@theta >= 0.01))
expect_true(all(x@theta <= 2))
expect_equal(x@lower, log(0.01))
expect_equal(x@upper, log(2))
## has missing values
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y[1:5] <- NA
spec <- Model(y ~ Poisson(mean ~ age + region))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
mean.y <- mean(y[6:20])
mean.expose <- mean(exposure[6:20])
shape <- ifelse(is.na(y), 1, y + 1)
rate <- ifelse(is.na(y), 1, exposure + 1)
theta <- rgamma(n = 20, shape = shape, rate = rate)
sigma <- runif(1)
theta <- array(theta, dim = dim(y), dimnames = dimnames(y))
betas <- loglm(~ age + region, theta)$param
theta <- as.numeric(theta)
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = NULL,
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priors)
expect_is(x, "PoissonVaryingUseExp")
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@nuSigma, new("DegreesFreedom", 7))
expect_identical(x@ASigma, new("Scale", 1))
expect_identical(x@tolerance, 1e-5)
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@nFailedPropTheta, new("Counter", 0L))
## saturated model
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age * region),
age:region ~ Exch(error = Error(scale = HalfT(df = 3, scale = 0.1, max = 0.6))))
x <- initialModel(spec, y = y, exposure = exposure)
expect_identical(x@ASigma, x@priorsBetas[[4]]@ATau)
expect_identical(x@sigmaMax, x@priorsBetas[[4]]@tauMax)
expect_identical(x@nuSigma, x@priorsBetas[[4]]@nuTau)
## must have exposure
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age * region))
expect_error(initialModel(spec, y = y, exposure = NULL),
"model 'y ~ Poisson\\(mean ~ age \\* region\\)' uses exposure, but no 'exposure' argument supplied")
})
test_that("initialModel creates object of class PoissonVaryingNotUseExp from valid inputs", {
initialModel <- demest:::initialModel
makeSigmaInitialPoisson <- demest:::makeSigmaInitialPoisson
## no missing values
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE),
priorSD = HalfT(df = 6, max = 10))
set.seed(1)
x <- initialModel(spec, y = y, exposure = NULL)
set.seed(1)
theta <- rgamma(n = 20,
shape = y + 1,
rate = 2)
sigma <- runif(1, max = sd(log(y+1)))
expect_is(x, "PoissonVaryingNotUseExp")
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@nuSigma, new("DegreesFreedom", 6))
expect_identical(x@ASigma, new("Scale", sd(log(y+1))))
expect_identical(x@sigmaMax, new("Scale", 10))
expect_identical(x@tolerance, 1e-5)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@nFailedPropTheta, new("Counter", 0L))
## intercept only
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ 1, useExpose = FALSE))
set.seed(1)
x <- initialModel(spec, y = y, exposure = NULL)
set.seed(1)
theta <- rgamma(n = 20,
shape = y + 1,
rate = 2)
sigma <- runif(1, max = sd(log(y+1)))
m <- mean(y)
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_identical(sapply(x@priorsBetas, class), "ExchFixed")
## has missing values
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y[1:5] <- NA
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE))
set.seed(1)
x <- initialModel(spec, y = y, exposure = NULL)
set.seed(1)
mean <- mean(y, na.rm = TRUE)
shape <- ifelse(is.na(y), 1, y + 1)
rate <- 2
theta <- rgamma(n = 20, shape = shape, rate = rate)
sigma <- runif(1, max = sd(log(y+1), na.rm = T))
expect_is(x, "PoissonVaryingNotUseExp")
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@nuSigma, new("DegreesFreedom", 7))
expect_identical(x@ASigma, new("Scale", sd(log(y+1), na.rm = TRUE)))
expect_identical(x@tolerance, 1e-5)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@nFailedPropTheta, new("Counter", 0L))
## saturated model
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age * region, useExpose = FALSE),
age:region ~ Exch(error = Error(scale = HalfT(df = 3, scale = 0.1, max = 0.6))))
x <- initialModel(spec, y = y, exposure = NULL)
expect_identical(x@ASigma, x@priorsBetas[[4]]@ATau)
expect_identical(x@sigmaMax, x@priorsBetas[[4]]@tauMax)
expect_identical(x@nuSigma, x@priorsBetas[[4]]@nuTau)
## must not have exposure
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age * region, useExpose = FALSE))
expect_error(initialModel(spec, y = y, exposure = exposure),
"exposure' argument supplied, but model 'y ~ Poisson\\(mean ~ age \\* region, useExpose = FALSE\\)' does not use exposure")
})
## PoissonBinomialMixture #############################################################
test_that("initialModel creates object of class PoissonBinomialMixture from valid inputs", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ PoissonBinomial(prob = 0.98))
x <- initialModel(spec, y = y, exposure = exposure)
expect_equal(x,
new("PoissonBinomialMixture",
call = call("Model", formula = y ~ PoissonBinomial(prob = 0.98)),
prob = 0.98,
metadataY = y@metadata))
})
## Round3 #############################################################
test_that("initialModel creates object of class Round3 from valid inputs", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- round3(y)
spec <- Model(y ~ Round3())
x <- initialModel(spec, y = y, exposure = exposure)
expect_equal(x,
new("Round3",
call = call("Model", formula = y ~ Round3()),
metadataY = y@metadata))
expect_error(initialModel(spec, y = y + 1L, exposure = exposure),
"using 'Round3' data model, but data contains values not divisible by 3")
})
## Exact #############################################################
test_that("initialModel creates object of class Exact from valid inputs", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Exact())
x <- initialModel(spec, y = y, exposure = exposure)
expect_equal(x,
new("Exact",
call = call("Model", formula = y ~ Exact()),
metadataY = y@metadata))
y[1] <- NA
expect_error(initialModel(spec, y = y, exposure = exposure),
"using 'Exact' data model, but data contains missing values")
})
## NormalFixed #############################################################
test_that("initialModel creates object of class NormalFixedUseExp from valid inputs", {
initialModel <- demest:::initialModel
mean <- Values(array(runif(n = 25),
dim = c(5, 5),
dimnames = list(age = 0:4, region = letters[1:5])))
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ NormalFixed(mean = mean, sd = 0.1))
ans.obtained <- initialModel(spec, y = y, exposure = exposure)
ans.expected <- new("NormalFixedUseExp",
call = call("Model", formula = y ~ NormalFixed(mean = mean, sd = 0.1)),
mean = new("ParameterVector", mean@.Data[1:20]),
sd = new("ScaleVec", rep(0.1, 20)),
metadataY = y@metadata,
meanAll = new("ParameterVector", mean@.Data),
sdAll = new("ScaleVec", rep(0.1, 25)),
metadataAll = mean@metadata)
expect_equal(ans.obtained, ans.expected)
})
test_that("initialModel throws appropriate errors with NormalFixedUseExp", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
mean.wrong <- Values(array(runif(n = 25),
dim = c(5, 5),
dimnames = list(age = 0:4, region = letters[2:6])))
spec <- Model(y ~ NormalFixed(mean = mean.wrong, sd = 0.1))
expect_error(initialModel(spec, y = y, exposure = exposure),
"'mean' from NormalFixed model not compatible with data :")
mean <- Values(array(runif(n = 25),
dim = c(5, 5),
dimnames = list(age = 0:4, region = letters[1:5])))
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ NormalFixed(mean = mean, sd = 0.1))
expect_error(initialModel(spec, y = y, exposure = NULL),
"model 'y ~ NormalFixed\\(mean = mean, sd = 0\\.1\\)' uses exposure, but no 'exposure' argument supplied")
})
test_that("initialModel creates object of class NormalFixedNotUseExp from valid inputs", {
initialModel <- demest:::initialModel
mean <- Values(array(runif(n = 25),
dim = c(5, 5),
dimnames = list(age = 0:4, region = letters[1:5])))
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ NormalFixed(mean = mean, sd = 0.1, useExpose = FALSE))
ans.obtained <- initialModel(spec, y = y, exposure = NULL)
ans.expected <- new("NormalFixedNotUseExp",
call = call("Model", formula = y ~ NormalFixed(mean = mean, sd = 0.1, useExpose = FALSE)),
mean = new("ParameterVector", mean@.Data[1:20]),
sd = new("ScaleVec", rep(0.1, 20)),
metadataY = y@metadata,
meanAll = new("ParameterVector", mean@.Data),
sdAll = new("ScaleVec", rep(0.1, 25)),
metadataAll = mean@metadata)
expect_equal(ans.obtained, ans.expected)
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ NormalFixed(mean = mean, sd = 0.1, useExpose = FALSE))
expect_error(initialModel(spec, y = y, exposure = exposure),
"'exposure' argument supplied, but model 'y ~ NormalFixed\\(mean = mean, sd = 0.1, useExpose = FALSE\\)' does not use exposure")
})
## TFixed #############################################################
test_that("initialModel creates object of class TFixedUseExp from valid inputs", {
initialModel <- demest:::initialModel
location <- Values(array(runif(n = 25),
dim = c(5, 5),
dimnames = list(age = 0:4, region = letters[1:5])))
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ TFixed(location = location, scale = 0.1, df = 8))
ans.obtained <- initialModel(spec, y = y, exposure = exposure)
ans.expected <- new("TFixedUseExp",
call = call("Model", formula = y ~ TFixed(location = location, scale = 0.1, df = 8)),
mean = new("ParameterVector", location@.Data[1:20]),
sd = new("ScaleVec", rep(0.1, 20)),
metadataY = y@metadata,
meanAll = new("ParameterVector", location@.Data),
sdAll = new("ScaleVec", rep(0.1, 25)),
nu = new("DegreesFreedom", 8),
metadataAll = location@metadata)
expect_equal(ans.obtained, ans.expected)
})
test_that("initialModel throws appropriate errors with TFixedUseExp", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
location.wrong <- Values(array(runif(n = 25),
dim = c(5, 5),
dimnames = list(age = 0:4, region = letters[2:6])))
spec <- Model(y ~ TFixed(location = location.wrong, scale = 0.1))
expect_error(initialModel(spec, y = y, exposure = exposure),
"'location' from TFixed model not compatible with data :")
location <- Values(array(runif(n = 25),
dim = c(5, 5),
dimnames = list(age = 0:4, region = letters[1:5])))
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ TFixed(location = location, scale = 0.1))
expect_error(initialModel(spec, y = y, exposure = NULL),
"model 'y ~ TFixed\\(location = location, scale = 0\\.1\\)' uses exposure, but no 'exposure' argument supplied")
})
test_that("initialModel creates object of class TFixedNotUseExp from valid inputs", {
initialModel <- demest:::initialModel
location <- Values(array(runif(n = 25),
dim = c(5, 5),
dimnames = list(age = 0:4, region = letters[1:5])))
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ TFixed(location = location, scale = 0.1, useExpose = FALSE))
ans.obtained <- initialModel(spec, y = y, exposure = NULL)
ans.expected <- new("TFixedNotUseExp",
call = call("Model", formula = y ~ TFixed(location = location, scale = 0.1, useExpose = FALSE)),
mean = new("ParameterVector", location@.Data[1:20]),
sd = new("ScaleVec", rep(0.1, 20)),
metadataY = y@metadata,
meanAll = new("ParameterVector", location@.Data),
sdAll = new("ScaleVec", rep(0.1, 25)),
nu = new("DegreesFreedom", 7),
metadataAll = location@metadata)
expect_equal(ans.obtained, ans.expected)
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ TFixed(location = location, scale = 0.1, useExpose = FALSE))
expect_error(initialModel(spec, y = y, exposure = exposure),
"'exposure' argument supplied, but model 'y ~ TFixed\\(location = location, scale = 0.1, useExpose = FALSE\\)' does not use exposure")
})
## LN2 ###############################################################################
test_that("initialModel creates object of class LN2 from valid inputs", {
initialModel <- demest:::initialModel
constraint <- Values(array(c(NA, -1L, 0L, 1L),
dim = c(2, 2),
dimnames = list(age = c("0-39", "40+"),
sex = c("Female", "Male"))))
y <- Counts(array(10L,
dim = c(2, 4, 3),
dimnames = c(list(sex = c("Female", "Male"),
age = c("0-19", "20-39", "40-59", "60+"),
time = c("2000", "2010", "2020")))))
exposure <- 2L * y
spec <- Model(y ~ LN2(constraint = constraint))
ans.obtained <- initialModel(spec,
y = y,
exposure = exposure)
expect_true(validObject(ans.obtained))
expect_true(ans.obtained@add1@.Data)
y[2,1:2,] <- 0L
sz <- Values(array(c(1L, 1L, 0L, 1L),
dim = c(2, 2),
dimnames = list(age = c("0-39", "40+"),
sex = c("Female", "Male"))))
concordances <- list(sex = Concordance(data.frame(from = c("F", "M", "Female", "Male"),
to = c("Female", "Male", "Female", "Male"))))
spec <- Model(y ~ LN2(constraint = constraint,
structuralZeros = sz,
concordances = concordances))
ans.obtained <- initialModel(spec,
y = y,
exposure = exposure)
expect_true(validObject(ans.obtained))
y[2,1:2,] <- 0L
sz <- Values(array(c(1L, 1L, 0L, 1L),
dim = c(2, 2),
dimnames = list(age = c("0-39", "40+"),
sex = c("Female", "Male"))))
concordances <- list(sex = Concordance(data.frame(from = c("F", "M", "Female", "Male"),
to = c("Female", "Male", "Female", "Male"))))
spec <- Model(y ~ LN2(constraint = constraint,
structuralZeros = sz,
concordances = concordances,
sd = 0.2))
ans.obtained <- initialModel(spec,
y = y,
exposure = exposure)
expect_true(validObject(ans.obtained))
expect_identical(ans.obtained@varsigma@.Data, 0.2)
sz <- Values(array(c(1L, 1L, 0L, 1L),
dim = c(2, 2),
dimnames = list(age = c("0-39", "40+"),
sex = c("Female", "Male"))))
spec <- Model(y ~ LN2(constraint = constraint,
structuralZeros = sz,
sd = InvChiSq(df = 5, scaleSq = 10)))
ans.obtained <- initialModel(spec,
y = y,
exposure = exposure)
expect_true(validObject(ans.obtained))
expect_false(ans.obtained@varsigmaLN2HasHalfT@.Data)
constraint <- Values(array(c(NA, -1L, 0L, 1L),
dim = c(2, 2),
dimnames = list(age = c("0-39", "40+"),
sex = c("Female", "Male"))))
y <- Counts(array(10L,
dim = c(2, 4, 3),
dimnames = c(list(sex = c("Female", "Male"),
age = c("0-19", "20-39", "40-59", "60+"),
time = c("2000", "2010", "2020")))))
exposure <- 2L * y
spec <- Model(y ~ LN2(constraint = constraint, add1 = FALSE))
ans.obtained <- initialModel(spec,
y = y,
exposure = exposure)
expect_true(validObject(ans.obtained))
expect_false(ans.obtained@add1@.Data)
})
## Aggregate #########################################################################
test_that("initialModel works with BinomialVarying and AgCertain", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])),
dimscales = c(age = "Intervals"))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])),
dimscales = c(age = "Intervals"))
## scalar
aggregate <- AgCertain(value = 0.4)
spec <- Model(y ~ Binomial(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "BinomialVaryingAgCertain")
## Values
value <- Values(array(c(0.2, 0.3, 0.4),
dim = 3,
dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgCertain(value = value)
spec <- Model(y ~ Binomial(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "BinomialVaryingAgCertain")
})
test_that("initialModel works with BinomialVarying and AgNormal", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgNormal(value = 0.4, sd = 0.1)
spec <- Model(y ~ Binomial(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "BinomialVaryingAgNormal")
## Values
value <- Values(array(c(0.2, 0.3, 0.4), dim = 3, dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgNormal(value = value, sd = sqrt(value))
spec <- Model(y ~ Binomial(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "BinomialVaryingAgNormal")
})
test_that("initialModel works with NormalVaryingVarsigmaKnown and AgCertain", {
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20, mean = 0, sd = 0.3),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgCertain(value = 0.4)
spec <- Model(y ~ Normal(mean ~ age + region, sd = 0.2),
aggregate = aggregate)
x <- initialModel(spec, y = y, weights = weights)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaKnownAgCertain")
## Values
value <- Values(array(c(0.2, 0.3, 0.4),
dim = 3,
dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgCertain(value = value)
spec <- Model(y ~ Normal(mean ~ age + region, sd = 0.2),
aggregate = aggregate)
x <- initialModel(spec, y = y, weights = weights)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaKnownAgCertain")
})
test_that("initialModel works with NormalVaryingVarsigmaKnown and AgNormal", {
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20, mean = 0, sd = 0.3),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgNormal(value = 0.4, sd = 0.1)
spec <- Model(y ~ Normal(mean ~ age + region, sd = 0.2),
aggregate = aggregate)
x <- initialModel(spec, y = y, weights = weights)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaKnownAgNormal")
## Values
value <- Values(array(c(0.2, 0.3, 0.4), dim = 3, dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgNormal(value = value, sd = sqrt(value))
spec <- Model(y ~ Normal(mean ~ age + region, sd = 0.2),
aggregate = aggregate)
x <- initialModel(spec, y = y, weights = weights)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaKnownAgNormal")
})
test_that("initialModel works with NormalVaryingVarsigmaUnknown and AgCertain", {
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20, mean = 0, sd = 0.3),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgCertain(value = 0.4)
spec <- Model(y ~ Normal(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, weights = weights)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaUnknownAgCertain")
## Values
value <- Values(array(c(0.2, 0.3, 0.4),
dim = 3,
dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgCertain(value = value)
spec <- Model(y ~ Normal(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, weights = weights)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaUnknownAgCertain")
})
test_that("initialModel works with NormalVaryingVarsigmaUnknown and AgNormal", {
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20, mean = 0, sd = 0.3),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgNormal(value = 0.4, sd = 0.1)
spec <- Model(y ~ Normal(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, weights = weights)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaUnknownAgNormal")
## Values
value <- Values(array(c(0.2, 0.3, 0.4), dim = 3, dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgNormal(value = value, sd = sqrt(value))
spec <- Model(y ~ Normal(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, weights = weights)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaUnknownAgNormal")
})
test_that("initialModel works with PoissonVaryingNotUseExp and AgCertain", {
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgCertain(value = 0.4)
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE),
upper = 3,
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = NULL)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingNotUseExpAgCertain")
## Values
value <- Values(array(c(0.2, 0.3, 0.4),
dim = 3,
dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgCertain(value = value)
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = NULL)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingNotUseExpAgCertain")
})
test_that("initialModel works with PoissonVaryingNotUseExp and AgNormal", {
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgNormal(value = 0.4, sd = 0.1)
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = NULL)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingNotUseExpAgNormal")
## Values
value <- Values(array(c(0.2, 0.3, 0.4), dim = 3, dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgNormal(value = value, sd = sqrt(value))
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = NULL)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingNotUseExpAgNormal")
})
test_that("initialModel works with PoissonVaryingUseExp and AgCertain", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.5 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgCertain(value = 0.4)
spec <- Model(y ~ Poisson(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingUseExpAgCertain")
## Values
value <- Values(array(c(0.2, 0.3, 0.4),
dim = 3,
dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgCertain(value = value)
spec <- Model(y ~ Poisson(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingUseExpAgCertain")
})
test_that("initialModel works with PoissonVaryingUseExp and AgNormal", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.5 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgNormal(value = 0.4, sd = 0.1)
spec <- Model(y ~ Poisson(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingUseExpAgNormal")
## Values
value <- Values(array(c(0.2, 0.3, 0.4), dim = 3, dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgNormal(value = value, sd = sqrt(value))
spec <- Model(y ~ Poisson(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingUseExpAgNormal")
})
test_that("initialModel works with PoissonVaryingNotUseExp and AgPoisson", {
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgPoisson(value = 0.4)
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = NULL)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingNotUseExpAgPoisson")
## Values
value <- Values(array(c(0.2, 0.3, 0.4), dim = 3, dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgPoisson(value = value)
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = NULL)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingNotUseExpAgPoisson")
})
test_that("initialModel works with PoissonVaryingUseExp and AgPoisson", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.5 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgPoisson(value = 0.4)
spec <- Model(y ~ Poisson(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingUseExpAgPoisson")
## Values
value <- Values(array(c(0.2, 0.3, 0.4), dim = 3, dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgPoisson(value = value)
spec <- Model(y ~ Poisson(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingUseExpAgPoisson")
})
test_that("initialModel works with PoissonVaryingUseExp and AgLife", {
initialModel <- demest:::initialModel
expose <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = c(0:3, "4+"), region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.5 * expose),
dim = c(5, 4),
dimnames = list(age = c(0:3, "4+"), region = letters[1:4])))
## scalar
aggregate <- AgLife(value = 20, sd = 0.1)
spec <- Model(y ~ Poisson(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = expose)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingUseExpAgLife")
## Values
value <- Values(array(c(0.2, 0.3, 0.4), dim = 3, dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgNormal(value = value, sd = sqrt(value))
spec <- Model(y ~ Poisson(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = expose)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingUseExpAgNormal")
})
## initialModelPredict ###############################################################
test_that("initialModelPredict works with BinomialVarying - without aggregate", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(as.integer(runif(n = 20, min = 5, max = 100)),
dim = c(5, 4),
dimnames = list(time = 2001:2005, region = 1:4)),
dimscales = c(time = "Intervals"))
weights.bench <- Counts(array(as.integer(runif(n = 20, min = 5, max = 100)),
dim = c(4, 4),
dimnames = list(time = 2006:2009, region = 1:4)),
dimscales = c(time = "Intervals"))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.7),
dim = c(5, 4),
dimnames = list(time = 2001:2005, region = 1:4)),
dimscales = c(time = "Intervals"))
spec <- Model(y ~ Binomial(mean ~ time + region))
mod.est <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(mod.est,
along = 1L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = NULL,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "BinomialVaryingPredict")
})
test_that("initialModelPredict works with BinomialVaryingPredict - with aggregate", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(as.integer(runif(n = 20, min = 5, max = 100)),
dim = c(5, 4),
dimnames = list(time = 2001:2005, region = 1:4)),
dimscales = c(time = "Intervals"))
weights.bench <- Counts(array(as.integer(runif(n = 20, min = 5, max = 100)),
dim = c(4, 4),
dimnames = list(time = 2006:2009, region = 1:4)),
dimscales = c(time = "Intervals"))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.7),
dim = c(5, 4),
dimnames = list(time = 2001:2005, region = 1:4)),
dimscales = c(time = "Intervals"))
## AgCertain with weights
aggregate <- AgCertain(value = 0.5, weights = weights.bench)
spec <- Model(y ~ Binomial(mean ~ time + region))
mod.est <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(mod.est,
along = 1L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = aggregate,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "BinomialVaryingPredictAgCertain")
## AgNormal with weights
aggregate <- AgNormal(value = 0.5, sd = 1, weights = weights.bench)
spec <- Model(y ~ Binomial(mean ~ time + region))
mod.est <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(mod.est,
along = 1L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = aggregate,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "BinomialVaryingPredictAgNormal")
})
test_that("initialModelPredict works with NormalVaryingVarsigmaUnknown - without aggregate", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20, mean = 0, sd = 0.3),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region))
x <- initialModel(spec, y = y, weights = weights)
set.seed(1)
x <- initialModelPredict(x,
along = 2L,
labels = c("a", "b", "c", "d"),
n = NULL,
offsetModel = 1L,
covariates = NULL,
aggregate = NULL,
lower = 1,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaUnknownPredict")
})
test_that("initialModelPredict works with NormalVaryingVarsigmaUnknownPredict - with aggregate", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(20),
dim = c(5, 4),
dimnames = list(time = 2001:2005, region = 1:4)),
dimscales = c(time = "Intervals"))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(time = 2001:2005, region = 1:4)),
dimscales = c(time = "Intervals"))
weights.ag <- Counts(array(1,
dim = c(4, 4),
dimnames = list(time = 2006:2009, region = 1:4)),
dimscales = c(time = "Intervals"))
## AgCertain
aggregate <- AgCertain(value = 0.5, weights = weights.ag)
spec <- Model(y ~ Normal(mean ~ time + region))
mod.est <- initialModel(spec, y = y, weights = weights)
x <- initialModelPredict(mod.est,
along = 1L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = aggregate,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaUnknownPredictAgCertain")
## AgNormal
aggregate <- AgNormal(value = 0.5, sd = 1, weights = weights.ag)
spec <- Model(y ~ Normal(mean ~ time + region))
mod.est <- initialModel(spec, y = y, weights = weights)
x <- initialModelPredict(mod.est,
along = 1L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = aggregate,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaUnknownPredictAgNormal")
})
test_that("initialModelPredict works with PoissonVaryingUseExp", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.5 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region))
x <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(x,
along = 2L,
labels = c("a", "b", "c", "d"),
n = NULL,
offsetModel = 1L,
covariates = NULL,
aggregate = NULL,
lower = NULL,
upper = 100000)
expect_is(x, "PoissonVaryingUseExpPredict")
})
test_that("initialModelPredict works with PoissonVaryingNotUseExpPredict - with aggregate", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
y <- Counts(array(rpois(20, lambda = 100),
dim = c(5, 4),
dimnames = list(time = 2001:2005, region = 1:4)),
dimscales = c(time = "Intervals"))
weights.bench <- Counts(array(1,
dim = c(4, 4),
dimnames = list(time = 2006:2009, region = 1:4)),
dimscales = c(time = "Intervals"))
## AgCertain with no weights
aggregate <- AgCertain(value = 0.5, weights = weights.bench)
spec <- Model(y ~ Poisson(mean ~ time + region, useExpose = FALSE))
mod.est <- initialModel(spec, y = y, exposure = NULL)
x <- initialModelPredict(mod.est,
along = 1L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = aggregate,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingNotUseExpPredictAgCertain")
expect_true(all(x@weightAg == 1))
## AgNormal
aggregate <- AgNormal(value = 0.5, sd = 1, weights = weights.bench)
spec <- Model(y ~ Poisson(mean ~ time + region, useExpose = FALSE))
mod.est <- initialModel(spec, y = y, exposure = NULL)
x <- initialModelPredict(mod.est,
along = 1L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = aggregate,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingNotUseExpPredictAgNormal")
})
test_that("initialModelPredict works with CMPVaryingUseExp", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.5 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ CMP(mean ~ age + region))
x <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(x,
along = 2L,
labels = c("a", "b", "c", "d"),
n = NULL,
offsetModel = 1L,
covariates = NULL,
aggregate = NULL,
lower = NULL,
upper = 100000)
expect_is(x, "CMPVaryingUseExpPredict")
})
test_that("initialModelPredict works with CMPVaryingUseExp", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.5 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ CMP(mean ~ age + region))
x <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(x,
along = 2L,
labels = c("a", "b", "c", "d"),
n = NULL,
offsetModel = 1L,
covariates = NULL,
aggregate = NULL,
lower = NULL,
upper = 100000)
expect_is(x, "CMPVaryingUseExpPredict")
})
test_that("initialModelPredict works with PoissonBinomial", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ PoissonBinomial(prob = 0.95))
model <- initialModel(spec, y = y, exposure = exposure)
ans <- initialModelPredict(model,
along = 1L,
labels = NULL,
n = 5L,
offsetModel = 1L)
expect_true(validObject(ans))
metadata.expected <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 5:9, region = letters[1:4])))@metadata
expect_identical(ans@metadataY, metadata.expected)
expect_identical(ans@iMethodModel, model@iMethodModel + 100L)
})
test_that("initialModelPredict works with Round3", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- round3(y)
spec <- Model(y ~ Round3())
model <- initialModel(spec, y = y, exposure = exposure)
ans <- initialModelPredict(model,
along = 1L,
labels = NULL,
n = 5L,
offsetModel = 1L)
expect_true(validObject(ans))
metadata.expected <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 5:9, region = letters[1:4])))@metadata
expect_identical(ans@metadataY, metadata.expected)
expect_identical(ans@iMethodModel, model@iMethodModel + 100L)
})
test_that("initialModelPredict works with Exact", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Exact())
model <- initialModel(spec, y = y, exposure = exposure)
ans <- initialModelPredict(model,
along = 1L,
labels = NULL,
n = 5L,
offsetModel = 1L)
expect_true(validObject(ans))
metadata.expected <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 5:9, region = letters[1:4])))@metadata
expect_identical(ans@metadataY, metadata.expected)
expect_identical(ans@iMethodModel, model@iMethodModel + 100L)
})
test_that("initialModelPredict works with NormFixedUseExp", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
mean <- Values(array(rnorm(n = 40),
dim = c(10, 4),
dimnames = list(age = 0:9, region = letters[1:4])))
spec <- Model(y ~ NormalFixed(mean = mean, sd = sqrt(abs(mean))))
model <- initialModel(spec, y = y, exposure = exposure)
ans <- initialModelPredict(model,
along = 1L,
labels = NULL,
n = 5L,
offsetModel = 1L)
expect_true(validObject(ans))
metadata.expected <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 5:9, region = letters[1:4])))@metadata
expect_identical(ans@metadataY, metadata.expected)
expect_identical(ans@iMethodModel, model@iMethodModel + 100L)
})
test_that("initialModelPredict works with NormFixedNotUseExp", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
mean <- Values(array(rnorm(n = 40),
dim = c(10, 4),
dimnames = list(age = 0:9, region = letters[1:4])))
spec <- Model(y ~ NormalFixed(mean = mean, sd = sqrt(abs(mean)), useExpose = FALSE))
model <- initialModel(spec, y = y, exposure = NULL)
ans <- initialModelPredict(model,
along = 1L,
labels = NULL,
n = 5L,
offsetModel = 1L)
expect_true(validObject(ans))
metadata.expected <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 5:9, region = letters[1:4])))@metadata
expect_identical(ans@metadataY, metadata.expected)
expect_identical(ans@iMethodModel, model@iMethodModel + 100L)
})
test_that("initialModelPredict works with NormFixedUseExp", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
location <- Values(array(rnorm(n = 40),
dim = c(10, 4),
dimnames = list(age = 0:9, region = letters[1:4])))
spec <- Model(y ~ TFixed(location = location, scale = sqrt(abs(location))))
model <- initialModel(spec, y = y, exposure = exposure)
ans <- initialModelPredict(model,
along = 1L,
labels = NULL,
n = 5L,
offsetModel = 1L)
expect_true(validObject(ans))
metadata.expected <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 5:9, region = letters[1:4])))@metadata
expect_identical(ans@metadataY, metadata.expected)
expect_identical(ans@iMethodModel, model@iMethodModel + 100L)
})
test_that("initialModelPredict works with NormFixedNotUseExp", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
location <- Values(array(rnorm(n = 40),
dim = c(10, 4),
dimnames = list(age = 0:9, region = letters[1:4])))
spec <- Model(y ~ TFixed(location = location, scale = sqrt(abs(location)), useExpose = FALSE))
model <- initialModel(spec, y = y, exposure = NULL)
ans <- initialModelPredict(model,
along = 1L,
labels = NULL,
n = 5L,
offsetModel = 1L)
expect_true(validObject(ans))
metadata.expected <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 5:9, region = letters[1:4])))@metadata
expect_identical(ans@metadataY, metadata.expected)
expect_identical(ans@iMethodModel, model@iMethodModel + 100L)
})
test_that("initialModelPredict creates object of class LN2Predict from valid inputs", {
initialModel <- demest:::initialModel
initialModelPredict <- demest:::initialModelPredict
## 'constraint' does not have time dimension
constraint <- Values(array(c(NA, -1L, 0L, 1L),
dim = c(2, 2),
dimnames = list(age = c("0-39", "40+"),
sex = c("Female", "Male"))))
y <- Counts(array(10L,
dim = c(2, 4, 3),
dimnames = c(list(sex = c("Female", "Male"),
age = c("0-19", "20-39", "40-59", "60+"),
time = c("2000", "2010", "2020")))))
exposure <- 2L * y
spec <- Model(y ~ LN2(constraint = constraint))
mod.est <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(mod.est,
along = 3L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = NULL,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "LN2Predict")
## 'constraint' has time dimension
constraint <- Values(array(c(NA, -1L),
dim = c(2, 7),
dimnames = list(age = c("0-39", "40+"),
time = seq(2000, 2060, 10))))
y <- Counts(array(10L,
dim = c(2, 4, 3),
dimnames = c(list(sex = c("Female", "Male"),
age = c("0-19", "20-39", "40-59", "60+"),
time = c("2000", "2010", "2020")))))
exposure <- 2L * y
spec <- Model(y ~ LN2(constraint = constraint))
mod.est <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(mod.est,
along = 3L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = NULL,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "LN2Predict")
expect_identical(as.character(limits(x@constraintLN2)$time), c("2030", "2060"))
expect_true(x@add1@.Data)
spec <- Model(y ~ LN2(constraint = constraint, add1 = FALSE))
mod.est <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(mod.est,
along = 3L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = NULL,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "LN2Predict")
expect_false(x@add1@.Data)
})
StatisticsNZ/demest documentation built on Nov. 2, 2023, 7:56 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.
context("Model-generators")
n.test <- 5
test.identity <- FALSE
test.extended <- FALSE
## addAg ########################################################################
## SpecAgPlaceholder
test_that("addAg works with BinomialVarying and SpecAgPlaceholder", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Binomial(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
ans.obtained <- addAg(model = model,
aggregate = new("SpecAgPlaceholder"))
ans.expected <- model
expect_identical(ans.obtained, ans.expected)
expect_is(ans.obtained, "BinomialVarying")
})
test_that("addAg works with NormalVaryingVarsigmaKnown and SpecAgPlaceholder", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region, sd = 0.5))
set.seed(1)
model <- initialModel(spec, y = y, weights = weights)
ans.obtained <- addAg(model = model,
aggregate = new("SpecAgPlaceholder"))
ans.expected <- model
expect_identical(ans.obtained, ans.expected)
expect_is(ans.obtained, "NormalVaryingVarsigmaKnown")
})
test_that("addAg works with NormalVaryingVarsigmaUnknown and SpecAgPlaceholder", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, weights = weights)
ans.obtained <- addAg(model = model,
aggregate = new("SpecAgPlaceholder"))
ans.expected <- model
expect_identical(ans.obtained, ans.expected)
expect_is(ans.obtained, "NormalVaryingVarsigmaUnknown")
})
test_that("addAg works with PoissonVaryingNotUseExp and SpecAgPlaceholder", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE))
set.seed(1)
model <- initialModel(spec, y = y, exposure = NULL)
ans.obtained <- addAg(model = model,
aggregate = new("SpecAgPlaceholder"))
ans.expected <- model
expect_identical(ans.obtained, ans.expected)
expect_is(ans.obtained, "PoissonVaryingNotUseExp")
})
test_that("addAg works with PoissonVaryingNotUseExp and SpecAgPlaceholder", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
exposure <- y + 1L
spec <- Model(y ~ Poisson(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
ans.obtained <- addAg(model = model,
aggregate = new("SpecAgPlaceholder"))
ans.expected <- model
expect_identical(ans.obtained, ans.expected)
expect_is(ans.obtained, "PoissonVaryingUseExp")
})
## SpecAgCertain
test_that("addAg works with BinomialVarying and SpecAgCertain", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Binomial(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
aggregate <- AgCertain(value = 0.3)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = exposure)
expect_is(ans.obtained, "BinomialVaryingAgCertain")
expect_identical(ans.obtained@slotsToExtract,
new("BinomialVaryingAgCertain")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("BinomialVaryingAgCertain")@iMethodModel)
})
test_that("addAg works with NormalVaryingVarsigmaKnown and SpecAgCertain", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region, sd = 0.5))
set.seed(1)
model <- initialModel(spec, y = y, weights = weights)
aggregate <- AgCertain(value = 0.3)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = weights)
expect_is(ans.obtained, "NormalVaryingVarsigmaKnownAgCertain")
expect_identical(ans.obtained@slotsToExtract,
new("NormalVaryingVarsigmaKnownAgCertain")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("NormalVaryingVarsigmaKnownAgCertain")@iMethodModel)
})
test_that("addAg works with NormalVaryingVarsigmaUnknown and SpecAgCertain", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, weights = weights)
aggregate <- AgCertain(value = 0.3)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = weights)
expect_is(ans.obtained, "NormalVaryingVarsigmaUnknownAgCertain")
expect_identical(ans.obtained@slotsToExtract,
new("NormalVaryingVarsigmaUnknownAgCertain")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("NormalVaryingVarsigmaUnknownAgCertain")@iMethodModel)
})
test_that("addAg works with PoissonVaryingNotUseExp and SpecAgCertain", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
defaultWeights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE))
set.seed(1)
model <- initialModel(spec, y = y, exposure = NULL)
aggregate <- AgCertain(value = 0.3)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = defaultWeights)
expect_is(ans.obtained, "PoissonVaryingNotUseExpAgCertain")
expect_identical(ans.obtained@slotsToExtract,
new("PoissonVaryingNotUseExpAgCertain")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("PoissonVaryingNotUseExpAgCertain")@iMethodModel)
})
test_that("addAg works with PoissonVaryingNotUseExp and SpecAgCertain", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
aggregate <- AgCertain(value = 0.3)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = exposure)
expect_is(ans.obtained, "PoissonVaryingUseExpAgCertain")
expect_identical(ans.obtained@slotsToExtract,
new("PoissonVaryingUseExpAgCertain")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("PoissonVaryingUseExpAgCertain")@iMethodModel)
})
## SpecAgNormal
test_that("addAg works with BinomialVarying and SpecAgNormal", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Binomial(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
aggregate <- AgNormal(value = 0.3, sd = 0.2)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = exposure)
expect_is(ans.obtained, "BinomialVaryingAgNormal")
expect_identical(ans.obtained@slotsToExtract,
new("BinomialVaryingAgNormal")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("BinomialVaryingAgNormal")@iMethodModel)
})
test_that("addAg works with NormalVaryingVarsigmaKnown and SpecAgNormal", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region, sd = 0.5))
set.seed(1)
model <- initialModel(spec, y = y, weights = weights)
aggregate <- AgNormal(value = 0.3, sd = 0.2)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = weights)
expect_is(ans.obtained, "NormalVaryingVarsigmaKnownAgNormal")
expect_identical(ans.obtained@slotsToExtract,
new("NormalVaryingVarsigmaKnownAgNormal")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("NormalVaryingVarsigmaKnownAgNormal")@iMethodModel)
})
test_that("addAg works with NormalVaryingVarsigmaUnknown and SpecAgNormal", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, weights = weights)
aggregate <- AgNormal(value = 0.3, sd = 0.2)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = weights)
expect_is(ans.obtained, "NormalVaryingVarsigmaUnknownAgNormal")
expect_identical(ans.obtained@slotsToExtract,
new("NormalVaryingVarsigmaUnknownAgNormal")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("NormalVaryingVarsigmaUnknownAgNormal")@iMethodModel)
})
test_that("addAg works with PoissonVaryingNotUseExp and SpecAgNormal", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
defaultWeights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE))
set.seed(1)
model <- initialModel(spec, y = y, exposure = NULL)
aggregate <- AgNormal(value = 0.3, sd = 0.2)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = defaultWeights)
expect_is(ans.obtained, "PoissonVaryingNotUseExpAgNormal")
expect_identical(ans.obtained@slotsToExtract,
new("PoissonVaryingNotUseExpAgNormal")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("PoissonVaryingNotUseExpAgNormal")@iMethodModel)
})
test_that("addAg works with PoissonVaryingUseExp and SpecAgNormal", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
aggregate <- AgNormal(value = 0.3, sd = 0.2)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = exposure)
expect_is(ans.obtained, "PoissonVaryingUseExpAgNormal")
expect_identical(ans.obtained@slotsToExtract,
new("PoissonVaryingUseExpAgNormal")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("PoissonVaryingUseExpAgNormal")@iMethodModel)
})
## SpecAgPoisson
test_that("addAg works with BinomialVarying and SpecAgPoisson", {
addAg <- demest:::addAg
expect_error(addAg(model = new("BinomialVarying"),
aggregate = new("SpecAgPoisson")),
"Poisson model for aggregate values can only be used with Poisson likelihood")
})
test_that("addAg works with NormalVaryingVarsigmaKnown and SpecAgPoisson", {
addAg <- demest:::addAg
expect_error(addAg(model = new("NormalVaryingVarsigmaKnown"),
aggregate = new("SpecAgPoisson")),
"Poisson model for aggregate values can only be used with Poisson likelihood")
})
test_that("addAg works with NormalVaryingVarsigmaUnknown and SpecAgPoisson", {
addAg <- demest:::addAg
expect_error(addAg(model = new("NormalVaryingVarsigmaUnknown"),
aggregate = new("SpecAgPoisson")),
"Poisson model for aggregate values can only be used with Poisson likelihood")
})
test_that("addAg works with PoissonVaryingNotUseExp and SpecAgPoisson", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
defaultWeights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE))
set.seed(1)
model <- initialModel(spec, y = y, exposure = NULL)
aggregate <- AgPoisson(value = 0.3)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = defaultWeights)
expect_is(ans.obtained, "PoissonVaryingNotUseExpAgPoisson")
expect_identical(ans.obtained@slotsToExtract,
new("PoissonVaryingNotUseExpAgPoisson")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("PoissonVaryingNotUseExpAgPoisson")@iMethodModel)
})
test_that("addAg works with PoissonVaryingNotUseExp and SpecAgPoisson", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
aggregate <- AgPoisson(value = 0.3)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = exposure)
expect_is(ans.obtained, "PoissonVaryingUseExpAgPoisson")
expect_identical(ans.obtained@slotsToExtract,
new("PoissonVaryingUseExpAgPoisson")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("PoissonVaryingUseExpAgPoisson")@iMethodModel)
})
## SpecAgFun
test_that("addAg works with BinomialVarying and SpecAgFun", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Binomial(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
aggregate <- AgFun(value = 0.3, sd = 0.2, FUN = function(x, weights) 1)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = exposure)
expect_is(ans.obtained, "BinomialVaryingAgFun")
expect_identical(ans.obtained@slotsToExtract,
new("BinomialVaryingAgFun")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("BinomialVaryingAgFun")@iMethodModel)
})
test_that("addAg works with NormalVaryingVarsigmaKnown and SpecAgFun", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region, sd = 0.5))
set.seed(1)
model <- initialModel(spec, y = y, weights = weights)
aggregate <- AgFun(value = 0.3, sd = 0.2, FUN = function(x, weights) 1)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = weights)
expect_is(ans.obtained, "NormalVaryingVarsigmaKnownAgFun")
expect_identical(ans.obtained@slotsToExtract,
new("NormalVaryingVarsigmaKnownAgFun")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("NormalVaryingVarsigmaKnownAgFun")@iMethodModel)
})
test_that("addAg works with NormalVaryingVarsigmaUnknown and SpecAgFun", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, weights = weights)
aggregate <- AgFun(value = 0.3, sd = 0.2, FUN = function(x, weights) 1)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = weights)
expect_is(ans.obtained, "NormalVaryingVarsigmaUnknownAgFun")
expect_identical(ans.obtained@slotsToExtract,
new("NormalVaryingVarsigmaUnknownAgFun")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("NormalVaryingVarsigmaUnknownAgFun")@iMethodModel)
})
test_that("addAg works with PoissonVaryingUseExp and SpecAgFun", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
y <- Counts(array(rbinom(n = 20, size = 20, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])),
dimscales = c(age = "Intervals"))
defaultWeights <- Counts(array(rbinom(n = 20, size = 20, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])),
dimscales = c(age = "Intervals"))
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE))
set.seed(1)
model <- initialModel(spec, y = y, exposure = NULL)
aggregate <- AgFun(value = 0.3, sd = 0.2, FUN = function(x, weights) 1)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = defaultWeights)
expect_is(ans.obtained, "PoissonVaryingNotUseExpAgFun")
expect_identical(ans.obtained@slotsToExtract,
new("PoissonVaryingNotUseExpAgFun")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("PoissonVaryingNotUseExpAgFun")@iMethodModel)
})
test_that("addAg works with PoissonVaryingUseExp and SpecAgFun", {
addAg <- demest:::addAg
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region))
set.seed(1)
model <- initialModel(spec, y = y, exposure = exposure)
aggregate <- AgFun(value = 0.3, sd = 0.2, FUN = function(x, weights) 1)
ans.obtained <- addAg(model = model,
aggregate = aggregate,
defaultWeights = exposure)
expect_is(ans.obtained, "PoissonVaryingUseExpAgFun")
expect_identical(ans.obtained@slotsToExtract,
new("PoissonVaryingUseExpAgFun")@slotsToExtract)
expect_identical(ans.obtained@iMethodModel,
new("PoissonVaryingUseExpAgFun")@iMethodModel)
})
## initialModel - Varying ##################################################################
test_that("initialModel creates object of class BinomialVarying from valid inputs", {
initialModel <- demest:::initialModel
makeLinearBetas <- demest:::makeLinearBetas
jitterBetas <- demest:::jitterBetas
makePriors <- demest:::makePriors
## main effect
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Binomial(mean ~ age + region))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
m <- sum(y) / sum(exposure)
nu <- m * (1 - m)
sigma <- new("Scale", runif(1, 0, nu))
theta <- rbeta(n = 20,
shape1 = m * (nu/sigma@.Data^2 - 1) + y,
shape2 = (1-m) * (nu/sigma@.Data^2 - 1) + exposure - y)
logit.theta <- array(log(theta / (1 - theta)), dim = dim(y), dimnames = dimnames(y))
betas <- makeLinearBetas(theta = logit.theta, formula = prob ~ age + region)
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = NULL,
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priors)
expect_identical(x@sigma, sigma)
expect_identical(x@theta, theta)
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@ASigma@.Data, 1)
expect_identical(x@nuSigma@.Data, 7)
expect_identical(x@sigmaMax@.Data, qhalft(p = 0.999, df = 7, scale = 1))
expect_identical(x@lower, -Inf)
expect_identical(x@upper, Inf)
expect_identical(x@tolerance, 1e-5)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@nFailedPropTheta, new("Counter", 0L))
## intercept only
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Binomial(mean ~ 1))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
m <- sum(y) / sum(exposure)
nu <- m * (1 - m)
sigma <- runif(1, 0, nu)
theta <- rbeta(n = 20,
shape1 = m * (nu/sigma^2 -1) + y,
shape2 = (1-m) * (nu/sigma^2 - 1) + exposure - y)
sigma <- new("Scale", sigma)
betas <- list("(Intercept)" = mean(log(theta / (1 - theta))))
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L),
y = y,
sY = NULL,
strucZeroArray = strucZeroArray)
betas <- unname(betas)
betas <- jitterBetas(betas, priors)
expect_identical(x@sigma, sigma)
expect_identical(x@theta, theta)
expect_identical(x@betas, betas)
expect_identical(sapply(x@priorsBetas, class), "ExchFixed")
## lower and upper specified
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Binomial(mean ~ age + region), lower = 0.4, upper = 0.8)
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(all(x@theta >= 0.4))
expect_true(all(x@theta <= 0.8))
expect_equal(x@lower, log(0.4 / 0.6))
expect_equal(x@upper, log(0.8 / 0.2))
## missing values
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y[1:4] <- NA
exposure[1:3] <- NA
spec <- Model(y ~ Binomial(mean ~ age + region))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
m <- sum(y[5:20]) / sum(exposure[5:20])
nu <- m * (1 - m)
sigma <- runif(1, 0, nu)
sigma <- new("Scale", sigma)
y[1:4] <- 0
exposure[1:4] <- 0
theta <- rbeta(n = 20,
shape1 = m * (nu/sigma@.Data^2 - 1) + y,
shape2 = (1-m) * (nu/sigma@.Data^2 - 1) + exposure - y)
expect_identical(x@sigma, sigma)
expect_identical(x@theta, theta)
## saturated model
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Binomial(mean ~ age * region),
age:region ~ Exch(error = Error(scale = HalfT(df = 3, scale = 0.1, max = 0.6))))
x <- initialModel(spec, y = y, exposure = exposure)
expect_identical(x@ASigma, x@priorsBetas[[4]]@ATau)
expect_identical(x@sigmaMax, x@priorsBetas[[4]]@tauMax)
expect_identical(x@nuSigma, x@priorsBetas[[4]]@nuTau)
})
test_that("initialModel methods for Varying model can cope with orig-dest dimtypes", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 25, lambda = 20),
dim = c(5, 5),
dimnames = list(reg_orig = c("a", "b", "c", "d", "e"),
reg_dest = c("a", "b", "c", "d", "e"))))
y <- Counts(array(rbinom(n = 25, size = exposure, prob = 0.5),
dim = c(5, 5),
dimnames = list(reg_orig = c("a", "b", "c", "d", "e"),
reg_dest = c("a", "b", "c", "d", "e"))))
spec <- Model(y ~ Binomial(mean ~ reg_orig + reg_dest))
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
})
test_that("initialModel creates object of class CMPVaryingUseExp from valid inputs", {
initialModel <- demest:::initialModel
makeSigmaInitialPoisson <- demest:::makeSigmaInitialPoisson
loglm <- MASS:::loglm
makePriors <- demest:::makePriors
jitterBetas <- demest:::jitterBetas
## no missing values
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ CMP(mean ~ age + region))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
theta <- rgamma(n = 20,
shape = 0.05 * mean(y) + 0.95 * y,
rate = 0.05 * mean(exposure) + 0.95 * exposure)
sigma <- runif(1, max = 1)
theta <- array(theta, dim = dim(y), dimnames = dimnames(y))
betas <- loglm(~ age + region, theta)$param
theta <- as.numeric(theta)
logNuCMP <- stats::rnorm(n = length(theta))
nuCMP <- exp(logNuCMP)
nuCMP[nuCMP < 0.5] <- 0.5
nuCMP[nuCMP > 2] <- 2
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = NULL,
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priorsBetas = priors)
expect_is(x, "CMPVaryingUseExp")
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@ASigma, new("Scale", 1))
expect_identical(x@nuSigma, new("DegreesFreedom", 7))
expect_identical(x@sigmaMax, new("Scale", qhalft(0.999, df = 7, scale = 1)))
expect_identical(x@tolerance, 1e-5)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@nFailedPropTheta, new("Counter", 0L))
expect_identical(x@nuCMP, new("ParameterVector", nuCMP))
expect_identical(x@meanLogNuCMP, new("Parameter", 0))
expect_identical(x@sdLogNuCMP, new("Scale", 1))
## intercept only
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ CMP(mean ~ 1))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
theta <- rgamma(n = 20,
shape = 0.05 * mean(y) + 0.95 * y,
rate = 0.05 * mean(exposure) + 0.95 * exposure)
sigma <- runif(1, max = 1)
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_identical(sapply(x@priorsBetas, class), "ExchFixed")
## specify upper and lower
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ CMP(mean ~ age + region), lower = 0.01, upper = 2)
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(all(x@theta >= 0.01))
expect_true(all(x@theta <= 2))
expect_equal(x@lower, log(0.01))
expect_equal(x@upper, log(2))
## has missing values
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y[1:5] <- NA
spec <- Model(y ~ CMP(mean ~ age + region))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
mean.y <- mean(y[6:20])
mean.expose <- mean(exposure[6:20])
shape <- 0.05 * mean.y + 0.95 * y
shape[1:5] <- mean.y
rate <- 0.05 * mean.expose + 0.95 * exposure
rate[1:5] <- mean.expose
theta <- rgamma(n = 20, shape = shape, rate = rate)
sigma <- runif(1)
theta <- array(theta, dim = dim(y), dimnames = dimnames(y))
betas <- loglm(~ age + region, theta)$param
theta <- as.numeric(theta)
nuCMP <- rnorm(n = length(theta), mean = 0, sd = 1)
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = NULL,
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priorsBetas = priors)
expect_is(x, "CMPVaryingUseExp")
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@nuSigma, new("DegreesFreedom", 7))
expect_identical(x@ASigma, new("Scale", 1))
expect_identical(x@tolerance, 1e-5)
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@nFailedPropTheta, new("Counter", 0L))
## saturated model
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ CMP(mean ~ age * region),
age:region ~ Exch(error = Error(scale = HalfT(df = 3, scale = 0.1, max = 0.6))))
x <- initialModel(spec, y = y, exposure = exposure)
expect_identical(x@ASigma, x@priorsBetas[[4]]@ATau)
expect_identical(x@sigmaMax, x@priorsBetas[[4]]@tauMax)
expect_identical(x@nuSigma, x@priorsBetas[[4]]@nuTau)
## must have exposure
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ CMP(mean ~ age * region))
expect_error(initialModel(spec, y = y, exposure = NULL),
"model 'y ~ CMP\\(mean ~ age \\* region\\)' uses exposure, but no 'exposure' argument supplied")
})
test_that("initialModel creates object of class NormalVaryingVarsigmaKnown from valid inputs", {
initialModel <- demest:::initialModel
makeLinearBetas <- demest:::makeLinearBetas
makePriors <- demest:::makePriors
jitterBetas <- demest:::jitterBetas
## y has no missing values
varsigma <- 1.3
y <- Counts(array(rnorm(n = 20, mean = 0, sd = varsigma),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region, sd = varsigma))
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
set.seed(1)
sigma <- runif(n = 1, max = sd(y))
theta <- rnorm(n = 20, mean = 0.5 * mean(y) + 0.5 * y, sd = sigma)
theta <- array(theta, dim = dim(y), dimnames = dimnames(y))
betas <- makeLinearBetas(theta = theta, formula = mean ~ age + region)
theta <- as.numeric(theta)
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = sd(y),
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priors)
expect_is(x, "NormalVaryingVarsigmaKnown")
expect_identical(x@theta, theta)
expect_identical(x@varsigma, new("Scale", varsigma))
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@nuSigma, new("DegreesFreedom", 7))
expect_identical(x@ASigma, new("Scale", sd(y)))
expect_identical(x@sigmaMax@.Data, qhalft(p = 0.999, df = 7, scale = sd(y)))
expect_identical(x@lower, -Inf)
expect_identical(x@upper, Inf)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@tolerance, 1e-5)
## intercept only
varsigma <- 2.1
y <- Counts(array(rnorm(n = 20, mean = 0, sd = varsigma),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ 1, sd = varsigma))
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
set.seed(1)
sigma <- runif(1, max = sd(y))
theta <- rnorm(n = 20, mean = 0.5 * mean(y) + 0.5 * y, sd = sigma)
betas <- list(mean(theta))
betas <- jitterBetas(betas)
expect_identical(x@theta, theta)
expect_identical(x@betas, betas)
expect_identical(sapply(x@priorsBetas, class), "ExchFixed")
## lower and upper specified
varsigma <- 1.3
y <- Counts(array(rnorm(n = 20, mean = 0, sd = varsigma),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region, sd = varsigma), lower = 0, upper = 3)
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
expect_true(all(x@theta > 0))
expect_true(all(x@theta < 3))
expect_equal(x@lower, 0)
expect_equal(x@upper, 3)
## y has two values
varsigma <- 2
y <- Counts(array(rnorm(n = 20, mean = 0, sd = varsigma),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y[-(1:2)] <- NA
weights[-(1:2)] <- NA
spec <- Model(y ~ Normal(mean ~ age + region, sd = 2))
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
set.seed(1)
sigma <- runif(1, max = sd(y, na.rm = T))
y.bar <- mean(y, na.rm = T)
mean <- rep(y.bar, 20)
mean[1:2] <- 0.5 * y.bar + 0.5 * y[1:2]
theta <- rnorm(20, mean = mean, sd = sigma)
theta <- array(theta, dim = dim(y), dimnames = dimnames(y))
betas <- makeLinearBetas(theta = theta, formula = mean ~ age + region)
theta <- as.numeric(theta)
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = sd(y, na.rm = T),
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priors)
expect_is(x, "NormalVaryingVarsigmaKnown")
expect_identical(x@theta, theta)
expect_identical(x@varsigma, new("Scale", varsigma))
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@nuSigma, new("DegreesFreedom", 7))
expect_identical(x@ASigma, new("Scale", sd(y, na.rm = T)))
expect_identical(x@lower, -Inf)
expect_identical(x@upper, Inf)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@tolerance, 1e-5)
## maxSigma specified
varsigma <- 1.3
y <- Counts(array(rnorm(n = 20, mean = 0, sd = varsigma),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region, sd = varsigma),
priorSD = HalfT(max = 1000))
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
expect_identical(x@sigmaMax@.Data, 1000)
## saturated model
y <- Values(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age * region, sd = 1),
age:region ~ Exch(error = Error(scale = HalfT(df = 3, scale = 0.1, max = 0.6))))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
x <- initialModel(spec, y = y, weights = weights)
expect_identical(x@ASigma, x@priorsBetas[[4]]@ATau)
expect_identical(x@sigmaMax, x@priorsBetas[[4]]@tauMax)
expect_identical(x@nuSigma, x@priorsBetas[[4]]@nuTau)
## varsigma is 0
y <- Values(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age * region, sd = 0),
age:region ~ Exch(error = Error(scale = HalfT(df = 3, scale = 0.1, max = 0.6))))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
x <- initialModel(spec, y = y, weights = weights)
expect_identical(x@varsigma@.Data, 0)
expect_identical(x@theta, as.numeric(y))
})
test_that("initialModel creates object of class NormalVaryingVarsigmaUnknown from valid inputs", {
initialModel <- demest:::initialModel
rinvchisq1 <- demest:::rinvchisq1
makeLinearBetas <- demest:::makeLinearBetas
makePriors <- demest:::makePriors
jitterBetas <- demest:::jitterBetas
## no missing values
y <- Counts(array(rnorm(n = 20, mean = 0),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region),
priorSD = HalfT(max = 20))
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
set.seed(1)
I <- length(y)
varsigma <- runif(1, max = sd(y))
sigma <- runif(1, max = sd(y))
theta <- rnorm(n = 20, mean = 0.5 * mean(y) + 0.5 * y, sd = sigma)
theta <- array(theta, dim = dim(y), dimnames = dimnames(y))
betas <- makeLinearBetas(theta = theta, formula = mean ~ age + region)
theta <- as.numeric(theta)
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = sd(y),
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priors)
expect_is(x, "NormalVaryingVarsigmaUnknown")
expect_identical(x@theta, theta)
expect_identical(x@varsigma, new("Scale", varsigma))
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@sigmaMax, new("Scale", 20))
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@lower, -Inf)
expect_identical(x@upper, Inf)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@tolerance, 1e-5)
## intercept only
y <- Counts(array(rnorm(n = 20, mean = 0, sd = varsigma),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ 1))
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
set.seed(1)
I <- length(y)
varsigma <- runif(1, max = sd(y))
sigma <- runif(1, max = sd(y))
theta <- rnorm(n = 20, mean = 0.5 * mean(y) + 0.5 * y, sd = sigma)
betas <- list(mean(theta))
betas <- jitterBetas(betas)
expect_identical(x@theta, theta)
expect_identical(x@betas, betas)
expect_identical(sapply(x@priorsBetas, class), "ExchFixed")
## lower and upper specified
y <- Counts(array(rnorm(n = 20, mean = 0, sd = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region), lower = 0, upper = 3)
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
expect_true(all(x@theta > 0))
expect_true(all(x@theta < 3))
expect_equal(x@lower, 0)
expect_equal(x@upper, 3)
## has missing values
y <- Counts(array(rnorm(n = 20, mean = 0),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y[1:10] <- NA
weights[1:5] <- NA
spec <- Model(y ~ Normal(mean ~ age + region))
set.seed(1)
x <- initialModel(spec, y = y, weights = weights)
set.seed(1)
varsigma <- runif(n = 1, max = sd(y, na.rm = T))
sigma <- runif(n = 1, max = sd(y, na.rm = T))
mean <- 0.5 * mean(y[11:20]) + 0.5 * y
mean[1:10] <- mean(y[11:20])
theta <- rnorm(n = 20, mean = mean, sd = sigma)
theta <- array(theta, dim = dim(y), dimnames = dimnames(y))
betas <- makeLinearBetas(theta = theta, formula = mean ~ age + region)
theta <- as.numeric(theta)
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = sd(y, na.rm = T),
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priors)
expect_is(x, "NormalVaryingVarsigmaUnknown")
expect_identical(x@theta, theta)
expect_identical(x@varsigma, new("Scale", varsigma))
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@lower, -Inf)
expect_identical(x@upper, Inf)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@tolerance, 1e-5)
## saturated model
y <- Values(array(rnorm(n = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age * region),
age:region ~ Exch(error = Error(scale = HalfT(df = 3, scale = 0.1, max = 0.6))))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
x <- initialModel(spec, y = y, weights = weights)
expect_identical(x@ASigma, x@priorsBetas[[4]]@ATau)
expect_identical(x@sigmaMax, x@priorsBetas[[4]]@tauMax)
expect_identical(x@nuSigma, x@priorsBetas[[4]]@nuTau)
})
test_that("initialModel creates object of class PoissonVaryingUseExp from valid inputs", {
initialModel <- demest:::initialModel
makeSigmaInitialPoisson <- demest:::makeSigmaInitialPoisson
loglm <- MASS:::loglm
makePriors <- demest:::makePriors
jitterBetas <- demest:::jitterBetas
## no missing values
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
theta <- rgamma(n = 20,
y + 1,
rate = exposure + 1)
sigma <- runif(1, max = 1)
theta <- array(theta, dim = dim(y), dimnames = dimnames(y))
betas <- loglm(~ age + region, theta)$param
theta <- as.numeric(theta)
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = NULL,
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priors)
expect_is(x, "PoissonVaryingUseExp")
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@ASigma, new("Scale", 1))
expect_identical(x@nuSigma, new("DegreesFreedom", 7))
expect_identical(x@sigmaMax, new("Scale", qhalft(0.999, df = 7, scale = 1)))
expect_identical(x@tolerance, 1e-5)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@nFailedPropTheta, new("Counter", 0L))
## intercept only
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ 1))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
theta <- rgamma(n = 20,
shape = y + 1,
rate = exposure + 1)
sigma <- runif(1, max = 1)
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_identical(sapply(x@priorsBetas, class), "ExchFixed")
## specify upper and lower
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region), lower = 0.01, upper = 2)
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(all(x@theta >= 0.01))
expect_true(all(x@theta <= 2))
expect_equal(x@lower, log(0.01))
expect_equal(x@upper, log(2))
## has missing values
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y[1:5] <- NA
spec <- Model(y ~ Poisson(mean ~ age + region))
set.seed(1)
x <- initialModel(spec, y = y, exposure = exposure)
set.seed(1)
mean.y <- mean(y[6:20])
mean.expose <- mean(exposure[6:20])
shape <- ifelse(is.na(y), 1, y + 1)
rate <- ifelse(is.na(y), 1, exposure + 1)
theta <- rgamma(n = 20, shape = shape, rate = rate)
sigma <- runif(1)
theta <- array(theta, dim = dim(y), dimnames = dimnames(y))
betas <- loglm(~ age + region, theta)$param
theta <- as.numeric(theta)
strucZeroArray <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
priors <- makePriors(betas = betas,
specs = spec@specsPriors,
namesSpecs = spec@namesSpecsPriors,
margins = list(0L, 1L, 2L),
y = y,
sY = NULL,
strucZeroArray = strucZeroArray)
betas <- unname(lapply(betas, as.numeric))
betas <- jitterBetas(betas, priors)
expect_is(x, "PoissonVaryingUseExp")
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_identical(x@betas, betas)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@nuSigma, new("DegreesFreedom", 7))
expect_identical(x@ASigma, new("Scale", 1))
expect_identical(x@tolerance, 1e-5)
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@nFailedPropTheta, new("Counter", 0L))
## saturated model
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age * region),
age:region ~ Exch(error = Error(scale = HalfT(df = 3, scale = 0.1, max = 0.6))))
x <- initialModel(spec, y = y, exposure = exposure)
expect_identical(x@ASigma, x@priorsBetas[[4]]@ATau)
expect_identical(x@sigmaMax, x@priorsBetas[[4]]@tauMax)
expect_identical(x@nuSigma, x@priorsBetas[[4]]@nuTau)
## must have exposure
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age * region))
expect_error(initialModel(spec, y = y, exposure = NULL),
"model 'y ~ Poisson\\(mean ~ age \\* region\\)' uses exposure, but no 'exposure' argument supplied")
})
test_that("initialModel creates object of class PoissonVaryingNotUseExp from valid inputs", {
initialModel <- demest:::initialModel
makeSigmaInitialPoisson <- demest:::makeSigmaInitialPoisson
## no missing values
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE),
priorSD = HalfT(df = 6, max = 10))
set.seed(1)
x <- initialModel(spec, y = y, exposure = NULL)
set.seed(1)
theta <- rgamma(n = 20,
shape = y + 1,
rate = 2)
sigma <- runif(1, max = sd(log(y+1)))
expect_is(x, "PoissonVaryingNotUseExp")
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@nuSigma, new("DegreesFreedom", 6))
expect_identical(x@ASigma, new("Scale", sd(log(y+1))))
expect_identical(x@sigmaMax, new("Scale", 10))
expect_identical(x@tolerance, 1e-5)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@nFailedPropTheta, new("Counter", 0L))
## intercept only
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ 1, useExpose = FALSE))
set.seed(1)
x <- initialModel(spec, y = y, exposure = NULL)
set.seed(1)
theta <- rgamma(n = 20,
shape = y + 1,
rate = 2)
sigma <- runif(1, max = sd(log(y+1)))
m <- mean(y)
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_identical(sapply(x@priorsBetas, class), "ExchFixed")
## has missing values
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y[1:5] <- NA
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE))
set.seed(1)
x <- initialModel(spec, y = y, exposure = NULL)
set.seed(1)
mean <- mean(y, na.rm = TRUE)
shape <- ifelse(is.na(y), 1, y + 1)
rate <- 2
theta <- rgamma(n = 20, shape = shape, rate = rate)
sigma <- runif(1, max = sd(log(y+1), na.rm = T))
expect_is(x, "PoissonVaryingNotUseExp")
expect_identical(x@sigma, new("Scale", sigma))
expect_identical(x@theta, theta)
expect_true(all(sapply(x@priorsBetas, is, "Prior")))
expect_is(x@priorsBetas[[1]], "ExchFixed")
expect_identical(x@nuSigma, new("DegreesFreedom", 7))
expect_identical(x@ASigma, new("Scale", sd(log(y+1), na.rm = TRUE)))
expect_identical(x@tolerance, 1e-5)
expect_identical(x@maxAttempt, 1000L)
expect_identical(x@nFailedPropTheta, new("Counter", 0L))
## saturated model
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age * region, useExpose = FALSE),
age:region ~ Exch(error = Error(scale = HalfT(df = 3, scale = 0.1, max = 0.6))))
x <- initialModel(spec, y = y, exposure = NULL)
expect_identical(x@ASigma, x@priorsBetas[[4]]@ATau)
expect_identical(x@sigmaMax, x@priorsBetas[[4]]@tauMax)
expect_identical(x@nuSigma, x@priorsBetas[[4]]@nuTau)
## must not have exposure
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age * region, useExpose = FALSE))
expect_error(initialModel(spec, y = y, exposure = exposure),
"exposure' argument supplied, but model 'y ~ Poisson\\(mean ~ age \\* region, useExpose = FALSE\\)' does not use exposure")
})
## PoissonBinomialMixture #############################################################
test_that("initialModel creates object of class PoissonBinomialMixture from valid inputs", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ PoissonBinomial(prob = 0.98))
x <- initialModel(spec, y = y, exposure = exposure)
expect_equal(x,
new("PoissonBinomialMixture",
call = call("Model", formula = y ~ PoissonBinomial(prob = 0.98)),
prob = 0.98,
metadataY = y@metadata))
})
## Round3 #############################################################
test_that("initialModel creates object of class Round3 from valid inputs", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- round3(y)
spec <- Model(y ~ Round3())
x <- initialModel(spec, y = y, exposure = exposure)
expect_equal(x,
new("Round3",
call = call("Model", formula = y ~ Round3()),
metadataY = y@metadata))
expect_error(initialModel(spec, y = y + 1L, exposure = exposure),
"using 'Round3' data model, but data contains values not divisible by 3")
})
## Exact #############################################################
test_that("initialModel creates object of class Exact from valid inputs", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Exact())
x <- initialModel(spec, y = y, exposure = exposure)
expect_equal(x,
new("Exact",
call = call("Model", formula = y ~ Exact()),
metadataY = y@metadata))
y[1] <- NA
expect_error(initialModel(spec, y = y, exposure = exposure),
"using 'Exact' data model, but data contains missing values")
})
## NormalFixed #############################################################
test_that("initialModel creates object of class NormalFixedUseExp from valid inputs", {
initialModel <- demest:::initialModel
mean <- Values(array(runif(n = 25),
dim = c(5, 5),
dimnames = list(age = 0:4, region = letters[1:5])))
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ NormalFixed(mean = mean, sd = 0.1))
ans.obtained <- initialModel(spec, y = y, exposure = exposure)
ans.expected <- new("NormalFixedUseExp",
call = call("Model", formula = y ~ NormalFixed(mean = mean, sd = 0.1)),
mean = new("ParameterVector", mean@.Data[1:20]),
sd = new("ScaleVec", rep(0.1, 20)),
metadataY = y@metadata,
meanAll = new("ParameterVector", mean@.Data),
sdAll = new("ScaleVec", rep(0.1, 25)),
metadataAll = mean@metadata)
expect_equal(ans.obtained, ans.expected)
})
test_that("initialModel throws appropriate errors with NormalFixedUseExp", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
mean.wrong <- Values(array(runif(n = 25),
dim = c(5, 5),
dimnames = list(age = 0:4, region = letters[2:6])))
spec <- Model(y ~ NormalFixed(mean = mean.wrong, sd = 0.1))
expect_error(initialModel(spec, y = y, exposure = exposure),
"'mean' from NormalFixed model not compatible with data :")
mean <- Values(array(runif(n = 25),
dim = c(5, 5),
dimnames = list(age = 0:4, region = letters[1:5])))
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ NormalFixed(mean = mean, sd = 0.1))
expect_error(initialModel(spec, y = y, exposure = NULL),
"model 'y ~ NormalFixed\\(mean = mean, sd = 0\\.1\\)' uses exposure, but no 'exposure' argument supplied")
})
test_that("initialModel creates object of class NormalFixedNotUseExp from valid inputs", {
initialModel <- demest:::initialModel
mean <- Values(array(runif(n = 25),
dim = c(5, 5),
dimnames = list(age = 0:4, region = letters[1:5])))
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ NormalFixed(mean = mean, sd = 0.1, useExpose = FALSE))
ans.obtained <- initialModel(spec, y = y, exposure = NULL)
ans.expected <- new("NormalFixedNotUseExp",
call = call("Model", formula = y ~ NormalFixed(mean = mean, sd = 0.1, useExpose = FALSE)),
mean = new("ParameterVector", mean@.Data[1:20]),
sd = new("ScaleVec", rep(0.1, 20)),
metadataY = y@metadata,
meanAll = new("ParameterVector", mean@.Data),
sdAll = new("ScaleVec", rep(0.1, 25)),
metadataAll = mean@metadata)
expect_equal(ans.obtained, ans.expected)
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ NormalFixed(mean = mean, sd = 0.1, useExpose = FALSE))
expect_error(initialModel(spec, y = y, exposure = exposure),
"'exposure' argument supplied, but model 'y ~ NormalFixed\\(mean = mean, sd = 0.1, useExpose = FALSE\\)' does not use exposure")
})
## TFixed #############################################################
test_that("initialModel creates object of class TFixedUseExp from valid inputs", {
initialModel <- demest:::initialModel
location <- Values(array(runif(n = 25),
dim = c(5, 5),
dimnames = list(age = 0:4, region = letters[1:5])))
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ TFixed(location = location, scale = 0.1, df = 8))
ans.obtained <- initialModel(spec, y = y, exposure = exposure)
ans.expected <- new("TFixedUseExp",
call = call("Model", formula = y ~ TFixed(location = location, scale = 0.1, df = 8)),
mean = new("ParameterVector", location@.Data[1:20]),
sd = new("ScaleVec", rep(0.1, 20)),
metadataY = y@metadata,
meanAll = new("ParameterVector", location@.Data),
sdAll = new("ScaleVec", rep(0.1, 25)),
nu = new("DegreesFreedom", 8),
metadataAll = location@metadata)
expect_equal(ans.obtained, ans.expected)
})
test_that("initialModel throws appropriate errors with TFixedUseExp", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
location.wrong <- Values(array(runif(n = 25),
dim = c(5, 5),
dimnames = list(age = 0:4, region = letters[2:6])))
spec <- Model(y ~ TFixed(location = location.wrong, scale = 0.1))
expect_error(initialModel(spec, y = y, exposure = exposure),
"'location' from TFixed model not compatible with data :")
location <- Values(array(runif(n = 25),
dim = c(5, 5),
dimnames = list(age = 0:4, region = letters[1:5])))
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ TFixed(location = location, scale = 0.1))
expect_error(initialModel(spec, y = y, exposure = NULL),
"model 'y ~ TFixed\\(location = location, scale = 0\\.1\\)' uses exposure, but no 'exposure' argument supplied")
})
test_that("initialModel creates object of class TFixedNotUseExp from valid inputs", {
initialModel <- demest:::initialModel
location <- Values(array(runif(n = 25),
dim = c(5, 5),
dimnames = list(age = 0:4, region = letters[1:5])))
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ TFixed(location = location, scale = 0.1, useExpose = FALSE))
ans.obtained <- initialModel(spec, y = y, exposure = NULL)
ans.expected <- new("TFixedNotUseExp",
call = call("Model", formula = y ~ TFixed(location = location, scale = 0.1, useExpose = FALSE)),
mean = new("ParameterVector", location@.Data[1:20]),
sd = new("ScaleVec", rep(0.1, 20)),
metadataY = y@metadata,
meanAll = new("ParameterVector", location@.Data),
sdAll = new("ScaleVec", rep(0.1, 25)),
nu = new("DegreesFreedom", 7),
metadataAll = location@metadata)
expect_equal(ans.obtained, ans.expected)
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.3 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ TFixed(location = location, scale = 0.1, useExpose = FALSE))
expect_error(initialModel(spec, y = y, exposure = exposure),
"'exposure' argument supplied, but model 'y ~ TFixed\\(location = location, scale = 0.1, useExpose = FALSE\\)' does not use exposure")
})
## LN2 ###############################################################################
test_that("initialModel creates object of class LN2 from valid inputs", {
initialModel <- demest:::initialModel
constraint <- Values(array(c(NA, -1L, 0L, 1L),
dim = c(2, 2),
dimnames = list(age = c("0-39", "40+"),
sex = c("Female", "Male"))))
y <- Counts(array(10L,
dim = c(2, 4, 3),
dimnames = c(list(sex = c("Female", "Male"),
age = c("0-19", "20-39", "40-59", "60+"),
time = c("2000", "2010", "2020")))))
exposure <- 2L * y
spec <- Model(y ~ LN2(constraint = constraint))
ans.obtained <- initialModel(spec,
y = y,
exposure = exposure)
expect_true(validObject(ans.obtained))
expect_true(ans.obtained@add1@.Data)
y[2,1:2,] <- 0L
sz <- Values(array(c(1L, 1L, 0L, 1L),
dim = c(2, 2),
dimnames = list(age = c("0-39", "40+"),
sex = c("Female", "Male"))))
concordances <- list(sex = Concordance(data.frame(from = c("F", "M", "Female", "Male"),
to = c("Female", "Male", "Female", "Male"))))
spec <- Model(y ~ LN2(constraint = constraint,
structuralZeros = sz,
concordances = concordances))
ans.obtained <- initialModel(spec,
y = y,
exposure = exposure)
expect_true(validObject(ans.obtained))
y[2,1:2,] <- 0L
sz <- Values(array(c(1L, 1L, 0L, 1L),
dim = c(2, 2),
dimnames = list(age = c("0-39", "40+"),
sex = c("Female", "Male"))))
concordances <- list(sex = Concordance(data.frame(from = c("F", "M", "Female", "Male"),
to = c("Female", "Male", "Female", "Male"))))
spec <- Model(y ~ LN2(constraint = constraint,
structuralZeros = sz,
concordances = concordances,
sd = 0.2))
ans.obtained <- initialModel(spec,
y = y,
exposure = exposure)
expect_true(validObject(ans.obtained))
expect_identical(ans.obtained@varsigma@.Data, 0.2)
sz <- Values(array(c(1L, 1L, 0L, 1L),
dim = c(2, 2),
dimnames = list(age = c("0-39", "40+"),
sex = c("Female", "Male"))))
spec <- Model(y ~ LN2(constraint = constraint,
structuralZeros = sz,
sd = InvChiSq(df = 5, scaleSq = 10)))
ans.obtained <- initialModel(spec,
y = y,
exposure = exposure)
expect_true(validObject(ans.obtained))
expect_false(ans.obtained@varsigmaLN2HasHalfT@.Data)
constraint <- Values(array(c(NA, -1L, 0L, 1L),
dim = c(2, 2),
dimnames = list(age = c("0-39", "40+"),
sex = c("Female", "Male"))))
y <- Counts(array(10L,
dim = c(2, 4, 3),
dimnames = c(list(sex = c("Female", "Male"),
age = c("0-19", "20-39", "40-59", "60+"),
time = c("2000", "2010", "2020")))))
exposure <- 2L * y
spec <- Model(y ~ LN2(constraint = constraint, add1 = FALSE))
ans.obtained <- initialModel(spec,
y = y,
exposure = exposure)
expect_true(validObject(ans.obtained))
expect_false(ans.obtained@add1@.Data)
})
## Aggregate #########################################################################
test_that("initialModel works with BinomialVarying and AgCertain", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])),
dimscales = c(age = "Intervals"))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])),
dimscales = c(age = "Intervals"))
## scalar
aggregate <- AgCertain(value = 0.4)
spec <- Model(y ~ Binomial(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "BinomialVaryingAgCertain")
## Values
value <- Values(array(c(0.2, 0.3, 0.4),
dim = 3,
dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgCertain(value = value)
spec <- Model(y ~ Binomial(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "BinomialVaryingAgCertain")
})
test_that("initialModel works with BinomialVarying and AgNormal", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgNormal(value = 0.4, sd = 0.1)
spec <- Model(y ~ Binomial(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "BinomialVaryingAgNormal")
## Values
value <- Values(array(c(0.2, 0.3, 0.4), dim = 3, dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgNormal(value = value, sd = sqrt(value))
spec <- Model(y ~ Binomial(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "BinomialVaryingAgNormal")
})
test_that("initialModel works with NormalVaryingVarsigmaKnown and AgCertain", {
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20, mean = 0, sd = 0.3),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgCertain(value = 0.4)
spec <- Model(y ~ Normal(mean ~ age + region, sd = 0.2),
aggregate = aggregate)
x <- initialModel(spec, y = y, weights = weights)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaKnownAgCertain")
## Values
value <- Values(array(c(0.2, 0.3, 0.4),
dim = 3,
dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgCertain(value = value)
spec <- Model(y ~ Normal(mean ~ age + region, sd = 0.2),
aggregate = aggregate)
x <- initialModel(spec, y = y, weights = weights)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaKnownAgCertain")
})
test_that("initialModel works with NormalVaryingVarsigmaKnown and AgNormal", {
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20, mean = 0, sd = 0.3),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgNormal(value = 0.4, sd = 0.1)
spec <- Model(y ~ Normal(mean ~ age + region, sd = 0.2),
aggregate = aggregate)
x <- initialModel(spec, y = y, weights = weights)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaKnownAgNormal")
## Values
value <- Values(array(c(0.2, 0.3, 0.4), dim = 3, dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgNormal(value = value, sd = sqrt(value))
spec <- Model(y ~ Normal(mean ~ age + region, sd = 0.2),
aggregate = aggregate)
x <- initialModel(spec, y = y, weights = weights)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaKnownAgNormal")
})
test_that("initialModel works with NormalVaryingVarsigmaUnknown and AgCertain", {
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20, mean = 0, sd = 0.3),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgCertain(value = 0.4)
spec <- Model(y ~ Normal(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, weights = weights)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaUnknownAgCertain")
## Values
value <- Values(array(c(0.2, 0.3, 0.4),
dim = 3,
dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgCertain(value = value)
spec <- Model(y ~ Normal(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, weights = weights)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaUnknownAgCertain")
})
test_that("initialModel works with NormalVaryingVarsigmaUnknown and AgNormal", {
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20, mean = 0, sd = 0.3),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgNormal(value = 0.4, sd = 0.1)
spec <- Model(y ~ Normal(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, weights = weights)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaUnknownAgNormal")
## Values
value <- Values(array(c(0.2, 0.3, 0.4), dim = 3, dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgNormal(value = value, sd = sqrt(value))
spec <- Model(y ~ Normal(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, weights = weights)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaUnknownAgNormal")
})
test_that("initialModel works with PoissonVaryingNotUseExp and AgCertain", {
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgCertain(value = 0.4)
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE),
upper = 3,
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = NULL)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingNotUseExpAgCertain")
## Values
value <- Values(array(c(0.2, 0.3, 0.4),
dim = 3,
dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgCertain(value = value)
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = NULL)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingNotUseExpAgCertain")
})
test_that("initialModel works with PoissonVaryingNotUseExp and AgNormal", {
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgNormal(value = 0.4, sd = 0.1)
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = NULL)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingNotUseExpAgNormal")
## Values
value <- Values(array(c(0.2, 0.3, 0.4), dim = 3, dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgNormal(value = value, sd = sqrt(value))
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = NULL)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingNotUseExpAgNormal")
})
test_that("initialModel works with PoissonVaryingUseExp and AgCertain", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.5 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgCertain(value = 0.4)
spec <- Model(y ~ Poisson(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingUseExpAgCertain")
## Values
value <- Values(array(c(0.2, 0.3, 0.4),
dim = 3,
dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgCertain(value = value)
spec <- Model(y ~ Poisson(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingUseExpAgCertain")
})
test_that("initialModel works with PoissonVaryingUseExp and AgNormal", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.5 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgNormal(value = 0.4, sd = 0.1)
spec <- Model(y ~ Poisson(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingUseExpAgNormal")
## Values
value <- Values(array(c(0.2, 0.3, 0.4), dim = 3, dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgNormal(value = value, sd = sqrt(value))
spec <- Model(y ~ Poisson(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingUseExpAgNormal")
})
test_that("initialModel works with PoissonVaryingNotUseExp and AgPoisson", {
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgPoisson(value = 0.4)
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = NULL)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingNotUseExpAgPoisson")
## Values
value <- Values(array(c(0.2, 0.3, 0.4), dim = 3, dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgPoisson(value = value)
spec <- Model(y ~ Poisson(mean ~ age + region, useExpose = FALSE),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = NULL)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingNotUseExpAgPoisson")
})
test_that("initialModel works with PoissonVaryingUseExp and AgPoisson", {
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.5 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
## scalar
aggregate <- AgPoisson(value = 0.4)
spec <- Model(y ~ Poisson(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingUseExpAgPoisson")
## Values
value <- Values(array(c(0.2, 0.3, 0.4), dim = 3, dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgPoisson(value = value)
spec <- Model(y ~ Poisson(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = exposure)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingUseExpAgPoisson")
})
test_that("initialModel works with PoissonVaryingUseExp and AgLife", {
initialModel <- demest:::initialModel
expose <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = c(0:3, "4+"), region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.5 * expose),
dim = c(5, 4),
dimnames = list(age = c(0:3, "4+"), region = letters[1:4])))
## scalar
aggregate <- AgLife(value = 20, sd = 0.1)
spec <- Model(y ~ Poisson(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = expose)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingUseExpAgLife")
## Values
value <- Values(array(c(0.2, 0.3, 0.4), dim = 3, dimnames = list(region = c("a", "b", "c"))))
aggregate <- AgNormal(value = value, sd = sqrt(value))
spec <- Model(y ~ Poisson(mean ~ age + region),
aggregate = aggregate)
x <- initialModel(spec, y = y, exposure = expose)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingUseExpAgNormal")
})
## initialModelPredict ###############################################################
test_that("initialModelPredict works with BinomialVarying - without aggregate", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(as.integer(runif(n = 20, min = 5, max = 100)),
dim = c(5, 4),
dimnames = list(time = 2001:2005, region = 1:4)),
dimscales = c(time = "Intervals"))
weights.bench <- Counts(array(as.integer(runif(n = 20, min = 5, max = 100)),
dim = c(4, 4),
dimnames = list(time = 2006:2009, region = 1:4)),
dimscales = c(time = "Intervals"))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.7),
dim = c(5, 4),
dimnames = list(time = 2001:2005, region = 1:4)),
dimscales = c(time = "Intervals"))
spec <- Model(y ~ Binomial(mean ~ time + region))
mod.est <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(mod.est,
along = 1L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = NULL,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "BinomialVaryingPredict")
})
test_that("initialModelPredict works with BinomialVaryingPredict - with aggregate", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(as.integer(runif(n = 20, min = 5, max = 100)),
dim = c(5, 4),
dimnames = list(time = 2001:2005, region = 1:4)),
dimscales = c(time = "Intervals"))
weights.bench <- Counts(array(as.integer(runif(n = 20, min = 5, max = 100)),
dim = c(4, 4),
dimnames = list(time = 2006:2009, region = 1:4)),
dimscales = c(time = "Intervals"))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.7),
dim = c(5, 4),
dimnames = list(time = 2001:2005, region = 1:4)),
dimscales = c(time = "Intervals"))
## AgCertain with weights
aggregate <- AgCertain(value = 0.5, weights = weights.bench)
spec <- Model(y ~ Binomial(mean ~ time + region))
mod.est <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(mod.est,
along = 1L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = aggregate,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "BinomialVaryingPredictAgCertain")
## AgNormal with weights
aggregate <- AgNormal(value = 0.5, sd = 1, weights = weights.bench)
spec <- Model(y ~ Binomial(mean ~ time + region))
mod.est <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(mod.est,
along = 1L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = aggregate,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "BinomialVaryingPredictAgNormal")
})
test_that("initialModelPredict works with NormalVaryingVarsigmaUnknown - without aggregate", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(n = 20, mean = 0, sd = 0.3),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
weights <- Counts(array(rbeta(n = 20, shape1 = 1, shape2 = 1),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Normal(mean ~ age + region))
x <- initialModel(spec, y = y, weights = weights)
set.seed(1)
x <- initialModelPredict(x,
along = 2L,
labels = c("a", "b", "c", "d"),
n = NULL,
offsetModel = 1L,
covariates = NULL,
aggregate = NULL,
lower = 1,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaUnknownPredict")
})
test_that("initialModelPredict works with NormalVaryingVarsigmaUnknownPredict - with aggregate", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
y <- Counts(array(rnorm(20),
dim = c(5, 4),
dimnames = list(time = 2001:2005, region = 1:4)),
dimscales = c(time = "Intervals"))
weights <- Counts(array(1,
dim = c(5, 4),
dimnames = list(time = 2001:2005, region = 1:4)),
dimscales = c(time = "Intervals"))
weights.ag <- Counts(array(1,
dim = c(4, 4),
dimnames = list(time = 2006:2009, region = 1:4)),
dimscales = c(time = "Intervals"))
## AgCertain
aggregate <- AgCertain(value = 0.5, weights = weights.ag)
spec <- Model(y ~ Normal(mean ~ time + region))
mod.est <- initialModel(spec, y = y, weights = weights)
x <- initialModelPredict(mod.est,
along = 1L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = aggregate,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaUnknownPredictAgCertain")
## AgNormal
aggregate <- AgNormal(value = 0.5, sd = 1, weights = weights.ag)
spec <- Model(y ~ Normal(mean ~ time + region))
mod.est <- initialModel(spec, y = y, weights = weights)
x <- initialModelPredict(mod.est,
along = 1L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = aggregate,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "NormalVaryingVarsigmaUnknownPredictAgNormal")
})
test_that("initialModelPredict works with PoissonVaryingUseExp", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.5 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Poisson(mean ~ age + region))
x <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(x,
along = 2L,
labels = c("a", "b", "c", "d"),
n = NULL,
offsetModel = 1L,
covariates = NULL,
aggregate = NULL,
lower = NULL,
upper = 100000)
expect_is(x, "PoissonVaryingUseExpPredict")
})
test_that("initialModelPredict works with PoissonVaryingNotUseExpPredict - with aggregate", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
y <- Counts(array(rpois(20, lambda = 100),
dim = c(5, 4),
dimnames = list(time = 2001:2005, region = 1:4)),
dimscales = c(time = "Intervals"))
weights.bench <- Counts(array(1,
dim = c(4, 4),
dimnames = list(time = 2006:2009, region = 1:4)),
dimscales = c(time = "Intervals"))
## AgCertain with no weights
aggregate <- AgCertain(value = 0.5, weights = weights.bench)
spec <- Model(y ~ Poisson(mean ~ time + region, useExpose = FALSE))
mod.est <- initialModel(spec, y = y, exposure = NULL)
x <- initialModelPredict(mod.est,
along = 1L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = aggregate,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingNotUseExpPredictAgCertain")
expect_true(all(x@weightAg == 1))
## AgNormal
aggregate <- AgNormal(value = 0.5, sd = 1, weights = weights.bench)
spec <- Model(y ~ Poisson(mean ~ time + region, useExpose = FALSE))
mod.est <- initialModel(spec, y = y, exposure = NULL)
x <- initialModelPredict(mod.est,
along = 1L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = aggregate,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "PoissonVaryingNotUseExpPredictAgNormal")
})
test_that("initialModelPredict works with CMPVaryingUseExp", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.5 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ CMP(mean ~ age + region))
x <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(x,
along = 2L,
labels = c("a", "b", "c", "d"),
n = NULL,
offsetModel = 1L,
covariates = NULL,
aggregate = NULL,
lower = NULL,
upper = 100000)
expect_is(x, "CMPVaryingUseExpPredict")
})
test_that("initialModelPredict works with CMPVaryingUseExp", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rpois(n = 20, lambda = 0.5 * exposure),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ CMP(mean ~ age + region))
x <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(x,
along = 2L,
labels = c("a", "b", "c", "d"),
n = NULL,
offsetModel = 1L,
covariates = NULL,
aggregate = NULL,
lower = NULL,
upper = 100000)
expect_is(x, "CMPVaryingUseExpPredict")
})
test_that("initialModelPredict works with PoissonBinomial", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ PoissonBinomial(prob = 0.95))
model <- initialModel(spec, y = y, exposure = exposure)
ans <- initialModelPredict(model,
along = 1L,
labels = NULL,
n = 5L,
offsetModel = 1L)
expect_true(validObject(ans))
metadata.expected <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 5:9, region = letters[1:4])))@metadata
expect_identical(ans@metadataY, metadata.expected)
expect_identical(ans@iMethodModel, model@iMethodModel + 100L)
})
test_that("initialModelPredict works with Round3", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- round3(y)
spec <- Model(y ~ Round3())
model <- initialModel(spec, y = y, exposure = exposure)
ans <- initialModelPredict(model,
along = 1L,
labels = NULL,
n = 5L,
offsetModel = 1L)
expect_true(validObject(ans))
metadata.expected <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 5:9, region = letters[1:4])))@metadata
expect_identical(ans@metadataY, metadata.expected)
expect_identical(ans@iMethodModel, model@iMethodModel + 100L)
})
test_that("initialModelPredict works with Exact", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
spec <- Model(y ~ Exact())
model <- initialModel(spec, y = y, exposure = exposure)
ans <- initialModelPredict(model,
along = 1L,
labels = NULL,
n = 5L,
offsetModel = 1L)
expect_true(validObject(ans))
metadata.expected <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 5:9, region = letters[1:4])))@metadata
expect_identical(ans@metadataY, metadata.expected)
expect_identical(ans@iMethodModel, model@iMethodModel + 100L)
})
test_that("initialModelPredict works with NormFixedUseExp", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
mean <- Values(array(rnorm(n = 40),
dim = c(10, 4),
dimnames = list(age = 0:9, region = letters[1:4])))
spec <- Model(y ~ NormalFixed(mean = mean, sd = sqrt(abs(mean))))
model <- initialModel(spec, y = y, exposure = exposure)
ans <- initialModelPredict(model,
along = 1L,
labels = NULL,
n = 5L,
offsetModel = 1L)
expect_true(validObject(ans))
metadata.expected <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 5:9, region = letters[1:4])))@metadata
expect_identical(ans@metadataY, metadata.expected)
expect_identical(ans@iMethodModel, model@iMethodModel + 100L)
})
test_that("initialModelPredict works with NormFixedNotUseExp", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
mean <- Values(array(rnorm(n = 40),
dim = c(10, 4),
dimnames = list(age = 0:9, region = letters[1:4])))
spec <- Model(y ~ NormalFixed(mean = mean, sd = sqrt(abs(mean)), useExpose = FALSE))
model <- initialModel(spec, y = y, exposure = NULL)
ans <- initialModelPredict(model,
along = 1L,
labels = NULL,
n = 5L,
offsetModel = 1L)
expect_true(validObject(ans))
metadata.expected <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 5:9, region = letters[1:4])))@metadata
expect_identical(ans@metadataY, metadata.expected)
expect_identical(ans@iMethodModel, model@iMethodModel + 100L)
})
test_that("initialModelPredict works with NormFixedUseExp", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
exposure <- Counts(array(rpois(n = 20, lambda = 20),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
y <- Counts(array(rbinom(n = 20, size = exposure, prob = 0.5),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
location <- Values(array(rnorm(n = 40),
dim = c(10, 4),
dimnames = list(age = 0:9, region = letters[1:4])))
spec <- Model(y ~ TFixed(location = location, scale = sqrt(abs(location))))
model <- initialModel(spec, y = y, exposure = exposure)
ans <- initialModelPredict(model,
along = 1L,
labels = NULL,
n = 5L,
offsetModel = 1L)
expect_true(validObject(ans))
metadata.expected <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 5:9, region = letters[1:4])))@metadata
expect_identical(ans@metadataY, metadata.expected)
expect_identical(ans@iMethodModel, model@iMethodModel + 100L)
})
test_that("initialModelPredict works with NormFixedNotUseExp", {
initialModelPredict <- demest:::initialModelPredict
initialModel <- demest:::initialModel
y <- Counts(array(rpois(n = 20, lambda = 10),
dim = c(5, 4),
dimnames = list(age = 0:4, region = letters[1:4])))
location <- Values(array(rnorm(n = 40),
dim = c(10, 4),
dimnames = list(age = 0:9, region = letters[1:4])))
spec <- Model(y ~ TFixed(location = location, scale = sqrt(abs(location)), useExpose = FALSE))
model <- initialModel(spec, y = y, exposure = NULL)
ans <- initialModelPredict(model,
along = 1L,
labels = NULL,
n = 5L,
offsetModel = 1L)
expect_true(validObject(ans))
metadata.expected <- Counts(array(1L,
dim = c(5, 4),
dimnames = list(age = 5:9, region = letters[1:4])))@metadata
expect_identical(ans@metadataY, metadata.expected)
expect_identical(ans@iMethodModel, model@iMethodModel + 100L)
})
test_that("initialModelPredict creates object of class LN2Predict from valid inputs", {
initialModel <- demest:::initialModel
initialModelPredict <- demest:::initialModelPredict
## 'constraint' does not have time dimension
constraint <- Values(array(c(NA, -1L, 0L, 1L),
dim = c(2, 2),
dimnames = list(age = c("0-39", "40+"),
sex = c("Female", "Male"))))
y <- Counts(array(10L,
dim = c(2, 4, 3),
dimnames = c(list(sex = c("Female", "Male"),
age = c("0-19", "20-39", "40-59", "60+"),
time = c("2000", "2010", "2020")))))
exposure <- 2L * y
spec <- Model(y ~ LN2(constraint = constraint))
mod.est <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(mod.est,
along = 3L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = NULL,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "LN2Predict")
## 'constraint' has time dimension
constraint <- Values(array(c(NA, -1L),
dim = c(2, 7),
dimnames = list(age = c("0-39", "40+"),
time = seq(2000, 2060, 10))))
y <- Counts(array(10L,
dim = c(2, 4, 3),
dimnames = c(list(sex = c("Female", "Male"),
age = c("0-19", "20-39", "40-59", "60+"),
time = c("2000", "2010", "2020")))))
exposure <- 2L * y
spec <- Model(y ~ LN2(constraint = constraint))
mod.est <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(mod.est,
along = 3L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = NULL,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "LN2Predict")
expect_identical(as.character(limits(x@constraintLN2)$time), c("2030", "2060"))
expect_true(x@add1@.Data)
spec <- Model(y ~ LN2(constraint = constraint, add1 = FALSE))
mod.est <- initialModel(spec, y = y, exposure = exposure)
x <- initialModelPredict(mod.est,
along = 3L,
labels = NULL,
n = 4L,
offsetModel = 1L,
covariates = NULL,
aggregate = NULL,
lower = NULL,
upper = NULL)
expect_true(validObject(x))
expect_is(x, "LN2Predict")
expect_false(x@add1@.Data)
})
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.