下面给出的是使用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)
如何计算相当于asreml
的predict()$avsed ^ 2
的lme4
中的vblup
和vblue
(差异的平均方差)?
我不太熟悉这种方差划分的东西,但我会尝试一下。
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%);我不知道这是否足以引起人们的关注。