inst/doc/tutorial.R

In agricolae: Statistical Procedures for Agricultural Research

### R code from vignette source 'tutorial.Rnw'

###################################################
### code chunk number 1: agricolae
###################################################
options(keep.source = TRUE, width = 60)
qp <- packageDescription("agricolae")


###################################################
### code chunk number 2: tutorial.Rnw:29-30
###################################################
rm(list=ls())


###################################################
### code chunk number 3: tutorial.Rnw:51-52
###################################################
library(agricolae)


###################################################
### code chunk number 4: tutorial.Rnw:58-60 (eval = FALSE)
###################################################
## help(package="agricolae")
## help(waller.test)


###################################################
### code chunk number 5: tutorial.Rnw:75-79
###################################################
detach(package:agricolae) # detach package agricole
library(agricolae) # Load the package to the memory
designs<-apropos("design")
designs[substr(designs,1,6)=="design"]


###################################################
### code chunk number 6: tutorial.Rnw:88-89
###################################################
library(agricolae) # Load the package to the memory: 


###################################################
### code chunk number 7: tutorial.Rnw:103-106
###################################################
A<-as.data.frame(data(package="agricolae")$results[,3:4])
A[,2]<-paste(substr(A[,2],1,35),"..",sep=".")
head(A)


###################################################
### code chunk number 8: tutorial.Rnw:117-120
###################################################
weight<-c( 68, 53, 69.5, 55, 71, 63, 76.5, 65.5, 69, 75, 76, 57, 70.5, 71.5, 56, 81.5,
           69, 59, 67.5, 61, 68, 59.5, 56.5, 73, 61, 72.5, 71.5, 59.5, 74.5, 63)
print(summary(weight)) 


###################################################
### code chunk number 9: f1
###################################################
oldpar<-par(mfrow=c(1,2),mar=c(4,4,0,1),cex=0.6)
h1<- graph.freq(weight,col=colors()[84],frequency=1,las=2,density=20,ylim=c(0,12),
                ylab="Frequency") 
x<-h1$breaks
h2<- plot(h1, frequency =2, axes= FALSE,ylim=c(0,0.4),xlab="weight",ylab="Relative (%)")
polygon.freq(h2, col=colors()[84], lwd=2, frequency =2)
axis(1,x,cex=0.6,las=2)
y<-seq(0,0.4,0.1)
axis(2, y,y*100,cex=0.6,las=1) 
par(oldpar)


###################################################
### code chunk number 10: tutorial.Rnw:144-145
###################################################
stat.freq(h1)


###################################################
### code chunk number 11: tutorial.Rnw:164-165
###################################################
print(summary(h1),row.names=FALSE)


###################################################
### code chunk number 12: tutorial.Rnw:170-174
###################################################
sturges.freq(weight)
inter.freq(h1)
join.freq(h1,1:3) -> h3
print(summary(h3))


###################################################
### code chunk number 13: f2
###################################################
oldpar<-par(mfrow=c(1,2),mar=c(4,4,0,1),cex=0.8)
plot(h3, frequency=2,col=colors()[84],ylim=c(0,0.6),axes=FALSE,xlab="weight",ylab="%",border=0)
y<-seq(0,0.6,0.2)
axis(2,y,y*100,las=2)
axis(1,h3$breaks)
normal.freq(h3,frequency=2,col=colors()[90])
ogive.freq(h3,col=colors()[84],xlab="weight")
par(oldpar)


###################################################
### code chunk number 14: f3
###################################################
oldpar<-par(mfrow=c(1,2),mar=c(4,3,2,1),cex=0.6)
h4<-hist(weight,xlab="Classes (h4)")
table.freq(h4)
# this is possible
# hh<-graph.freq(h4,plot=FALSE)
# summary(hh)
# new class
classes <- c(0, 10, 20, 30, 40, 50)
freq <- c(3, 8, 15, 18, 6)
h5 <- graph.freq(classes,counts=freq, xlab="Classes (h5)",main="Histogram grouped data")
par(oldpar)


###################################################
### code chunk number 15: tutorial.Rnw:222-223
###################################################
print(summary(h5),row.names=FALSE)


###################################################
### code chunk number 16: tutorial.Rnw:267-273
###################################################
str(design.crd)
trt <- c("A", "B", "C")
repeticion <- c(4, 3, 4)
outdesign <- design.crd(trt,r=repeticion,seed=777,serie=0)
book1 <- outdesign$book
head(book1)


###################################################
### code chunk number 17: tutorial.Rnw:282-290
###################################################
str(design.rcbd)
trt <- c("A", "B", "C","D","E")
repeticion <- 4
outdesign <- design.rcbd(trt,r=repeticion, seed=-513, serie=2)
# book2 <- outdesign$book
book2<- zigzag(outdesign) # zigzag numeration
print(outdesign$sketch)
print(matrix(book2[,1],byrow = TRUE, ncol = 5))


###################################################
### code chunk number 18: tutorial.Rnw:296-300
###################################################
str(design.lsd)
trt <- c("A", "B", "C", "D")
outdesign <- design.lsd(trt, seed=543, serie=2)
print(outdesign$sketch)


###################################################
### code chunk number 19: tutorial.Rnw:304-306
###################################################
book <- zigzag(outdesign)
print(matrix(book[,1],byrow = TRUE, ncol = 4))


###################################################
### code chunk number 20: tutorial.Rnw:314-319
###################################################
str(design.graeco)
trt1 <- c("A", "B", "C", "D")
trt2 <- 1:4
outdesign <- design.graeco(trt1,trt2, seed=543, serie=2)
print(outdesign$sketch)


