Fix up VarCorr

[pitakakariki]

Pivoting bug. Look at safe_chol in lme4 codebase?

[pitakakariki]

Add lots of tests and maybe error in VarCorr<- if VarCorr doesn't match afterwards.

[pitakakariki]

> V2
          [,1]      [,2]      [,3]
[1,] 0.1583598 0.1166801 0.1761762
[2,] 0.1166801 0.3506677 0.2832665
[3,] 0.1761762 0.2832665 0.3423178
> model <- makeLmer(y ~ x + x2 + (x + x2|g),
+                   fixef=b, VarCorr=V2,
+                   sigma=s, data=X)
> print(model)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ x + x2 + (x + x2 | g)
   Data: X
REML criterion at convergence: 6108.123
Random effects:
 Groups   Name        Std.Dev. Corr     
 g        (Intercept) 0.3457           
          x           0.5922   0.00     
          x2          0.5851   0.40 0.82
 Residual             1.1076           
Number of obs: 1920, groups:  g, 32
Fixed Effects:
(Intercept)            x           x2 
    -0.0166       0.1216       0.1942

[pitakakariki]

While you're at it make some more documentation for VarCorr.

[pitakakariki]

Note that VarCorr<- can also update sigma so that should be properly documented.

[hulahoopmagic]

Hi, I was wondering whether there is still a bug in how makeLmer uses the variance-covariance matrix? In the example you posted above, I don't feel like the variances and covariances in V2 correspond to the variances and covariances in model? For instance, shouldn't be the value of the first column and first row (0.1583598) be the squared standard deviation of the by-g random intercept (and hence be 0.1195085)? Further, the value of the second column and first row and first column and second row (0.1166801) should be zero if there is indeed a correlation of zero of the by-g random intercept and by-g random slope of x, right?

[pitakakariki]

Yes, the 30 April post is an example of the bug. Still working out how to fix it, but I have some workarounds if you're being affected by the bug.

[hulahoopmagic]

Ah I see! Yes, it seems like I am affected by this bug. I tried to figure out what’s wrong with my variance-covariance matrix, but probably it’s nothing wrong with it :) Hence, I am very interested in your workarounds!

[pitakakariki]

Where did yours come from? Is it from a previously fit lme4 model, or have you specified it by hand?

[hulahoopmagic]

I actually tried both, one that I made up and another one that I extracted from a lme4 model :)

[pitakakariki]

If you already have one in another model, that's the easiest option because you can just copy the theta parameter across:

newmodel@theta <- oldmodel@theta

If your variance-covariance matrix V is positive-definite then this should work too:

newmodel@theta <- t(chol(V))[lower.tri(V, diag=TRUE)]

[hulahoopmagic]

Wow, cool! Thanks so much, that works just fine :-)