1. Home
2. Videos: ASReml-R

Videos: ASReml-R

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.

Split-plot design

A split-plot design is recommended when a factor that you want to analyse can’t be changed as easily as your other factors.

R script used in this video

library(asreml)

# load the ‘oats’ data set – an internal data set in the ASReml package
data(oats)

# display the data set using View function
View(oats)

# use str function to investigate structure
str(oats)

# Split Plot Design analysis model
oats.asr <- asreml(fixed = yield ~ Variety + Nitrogen + Variety:Nitrogen,
random = ~ Blocks + Blocks:Wplots,
data = oats)

# look at the residuals
plot(oats.asr)

# summarise the variance components
summary(oats.asr)\$varcomp

# fixed(or random) effects
summary(oats.asr, coef = TRUE)\$coef.fixed
summary(oats.asr, coef = TRUE)\$coef.random

# denDF = “none” and ssType = “incremental” (the defaults)
wald(oats.asr, denDF = “default”)

#Prediction
oatsN.pv <- predict(oats.asr, classify = “Nitrogen”, sed = T)

str(oatsN.pv)

# predicted means
oatsN.pv\$pvals

# full matrix of SEDs
oatsN.pv\$sed

# minimum, mean and maximum SED
oatsN.pv\$avsed

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.

R script used in this video

library(asreml)

#check internal data in asreml
data(package=”asreml”)

#data input
data(“harvey”)
View(harvey)

#check the pedigree data

# 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

#the first 15 E-BLUP estimates

#######################
#2) sire model

adg1.asr <- asreml(y3 ~ Line, random = ~Sire, data = harvey)

############################################################
#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)

#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)

########################
# 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)

#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)

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.

R script used in this video

library(asreml)

data(“shf”)

#print the first few rows

#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

Updated on May 15, 2019