###################################################
### code chunk number 21: tutorial.Rnw:324-326
###################################################
book <- zigzag(outdesign)
print(matrix(book[,1],byrow = TRUE, ncol = 4))


###################################################
### code chunk number 22: tutorial.Rnw:333-341
###################################################
str(design.youden)
varieties<-c("perricholi","yungay","maria bonita","tomasa")
r<-3
outdesign <-design.youden(varieties,r,serie=2,seed=23)
print(outdesign$sketch)
book <- outdesign$book
print(book) # field book.
print(matrix(as.numeric(book[,1]),byrow = TRUE, ncol = r))


###################################################
### code chunk number 23: tutorial.Rnw:346-348
###################################################
book <- zigzag(outdesign)
print(matrix(as.numeric(book[,1]),byrow = TRUE, ncol = r))


###################################################
### code chunk number 24: tutorial.Rnw:356-363
###################################################
str(design.bib)
trt <- c("A", "B", "C", "D", "E" )
k <- 4
outdesign <- design.bib(trt,k, seed=543, serie=2)
book5 <- outdesign$book
outdesign$statistics
outdesign$parameters


###################################################
### code chunk number 25: tutorial.Rnw:368-369
###################################################
outdesign$sketch


###################################################
### code chunk number 26: tutorial.Rnw:375-377
###################################################
book <- zigzag(outdesign)
matrix(book[,1],byrow = TRUE, ncol = 4)


###################################################
### code chunk number 27: tutorial.Rnw:384-390
###################################################
str(design.cyclic)
trt <- c("A", "B", "C", "D", "E", "F" )
outdesign <- design.cyclic(trt,k=3, r=6, seed=543, serie=2)
book6 <- outdesign$book
outdesign$sketch[[1]]
outdesign$sketch[[2]]


###################################################
### code chunk number 28: tutorial.Rnw:396-400
###################################################
book <- zigzag(outdesign)
array(book$plots,c(3,6,2))->X
t(X[,,1])
t(X[,,2])


###################################################
### code chunk number 29: tutorial.Rnw:406-407
###################################################
str(design.lattice)


###################################################
### code chunk number 30: tutorial.Rnw:413-420
###################################################
trt<-letters[1:9]
outdesign <-design.lattice(trt, r = 3, serie = 2, seed = 33, 
    kinds =  "Super-Duper")
book7 <- outdesign$book
outdesign$parameters
outdesign$sketch
head(book7)


###################################################
### code chunk number 31: tutorial.Rnw:425-430
###################################################
book <- zigzag(outdesign)
array(book$plots,c(3,3,3)) -> X
t(X[,,1])
t(X[,,2])
t(X[,,3])


###################################################
### code chunk number 32: tutorial.Rnw:436-446
###################################################
str(design.alpha)
trt <- letters[1:15]
outdesign <- design.alpha(trt,k=3,r=2,seed=543)
book8 <- outdesign$book
outdesign$statistics
outdesign$sketch
# codification of the plots
A<-array(book8[,1], c(3,5,2))
t(A[,,1])
t(A[,,2])


###################################################
### code chunk number 33: tutorial.Rnw:451-455
###################################################
book <- zigzag(outdesign)
A<-array(book[,1], c(3,5,2))
t(A[,,1])
t(A[,,2])


###################################################
### code chunk number 34: tutorial.Rnw:461-468
###################################################
str(design.dau)
rm(list=ls())
trt1 <- c("A", "B", "C", "D")
trt2 <- c("t","u","v","w","x","y","z")
outdesign <- design.dau(trt1, trt2, r=5, seed=543, serie=2)
book9 <- outdesign$book
with(book9,by(trt, block,as.character))


###################################################
### code chunk number 35: tutorial.Rnw:473-476
###################################################
book <- zigzag(outdesign)
with(book,by(plots, block, as.character))
head(book)


###################################################
### code chunk number 36: tutorial.Rnw:484-485
###################################################
str(design.split)


###################################################
### code chunk number 37: tutorial.Rnw:490-499
###################################################
trt1<-c("A","B","C","D")
trt2<-c("a","b","c")
outdesign <-design.split(trt1,trt2,r=3,serie=2,seed=543)
book10 <- outdesign$book
head(book10)
p<-book10$trt1[seq(1,36,3)]
q<-NULL
for(i in 1:12) 
q <- c(q,paste(book10$trt2[3*(i-1)+1],book10$trt2[3*(i-1)+2], book10$trt2[3*(i-1)+3]))


###################################################
### code chunk number 38: tutorial.Rnw:504-505
###################################################
print(t(matrix(p,c(4,3))))


###################################################
### code chunk number 39: tutorial.Rnw:510-511
###################################################
print(t(matrix(q,c(4,3))))


###################################################
### code chunk number 40: tutorial.Rnw:516-518
###################################################
book <- zigzag(outdesign)
head(book,5)


###################################################
### code chunk number 41: tutorial.Rnw:524-525
###################################################
str(design.strip)


###################################################
### code chunk number 42: tutorial.Rnw:530-542
###################################################
trt1<-c("A","B","C","D")
trt2<-c("a","b","c")
outdesign <-design.strip(trt1,trt2,r=3,serie=2,seed=543)
book11 <- outdesign$book
head(book11)
t3<-paste(book11$trt1, book11$trt2)
B1<-t(matrix(t3[1:12],c(4,3)))
B2<-t(matrix(t3[13:24],c(3,4)))
B3<-t(matrix(t3[25:36],c(3,4)))
print(B1)
print(B2)
print(B3)


