Prior for covariance (or precision) matrix

20 views
Skip to first unread message

Miranda Fix

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

Thanks in advance for your help!!
-Miranda

Daniel Gibson

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

David Pleydell

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


Miranda Fix

unread,
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