I'm running the below Bayesian mixing model in R using the rjags package, but I am having difficultly in specifying the scale matrix for the Wishart distribution. Essentially, I want Sigma.inv to be a covariance matrix that accounts for instrumental measurement error as determined from repeat analysis of a control sample. I have chosen the Wishart distribution as the Sigma.inv prior, with the scale matrix taking the form of the covariance matrix of the control samples.However, this results in the model failing to converge, even after 1 million runs with 500,000 burn-in and 125 thinning length, and tiny errors are generated for the estimates of p.
So my question is should the scale matrix be in a different form? Or is there perhaps an alternative to the Wishart distribution that would be more appropriate in this case? I should add that the model works quite well if a use the identity matrix to scale the Wishart distribution, but of course this doesn't represent the measurement error that I want to incorporate.
Any help will be most appreciated.
model
{
# Likelihood
for(i in 1:N) {
Y[i,1:J] ~ dmnorm(mu[i,],Sigma.inv)
mu[i,1:J] <- p[i,1:K]%*%s[1:K,1:J,i]
}
# Priors on s
for(i in 1:N) {
for(k in 1:K) {
s[k,1:J,i] ~ dmnorm(mu.s[,k],Sigma.s.inv[,,k])
}
}
# Prior on p
for(i in 1:N) {
p[i,1:K] <- expphiV[i,1:K]/sum(expphiV[i,1:K])
phiV[i,1:K] <- phi[i,1:(K-1)]%*%tV
for(k in 1:K) {
expphiV[i,k] <- exp(phiV[i,k])
}
}
# Prior on phi
for(i in 1:N) {
for(k in 1:(K-1)) {
phi[i,k] ~ dnorm(0,kappa.inv[k])
}
philast[i] <- -sum(phi[i,1:(K-1)])
}
for(k in 1:(K-1)) {
kappa.inv[k] ~ dgamma(2,1)
kappa[k] <- 1/kappa.inv[k]
}
# Prior on Sigma
Sigma.inv ~ dwish(R.Sigma,k.Sigma)
}
#
R Code to run model
#
# Y values
Y<-matrix(c(5.8565, 18.5979, 0.0045, 5.7579, 1.1401, 0.4945, 0.2122, 0.2613, 0.3984,
14.0696, 10.6925, 0.0084, 8.5126, 2.1779, 1.1406, 0.4051, 0.3472, 0.5971,
13.7417, 10.0652, 0.0084, 8.2977, 2.2172, 1.1204, 0.3929, 0.3583, 0.6143,
11.7867, 15.1036, 0.0072, 8.3222, 1.9592, 0.9697, 0.3454, 0.3067, 0.5695),
byrow=TRUE, ncol=9)
colnames(Y)<-c("Al","Ca","Ce","Fe","K","Mg","Na","P","Ti")
# mu.s values
mu.s<-matrix(c(6.97463333, 6.657444444, 10.560500000, 13.944108333,
35.47273333, 13.221666667, 6.794500000, 4.149395833,
0.00358000, 0.004533333, 0.008404167, 0.008891667,
5.04523333, 7.128111111, 6.149708333, 6.932004167,
1.19620000, 1.121222222, 2.086083333, 2.446995833,
0.61910000, 0.485333333, 1.010833333, 0.880383333,
0.19653333, 0.199222222, 0.491125000, 0.406904167,
0.07113333, 0.312111111, 0.339916667, 0.298254167,
0.45840000, 0.307222222, 0.616708333, 0.658416667),
byrow=TRUE, ncol=4)
rownames(mu.s)<-c("Al","Ca","Ce","Fe","K","Mg","Na","P","Ti")
colnames(mu.s)<-c("CB","FD","RV","TS")
# Covariance matrix accounting for measurement error of control samples
Sigma.inv<-structure(c(7.516889e-03, 4.719111e-04, 5.117778e-06, 4.896822e-03, 2.683444e-03, 5.368000e-04, 4.266667e-05, -8.364444e-05, 4.630667e-04,
4.719111e-04, 8.097005e-05, 1.132174e-06, 8.006300e-04, 3.065903e-04, 4.421449e-05, -4.869565e-06, -5.677295e-06, 7.178551e-05,
5.117778e-06, 1.132174e-06, 9.927053e-08, 1.090671e-05, 4.145217e-06, 2.788406e-07, -4.459903e-07, -2.015942e-07, 9.478261e-07,
4.896822e-03, 8.006300e-04, 1.090671e-05, 9.350603e-03, 3.260499e-03, 4.947188e-04, -3.441932e-05, -4.135604e-05, 7.779478e-04,
2.683444e-03, 3.065903e-04, 4.145217e-06, 3.260499e-03, 1.390332e-03, 2.255014e-04, -9.909179e-06, -3.122995e-05, 2.906319e-04,
5.368000e-04, 4.421449e-05, 2.788406e-07, 4.947188e-04, 2.255014e-04, 7.714783e-05, 1.326377e-05, -6.255072e-06, 4.525217e-05,
4.266667e-05, -4.869565e-06, -4.459903e-07, -3.441932e-05, -9.909179e-06, 1.326377e-05, 1.484058e-05, 2.859903e-06, -3.663768e-06,
-8.364444e-05, -5.677295e-06, -2.015942e-07, -4.135604e-05, -3.122995e-05, -6.255072e-06, 2.859903e-06, 4.060386e-06, -5.611594e-06,
4.630667e-04, 7.178551e-05, 9.478261e-07, 7.779478e-04, 2.906319e-04, 4.525217e-05, -3.663768e-06, -5.611594e-06, 7.074783e-05), dim=c(9,9))
Sigma.inv<-as.positive.definite(Sigma.inv)
# Sigma.s.inv values - a precision matrix
Sigma.s.inv<-structure(c(53.2149530065154, 7.78555807206446, -23455.5492262134,
3.98226035437711, 181.147593576025, -554.144707627592, 810.959913329385,
-222.049892052889, -676.160229178196, 7.78555807206446, 1.66555095380872,
-2875.24335449319, 0.59057208979672, 35.2056764704436, -90.5513686454713,
122.419754836375, -24.1464188840069, -110.985869289511, -23455.5492262134,
-2875.24335449319, 26623099.2220065, -551.965187580738, -109728.586986558,
276663.105767497, -386431.378216831, -82170.2904382779, 259038.414073531,
3.98226035437711, 0.59057208979672, -551.965187580738, 0.918694682125983,
9.31656679803296, -36.0723479510435, 38.501457515572, -39.2951967395694,
-43.3516097843844, 181.147593576025, 35.2056764704436, -109728.586986558,
9.31656679803296, 1046.08776695798, -2268.37491730957, 2936.52705816089,
60.3905421599093, -2791.52606704426, -554.144707627592, -90.5513686454713,
276663.105767497, -36.0723479510435, -2268.37491730957, 6336.11679874626,
-8941.85134505843, 1263.96712724819, 7151.07272932156, 810.959913329385,
122.419754836375, -386431.378216831, 38.501457515572, 2936.52705816089,
-8941.85134505843, 14092.5591797243, -2587.63879760624, -10498.261488095,
-222.049892052889, -24.1464188840069, -82170.2904382779, -39.2951967395694,
60.3905421599093, 1263.96712724819, -2587.63879760624, 5072.16198995712,
2677.23974494659, -676.160229178196, -110.985869289511, 259038.414073531,
-43.3516097843844, -2791.52606704426, 7151.07272932156, -10498.261488095,
2677.23974494659, 10377.5013426694, 39071331365.1695, -857206689.09299,
13250725351142.8, 4293968882.52526, -138526471018.157, -346806050176.474,
-185731096148.676, -40154997287.1756, -13832086879.5978, -857206689.09299,
45035996.2737049, -291049796239.663, -94316253.2767046, 3042708991.03323,
7617532441.71197, 4079550081.74548, 881997284.581011, 303819297.524648,
13250725351142.8, -291049796239.663, 4499056063320832, 1457943322445.14,
-47034282019304.7, -117752032879430, -63061801053651.4, -13633939080436.7,
-4696447331898.41, 4293968882.52526, -94316253.2767046, 1457943322445.14,
513861374.526114, -15241712287.5009, -38158180147.8331, -20435516111.352,
-4418151354.45332, -1521909040.23616, -138526471018.157, 3042708991.03323,
-47034282019304.7, -15241712287.5009, 491753442721.598, 1231009860799.41,
659264193054.762, 142532685331.566, 49097861284.5838, -346806050176.474,
7617532441.71197, -117752032879430, -38158180147.8331, 1231009860799.41,
3081922920585.39, 1650491845606.43, 356835753178.233, 122918278499.906,
-185731096148.676, 4079550081.74548, -63061801053651.4, -20435516111.352,
659264193054.762, 1650491845606.43, 883961734443.138, 191102478025.129,
65828570726.718, -40154997287.1756, 881997284.581011, -13633939080436.7,
-4418151354.45332, 142532685331.566, 356835753178.233, 191102478025.129,
41361275161.6651, 14232113705.0418, -13832086879.5978, 303819297.524648,
-4696447331898.41, -1521909040.23616, 49097861284.5838, 122918278499.906,
65828570726.718, 14232113705.0418, 4947128692.00388, 123.244158717727,
15.1472733328015, -29469.8213424378, -162.933686694468, -830.16396214664,
94.7853096969282, -68.8560246643179, 741.612274164014, 2264.45755928194,
15.1472733328015, 5.35479274130578, -399.822164505876, -26.2908447667995,
-48.3456870303433, -20.9894871349876, -7.0699744151633, 67.0416328725067,
116.600590769204, -29469.8213424378, -399.822164505876, 25454078.4009846,
34253.206150171, 238183.903396863, -33771.2583884048, 38523.5333162132,
-288045.203177399, -1016390.47602203, -162.933686694468, -26.2908447667995,
34253.206150171, 237.878595765458, 986.561801969833, -68.5166440497585,
209.1509212556, -915.217874817308, -2932.38505457965, -830.16396214664,
-48.3456870303433, 238183.903396863, 986.561801969833, 6669.29110276054,
-1078.91990285319, 352.917095152379, -5609.73421176863, -17946.4520240348,
94.7853096969282, -20.9894871349876, -33771.2583884048, -68.5166440497585,
-1078.91990285319, 1052.57768971324, -913.546632402807, 484.410845136349,
2431.79560931122, -68.8560246643179, -7.0699744151633, 38523.5333162132,
209.1509212556, 352.917095152379, -913.546632402807, 2990.61062774375,
525.038417043516, -3685.88997836577, 741.612274164014, 67.0416328725067,
-288045.203177399, -915.217874817308, -5609.73421176863, 484.410845136349,
525.038417043516, 7256.45767695138, 17207.7716622753, 2264.45755928194,
116.600590769204, -1016390.47602203, -2932.38505457965, -17946.4520240348,
2431.79560931122, -3685.88997836577, 17207.7716622753, 62661.2158837379,
7.0546793436162, -0.0754213178285816, -2357.02875362113, 2.50749046482452,
12.5166502527345, -86.6941792229212, 1.66923366612312, 7.93904545083449,
-126.082128066627, -0.0754213178285816, 0.618360785191319, 1538.6589391227,
-0.31799222497181, -1.26341570245756, 2.43556703274056, 18.6459686649981,
-1.52020047556644, -1.33047266147014, -2357.02875362113, 1538.6589391227,
7450231.40845542, -1467.47518565875, -8484.65408805661, 38482.6779430954,
36677.1639056923, 5037.57654299884, -16285.3259202736, 2.50749046482452,
-0.31799222497181, -1467.47518565875, 1.63035346006115, 1.13929112479727,
-17.7642468677628, 2.98912100156106, 7.74881374361886, -54.902708360705,
12.5166502527345, -1.26341570245756, -8484.65408805661, 1.13929112479727,
106.706180072993, -244.571429668905, 24.7937851688084, -114.000942411739,
-690.553652091766, -86.6941792229212, 2.43556703274056, 38482.6779430954,
-17.7642468677628, -244.571429668905, 1534.27886445762, 549.368288309068,
64.0901199201054, 980.291164452268, 1.66923366612312, 18.6459686649981,
36677.1639056923, 2.98912100156106, 24.7937851688084, 549.368288309068,
1781.88219182038, -107.322815612669, -1733.09607995999, 7.93904545083449,
-1.52020047556644, 5037.57654299884, 7.74881374361886, -114.000942411739,
64.0901199201054, -107.322815612669, 507.111649502059, 661.404359682356,
-126.082128066627, -1.33047266147014, -16285.3259202736, -54.902708360705,
-690.553652091766, 980.291164452268, -1733.09607995999, 661.404359682356,
9442.32984843676), .Dim = c(9L, 9L, 4L))
#N
N<-nrow(Y)
#K
K<-4
#J
J<-9
myjagsfile<-"JAGSmodel.txt" # the above jags model
tV = t(ilrBase(D=K))
mydata <-list(N=N, J=J, K=K, Y=Y, mu.s=mu.s, Sigma.s.inv=Sigma.s.inv, R.Sigma=Sigma.inv, k.Sigma=J, tV=tV)
mypars <-c("Sigma.inv","p","kappa","phi", "s")
myiterations<-25000
myburnin<-12500
mythin<-175
set.seed(5000)
run.mcmc <- jags(data=mydata, inits=NULL, parameters.to.save=mypars, model.file=myjagsfile, n.chains=2,
n.iter=myiterations, n.burnin=myburnin, n.thin=mythin, DIC=TRUE)