`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를 사용하여 위를 복제하는 것을 시도하고있다.

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

방법 asreml의 predict()$pred$avsed^2 등가 lme4에서 vblup 및 vblue (차의 평균 편차)을 계산하는 방법?

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

이들 모두가 동일하다면, 임의의 두 차이의 평균 분산은 단지 두 배의 분산 (Var (xy) = Var (x) + Var (y)) 일 것입니다. BLUPs는 모두 독립적입니다.) 나는 이것을 사용하겠습니다.

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 : 고정 효과로 끼워 gen과 모델의

vblup <- 2*mean(vv1) 

, 이제 (제 수준에서 예상되는 값의 차이 임) 유전형에 관한 파라미터의 차이를 추출하도록

vblue <- mean(vv2) 

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

,617,451 (이러한 차이는 이미 분산되어 있기 때문에, 하지 두 값) '은 이들 값을 이용 걸릴 것있어

H^2 추정치가 올바르게 보이지만 H^2c 추정치는 약간 다릅니다 (0.797 vs. 0.809, 1.5 % 상대 차이). 그게 큰 관심사인지 아닌지 나는 모른다.