###################################################
### code chunk number 43: tutorial.Rnw:547-553
###################################################
book <- zigzag(outdesign)
head(book)
array(book$plots,c(3,4,3))->X
t(X[,,1])
t(X[,,2])
t(X[,,3])


###################################################
### code chunk number 44: tutorial.Rnw:559-560
###################################################
str(design.ab)


###################################################
### code chunk number 45: tutorial.Rnw:563-564
###################################################
trt <- c (4,2,3) # three factors with  4,2 and 3 levels.


###################################################
### code chunk number 46: tutorial.Rnw:568-572
###################################################
trt<-c(3,2) # factorial 3x2
outdesign <-design.ab(trt, r=3, serie=2)
book12 <- outdesign$book
head(book12) # print of the field book


###################################################
### code chunk number 47: tutorial.Rnw:577-579
###################################################
book <- zigzag(outdesign)
head(book)


###################################################
### code chunk number 48: tutorial.Rnw:584-589
###################################################
trt<-c(2,2,2)
crd<-design.ab(trt, r=5, serie=2,design="crd")
names(crd)
crd$parameters
head(crd$book)


###################################################
### code chunk number 49: tutorial.Rnw:595-596
###################################################
str(design.mat)


###################################################
### code chunk number 50: tutorial.Rnw:602-606
###################################################
data(sweetpotato)
X <- design.mat(sweetpotato,1)
n<-nrow(X)
print(X)


###################################################
### code chunk number 51: tutorial.Rnw:609-619
###################################################
X<- rbind(X,c(0,3,3,3,3))
y<-sweetpotato[,2]
Y<- c(y,0)
XX<-t(X)%*%X
XY<-t(X)%*%Y
YY<-t(Y)%*%Y
beta<-solve(XX,XY)
e<-Y-X%*%beta
tau<- beta[-1,]
print(tau) #Efectos de tratamientos


###################################################
### code chunk number 52: tutorial.Rnw:622-628
###################################################
SCtot<-t(Y)%*%Y 
SCb<-t(beta)%*%XX%*%beta
SCe<-t(e)%*%e
SS<-rbind(SCtot,SCb,SCe)
dimnames(SS)<-list(c("Total","model","Error"),"SC")
print(SS)


###################################################
### code chunk number 53: tutorial.Rnw:631-640
###################################################
trt<-tau
t<-length(trt)
XX1<-XX[2:5,2:5]
SCtot<-t(Y)%*%Y-beta[1,1]%*%XX[1,1]%*%beta[1,1] ; GLtotal<- n-1
SCt<-trt%*%XX1%*%trt ; GLt <- t-1; GLe<-GLtotal - GLt
SC<-cbind(c(GLtotal,GLt,GLe),c(SCtot,SCt,SCe))
MS<-cbind(SC,SC[,2]/SC[,1])
dimnames(MS)<-list(c("Total","trt","Error"),c("Gl","SC","CM"))
print(MS)


###################################################
### code chunk number 54: tutorial.Rnw:651-655
###################################################
data(sweetpotato)
model<-aov(yield~virus, data=sweetpotato)
cv.model(model)
with(sweetpotato,mean(yield))


###################################################
### code chunk number 55: tutorial.Rnw:660-662
###################################################
df<-df.residual(model)
MSerror<-deviance(model)/df


###################################################
### code chunk number 56: tutorial.Rnw:669-671
###################################################
# comparison <- LSD.test(yield,virus,df,MSerror)
LSD.test(model, "virus",console=TRUE)


###################################################
### code chunk number 57: tutorial.Rnw:676-678
###################################################
# comparison <- LSD.test(yield, virus,df, MSerror, group=FALSE)
outLSD <-LSD.test(model, "virus", group=FALSE,console=TRUE)


###################################################
### code chunk number 58: tutorial.Rnw:685-687
###################################################
options(digits=2)
print(outLSD)


###################################################
### code chunk number 59: tutorial.Rnw:693-698
###################################################
LSD.test(model, "virus", group=FALSE, p.adj= "bon",console=TRUE)
out<-LSD.test(model, "virus", group=TRUE, p.adj= "holm")
print(out$group)
out<-LSD.test(model, "virus", group=FALSE, p.adj= "holm")
print(out$comparison)


###################################################
### code chunk number 60: tutorial.Rnw:709-710
###################################################
duncan.test(model, "virus",console=TRUE)


###################################################
### code chunk number 61: tutorial.Rnw:717-719
###################################################
# SNK.test(model, "virus", alpha=0.05,console=TRUE)
SNK.test(model, "virus", group=FALSE,console=TRUE)


###################################################
### code chunk number 62: tutorial.Rnw:726-728
###################################################
# REGW.test(model, "virus", alpha=0.05,console=TRUE)
REGW.test(model, "virus", group=FALSE,console=TRUE)


###################################################
### code chunk number 63: tutorial.Rnw:735-737
###################################################
outHSD<- HSD.test(model, "virus",console=TRUE)
outHSD


###################################################
### code chunk number 64: tutorial.Rnw:744-750
###################################################
A<-sweetpotato[-c(4,5,7),]
modelUnbalanced <- aov(yield ~ virus, data=A)
outUn <-HSD.test(modelUnbalanced, "virus",group=FALSE, unbalanced = TRUE)
print(outUn$comparison)
outUn <-HSD.test(modelUnbalanced, "virus",group=TRUE, unbalanced = TRUE)
print(outUn$groups)


###################################################
### code chunk number 65: tutorial.Rnw:757-761
###################################################
outUn <-HSD.test(modelUnbalanced, "virus",group=FALSE)
print(outUn$comparison[,1:2])
outUn <-HSD.test(modelUnbalanced, "virus",group=TRUE)
print(outUn$groups)


