Hi everyone! I am new to Stan, and I am trying to estimate variance components in a hierarchical linear model. I have some example data, and I have the outputs of two different statistical software packages, the results of which I am trying to match. So far, I have been terribly unsuccessful.
I am trying to fit the model y_{ijk} \sim \mu + \alpha_i + \beta_j + \alpha\beta_{i,j} + \epsilon, where \alpha_i \sim \mathcal{N}(0,\sigma^2_\alpha), \beta_i \sim \mathcal{N}(0,\sigma^2_\beta), \alpha \beta _i \sim \mathcal{N}(0,\sigma^2_{\alpha\beta}), \epsilon_{ijk} \sim \mathcal{N}(0,\sigma^2_e)
Both of the other statistical packages I am using give both REML and “Bayesian” estimates of these parameters, as shown below. I am attempting to recreate – or at least get close to – the Bayesian estimates. (However, getting the REML estimates would be very beneficial to my understanding)
Software 1:
REML:
sig1 = .0514553
sig2 = 1.0864466
sig3 = 0.0000000
sig4 = .0399733
Bayesian:
sig1 = .024646
sig2 = .884780
sig3 = .00201177
sig4 = .0382738
Software 2:
REML:
sig1 = .05146
sig2 = 1.0864
sig3 = 0
sig4 = .03997
Bayesian:
sig1 = .024795
sig2 = .89201
sig3 = .002013
sig4 = .03867
I have seen some documentation about how to use covariance matrices, but it appears that such matrices must be positive definite, and the software I am trying to replicate reports that “The Estimated G Matrix is not positive definite.”
How can I estimate these components of variance?
Stan Code:
data {
int LENGTH;
vector[LENGTH] x1;
vector[LENGTH] x2;
vector[LENGTH] Y;
}
parameters {
real mu;
real a;
real b;
real ab;
real eps;
real<lower=0> sig_a;
real<lower=0> sig_b;
real<lower=0> sig_ab;
real<lower=0> sig_e;
}
model {
Y ~ normal(mu + a*x1 + b*x2 + ab * (x1 .* x2) + eps, sig_a + sig_b + sig_ab + sig_e);
a ~ normal(0,sig_a);
b ~ normal(0,sig_b);
ab ~ normal(0,sig_ab);
eps ~ normal(0,sig_e);
}
Data:
x1,x2,y
1,1,0.29
1,1,0.41
1,1,0.64
2,1,0.08
2,1,0.25
2,1,0.07
3,1,0.04
3,1,-0.11
3,1,-0.15
1,2,-0.56
1,2,-0.68
1,2,-0.58
2,2,-0.47
2,2,-1.22
2,2,-0.68
3,2,-1.38
3,2,-1.13
3,2,-0.96
1,3,1.34
1,3,1.17
1,3,1.27
2,3,1.19
2,3,0.94
2,3,1.34
3,3,0.88
3,3,1.09
3,3,0.67
1,4,0.47
1,4,0.5
1,4,0.64
2,4,0.01
2,4,1.03
2,4,0.2
3,4,0.14
3,4,0.2
3,4,0.11
1,5,-0.8
1,5,-0.92
1,5,-0.84
2,5,-0.56
2,5,-1.2
2,5,-1.28
3,5,-1.46
3,5,-1.07
3,5,-1.45
1,6,0.02
1,6,-0.11
1,6,-0.21
2,6,-0.2
2,6,0.22
2,6,0.06
3,6,-0.29
3,6,-0.67
3,6,-0.49
1,7,0.59
1,7,0.75
1,7,0.66
2,7,0.47
2,7,0.55
2,7,0.83
3,7,0.02
3,7,0.01
3,7,0.21
1,8,-0.31
1,8,-0.2
1,8,-0.17
2,8,-0.63
2,8,0.08
2,8,-0.34
3,8,-0.46
3,8,-0.56
3,8,-0.49
1,9,2.26
1,9,1.99
1,9,2.01
2,9,1.8
2,9,2.12
2,9,2.19
3,9,1.77
3,9,1.45
3,9,1.87
1,10,-1.36
1,10,-1.25
1,10,-1.31
2,10,-1.68
2,10,-1.62
2,10,-1.5
3,10,-1.49
3,10,-1.77
3,10,-2.16