Mittlere Varianz eines Unterschieds von BLAUEN oder BLUPs in "lme4"

Question

Unten ist der Code für die Analyse eines auflösbaren Alpha-Designs (Alpha-Gitter-Design) unter Verwendung des R Pakets asreml angegeben.

# load the data 
library(agridat) 
data(john.alpha) 
dat <- john.alpha 

# load asreml 
library(asreml) 

# model1 - random `gen` 
#---------------------- 
# fitting the model 
model1 <- asreml(yield ~ 1 + rep, data=dat, random=~ gen + rep:block) 
# variance due to `gen` 
sg2 <- summary(model1)$varcomp[1,'component'] 
# mean variance of a difference of two BLUPs 
vblup <- predict(model1 , classify="gen")$pred$avsed^2 

# model2 - fixed `gen` 
#---------------------- 
model2 <- asreml(yield ~ 1 + gen + rep, data=dat, random = ~ rep:block) 
# mean variance of a difference of two adjusted treatment means (BLUE) 
vblue <- predict(model2 , classify="gen")$pred$avsed^2 

# H^2 = .803 
sg2/(sg2 + vblue/2) 
# H^2c = .809 
1-(vblup/2/sg2) 

Ich versuche, die oben mit dem R Paket lme4 zu replizieren.

# model1 - random `gen` 
#---------------------- 
# fitting the model 
model1 <- lmer(yield ~ 1 + (1|gen) + rep + (1|rep:block), dat) 
# variance due to `gen` 
varcomp <- VarCorr(model1) 
varcomp <- data.frame(print(varcomp, comp = "Variance")) 
sg2 <- varcomp[varcomp$grp == "gen",]$vcov 

# model2 - fixed `gen` 
#---------------------- 
model2 <- lmer(yield ~ 1 + gen + rep + (1|rep:block), dat) 

Wie die vblup und vblue (mittlere Varianz der Differenz) in lme4 entspricht predict()$pred$avsed^2 von asreml berechnen?

Answer 1

Ich bin nicht vertraut mit dieser Varianz Partitionierung, aber ich werde eine Aufnahme machen.

library(lme4) 
model1 <- lmer(yield ~ 1 + rep + (1|gen) + (1|rep:block), john.alpha) 
model2 <- update(model1, . ~ . + gen - (1|gen)) 

## variance due to `gen` 
sg2 <- c(VarCorr(model1)[["gen"]]) ## 0.142902 

Erhalten bedingte Abweichungen von BLUPs:

rr1 <- ranef(model1,condVar=TRUE) 
vv1 <- attr(rr$gen,"postVar") 
str(vv1) 
## num [1, 1, 1:24] 0.0289 0.0289 0.0289 0.0289 0.0289 ... 

Dies ist ein 1x1x24 Array (effektiv nur ein Vektor der Varianzen, wir verwenden c() kollabieren könnten, wenn wir nach Bedarf). Sie sind nicht alle gleich, aber sie sind ziemlich nah dran ... Ich weiß nicht, ob sie sollte alle identisch sein (und dies ist eine Abrundungs ​​Ausgabe)

(uv <- unique(vv1)) 
## [1] 0.02887451 0.02885887 0.02885887 

Die relative Abweichung beträgt ca. 5.4e-4 ...

Wenn diese alle gleich wären, dann wäre die mittlere Varianz einer Differenz von zwei nur das Doppelte der Varianz (Var (xy) = Var (x) + Var (y); Aufbau der BLUPs sind alle unabhängig). Ich werde weitermachen und das benutzen.

vblup <- 2*mean(vv1) 

Für das Modell mit gen als festen Effekt ausgestattet, lassen Sie sich die Varianzen des Parameters Genotypen in Bezug extrahiert (die von der ersten Ebene Unterschiede in dem erwarteten Wert sind):

vv2 <- diag(vcov(model2))[-(1:3)] 
summary(vv2) 
## 
## Min. 1st Qu. Median Mean 3rd Qu. Max. 
## 0.06631 0.06678 0.07189 0.07013 0.07246 0.07286 

I bin dabei, die Mittel dieser Werte zu übernehmen (nicht das doppelte der Werte, da diese bereits die Varianzen der Unterschiede)

vblue <- mean(vv2) 

sg2/(sg2+vblue/2) ## 0.8029779 
1-(vblup/2/sg2)  ## 0.7979965 

Die Schätzung H^2 sieht gut aus, aber die H^2c Schätzung ist ein wenig anders (0,797 vs. 0,809, eine 1,5% relative Differenz); Ich weiß nicht, ob das groß genug ist, um es zu beunruhigen oder nicht.