###################################################
### code chunk number 66: tutorial.Rnw:768-771
###################################################
# variance analysis:
anova(model)
with(sweetpotato,waller.test(yield,virus,df,MSerror,Fc= 17.345, group=FALSE,console=TRUE))


###################################################
### code chunk number 67: tutorial.Rnw:775-776
###################################################
outWaller <- waller.test(model, "virus", group=FALSE,console=FALSE)


###################################################
### code chunk number 68: tutorial.Rnw:781-783
###################################################
names(outWaller)
print(outWaller$comparison)


###################################################
### code chunk number 69: tutorial.Rnw:788-789
###################################################
outWaller$statistics


###################################################
### code chunk number 70: tutorial.Rnw:796-799
###################################################
# analysis of variance: 
scheffe.test(model,"virus", group=TRUE,console=TRUE,
main="Yield of sweetpotato\nDealt with different virus")


###################################################
### code chunk number 71: tutorial.Rnw:805-806
###################################################
outScheffe <- scheffe.test(model,"virus", group=FALSE, console=TRUE)


###################################################
### code chunk number 72: tutorial.Rnw:813-815
###################################################
# modelABC <-aov (y ~ A * B * C, data)
# compare <-LSD.test (modelABC, c ("A", "B", "C"),console=TRUE)


###################################################
### code chunk number 73: tutorial.Rnw:822-834
###################################################
yield <-scan (text =
 "6 7 9 13 16 20 8 8 9
  7 8 8 12 17 18 10 9 12
  9 9 9 14 18 21 11 12 11
  8 10 10 15 16 22 9 9 9 "
 )
block <-gl (4, 9)
clone <-rep (gl (3, 3, labels = c ("c1", "c2", "c3")), 4)
nitrogen <-rep (gl (3, 1, labels = c ("n1", "n2", "n3")), 12)
A <-data.frame (block, clone, nitrogen, yield)
head (A)
outAOV <-aov (yield ~ block + clone * nitrogen, data = A)


###################################################
### code chunk number 74: tutorial.Rnw:837-840
###################################################
anova (outAOV)
outFactorial <-LSD.test (outAOV, c("clone", "nitrogen"), 
main = "Yield ~ block + nitrogen + clone + clone:nitrogen",console=TRUE)


###################################################
### code chunk number 75: tutorial.Rnw:843-849 (eval = FALSE)
###################################################
## oldpar<-par(mar=c(3,3,2,0))
## pic1<-bar.err(outFactorial$means,variation="range",ylim=c(5,25), bar=FALSE,col=0,las=1)
## points(pic1$index,pic1$means,pch=18,cex=1.5,col="blue")
## axis(1,pic1$index,labels=FALSE)
## title(main="average and range\nclon:nitrogen")
## par(oldpar)


###################################################
### code chunk number 76: tutorial.Rnw:856-865
###################################################
# Example linear estimation and design of experiments. (Joshi)
# Institute of Social Sciences Agra, India
# 6 varieties of wheat in 10 blocks of 3 plots each.
block<-gl(10,3)
variety<-c(1,2,3,1,2,4,1,3,5,1,4,6,1,5,6,2,3,6,2,4,5,2,5,6,3,4,5,3, 4,6)
Y<-c(69,54,50,77,65,38,72,45,54,63,60,39,70,65,54,65,68,67,57,60,62, 
59,65,63,75,62,61,59,55,56)
head(cbind(block,variety,Y))
BIB.test(block, variety, Y,console=TRUE)


###################################################
### code chunk number 77: tutorial.Rnw:870-873
###################################################
out <-BIB.test(block, trt=variety, Y, test="tukey", group=FALSE, console=TRUE)
names(out)
rm(block,variety)


###################################################
### code chunk number 78: tutorial.Rnw:886-891
###################################################
# alpha design 
Genotype<-paste("geno",1:30,sep="")
r<-2
k<-3
plan<-design.alpha(Genotype,k,r,seed=5)


###################################################
### code chunk number 79: tutorial.Rnw:897-899
###################################################
yield <-c(5,2,7,6,4,9,7,6,7,9,6,2,1,1,3,2,4,6,7,9,8,7,6,4,3,2,2,1,1, 
          2,1,1,2,4,5,6,7,8,6,5,4,3,1,1,2,5,4,2,7,6,6,5,6,4,5,7,6,5,5,4)


###################################################
### code chunk number 80: tutorial.Rnw:903-907
###################################################
data<-data.frame(plan$book,yield)
# The analysis: 
modelPBIB <- with(data,PBIB.test(block, Genotype, replication, yield, k=3, 
                                 group=TRUE, console=TRUE))


###################################################
### code chunk number 81: tutorial.Rnw:942-948
###################################################
trt<-c(1,8,5,5,2,9,2,7,3,6,4,9,4,6,9,8,7,6,1,5,8,3,2,7,3,7,2,1,3,4,6,4,9,5,8,1)
yield<-c(48.76,10.83,12.54,11.07,22,47.43,27.67,30,13.78,37,42.37,39,14.46,30.69,42.01,
22,42.8,28.28,50,24,24,15.42,30,23.8,19.68,31,23,41,12.9,49.95,25,45.57,30,20,18,43.81)
sqr<-rep(gl(4,3),3)
block<-rep(1:12,3)
modelLattice<-PBIB.test(block,trt,sqr,yield,k=3,console=TRUE, method="VC")


