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