Skip to content

Commit 644bcc4

Browse files
johnrbryantjohnrbryant
committed
tests passing
1 parent d17ad3f commit 644bcc4

24 files changed

+1040
-128
lines changed

‎R/AllClasses-14-Model-hierarchy-spec.R‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -187,7 +187,8 @@ setClass("SpecLikelihoodLN2",
187187
"SpecAVarsigmaMixin",
188188
"SpecVarsigmaMaxMixin",
189189
"StructuralZerosMixin"),
190-
slots = c(constraintLN2 = "Values",
190+
slots = c(add1 = "LogicalFlag",
191+
constraintLN2 = "Values",
191192
concordances = "list",
192193
updateVarsigmaLN2 = "LogicalFlag",
193194
varsigmaLN2 = "Scale",
@@ -478,7 +479,8 @@ setClass("SpecLN2",
478479
"SpecSeriesMixin",
479480
"StructuralZerosMixin"),
480481
prototype = prototype(useExpose = methods::new("LogicalFlag", TRUE)),
481-
slots = c(constraintLN2 = "Values",
482+
slots = c(add1 = "LogicalFlag",
483+
constraintLN2 = "Values",
482484
concordances = "list",
483485
updateVarsigmaLN2 = "LogicalFlag",
484486
varsigmaLN2 = "Scale",

‎R/AllClasses-15-Model-hierarchy.R‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -402,7 +402,8 @@ setClass("TFixedUseExp",
402402

403403
## NO_TESTS
404404
setClass("LN2",
405-
slots = c(alphaLN2 = "ParameterVector",
405+
slots = c(add1 = "LogicalFlag",
406+
alphaLN2 = "ParameterVector",
406407
constraintLN2 = "Values",
407408
constraintAllLN2 = "Values",
408409
nCellBeforeLN2 = "integer",

‎R/AllClasses-19-Skeleton.R‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -608,7 +608,8 @@ setClass("SkeletonMissingDatasetNormalFixedUseExp",
608608
## HAS_TESTS
609609
## HAS_FETCH
610610
setClass("SkeletonMissingDatasetLN2",
611-
slots = c(offsetsAlphaLN2 = "Offsets",
611+
slots = c(add1 = "LogicalFlag",
612+
offsetsAlphaLN2 = "Offsets",
612613
offsetsVarsigmaLN2 = "Offsets",
613614
transformLN2 = "CollapseTransformExtra",
614615
updateVarsigmaLN2 = "LogicalFlag",

‎R/Model-generators.R‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1526,6 +1526,7 @@ setMethod("initialModel",
15261526
function(object, y, exposure) {
15271527
A.varsigma <- object@AVarsigma@.Data
15281528
A.sigma <- object@ASigma@.Data
1529+
add1 <- object@add1
15291530
call <- object@call
15301531
concordances <- object@concordances
15311532
constraint.all <- object@constraintLN2
@@ -1600,7 +1601,10 @@ setMethod("initialModel",
16001601
sigma.max <- methods::new("Scale", sigma.max)
16011602
## randomly draw 'alpha', using residuals and constraint
16021603
alpha <- numeric(length = length(constraint))
1603-
resid <- log1p(y) - log1p(exposure)
1604+
if (add1@.Data)
1605+
resid <- log1p(y) - log1p(exposure)
1606+
else
1607+
resid <- log(y) - log(exposure)
16041608
resid[is.na(resid)] <- 0
16051609
resid <- dembase::collapse(resid,
16061610
transform = transform)
@@ -1630,6 +1634,7 @@ setMethod("initialModel",
16301634
n.cell.before <- rep(0L, times = length(alpha)) # temporary value
16311635
cellInLik <- rep(TRUE, times = length(y)) # temporary value
16321636
model <- methods::new("LN2",
1637+
add1 = add1,
16331638
alphaLN2 = alpha,
16341639
ASigma = A.sigma,
16351640
AVarsigma = A.varsigma,

‎R/Skeleton-generator.R‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -371,6 +371,7 @@ setMethod("SkeletonMissingDataset",
371371
offsets.component <- methods::new("Offsets",
372372
c(skeletonComponent@first, skeletonComponent@last))
373373
methods::new("SkeletonMissingDatasetLN2",
374+
add1 = model@add1,
374375
offsetsAlphaLN2 = offsets.alpha,
375376
offsetsVarsigmaLN2 = offsets.varsigma,
376377
offsetsComponent = offsets.component,

‎R/Skeleton-methods.R‎

Lines changed: 70 additions & 69 deletions
Original file line numberDiff line numberDiff line change
@@ -942,78 +942,79 @@ setMethod("fetchResults",
942942
function(object, nameObject, filename, iterations,
943943
nIteration, lengthIter,
944944
impute = FALSE) {
945-
data <- object@data
946-
if (impute) {
947-
offsets.alpha <- object@offsetsAlphaLN2
948-
offsets.varsigma <- object@offsetsVarsigmaLN2
949-
offsets.comp <- object@offsetsComponent
950-
transform.alpha <- object@transformLN2
951-
transform.comp <- object@transformComponent
952-
struc.zero.array <- object@strucZeroArray
953-
update.varsigma <- object@updateVarsigmaLN2@.Data
954-
if (is.null(iterations))
955-
iterations <- seq_len(nIteration)
956-
metadata <- data@metadata
957-
n.data <- prod(dim(metadata))
958-
metadata <- dembase::addIterationsToMetadata(metadata,
959-
iterations = iterations)
960-
n.iter <- length(iterations)
961-
transform.alpha <- addIterationsToTransform(transform.alpha,
962-
nIter = n.iter)
963-
transform.comp <- addIterationsToTransform(transform.comp,
964-
nIter = n.iter)
965-
transform.alpha <- methods::new("ExtendTransform",
966-
indices = transform.alpha@indices,
967-
dims = transform.alpha@dims,
968-
dimBefore = transform.alpha@dimAfter,
969-
dimAfter = transform.alpha@dimBefore)
970-
.Data <- array(data@.Data,
971-
dim = dim(metadata),
972-
dimnames = dimnames(metadata))
973-
alpha <- getDataFromFile(filename = filename,
974-
first = offsets.alpha[1L],
975-
last = offsets.alpha[2L],
976-
lengthIter = lengthIter,
977-
iterations = iterations)
978-
if (update.varsigma)
979-
varsigma <- getDataFromFile(filename = filename,
980-
first = offsets.varsigma[1L],
981-
last = offsets.varsigma[2L],
945+
data <- object@data
946+
if (impute) {
947+
add1 <- object@add1@.Data
948+
offsets.alpha <- object@offsetsAlphaLN2
949+
offsets.varsigma <- object@offsetsVarsigmaLN2
950+
offsets.comp <- object@offsetsComponent
951+
transform.alpha <- object@transformLN2
952+
transform.comp <- object@transformComponent
953+
struc.zero.array <- object@strucZeroArray
954+
update.varsigma <- object@updateVarsigmaLN2@.Data
955+
if (is.null(iterations))
956+
iterations <- seq_len(nIteration)
957+
metadata <- data@metadata
958+
n.data <- prod(dim(metadata))
959+
metadata <- dembase::addIterationsToMetadata(metadata,
960+
iterations = iterations)
961+
n.iter <- length(iterations)
962+
transform.alpha <- addIterationsToTransform(transform.alpha,
963+
nIter = n.iter)
964+
transform.comp <- addIterationsToTransform(transform.comp,
965+
nIter = n.iter)
966+
transform.alpha <- methods::new("ExtendTransform",
967+
indices = transform.alpha@indices,
968+
dims = transform.alpha@dims,
969+
dimBefore = transform.alpha@dimAfter,
970+
dimAfter = transform.alpha@dimBefore)
971+
.Data <- array(data@.Data,
972+
dim = dim(metadata),
973+
dimnames = dimnames(metadata))
974+
alpha <- getDataFromFile(filename = filename,
975+
first = offsets.alpha[1L],
976+
last = offsets.alpha[2L],
977+
lengthIter = lengthIter,
978+
iterations = iterations)
979+
if (update.varsigma)
980+
varsigma <- getDataFromFile(filename = filename,
981+
first = offsets.varsigma[1L],
982+
last = offsets.varsigma[2L],
983+
lengthIter = lengthIter,
984+
iterations = iterations)
985+
else
986+
varsigma <- rep(object@varsigma@.Data,
987+
times = n.iter)
988+
exposure <- getDataFromFile(filename = filename,
989+
first = offsets.comp[1L],
990+
last = offsets.comp[2L],
982991
lengthIter = lengthIter,
983992
iterations = iterations)
993+
alpha <- array(alpha,
994+
dim = transform.alpha@dimBefore)
995+
varsigma <- rep(varsigma, each = n.data)
996+
exposure <- array(exposure,
997+
dim = transform.comp@dimBefore)
998+
alpha <- dembase::extend(alpha,
999+
transform = transform.alpha)
1000+
exposure <- dembase::collapse(exposure,
1001+
transform = transform.comp)
1002+
is.struc.zero <- struc.zero.array@.Data == 0L
1003+
is.struc.zero <- array(is.struc.zero,
1004+
dim = dim(.Data))
1005+
is.missing <- is.na(.Data)
1006+
n.impute <- sum(is.missing)
1007+
mean <- log(exposure[is.missing] + add1) + alpha[is.missing]
1008+
sd <- varsigma[is.missing]
1009+
imputed <- exp(stats::rnorm(n = n.impute, mean = mean, sd = sd)) - add1
1010+
.Data[is.missing] <- imputed
1011+
.Data[is.struc.zero] <- 0L
1012+
methods::new("Counts",
1013+
.Data = .Data,
1014+
metadata = metadata)
1015+
}
9841016
else
985-
varsigma <- rep(object@varsigma@.Data,
986-
times = n.iter)
987-
exposure <- getDataFromFile(filename = filename,
988-
first = offsets.comp[1L],
989-
last = offsets.comp[2L],
990-
lengthIter = lengthIter,
991-
iterations = iterations)
992-
alpha <- array(alpha,
993-
dim = transform.alpha@dimBefore)
994-
varsigma <- rep(varsigma, each = n.data)
995-
exposure <- array(exposure,
996-
dim = transform.comp@dimBefore)
997-
alpha <- dembase::extend(alpha,
998-
transform = transform.alpha)
999-
exposure <- dembase::collapse(exposure,
1000-
transform = transform.comp)
1001-
is.struc.zero <- struc.zero.array@.Data == 0L
1002-
is.struc.zero <- array(is.struc.zero,
1003-
dim = dim(.Data))
1004-
is.missing <- is.na(.Data)
1005-
n.impute <- sum(is.missing)
1006-
mean <- log(exposure[is.missing] + 1) + alpha[is.missing]
1007-
sd <- varsigma[is.missing]
1008-
imputed <- exp(stats::rnorm(n = n.impute, mean = mean, sd = sd)) - 1
1009-
.Data[is.missing] <- imputed
1010-
.Data[is.struc.zero] <- 0L
1011-
methods::new("Counts",
1012-
.Data = .Data,
1013-
metadata = metadata)
1014-
}
1015-
else
1016-
data
1017+
data
10171018
})
10181019

10191020

‎R/SpecModel-generators.R‎

Lines changed: 34 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -240,13 +240,22 @@ Binomial <- function(formula, structuralZeros = NULL) {
240240
## HAS_TESTS
241241
#' Log-normal model with two levels
242242
#'
243-
#' THIS FUNCTION IS STILL EXPERIMENTAL.
243+
#' A hierarchical model for datasets with systematic
244+
#' biases, where these biases are modelled as proportional
245+
#' differences.
244246
#'
245-
#' A simple hierarchical model in which the log of y + 1
246-
#' equals the log of exposure + 1, plus a bias term.
247-
#' The bias term can be subject to constraints.
247+
#' When \code{add1} is \code{TRUE} (the default), the
248+
#' model assumes
249+
#'
250+
#' \deqn{y[i] + 1 \sim N(exposure[i] + 1 + alpha[j[i]], sd^2)}
251+
#'
252+
#' and when \code{add1} is \code{FALSE}, the model
253+
#' assumes
248254
#'
249-
#' The bias term has the same level of detail as the constraints,
255+
#' \deqn{y[i] \sim N(exposure[i] + alpha[j[i]], sd^2)}.
256+
#'
257+
#' The \code{alpha} term measures bias. It has
258+
#' the same level of detail as the constraints,
250259
#' which will almost always be less than the level of detail
251260
#' of \code{y}.
252261
#'
@@ -257,6 +266,9 @@ Binomial <- function(formula, structuralZeros = NULL) {
257266
#' implies that it must be positive, and a NA
258267
#' implies that there is no constraint.
259268
#'
269+
#' Adding 1 to the response and exposure is a way of dealing with
270+
#' zeros in the response.
271+
#'
260272
#' @inheritParams likelihood
261273
#' @param constraint An object of class \code{\link[dembase:Values-class]{Values}},
262274
#' with values 0, -1, 1, and NA.
@@ -266,10 +278,18 @@ Binomial <- function(formula, structuralZeros = NULL) {
266278
#' fixed value or a prior distribution. The prior distribution
267279
#' is specified via \code{\link{HalfT}} or \code{\link{InvChiSq}}.
268280
#' Defaults to a half-t distribution.
269-
#'
281+
#' @param add1 Whether to add 1 to the response and exposure.
282+
#'
283+
#' @examples
284+
#' constraint <- ValuesOne(c(NA, NA, -1),
285+
#' labels = c("0-19", "20-64", "65"),
286+
#' name = "age")
287+
#' spec <- LN2(constraint)
288+
#' spec
270289
#' @export
271290
LN2 <- function(constraint, structuralZeros = NULL,
272-
concordances = list(), sd = HalfT()) {
291+
concordances = list(), sd = HalfT(),
292+
add1 = TRUE) {
273293
## constraint
274294
if (!methods::is(constraint, "Values"))
275295
stop(gettextf("'%s' has class \"%s\"",
@@ -328,9 +348,14 @@ LN2 <- function(constraint, structuralZeros = NULL,
328348
varsigmaLN2 <- methods::new("Scale", varsigmaLN2)
329349
updateVarsigmaLN2 <- methods::new("LogicalFlag", updateVarsigmaLN2)
330350
varsigmaLN2HasHalfT <- methods::new("LogicalFlag", sd.is.half.t)
351+
## add1
352+
checkLogical(x = add1,
353+
name = "add1")
354+
add1 <- new("LogicalFlag", add1)
331355
## return
332356
methods::new("SpecLikelihoodLN2",
333357
AVarsigma = AVarsigma,
358+
add1 = add1,
334359
concordances = concordances,
335360
multVarsigma = multVarsigma,
336361
constraintLN2 = constraint,
@@ -1357,6 +1382,7 @@ setMethod("SpecModel",
13571382
priorSD, jump,
13581383
series, aggregate) {
13591384
A.varsigma <- specInner@AVarsigma
1385+
add1 <- specInner@add1
13601386
concordances <- specInner@concordances
13611387
constraintLN2 <- specInner@constraintLN2
13621388
mult.varsigma <- specInner@multVarsigma
@@ -1413,6 +1439,7 @@ setMethod("SpecModel",
14131439
methods::new("SpecLN2",
14141440
ASigma = A.sigma,
14151441
AVarsigma = A.varsigma,
1442+
add1 = add1,
14161443
call = call,
14171444
concordances = concordances,
14181445
multSigma = mult.sigma,

‎R/helper-functions.R‎

Lines changed: 9 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2313,10 +2313,17 @@ logLikelihood_LN2 <- function(model, count, dataset, i, useC = FALSE) {
23132313
alpha <- model@alphaLN2@.Data
23142314
transform <- model@transformLN2
23152315
sd <- model@varsigma@.Data
2316-
x <- log1p(dataset[[i]])
2316+
add1 <- model@add1@.Data
23172317
j <- dembase::getIAfter(i = i,
23182318
transform = transform)
2319-
mean <- log1p(count) + alpha[j]
2319+
if (add1) {
2320+
x <- log1p(dataset[[i]])
2321+
mean <- log1p(count) + alpha[j]
2322+
}
2323+
else {
2324+
x <- log(dataset[[i]])
2325+
mean <- log(count) + alpha[j]
2326+
}
23202327
stats::dnorm(x = x,
23212328
mean = mean,
23222329
sd = sd,

‎R/helper-printing.R‎

Lines changed: 23 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -360,7 +360,11 @@ printLN2LikEqns <- function(object) {
360360
nuVarsigma <- object@nuVarsigma@.Data
361361
AVarsigma <- object@AVarsigma@.Data
362362
varsigmaMax <- object@varsigmaMax@.Data
363-
cat(" log(y[i] + 1) ~ N(log(exposure[i] + 1) + alpha[j[i]], ")
363+
add1 <- object@add1@.Data
364+
if (add1)
365+
cat(" log(y[i] + 1) ~ N(log(exposure[i] + 1) + alpha[j[i]], ")
366+
else
367+
cat(" log(y[i]) ~ N(log(exposure[i]) + alpha[j[i]], ")
364368
if (updateVarsigma)
365369
cat("sdData^2)\n")
366370
else
@@ -390,6 +394,7 @@ printLN2ModEqns <- function(object) {
390394
AVarsigma <- object@AVarsigma@.Data
391395
varsigmaMax <- object@varsigmaMax
392396
nuVarsigma <- object@nuVarsigma
397+
add1 <- object@add1@.Data
393398
name.y <- deparse(call$formula[[2L]])
394399
if (is.null(series)) {
395400
if (identical(name.y, "y"))
@@ -401,8 +406,14 @@ printLN2ModEqns <- function(object) {
401406
exposure <- series
402407
n.spaces <- max(5L - nchar(name.y), 0L)
403408
spaces <- rep(" ", n.spaces)
404-
response <- sprintf("%slog(%s + 1)", spaces, name.y)
405-
cat(name.y, response, " ~ N(log(", exposure, "[i] + 1) + alpha[j[i]], ", sep = "")
409+
if (add1)
410+
response <- sprintf("%slog(%s + 1)", spaces, name.y)
411+
else
412+
response <- sprintf("%slog(%s)", spaces, name.y)
413+
if (add1)
414+
cat(name.y, response, " ~ N(log(", exposure, "[i] + 1) + alpha[j[i]], ", sep = "")
415+
else
416+
cat(name.y, response, " ~ N(log(", exposure, "[i]) + alpha[j[i]], ", sep = "")
406417
if (updateVarsigma)
407418
cat("sdData^2)\n")
408419
else
@@ -447,6 +458,7 @@ printLN2SpecEqns <- function(object) {
447458
varsigmaMax <- object@varsigmaMax
448459
nuVarsigma <- object@nuVarsigma
449460
nameY <- object@nameY
461+
add1 <- object@add1@.Data
450462
has.series <- !is.na(series)
451463
name.y <- deparse(call$formula[[2L]])
452464
if (has.series)
@@ -455,8 +467,14 @@ printLN2SpecEqns <- function(object) {
455467
exposure <- "exposure"
456468
n.spaces <- max(5L - nchar(name.y), 0L)
457469
spaces <- paste(rep(" ", n.spaces), collapse = "")
458-
response <- sprintf("%slog(%s + 1)", spaces, name.y)
459-
cat(response, " ~ N(log(", exposure, "[i] + 1) + alpha[j[i]], ", sep = "")
470+
if (add1)
471+
response <- sprintf("%slog(%s + 1)", spaces, name.y)
472+
else
473+
response <- sprintf("%slog(%s)", spaces, name.y)
474+
if (add1)
475+
cat(response, " ~ N(log(", exposure, "[i] + 1) + alpha[j[i]], ", sep = "")
476+
else
477+
cat(response, " ~ N(log(", exposure, "[i]) + alpha[j[i]], ", sep = "")
460478
if (updateVarsigma)
461479
cat("sdData^2)\n")
462480
else

0 commit comments

Comments
 (0)