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.
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]
> 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
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