###################################################
### code chunk number 82: tutorial.Rnw:963-968
###################################################
block<-c(rep("I",7),rep("II",6),rep("III",7))
trt<-c("A","B","C","D","g","k","l","A","B","C","D","e","i","A","B", "C",
"D","f","h","j")
yield<-c(83,77,78,78,70,75,74,79,81,81,91,79,78,92,79,87,81,89,96, 82)
head(data.frame(block, trt, yield))


###################################################
### code chunk number 83: tutorial.Rnw:973-974
###################################################
by(trt,block,as.character)


###################################################
### code chunk number 84: tutorial.Rnw:979-980
###################################################
by(yield,block,as.character)


###################################################
### code chunk number 85: tutorial.Rnw:985-988
###################################################
modelDAU<- DAU.test(block,trt,yield,method="lsd",console=TRUE)
options(digits = 2)
modelDAU$means


###################################################
### code chunk number 86: tutorial.Rnw:991-993
###################################################
modelDAU<- DAU.test(block,trt,yield,method="lsd",group=FALSE,console=FALSE)
head(modelDAU$comparison,8)


###################################################
### code chunk number 87: tutorial.Rnw:1015-1017
###################################################
data(corn)
str(corn)


###################################################
### code chunk number 88: tutorial.Rnw:1024-1025
###################################################
str(kruskal)


###################################################
### code chunk number 89: tutorial.Rnw:1030-1031
###################################################
outKruskal<-with(corn,kruskal(observation,method,group=TRUE, main="corn", console=TRUE))


###################################################
### code chunk number 90: tutorial.Rnw:1040-1044
###################################################
out<-with(corn,kruskal(observation,method,group=TRUE, main="corn", p.adj="holm"))
print(out$group)
out<-with(corn,kruskal(observation,method,group=FALSE, main="corn", p.adj="holm"))
print(out$comparison)


###################################################
### code chunk number 91: tutorial.Rnw:1050-1051
###################################################
str(friedman)


###################################################
### code chunk number 92: tutorial.Rnw:1056-1059
###################################################
data(grass)
out<-with(grass,friedman(judge,trt, evaluation,alpha=0.05, group=FALSE,
main="Data of the book of Conover",console=TRUE))


###################################################
### code chunk number 93: tutorial.Rnw:1065-1066
###################################################
str(waerden.test)


###################################################
### code chunk number 94: tutorial.Rnw:1071-1073
###################################################
data(sweetpotato)
outWaerden<-with(sweetpotato,waerden.test(yield,virus,alpha=0.01,group=TRUE,console=TRUE))


###################################################
### code chunk number 95: tutorial.Rnw:1077-1080
###################################################

names(outWaerden)



###################################################
### code chunk number 96: tutorial.Rnw:1083-1084
###################################################
out<-with(sweetpotato,waerden.test(yield,virus,group=FALSE,console=TRUE))


###################################################
### code chunk number 97: tutorial.Rnw:1092-1093
###################################################
str(Median.test)


###################################################
### code chunk number 98: tutorial.Rnw:1096-1097
###################################################
str(Median.test)


###################################################
### code chunk number 99: tutorial.Rnw:1102-1107
###################################################
data(sweetpotato)
outMedian<-with(sweetpotato,Median.test(yield,virus,console=TRUE))
names(outMedian)
outMedian$statistics
outMedian$medians


###################################################
### code chunk number 100: f14
###################################################
oldpar<-par(mfrow=c(2,2),mar=c(3,3,1,1),cex=0.8)
# Graphics
bar.group(outMedian$groups,ylim=c(0,50))
bar.group(outMedian$groups,xlim=c(0,50),horiz = TRUE)
plot(outMedian)
plot(outMedian,variation="IQR",horiz = TRUE)
par(oldpar)


###################################################
### code chunk number 101: tutorial.Rnw:1131-1132
###################################################
str(durbin.test)


###################################################
### code chunk number 102: tutorial.Rnw:1137-1147
###################################################
days <-gl(7,3)
chemical<-c("A","B","D","A","C","E","C","D","G","A","F","G", "B","C","F", 
"B","E","G","D","E","F")
toxic<-c(0.465,0.343,0.396,0.602,0.873,0.634,0.875,0.325,0.330, 0.423,0.987,0.426,
0.652,1.142,0.989,0.536,0.409,0.309, 0.609,0.417,0.931)
head(data.frame(days,chemical,toxic))
out<-durbin.test(days,chemical,toxic,group=FALSE,console=TRUE,
main="Logarithm of the toxic dose")
names(out)
out$statistics


###################################################
### code chunk number 103: tutorial.Rnw:1160-1164
###################################################
# model <-aov (yield ~ fertilizer, data = field) 
# out <-LSD.test (model, "fertilizer", group = TRUE) 
# bar.group (out $ group)
str(bar.group)


###################################################
### code chunk number 104: tutorial.Rnw:1172-1176
###################################################
# model <-aov (yield ~ fertilizer, data = field) 
# out <-LSD.test (model, "fertilizer", group = TRUE) 
# bar.err(out$means)
str(bar.err)


###################################################
### code chunk number 105: f4
###################################################
oldpar<-par(mfrow=c(2,2),mar=c(3,3,2,1),cex=0.7)
c1<-colors()[480]; c2=colors()[65] 
bar.err(outHSD$means, variation="range",ylim=c(0,50),col=c1,las=1)
bar.err(outHSD$means, variation="IQR",horiz=TRUE, xlim=c(0,50),col=c2,las=1)
plot(outHSD, variation="range",las=1)
plot(outHSD, horiz=TRUE, variation="SD",las=1)
par(oldpar)


