Extract summary fit stan model

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!

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.

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

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.

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

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.