These tutorials use the real-life data sets included in the ASReml-R package to help you explore and learn how to use the software.
Specifying variance structures using genetic and other relationships
Learn how to include known variance structures, such as pedigree information, into a linear mixed model analysis.
#load asreml package
library(asreml)
#check internal data in asreml
data(package=”asreml”)
#data input
data(“harvey”)
View(harvey)
#check the pedigree data
head(harvey.ped,20)
# generate A inverse and check it
harvey.ainv <- ainverse(harvey.ped)
# check the attribute rowNames
attr(harvey.ainv, “rowNames”)
harvey.ainv[1:20, ]
# ASReml-R3
# harvey.ainv <- asreml.Ainverse(harvey.ped)$ginv
# attr(harvey.ainv,”rowNames”)
########################
#1) fitting animal model
adg0.asr <- asreml(y3 ~ Line, random = ~vm(Calf, harvey.ainv), data = harvey)
#in Version 3
# adg.asr <- asreml(y3 ~ Line,
# random = ~ ped(Calf, var=T),
# ginverse=list(Calf=harvey.ainv), data=harvey)
#variance component
summary(adg0.asr)$varcomp
#the first 15 E-BLUP estimates
head(adg0.asr$coef$random, 15)
#######################
#2) sire model
adg1.asr <- asreml(y3 ~ Line, random = ~Sire, data = harvey)
summary(adg1.asr)$varcomp
head(adg1.asr$coef$random, 9)
############################################################
#3)a. specified a general inverse matrix for Sire as 0.5I9
# define a general inverse matrix for Sire as 0.5I9
sire.giv <- data.frame(row = seq(1, 9), column = seq(1, 9), value = rep(0.5, 9))
sire.giv
#attribute INVERSE set to TRUE
attr(sire.giv, “INVERSE”) <- TRUE
#rowNames attribute
attr(sire.giv, “rowNames”) <- paste(“Sire_”, seq(1, 9), sep = “”)
#asreml.options(gammaPar = TRUE)
adg2.asr <- asreml(y3 ~ Line, random = ~vm(Sire, sire.giv), data = harvey)
summary(adg2.asr)$varcomp
head(adg2.asr$coef$random, 9)
#in Version 3
# adg2.asr <- asreml(y3 ~ Line,
# random = ~ giv(Sire, var=T),
# ginverse=list(Sire=sire.giv),
# data=harvey)
##########################################################
#3)b. alternative define a general matrix for Sire as 2I9
sire.mat <- data.frame(row = seq(1, 9), column = seq(1, 9), value = rep(2, 9))
sire.mat <- as.matrix(sire.mat)
attr(sire.mat, “rowNames”) <- paste(“Sire_”, seq(1, 9), sep = “”)
adg3.asr <- asreml(y3 ~ Line, random = ~vm(Sire, sire.mat), data = harvey)
summary(adg3.asr)$varcomp
########################
# 4)source is not p.d.
# a) source is p.s.d.
# the relationship matrix is read in as a vector
ped.vec <- c(2, 0, 2, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0,
0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
attr(ped.vec, “rowNames”) <- paste(“Sire_”, seq(1, 10), sep = “”)
adg4.asr <- asreml(y3 ~ Line, random = ~vm(Sire, ped.vec, “PSD”), data = harvey)
summary(adg4.asr)$varcomp
#b) source is p.z.n.d.
spped.mat <- data.frame(Row = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
Column = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 9),
Value = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1))
spped.mat <- as.matrix(spped.mat)
spped.mat
attr(spped.mat, “rowNames”) <- paste(“Sire_”, seq(1, 10), sep = “”)
adg5.asr <- asreml(y3 ~ Line, random = ~vm(Sire, spped.mat, “NSD”), data = harvey)
summary(adg5.asr)$varcomp
Spatial analysis of a field trial
A spatial analysis can improve the estimation of treatment effects by modelling more accurately the spatial distribution of error effects in a field trial.
#load the asreml package
library(asreml)
#load the data
data(“shf”)
#print the first few rows
head(shf)
#switch to gamma parameterization for estimation
asreml.options(gammaPar = TRUE)
#MODEL 1
#fit the separable autoregressive error model
barley1.asr <- asreml(yield ~ Variety, residual = ~ar1v(Row):ar1(Column), data = shf)
#REML log-likelihood, random components and Wald statistics from the fit
barley1.asr$loglik
summary(barley1.asr)$varcomp
wald(barley1.asr)
#plot the variogram
plot(varioGram(barley1.asr))
#MODEL 2
# an extension to this model includes a measurement error or nugget effect term
barley2.asr <- asreml(yield ~ Variety, random = ~idv(units), residual =
~ar1v(Row):ar1(Column), data = shf)
#REML log-likelihood, random components and Wald statistics from the fit
barley2.asr$loglik
summary(barley2.asr)$varcomp
wald(barley2.asr)
#log-likelihood ratio test to compare the two spatial models
lrt (barley1.asr, barley2.asr)
#MODEL 3
#incomplete block analysis (with recovery of inter-block information)
barley3.asr <- asreml(yield ~ Variety, random = ~Rep + RowBlk + ColBlk, data = shf)
#REML log-likelihood, random components and Wald statistics from the fit
barley3.asr$loglik
summary(barley3.asr)$varcomp
wald(barley3.asr)
#compare the 3 models using BIC
summary(barley3.asr)$bic
summary(barley2.asr)$bic
summary(barley1.asr)$bic
#PREDICTION
#predict the Variety means
barley2.pv <- predict(barley2.asr, classify = “Variety”)
#Variety means with standard errors
barley2.pv$pvals
#average standard error of difference between predicted means
barley2.pv$avsed