1
$\begingroup$

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: In red the distribution of the residuals obtained from the first method and in black from the second method.

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

$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$
    – dipetkov
    Commented Aug 29, 2023 at 10:03
  • 1
    $\begingroup$ 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.) $\endgroup$
    – dipetkov
    Commented Aug 29, 2023 at 17:31
  • 1
    $\begingroup$ 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. $\endgroup$
    – dipetkov
    Commented Aug 29, 2023 at 17:32

0

Browse other questions tagged or ask your own question.