3

R下面给出的是使用包分析可解析的 alpha 设计（alpha lattice design）的代码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)

我正在尝试使用Rpackage复制上述内容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 ^ 2asreml

4

1 回答 1

2

我对这种方差分区的东西不太熟悉，但我会试一试。

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(xy) = 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%）；我不知道这是否足以引起关注。

于 2016-08-01T18:37:53.177 回答