Hello,

I am trying to recreate certain L Cholesky factors used in the calculation
of the profiled deviance but I seem to be unable to get exact results, even
when reusing lme4-provided intermediate results. I want to ask if that is
due only to different numerical procedure used or is actually because of
something different happening.
I first came across this problem with relative simple examples:

fm01 <- lmer(Yield ~ 1 + (1|Batch),Dyestuff, verbose=F)

L_fm01_native <- fm01@L
#Implement Cholesky factorization as described in "Penalized least squares
versus generalized least squares representations of linear mixed models"
Eq. 9
Lambda_fm01 <- with(expand(fm01), S %*% T)
Ut_fm01 <- crossprod(Lambda_fm01, fm01@Zt)
L_fm01_external <- Cholesky(tcrossprod(Ut_fm01), LDL = FALSE, Imult = 1)

where the difference
max(abs(L_fm01_external@x - L_fm01_native@x))
[1] 0.0001942293
is small enough compared to the values of x to probably assume it is due to
using difference numerical packages. (Permutation information is the same)

Nevertheless slightly more involved examples:

fm03 <- lmer(strength ~ 1 + (1|batch) + (0+ batch|sample), Pastes,
verbose=F)

L_fm03_native <- fm03@L
Lambda_fm03 <- with(expand(fm03), S %*% T)
Ut_fm03 <- crossprod(Lambda_fm03, fm03@Zt)
L_fm03_external <- Cholesky(tcrossprod(Ut_fm03), LDL = FALSE, Imult = 1)

yield higher differences
max(abs(L_fm03_external@x - L_fm03_native@x))
[1] 0.7636072
that are quite substantial, this being easily seen even by just looking at
the vectors x.

tail(L_fm03_native@x)
[1]  1.0696245 -0.1715797 -0.8431633  1.0974082  0.9148627  3.6455565
tail(L_fm03_external@x)
[1]  1.00046651 -0.01618789 -0.07955603  1.13172097  1.21935069  4.10916295

Is there an obvious way why this should happen?
I was not worried with the first example but the differences in Lambda
estimation in fm03 do result in noticeable differences in final estimates.
In "Computational Methods for mixed models" it is mentioned (right after
Eq.13) that "the factor L is only determined up to changes in the sign of
the diagonal elements", but aside this I can not justify this difference
(and I am convinced if that is the reason anyway).

I am using package versions lme4_0.999999-0 and Matrix_1.0-11 respectively.

Thank you very much in advance.

Pantelis Hadjipantelis

	[[alternative HTML version deleted]]

