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

[H. Skaug]

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
> Users at admb-project.org
> http://lists.admb-project.org/mailman/listinfo/users