model$calculate("node") unexpectedly produces a positive log prob value
gesta...@gmail.com
I have a survival model where I specify an exponentially decreasing age effect on the log hazard (age_effect = b1 * exp(-b2 * age)). When I build the model and do model$calculate(), I get a sensible, finite, negative value. But when I do model$calculate("b2"), I unexpectedly get a positive value, and I am not sure how that is possible.
If I reparameterize the model, and define b3 = b2 * 100, and age_effect = b1 * exp(-b3/100 * age), and then build the model and do model$calculate("b3"), I now get a negative value, as expected.
In the first case, how is it possible that I am getting a positive value, and is it indicative of a problem in the model? If I go ahead and run the model, it produces estimates that make sense and seem reasonable.
Why does simple rescaling of the exponential rate parameter (and appropriately adjusting the prior) seem to correct the problem (if it is a problem)?
Thanks,
Glenn
Below is a simple working example to illustrate the issue. The use of a custom nimble function in the model is because in application I use a much more complicated likelihood that is not possible to specify without use of custom distributions.
library(nimble)
#-------------------------------------
# generate data
#-------------------------------------
b0 <- -6
b1 <- 3
b2 <- 0.1
age_effect <- b1 * exp(-b2 * 0:1500)
haz <- exp(b0 + age_effect)
p_death <- 1- exp(-haz)
s <- rep(0, 1000)
for(i in 1:1000){
event <- rbinom(1500, 1, p_death)
s[i] <- ifelse(max(event) == 0, 5000, min(which(event == 1)) + 1)
}
r <- ifelse(s > 19, 19, s-1)
s <- ifelse(s > 19, Inf, s)
z <- ifelse(is.infinite(s), r, s)
#-------------------------------------
# nimble function for data likelihood
#-------------------------------------
dLik <- nimble::nimbleFunction(
run = function( ### argument type declarations
x = integer(0),
n_samples = integer(0), # number of samples in dataset
r = double(1), # r, last known alive
s = double(1), # s, first known dead
z = double(1), # z, terminal interval (r or s)
beta0s = double(0),
age_effect_surv = double(1),
log = double()) {
# start the loop through individuals
sumllik <- 0
lam <- exp(beta0s + age_effect_surv)
for (i in 1:n_samples) {
llive_r <- 0
ldead <- 1
if(1 < r[i]){
llive_r <- -sum(lam[1:(r[i] - 1)])
}
if(s[i] != Inf){
ldead <- 1-exp(-sum(lam[r[i]:(s[i] - 1)]))
}
llik <- llive_r + log(ldead)
sumllik <- sumllik + llik
}#end n_samples
returnType(double(0))
if (log) {
return(sumllik)
} else {
return(exp(sumllik))
} ## return log-likelihood
})
nimble::registerDistributions(list(
dLik = list(
BUGSdist = "dLik(n_samples, r, s, z, beta0s, age_effect_surv)",
types = c(
"value = integer(0)",
"n_samples = integer(0)",
"r = double(1)",
"s = double(1)",
"z = double(1)",
"beta0s = double(0)",
"age_effect_surv = double(1)",
"log = double(0)"
),
discrete = TRUE
)
))
assign("dLik", dLik, envir = .GlobalEnv)
#-------------------------------------
# nimble model
#-------------------------------------
modelcode <- nimbleCode({
# priors
beta0s ~ dnorm(-6, sd = 2)
b_age1 ~ dgamma(30, 10)
b_age2 ~ dgamma(5, 50)
for (i in 1:19) {
age_effect_surv[i] <- b_age1 * exp(-b_age2 *(i-1))
}
# the data likelihood function
y_surv ~ dLik(
n_samples = n_samples,
r = r[1:n_samples],
s = s[1:n_samples],
z = z[1:n_samples],
beta0s = beta0s,
age_effect_surv = age_effect_surv[1:19]
)
})#end model statement
#-------------------------------------
# set up and run the model
#-------------------------------------
nimData <- list(
y_surv = 1,
r = r,
s = s,
z = z
)
nimConsts <- list(
n_samples = length(s)
)
initsFun <- function()list(
beta0s = rnorm(1, -5.5, 0.1),
b_age1 = rgamma(1, 30,10),
b_age2 = rgamma(1, 5, 50)
)
init <- initsFun()
parameters <- c("beta0s", "b_age1", "b_age2")
ni <- 120000
nb <- .2
bin <- ni * nb
nt <- 100
nc <- 3
Rmodel <- nimbleModel(code = modelcode,
name = "modelcode",
data = nimData,
constants = nimConsts,
inits = init)
Rmodel$calculate()
Rmodel$calculate("b_age2")
# Cnim <- compileNimble(Rmodel)
#
# confMCMC <- configureMCMC(Rmodel,
# monitors = parameters,
# thin = nt,
# useConjugacy = FALSE)
# nimMCMC <- buildMCMC(confMCMC)
# CnimMCMC <- compileNimble(nimMCMC,
# project = Rmodel)
#
# out <- runMCMC(CnimMCMC,
# niter = ni,
# nburnin = bin,
# nchains = nc,
# samplesAsCodaMCMC = TRUE
# )
#-------------------------------------
# reparameterize b_age2
#-------------------------------------
modelcode2 <- nimbleCode({
# priors
beta0s ~ dnorm(-6, sd = 2)
b_age1 ~ dgamma(30, 10)
b_age2 ~ dgamma(10, 1)
for (i in 1:19) {
age_effect_surv[i] <- b_age1 * exp(-b_age2/100 *(i-1))
}
# the data likelihood function
y_surv ~ dLik(
n_samples = n_samples,
r = r[1:n_samples],
s = s[1:n_samples],
z = z[1:n_samples],
beta0s = beta0s,
age_effect_surv = age_effect_surv[1:19]
)
})#end model statement
initsFun2 <- function()list(
beta0s = rnorm(1, -5.5, 0.1),
b_age1 = rgamma(1, 30,10),
b_age2 = rgamma(1, 10, 1)
)
init2 <- initsFun2()
Rmodel2 <- nimbleModel(code = modelcode2,
name = "modelcode2",
data = nimData,
constants = nimConsts,
inits = init2)
Rmodel2$calculate()
Rmodel2$calculate("b_age2")
Perry de Valpine
--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/nimble-users/b06c4999-337d-417d-ae8c-488b965ab53dn%40googlegroups.com.