# [ADMB Users] difference between ADMB-RE and R/mgcv in SEs for smoother coefficients in a GAM fitted by maximum likelihood

H. Skaug hskaug at gmail.com
Thu Nov 29 13:10:53 PST 2012

```I think there may be a problem with the correlation
between fixed and random effects. I modified Dave's
example a bit to cast it in probabilistic terms where
I know what the covariance matrix should be:

********ADMB program with explanation of hierarchical model:
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 N(0,1)

******* Correlation matrix where the correlation has the wrong  sign
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

Comment: it is counter intuitive that x and u should
be negatively correlated.

***** Calculations in R to find the covarinace 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]  #
Conditional covariance (x,u) given 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

Conclusion: sd's match ADMB, but corrrelation has oposite sign.

Hans

On Thu, Nov 29, 2012 at 6:53 PM, dave fournier <davef at otter-rsch.com> wrote:
> A very simple example shows how the method works.
>
>
> DATA_SECTION
>
> PARAMETER_SECTION
>   init_number x
>   random_effects_vector u(1,1)
>   objective_function_value f
>   sdreport_number sd1
>   sdreport_number sd2
> PROCEDURE_SECTION
>
>   f= square(x-10) + 0.5*square(u(1)) + square(u(1)-x);
>   sd1=u(1)-x;
>   sd2=x-u(1);
>
>
> The logarithm of the determinant of the hessian = 0.980829
>  index   name    value      std.dev       1       2       3
>      1   x    7.5002e+00 6.1237e-01   1.0000
>      2   u    5.0001e+00 5.7735e-01   0.0000  1.0000
>      3   sd1 -2.5001e+00 8.4163e-01  -0.7276  0.6860  1.0000
>      4   sd2  2.5001e+00 8.4163e-01   0.7276 -0.6860 -1.0000  1.0000
>
> I guess the question is whether this is reasonable for sd1. It seems a bit
> large to me.
>
> _______________________________________________
> Users mailing list