# Prior for covariance (or precision) matrix

20 views

### Miranda Fix

Apr 30, 2021, 2:27:23 PMApr 30
to nimble-users
Hi there,

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)

-Miranda

### Daniel Gibson

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]
}
}
for(k in 1:n.ind){
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.

### David Pleydell

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