###################################################
### code chunk number 106: tutorial.Rnw:1204-1219 (eval = FALSE)
###################################################
## oldpar<-par(mfrow=c(2,2),cex=0.7,mar=c(3.5,1.5,3,1))
## C1<-bar.err(modelPBIB$means[1:7, ], ylim=c(0,9), col=0, main="C1", 
## variation="range",border=3,las=2)
## C2<-bar.err(modelPBIB$means[8:15,], ylim=c(0,9), col=0, main="C2", 
## variation="range", border =4,las=2)
## # Others graphic
## C3<-bar.err(modelPBIB$means[16:22,], ylim=c(0,9), col=0, main="C3", 
## variation="range",border =2,las=2)
## C4<-bar.err(modelPBIB$means[23:30,], ylim=c(0,9), col=0, main="C4", 
## variation="range", border =6,las=2)
## # Lattice graphics
## par(oldpar)
## oldpar<-par(mar=c(2.5,2.5,1,0),cex=0.6)
## bar.group(modelLattice$group,ylim=c(0,55),density=10,las=1)
## par(oldpar)


###################################################
### code chunk number 107: f13
###################################################
# model : yield ~ virus
# Important group=TRUE
oldpar<-par(mfrow=c(1,2),mar=c(3,3,1,1),cex=0.8)
x<-duncan.test(model, "virus", group=TRUE)
plot(x,las=1)
plot(x,variation="IQR",horiz=TRUE,las=1)
par(oldpar)


###################################################
### code chunk number 108: f5
###################################################
# function (x, main = NULL, color1 = "red", color2 = "blue", 
#    color3 = "black", cex.axis = 0.8, las = 1, pch = 20, 
#    bty = "l", cex = 0.8, lwd = 1, xlab = "", ylab = "", 
#   ...)
# model : yield ~ virus
# Important group=FALSE
x<-HSD.test(model, "virus", group=FALSE)
diffograph(x,cex.axis=0.9,xlab="Yield",ylab="Yield",cex=0.9)


###################################################
### code chunk number 109: tutorial.Rnw:1272-1276
###################################################
options(digit=2)
f <- system.file("external/dataStb.csv", package="agricolae")
dataStb<-read.csv(f)
stability.par(dataStb, rep=4, MSerror=1.8, alpha=0.1, main="Genotype",console=TRUE)


###################################################
### code chunk number 110: tutorial.Rnw:1283-1286
###################################################
output <- stability.par(dataStb, rep=4, MSerror=2)
names(output)
print(output$stability)


###################################################
### code chunk number 111: tutorial.Rnw:1291-1295
###################################################
data5<-dataStb[,1:5]
altitude<-c(1200, 1300, 800, 1600, 2400)
stability <- stability.par(data5,rep=4,MSerror=2, cova=TRUE, name.cov= "altitude",
file.cov=altitude)


###################################################
### code chunk number 112: tutorial.Rnw:1302-1306
###################################################
data <- data.frame(name=row.names(dataStb), dataStb)
output<-stability.nonpar(data, "YIELD", ranking=TRUE)
names(output)
output$statistics


###################################################
### code chunk number 113: tutorial.Rnw:1323-1324
###################################################
  str(AMMI)


###################################################
### code chunk number 114: tutorial.Rnw:1327-1328
###################################################
  str(plot.AMMI)


###################################################
### code chunk number 115: tutorial.Rnw:1331-1339
###################################################
data(plrv)
model<-with(plrv,AMMI(Locality, Genotype, Rep, Yield, console=FALSE))
names(model)
model$ANOVA
model$analysis
pc <- model$analysis[, 1]
pc12<-sum(pc[1:2])
pc123<-sum(pc[1:3])


###################################################
### code chunk number 116: f6
###################################################
oldpar<-par(cex=0.4,mar=c(4,4,1,2))
plot(model,las=1,xlim=c(-5,6))
par(oldpar)


###################################################
### code chunk number 117: tutorial.Rnw:1357-1364
###################################################
data(plrv)
model<- with(plrv,AMMI(Locality, Genotype, Rep, Yield, console=FALSE))
index<-index.AMMI(model)
# Crops with improved stability according AMMI.
print(index[order(index[,3]),])
# Crops with better response and improved stability according AMMI.
print(index[order(index[,4]),])


###################################################
### code chunk number 118: f7
###################################################
oldpar<-par(cex=0.6,mar=c(3,3,2,1))
data(pamCIP)
rownames(pamCIP)<-substr(rownames(pamCIP),1,6)
output<-consensus(pamCIP,distance="binary", method="complete", nboot=5)
par(oldpar)


###################################################
### code chunk number 119: f8
###################################################
oldpar<-par(cex=0.6,mar=c(3,3,1.5,1))
out1<- hcut(output,h=0.4,group=8,type="t",edgePar = list(lty=1:2, col=colors()[c(42,84)]),
main="group 8" ,col.text="blue",cex.text=1,las=1)
par(oldpar)


###################################################
### code chunk number 120: tutorial.Rnw:1402-1403
###################################################
names(output)


###################################################
### code chunk number 121: tutorial.Rnw:1409-1412
###################################################
dend <- as.dendrogram(output$dendrogram)
data <- output$table.dend
head(output$table.dend)


###################################################
### code chunk number 122: tutorial.Rnw:1414-1418 (eval = FALSE)
###################################################
## oldpar<-par(mar=c(3,3,1,1),cex=0.6)
## plot(dend,type="r",edgePar = list(lty=1:2, col=colors()[c(42,84)]) ,las=1)
## text(data[,3],data[,4],data[,5],col="blue",cex=1)
## par(oldpar)


