Hi all,
After running the model with the code specified below, the code takes too long, showing the following message:
I am running the code in R, and I specified the data in the following way:
*herpesdata <- data14NOMISSING[order(Time), ]
nobs = nrow(herpesdata) - 60 # number of observations in times = 0
timescero <- filter(herpesdata,Time==0)
nrow(timescero)
t0 = 0.0
y0 = as.array(c(1756, 654), dim = 2)
ts = (herpesdata$Time/218) + 1
tsg0 = filter(herpesdata,Time>0) %>% select(Time)
tsg0 <- as_vector(tsg0)
tsg0 <- as.vector(tsg0)
indivs = length(unique(herpesdata$ID))
antibodies = filter(herpesdata,Time>0) %>% select(AB)
antibodies <- as_vector(antibodies)
antibodies <- as.vector(antibodies)
subj = herpesdata$ID
subjlength = length(unique(subj))
datalistg0 = list(nobs = nobs, t0 = t0, y0 = y0,
ts = tsg0, indivs = indivs, antib = antibodies)
z1 <- filter(herpesdata,Time==0) %>% select(AB)
z1 <- as_vector(z1)
z1 <- as.vector(z1)
z2 <- rnorm(indivs, 654, 0.01)
zinitials <- cbind(z1, z2)
zinitials2 <- rbind(z1, z2)
init = function(){
list(pAB0= -2, uAB = 0.053, uASC0 = -2, sigma = 0.01,
sigmapAB = 0.01, sigmauASC = 0.01, rpAB = rnorm(indivs, 0.001, 0.01),
ruASC = rnorm(indivs, 0.001, 0.01), z_init = zinitials)
}*
Any ideas about what might be going wrong (if there is something), that can be slowing down when I run the model?
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;
real ts[nobs];
int <lower=1> indivs;
real <lower=0> antib[nobs];
//real <lower=1, upper=indivs> subj[nobs];
}
parameters {
real pAB0;
real <lower=0> uAB;
real uASC0;
real <lower=0> sigma;
real <lower=0> sigmapAB;
real <lower=0> sigmauASC;
real rpAB[indivs];
real ruASC[indivs];
vector<lower=0>[2] z_init[indivs];
}
model {
vector[2] yhat[nobs];
//prior distributions
pAB0 ~ normal(-2, 0.5);
uASC0 ~ normal(-2, 0.5);
uAB ~ normal(0, 0.001);
sigmapAB ~ inv_gamma(0.01, 0.01);
sigmauASC ~ inv_gamma(0.01, 0.01);
sigma ~ inv_gamma(0.01, 0.01);
for (j in 1:indivs) {
real pAB = exp(pAB0)*exp(rpAB[j]);
real uABt = uAB;
real uASC = exp(uASC0)*exp(ruASC[j]);
vector[2] zinitials;
zinitials[1] = z_init[j,1];
zinitials[2] = z_init[j,2];
yhat = ode_rk45(model3, zinitials, t0, ts, pAB, uABt, uASC);
//likelihood
for (i in 1:nobs) {
antib[i] ~ normal(yhat[i,2], sigma);
// where mu of the r.e. is 0 or -sigmapAB/2
rpAB[j] ~ normal(0.01, sigmapAB);
ruASC[j] ~ normal(0.01, sigmauASC);
z_init[j,1] ~ normal(1764, 0.001);
z_init[j,2] ~ normal(654, 0.001);
}
}
}
generated quantities {
real z_pred[nobs];
for (t in 1:nobs){
z_pred[t] = normal_rng(log( antib[t] ), sigma);
}
}
Edited by @jsocolar for syntax highlighting.