Personal tools
You are here: Home Math Covariance in RE models

Covariance in RE models

In RE (random effects) models there are two types of parameters, x (parameter in the ordinary sense) and u (random effect). This documents tries to shed some light on how the variance of x and u are calculated. The covariance matrix of the x vector is based on the (marginal) likelihood, obtained via the Laplace approximation, and corresponds to the way covariance matrices are calculated in non-RE models in ADMB.

Theory

Here is what the user manual says about the variance of u (the manual talks about "theta" instead of "x"): usermanual.pdf A few words can be added to this.

The formula is based on the Law of total variance: http://en.wikipedia.org/wiki/Law_of_total_variance

In our context this says:

Var(u) = Ex[var(u|x)] + varx(E(u|theta)

The expectation "Ex" is obtained simply by inserting the point estimate of x into "var(u|x)". The second term is based on the "delta method" which is used elsewhere in ADMB, in combination with the covariance matrix of x (described above). Everything in these calculations are conditional on "data".

Example: simple hierarchical model

Consider the following simple Gaussian hierarchical model:

``` Prior on x: x = e1
u|x: u = x + e2
y|u: y = u + e3
// where e1, e2, e3 are all distributed N(0,1)```

R code for the covariance matrix

```S = matrix(0,3,3,row=c("x","u","y"),col=c("x","u","y"))
S[,]=1
S[2:3,2:3]=2
S[3,3]=3
S12_3 = S[1:2,1:2] - S[1:2,3]%*%solve(S[3,3])%*%S[3,1:2]```
We are interested in the conditional variance of x and u given data (y).
```> sqrt(diag(S12_3))
x         u
0.8164966 0.8164966
> cov2cor(S12_3)
x   u
x 1.0 0.5
u 0.5 1.0```

An implementation of this model in ADMB is:
``` DATA_SECTION
number y
!! y=10.0;
PARAMETER_SECTION
init_number x
random_effects_vector u(1,1)
objective_function_value f
PROCEDURE_SECTION
f = 0.0;
f -= -0.5*square(x);	 	// Prior on x: x = e1
f -= -0.5*square(u(1)-x);    // u|x: u = x + e2
f -= -0.5*square(y-u(1));    // y|u: y = u + e3
// where e1, e2, e3 are all distributed N(0,1); standard normal
GLOBALS_SECTION
#include "getbigs.cpp"```

NOTE: as of Nov 29 2012 you need to include "getbigs.cpp" due to a recently discovered bug in ADMB.

The result you get when you run ADMB matches those from R:

```D:\tmp\tmp>more simple_variance.cor
The logarithm of the determinant of the hessian = 0.405465
index   name    value      std dev       1
1   x 3.3335e+000 8.1650e-001   1.0000
2   u 6.6668e+000 8.1650e-001   0.5000  1.0000```