r - `lme4` 中 BLUE 或 BLUP 差异的平均方差

Question

下面给出的是使用 R 分析可解析 alpha 设计(alpha 点阵设计)的代码。套餐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")$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")$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)

如何计算vblup和vblue (差异的平均方差)lme4相当于predict()$avsed ^ 2的asreml ？

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

获取 BLUP 的条件方差:

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(x-y) = Var(x)+Var(y)；通过构造，BLUP 都是独立的)。我将继续使用它。

vblup <- 2*mean(vv1)

对于将gen拟合为固定效应的模型，让我们提取与基因型相关的参数的方差(即与第一级的期望值的差异):

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

我将采用这些值的平均值(不是将值加倍，因为这些已经是差异的方差)

vblue <- mean(vv2)

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

H^2 估计值看起来正确，但 H^2c 估计值略有不同(0.797 与 0.809，相对差异为 1.5%)；我不知道这是否足够大到值得关注。