rm(list=ls()) # Removes all variables in memory
setwd("C:/VSNC/tuiwen")
alfalfa<-read.table("ALFALFA.TXT",header=TRUE,na.string='*')
head(alfalfa)
## ID Variety Block Resp
## 1 1 A 1 2.17
## 2 2 B 1 1.58
## 3 3 C 1 2.29
## 4 4 D 1 2.23
## 5 5 A 2 1.88
## 6 6 B 2 1.26
summary(alfalfa)
## ID Variety Block Resp
## Min. : 1.00 A : 6 Min. :1.0 Min. :0.880
## 1st Qu.:18.75 B : 6 1st Qu.:2.0 1st Qu.:1.310
## Median :36.50 C : 6 Median :3.5 Median :1.585
## Mean :36.50 D : 6 Mean :3.5 Mean :1.597
## 3rd Qu.:54.25 E : 6 3rd Qu.:5.0 3rd Qu.:1.820
## Max. :72.00 F : 6 Max. :6.0 Max. :2.340
## (Other):36
str(alfalfa)
## 'data.frame': 72 obs. of 4 variables:
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Variety: Factor w/ 12 levels "A","B","C","D",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ Block : int 1 1 1 1 2 2 2 2 3 3 ...
## $ Resp : num 2.17 1.58 2.29 2.23 1.88 1.26 1.6 2.01 1.62 1.22 ...
# Creating factors
alfalfa$Variety<-as.factor(alfalfa$Variety)
alfalfa$Block<-as.factor(alfalfa$Block)
str(alfalfa)
## 'data.frame': 72 obs. of 4 variables:
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Variety: Factor w/ 12 levels "A","B","C","D",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ Block : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ Resp : num 2.17 1.58 2.29 2.23 1.88 1.26 1.6 2.01 1.62 1.22 ...
# Analysis using ASReml
library(Matrix)
library(asreml)
## Warning: package 'asreml' was built under R version 3.5.1
# fit mixed model- Variety random
# Full model Base
alf.r<-asreml(fixed=Resp~Block,
random=~Variety,
data=alfalfa)
## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Mon Jul 1 11:09:47 2019
## LogLik Sigma2 DF wall cpu
## 1 48.7345 0.0619742 66 11:09:47 0.0
## 2 50.0211 0.0573191 66 11:09:47 0.0
## 3 51.1501 0.0525523 66 11:09:47 0.0
## 4 51.6975 0.0487491 66 11:09:47 0.0
## 5 51.7366 0.0477512 66 11:09:47 0.0
## 6 51.7370 0.0476536 66 11:09:47 0.0
plot(alf.r)
summary(alf.r)$varcomp
## component std.error z.ratio bound %ch
## Variety 0.02768043 0.015262120 1.813669 P 0
## units!R 0.04765360 0.009087279 5.243989 P 0
wald(alf.r)
## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: Resp[0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 40.936 859.04 < 2.2e-16 ***
## Block 5 4.150 87.08 < 2.2e-16 ***
## residual (MS) 0.048
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Heritability
(H2<-vpredict(alf.r, H2~V1/(V1+V2)))
## Estimate SE
## H2 0.3674359 0.1405692
# Heritability
H2<-0.02768043/(0.02768043+0.04765360) #Broad-sense Herita.
H2
## [1] 0.3674359
# Looking at significance variance components
rest0<-asreml(fixed=Resp~Block,data=alfalfa)
## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Mon Jul 1 11:09:50 2019
## LogLik Sigma2 DF wall cpu
## 1 44.8781 0.0753324 66 11:09:50 0.0
summary(rest0)$varcomp
## component std.error z.ratio bound %ch
## units!R 0.07533245 0.0131137 5.744563 P 0
# LRT - Model significant testing
lrt(alf.r,rest0)
## Likelihood ratio test(s) assuming nested random models.
## (See Self & Liang, 1987)
##
## df LR-statistic Pr(Chisq)
## alf.r/rest0 1 13.718 0.0001062 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(alf.r)$aic
## [1] -99.47393
## attr(,"parameters")
## [1] 2
summary(alf.r)$bic
## [1] -95.09462
## attr(,"parameters")
## [1] 2
summary(alf.r)$loglik
## [1] 51.73696