`lme4`中BLUEs或BLUPs差异的均方差

Question

下面是使用 R 包 asreml 进行可解析 alpha 设计（alpha 格点设计）分析的代码。

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

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

- Crops

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

获取最优线性无偏预测值的条件方差：

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）;按照构造，BLUPs都是独立的）。 我将使用它。

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 vs 0.809，相对差异约为1.5%）; 我不知道这是否足够引起关注。

- Ben Bolker