data {
//number of observations
int<lower=0> N_obs;
//number of groups
int<lower=0> Ngroups;
//number of columns in the design matrix
int K;
//the column containing the groups
int<lower=1> group[N_obs];
//number of observations per plate
vector[N_obs] n_t;
// response
vector[N_obs] y_obs;
//timevar
//vector[N_obs] timevar;
//design matrix
matrix [N_obs, K] X_obs;
//vector of S2
vector[N_obs] S2;
}
//transfored data
transformed data {
//observed data transformation
matrix[N_obs, K] Q_ast;
matrix[K, K] R_ast;
matrix[K, K] R_ast_inverse;
// thin and scale the QR decomposition
Q_ast = qr_thin_Q(X_obs) * sqrt(N_obs - 1);
R_ast = qr_thin_R(X_obs) / sqrt(N_obs - 1);
R_ast_inverse = inverse(R_ast);
}
// The parameters accepted by the model. Our model
parameters {
// regression coefficient for mean structure
real alpha; //intercept
vector[K] theta;
//vraiance
real<lower=0> sigma;
real<lower=0> nu;
real<lower=0> sigmaint; //variance for random intercept
vector[Ngroups] etaint; //reparameterization for the random effects
}
transformed parameters {
vector[Ngroups] ranint; //random intercepts
vector[N_obs] yhat;
//computing the random intercept
ranint = sigmaint * etaint;
//mean function
for (i in 1:N_obs)
yhat[i] = Q_ast[i] * theta + ranint[group[i]] + alpha;
}
// The model to be estimated.
model {
vector[N_obs] b2;
// priors
alpha ~ cauchy(0, 10); //prior for the intercept
sigma ~ scaled_inv_chi_square(nu, 0.4);
nu ~ gamma(2, 0.1);
etaint ~ normal(0, 1);
//constrain such that (n - 1)S^2 / sigma^2 is chi^2(n - 1)
for(i in 1:N_obs)
b2[i] = ((n_t[i] - 1) * S2[i]) / pow(sigma, 2);
//jacobian such that we sample from the target distribution for b2
for(i in 1:N_obs)
b2[i] ~ chi_square(n_t[i] - 1); //distribution for b2
target += log(n_t - 1) + log(S2) - 2*log(pow(sigma, 2));
// likelihood
for (n in 1:N_obs)
y_obs[n] ~ normal(yhat[n], sigma/pow(n_t[n], 0.5));
}
generated quantities {
vector[K] beta;
vector[N_obs] log_lik;
beta = R_ast_inverse * theta; // coefficients on x
//log-likelihood
for (n in 1:N_obs)
log_lik[n] = normal_lpdf(y_obs[n]| yhat[n], sigma/pow(n_t[n], 0.5));
}
I am building a model for \bar{y} \sim N(\mu, \frac{\sigma}{\sqrt{n}}), such that \frac{(n - 1) S^2}{\sigma^2} \sim \chi^2. I did this because \sigma^2 was of interest and the data contained S^2 and n - 1 for each row. The model converged without any problems, but whenever I try to apply the loo function, I get NaN for the relative efficiency of 3 observations, so I am unable to compute loo. Note, I am trying to compare this model with a model where I don’t estimate \sigma^2, but the standard error directly i.e. \bar{y} \sim N(\mu, \sigma_{ \bar{y} }).