4

Given below is the code for analysis of a resolvable alpha design (alpha lattice design) using the R package 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)

I am trying to replicate the above using the R package 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)

How to compute the vblup and vblue (mean variance of difference) in lme4 equivalent to predict()$avsed ^ 2 of asreml ?

Kevin Wright
  • 2,397
  • 22
  • 29
Crops
  • 5,024
  • 5
  • 38
  • 65

1 Answers1

4

I'm not that familiar with this variance partitioning stuff, but I'll take a shot.

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

Get conditional variances of 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 ...

This is a 1x1x24 array (effectively just a vector of variances; we could collapse using c() if we needed to). They're not all the same, but they're pretty close ... I don't know whether they should all be identical (and this is a roundoff issue)

(uv <- unique(vv1))
## [1] 0.02887451 0.02885887 0.02885887

The relative variation is approximately 5.4e-4 ...

If these were all the same then the mean variance of a difference of any two would be just twice the variance (Var(x-y) = Var(x)+Var(y); by construction the BLUPs are all independent). I'm going to go ahead and use this.

vblup <- 2*mean(vv1)

For the model with gen fitted as a fixed effect, let's extract the variances of the parameters relating to genotypes (which are differences in the expected value from the first level):

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'm going to take the means of these values (not double the values, since these are already the variances of differences)

vblue <- mean(vv2)

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

The H^2 estimate looks right on, but the H^2c estimate is a little different (0.797 vs. 0.809, a 1.5% relative difference); I don't know if that is big enough to be of concern or not.

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453