Wie Kovarianzmatrix für zufällige Effekte (BLUPs/bedingten Modi) von lme4

Question

Also, ich ein lineares gemischtes Modell mit zwei zufälligen Abschnitte in R ausgerüstet haben:

Y = X beta + Z b + e_i, 

wo b ~ MVN (0, Sigma); X und Z sind die Fixed- und Random-Effects-Modellmatrizen bzw. beta und b die Fixed-Effect-Parameter und Random-Effects-BLUPs/Conditional-Modes.

Ich möchte die zugrunde liegende Kovarianz-Matrix von b in die Hände bekommen, was in lme4 Paket nicht so trivial zu sein scheint. Sie können nur die Varianzen durch VarCorr erhalten, nicht die tatsächliche Korrelationsmatrix.

Gemäß one of the package vignettes (Seite 2) können Sie die Kovarianz von Beta berechnen: e_i * lambda * t(lambda). Und all diese Komponenten können Sie aus der Ausgabe von lme4 extrahieren.

Ich fragte mich, ob das der richtige Weg ist? Oder haben Sie andere Vorschläge?

Answer 1

Von ?ranef:

If ‚condvar‘ ist ‚TRUE‘ jeder der Datenrahmen ein Attribut hat genannte ‚postVar‘ ", die eine dreidimensionale Anordnung mit symmetrischen Flächen ist; Jede Seite enthält die Varianz-Kovarianz-Matrix für eine bestimmte Ebene des Gruppierungsfaktors. (Der Name dieses Attributs ist ein historisches Artefakt und kann in Zukunft irgend Punkt zu ‚condvar‘ geändert werden.)

ein Beispiel ein:

library(lme4) 
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy) 
rr <- ranef(fm1,condVar=TRUE) 

die Varianz Get -covariance Matrix unter den b Werte für den Intercept

pv <- attr(rr[[1]],"postVar") 
str(pv) 
##num [1:2, 1:2, 1:18] 145.71 -21.44 -21.44 5.31 145.71 ... 

das ist also ein 2x2x18 Array; Jede Schicht ist die Varianz-Kovarianz-Matrix unter dem bedingten Abschnitt und der Steigung für ein bestimmtes Subjekt (per Definition sind die Intercepts und die Steigungen für jedes Subjekt unabhängig von den Interzepten und Abhängen für alle anderen Subjekte).

Um dies zu einem Varianz-Kovarianzmatrix zu konvertieren (siehe getMethod("image",sig="dgTMatrix") ...)

library(Matrix) 
vc <- bdiag( ## make a block-diagonal matrix 
      lapply(
       ## split 3d array into a list of sub-matrices 
       split(pv,slice.index(pv,3)), 
        ## ... put them back into 2x2 matrices 
         matrix,2)) 
image(vc,sub="",xlab="",ylab="",useRaster=TRUE)