Perpexing errors in vectorised GLM

[Bob O'Hara]

I'm trying to implement a weighted Bayesian lasso in Nimble for a binary response (the long version: I want a Bayesian version of MaxEnt so I can do some uncertainty propagation). I'm trying to adapt the final method here: <https://r-nimble.org/examples/linpred.html>, but I'm hitting error messages that I understand but that (to me) don't make sense.

Reproducible code is below. With the weighted function, I'm told "Error in chol.default(model$omega[1:11, 1:11]) : the leading minor of order 1 is not positive". If I set omega to be an identity matrix this doesn't happen, but if I make it a diagonal matrix with dgamma(N,N) (N large) priors I do get this error. This does not seem to be in the part of the model where the weights are, though.

I tried to make a minimal reproducible example by removing the weights, and got a different error: "bad parameters for distribution dbern_vec(p = prob[1:n]). (No available re-parameterization found.)". This doesn't seem to have anything to do with the weights: my guess is that with weights omega fails before we hit this error.

I'm sure this all does make sense, but any help would be gratefully received.

Bob

library(disdat)
library(nimble)

dbern_vec_wt <- nimbleFunction( ## Define the distribution
  run = function(x = double(1), p = double(1), wt = double(1), log = integer(0, default = 0)) { #
    returnType(double(0))
    logProbs <- wt*dbinom(x=x, size=1, prob=p, log = TRUE)
    logProb <- sum(logProbs)
    #    logProb <- sum(dbinom(x=x, size=1, prob=p, log = TRUE))
    if(log) return(logProb)
    else return(exp(logProb))
  })

dbern_vec <- nimbleFunction( ## Define the distribution
  run = function(x = double(1), p = double(1), log = integer(0, default = 0)) {
    returnType(double(0))
    logProbs <- dbinom(x=x, size=1, prob=p, log = TRUE)
    logProb <- sum(logProbs)
    #    logProb <- sum(dbinom(x=x, size=1, prob=p, log = TRUE))
    if(log) return(logProb)
    else return(exp(logProb))
  })

rbern_vec <- nimbleFunction( ## Define a simulation function, optionally.
  run = function(n = integer(0), p = double(1)) {
    returnType(double(1))
    if(n != 1) print("rbern_vec only allows n = 1; using n = 1.")
    smp <- rbinom(n, p, size=1)
    return(smp)
  })

Weightedcode <- nimbleCode({
  for(j in 1:k) {
    ratetau[j] <- MeanLambda
    tau2[j] ~ dexp(rate=ratetau[j])
  }
  omega[1:k, 1:k] <- diag(1/tau2[1:k])
  alpha ~ dnorm(0, sd = 1)
  beta[1:k] ~ dmnorm(zeros[1:k], omega[1:k, 1:k])
  linpred[1:n] <- alpha + X[1:n, 1:k] %*% beta[1:k]
  logit(prob[1:n]) <- linpred[1:n]
  y[1:n] ~ dbern_vec_wt(prob[1:n], W[1:n])
})

Code <- nimbleCode({
  for(j in 1:k) {
    ratetau[j] <- MeanLambda
    tau2[j] ~ dexp(rate=ratetau[j])
  }
  omega[1:k, 1:k] <- diag(1/tau2[1:k])
  alpha ~ dnorm(0, sd = 1)
  beta[1:k] ~ dmnorm(zeros[1:k], omega[1:k, 1:k])
  linpred[1:n] <- alpha + X[1:n, 1:k] %*% beta[1:k]
  logit(prob[1:n]) <- linpred[1:n]
  y[1:n] ~ dbern_vec(prob[1:n])
})


# Use one species from Canadian data
CAN <- disData("CAN")
EnvPO <- c("alt", "asp2", "ontprec", "ontprec4", "ontprecsd", "ontslp",
           "onttemp", "onttempsd", "onttmin4", "ontveg", "watdist")
Species <- "can02"

SpData <- rbind(CAN$po[CAN$po$spid==Species,c("occ", EnvPO)], CAN$bg[,c("occ", EnvPO)])
SpData[,EnvPO] <- apply(SpData[,EnvPO], 2, scale) # scale covariates

mm <- as.matrix(SpData[,EnvPO])
dd <- list(y = SpData$occ)

Params <- c("alpha", "beta")

Const <- list(n = nrow(mm), k = ncol(mm), MeanLambda = 1)
Data <- list(y = SpData$occ, X = mm, W=rep(1,nrow(mm)))
Inits <- list(alpha = 0, beta = rnorm(Const$k, 0, sd=0.01))


Model_wt <- nimbleModel(code = Weightedcode, data = Data, constants = Const, inits = Inits)

# Gives the following error:
# Error in chol.default(model$omega[1:11, 1:11]) :
#  the leading minor of order 1 is not positive

Model <- nimbleModel(code = Code, data = Data, constants = Const, inits = Inits)

# Gives the following error:
# Error: bad parameters for distribution dbern_vec(p = prob[1:n]). (No available re-parameterization found.)

[Daniel Turek]

Bob, good questions.  The error:

# Error in chol.default(model$omega[1:11, 1:11]) :
#  the leading minor of order 1 is not positive

is just a problem with initial values.  Since the dmnorm distribution (under the hood) finds the cholesky factorization of its precision/covariance argument, it needs a valid initial value for this matrix.  Your omega matrix is determined by the values in tau2, so what you need to do is provide valid (positive) initial values for the tau2 vector, and that will provide a valid initial value for the omega matrix, and that error goes away.  You can try, for example:

Inits <- list(alpha = 0, beta = rnorm(Const$k, 0, sd=0.01),
              tau2 = rep(1, Const$k))

For your second example (using nimbleModel(code = Code,...) ), I didn't get the same error you did.  In fact, when using the valid initial values (shown above), I didn't get any error at all from this line:

Model <- nimbleModel(code = Code, data = Data, constants = Const, inits = Inits)

Maybe I made a mistake running your code.  Are you able to try again, and maybe send another reproducible example of this error?

Cheers!

Daniel

--
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/0dc7eeb0-cd39-4dd4-8eda-fa0f0c6c28fbn%40googlegroups.com.

[Bob O'Hara]

Wonderful, thanks. 

Restarting R solved the second problem, so that was at my end. 

Bob