Hi all,
I have no idea why the following code have the issues, if anyone can help with this, i really appreciate it.
dat<-read.table(header = TRUE, text = "
id[] dose[] n[] y[,1] y[,2] y[,3] y[,4] y[,5] y[,6]
1 10 3 -0.43 -0.64 0 NA NA NA
2 10 3 0.35 -2 -0.92 NA NA NA
3 10 3 0.07 -0.07 0.15 NA NA NA
4 10 3 0.29 -1.42 0.14 NA NA NA
5 50 3 0.57 -2.57 -1.07 NA NA NA
6 10 3 0.08 0.43 -1.07 NA NA NA
7 50 3 -0.36 -0.14 0.64 NA NA NA
8 10 3 -0.21 1.28 -1.5 NA NA NA
9 10 3 0.68 -0.25 -0.09 NA NA NA
10 10 3 -0.64 1 -0.29 NA NA NA
11 10 4 0.05 -0.22 0.57 0.36 NA NA
12 25 6 0.64 1.08 -0.36 0.79 -0.64 1.5
13 5 3 -0.29 1.08 1.5 NA NA NA
14 20 4 4.29 3.15 0.78 4.49 NA NA
15 20 3 0.64 0 0.5 NA NA NA
16 10 4 1.22 1.07 -0.08 0.5 NA NA
17 25 4 -0.08 0.86 1.07 1.15 NA NA
18 10 3 1.85 0.78 0.54 NA NA NA
19 50 4 0.86 1.43 0.65 1.86 NA NA
20 20 3 0.86 1.39 1.85 NA NA NA
21 10 3 1.35 1.36 0.79 NA NA NA
22 10 3 0.47 0.93 0.86 NA NA NA
23 25 3 0.42 0.71 0.66 NA NA NA
")
names(dat) <- gsub(".", "", names(dat), fixed = TRUE)
y1 = as.array(as.numeric(na.exclude (dat$y1)))
y2 = as.array( as.numeric(na.exclude (dat$y2)))
y3 = as.array(as.numeric(na.exclude (dat$y3)))
y4 = as.array(as.numeric(na.exclude (dat$y4)))
y5 = as.array(as.numeric(na.exclude (dat$y5)))
y6 = as.array(as.numeric(na.omit (dat$y6)))
y=as.array(c(y1,y2,y3,y4,y5,y6))
s=dat$n
options(mc.cores = parallel::detectCores())
dataList = list(
Ntotal = length(y),
s = s,
K= length(s),
y= y
)
modelString1 = "
data{
int<lower=0> Ntotal;
int<lower=0> K;
vector[Ntotal] y;
int s[K];
}
parameters {
real theta[K];
real<lower=-10,upper=10> mutheta;
real<lower=0> sigma2[K];
real<lower=-10,upper=10> musigma;
real<lower=0,upper=10> tausigma;
real<lower=0,upper=50> tautheta;
}
model {
int pos;
pos = 1;
for (k in 1:K) {
segment(y, pos,s[k]) ~ normal(theta[k],sqrt(sigma2[k]));
pos = pos + s[k];
theta[k] ~ normal(mutheta, tautheta) ;
sigma2[k] ~ lognormal(musigma,tausigma);
}
mutheta ~ uniform(-10,10);
musigma ~ uniform(-10,10);
tausigma ~ uniform(0,10);
tautheta ~ uniform (0,50);
}
generated quantities {
real<lower=0,upper=1> Ppop ;
real<lower=0,upper=1> Ppopdiff;
real<lower=0> meanvar;
Ppop = mutheta>0;
Ppopdiff = mutheta>0.5;
meanvar = exp(-musigma+tausigma/2);
}
"
parameters = c( "mutheta", "Ppop", "Ppopdiff", "tautheta","tausigma","meanvar")
writeLines( modelString1 , con="TEMPmodel.stan" )
stanModel <- stan(
data = dataList ,
file = "TEMPmodel.stan",
chains = 4,
warmup = 5000,
pars = parameters ,
iter = 25000,
control = list(adapt_delta = 0.99))
Warnings are as follows:
Warning messages:
2: There were 99 divergent transitions after warmup. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
3: Examine the pairs() plot to diagnose sampling problems