###################################################
### code chunk number 123: tutorial.Rnw:1425-1429
###################################################
data(soil)
# set.seed(9473)
simulated <- montecarlo(soil$pH,1000)
h<-graph.freq(simulated,nclass=7,plot=FALSE)


###################################################
### code chunk number 124: f9
###################################################
oldpar<-par(mar=c(2,0,2,1),cex=0.6)
plot(density(soil$pH),axes=FALSE,main="pH density of the soil\ncon Ralstonia",xlab="",lwd=4)
lines(density(simulated), col="blue", lty=4,lwd=4)
axis(1,0:12)
legend("topright",c("Original","Simulated"),lty=c(1,4),col=c("black", "blue"), lwd=4)
par(oldpar)


###################################################
### code chunk number 125: tutorial.Rnw:1445-1446
###################################################
round(table.freq(h),2)


###################################################
### code chunk number 126: tutorial.Rnw:1449-1450
###################################################
summary(soil$pH)


###################################################
### code chunk number 127: tutorial.Rnw:1453-1454
###################################################
summary(simulated)


###################################################
### code chunk number 128: tutorial.Rnw:1460-1467
###################################################
data(potato)
potato[,1]<-as.factor(potato[,1])
potato[,2]<-as.factor(potato[,2])
model<-"cutting~variety + date + variety:date"
analysis<-resampling.model(model, potato, k=100)
Xsol<-as.matrix(round(analysis$solution,2))
print(Xsol,na.print = "")


###################################################
### code chunk number 129: tutorial.Rnw:1476-1477
###################################################
simModel <- simulation.model(model, potato, k=100,console=TRUE)


###################################################
### code chunk number 130: tutorial.Rnw:1479-1480
###################################################
ab<-simModel$simulation[3,3]


###################################################
### code chunk number 131: tutorial.Rnw:1487-1493
###################################################
corr.x<- matrix(c(1,0.5,0.5,1),c(2,2))
corr.y<- rbind(0.6,0.7)
names<-c("X1","X2")
dimnames(corr.x)<-list(names,names)
dimnames(corr.y)<-list(names,"Y")
output<-path.analysis(corr.x,corr.y)


###################################################
### code chunk number 132: tutorial.Rnw:1496-1497
###################################################
output


###################################################
### code chunk number 133: tutorial.Rnw:1532-1540
###################################################
rm(list=ls())
options(digits = 2)
data(heterosis)
str(heterosis)
site2<-subset(heterosis,heterosis[,1]==2)
site2<-subset(site2[,c(2,5,6,8)],site2[,4]!="Control")
output1<-with(site2,lineXtester(Replication, Female, Male, v2))
options(digits = 7)


###################################################
### code chunk number 134: f10
###################################################
oldpar<-par(mar=c(3,3,4,1),cex=0.7)
data(rice)
table<-index.smith(rice, col="blue",
 main="Interaction between the CV and the plot size",type="l",xlab="Size")
par(oldpar)


###################################################
### code chunk number 135: tutorial.Rnw:1560-1562
###################################################
uniformity <- data.frame(table$uniformity)
head(uniformity)


###################################################
### code chunk number 136: tutorial.Rnw:1571-1574
###################################################
data(paracsho)
species <- paracsho[79:87,4:6]
species


###################################################
### code chunk number 137: tutorial.Rnw:1579-1580
###################################################
output <- index.bio(species[,3],method="Shannon",level=95,nboot=200)


###################################################
### code chunk number 138: tutorial.Rnw:1588-1591
###################################################
data(soil)
correlation(soil[,2:4],method="pearson")
with(soil,correlation(pH,soil[,3:4],method="pearson"))


###################################################
### code chunk number 139: tutorial.Rnw:1602-1604
###################################################
data(RioChillon)
with(RioChillon$babies,tapply.stat(yield,farmer,function(x) max(x)-min(x)))


###################################################
### code chunk number 140: tutorial.Rnw:1623-1626
###################################################
data(sweetpotato)
model <- model<-aov(yield ~ virus, data=sweetpotato)
cv.model(model)


###################################################
### code chunk number 141: tutorial.Rnw:1634-1635
###################################################
x<-c(3,4,5,2,3,4,5,6,4,NA,7)


###################################################
### code chunk number 142: tutorial.Rnw:1640-1641
###################################################
skewness(x)


###################################################
### code chunk number 143: tutorial.Rnw:1646-1647
###################################################
kurtosis(x)


###################################################
### code chunk number 144: tutorial.Rnw:1654-1663
###################################################
q<-5
f<-15
K<-seq(10,1000,100)
n<-length(K)
y<-rep(0,3*n)
dim(y)<-c(n,3)
for(i in 1:n) y[i,1]<-waller(K[i],q,f,Fc=2)
for(i in 1:n) y[i,2]<-waller(K[i],q,f,Fc=4)
for(i in 1:n) y[i,3]<-waller(K[i],q,f,Fc=8)


###################################################
### code chunk number 145: tutorial.Rnw:1669-1677 (eval = FALSE)
###################################################
## oldpar<-par(mar=c(3,3,4,1),cex=0.7)
## plot(K,y[,1],type="l",col="blue",ylab="waller",bty="l")
## lines(K,y[,2],type="l",col="brown",lty=2,lwd=2)
## lines(K,y[,3],type="l",col="green",lty=4,lwd=2)
## legend("topleft",c("2","4","8"),col=c("blue","brown","green"),lty=c(1,8,20),
## lwd=2,title="Fc")
## title(main="Waller in function of K")
## par(oldpar)


