Subject-wise gradient and Hessian

Dear all,

I just use these simple codes to obtain gradient and Hessian of the fitted model to Orthodont data (for example) for different parameters,

fm1 <- lmer(scale(distance)~Sex*I(scale(age))+(I(scale(age))|Subject), data=Orthodont)
fm1 at optinfo$derivs$gradient

[1] 1.578826e-07 1.200391e-06 4.504841e-08

fm1 at optinfo$derivs$Hessian

[,1]       [,2]       [,3]
[1,]  30.766258  -4.868289 -13.459066
[2,]  -4.868289 128.743919   1.404892
[3,] -13.459066   1.404892  58.653149

So, it seems I have seven parameters, 4 fixed effects, 2 random effects, and one error variance, he so why I get a 3 \times 3 Hessian and 3 \times 1 gradient?

I also tried to use:

fm1_update <- update(fm1,devFunOnly=TRUE)
params <- unlist(getME(fm1,c("theta","beta")))
grad(fm1_update,params)

(7!=3)
Error: theta size mismatch

But you see the error I got, it only works fine when I just include "theta" and not the fixed effects "beta"

I wonder if anyone have an idea about how I can obtain gradient and Hessian for all the parameters in the model.

Afterward, I would try to find a way to calculate these quantities by subject. That would be probably possible by fitting a new model for each subject, but starting from fitted values and just for one iteration. But that is the next step after I know how I can obtain gradient and Hessian for each parameter at the first place.

Thanks!
Vahid.