Two different methods to plot residuals with rstan but two different distributions [closed]

The following is the first exercice of chapter 5 from the Book 'Bayesian Statistical Modeling with Stan, R and Python' of Kentaro Matsuura, 2023.

I am fitting a bayesian linear model in rstan and I tried to plot the residuals with two different methods but I do not obtain the same result and I do not understand why. For some context, the model is the following:

$$Y = b_1 + b_2 \cdot Sex + b_3 \cdot Income + \epsilon,$$ where $$Y$$ is the proportion of purchase during a year, Sex is the sex of the person who bought something, and Income is the income of this same person. $$b_1,b_2,b_3$$ are the coefficients of the linear model and $$\epsilon \sim \mathcal{N}(0,\sigma)$$.

For the first method, I computed the residuals $$\epsilon$$ from the fitted mean values and $$Y$$ values of the data: $$\epsilon_i = Y_i - \mu_i$$. The corresponding STAN code is:

// The input data is a vector 'y' of length 'N'.
data {
int<lower=0> N;
vector<lower = 0, upper = 1>[N] Y;
vector<lower = 0, upper = 1>[N] Sex;
vector<lower = 0>[N] Income;
}

// The parameters accepted by the model. Our model
// accepts two parameters 'mu' and 'sigma'.
parameters {
vector[3] b;
real<lower=0> sigma;
}

transformed parameters {
vector[N] mu = b[1] + b[2]*Sex[1:N] + b[3]*Income[1:N];
}

model {
Y[1:N] ~ normal(mu[1:N], sigma);
}


and the R-code

d <- read.csv("MatsuuraKenaro_Book/Chapter 5/data-shopping-1.csv")

data <- list(N = length(d[,1]), Y = d$$Y, Sex = d$$Sex, Income = d$Income) model <- cmdstan_model('Exercise5.1.stan') fit <- model$sample(
data=data, seed=123,
chains=4,
iter_warmup=500, #depend on trace plot; value between 200-1000
iter_sampling=2000, #for trying the model : 300-1000 iterations, for final model 2000-10000 iterations
parallel_chains=4,
save_warmup=TRUE
)

d.ms <- as.data.frame(fit$draws()) mu <- d.ms[,21:220] epsilon <- as.matrix(d$Y - mu)
hist(epsilon)


For the second method, I tried to generate $$\epsilon$$ directly with the STAN code. Here is my code snippets:

// The input data is a vector 'y' of length 'N'.
data {
int<lower=0> N;
vector<lower = 0, upper = 1>[N] Y;
vector<lower = 0, upper = 1>[N] Sex;
vector<lower = 0>[N] Income;
}

// The parameters accepted by the model. Our model
// accepts two parameters 'mu' and 'sigma'.
parameters {
vector[3] b;
real<lower=0> sigma;
}

transformed parameters {
vector[N] mu = b[1] + b[2]*Sex[1:N] + b[3]*Income[1:N];
}

model {
Y[1:N] ~ normal(mu[1:N], sigma);
}

generated quantities {
vector[N] epsilon = Y[1:N] - mu[1:N];
}


and the R-code

model2 <- cmdstan_model('Exercise5.2.stan')

fit2 <- model$sample( data=data, seed=123, chains=4, iter_warmup=500, #depend on trace plot; value between 200-1000 iter_sampling=2000, #for trying the model : 300-1000 iterations, for final model 2000-10000 iterations parallel_chains=4, save_warmup=TRUE ) d.ms2 <- as.data.frame(fit2$draws())
epsilon2 <- as.matrix(d.ms2[,421:620])
hist(epsilon2)


The resulting histogramm is shown below:

If you want to have a look at the data, here is the link to the dataset :dataset.

• Hello. Since you've provided the book title, it was easy to find that the book "Bayesian Statistical Modeling with Stan, R, and Python" has a GitHub repo with all the data here. Please update your question to replace the image of the data frame with a link to the repo and the data file. Commented Aug 29, 2023 at 10:03
• This is much more of a programming question than a statistics question. It might be better ask it on StackOverflow to get more help. I suspect the issue is with the first method. (The red histogram "looks" wrong because the estimate of sigma is 0.04 and the red histogram is way wider than that.) Commented Aug 29, 2023 at 17:31
• So fit$draws() is an array with size (iter_sampling, chains, N) = (2000, 4, 50) while d$Y is numeric with length N = 50. You flatten the array with as.data.frame and then subtract d$Y - mu. I guess that there's an index mismatch between d$Y and mu. You need to subtract d\$Y[i] from the corresponding draws[,,i]. Even better than learning how R flattens arrays may be learn how to use the posterior package to manipulate posterior draws. Commented Aug 29, 2023 at 17:32