Several parameters of our model fail to converge in WinBUGS

Question

I established a Beyesian model, which including some equations from the literature. I want to simultaneously simulate ten parameters (P, R, Rear, K, R20, mu, a, b, tau, and sigma) according to the monitoring data. The model is syntactically correct, loaded data, and compiled well. I updated 1000000. It can give the results of each parameter.

But three parameters (mu, a and b) of our model seem always to fail to converge in WinBUGS. When I only simulated these three parameters, these also failed to converge.

The problem is that I don't know how to adjust that. I tried to change the pior distribution of some paramters many times. But no use. I also don't know how to estimate the model, because the DIC results shows "undefined real results".

I searched our website and find that there was also a simliar question.

I got the following answer: "WinBUGS is several years old, and it is 8 years since its last update! I think you should forget it, because there are several better alternatives. You may try virtually the same code in Jags or Stan, in which both are usable in R via rJags and RStan. Stan is specially important, because it uses MCMC which converges in many situations that WinBUGS does not. " WinBUGS Weibull Network Meta-Analysis

However, I am also new to Jags or Stan. I don't know how to change it to the STAN code.

I would appreciate anyone suggestions about how to adjust this (the model structure, the equation, pior distribution, initial values, etc...) to get it to run properly. In addition, I am tring to do the analyses using R.

Thanks in advance.

Here is my model

    model; 
{ 
for( i in 1 : N ) {
lambda[i] ~ dnorm(u[i],tau)  }

for( i in 1 : N ) {
u[i] <- P[i] - R[i] + Rear[i]
P[i] <- mu * IZ[i]  
IZ[i] <- I[i] * exp(e[i])  
e[i] <- (-1) * Kd[i]
Kd[i] <- a + b * log(Tur[i])

R[i] <- R20 * pow(1.047,T[i] - 20)
Rear[i] <- K * (Os[i] - O[i])
Os[i] <- 14.62 - 0.3671 * T[i] + (0.004497 * T[i]) * T[i] - 0.966 * Sa[i] + (0.00205 * Sa[i]) * T[i] + (2.739E-4 * Sa[i]) * Sa[i]

T[i] ~ dnorm( 0.0,1.0E-6)
I[i] ~ dnorm( 0.0,1.0E-6) 
Tur[i] ~ dnorm(10, 0.1)I(0,)  
Sa[i] ~ dnorm( 0.0,1.0E-6)
O[i] ~ dnorm( 0.0,1.0E-6)  
  }


tau ~ dgamma(0.001,0.001)
R20 ~ dnorm( 0.0,1.0E-6)
K ~ dnorm( 0.0,1.0E-6)
mu ~ dnorm( 0.0,1.0E-6)
b ~ dnorm( 0.0,1.0E-6)
a ~ dnorm( 0.0,1.0E-6)
sigma <- 1 / sqrt(tau) 
} 

Here is my data

list(lambda=c(0.3, 0, 0.03, 0.12,   0.13,   0.12,   0.03,   0.27,   0.29,   0.02,   0.2,    0.25,   0.26,   0.15,   0.16,   0.74,   0.3,    0.4,    0.4,    0.28,   0.15,   -0.15,  -0.07,  0.02,   -0.13,  -0.3,   -0.26,  -0.36,  -0.26,  -0.28,  -0.26,  -0.32,  -0.18,  -0.29,  -0.27,  -0.09,  -0.32,  -0.21,  -0.18,  -0.16,  -0.23,  -0.18,  -0.16,  -0.13,  -0.18,  -0.07,  -0.15,  -0.11,  -0.03,  0.01,   0.03,   0.12,   -0.07,  0.12,   0.08,   0.18,   0.24,   0.3,    0.08,   0.35,   0.27,   0.29,   0.02,   0.2,    0.25,   0.26,   0.24,   0.2,    0.2,    -0.09,  -0.2,   -0.16,  -0.16,  -0.08,  -0.06,  -0.17,  -0.32,  -0.14,  -0.18,  -0.32,  -0.16,  -0.28,  -0.04,  -0.27,  -0.15,  -0.06,  -0.15,  -0.11,  -0.15,  -0.11,  -0.03,  -0.25,  -0.02,  -0.06,  -0.16,  -0.08,  0.08),
  T=c(27.45,27.56,27.33,27.39,27.61,27.61,28,28.4,28.78,29,29.5,29.94,30,30.22,30.56,31,31.4,31.8,31.9,30.82,30.15,30.25,30.25,29.86,29.94,29.76,29.52,29.26,29.02,28.8,28.59,28.33,28.09,27.85,27.63,27.46,27.31,27.17,27.04,26.91,26.77,26.63,26.52,26.41,26.29,26.13,26.01,25.88,25.79,25.79,25.77,25.75,25.76,25.86,25.8,25.85,26.01,26.15,26.27,26.49,26.74,27,27.56,27.87,27.83,28.02,28.36,28.67,29.11,29.45,29.55,29.35,28.67,28.52,28.05,27.98,27.91,27.55,27.26,26.99,26.72,26.46,26.21,25.98,25.71,25.44,25.23,25.01,24.79,24.6,24.41,24.23,24.05,23.88,23.73,23.59,23.48),
I=c(24,52,141,138,179,193,148,306,210,423,306,258,235,385,157,382,598,517,427,603,279,196,235,124,132,87,42,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,7,20,24,52,141,138,179,193,148,306,210,423,306,258,235,385,157,382,598,517,427,603,279,196,235,124,132,87,42,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,7,20,59,90,109,112,86,143,206,248,244,260,384,361,460,563,398,204,510,614,608,588,523,430,332,246,161,73,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,22,90),
 Sa=c(2.67,2.67,2.68,2.68,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.68,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.68,2.68,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.67,2.68,2.67,2.67,2.67,2.68,2.67,2.67,2.67,2.68,2.67,2.67,2.67,2.68,2.67,2.68,2.67,2.68,2.67,2.68,2.71,2.68,2.6,2.68,2.68,2.68,2.68,2.69,2.69,2.69,2.69,2.68,2.68,2.69,2.68,2.7,2.69,2.69,2.695,2.697,2.699,2.701,2.703,2.69,2.7,2.7,2.69,2.7,2.7,2.69,2.7,2.69,2.7,2.7,2.7,2.7,2.7,2.7,2.7,2.7,2.69,2.69,2.7,2.7),
 O=c(6.23,6.53,6.53,6.56,6.68,6.81,6.93,6.96,7.23,7.52,7.54,7.74,7.99,8.25,8.4,8.56,9.3,9.6,10,10.4,10.68,10.83,10.68,10.61,10.63,10.5,10.2,9.94,9.58,9.32,9.04,8.78,8.46,8.28,7.99,7.72,7.63,7.31,7.1,6.92,6.76,6.53,6.35,6.19,6.06,5.88,5.81,5.66,5.55,5.52,5.53,5.56,5.68,5.61,5.73,5.81,5.99,6.23,6.53,6.61,6.96,7.23,7.52,7.54,7.74,7.99,8.25,8.49,8.69,8.89,8.8,8.6,8.44,8.28,8.2,8.14,7.97,7.65,7.51,7.33,7.01,6.85,6.57,6.53,6.26,6.11,6.05,5.9,5.79,5.64,5.53,5.5,5.25,5.23,5.17,5.01,4.93),   

N=97)

