20 views

Skip to first unread message

Apr 30, 2021, 2:27:23 PMApr 30

to nimble-users

Hi there,

for( j in 1:Nstands ){

a_0[j, 1:2] ~ dmnorm(a_0_mean[1:2], prec = Omega[1:2,1:2])

}

a_0_mean[1:2] ~ dmnorm(mu0[1:2], prec = S2[1:2,1:2])

Omega[1:2,1:2] ~ dwish(R2[1:2,1:2], 2)

Sigma_a_0[1:2,1:2] <- inverse(Omega[1:2,1:2])

sig_pre <- sqrt(Sigma_a_0[1,1])

sig_post <- sqrt(Sigma_a_0[2,2])

cor_pp <- Sigma_a_0[1,2]/(sig_pre*sig_post)

R2 = matrix(c(1,0,0,1),nrow=2)

I am trying to code up a model in NIMBLE with a bivariate random effect using dmnorm. Following some examples in BUGS, I gave the precision matrix a wishart prior using dwish. However, it seems like the results are pretty sensitive to the parameters of this wishart prior. I have also heard that there are issues with using a wishart prior, but I couldn't find any examples of a different way of specifying a prior for a covariance (or precision) matrix in NIMBLE.

Does anyone have any suggestions? Either on how to choose the parameters for the wishart prior, or how to specify a different (better?) prior for a covariance (or precision) matrix?

Below are some relevant lines of code. I do actually care about getting good posterior estimates for sig_pre, sig_post, and cor_pp.

for( j in 1:Nstands ){

a_0[j, 1:2] ~ dmnorm(a_0_mean[1:2], prec = Omega[1:2,1:2])

}

a_0_mean[1:2] ~ dmnorm(mu0[1:2], prec = S2[1:2,1:2])

Omega[1:2,1:2] ~ dwish(R2[1:2,1:2], 2)

Sigma_a_0[1:2,1:2] <- inverse(Omega[1:2,1:2])

sig_pre <- sqrt(Sigma_a_0[1,1])

sig_post <- sqrt(Sigma_a_0[2,2])

cor_pp <- Sigma_a_0[1,2]/(sig_pre*sig_post)

I am currently specifying the following constants (for lack of other ideas), but changing R2 can change the results quite a bit.

S2 = matrix(c(1,0,0,1),nrow=2)/9,R2 = matrix(c(1,0,0,1),nrow=2)

Thanks in advance for your help!!

-Miranda

Apr 30, 2021, 2:48:02 PMApr 30

to Miranda Fix, nimble-users

Hi Miranda,

I may have shared something similar to this on the listserv before, but this is my preferred way at the moment as it can handle any number of correlated features.

Here, n = the number of co-varying features (your case, 2), and n.ind = the number of observations. It is specified to create a variance-covariance matrix as opposed to the precision matrix, so it just needs users to provide a prior for the correlation (rho), which is between -1 and 1, and some relatively diffuse prior for each variance term (sig). I used a half-Cauchy, but you can definitely switch this to a uniform distribution that resembles the potential variance in your system.

for(m in 1:n){

mu[m] <- 0

sig[m] ~ T(dt(0, pow(2.5,-2), 1),0,)

for(i in 1:n){

Sigma[m,i] <- rho[m,i] * sig[m] * sig[i]

}

}

for(mm in 2:n){

for(ii in 1:(mm-1)){

rho[mm,ii] ~ dunif(-1,1)

rho[ii,mm] <- rho[mm,ii]

}

}

mu[m] <- 0

sig[m] ~ T(dt(0, pow(2.5,-2), 1),0,)

for(i in 1:n){

Sigma[m,i] <- rho[m,i] * sig[m] * sig[i]

}

}

for(mm in 2:n){

for(ii in 1:(mm-1)){

rho[mm,ii] ~ dunif(-1,1)

rho[ii,mm] <- rho[mm,ii]

}

}

for(k in 1:n.ind){

eps[1:n,k] ~ dmnorm(mu[1:n], cov = Sigma[1:n,1:n])

}

eps[1:n,k] ~ dmnorm(mu[1:n], cov = Sigma[1:n,1:n])

}

Best

-Dan

--

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 on the web visit https://groups.google.com/d/msgid/nimble-users/98b2bdb5-2e5f-489d-a760-88ba5bca1001n%40googlegroups.com.

Apr 30, 2021, 5:20:21 PMApr 30

to nimble-users

Hi Miranda

An alternative specification for a correlation matrix is the so called LKJ prior

You can hard-wire it directly into your BUGs code using the the code below, which is adapted from here

Cheers

David

p.s. the code uses a transform of the beta distribution on to (-Inf, Inf) - this is purely to avoid Metropolis Hastings wasting time making proposals outside [0,1].

## ############################################

## Model correlation matrix based on cholesky onion method ##

## ############################################

eta <- 1 # Must be one or larger for all Beta distributions to be proper

alpha[1] <- eta + (K - 2)/2

logit(corY[1]) ~ dLogitBeta(alpha[1], alpha[1])

u[1,1] <- 1

u[1,2] <- 2 * corY[1] - 1

u[2,2] <- 2*sqrt(corY[1] - corY[1]^2)

u[2:K,1] <- 0

for (m in 2:(K-1)) {

## Draw beta random variable

alpha[m] <- alpha[m-1] - 0.5

logit(corY[m]) ~ dLogitBeta(m / 2, alpha[m])

## Draw uniformly on a hypersphere

for (jj in 1:m) {

corZ[m, jj] ~ dnorm(mean=0, sd=1)

}

scZ[m, 1:m] <- corZ[m, 1:m] / sqrt(inprod(corZ[m, 1:m], corZ[m, 1:m]))

u[1:m,m+1] <- sqrt(corY[m]) * scZ[m,1:m]

u[m+1,m+1] <- sqrt(1 - corY[m])

u[(m+1):K,m] <- 0

}

## ################################################

## A DISTRIBUTION GIVING THE LOGIT OF A BETA DISTRIBUTION ##

## ################################################

dLogitBeta <- nimbleFunction (

## Returns density of x where

## y ~ Beta(a1,a2)

## x = logit(y)

run = function(x = double(0),

shape1=double(0, default=1.0),

shape2=double(0, default=1.0),

log = integer(0, default=0)) {

returnType(double(0))

y = ilogit(x)

logProbX = log(y) + log(1 - y) + dbeta(y, shape1=shape1, shape2=shape2, log=TRUE) ## Via change of variables

if (log)

return(logProbX)

return(exp(logProbX))

}

)

rLogitBeta <- nimbleFunction (

## Generates y ~ Beta(a1,a2)

## Returns x = logit(y)

run = function(n = integer(0, default=1),

shape1 = double(0, default=1.0),

shape2 = double(0, default=1.0)) {

returnType(double(0))

if(n != 1)

nimPrint("Warning: rLogitBeta only allows n = 1; Using n = 1.\n")

y <- rbeta(1, shape1=shape1, shape2=shape2)

x <- logit(y)

return(x)

}

)

registerDistributions(list(dLogitBeta = list(

BUGSdist = "dLogitBeta(shape1, shape2)",

discrete = FALSE,

types = c("value=double(0)"), ## , "para=double(0)"

pqAvail = FALSE)))

May 3, 2021, 12:02:55 PMMay 3

to nimble-users

Thanks so much Dan & David!

I will try these out.

Best,

Miranda

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu