Hi all,
I am running over and over into the same issue:
Even if I add restrictions to my parameters, change some prior distributions and set a list of initial values when I run the model in R as follows, I get the same issue:
nobs = nrow(data14NOMISSING)
t0 = 0.0
zAB0 ← data14NOMISSING$AB0
ts ← data14NOMISSING$Time
indivs = length(unique(data14NOMISSING$ID))
indivs
zASC0 ← rnorm(nobs, 654, 0.5)
antibodies = data14NOMISSING$AB
subj = data14NOMISSING$ID
zinitials = cbind(zAB0, zASC0)
datalist2 ← list(nobs = nobs, t0 = t0, y0 = zinitials,
ts = ts, indivs = indivs, antib = antibodies, subj = subj)
init = function(){
list(pAB0= -2, uAB = 0.053, uASC0 = -2, AB0 = 1766, ASC0 = 654, sigma = 1,
sigmapAB = 0.05, sigmauASC = 0.05, sigmaAB0 = 0.05,
sigmaASC0 = 0.05, rpAB = runif(nobs, 1, 1.5), ruASC = runif(nobs, 1, 1.5),
rAB0 = runif(nobs, 1, 1.5), rASC0 = runif(nobs, 1, 1.5))
}
Modelnosorted ← stan_model(“C:/Users/IGarciaFogeda/Documents/R/Rstan/HIERARCHICAL/Longitudinal/longitudinalmodel.stan”)
fit_modelnosorted<- sampling(Modelnosorted,
data = datalist2, init = init,
iter = 2000,
chains = 1,
seed = 0, control = list (stepsize = 0.1))
Where does the 0 of the log probability come from?
functions {
vector model3(real t, vector y, real pAB, real uAB, real uASC) {
vector[2] dydt;
dydt[1] = pAB*y[2]-uAB*y[1];
dydt[2] = -uASC*y[2];
return dydt;
}
}
data {
int <lower=1> nobs;
real t0;
vector[2] y0[nobs];
real ts[nobs];
int <lower=1> indivs;
real <lower=0> antib[nobs];
int <lower=1> subj[nobs];
}
parameters {
real pAB0;
real <lower=0> uAB;
real uASC0;
real <lower=5> AB0;
real <lower=2> ASC0;
real <lower=0> sigma;
real <lower=0> sigmapAB;
real <lower=0> sigmauASC;
real <lower=0> sigmaAB0;
real <lower=0> sigmaASC0;
real <lower=0> rpAB[nobs];
real <lower=0> ruASC[nobs];
real <lower=0> rAB0[nobs];
real <lower=0> rASC0[nobs];
}
model {
vector[2] yhat[nobs];
vector[2] yhatmu[nobs];
real eval_time[1];
//prior distributions
pAB0 ~ normal(0, 2.5);
uASC0 ~ normal(0, 2.5);
uAB ~ uniform(0, 1);
AB0 ~ uniform(5,10);
ASC0 ~ uniform(2,4);
sigmapAB ~ inv_gamma(0.01, 0.01);
sigmauASC ~ inv_gamma(0.01, 0.01);
sigmaAB0 ~ inv_gamma(0.01, 0.01);
sigmaASC0 ~ inv_gamma(0.01, 0.01);
sigma ~ inv_gamma(0.01, 0.01);
for (j in 1:nobs) {
real pAB = exp(pAB0)*exp(rpAB[subj[j]]);
real uABt = uAB;
real uASC = exp(uASC0)*exp(ruASC[subj[j]]);
vector[2] zinitials;
zinitials[1] = AB0 + rAB0[subj[j]];
//are we sure this is not yielding to negative values? is it hierarchical?
zinitials[2] = ASC0 + rASC0[subj[j]];
// where mu of the r.e. is 0 or -sigmapAB/2
rpAB[subj[j]] ~ normal(-sigmapAB/2, sigmapAB);
ruASC[subj[j]] ~ normal(-sigmauASC/2, sigmauASC);
rAB0[subj[j]] ~ normal(-sigmaAB0/2, sigmaAB0);
rASC0[subj[j]] ~ normal(-sigmaASC0/2, sigmaASC0);
eval_time[1] = ts[j];
if (eval_time[1] == 0){
yhatmu[j,1] = zinitials[1] - pow(sigma, 2.0)/2;
yhatmu[j,2] = zinitials[2] - pow(sigma, 2.0)/2;}
else{
yhat = ode_rk45(model3, zinitials, t0, eval_time, pAB, uABt, uASC);
yhatmu[j,2] = yhat[1,1] - pow(sigma, 2.0)/2;
}
//likelihood
antib[j] ~ lognormal(yhatmu[j,2], sigma);
}
}
generated quantities {
real z_pred[nobs];
for (t in 1:nobs){
z_pred[t] = lognormal_rng(antib[t], sigma);
}
}