list(mu=0.5, a=1, b=2, K=0.5, R20=1, tau=1) 

Part of the Results

ode statistics 
node     mean    sd  MC error   2.5%    median  97.5%   start   sample
mu  713.5   620.3   17.57   0.07538 590.0   2162.0  4001    996000   

node     mean    sd  MC error   2.5%    median  97.5%   start   sample
a   8.924   2.585   0.08136 0.6363  9.707   11.45   4001    996000

the time series of iteration

enter image description here enter image description here

A quick summary of the variables entered into the model:

T : water temperature
I: solar irradiance at water surface
Sa: water salinity
O: disolved oxygen

(P, R, Rear, K, R20, mu, a, b, tau, and sigma)
P: primary prodution 
R: ecosystem respiration
Rear: molecular diffusivity of dissolved oxygen across the air–water interface
K: gas transfer velocity parmeter
R20: respiration at T= 20℃
mu: the slope of the photosynthesis–irradiance relationship at low light conditions
a and b: regression coefficients between turbidity and irradiance attenuation coefficient
tau: the precision?
sigma: the standard error?

Answer 1

It would be written something like this in the Stan language, I think

data {
  int<lower=1> N;
  vector[N] lambda;
  vector[N] T;
  vector[N] I;
  vector[N] Sa;
  vector[N] O;
}
parameters {
  vector<lower=0>[N] Tur;
  real<lower=0> sigma;
  real R20;
  real K;
  real mu;
  real a;
  real b;
}
model {
  vector[N] u;
  for (i in 1:N) {
    real Os = 14.62 - 0.3671 * T[i] + (0.004497 * T[i]) * T[i] - 0.966 * Sa[i] + (0.00205 * Sa[i]) * T[i] + (2.739E-4 * Sa[i]) * Sa[i];
    real Rear = K * (Os - O[i]);
    real R = R20 * 1.047^(T[i] - 20);
    real Kd = a + b * log(Tur[i]);
    real e = -1 * Kd;
    real IZ = I[i] * exp(e);
    real P = mu * IZ;
    u[i] = P - R + Rear;
  }
  target += normal_lpdf(lambda | u, sigma);

  target += normal_lpdf(Tur | 10, 3.162278); // do not need to truncate
  // other priors in the form of target += distribution_lpdf( | , );
}

However, there are several confusing things in your BUGS model. First, you have sampling statements for T, I, Sa, and O even though those are apparently data and the sampling statements only involve constants. Constants are irrelevant to the parameter proposals in Stan or which proposals get accepted, so you can just omit all that. If you need random draws of them in order to predict in some larger population, you can do that in Stan with a generated quantities block at the end of your Stan program or you can do it in R afterward.

Second, these priors you are using --- while common in BUGS models --- are antithetical to inference. For a model like this (and many other models), you are going to have to express what you actually believe about these parameters, rather than saying vacuous things like there is a two-thirds chance that mu is between -1000 and 1000. There are many recommended priors for Stan. Also, keep in mind that the parameterization of the normal distribution in Stan is in terms of an expectation and a standard deviation (rather than a precision).

See also Appendix B of the Stan user manual and the FAQ. If you use much better priors, it may work, but pay attention to the warning messages you will probably get the first time you try it (that indicate the posterior distribution you have defined is not conducive to sampling from).