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.
Random coefficient regression using ASReml-R 4
Random coefficient regression can be used when we want to explore the relationship between a response variable (y) and a continuous explanatory variable (x) and we have repeated measurements of x and y on individual subjects.
Data source: RatGrowth.txt
https://vsnidownloads.s3.eu-west-2.amazonaws.com/130fd7da93a9ad8e9525104d753b025c3d62a178/RatGrowth.txt
R script
##################################
# RANDOM COEFFICIENT REGRESSION
# Using Rat Growth Example
##################################
# Description:
# These data originated from an experiment to study the effect of drugs on the growth rates of rats.
# There were three treatment groups (control, thyroxin and thiouracil) and the weight of the rats
# were measured weekly.
# Reference:
# Box, G. (1950). Problems in the Analysis of Growth and Wear Curves. Biometrics 6(4), 362-389.
# Set working directory
setwd(“C:/Users/jane/Desktop/RatGrowth”)
# Load data and explore its structure
RatGrowth<-read.table(“RatGrowth.txt”,header=TRUE)
head(RatGrowth)
tail(RatGrowth)
str(RatGrowth)
# Convert rat to a factor
RatGrowth$rat<-as.factor(RatGrowth$rat)
str(RatGrowth)
# Plot the data
library(lattice)
xyplot(weight~time|drug,data=RatGrowth,type=”l”,groups=rat)
# Calculate time centered and time centered squared
RatGrowth$timec<-RatGrowth$time-mean(RatGrowth$time)
RatGrowth$timecSQ<-RatGrowth$timec^2
head(RatGrowth)
# Data frame with thiouracil
RatGrowth_thiouracil<-RatGrowth[RatGrowth$drug %in% “thiouracil”,]
# Loading asreml library
library(asreml)
# Random coefficient regression model for thiouracil drug
RCModel_thiouracil<-asreml(fixed=weight~timec+timecSQ,
random=~str(~rat+timec:rat+timecSQ:rat,~corgh(3):id(rat)),
data=RatGrowth_thiouracil,maxit=20)
plot(RCModel_thiouracil) # residual plots
summary(RCModel_thiouracil)$varcomp # variance components
summary(RCModel_thiouracil,coef=TRUE)$coef.fixed # fixed effects
summary(RCModel_thiouracil,coef=TRUE)$coef.random # random effects
wald.asreml(RCModel_thiouracil,denDF=’algebraic’)$Wald # Wald tests
# Random coefficient regression model for all data
RCModel<-asreml(fixed=weight~timec+timecSQ+drug+drug:timec+drug:timecSQ,
random=~str(~rat+timec:rat+timecSQ:rat,~corgh(3):id(rat)),
data=RatGrowth)
plot(RCModel) # residual plots
summary(RCModel)$varcomp # variance components
summary(RCModel,coef=TRUE)$coef.fixed # fixed effects
summary(RCModel,coef=TRUE)$coef.random # random effects
wald.asreml(RCModel,denDF=’algebraic’)$Wald # Wald tests
# Plot of the predicted means
# Note: The values of timec and timecSQ change in parallel.
predTime<-predict(RCModel,classify=”drug:timec:timecSQ”,
levels=list(drug=c( 1,1,1,1,1, 2,2,2,2,2, 3,3,3,3,3),
timec=c( -2,-1,0,1,2, -2,-1,0,1,2, -2,-1,0,1,2),
timecSQ=c(4,1,0,1,4, 4,1,0,1,4, 4,1,0,1,4)),
parallel=TRUE,sed=FALSE)$pvals
head(predTime)
tail(predTime)
xyplot(predicted.value~(timec+mean(RatGrowth$time)),groups=drug,data=predTime,
type=”l”,lwd=2,auto.key=TRUE)
Creating and fitting a mixed effects model in ASReml-R4
GBLUP using ASReml-R 4
Genomic Best Linear Unbiased Prediction, or GBLUP, is a genomic selection method that uses genetic relationships derived from molecular markers to estimate the genetic merit of an individual.
GBLUP_Apple.R
https://downloads.vsni.digital/642012ca469933c905f8f0741c04057f902008bb/GBLUP_Apple.R
G-Inv.txt
https://downloads.vsni.digital/8089e1718782adb8dc5d8b4cd1128354107c12ee/G-Inv.txt
PhenotApple.txt
https://downloads.vsni.digital/3bff07343d91497f08ea5260eac1f1ca19a8fd2b/PhenotApple.txt
GenotApple.txt
https://downloads.vsni.digital/fe712a47e75d2aa1eb6b94ac6cf44434e1dace74/GenotApple.txt
Repeated measures
A repeated measures analysis is needed when a trait is measured on the same individual over time.
# load the asreml package
library(asreml)
# print the first and last few rows
head(grassUV)
tail(grassUV)
library(ggplot2)
ggplot(data = grassUV,
aes(x = Time, y = y, group = Plant, colour=Plant)) +
geom_line() + facet_grid(. ~ Tmt)
grass
####### Uniform Model
# Univariate
grassUV.uniform <- asreml(y ~ Tmt + Time + Tmt:Time,
residual = ~Plant:cor(Time),
data =grassUV)
summary(grassUV.uniform)$bic
# Multivariate
grass.uniform <- asreml(cbind(y1, y3, y5, y7, y10) ~ Tmt + trait + Tmt:trait,
residual = ~Plant:cor(trait),
data = grass)
summary(grass.uniform)$bic
####### Power Model
# Convert to Univariate
grassConvert<-reshape(grass, direction=”long”, varying=3:7, sep = “”)
grassConvert<-grassConvert[order(grassConvert$Plant,grassConvert$time),]
grassConvert$time<-factor(grassConvert$time)
# Univariate
grassUV.power <- asreml(y ~ Tmt + Time + Tmt:Time,
residual = ~Plant:exp(Time),
data =grassUV)
summary(grassUV.power)$bic
plot(grassUV$Time,grassUV.power$resid)
####### Heterogeneous Power Model
grassUV.hetpower <- asreml(y ~ Tmt + Time + Tmt:Time,
residual = ~Plant:exph(Time),
data =grassUV)
summary(grassUV.hetpower)$bic
####### Antedependence Model
# Univariate
grassUV.ante <- asreml(y ~ Tmt + Time + Tmt:Time,
residual = ~Plant:ante(Time),
data =grassUV)
summary(grassUV.ante)$bic
# Multivariate
grass.ante <- asreml(cbind(y1, y3, y5, y7, y10) ~ trait + Tmt + trait:Tmt,
residual = ~Plant:ante(trait),
data = grass)
summary(grass.ante)$bic
####### Unstructured Model
# Univariate
grassUV.us <- asreml(y ~ Tmt + Time + Tmt:Time,
residual = ~Plant:us(Time),
data =grassUV)
summary(grassUV.us)$bic
# Multivariate
grass.us <- asreml(cbind(y1, y3, y5, y7, y10) ~ trait + Tmt + trait:Tmt,
residual = ~Plant:us(trait),
data = grass)
summary(grass.us)$bic
####### BICs
summary(grassUV.uniform)$bic
summary(grassUV.power)$bic
summary(grassUV.hetpower)$bic
summary(grassUV.ante)$bic
summary(grassUV.us)$bic
####### Wald tests for chosen model
# Univariate
wald(grassUV.ante)
# Multivariate
wald(grass.ante)
####### Predicted treatment by time means
# Univariate
predict(grassUV.ante, classify = “Tmt:Time”)
# Multivariate
predict(grass.ante, classify = “Tmt:trait”)
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.
# load the asreml package
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.
#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