Skip to content

Commit a65d77a

Browse files
kenkellnerkenkellner
authored
LINPRED and LINPRED_PRIORS automatically generate inits (#31)
* LINPRED and LINPRED_PRIORS auto-generate default inits * Add tests and fix some bugs
1 parent 76e50bc commit a65d77a

File tree

2 files changed

+189
-7
lines changed

2 files changed

+189
-7
lines changed

‎R/LINPRED.R‎

Lines changed: 82 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -412,11 +412,75 @@ fixTerms <- function(trms, formula_info){
412412
check_terms <- strsplit(trms, ":")
413413
reorder_terms <- sapply(check_terms, function(x, ref){
414414
ind <- which(sapply(ref, function(z) identical(sort(z), sort(x))))
415-
formula_info$terms[ind]
415+
if(length(ind) > 0) return(formula_info$terms[ind])
416+
paste(x, collapse=":") # return original if not found
416417
}, ref=formula_info$terms_ref)
417418
reorder_terms
418419
}
419420

421+
# Here 'formula' should be the output of processFormula
422+
makeDefaultCoefficientInits <- function(form_info, orig_formula,
423+
data, coefPrefix, modMatNames){
424+
425+
# Get structure of coefficients
426+
inits <- makeEmptyParameterStructure(removeBracketsFromFormula(form_info$formula), data)
427+
# Update names to match code
428+
names(inits) <- paste0(safeDeparse(coefPrefix), names(inits))
429+
names(inits) <- gsub(":", "_", names(inits), fixed=TRUE)
430+
431+
# Repeat with original formula to account for centerVar
432+
new_form <- removeBracketsFromFormula(as.formula(orig_formula))
433+
new_form <- reformulas::nobars(new_form)
434+
# make sure structure matches
435+
trms <- stats::terms(new_form)
436+
has_int <- ifelse(attr(trms, "intercept") == 1, "1", "0")
437+
trms <- attr(trms, "term.labels")
438+
trms <- c(has_int, trms)
439+
trms <- fixTerms(trms, form_info)
440+
new_form <- stats::reformulate(trms)
441+
442+
inits2 <- makeEmptyParameterStructure(new_form, data)
443+
if(length(inits2) > 0){
444+
# Update names to match code
445+
names(inits2) <- paste0(safeDeparse(coefPrefix), names(inits2))
446+
names(inits2) <- gsub(":", "_", names(inits2), fixed=TRUE)
447+
if(is.null(inits)) inits <- list()
448+
inits <- utils::modifyList(inits, inits2)
449+
}
450+
451+
# Cleanup structure
452+
inits <- lapply(inits, function(x){
453+
if(length(x) == 0){
454+
return(0)
455+
}
456+
drop(x)
457+
})
458+
459+
# If modMatNames, we also need inits for the re-named parameters
460+
if(modMatNames){
461+
nms <- makeParameterStructureModMatNames(new_form, data)
462+
nms <- unlist(nms)
463+
nms <- nms[nms != "0"]
464+
nms <- paste0(coefPrefix, nms)
465+
modmat_inits <- as.list(rep(0, length(nms)))
466+
names(modmat_inits) <- nms
467+
inits <- utils::modifyList(inits, modmat_inits)
468+
}
469+
inits
470+
}
471+
472+
makeDefaultSDInits <- function(formula, modelInfo, formula_info, sdPrefix){
473+
bars <- reformulas::findbars(formula)
474+
sd_names <- lapply(bars, getHyperpriorNames, modelInfo=modelInfo,
475+
formula_info=formula_info, prefix=sdPrefix)
476+
sd_names <- unlist(sd_names)
477+
sd_names <- sapply(sd_names, safeDeparse)
478+
sd_names
479+
inits <- as.list(rep(1, length(sd_names)))
480+
names(inits) <- sd_names
481+
inits
482+
}
483+
420484
#' Macro to build code for linear predictor from R formula
421485
#'
422486
#' Converts an R formula into corresponding code for a linear predictor in NIMBLE model code.
@@ -496,6 +560,13 @@ function(stoch, LHS, formula, link=NULL, coefPrefix=quote(beta_),
496560
RHS <- as.call(list(quote(nimbleMacros::FORLOOP), RHS))
497561
# Combine LHS and RHS
498562
code <- substitute(LHS <- RHS, list(LHS = LHS, RHS = RHS))
563+
# Add default inits
564+
inits <- makeDefaultCoefficientInits(form_info, formula, dat, coefPrefix,
565+
modMatNames)
566+
if(length(inits) > 0){
567+
if(is.null(modelInfo$inits)) modelInfo$inits <- list()
568+
modelInfo$inits <- utils::modifyList(modelInfo$inits, inits)
569+
}
499570

500571
# Add code for priors to output if needed
501572
if(!is.null(priorSpecs)){
@@ -679,6 +750,16 @@ function(form, coefPrefix=quote(beta_), sdPrefix=NULL, priorSpecs=setPriors(),
679750
out <- out[!sapply(out, is.null)]
680751
out <- embedLinesInCurlyBrackets(out)
681752
out <- removeExtraBrackets(out)
753+
754+
# Add default inits
755+
inits <- makeDefaultCoefficientInits(form_info, form, dat, coefPrefix,
756+
modMatNames)
757+
sd_inits <- makeDefaultSDInits(form, modelInfo, form_info, sdPrefix)
758+
inits <- utils::modifyList(inits, sd_inits)
759+
if(length(inits) > 0){
760+
if(is.null(modelInfo$inits)) modelInfo$inits <- list()
761+
modelInfo$inits <- utils::modifyList(modelInfo$inits, inits)
762+
}
682763

683764
list(code=out, modelInfo=modelInfo)
684765
},

‎tests/testthat/test_LINPRED.R‎

Lines changed: 107 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -385,7 +385,8 @@ test_that("LINPRED", {
385385
set.seed(123)
386386
modInfo <- list(constants=list(y = rnorm(10), x=factor(sample(letters[1:3], 10, replace=T)),
387387
x2=factor(sample(letters[4:5], 10, replace=T)),
388-
x3=round(rnorm(10),3)))
388+
x3=round(rnorm(10),3)),
389+
inits=list(beta_Intercept=0))
389390

390391
code <- quote(y[1:n] <- LINPRED(~1, priorSpecs=NULL))
391392

@@ -400,10 +401,17 @@ test_that("LINPRED", {
400401
)
401402

402403
code <- quote(y[1:n] ~ LINPRED(~x + x3, priorSpecs=NULL))
404+
out <- LINPRED$process(code, modInfo, NULL)
403405
expect_equal(
404-
LINPRED$process(code, modInfo, NULL)$code,
406+
out$code,
405407
quote(y[1:n] <- nimbleMacros::FORLOOP(beta_Intercept + beta_x[x[1:n]] + beta_x3 * x3[1:n]))
406408
)
409+
expect_equal(
410+
out$modelInfo$inits,
411+
list(beta_Intercept = 0,
412+
beta_x = structure(c(xa = 0, xb = 0, xc = 0), dim = 3L, dimnames = list(c("xa", "xb", "xc"))),
413+
beta_x3 = c(x3 = 0))
414+
)
407415

408416
code <- quote(y[1:n] ~ LINPRED(~x, priorSpecs=NULL, coefPrefix=alpha_))
409417
expect_equal(
@@ -438,18 +446,32 @@ test_that("LINPRED", {
438446
)
439447

440448
code <- quote(y[1:n] ~ LINPRED(~x + (x3|x2), priorSpecs=NULL))
449+
out <- LINPRED$process(code, modInfo, NULL)
441450
expect_equal(
442-
LINPRED$process(code, modInfo, NULL)$code,
451+
out$code,
443452
quote(y[1:n] <- nimbleMacros::FORLOOP(beta_Intercept + beta_x[x[1:n]] +
444453
beta_x2[x2[1:n]] + beta_x2_x3[x2[1:n]] * x3[1:n]))
445454
)
455+
expect_equal(
456+
out$modelInfo$inits,
457+
list(beta_Intercept = 0, beta_x = structure(c(xa = 0, xb = 0, xc = 0), dim = 3L, dimnames = list(c("xa", "xb", "xc"))),
458+
beta_x2 = structure(c(x2d = 0, x2e = 0), dim = 2L, dimnames = list(c("x2d", "x2e"))),
459+
beta_x2_x3 = c(x2d = 0, x2e = 0))
460+
)
446461

447462
# With modMatNames = TRUE (code should be unchanged)
448463
code <- quote(y[1:n] <- LINPRED(~x2, link=log, priorSpecs=NULL, modMatNames=TRUE))
464+
out <- LINPRED$process(code, modInfo, NULL)
449465
expect_equal(
450-
LINPRED$process(code, modInfo, NULL)$code,
466+
out$code,
451467
quote(log(y[1:n]) <- nimbleMacros::FORLOOP(beta_Intercept + beta_x2[x2[1:n]]))
452468
)
469+
expect_equal(
470+
out$modelInfo$inits,
471+
list(beta_Intercept = 0,
472+
beta_x2 = structure(c(x2d = 0, x2e = 0), dim = 2L, dimnames = list(c("x2d", "x2e"))),
473+
beta_x2e = 0)
474+
)
453475

454476
code2 <- quote(y[1:n] <- LINPRED(~x2, link=log, priorSpecs=NULL, modMatNames=FALSE))
455477
expect_equal(
@@ -479,6 +501,11 @@ test_that("LINPRED with random effect", {
479501
out$modelInfo$constants,
480502
modInfo$constants
481503
)
504+
expect_equal(
505+
out$modelInfo$inits,
506+
list(beta_Intercept = 0, beta_x3 = c(x3 = 0),
507+
beta_x = structure(c(xa = 0, xb = 0, xc = 0), dim = 3L, dimnames = list(c("xa", "xb", "xc"))))
508+
)
482509

483510
# With subtracted intercepts
484511
code <- quote(y[1:n] ~ LINPRED(~-1 + x, priorSpecs=NULL))
@@ -487,17 +514,26 @@ test_that("LINPRED with random effect", {
487514
out$code,
488515
quote(y[1:n] <- nimbleMacros::FORLOOP(beta_x[x[1:n]]))
489516
)
517+
expect_equal(
518+
out$modelInfo$inits,
519+
list(beta_x = structure(c(xa = 0, xb = 0, xc = 0), dim = 3L, dimnames = list(c("xa", "xb", "xc"))))
520+
)
490521

491522
code2 <- quote(y[1:n] ~ LINPRED(~-1 + (1|x), priorSpecs=NULL))
492523
out2 <- LINPRED$process(code2, modInfo, NULL)
493524
expect_equal(out2$code, out$code)
525+
expect_equal(out2$modelInfo$inits, out$modelInfo$inits)
494526

495527
code3 <- quote(y[1:n] ~ LINPRED(~-1 + (-1+x3|x), priorSpecs=NULL))
496528
out3 <- LINPRED$process(code3, modInfo, NULL)
497529
expect_equal(
498530
out3$code,
499531
quote(y[1:n] <- nimbleMacros::FORLOOP(beta_x3_x[x[1:n]] * x3[1:n]))
500532
)
533+
expect_equal(
534+
out3$modelInfo$inits,
535+
list(beta_x3_x = c(xa = 0, xb = 0, xc = 0))
536+
)
501537

502538

503539
code3 <- quote(LINPRED_PRIORS(~-1 + (-1+x3|x), priorSpecs=setPriors()))
@@ -539,6 +575,13 @@ test_that("LINPRED with random effect", {
539575
beta_x2_x[x2[1:n], x[1:n]]))
540576
)
541577

578+
expect_equal(
579+
out5$modelInfo$inits,
580+
list(beta_x2 = structure(c(x2d = 0, x2e = 0), dim = 2L, dimnames = list(c("x2d", "x2e"))),
581+
beta_x = structure(c(xa = 0, xb = 0, xc = 0), dim = 3L, dimnames = list(c("xa", "xb", "xc"))),
582+
beta_x2_x = structure(c(0, 0, 0, 0, 0, 0), dim = 2:3, dimnames = list(c("x2d", "x2e"), c("xa", "xb", "xc"))))
583+
)
584+
542585
code6 <- quote(LINPRED_PRIORS(~x*x2 - 1 - x + (1|x) ))
543586
out6 <- LINPRED_PRIORS$process(code6, modInfo, NULL)
544587
expect_equal(
@@ -556,7 +599,13 @@ test_that("LINPRED with random effect", {
556599
beta_x[1:3] ~ nimbleMacros::FORLOOP(dnorm(0, sd = sd_x))
557600
})
558601
)
559-
602+
expect_equal(
603+
out6$modelInfo$inits,
604+
list(beta_x2 = structure(c(x2d = 0, x2e = 0), dim = 2L, dimnames = list(c("x2d", "x2e"))),
605+
beta_x = structure(c(xa = 0, xb = 0, xc = 0), dim = 3L, dimnames = list(c("xa", "xb", "xc"))),
606+
beta_x2_x = structure(c(0, 0, 0, 0, 0, 0), dim = 2:3, dimnames = list(c("x2d", "x2e"), c("xa", "xb", "xc"))), sd_x = 1)
607+
)
608+
560609
# Generate error when trying to get random slope for factor
561610
code6 <- quote(y[1:n] ~ LINPRED(~ (x2|x), priorSpecs=NULL ))
562611
expect_error(LINPRED$process(code6, modInfo, NULL))
@@ -581,30 +630,77 @@ test_that("LINPRED with 'centered' random effect", {
581630
out$modelInfo$constants,
582631
modInfo$constants
583632
)
633+
expect_equal(
634+
out$modelInfo$inits,
635+
list(beta_x3 = c(x3 = 0),
636+
beta_x = structure(c(xa = 0, xb = 0, xc = 0), dim = 3L, dimnames = list(c("xa", "xb", "xc"))),
637+
beta_Intercept = 0)
638+
)
584639

640+
code2 <- quote(LINPRED_PRIORS(~x3 + (1|x), centerVar=x))
641+
out2 <- LINPRED_PRIORS$process(code2, modInfo, NULL)
642+
expect_equal(out2$modelInfo$inits, c(out$modelInfo$inits, list(sd_x=1)))
643+
585644
code <- quote(y[1:n] ~ LINPRED(~x3 + (x3|x), priorSpecs=NULL, centerVar=x))
586645

587646
out <- LINPRED$process(code, modInfo, NULL)
588647
expect_equal(
589648
out$code,
590649
quote(y[1:n] <- nimbleMacros::FORLOOP(beta_x[x[1:n]] + beta_x_x3[x[1:n]] * x3[1:n]))
591650
)
592-
651+
expect_equal(
652+
out$modelInfo$inits,
653+
list(beta_x = structure(c(xa = 0, xb = 0, xc = 0), dim = 3L, dimnames = list(c("xa", "xb", "xc"))),
654+
beta_x_x3 = c(xa = 0, xb = 0, xc = 0),
655+
beta_Intercept = 0, beta_x3 = c(x3 = 0))
656+
)
657+
593658
code <- quote(y[1:n] ~ LINPRED(~x3 + (x3|x) + (1|x2), priorSpecs=NULL, centerVar=x))
594659

595660
out <- LINPRED$process(code, modInfo, NULL)
596661
expect_equal(
597662
out$code,
598663
quote(y[1:n] <- nimbleMacros::FORLOOP(beta_x[x[1:n]] + beta_x2[x2[1:n]] + beta_x_x3[x[1:n]] * x3[1:n]))
599664
)
665+
expect_equal(
666+
out$modelInfo$inits,
667+
list(beta_x = structure(c(xa = 0, xb = 0, xc = 0), dim = 3L, dimnames = list(c("xa", "xb", "xc"))),
668+
beta_x2 = structure(c(x2d = 0, x2e = 0), dim = 2L, dimnames = list(c("x2d", "x2e"))),
669+
beta_x_x3 = c(xa = 0, xb = 0, xc = 0),
670+
beta_Intercept = 0, beta_x3 = c(x3 = 0))
671+
)
672+
673+
code2 <- quote(LINPRED_PRIORS(~x3+(x3|x)+(1|x2), centerVar=x))
674+
out2 <- LINPRED_PRIORS$process(code2, modInfo,NULL)
675+
expect_equal(out2$modelInfo$inits,
676+
c(out$modelInfo$inits, list(sd_x=1, sd_x_x3=1, sd_x2=1)))
600677

601678
code <- quote(y[1:n] ~ LINPRED(~(x3|x), priorSpecs=NULL, centerVar=x))
602679
out <- LINPRED$process(code, modInfo, NULL)
603680
expect_equal(
604681
out$code,
605682
quote(y[1:n] <- nimbleMacros::FORLOOP(beta_x[x[1:n]] + beta_x_x3[x[1:n]] * x3[1:n]))
606683
)
684+
expect_equal(
685+
out$modelInfo$inits,
686+
list(beta_x = structure(c(xa = 0, xb = 0, xc = 0), dim = 3L, dimnames = list(c("xa", "xb", "xc"))),
687+
beta_x_x3 = c(xa = 0, xb = 0, xc = 0),
688+
beta_Intercept = 0)
689+
)
690+
691+
# NOTE: The below code does not provide a prior for beta_x3, which it probably should
692+
# you can fix this by adding x3 manually to the fixed part of the formula
693+
# not sure if this is a bug exactly? Probably need to check it somehow
694+
code2 <- quote(LINPRED_PRIORS(~(x3|x), centerVar=x))
695+
out2 <- LINPRED_PRIORS$process(code2, modInfo, NULL)
607696

697+
#test_code <- nimbleCode({
698+
# mu[1:n] <- LINPRED(~(x3|x), centerVar=x)
699+
# y[1:n] ~ FORLOOP(dnorm(mu[1:n], sd = sd_res))
700+
# sd_res ~ dunif(0, 100)
701+
#})
702+
#mod <- nimbleModel(test_code, constants=c(modInfo$constants, list(n=10, y=rnorm(10))))
703+
#samples <- nimbleMCMC(mod, niter=10, nburnin=0)
608704
})
609705

610706
test_that("LINPRED with factor array covariate", {
@@ -619,6 +715,11 @@ test_that("LINPRED with factor array covariate", {
619715
quote(y[1:M, 1:J] <- nimbleMacros::FORLOOP(beta_Intercept + beta_x[x[1:M, 1:J]]))
620716
)
621717
expect_equal(dim(out$modelInfo$constants$x), c(3,4))
718+
expect_equal(
719+
out$modelInfo$inits,
720+
list(beta_Intercept = 0,
721+
beta_x = structure(c(xa = 0, xb = 0, xc = 0), dim = 3L, dimnames = list(c("xa", "xb", "xc"))))
722+
)
622723

623724
p <- nimble:::codeProcessModelMacros(code, modInfo, environment())
624725
expect_true(is.numeric(p$modelInfo$constants$x))

0 commit comments

Comments
 (0)