model$calculate("node") produces positive result

[gesta...@gmail.com]

I have a survival model with an age effect that is exponentially decreasing on the log hazard scale, like age_effect = b1 * exp(-b2 * age), so double-exponential on the hazard scale. When I do model$calculate(), I get a reasonable and expected finite negative number. But when I do model$calculate("b2"), I get a positive number, and I do not understand how this is possible. When I go ahead and run the model, it does seem to produce sensible estimates.

When I reparameterize the model such that age_effect = b1 * exp(-b3/100 * age), where b3 = b2 *100, then model$calulate("b3") produces a negative number, as expected.

In the first case, why do I get a non-negative log prob, and is that indicative of a problem with the model?

Secondly, why does simply rescaling the rate parameter (and appropriately adjusting the prior) change the outcome?

Thanks,

Glenn

Below is a simple working example to illustrate. The nimple function included in the model is because in application I specify much more complicated likelihoods that would not be possible without cunstom nimble functions.

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")

[Chris Paciorek]

This looks to have been a duplicated post that got caught in our spam filter for some time.

--
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/b059292a-ef39-473e-b510-691ea8bdf12fn%40googlegroups.com.