###################################################
### code chunk number 146: tutorial.Rnw:1682-1698
###################################################
K<-100
Fc<-1.2
q<-c(seq(6,20,1),30,40,100)
f<-c(seq(4,20,2),24,30)
n<-length(q)
m<-length(f)
W.D <-rep(0,n*m)
dim(W.D)<-c(n,m)
for (i in 1:n) {
for (j in 1:m) {
W.D[i,j]<-waller(K, q[i], f[j], Fc)
}}
W.D<-round(W.D,2)
dimnames(W.D)<-list(q,f)
cat("table: Waller Duncan k=100, F=1.2")
print(W.D)


###################################################
### code chunk number 147: tutorial.Rnw:1705-1710
###################################################
days<-c(7,14,21,28,35,42)
evaluation<-data.frame(E1=10,E2=40,E3=50,E4=70,E5=80,E6=90)
print(evaluation)
absolute1 <-audpc(evaluation,days)
relative1 <-round(audpc(evaluation,days,"relative"),2)


###################################################
### code chunk number 148: tutorial.Rnw:1717-1719
###################################################
absolute2 <-audps(evaluation,days)
relative2 <-round(audps(evaluation,days,"relative"),2)


###################################################
### code chunk number 149: f11
###################################################
oldpar<-par(mfrow=c(1,2),mar=c(3,3,1,1),cex=0.7)
plot(days, evaluation,type="h",ylim=c(0,100),axes=FALSE,col= colors()[42],xlab="Days", ylab="Evaluation")
lines(days,evaluation,col= colors()[42])
axis(1,days)
axis(2,seq(0,100,20),las=2)
rect(7,0,42,100)
text(15,80,substitute(paste("Audpc Abs.=",A1),list(A1=absolute1)))
text(15,70,substitute(paste("Audpc Rel.=",A2),list(A2=relative1)))
x<-seq(3.5,45.5,7)
evaluation<-c(evaluation,evaluation[6])
plot(x, evaluation,type="s",ylim=c(0,100),axes=FALSE,col= colors()[42],xlab="Days", ylab="Evaluation")
points(x,evaluation,type="h",col= colors()[42])
points(days,evaluation[-6],pch=16,col=4)
text(days,5,days,col=4,cex=1)
axis(1,x,pos=0)
axis(2,seq(0,100,20),las=2)
rect(3.5,0,45.5,100)
text(13,80,substitute(paste("Audps Abs.=",A1),list(A1=absolute2)))
text(13,70,substitute(paste("Audps Rel.=",A2),list(A2=relative2)))


###################################################
### code chunk number 150: tutorial.Rnw:1759-1765
###################################################
data(potato)
potato[,1]<-as.factor(potato[,1])
model<-lm(cutting ~ date + variety,potato)
df<-df.residual(model)
MSerror<-deviance(model)/df
analysis<-with(potato,nonadditivity(cutting, date, variety, df, MSerror))


###################################################
### code chunk number 151: tutorial.Rnw:1774-1801
###################################################
options(digits=2)
f <- system.file("external/weather.csv", package="agricolae")
weather <- read.csv(f,header=FALSE)
f <- system.file("external/severity.csv", package="agricolae")
severity <- read.csv(f)
weather[,1]<-as.Date(weather[,1],format = "%m/%d/%Y")
# Parameters dates
dates<-c("2000-03-25","2000-04-09","2000-04-12","2000-04-16","2000-04-22")
dates<-as.Date(dates)
EmergDate <- as.Date("2000/01/19")
EndEpidDate <- as.Date("2000-04-22")
dates<-as.Date(dates)
NoReadingsH<- 1
RHthreshold <- 90
WS<-weatherSeverity(weather,severity,dates,EmergDate,EndEpidDate,
NoReadingsH,RHthreshold)
# Parameters to Lateblight function
InocDate<-"2000-03-18"
LGR <- 0.00410
IniSpor <- 0
SR <- 292000000
IE <- 1.0
LP <- 2.82
InMicCol <- 9
Cultivar <- "NICOLA"
ApplSys <- "NOFUNGICIDE"
main<-"Cultivar: NICOLA"


###################################################
### code chunk number 152: f12
###################################################
oldpar<-par(mar=c(3,3,4,1),cex=0.7)
#--------------------------
model<-lateblight(WS, Cultivar,ApplSys, InocDate, LGR,IniSpor,SR,IE, 
LP,MatTime='LATESEASON',InMicCol,main=main,type="l",xlim=c(65,95),lwd=1.5,
xlab="Time (days after emergence)", ylab="Severity (Percentage)")
par(oldpar)


###################################################
### code chunk number 153: tutorial.Rnw:1821-1824
###################################################
head(model$Gfile)
str(model$Ofile)
head(model$Ofile[,1:7])


###################################################
### code chunk number 154: tutorial.Rnw:1828-1834
###################################################
x<- model$Ofile$nday
y<- model$Ofile$SimSeverity
w<- model$Gfile$nday
z<- model$Gfile$MeanSeverity
Min<-model$Gfile$MinObs
Max<-model$Gfile$MaxObs


###################################################
### code chunk number 155: tutorial.Rnw:1837-1845 (eval = FALSE)
###################################################
## oldpar<-par(mar=c(3,2.5,1,1),cex=0.7)
## plot(x,y,type="l",xlim=c(65,95),lwd=1.5,xlab="Time (days after emergence)",
## ylab="Severity (Percentage)")
## points(w,z,col="red",cex=1,pch=19); npoints <- length(w)
## for ( i in 1:npoints)segments(w[i],Min[i],w[i],Max[i],lwd=1.5,col="red")
## legend("topleft",c("Disease progress curves","Weather-Severity"),
## title="Description",lty=1,pch=c(3,19),col=c("black","red"))
## par(oldpar)