Prior for covariance (or precision) matrix

[Miranda Fix]

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]

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]

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]

Thanks so much Dan & David!

I will try these out.

Best,

Miranda