View source for BLUP

[[:Category:R]]

==r-BLUP==

# computation of BLUP's, using seedling dry weight data as example

library(lme4)

d <- read.table('SeedlingDryWeight2.txt',as.is=T, header=T)

# create factors
d$g <- as.factor(d$Genotype)
d$t <- as.factor(d$Tray)
d$s <- as.factor(d$Seedling)

# fit the mixed effects model model
temp <- lmer(SeedlingWeight~g+(1|t),data=d)

# BLUPs
blup <- ranef(temp)

# look at ranef() output to see its structure
str(blup)

# there are options to ranef(drop=T) to simplify structure of output

# compare BLUP to genotype averages
gmean <- tapply(d$SeedlingWeight,d$t,mean)

plot(gmean,blup$t$'(Intercept)')

# BLUPs + fixed effects
#  there is no predict() for an lmer object
#  these are the equivalents

betau <- coef(temp)
# values are beta + u for each row

trayX <- cbind(rep(1,8), rep(c(0,1),c(4,4)))
# model.matrix with one row per tray
pred.tray <- apply(betau$t*trayX,1,sum)

plot(gmean, pred.tray)