Extract summary fit stan model

Modeling
1

Hi all, I’m trying extract mean of Y_pred Posterior predictive distribution

cat(file = “VB.stan”, "
data {
int N;
vector[N] dt;
vector[N] L;
vector[N] A;
}

parameters {
real<lower = 0> Linf;
real<lower = 0> k;
real<lower = 0> tau;
}

transformed parameters {
vector[N] mu=Linf * (1 - exp(-k * (A + dt)));
real sigma;
sigma = 1 / sqrt(tau); 
}

model{
Linf ~ normal(50,10)T[0,];
k ~ normal(0, 0.1)T[0,];
tau ~ gamma(.0001, .0001);
L ~ normal(mu, sigma);
}
generated quantities {
real Y_mean[N]; 
real Y_pred[N];

for (i in 1:N) {
Y_mean[i] = Linf * (1 - exp(-k * (A[i] + dt[i])));   
Y_pred[i] = normal_rng(Y_mean[i], sigma);
}
}
")

Y_pred_mean <- apply(Y_pred$Y_pred, 2, mean)
plot(L ~ A, xlab=“Longitud total (mm)”, ylab=“Edad estimada (años)”)
lines(A, Y_pred_mean)

Gave me multiple curves of Y_pred_mean,

Y_pred_mean

I want to know how to extrac only the mean of all iterations for 3 chains, thanks!

2

Apply the get_posterior_mean function to the stanfit object and filter down to the rows that pertain to Y_pred or apply the extract function to the stanfit object with the argument pars = Y_pred and take the means over the iterations.

3

I can’t, this is Y_pred Posterior predictive distribution extract

Y_pred ← extract(fitstan, “Y_pred”)

str(Y_pred)
List of 1
Y_pred: num [1:1500, 1:2163] 28.2 21.7 33.4 37.4 25.6 ... ..- attr(*, "dimnames")=List of 2 .. .. iterations: NULL
… …$ : NULL
When I extract the mean, gave me:
Y_pred<-get_posterior_mean(fitstan,pars=“Y_pred”)
str(Y_pred)
num [1:2163, 1:4] 32.9 34.8 29.9 33.8 35.6 …

  • attr(*, “dimnames”)=List of 2
    … : chr [1:2163] "Y_pred[1]" "Y_pred[2]" "Y_pred[3]" "Y_pred[4]" ... .. : chr [1:4] “mean-chain:1” “mean-chain:2” “mean-chain:3” “mean-all chains”

I don’t know how apply mean

4

Those are the means, one column for each of the three chains and then the last column is the mean of all chains. The rows are named by parameter.

5

Exactly when I apply this, gave me multiple line in my plot and not only one line predict
lines(A, Y_pred[,4])

maybe my model is wrong

6

Well, lines is going to be weird unless A is in ascending order because it tries to connect consecutive dots. Even then the line segments are going to be noisy. I would just plot the points.