「lme4」のBLUEまたはBLUPの差の平均分散

Question

Rパッケージasremlを使用して解決可能なアルファデザイン（アルファ格子デザイン）を分析するコードを以下に示します。

# 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) 

私はRパッケージlme4を使用して、上記複製しようとしています。 asremlのpredict()$pred$avsed^2に相当lme4にvblupとvblueを（差の分散を意味する）を計算する方法

# 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) 

？

Answer 1

私はこの分散パーティショニングのことに慣れていませんが、私はショットを撮ります。

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 

は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 ... 

これは1x1x24配列である（分散の効果だけベクトルを、我々は我々がする必要がある場合c()を使用して崩壊する可能性があり）。相対的な変動はおよそ私は、彼らががすべて同一である必要があります（これは丸め問題である）かどうかわからない...

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

彼らはすべて同じではないですが、彼らはかなり近いです5.4e-4 ...

これらがすべて同じ場合、任意の2つの差の平均分散は、（Var（xy）= Var（x）+ Var（y））の2倍になります。 BLUPはすべて独立しています）。私は先に進んでこれを使うつもりです。固定効果として嵌め込まgen内蔵モデル

vblup <- 2*mean(vv1) 

、のは、（最初​​のレベルからの期待値の差である）の遺伝子型に関連するパラメータの差異を抽出してみましょう：

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 「

vblue <- mean(vv2) 

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

（これらはすでに違いの分散であるため、 ないダブル値）これらの値の平均を取るつもりメートル

H^2の見積もりが正しいと思われますが、H^2cの見積もりは少し異なります（0.797 vs. 0.809、相対差1.5％）。それが懸念されるほど大きいかどうかはわかりません。