| Title: | Statistical Procedures for Agricultural Research |
|---|---|
| Description: | Original idea was presented in the thesis "A statistical analysis tool for agricultural research" to obtain the degree of Master on science, National Engineering University (UNI), Lima-Peru. Some experimental data for the examples come from the CIP and others research. Agricolae offers extensive functionality on experimental design especially for agricultural and plant breeding experiments, which can also be useful for other purposes. It supports planning of lattice, Alpha, Cyclic, Complete Block, Latin Square, Graeco-Latin Squares, augmented block, factorial, split and strip plot designs. There are also various analysis facilities for experimental data, e.g. treatment comparison procedures and several non-parametric tests comparison, biodiversity indexes and consensus cluster. |
| Authors: | Felipe de Mendiburu |
| Maintainer: | Felipe de Mendiburu <[email protected]> |
| License: | GPL |
| Version: | 1.3-7 |
| Built: | 2024-12-15 07:33:18 UTC |
| Source: | CRAN |
This package contains functionality for the Statistical Analysis of experimental designs applied specially for field experiments in agriculture and plant breeding.
| Package: | agricolae |
| Type: | Package |
| Version: | 1.3-7 |
| Date: | 2023-10-22 |
| License: | GPL |
Planning of field experiments: lattice, factorial, RCBD, CRD, Latin Square, Youden, Graeco, BIB, Alpha design, Cyclic, augmented block, split and strip plot Designs. Comparison of multi-location trials: AMMI, Index AMMI Stability (biplot, triplot), comparison between treatments: LSD, Bonferroni and other p-adjust, HSD, Waller, Student Newman Keuls SNK, Duncan, REGW, Scheffe; Non parametric tests: Kruskal, Friedman, Durbin, Van Der Waerden, Median. Analysis of genetic experiments: North Carolina designs, LinexTester, Balanced Incomplete Block, Strip plot, Split-Plot, Partially Balanced Incomplete Block, analysis Mother and baby trials (see data RioChillon). Resampling and simulation: resampling.model, simulation.model, montecarlo, lateblight Simulator for potato. Ecology: Biodiversity Index, Path Analysis. Soil Uniformity: Smith's Index. Cluster Analysis: Consensus Cluster. Descriptive statistics utilities: *.freq
Felipe de Mendiburu Statistical Engineer Master in Systems Engineering Professor of Applied Statistics
Maintainer: Felipe de Mendiburu <[email protected]>
De Mendiburu, Felipe (2009). Una herramienta de analisis estadistico para la investigacion agricola. Tesis. Universidad Nacional de Ingenieria (UNI-PERU).
Universidad Nacional Agraria La Molina, Lima-PERU. Facultad de Economia y Planificacion Departamento Academico de Estadistica e Informatica
Additive Main Effects and Multiplicative Interaction Models (AMMI) are widely used to analyze main effects and genotype by environment (GEN, ENV) interactions in multilocation variety trials. Furthermore, this function generates data to biplot, triplot graphs and analysis.
AMMI(ENV, GEN, REP, Y, MSE = 0,console=FALSE,PC=FALSE)AMMI(ENV, GEN, REP, Y, MSE = 0,console=FALSE,PC=FALSE)
ENV |
Environment |
GEN |
Genotype |
REP |
Replication |
Y |
Response |
MSE |
Mean Square Error |
console |
ouput TRUE or FALSE |
PC |
Principal components ouput TRUE or FALSE |
additional graphics see help(plot.AMMI).
ANOVA |
analysis of variance general |
genXenv |
class by, genopyte and environment |
analysis |
analysis of variance principal components |
means |
average genotype and environment |
biplot |
data to produce graphics |
PC |
class princomp |
F. de Mendiburu
Crossa, J. 1990. Statistical analysis of multilocation trials. Advances in Agronomy 44:55-85
# Full replications library(agricolae) # Example 1 data(plrv) model<- with(plrv,AMMI(Locality, Genotype, Rep, Yield, console=FALSE)) model$ANOVA # see help(plot.AMMI) # biplot plot(model) # biplot PC1 vs Yield plot(model, first=0,second=1, number=TRUE) # Example 2 data(CIC) data1<-CIC$comas[,c(1,6,7,17,18)] data2<-CIC$oxapampa[,c(1,6,7,19,20)] cic <- rbind(data1,data2) model<-with(cic,AMMI(Locality, Genotype, Rep, relative)) model$ANOVA plot(model,0,1,angle=20,ecol="brown") # Example 3 # Only means. Mean square error is well-known. data(sinRepAmmi) REP <- 3 MSerror <- 93.24224 #startgraph model<-with(sinRepAmmi,AMMI(ENV, GEN, REP, YLD, MSerror,PC=TRUE)) # print anova print(model$ANOVA,na.print = "") # Biplot with the one restored observed. plot(model,0,1) # with principal components model$PC is class "princomp" pc<- model$PC pc$loadings summary(pc) biplot(pc) # Principal components by means of the covariance similar AMMI # It is to compare results with AMMI cova<-cov(model$genXenv) values<-eigen(cova) total<-sum(values$values) round(values$values*100/total,2) # AMMI: 64.81 18.58 13.50 3.11 0.00# Full replications library(agricolae) # Example 1 data(plrv) model<- with(plrv,AMMI(Locality, Genotype, Rep, Yield, console=FALSE)) model$ANOVA # see help(plot.AMMI) # biplot plot(model) # biplot PC1 vs Yield plot(model, first=0,second=1, number=TRUE) # Example 2 data(CIC) data1<-CIC$comas[,c(1,6,7,17,18)] data2<-CIC$oxapampa[,c(1,6,7,19,20)] cic <- rbind(data1,data2) model<-with(cic,AMMI(Locality, Genotype, Rep, relative)) model$ANOVA plot(model,0,1,angle=20,ecol="brown") # Example 3 # Only means. Mean square error is well-known. data(sinRepAmmi) REP <- 3 MSerror <- 93.24224 #startgraph model<-with(sinRepAmmi,AMMI(ENV, GEN, REP, YLD, MSerror,PC=TRUE)) # print anova print(model$ANOVA,na.print = "") # Biplot with the one restored observed. plot(model,0,1) # with principal components model$PC is class "princomp" pc<- model$PC pc$loadings summary(pc) biplot(pc) # Principal components by means of the covariance similar AMMI # It is to compare results with AMMI cova<-cov(model$genXenv) values<-eigen(cova) total<-sum(values$values) round(values$values*100/total,2) # AMMI: 64.81 18.58 13.50 3.11 0.00
Draws a polygon or a circumference around the center of the Biplot with a proportional radio at the longest distance of the genotype.
AMMI.contour(model, distance, shape, ...)AMMI.contour(model, distance, shape, ...)
model |
Object |
distance |
Circumference radius >0 and <=1 |
shape |
Numerical, relating to the shape of the polygon outline. |
... |
Parameters corresponding to the R lines function |
First, it is necessary to execute the AMMI function. It is only valid for the BIPLOT function but not for the TRIPLOT one.
Genotypes within and outside the area.
distance |
Distance from genotype to origin (0,0) |
Complement graphics AMMI
Felipe de Mendiburu
library(agricolae) # see AMMI. data(sinRepAmmi) Environment <- sinRepAmmi$ENV Genotype <- sinRepAmmi$GEN Yield <- sinRepAmmi$YLD REP <- 3 MSerror <- 93.24224 model<-AMMI(Environment, Genotype, REP, Yield, MSerror) plot(model) AMMI.contour(model,distance=0.7,shape=8,col="red",lwd=2,lty=5)library(agricolae) # see AMMI. data(sinRepAmmi) Environment <- sinRepAmmi$ENV Genotype <- sinRepAmmi$GEN Yield <- sinRepAmmi$YLD REP <- 3 MSerror <- 93.24224 model<-AMMI(Environment, Genotype, REP, Yield, MSerror) plot(model) AMMI.contour(model,distance=0.7,shape=8,col="red",lwd=2,lty=5)
Area Under Disease Progress Curve. The AUDPC measures the disease throughout a period. The AUDPC is the area that is determined by the sum of trapezes under the curve.
audpc(evaluation, dates, type = "absolute")audpc(evaluation, dates, type = "absolute")
evaluation |
Table of data of the evaluations: Data frame |
dates |
Vector of dates corresponding to each evaluation |
type |
relative, absolute |
AUDPC. For the illustration one considers three evaluations (14, 21 and 28 days) and percentage of damage in the plant 40, 80 and 90 (interval between dates of evaluation 7 days). AUDPC = 1045. The evaluations can be at different interval.
Vector with relative or absolute audpc.
Felipe de Mendiburu
Campbell, C. L., L. V. Madden. (1990): Introduction to Plant Disease Epidemiology. John Wiley & Sons, New York City.
library(agricolae) dates<-c(14,21,28) # days # example 1: evaluation - vector evaluation<-c(40,80,90) audpc(evaluation,dates) # example 2: evaluation: dataframe nrow=1 evaluation<-data.frame(E1=40,E2=80,E3=90) # percentages plot(dates,evaluation,type="h",ylim=c(0,100),col="red",axes=FALSE) title(cex.main=0.8,main="Absolute or Relative AUDPC\nTotal area = 100*(28-14)=1400") lines(dates,evaluation,col="red") text(dates,evaluation+5,evaluation) text(18,20,"A = (21-14)*(80+40)/2") text(25,60,"B = (28-21)*(90+80)/2") text(25,40,"audpc = A+B = 1015") text(24.5,33,"relative = audpc/area = 0.725") abline(h=0) axis(1,dates) axis(2,seq(0,100,5),las=2) lines(rbind(c(14,40),c(14,100)),lty=8,col="green") lines(rbind(c(14,100),c(28,100)),lty=8,col="green") lines(rbind(c(28,90),c(28,100)),lty=8,col="green") # It calculates audpc absolute absolute<-audpc(evaluation,dates,type="absolute") print(absolute) rm(evaluation, dates, absolute) # example 3: evaluation dataframe nrow>1 data(disease) dates<-c(1,2,3) # week evaluation<-disease[,c(4,5,6)] # It calculates audpc relative index <-audpc(evaluation, dates, type = "relative") # Correlation between the yield and audpc correlation(disease$yield, index, method="kendall") # example 4: days infile data(CIC) comas <- CIC$comas oxapampa <- CIC$oxapampa dcomas <- names(comas)[9:16] days<- as.numeric(substr(dcomas,2,3)) AUDPC<- audpc(comas[,9:16],days) relative<-audpc(comas[,9:16],days,type = "relative") h1<-graph.freq(AUDPC,border="red",density=4,col="blue") table.freq(h1) h2<-graph.freq(relative,border="red",density=4,col="blue", frequency=2, ylab="relative frequency")library(agricolae) dates<-c(14,21,28) # days # example 1: evaluation - vector evaluation<-c(40,80,90) audpc(evaluation,dates) # example 2: evaluation: dataframe nrow=1 evaluation<-data.frame(E1=40,E2=80,E3=90) # percentages plot(dates,evaluation,type="h",ylim=c(0,100),col="red",axes=FALSE) title(cex.main=0.8,main="Absolute or Relative AUDPC\nTotal area = 100*(28-14)=1400") lines(dates,evaluation,col="red") text(dates,evaluation+5,evaluation) text(18,20,"A = (21-14)*(80+40)/2") text(25,60,"B = (28-21)*(90+80)/2") text(25,40,"audpc = A+B = 1015") text(24.5,33,"relative = audpc/area = 0.725") abline(h=0) axis(1,dates) axis(2,seq(0,100,5),las=2) lines(rbind(c(14,40),c(14,100)),lty=8,col="green") lines(rbind(c(14,100),c(28,100)),lty=8,col="green") lines(rbind(c(28,90),c(28,100)),lty=8,col="green") # It calculates audpc absolute absolute<-audpc(evaluation,dates,type="absolute") print(absolute) rm(evaluation, dates, absolute) # example 3: evaluation dataframe nrow>1 data(disease) dates<-c(1,2,3) # week evaluation<-disease[,c(4,5,6)] # It calculates audpc relative index <-audpc(evaluation, dates, type = "relative") # Correlation between the yield and audpc correlation(disease$yield, index, method="kendall") # example 4: days infile data(CIC) comas <- CIC$comas oxapampa <- CIC$oxapampa dcomas <- names(comas)[9:16] days<- as.numeric(substr(dcomas,2,3)) AUDPC<- audpc(comas[,9:16],days) relative<-audpc(comas[,9:16],days,type = "relative") h1<-graph.freq(AUDPC,border="red",density=4,col="blue") table.freq(h1) h2<-graph.freq(relative,border="red",density=4,col="blue", frequency=2, ylab="relative frequency")
A better estimate of disease progress is the area under the disease progress stairs (AUDPS). The AUDPS approach improves the estimation of disease progress by giving a weight closer to optimal to the first and last observations.
audps(evaluation, dates, type = "absolute")audps(evaluation, dates, type = "absolute")
evaluation |
Table of data of the evaluations: Data frame |
dates |
Vector of dates corresponding to each evaluation |
type |
relative, absolute |
AUDPS. For the illustration one considers three evaluations (14, 21 and 28 days) and percentage of damage in the plant 40, 80 and 90 (interval between dates of evaluation 7 days). AUDPS = 1470. The evaluations can be at different interval. AUDPS= sum( rectangle area by interval in times evaluation ) see example.
Vector with relative or absolute audps.
Felipe de Mendiburu
Ivan Simko, and Hans-Peter Piepho, (2012). The area under the disease progress stairs: Calculation, advantage, and application. Phytopathology 102:381- 389.
library(agricolae) dates<-c(14,21,28) # days # example 1: evaluation - vector evaluation<-c(40,80,90) audps(evaluation,dates) audps(evaluation,dates,"relative") x<-seq(10.5,31.5,7) y<-c(40,80,90,90) plot(x,y,"s",ylim=c(0,100),xlim=c(10,32),axes=FALSE,col="red" ,ylab="",xlab="") title(cex.main=0.8,main="Absolute or Relative AUDPS\nTotal area=(31.5-10.5)*100=2100", ylab="evaluation",xlab="dates" ) points(x,y,type="h") z<-c(14,21,28) points(z,y[-3],col="blue",lty=2,pch=19) points(z,y[-3],col="blue",lty=2,pch=19) axis(1,x,pos=0) axis(2,c(0,40,80,90,100),las=2) text(dates,evaluation+5,dates,col="blue") text(14,20,"A = (17.5-10.5)*40",cex=0.8) text(21,40,"B = (24.5-17.5)*80",cex=0.8) text(28,60,"C = (31.5-24.5)*90",cex=0.8) text(14,95,"audps = A+B+C = 1470") text(14,90,"relative = audps/area = 0.7") # It calculates audpc absolute absolute<-audps(evaluation,dates,type="absolute") print(absolute) rm(evaluation, dates, absolute)library(agricolae) dates<-c(14,21,28) # days # example 1: evaluation - vector evaluation<-c(40,80,90) audps(evaluation,dates) audps(evaluation,dates,"relative") x<-seq(10.5,31.5,7) y<-c(40,80,90,90) plot(x,y,"s",ylim=c(0,100),xlim=c(10,32),axes=FALSE,col="red" ,ylab="",xlab="") title(cex.main=0.8,main="Absolute or Relative AUDPS\nTotal area=(31.5-10.5)*100=2100", ylab="evaluation",xlab="dates" ) points(x,y,type="h") z<-c(14,21,28) points(z,y[-3],col="blue",lty=2,pch=19) points(z,y[-3],col="blue",lty=2,pch=19) axis(1,x,pos=0) axis(2,c(0,40,80,90,100),las=2) text(dates,evaluation+5,dates,col="blue") text(14,20,"A = (17.5-10.5)*40",cex=0.8) text(21,40,"B = (24.5-17.5)*80",cex=0.8) text(28,60,"C = (31.5-24.5)*90",cex=0.8) text(14,95,"audps = A+B+C = 1470") text(14,90,"relative = audps/area = 0.7") # It calculates audpc absolute absolute<-audps(evaluation,dates,type="absolute") print(absolute) rm(evaluation, dates, absolute)
It plots bars of the averages of treatments and standard error or standard deviance. It uses the objects generated by a procedure of comparison like LSD, HSD, Kruskal and Waller-Duncan.
bar.err(x,variation=c("SE","SD","range","IQR"),horiz=FALSE, bar=TRUE,...)bar.err(x,variation=c("SE","SD","range","IQR"),horiz=FALSE, bar=TRUE,...)
x |
object means of the comparisons the LSD.test, HSD.test,...,etc |
variation |
SE=standard error, range=Max-Min or IQR=interquartil range |
horiz |
Horizontal or vertical bars |
bar |
paint bar |
... |
Parameters of the function barplot() |
x: data frame formed by 5 columns: name of the bars, height, level out: LSD.test, HSD, waller.test, scheffe.test, duncan.test, SNK.test, friedman, kruskal, waerden.test and Median.test.
A list with numeric vectors giving the coordinates of all the bar midpoints drawn.
x |
eje-1 coordinate |
height |
eje-2 coordinate by group |
Felipe de Mendiburu
LSD.test, HSD.test,
waller.test, kruskal, bar.group
library(agricolae) data(sweetpotato) model<-aov(yield~virus,data=sweetpotato) out <- waller.test(model,"virus", console=TRUE, main="Yield of sweetpotato\ndealt with different virus") oldpar<-par(mfrow=c(2,2),cex=1) bar.err(out$means,variation="range",horiz=TRUE,xlim=c(0,45),angle=125,density=6, main="range") bar.err(out$means,variation="SD",ylim=c(0,45),col=colors()[30], main="Standard deviation",density=8) bar.err(out$means,variation="SE",horiz=TRUE,xlim=c(0,45),density=8, col="brown",main="Standard error") bar.err(out$means,variation="range",ylim=c(0,45),bar=FALSE,col="green", main="range") par(oldpar) oldpar<-par(mfrow=c(1,2),cex=1) bar.err(out$means,variation="range",ylim=c(0,45),bar=FALSE,col=0) abline(h=0) # horiz = TRUE bar.err(out$means,variation="SE",horiz=TRUE,xlim=c(0,45),bar=FALSE,col=0) #startgraph par(oldpar) #endgraphlibrary(agricolae) data(sweetpotato) model<-aov(yield~virus,data=sweetpotato) out <- waller.test(model,"virus", console=TRUE, main="Yield of sweetpotato\ndealt with different virus") oldpar<-par(mfrow=c(2,2),cex=1) bar.err(out$means,variation="range",horiz=TRUE,xlim=c(0,45),angle=125,density=6, main="range") bar.err(out$means,variation="SD",ylim=c(0,45),col=colors()[30], main="Standard deviation",density=8) bar.err(out$means,variation="SE",horiz=TRUE,xlim=c(0,45),density=8, col="brown",main="Standard error") bar.err(out$means,variation="range",ylim=c(0,45),bar=FALSE,col="green", main="range") par(oldpar) oldpar<-par(mfrow=c(1,2),cex=1) bar.err(out$means,variation="range",ylim=c(0,45),bar=FALSE,col=0) abline(h=0) # horiz = TRUE bar.err(out$means,variation="SE",horiz=TRUE,xlim=c(0,45),bar=FALSE,col=0) #startgraph par(oldpar) #endgraph
It plots bars of the averages of treatments to compare. It uses the objects generated by a procedure of comparison like LSD, HSD, Kruskall, Waller-Duncan, Friedman or Durbin. It can also display the 'average' value over each bar in a bar chart.
bar.group(x,horiz=FALSE, decreasing=TRUE, ...)bar.group(x,horiz=FALSE, decreasing=TRUE, ...)
x |
Object created by a test of comparison |
horiz |
Horizontal or vertical bars |
decreasing |
Logical, decreasing order of the mean |
... |
Parameters of the function barplot() |
x: data frame formed by 5 columns: name of the bars, height and level of the bar.
A list with numeric vectors giving the coordinates of all the bar midpoints drawn.
x |
eje-1 coordinate |
height |
eje-2 coordinate by group |
Felipe de Meniburu
LSD.test, HSD.test, kruskal
, friedman, durbin.test, waller.test
, plot.group
# Example 1 library(agricolae) data(sweetpotato) model<-aov(yield~virus,data=sweetpotato) comparison<- LSD.test(model,"virus",alpha=0.01,group=TRUE) print(comparison$groups) oldpar<-par(cex=1.5) bar.group(comparison$groups,horiz=TRUE,density=8,col="blue",border="red", xlim=c(0,50),las=1) title(cex.main=0.8,main="Comparison between\ntreatment means",xlab="Yield",ylab="Virus") # Example 2 library(agricolae) x <- 1:4 y <- c(0.29, 0.44, 0.09, 0.49) xy <- data.frame(x,y,y) par(oldpar) oldpar<-par(cex=1.5) bar.group(xy,density=30,angle=90,col="brown",border=FALSE,ylim=c(0,0.6),lwd=2,las=1) par(oldpar)# Example 1 library(agricolae) data(sweetpotato) model<-aov(yield~virus,data=sweetpotato) comparison<- LSD.test(model,"virus",alpha=0.01,group=TRUE) print(comparison$groups) oldpar<-par(cex=1.5) bar.group(comparison$groups,horiz=TRUE,density=8,col="blue",border="red", xlim=c(0,50),las=1) title(cex.main=0.8,main="Comparison between\ntreatment means",xlab="Yield",ylab="Virus") # Example 2 library(agricolae) x <- 1:4 y <- c(0.29, 0.44, 0.09, 0.49) xy <- data.frame(x,y,y) par(oldpar) oldpar<-par(cex=1.5) bar.group(xy,density=30,angle=90,col="brown",border=FALSE,ylim=c(0,0.6),lwd=2,las=1) par(oldpar)
Analysis of variance BIB and comparison mean adjusted.
BIB.test(block, trt, y, test = c("lsd","tukey","duncan","waller","snk"), alpha = 0.05, group = TRUE,console=FALSE)BIB.test(block, trt, y, test = c("lsd","tukey","duncan","waller","snk"), alpha = 0.05, group = TRUE,console=FALSE)
block |
blocks |
trt |
Treatment |
y |
Response |
test |
Comparison treatments |
alpha |
Significant test |
group |
logical |
console |
logical, print output |
Test of comparison treatment. lsd: Least significant difference. tukey: Honestly significant differente. duncan: Duncan's new multiple range test waller: Waller-Duncan test. snk: Student-Newman-Keuls (SNK)
parameters |
Design parameters |
statistics |
Statistics of the model |
comparison |
Comparison between treatments |
means |
Adjusted mean and statistics summary |
groups |
Grouping of treatments |
F. de Mendiburu
Design of Experiments. Robert O. Kuehl. 2nd ed., Duxbury, 2000 Linear Estimation and Design of Experiments. D.D. Joshi. WILEY EASTERN LIMITED 1987, New Delhi, India. Introduction to experimental statistics. Ching Chun Li McGraw - Hill Book Company, Inc. New York. 1964
DAU.test, duncan.test, durbin.test,
friedman, HSD.test, kruskal,
LSD.test, Median.test, PBIB.test,
REGW.test, scheffe.test, SNK.test,
waerden.test, waller.test, plot.group
library(agricolae) # Example Design of Experiments. Robert O. Kuehl. 2da. Edicion. 2001 run<-gl(10,3) psi<-c(250,325,475,250,475,550,325,400,550,400,475,550,325,475,550, 250,400,475,250,325,400,250,400,550,250,325,550,325,400,475) monovinyl<-c(16,18,32,19,46,45,26,39,61,21,35,55,19,47,48,20,33,31,13,13,34,21, 30,52,24,10,50,24,31,37) out<-BIB.test(run,psi,monovinyl,test="waller",group=FALSE) print(out) bar.err(out$means,variation="range",ylim=c(0,60),bar=FALSE,col=0) out<-BIB.test(run,psi,monovinyl,test="waller",group=TRUE) out<-BIB.test(run,psi,monovinyl,test="tukey",group=TRUE,console=TRUE) out<-BIB.test(run,psi,monovinyl,test="tukey",group=FALSE,console=TRUE) rm(run,psi,monovinyl,out) # Example linear estimation and design of experiments. D.D. Joshi. 1987 # Professor of Statistics, Institute of Social Sciences Agra, India # 6 varieties of wheat crop in a BIB whit 10 blocks of 3 plots each. y <-c(69,77,72,63,70,54,65,65,57,59,50,45,68,75,59,38,60,60,62, 55,54,65,62,65,61,39,54,67,63,56) varieties<-gl(6,5) block <- c(1,2,3,4,5,1,2,6,7,8,1,3,6,9,10,2,4,7,9,10,3,5,7,8,9,4,5,6,8,10) BIB.test(block, varieties, y) # Example Introduction to experimental statistics. Ching Chun Li. 1964 # pag. 395 table. 27.2 # 7 trt, k=3 and b=7. y <-c(10,15,11,4,12,15,5,14,10,14,19,19,8,10,17,6,11,12,5,14,21) block<-gl(7,3) trt <- c(1,2,4,2,3,5,3,4,6,4,5,7,1,5,6,2,6,7,1,3,7) out<-BIB.test(block, trt, y, test="duncan") bar.group(out$groups,col="blue",density=4,ylim=c(0,max(y))) rm(y,block,trt,out)library(agricolae) # Example Design of Experiments. Robert O. Kuehl. 2da. Edicion. 2001 run<-gl(10,3) psi<-c(250,325,475,250,475,550,325,400,550,400,475,550,325,475,550, 250,400,475,250,325,400,250,400,550,250,325,550,325,400,475) monovinyl<-c(16,18,32,19,46,45,26,39,61,21,35,55,19,47,48,20,33,31,13,13,34,21, 30,52,24,10,50,24,31,37) out<-BIB.test(run,psi,monovinyl,test="waller",group=FALSE) print(out) bar.err(out$means,variation="range",ylim=c(0,60),bar=FALSE,col=0) out<-BIB.test(run,psi,monovinyl,test="waller",group=TRUE) out<-BIB.test(run,psi,monovinyl,test="tukey",group=TRUE,console=TRUE) out<-BIB.test(run,psi,monovinyl,test="tukey",group=FALSE,console=TRUE) rm(run,psi,monovinyl,out) # Example linear estimation and design of experiments. D.D. Joshi. 1987 # Professor of Statistics, Institute of Social Sciences Agra, India # 6 varieties of wheat crop in a BIB whit 10 blocks of 3 plots each. y <-c(69,77,72,63,70,54,65,65,57,59,50,45,68,75,59,38,60,60,62, 55,54,65,62,65,61,39,54,67,63,56) varieties<-gl(6,5) block <- c(1,2,3,4,5,1,2,6,7,8,1,3,6,9,10,2,4,7,9,10,3,5,7,8,9,4,5,6,8,10) BIB.test(block, varieties, y) # Example Introduction to experimental statistics. Ching Chun Li. 1964 # pag. 395 table. 27.2 # 7 trt, k=3 and b=7. y <-c(10,15,11,4,12,15,5,14,10,14,19,19,8,10,17,6,11,12,5,14,21) block<-gl(7,3) trt <- c(1,2,4,2,3,5,3,4,6,4,5,7,1,5,6,2,6,7,1,3,7) out<-BIB.test(block, trt, y, test="duncan") bar.group(out$groups,col="blue",density=4,ylim=c(0,max(y))) rm(y,block,trt,out)
Statistic analysis of the Carolina I, II and III genetic designs.
carolina(model,data)carolina(model,data)
model |
Constant |
data |
Data frame |
model = 1,2 and 3 is I, II and III see carolina1,2 and 3.
model |
model analysis (I, II or III) of caroline design |
and variance and additive variance of male, female and male.female interaction.
Felipe de Mendiburu
Biometrical Methods in Quantitative Genetic Analysis, Singh, Chaudhary. 1979
library(agricolae) data(DC) carolina1 <- DC$carolina1 # str(carolina1) output<-carolina(model=1,carolina1) output[][-1] carolina2 <- DC$carolina2 # str(carolina2) majes<-subset(carolina2,carolina2[,1]==1) majes<-majes[,c(2,5,4,3,6:8)] output<-carolina(model=2,majes[,c(1:4,6)]) output[][-1] carolina3 <- DC$carolina3 # str(carolina3) output<-carolina(model=3,carolina3) output[][-1]library(agricolae) data(DC) carolina1 <- DC$carolina1 # str(carolina1) output<-carolina(model=1,carolina1) output[][-1] carolina2 <- DC$carolina2 # str(carolina2) majes<-subset(carolina2,carolina2[,1]==1) majes<-majes[,c(2,5,4,3,6:8)] output<-carolina(model=2,majes[,c(1:4,6)]) output[][-1] carolina3 <- DC$carolina3 # str(carolina3) output<-carolina(model=3,carolina3) output[][-1]
Incidents and performance of healthy tubers and rotten potato field infested with naturally Ralstonia solanacearum Race 3/Bv 2A, after application of inorganic amendments and a rotation crop in Carhuaz Peru, 2006.
data(Chz2006)data(Chz2006)
The format is: List of 2
amendmenta factor
cropa factor
blocka numeric vector, replications
planta numeric vector, number plant
wilt_percenta numeric vector, wilt percentage at 60 days
healtha numeric vector, kg/8m2
rota numeric vector, kg/8m2
Application of inorganic amendment and crop rotation to control bacterial wilt of the potato (MBP).
Experimental field, 2006. Data Kindly provided by Pedro Aley.
International Potato Center. CIP - Lima Peru.
library(agricolae) data(Chz2006) str(Chz2006) wilt<-Chz2006$wilt yield<-Chz2006$yield means <- tapply.stat(wilt[,5],wilt[,1:3],function(x) mean(x,na.rm=TRUE)) names(means)[4]<-"wilt_percent" model <- aov(wilt_percent ~ block + crop, means) anova(model) cv.model(model) yield<-yield[order(paste(yield[,1],yield[,2],yield[,3])),] correlation(means[,4],yield[,4],method="spearman")library(agricolae) data(Chz2006) str(Chz2006) wilt<-Chz2006$wilt yield<-Chz2006$yield means <- tapply.stat(wilt[,5],wilt[,1:3],function(x) mean(x,na.rm=TRUE)) names(means)[4]<-"wilt_percent" model <- aov(wilt_percent ~ block + crop, means) anova(model) cv.model(model) yield<-yield[order(paste(yield[,1],yield[,2],yield[,3])),] correlation(means[,4],yield[,4],method="spearman")
A study of Phytophthora infestans in the potato plant in the localities of Comas and Oxapampa in Peru, 2005.
data(CIC)data(CIC)
The format is: List of 2 (comas, oxapampa)
Localitya factor with levels Comas Oxapampa
Genotypea factor
Repa numeric vector, replications
E9a numeric vector, infestans percentaje to 9 days
AUDPCa numeric vector: the area under the disease-progress curve
Relativea numeric vector, relative area
comas: temperature=59.9 Fahrenheit, relative humidity=83.3 oxapampa: temperature=64.8 Fahrenheit, relative humidity=86.2 AUDPC and relative see function audpc(). help(audpc) Exx: Evaluation in percentaje, xx is days. ORD1, ORD2, SBLK and row are references location of the plot in the field.
Experimental field, 2004-2005. Data Kindly provided by Matilde Orrillo.
International Potato Center. CIP - Lima Peru.
library(agricolae) data(CIC) CIC$comas CIC$oxapampalibrary(agricolae) data(CIC) CIC$comas CIC$oxapampa
An evaluation over a time period.
data(clay)data(clay)
A data frame with 69 observations on the following 3 variables.
per.claya numeric vector
daysa numeric vector
ralstoniaa numeric vector
Experimental field.
International Potato Center. CIP - Lima Peru.
library(agricolae) data(clay) str(clay)library(agricolae) data(clay) str(clay)
Fifty-three potato varieties developed by the breeding program of the International Potato Center and released in different countries around the world were evaluated for their resistance to late blight in two locations in Peru.
data(ComasOxapampa)data(ComasOxapampa)
A data frame with 168 observations on the following 4 variables.
cultivara factor with 56 levels
replicationa factor with 3 levels
comasa numeric vector
oxapampaa numeric vector
The experimental design was a randomized complete block design with 3 replications of 15 apical stem cuttings in Oxapampa and 10 tubers in Mariscal Castilla. Plots were 11.9 x 18.5 m in size with 30 cm in-row and 0.9 m between-row spacings. Spreader rows around plots were used at each site. Mancozeb was applied weekly until 30 days after transplanting or planting, after which the plants were left to natural infection. Due to climatic conditions not conductive to the disease in Oxapampa, inoculum was enhanced with local isolate (POX 067, with virulence R1, 2, 3, 4, 5, 6, 7, 10, 11) at a concentration of 5000-sporangia/ ml at 49 days after planting. Percentage of foliar infection was estimated visually every 3 days for 8 times in Oxapampa and every 7 days for 12 times in Comas, then values were converted to the relative area under the diseases progress curve (rAUPDC). rAUDPC rankings were analyzed for phenotypic stability with nonparametric measures.
Experimental field, 2002. Data Kindly provided by Wilmer Perez.
International Potato Center. CIP - Lima Peru.
library(agricolae) data(ComasOxapampa) # Oxapampa (10 35 31 S latitude, 75 23 0 E longitude, 1813 m.a.s.l ) # Comas, Mariscal Castilla (11 42 54 S latitude, 75 04 45 E longitude, 2800 m.a.s.l,) # cultivars LBr-40 (resistant), Cruza 148 (moderately resistant) and Pimpernell (susceptible) str(ComasOxapampa) means <- tapply.stat(ComasOxapampa[,3:4],ComasOxapampa$cultivar,mean) correlation(means$comas,means$oxapampa, method="kendall")library(agricolae) data(ComasOxapampa) # Oxapampa (10 35 31 S latitude, 75 23 0 E longitude, 1813 m.a.s.l ) # Comas, Mariscal Castilla (11 42 54 S latitude, 75 04 45 E longitude, 2800 m.a.s.l,) # cultivars LBr-40 (resistant), Cruza 148 (moderately resistant) and Pimpernell (susceptible) str(ComasOxapampa) means <- tapply.stat(ComasOxapampa[,3:4],ComasOxapampa$cultivar,mean) correlation(means$comas,means$oxapampa, method="kendall")
The criterion of the consensus is to produce many trees by means of boostrap and to such calculate the relative frequency with members of the clusters.
consensus(data,distance=c("binary","euclidean","maximum","manhattan", "canberra", "minkowski", "gower","chisq"),method=c("complete","ward","single","average", "mcquitty","median", "centroid"),nboot=500,duplicate=TRUE,cex.text=1, col.text="red", ...)consensus(data,distance=c("binary","euclidean","maximum","manhattan", "canberra", "minkowski", "gower","chisq"),method=c("complete","ward","single","average", "mcquitty","median", "centroid"),nboot=500,duplicate=TRUE,cex.text=1, col.text="red", ...)
data |
data frame |
distance |
method distance, see dist() |
method |
method cluster, see hclust() |
nboot |
The number of bootstrap samples desired. |
duplicate |
control is TRUE other case is FALSE |
cex.text |
size text on percentage consensus |
col.text |
color text on percentage consensus |
... |
parameters of the plot dendrogram |
distance: "euclidean", "maximum", "manhattan", "canberra", "binary", "minkowski", "gower", "chisq". Method: "ward", "single", "complete", "average", "mcquitty", "median", "centroid". see functions: dist(), hclust() and daisy() of cluster.
table.dend |
The groups and consensus percentage |
dendrogram |
The class object is hclust, dendrogram plot |
duplicate |
Homonymous elements |
F. de Mendiburu
An Introduction to the Boostrap. Bradley Efron and Robert J. Tibshirani. 1993. Chapman and Hall/CRC
library(agricolae) data(pamCIP) # only code rownames(pamCIP)<-substr(rownames(pamCIP),1,6) output<-consensus( pamCIP,distance="binary", method="complete",nboot=5) # Order consensus Groups<-output$table.dend[,c(6,5)] Groups<-Groups[order(Groups[,2],decreasing=TRUE),] print(Groups) ## Identification of the codes with the numbers. cbind(output$dendrogram$labels) ## To reproduce dendrogram dend<-output$dendrogram data<-output$table.dend plot(dend) text(data[,3],data[,4],data[,5]) # Other examples # classical dendrogram dend<-as.dendrogram(output$dendrogram) plot(dend,type="r",edgePar = list(lty=1:2, col=2:1)) text(data[,3],data[,4],data[,5],col="blue",cex=1) plot(dend,type="t",edgePar = list(lty=1:2, col=2:1)) text(data[,3],data[,4],data[,5],col="blue",cex=1) ## Without the control of duplicates output<-consensus( pamCIP,duplicate=FALSE,nboot=5) ## using distance gower, require cluster package. # output<-consensus( pamCIP,distance="gower", method="complete",nboot=5)library(agricolae) data(pamCIP) # only code rownames(pamCIP)<-substr(rownames(pamCIP),1,6) output<-consensus( pamCIP,distance="binary", method="complete",nboot=5) # Order consensus Groups<-output$table.dend[,c(6,5)] Groups<-Groups[order(Groups[,2],decreasing=TRUE),] print(Groups) ## Identification of the codes with the numbers. cbind(output$dendrogram$labels) ## To reproduce dendrogram dend<-output$dendrogram data<-output$table.dend plot(dend) text(data[,3],data[,4],data[,5]) # Other examples # classical dendrogram dend<-as.dendrogram(output$dendrogram) plot(dend,type="r",edgePar = list(lty=1:2, col=2:1)) text(data[,3],data[,4],data[,5],col="blue",cex=1) plot(dend,type="t",edgePar = list(lty=1:2, col=2:1)) text(data[,3],data[,4],data[,5],col="blue",cex=1) ## Without the control of duplicates output<-consensus( pamCIP,duplicate=FALSE,nboot=5) ## using distance gower, require cluster package. # output<-consensus( pamCIP,distance="gower", method="complete",nboot=5)
Data from a completely randomized design where four different methods of growing corn resulted in various yields per acre on various plots of ground where the four methods were tried. Ordinarily, only one statistical analysis is used, but here we will use the kuskal-wallis test so that a rough comparison may be made with the mediasn test.
data(corn)data(corn)
A data frame with 34 observations on the following 3 variables.
methoda numeric vector
observationa numeric vector
rxa numeric vector
The observations are ranked from the smallest, 77, of rank 1 to the largest 101, of rank N=34. Ties values receive the averarge rank.
Book: Practical Nonparametric Statistics.
Practical Nonparametrics Statistics. W.J. Conover. Third Edition, 1999.
data(corn) str(corn)data(corn) str(corn)
An exact correlation for ties or without ties. Methods of Kendall, Spearman and Pearson.
correl(x, y, method = "pearson",alternative="two.sided")correl(x, y, method = "pearson",alternative="two.sided")
x |
Vector |
y |
Vector |
method |
"pearson", "kendall", "spearman" |
alternative |
"two.sided", "less", "greater" |
The correlation of x,y vector with the statistical value and its probability
Felipe de Mendiburu
Numerical Recipes in C. Second Edition.
library(agricolae) data(soil) with(soil,correl(pH,clay,method="kendall")) with(soil,correl(pH,clay,method="spearman")) with(soil,correl(pH,clay,method="pearson"))library(agricolae) data(soil) with(soil,correl(pH,clay,method="kendall")) with(soil,correl(pH,clay,method="spearman")) with(soil,correl(pH,clay,method="pearson"))
It obtains the coefficients of correlation and p-value between all the variables of a data table. The methods to apply are Pearson, Spearman , Kendall and lin's concordance index. In case of not specifying the method, the Pearson method will be used. The results are similar to SAS.
correlation(x,y=NULL, method = c("pearson", "kendall", "spearman", "lin") ,alternative="two.sided")correlation(x,y=NULL, method = c("pearson", "kendall", "spearman", "lin") ,alternative="two.sided")
x |
table, matrix or vector |
y |
table, matrix or vector |
method |
"pearson", "kendall", "spearman", "lin" |
alternative |
"two.sided", "less", "greater" |
Parameters equal to function cor()
The correlation matrix with its probability
Felipe de Mendiburu
Lin LI. A concordance correlation coefficient to evaluate reproducibility. Biometrics. 1989; 45, 255-268.
library(agricolae) data(soil) # example 1 analysis<-correlation(soil[,2:8],method="pearson") analysis # Example 2: correlation between pH, variable 2 and other elements from soil. analysis<-with(soil,correlation(pH,soil[,3:8],method="pearson",alternative="less")) analysis # Example 3: correlation between pH and clay method kendall. with(soil,correlation(pH,clay,method="kendall", alternative="two.sided"))library(agricolae) data(soil) # example 1 analysis<-correlation(soil[,2:8],method="pearson") analysis # Example 2: correlation between pH, variable 2 and other elements from soil. analysis<-with(soil,correlation(pH,soil[,3:8],method="pearson",alternative="less")) analysis # Example 3: correlation between pH and clay method kendall. with(soil,correlation(pH,clay,method="kendall", alternative="two.sided"))
Data of cotton collected in experiments of two localities in Lima and Pisco, Peru.
data(cotton)data(cotton)
A data frame with 96 observations on the following 5 variables.
sitea factor with levels Lima Pisco
blocka factor with levels I II III IV V VI
lineagea numeric vector
epocaa numeric vector
yielda numeric vector
Book spanish: Metodos estadisticos para la investigacion. Autor: Calzada Benza Universidad Nacional Agraria - La Molina - Peru..
Book spanish: Metodos estadisticos para la investigacion. Autor: Calzada Benza Universidad Nacional Agraria - La Molina - Peru.
library(agricolae) data(cotton) str(cotton)library(agricolae) data(cotton) str(cotton)
It obtains the coefficient of variation of the experiment obtained by models lm() or aov()
cv.model(x)cv.model(x)
x |
object of model lm() or AOV() |
sqrt(MSerror)*100/mean(x)
Returns the coefficient of variation of the experiment according to the applied statistical model
Felipe de Mendiburu
LSD.test, HSD.test,
waller.test
# see examples from LSD , Waller-Duncan or HSD and complete with it: library(agricolae) # not run # cv<-cv.model(model)# see examples from LSD , Waller-Duncan or HSD and complete with it: library(agricolae) # not run # cv<-cv.model(model)
This process consists of finding the coefficient of the distances of similarity of binary tables (1 and 0) as used for scoring molecular marker data for presence and absence of PCR amplification products.
cv.similarity(A)cv.similarity(A)
A |
matrix of binary data |
Returns the coefficient of variation of the similarity model
Felipe de Mendiburu
# molecular markers. library(agricolae) data(markers) cv<-cv.similarity(markers)# molecular markers. library(agricolae) data(markers) cv<-cv.similarity(markers)
Analysis of variance Augmented block and comparison mean adjusted.
DAU.test(block, trt, y, method = c("lsd","tukey"),alpha=0.05,group=TRUE,console=FALSE)DAU.test(block, trt, y, method = c("lsd","tukey"),alpha=0.05,group=TRUE,console=FALSE)
block |
blocks |
trt |
Treatment |
y |
Response |
method |
Comparison treatments |
alpha |
Significant test |
group |
TRUE or FALSE |
console |
logical, print output |
Method of comparison treatment. lsd: Least significant difference. tukey: Honestly significant differente. The controls can have different repetitions, at least two, do not use missing data.
means |
Statistical summary of the study variable |
parameters |
Design parameters |
statistics |
Statistics of the model |
comparison |
Comparison between treatments |
groups |
Formation of treatment groups |
SE.difference |
Standard error of: |
vartau |
Variance-covariance matrix of the difference in treatments |
F. de Mendiburu
Federer, W. T. (1956). Augmented (or hoonuiaku) designs. Hawaiian Planters, Record LV(2):191-208.
BIB.test, duncan.test, durbin.test,
friedman, HSD.test, kruskal,
LSD.test, Median.test, PBIB.test,
REGW.test, scheffe.test, SNK.test,
waerden.test, waller.test, plot.group
library(agricolae) 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) out<- DAU.test(block,trt,yield,method="lsd", group=TRUE) print(out$groups) plot(out)library(agricolae) 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) out<- DAU.test(block,trt,yield,method="lsd", group=TRUE) print(out$groups) plot(out)
Data for the analysis of carolina I, II and III genetic design
data(DC)data(DC)
DC is list, 3 data.frame: carolina1(72 obs, 6 var), carolina2(300 obs, 9 var) and carolina3(64 obs, 5 var).
Carolina1: Data for the analysis of Carolina I Genetic design. In this design F2 or any advanced generation maintained by random mating, produced from cross between two pure-lines, is taken as base population. From the population an individual is randomly selected and used as a male. A set of 4 randomly selected plans are used as females and are mated to the above male. Thus a set of 4 full-sib families are produced. This is denoted as a male group. Similarly, a large number of male groups are produced. No female is used for any second mating. four male groups (16 female groups) from a set.
Carolina2: Data for the analysis of Carolina II Genetic design. Both paternal and maternal half-sibs are produced in this design. From an F2 population, n1 males and n2 females are randomly selected and each male is crossed to each of the females. Thus n1 x n2 progenies are produced whitch are analysed in a suitably laid experiment.
Carolina3: Data for the analysis of Carolina III genetic design. The F2 population is produced by crossing two inbreds, say L1 and L2. The material for estimation of genetic parameters is produced by back crossing randomly selected F2 individuals (using as males) to each of the inbreds (used as females).
Biometrical Methods in Quantitative Genetic Analysis, Singh, Chaudhary. 1979.
Biometrical Methods in Quantitative Genetic Analysis, Singh, Chaudhary. 1979.
data(DC) names(DC) str(DC$carolina1) str(DC$carolina2) str(DC$carolina3)data(DC) names(DC) str(DC$carolina1) str(DC$carolina2) str(DC$carolina3)
In many situations it is required to omit the rows or columns less or greater with NA of the matrix.
delete.na(x, alternative=c("less", "greater") )delete.na(x, alternative=c("less", "greater") )
x |
matrix with NA |
alternative |
"less" or "greater" |
x |
matrix |
Felipe de Mendiburu
library(agricolae) x<-c(2,5,3,7,5,NA,8,0,4,3,NA,NA) dim(x)<-c(4,3) x # [,1] [,2] [,3] #[1,] 2 5 4 #[2,] 5 NA 3 #[3,] 3 8 NA #[4,] 7 0 NA delete.na(x,"less") # [,1] #[1,] 2 #[2,] 5 #[3,] 3 #[4,] 7 delete.na(x,"greater") # [,1] [,2] [,3] #[1,] 2 5 4library(agricolae) x<-c(2,5,3,7,5,NA,8,0,4,3,NA,NA) dim(x)<-c(4,3) x # [,1] [,2] [,3] #[1,] 2 5 4 #[2,] 5 NA 3 #[3,] 3 8 NA #[4,] 7 0 NA delete.na(x,"less") # [,1] #[1,] 2 #[2,] 5 #[3,] 3 #[4,] 7 delete.na(x,"greater") # [,1] [,2] [,3] #[1,] 2 5 4
It generates a design of blocks, randomize and latin square for combined n. factors uses the methods of number generation in R. The seed is by set.seed(seed, kinds).
design.ab(trt, r, serie = 2, design=c("rcbd","crd","lsd"), seed = 0, kinds = "Super-Duper",first=TRUE,randomization=TRUE)design.ab(trt, r, serie = 2, design=c("rcbd","crd","lsd"), seed = 0, kinds = "Super-Duper",first=TRUE,randomization=TRUE)
trt |
n levels factors |
r |
Replications or Blocks |
serie |
number plot, 1: 11,12; 2: 101,102; 3: 1001,1002 |
design |
type |
seed |
Seed |
kinds |
Method for to randomize |
first |
TRUE or FALSE - randomize rep 1 |
randomization |
TRUE or FALSE - randomize |
kinds <- c("Wichmann-Hill", "Marsaglia-Multicarry", "Super-Duper", "Mersenne-Twister", "Knuth-TAOCP", "user-supplied", "Knuth-TAOCP-2002", "default" )
parameters |
Design parameters |
book |
Fieldbook |
Felipe de Mendiburu
Introduction to Experimental Statistics. Ching Chun Li. McGraw-Hill Book Company, INC, New. York, 1964
design.split, design.alpha,design.bib,
design.crd , design.cyclic , design.dau ,
design.graeco, design.lattice, design.lsd,
design.rcbd, design.strip
# factorial 3 x 2 with 3 blocks library(agricolae) trt<-c(3,2) # factorial 3x2 outdesign <-design.ab(trt, r=3, serie=2) book<-outdesign$book head(book,10) # print of the field book # factorial 2 x 2 x 2 with 5 replications in completely randomized design. trt<-c(2,2,2) outdesign<-design.ab(trt, r=5, serie=2,design="crd") book<-outdesign$book print(book) # factorial 3 x 3 in latin square design. trt <-c(3,3) outdesign<-design.ab(trt, serie=2, design="lsd") book<-outdesign$book print(book)# factorial 3 x 2 with 3 blocks library(agricolae) trt<-c(3,2) # factorial 3x2 outdesign <-design.ab(trt, r=3, serie=2) book<-outdesign$book head(book,10) # print of the field book # factorial 2 x 2 x 2 with 5 replications in completely randomized design. trt<-c(2,2,2) outdesign<-design.ab(trt, r=5, serie=2,design="crd") book<-outdesign$book print(book) # factorial 3 x 3 in latin square design. trt <-c(3,3) outdesign<-design.ab(trt, serie=2, design="lsd") book<-outdesign$book print(book)
Generates an alpha designs starting from the alpha design fixing under the series formulated by Patterson and Williams. These designs are generated by the alpha arrangements. They are similar to the lattice designs, but the tables are rectangular s by k (with s blocks and k<s columns. The number of treatments should be equal to s*k and all the experimental units r*s*k (r replications).
design.alpha(trt, k, r, serie = 2, seed = 0, kinds = "Super-Duper",randomization=TRUE)design.alpha(trt, k, r, serie = 2, seed = 0, kinds = "Super-Duper",randomization=TRUE)
trt |
Treatments |
k |
size block |
r |
Replications |
serie |
number plot, 1: 11,12; 2: 101,102; 3: 1001,1002 |
seed |
seed |
kinds |
method for to randomize |
randomization |
TRUE or FALSE - randomize |
Parameters for the alpha design: I. r=2, k <= s; II. r=3, s odd, k <= s; III.r=3, s even, k <= s-1; IV. r=4, s odd but not a multiple of 3, k<=s
r= replications s=number of blocks k=size of block Number of treatment is equal to k*s
parameters |
Design parameters |
statistics |
Design statistics |
sketch |
Design sketch |
book |
Fieldbook |
Felipe de Mendiburu
H.D. Patterson and E.R. Williams. Biometrika (1976) A new class of resolvable incomplete block designs. printed in Great Britain. Online: http://biomet.oxfordjournals.org/cgi/content/abstract/63/1/83
design.ab, design.split,design.bib,
design.crd , design.cyclic , design.dau ,
design.graeco, design.lattice, design.lsd,
design.rcbd, design.strip
library(agricolae) #Example one trt<-1:30 t <- length(trt) # size block k k<-3 # Blocks s s<-t/k # replications r r <- 2 outdesign<- design.alpha(trt,k,r,serie=2) book<-outdesign$book plots<-book[,1] dim(plots)<-c(k,s,r) for (i in 1:r) print(t(plots[,,i])) outdesign$sketch # Example two trt<-letters[1:12] t <- length(trt) k<-3 r<-3 s<-t/k outdesign<- design.alpha(trt,k,r,serie=2) book<-outdesign$book plots<-book[,1] dim(plots)<-c(k,s,r) for (i in 1:r) print(t(plots[,,i])) outdesign$sketchlibrary(agricolae) #Example one trt<-1:30 t <- length(trt) # size block k k<-3 # Blocks s s<-t/k # replications r r <- 2 outdesign<- design.alpha(trt,k,r,serie=2) book<-outdesign$book plots<-book[,1] dim(plots)<-c(k,s,r) for (i in 1:r) print(t(plots[,,i])) outdesign$sketch # Example two trt<-letters[1:12] t <- length(trt) k<-3 r<-3 s<-t/k outdesign<- design.alpha(trt,k,r,serie=2) book<-outdesign$book plots<-book[,1] dim(plots)<-c(k,s,r) for (i in 1:r) print(t(plots[,,i])) outdesign$sketch
Creates Randomized Balanced Incomplete Block Design. "Random" uses the methods of number generation in R. The seed is by set.seed(seed, kinds).
design.bib(trt, k, r=NULL, serie = 2, seed = 0, kinds = "Super-Duper", maxRep=20,randomization=TRUE)design.bib(trt, k, r=NULL, serie = 2, seed = 0, kinds = "Super-Duper", maxRep=20,randomization=TRUE)
trt |
Treatments |
k |
size block |
r |
Replications |
serie |
number plot, 1: 11,12; 2: 101,102; 3: 1001,1002 |
seed |
seed |
kinds |
method for to randomize |
maxRep |
repetition maximum |
randomization |
TRUE or FALSE - randomize |
The package AlgDesign is necessary.
if r = NULL, then it calculates the value of r smaller for k defined. In the case of r = value, then the possible values for "r" is calculated
K is the smallest integer number of treatments and both values are consistent in design.
kinds <- c("Wichmann-Hill", "Marsaglia-Multicarry", "Super-Duper", "Mersenne-Twister", "Knuth-TAOCP", "user-supplied", "Knuth-TAOCP-2002", "default" )
parameters |
Design parameters |
statistics |
Design statistics |
sketch |
Design sketch |
book |
Fieldbook |
Felipe de Mendiburu
1. Experimental design. Cochran and Cox. Second edition. Wiley Classics Library Edition published 1992
2. Optimal Experimental Design with R. Dieter Rasch, Jurgen Pilz, Rob Verdooren and Albrecht Gebhardt. 2011 by Taylor and Francis Group, LLC CRC Press is an imprint of Taylor and Francis Group, an Informa business.
3. Design of Experiments. Robert O. Kuehl. 2nd ed., Duxbury, 2000.
design.ab, design.alpha,design.split,
design.crd , design.cyclic , design.dau ,
design.graeco, design.lattice, design.lsd,
design.rcbd, design.strip
library(agricolae) # 4 treatments and k=3 size block trt<-c("A","B","C","D") k<-3 outdesign<-design.bib(trt,k,serie=2,seed =41,kinds ="Super-Duper") # seed = 41 print(outdesign$parameters) book<-outdesign$book plots <-as.numeric(book[,1]) matrix(plots,byrow=TRUE,ncol=k) print(outdesign$sketch) # write in hard disk # write.csv(book,"book.csv", row.names=FALSE) # file.show("book.csv")library(agricolae) # 4 treatments and k=3 size block trt<-c("A","B","C","D") k<-3 outdesign<-design.bib(trt,k,serie=2,seed =41,kinds ="Super-Duper") # seed = 41 print(outdesign$parameters) book<-outdesign$book plots <-as.numeric(book[,1]) matrix(plots,byrow=TRUE,ncol=k) print(outdesign$sketch) # write in hard disk # write.csv(book,"book.csv", row.names=FALSE) # file.show("book.csv")
It generates completely a randomized design with equal or different repetition. "Random" uses the methods of number generation in R. The seed is by set.seed(seed, kinds).
design.crd(trt, r, serie = 2, seed = 0, kinds = "Super-Duper",randomization=TRUE)design.crd(trt, r, serie = 2, seed = 0, kinds = "Super-Duper",randomization=TRUE)
trt |
Treatments |
r |
Replications |
serie |
number plot, 1: 11,12; 2: 101,102; 3: 1001,1002 |
seed |
seed |
kinds |
method for to randomize |
randomization |
TRUE or FALSE - randomize |
kinds <- c("Wichmann-Hill", "Marsaglia-Multicarry", "Super-Duper", "Mersenne-Twister", "Knuth-TAOCP", "user-supplied", "Knuth-TAOCP-2002", "default" )
parameters |
Design parameters |
book |
Fieldbook |
Felipe de Mendiburu
Introduction to Experimental Statistics. Ching Chun Li. McGraw-Hill Book Company, INC, New. York, 1964
design.ab, design.alpha,design.bib,
design.split , design.cyclic , design.dau ,
design.graeco, design.lattice, design.lsd,
design.rcbd, design.strip
library(agricolae) trt <-c("CIP-101","CIP-201","CIP-301","CIP-401","CIP-501") r <-c(4,3,5,4,3) # seed = 12543 outdesign1 <-design.crd(trt,r,serie=2,2543,"Mersenne-Twister") book1<-outdesign1 # no seed outdesign2 <-design.crd(trt,r,serie=3) print(outdesign2$parameters) book2<-outdesign2 # write to hard disk # write.table(book1,"crd.txt", row.names=FALSE, sep="\t") # file.show("crd.txt")library(agricolae) trt <-c("CIP-101","CIP-201","CIP-301","CIP-401","CIP-501") r <-c(4,3,5,4,3) # seed = 12543 outdesign1 <-design.crd(trt,r,serie=2,2543,"Mersenne-Twister") book1<-outdesign1 # no seed outdesign2 <-design.crd(trt,r,serie=3) print(outdesign2$parameters) book2<-outdesign2 # write to hard disk # write.table(book1,"crd.txt", row.names=FALSE, sep="\t") # file.show("crd.txt")
The cyclic design is a incomplete blocks designs, it is generated from a incomplete block initial of the size k, the plan is generated and randomized. The efficient and robust cyclic designs for 6 to 30 treatments, replications <= 10.
design.cyclic(trt, k, r, serie = 2, rowcol = FALSE, seed = 0, kinds = "Super-Duper" ,randomization=TRUE)design.cyclic(trt, k, r, serie = 2, rowcol = FALSE, seed = 0, kinds = "Super-Duper" ,randomization=TRUE)
trt |
vector treatments |
k |
block size |
r |
Replications |
serie |
number plot, 1: 11,12; 2: 101,102; 3: 1001,1002 |
rowcol |
TRUE: row-column design |
seed |
init seed random |
kinds |
random method |
randomization |
TRUE or FALSE - randomize |
Number o treatment 6 to 30. (r) Replication 2 to 10. (k) size of block 2 to 10. replication = i*k, "i" is value integer.
parameters |
Design parameters |
sketch |
Design sketch |
book |
Fieldbook |
Felipe de Mendiburu
Kuehl, Robert(2000), Design of Experiments. 2nd ed., Duxbury. John, J.A. (1981) Efficient Cyclic Design. J. R. Statist. Soc. B, 43, No. 1, pp, 76-80.
design.ab, design.alpha,design.bib,
design.crd , design.split , design.dau ,
design.graeco, design.lattice, design.lsd,
design.rcbd, design.strip
library(agricolae) trt<-letters[1:8] # block size = 2, replication = 6 outdesign1 <- design.cyclic(trt,k=2, r=6,serie=2) names(outdesign1) # groups 1,2,3 outdesign1$sketch[[1]] outdesign1$sketch[[2]] outdesign1$sketch[[3]] outdesign1$book # row-column design outdesign2<- design.cyclic(trt,k=2, r=6, serie=2, rowcol=TRUE) outdesign2$sketchlibrary(agricolae) trt<-letters[1:8] # block size = 2, replication = 6 outdesign1 <- design.cyclic(trt,k=2, r=6,serie=2) names(outdesign1) # groups 1,2,3 outdesign1$sketch[[1]] outdesign1$sketch[[2]] outdesign1$sketch[[3]] outdesign1$book # row-column design outdesign2<- design.cyclic(trt,k=2, r=6, serie=2, rowcol=TRUE) outdesign2$sketch
These are designs for two types of treatments: the control treatments (common) and the increased treatments. The common treatments are applied in complete randomized blocks, and the increased treatments, at random. Each treatment should be applied in any block once only. It is understood that the common treatments are of a greater interest; the standard error of the difference is much smaller than when between two increased ones in different blocks.
design.dau(trt1, trt2, r, serie = 2, seed = 0, kinds = "Super-Duper", name="trt" ,randomization=TRUE)design.dau(trt1, trt2, r, serie = 2, seed = 0, kinds = "Super-Duper", name="trt" ,randomization=TRUE)
trt1 |
checks |
trt2 |
new |
r |
Replications or blocks |
serie |
number plot, 1: 11,12; 2: 101,102; 3: 1001,1002 |
seed |
seed |
kinds |
method for to randomize |
name |
name of treatments |
randomization |
TRUE or FALSE - randomize |
kinds <- c("Wichmann-Hill", "Marsaglia-Multicarry", "Super-Duper", "Mersenne-Twister", "Knuth-TAOCP", "user-supplied", "Knuth-TAOCP-2002", "default" )
parameters |
Design parameters |
book |
Fieldbook |
Felipe de Mendiburu
1. Augmented (or Hoonuiaku) Design. Federer, W.T. (1956), Hawaii Plr. rec., 55: 191-208. 2. In Augmented Designs. Federer, W.T and Raghavarao, D. (1975). Bometrics, vol. 31, No. 1 (mar.., 1975), pp. 29-35
design.ab, design.alpha,design.bib,
design.crd , design.cyclic , design.split ,
design.graeco, design.lattice, design.lsd,
design.rcbd, design.strip
library(agricolae) # 4 treatments and 5 blocks T1<-c("A","B","C","D") T2<-letters[20:26] outdesign <-design.dau(T1,T2, r=5,serie=2) # field book book<-outdesign$book by(book,book[2],function(x) paste(x[,1],"-",as.character(x[,3]))) # write in hard disk # write.table(book,"dau.txt", row.names=FALSE, sep="\t") # file.show("dau.txt") # Augmented designs in Completely Randomized Design trt<-c(T1,T2) r<-c(4,4,4,4,1,1,1,1,1,1,1) outdesign <- design.crd(trt,r) outdesign$booklibrary(agricolae) # 4 treatments and 5 blocks T1<-c("A","B","C","D") T2<-letters[20:26] outdesign <-design.dau(T1,T2, r=5,serie=2) # field book book<-outdesign$book by(book,book[2],function(x) paste(x[,1],"-",as.character(x[,3]))) # write in hard disk # write.table(book,"dau.txt", row.names=FALSE, sep="\t") # file.show("dau.txt") # Augmented designs in Completely Randomized Design trt<-c(T1,T2) r<-c(4,4,4,4,1,1,1,1,1,1,1) outdesign <- design.crd(trt,r) outdesign$book
A graeco - latin square is a KxK pattern that permits the study of k treatments simultaneously with three different blocking variables, each at k levels.
The function is only for squares of the odd numbers and even numbers (4, 8, 10 and 12)
design.graeco(trt1, trt2, serie = 2, seed = 0, kinds = "Super-Duper",randomization=TRUE)design.graeco(trt1, trt2, serie = 2, seed = 0, kinds = "Super-Duper",randomization=TRUE)
trt1 |
Treatments |
trt2 |
Treatments |
serie |
number plot, 1: 11,12; 2: 101,102; 3: 1001,1002 |
seed |
seed |
kinds |
method for to randomize |
randomization |
TRUE or FALSE - randomize |
kinds <- c("Wichmann-Hill", "Marsaglia-Multicarry", "Super-Duper", "Mersenne-Twister", "Knuth-TAOCP", "user-supplied", "Knuth-TAOCP-2002", "default" )
parameters |
Design parameters |
book |
Fieldbook |
Felipe de Mendiburu
1. Statistics for Experimenters Design, Innovation, and Discovery Second Edition. George E. P. Box. Wiley-Interscience. 2005.
2. Experimental design. Cochran and Cox. Second edition. Wiley Classics Library Edition published 1992.
design.ab, design.alpha,design.bib,
design.crd , design.cyclic , design.dau ,
design.split, design.lattice, design.lsd,
design.rcbd, design.strip
library(agricolae) T1<-c("a","b","c","d") T2<-c("v","w","x","y") outdesign <- design.graeco(T1,T2,serie=1) graeco<-outdesign$book plots <-as.numeric(graeco[,1]) print(outdesign$sketch) print(matrix(plots,byrow=TRUE,ncol=4)) # 10 x 10 T1 <- letters[1:10] T2 <- 1:10 outdesign <- design.graeco(T1,T2,serie=2) print(outdesign$sketch)library(agricolae) T1<-c("a","b","c","d") T2<-c("v","w","x","y") outdesign <- design.graeco(T1,T2,serie=1) graeco<-outdesign$book plots <-as.numeric(graeco[,1]) print(outdesign$sketch) print(matrix(plots,byrow=TRUE,ncol=4)) # 10 x 10 T1 <- letters[1:10] T2 <- 1:10 outdesign <- design.graeco(T1,T2,serie=2) print(outdesign$sketch)
SIMPLE and TRIPLE lattice designs. It randomizes treatments in k x k lattice.
design.lattice(trt, r=3, serie = 2, seed = 0, kinds = "Super-Duper",randomization=TRUE)design.lattice(trt, r=3, serie = 2, seed = 0, kinds = "Super-Duper",randomization=TRUE)
trt |
treatments |
r |
r=2(simple) or r=3(triple) lattice |
serie |
number plot, 1: 11,12; 2: 101,102; 3: 1001,1002 |
seed |
seed |
kinds |
method for to randomize |
randomization |
TRUE or FALSE - randomize |
kinds <- c("Wichmann-Hill", "Marsaglia-Multicarry", "Super-Duper", "Mersenne-Twister", "Knuth-TAOCP", "user-supplied", "Knuth-TAOCP-2002", "default" )
parameters |
Design parameters |
statistics |
Design statistics |
sketch |
Design sketch |
book |
Fieldbook |
Felipe de Mendiburu
FIELD PLOT TECHNIQUE. Erwin L. LeCLERG. 2nd ed., 1962, Burgess Publishing Company, Minnesota
design.ab, design.alpha,design.bib,
design.crd , design.cyclic , design.dau ,
design.graeco, design.split, design.lsd,
design.rcbd, design.strip
library(agricolae) # triple lattice trt<-LETTERS[1:9] outdesign<-design.lattice(trt,r=3,serie=2) # triple lattice design ( 9 trt) # simple lattice trt<-1:100 outdesign<-design.lattice(trt,r=2,serie=3) # simple lattice design, 10x10library(agricolae) # triple lattice trt<-LETTERS[1:9] outdesign<-design.lattice(trt,r=3,serie=2) # triple lattice design ( 9 trt) # simple lattice trt<-1:100 outdesign<-design.lattice(trt,r=2,serie=3) # simple lattice design, 10x10
It generates Latin Square Design. "Random" uses the methods of number generation in R. The seed is by set.seed(seed, kinds).
design.lsd(trt, serie = 2, seed = 0, kinds = "Super-Duper",first=TRUE,randomization=TRUE)design.lsd(trt, serie = 2, seed = 0, kinds = "Super-Duper",first=TRUE,randomization=TRUE)
trt |
Treatments |
serie |
number plot, 1: 11,12; 2: 101,102; 3: 1001,1002 |
seed |
seed |
kinds |
method for to randomize |
first |
TRUE or FALSE - randomize rep 1 |
randomization |
TRUE or FALSE - randomize |
kinds <- c("Wichmann-Hill", "Marsaglia-Multicarry", "Super-Duper", "Mersenne-Twister", "Knuth-TAOCP", "user-supplied", "Knuth-TAOCP-2002", "default" )
parameters |
Design parameters |
book |
Fieldbook |
Felipe de Mendiburu
Introduction to Experimental Statistics. Ching Chun Li. McGraw-Hill Book Company, INC, New. York, 1969
design.ab, design.alpha,design.bib,
design.crd , design.cyclic , design.dau ,
design.graeco, design.lattice, design.split,
design.rcbd, design.strip
library(agricolae) varieties<-c("perricholi","yungay","maria bonita","tomasa") outdesign <-design.lsd(varieties,serie=2,seed=23) lsd <- outdesign$book print(outdesign$sketch) print(lsd) # field book. plots <-as.numeric(lsd[,1]) print(matrix(plots,byrow = TRUE, ncol = 4)) # Write on hard disk. # write.table(lsd,"lsd.txt", row.names=FALSE, sep="\t") # file.show("lsd.txt")library(agricolae) varieties<-c("perricholi","yungay","maria bonita","tomasa") outdesign <-design.lsd(varieties,serie=2,seed=23) lsd <- outdesign$book print(outdesign$sketch) print(lsd) # field book. plots <-as.numeric(lsd[,1]) print(matrix(plots,byrow = TRUE, ncol = 4)) # Write on hard disk. # write.table(lsd,"lsd.txt", row.names=FALSE, sep="\t") # file.show("lsd.txt")
Generate the design matrix from the fieldbook generated by an experimental plan or a dataframe for analysis.
design.mat(book, locations)design.mat(book, locations)
book |
data frame or matrix, field book |
locations |
numeric, column position of the field book |
X is matrix design.
Felipe de Mendiburu
design.ab, design.alpha,design.bib,
design.crd , design.cyclic , design.split ,
design.graeco, design.lattice, design.lsd,
design.rcbd, design.strip, design.dau
# dataframe: data analysis library(agricolae) data(sweetpotato) X<-design.mat(sweetpotato,1) print(X) # fieldbook: RCBD design trt <- LETTERS[1:4] r<-3 plan<-design.rcbd(trt,r,seed=11) X<-design.mat(plan$book,2:3) print(X)# dataframe: data analysis library(agricolae) data(sweetpotato) X<-design.mat(sweetpotato,1) print(X) # fieldbook: RCBD design trt <- LETTERS[1:4] r<-3 plan<-design.rcbd(trt,r,seed=11) X<-design.mat(plan$book,2:3) print(X)
It generates Randomized Complete Block Design. "Random" uses the methods of number generation in R. The seed is by set.seed(seed, kinds).
design.rcbd(trt, r, serie = 2, seed = 0, kinds = "Super-Duper", first=TRUE, continue=FALSE,randomization=TRUE )design.rcbd(trt, r, serie = 2, seed = 0, kinds = "Super-Duper", first=TRUE, continue=FALSE,randomization=TRUE )
trt |
Treatments |
r |
Replications or blocks |
serie |
number plot, 1: 11,12; 2: 101,102; 3: 1001,1002 |
seed |
seed |
kinds |
method for to randomize |
first |
TRUE or FALSE - randomize rep 1 |
continue |
TRUE or FALSE, continuous numbering of plot |
randomization |
TRUE or FALSE - randomize |
kinds <- c("Wichmann-Hill", "Marsaglia-Multicarry", "Super-Duper", "Mersenne-Twister", "Knuth-TAOCP", "user-supplied", "Knuth-TAOCP-2002", "default" )
parameters |
Design parameters |
sketch |
Design sketch |
book |
Fieldbook |
Felipe de Mendiburu
Introduction to Experimental Statistics. Ching Chun Li. McGraw-Hill Book Company, INC, New. York, 1964
design.ab, design.alpha,design.bib,
design.crd , design.cyclic , design.dau ,
design.graeco, design.lattice, design.lsd,
design.split, design.strip
library(agricolae) # 5 treatments and 6 blocks trt<-c("A","B","C","D","E") outdesign <-design.rcbd(trt,6,serie=2,986,"Wichmann-Hill") # seed = 986 book <-outdesign$book # field book # write in hard disk # write.table(book,"rcbd.txt", row.names=FALSE, sep="\t") # file.show("rcbd.txt") # Plots in field model ZIGZAG fieldbook <- zigzag(outdesign) print(outdesign$sketch) print(matrix(fieldbook[,1],byrow=TRUE,ncol=5)) # continuous numbering of plot outdesign <-design.rcbd(trt,6,serie=0,continue=TRUE) head(outdesign$book)library(agricolae) # 5 treatments and 6 blocks trt<-c("A","B","C","D","E") outdesign <-design.rcbd(trt,6,serie=2,986,"Wichmann-Hill") # seed = 986 book <-outdesign$book # field book # write in hard disk # write.table(book,"rcbd.txt", row.names=FALSE, sep="\t") # file.show("rcbd.txt") # Plots in field model ZIGZAG fieldbook <- zigzag(outdesign) print(outdesign$sketch) print(matrix(fieldbook[,1],byrow=TRUE,ncol=5)) # continuous numbering of plot outdesign <-design.rcbd(trt,6,serie=0,continue=TRUE) head(outdesign$book)
It generates split plot design. "Random" uses the methods of number generation in R. The seed is by set.seed(seed, kinds).
design.split(trt1, trt2,r=NULL, design=c("rcbd","crd","lsd"),serie = 2, seed = 0, kinds = "Super-Duper", first=TRUE,randomization=TRUE)design.split(trt1, trt2,r=NULL, design=c("rcbd","crd","lsd"),serie = 2, seed = 0, kinds = "Super-Duper", first=TRUE,randomization=TRUE)
trt1 |
Treatments in Plots |
trt2 |
Treatments in Subplots |
r |
Replications or blocks |
design |
Experimental design |
serie |
number plot, 1: 11,12; 2: 101,102; 3: 1001,1002 |
seed |
seed |
kinds |
method for to randomize |
first |
TRUE or FALSE - randomize rep 1 |
randomization |
TRUE or FALSE - randomize |
kinds <- c("Wichmann-Hill", "Marsaglia-Multicarry", "Super-Duper", "Mersenne-Twister", "Knuth-TAOCP", "user-supplied", "Knuth-TAOCP-2002", "default" )
parameters |
Design parameters |
book |
Fieldbook |
Felipe de Mendiburu
Statistical Procedures for Agricultural Research. Kwanchai A. Gomez, Arturo A. Gomez. John Wiley & Sons, new York, 1984
design.ab, design.alpha,design.bib,
design.crd , design.cyclic , design.dau ,
design.graeco, design.lattice, design.lsd,
design.rcbd, design.strip
library(agricolae) # 4 treatments and 5 blocks in split-plot t1<-c("A","B","C","D") t2<-c(1,2,3) outdesign <-design.split(t1,t2,r=3,serie=2,seed=45,kinds ="Super-Duper")#seed=45 book<-outdesign$book# field book # write in hard disk # write.table(book,"book.txt", row.names=FALSE, sep="\t") # file.show("book.txt")library(agricolae) # 4 treatments and 5 blocks in split-plot t1<-c("A","B","C","D") t2<-c(1,2,3) outdesign <-design.split(t1,t2,r=3,serie=2,seed=45,kinds ="Super-Duper")#seed=45 book<-outdesign$book# field book # write in hard disk # write.table(book,"book.txt", row.names=FALSE, sep="\t") # file.show("book.txt")
It generates strip plot design. "Random" uses the methods of number generation in R. The seed is by set.seed(seed, kinds).
design.strip(trt1, trt2,r, serie = 2, seed = 0, kinds = "Super-Duper",randomization=TRUE)design.strip(trt1, trt2,r, serie = 2, seed = 0, kinds = "Super-Duper",randomization=TRUE)
trt1 |
Row treatments |
trt2 |
column treatments |
r |
Replications |
serie |
number plot, 1: 11,12; 2: 101,102; 3: 1001,1002 |
seed |
seed |
kinds |
method for to randomize |
randomization |
TRUE or FALSE - randomize |
kinds <- c("Wichmann-Hill", "Marsaglia-Multicarry", "Super-Duper", "Mersenne-Twister", "Knuth-TAOCP", "user-supplied", "Knuth-TAOCP-2002", "default" )
parameters |
Design parameters |
book |
Fieldbook |
Felipe de Mendiburu
Statistical Procedures for Agricultural Research. Kwanchai A. Gomez, Arturo A. Gomez. John Wiley & Sons, new York, 1984
design.ab, design.alpha,design.bib,
design.crd , design.cyclic , design.dau ,
design.graeco, design.lattice, design.lsd,
design.rcbd, design.split
library(agricolae) # 4 and 3 treatments and 3 blocks in strip-plot t1<-c("A","B","C","D") t2<-c(1,2,3) r<-3 outdesign <-design.strip(t1,t2,r, serie=2,seed=45,kinds ="Super-Duper") # seed = 45 book <-outdesign$book # field book # write in hard disk # write.table(book,"book.txt", row.names=FALSE, sep="\t") # file.show("book.txt")library(agricolae) # 4 and 3 treatments and 3 blocks in strip-plot t1<-c("A","B","C","D") t2<-c(1,2,3) r<-3 outdesign <-design.strip(t1,t2,r, serie=2,seed=45,kinds ="Super-Duper") # seed = 45 book <-outdesign$book # field book # write in hard disk # write.table(book,"book.txt", row.names=FALSE, sep="\t") # file.show("book.txt")
Such designs are referred to as Youden squares since they were introduced by Youden (1937) after Yates (1936) considered the special case of column equal to number treatment minus 1. "Random" uses the methods of number generation in R. The seed is by set.seed(seed, kinds).
design.youden(trt, r, serie = 2, seed = 0, kinds = "Super-Duper",first=TRUE ,randomization=TRUE)design.youden(trt, r, serie = 2, seed = 0, kinds = "Super-Duper",first=TRUE ,randomization=TRUE)
trt |
Treatments |
r |
Replications or number of columns |
serie |
number plot, 1: 11,12; 2: 101,102; 3: 1001,1002 |
seed |
seed |
kinds |
method for to randomize |
first |
TRUE or FALSE - randomize rep 1 |
randomization |
TRUE or FALSE - randomize |
kinds <- c("Wichmann-Hill", "Marsaglia-Multicarry", "Super-Duper", "Mersenne-Twister", "Knuth-TAOCP", "user-supplied", "Knuth-TAOCP-2002", "default" )
parameters |
Design parameters |
sketch |
Design sketch |
book |
Fieldbook |
Felipe de Mendiburu
Design and Analysis of experiment. Hinkelmann, Klaus and Kempthorne, Oscar. Wiley-Interscience. Copyright (2008) by John Wiley and Sons. Inc., Hoboken, new Yersy
design.ab, design.alpha,design.bib,
design.crd , design.cyclic , design.dau ,
design.graeco, design.lattice, design.split,
design.rcbd, design.strip, design.lsd
library(agricolae) varieties<-c("perricholi","yungay","maria bonita","tomasa") r<-3 outdesign <-design.youden(varieties,r,serie=2,seed=23) youden <- outdesign$book print(outdesign$sketch) plots <-as.numeric(youden[,1]) print(matrix(plots,byrow=TRUE,ncol=r)) print(youden) # field book. # Write on hard disk. # write.table(youden,"youden.txt", row.names=FALSE, sep="\t") # file.show("youden.txt")library(agricolae) varieties<-c("perricholi","yungay","maria bonita","tomasa") r<-3 outdesign <-design.youden(varieties,r,serie=2,seed=23) youden <- outdesign$book print(outdesign$sketch) plots <-as.numeric(youden[,1]) print(matrix(plots,byrow=TRUE,ncol=r)) print(youden) # field book. # Write on hard disk. # write.table(youden,"youden.txt", row.names=FALSE, sep="\t") # file.show("youden.txt")
It plots bars of the averages of treatments to compare. It uses the objects generated by a procedure of comparison like LSD (Fisher), duncan, Tukey (HSD), Student Newman Keul (SNK), Scheffe, Ryan, Einot and Gabriel and Welsch (REGW), Kruskal Wallis, Friedman and Waerden.
diffograph(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="",...)diffograph(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="",...)
x |
Object created by a test of comparison, group=FALSE |
main |
The main title (on top) |
color1 |
non significant color |
color2 |
significant color |
color3 |
center line color |
cex.axis |
parameters of the plot() |
las |
parameters of the plot() |
pch |
parameters of the plot() |
bty |
parameters of the plot() |
cex |
parameters of the plot() |
lwd |
parameters of the plot() |
xlab |
parameters of the plot() |
ylab |
parameters of the plot() |
... |
Other parameters of the function plot() |
The graph.diff function should be used for functions: LSD, duncan, SNK, scheffe, REGW, HSD, kruskal, friedman and waerden test.
x |
list, object comparison test |
Felipe de Mendiburu
Multiple comparisons theory and methods. Departament of statistics the Ohio State University. USA, 1996. Jason C. Hsu. Chapman Hall/CRC
LSD.test, HSD.test, duncan.test, SNK.test,
scheffe.test, REGW.test, kruskal,friedman,
waerden.test
# Example 1 library(agricolae) data(sweetpotato) model<-aov(yield~virus,data=sweetpotato) x<- LSD.test(model,"virus",alpha=0.01,group=FALSE) diffograph(x,cex.axis=0.8,xlab="Yield",ylab="") # Example 2 x<- REGW.test(model,"virus",alpha=0.01,group=FALSE) diffograph(x,cex.axis=0.6,xlab="Yield",ylab="",color1="brown",color2="green")# Example 1 library(agricolae) data(sweetpotato) model<-aov(yield~virus,data=sweetpotato) x<- LSD.test(model,"virus",alpha=0.01,group=FALSE) diffograph(x,cex.axis=0.8,xlab="Yield",ylab="") # Example 2 x<- REGW.test(model,"virus",alpha=0.01,group=FALSE) diffograph(x,cex.axis=0.6,xlab="Yield",ylab="",color1="brown",color2="green")
Three evaluations over time and the potato yield when applying several treatments.
data(disease)data(disease)
A data frame with 21 observations on the following 7 variables.
plotsa numeric vector
repa numeric vector
trta factor with levels T0 T1 T2 T3 T4 T5 T6
E2a numeric vector
E5a numeric vector
E7a numeric vector
yielda numeric vector
Experimental data.
International Potato Center. CIP - Lima Peru.
library(agricolae) data(disease) str(disease)library(agricolae) data(disease) str(disease)
This test is adapted from the Newman-Keuls method. Duncan's test does not control family wise error rate at the specified alpha level. It has more power than the other post tests, but only because it doesn't control the error rate properly. The Experimentwise Error Rate at: 1-(1-alpha)^(a-1); where "a" is the number of means and is the Per-Comparison Error Rate. Duncan's procedure is only very slightly more conservative than LSD. The level by alpha default is 0.05.
duncan.test(y, trt, DFerror, MSerror, alpha = 0.05, group=TRUE, main = NULL,console=FALSE)duncan.test(y, trt, DFerror, MSerror, alpha = 0.05, group=TRUE, main = NULL,console=FALSE)
y |
model(aov or lm) or answer of the experimental unit |
trt |
Constant( only y=model) or vector treatment applied to each experimental unit |
DFerror |
Degree free |
MSerror |
Mean Square Error |
alpha |
Significant level |
group |
TRUE or FALSE |
main |
Title |
console |
logical, print output |
It is necessary first makes a analysis of variance.
if y = model, then to apply the instruction:
duncan.test(model, "trt", alpha = 0.05, group = TRUE, main = NULL, console = FALSE)
where the model class is aov or lm.
statistics |
Statistics of the model |
parameters |
Design parameters |
duncan |
Critical Range Table |
means |
Statistical summary of the study variable |
comparison |
Comparison between treatments |
groups |
Formation of treatment groups |
Felipe de Mendiburu
1. Principles and procedures of statistics a biometrical approach
Steel & Torry & Dickey. Third Edition 1997
2. Multiple comparisons theory and methods. Departament of statistics
the Ohio State University. USA, 1996. Jason C. Hsu. Chapman Hall/CRC.
BIB.test, DAU.test, durbin.test,
friedman, HSD.test, kruskal,
LSD.test, Median.test, PBIB.test,
REGW.test, scheffe.test, SNK.test,
waerden.test, waller.test, plot.group
library(agricolae) data(sweetpotato) model<-aov(yield~virus,data=sweetpotato) out <- duncan.test(model,"virus", main="Yield of sweetpotato. Dealt with different virus") plot(out,variation="IQR") duncan.test(model,"virus",alpha=0.01,console=TRUE) # version old duncan.test() df<-df.residual(model) MSerror<-deviance(model)/df out <- with(sweetpotato,duncan.test(yield,virus,df,MSerror, group=TRUE)) plot(out,horiz=TRUE,las=1) print(out$groups)library(agricolae) data(sweetpotato) model<-aov(yield~virus,data=sweetpotato) out <- duncan.test(model,"virus", main="Yield of sweetpotato. Dealt with different virus") plot(out,variation="IQR") duncan.test(model,"virus",alpha=0.01,console=TRUE) # version old duncan.test() df<-df.residual(model) MSerror<-deviance(model)/df out <- with(sweetpotato,duncan.test(yield,virus,df,MSerror, group=TRUE)) plot(out,horiz=TRUE,las=1) print(out$groups)
A multiple comparison of the Durbin test for the balanced incomplete blocks for sensorial or categorical evaluation. It forms groups according to the demanded ones for level of significance (alpha); by default, 0.05.
durbin.test(judge, trt, evaluation, alpha = 0.05, group =TRUE, main = NULL, console=FALSE)durbin.test(judge, trt, evaluation, alpha = 0.05, group =TRUE, main = NULL, console=FALSE)
judge |
Identification of the judge in the evaluation |
trt |
Treatments |
evaluation |
variable |
alpha |
level of significant |
group |
TRUE or FALSE |
main |
Title |
console |
logical, print output |
The post hoc test is using the criterium Fisher's least significant difference.
statistics |
Statistics of the model |
parameters |
Design parameters |
means |
Statistical summary of the study variable |
rank |
rank table of the study variable |
comparison |
Comparison between treatments |
groups |
Formation of treatment groups |
Felipe de Mendiburu
Practical Nonparametrics Statistics. W.J. Conover, 1999 Nonparametric Statistical Methods. Myles Hollander and Douglas A. Wofe, 1999
BIB.test, DAU.test, duncan.test,
friedman, HSD.test, kruskal,
LSD.test, Median.test, PBIB.test,
REGW.test, scheffe.test, SNK.test,
waerden.test, waller.test, plot.group
library(agricolae) # Example 1. Conover, pag 391 person<-gl(7,3) variety<-c(1,2,4,2,3,5,3,4,6,4,5,7,1,5,6,2,6,7,1,3,7) preference<-c(2,3,1,3,1,2,2,1,3,1,2,3,3,1,2,3,1,2,3,1,2) out<-durbin.test(person,variety,preference,group=TRUE,console=TRUE, main="Seven varieties of ice cream manufacturer") #startgraph bar.group(out$groups,horiz=TRUE,xlim=c(0,10),density=4,las=1) #endgraph # Example 2. Myles Hollander, pag 311 # Source: W. Moore and C.I. Bliss. 1942 day<-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) out<-durbin.test(day,chemical,toxic,group=TRUE,console=TRUE, main="Logarithm of Toxic Dosages") plot(out)library(agricolae) # Example 1. Conover, pag 391 person<-gl(7,3) variety<-c(1,2,4,2,3,5,3,4,6,4,5,7,1,5,6,2,6,7,1,3,7) preference<-c(2,3,1,3,1,2,2,1,3,1,2,3,3,1,2,3,1,2,3,1,2) out<-durbin.test(person,variety,preference,group=TRUE,console=TRUE, main="Seven varieties of ice cream manufacturer") #startgraph bar.group(out$groups,horiz=TRUE,xlim=c(0,10),density=4,las=1) #endgraph # Example 2. Myles Hollander, pag 311 # Source: W. Moore and C.I. Bliss. 1942 day<-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) out<-durbin.test(day,chemical,toxic,group=TRUE,console=TRUE, main="Logarithm of Toxic Dosages") plot(out)
The data consist of b-blocks mutually independent k-variate random variables Xij, i=1,..,b; j=1,..k. The random variable X is in block i and is associated with treatment j. It makes the multiple comparison of the Friedman test with or without ties. A first result is obtained by friedman.test of R.
friedman(judge,trt,evaluation,alpha=0.05,group=TRUE,main=NULL,console=FALSE)friedman(judge,trt,evaluation,alpha=0.05,group=TRUE,main=NULL,console=FALSE)
judge |
Identification of the judge in the evaluation |
trt |
Treatment |
evaluation |
Variable |
alpha |
Significant test |
group |
TRUE or FALSE |
main |
Title |
console |
logical, print output |
The post hoc friedman test is using the criterium Fisher's least significant difference (LSD)
statistics |
Statistics of the model |
parameters |
Design parameters |
means |
Statistical summary of the study variable |
comparison |
Comparison between treatments |
groups |
Formation of treatment groups |
Felipe de Mendiburu
Practical Nonparametrics Statistics. W.J. Conover, 1999
BIB.test, DAU.test, duncan.test,
durbin.test, HSD.test, kruskal,
LSD.test, Median.test, PBIB.test,
REGW.test, scheffe.test, SNK.test,
waerden.test, waller.test, plot.group
library(agricolae) data(grass) out<-with(grass,friedman(judge,trt, evaluation,alpha=0.05, group=TRUE,console=TRUE, main="Data of the book of Conover")) #startgraph plot(out,variation="IQR") #endgraphlibrary(agricolae) data(grass) out<-with(grass,friedman(judge,trt, evaluation,alpha=0.05, group=TRUE,console=TRUE, main="Data of the book of Conover")) #startgraph plot(out,variation="IQR") #endgraph
Data of frijol under 4 technologies for the homogeneity of regression study. Yield of Frijol in kg/ha in clean and dry grain.
Tecnnologies: 20-40-20 kg/ha. N. P2O5 and K2O + 2 t/ha of gallinaza. 40-80-40 kg/ha. N. P2O5 and K2O + 2 t/ha of gallinaza. 60-120-60 kg/ha. N. P2O5 and K2O + 2 t/ha of gallinaza. 40-80-40 kg/ha. N. P2O5 and K2O + 4 t/ha of gallinaza.
data(frijol)data(frijol)
A data frame with 84 observations on the following 3 variables.
technologya factor with levels a b c d
productiona numeric vector
indexa numeric vector
Oriente antioqueno (1972) (ICA.- Orlando Martinez W.) Colombia.
library(agricolae) data(frijol) str(frijol)library(agricolae) data(frijol) str(frijol)
50 genotypes and 5 environments.
data(genxenv)data(genxenv)
A data frame with 250 observations on the following 3 variables.
ENVa numeric vector
GENa numeric vector
YLDa numeric vector
International Potato Center. CIP - Lima Peru.
International Potato Center. CIP - Lima Peru.
library(agricolae) data(genxenv) str(genxenv)library(agricolae) data(genxenv) str(genxenv)
A measurement of the Glycoalkaloids by two methods: HPLC and spectrophotometer.
data(Glycoalkaloids)data(Glycoalkaloids)
A data frame with 25 observations on the following 2 variables.
HPLCa numeric vector
spectrophotometera numeric vector
International Potato Center. CIP - Lima Peru.
International Potato Center. CIP - Lima Peru.
library(agricolae) data(Glycoalkaloids) str(Glycoalkaloids)library(agricolae) data(Glycoalkaloids) str(Glycoalkaloids)
In many situations it has intervals of class defined with its respective frequencies. By means of this function, the graphic of frequency is obtained and it is possible to superpose the normal distribution, polygon of frequency, Ojiva and to construct the table of complete frequency.
graph.freq(x, breaks=NULL,counts=NULL,frequency=1, plot=TRUE, nclass=NULL, xlab="",ylab="",axes = "",las=1,...)graph.freq(x, breaks=NULL,counts=NULL,frequency=1, plot=TRUE, nclass=NULL, xlab="",ylab="",axes = "",las=1,...)
x |
a vector of values, a object hist(), graph.freq() |
counts |
frequency and x is class intervals |
breaks |
a vector giving the breakpoints between histogram cells |
frequency |
1=counts, 2=relative, 3=density |
plot |
logic |
nclass |
number of classes |
xlab |
x labels |
ylab |
y labels |
las |
values 0,1,2 and 3 are the axis styles. see plot() |
axes |
TRUE or FALSE |
... |
other parameters of plot |
breaks |
a vector giving the breakpoints between histogram cells |
counts |
frequency and x is class intervals |
mids |
center point in class |
relative |
Relative frequency, height |
density |
Density frequency, height |
Felipe de Mendiburu
polygon.freq, table.freq,
stat.freq,inter.freq,sturges.freq,
join.freq,ogive.freq, normal.freq
library(agricolae) data(genxenv) yield <- subset(genxenv$YLD,genxenv$ENV==2) yield <- round(yield,1) h<- graph.freq(yield,axes=FALSE, frequency=1, ylab="frequency",col="yellow") axis(1,h$breaks) axis(2,seq(0,20,0.1)) # To reproduce histogram. h1 <- graph.freq(h, col="blue", frequency=2,border="red", density=8,axes=FALSE, xlab="YIELD",ylab="relative") axis(1,h$breaks) axis(2,seq(0,.4,0.1)) # summary, only frecuency limits <-seq(10,40,5) frequencies <-c(2,6,8,7,3,4) #startgraph h<-graph.freq(limits,counts=frequencies,col="bisque",xlab="Classes") polygon.freq(h,col="red") title( main="Histogram and polygon of frequency", ylab="frequency") #endgraph # Statistics measures<-stat.freq(h) print(measures) # frequency table full round(table.freq(h),2) #startgraph # ogive ogive.freq(h,col="red",type="b",ylab="Accumulated relative frequency", xlab="Variable") # only .frequency polygon h<-graph.freq(limits,counts=frequencies,border=FALSE,col=NULL,xlab=" ",ylab="") title( main="Polygon of frequency", xlab="Variable", ylab="Frecuency") polygon.freq(h,col="blue") grid(col="brown") #endgraph # Draw curve for Histogram h<- graph.freq(yield,axes=FALSE, frequency=3, ylab="f(yield)",col="yellow") axis(1,h$breaks) axis(2,seq(0,0.18,0.03),las=2) lines(density(yield), col = "red", lwd = 2) title("Draw curve for Histogram")library(agricolae) data(genxenv) yield <- subset(genxenv$YLD,genxenv$ENV==2) yield <- round(yield,1) h<- graph.freq(yield,axes=FALSE, frequency=1, ylab="frequency",col="yellow") axis(1,h$breaks) axis(2,seq(0,20,0.1)) # To reproduce histogram. h1 <- graph.freq(h, col="blue", frequency=2,border="red", density=8,axes=FALSE, xlab="YIELD",ylab="relative") axis(1,h$breaks) axis(2,seq(0,.4,0.1)) # summary, only frecuency limits <-seq(10,40,5) frequencies <-c(2,6,8,7,3,4) #startgraph h<-graph.freq(limits,counts=frequencies,col="bisque",xlab="Classes") polygon.freq(h,col="red") title( main="Histogram and polygon of frequency", ylab="frequency") #endgraph # Statistics measures<-stat.freq(h) print(measures) # frequency table full round(table.freq(h),2) #startgraph # ogive ogive.freq(h,col="red",type="b",ylab="Accumulated relative frequency", xlab="Variable") # only .frequency polygon h<-graph.freq(limits,counts=frequencies,border=FALSE,col=NULL,xlab=" ",ylab="") title( main="Polygon of frequency", xlab="Variable", ylab="Frecuency") polygon.freq(h,col="blue") grid(col="brown") #endgraph # Draw curve for Histogram h<- graph.freq(yield,axes=FALSE, frequency=3, ylab="f(yield)",col="yellow") axis(1,h$breaks) axis(2,seq(0,0.18,0.03),las=2) lines(density(yield), col = "red", lwd = 2) title("Draw curve for Histogram")
Twelve homeowners are selected randomly to participate in an experiment with a plant nursery. Each homeowner is asked to select four fairly identical areas in his yard and to plant four different types of grasses, one in each area.
data(grass)data(grass)
A data frame with 48 observations on the following 3 variables.
judgea numeric vector
trta factor with levels t1 t2 t3 t4
evaluationa numeric vector
Each of the 12 blocks consists of four fairly identical plots of land, each receiving care of approximately the same degree of skill because the four plots are presumably cared for by the same homeowern.
Book: Practical Nonparametrics Statistics, pag 372.
Practical Nonparametrics Statistics. W.J. Conover, 1999
data(grass) str(grass)data(grass) str(grass)
Potato minituber production in greenhouse, three sets of data in potato varieties with different methods: hydroponics, Aeroponic, Pots and Plant beds, the unit is in grams and the number of tubers in units,
data(greenhouse)data(greenhouse)
greenhouse is list, three tables: greenhouse1(480 obs, 5 var), yield for plant, unit is grams. greenhouse2(48 obs, 5 var), Yields of 10 plants by experimental unit(grams). planting date(April 24, 2004) and harvest date(July 16, 2004) and greenhouse3(480 obs, 5 var), Tubers by plants.
International Potato Center(CIP). Lima-Peru. Data Kindly provided by Carlos Chuquillanqui.
- Produccion de semila de papa por hidroponia tecnica de flujo continuo de una pelicula de solucion nutritiva (nft) Carlos Chuquillanqui(CIP), Jorge Tenorio(CIP) and L. F. Salazar(Agdia Inc). AGROENFOQUE Lima-Peru (2004) - Potato Minituber Production Using Aeroponics: Effect of Plant Density and Harvesting Intervals American Journal of Potato Research, Jan/Feb 2006 by Farran, Imma, Mingo-Castel, Angel M
library(agricolae) data(greenhouse) greenhouse1 <- greenhouse$greenhouse1 greenhouse2 <- greenhouse$greenhouse2 greenhouse3 <- greenhouse$greenhouse3library(agricolae) data(greenhouse) greenhouse1 <- greenhouse$greenhouse1 greenhouse2 <- greenhouse$greenhouse2 greenhouse3 <- greenhouse$greenhouse3
Data growth of pijuayo trees in several localities.
data(growth)data(growth)
A data frame with 30 observations on the following 3 variables.
placea factor with levels L1 L2
slimea numeric vector
heighta numeric vector
Experimental data (Pucallpa - Peru)
ICRAF lima Peru.
library(agricolae) data(growth) str(growth)library(agricolae) data(growth) str(growth)
Published data. Haynes. Mean area under the disease progress curve (AUDPC) for each of 16 potato clones evaluated at eight sites across the United States in 1996
data(haynes)data(haynes)
A data frame with 16 observations on the following 9 variables.
clonea factor with levels A84118-3 AO80432-1 AO84275-3 AWN86514-2 B0692-4 B0718-3 B0749-2F B0767-2 Bertita Bzura C0083008-1 Elba Greta Krantz Libertas Stobrawa
FLa numeric vector
MIa numeric vector
MEa numeric vector
MNa numeric vector
NDa numeric vector
NYa numeric vector
PAa numeric vector
WIa numeric vector
Haynes K G, Lambert D H, Christ B J, Weingartner D P, Douches D S, Backlund J E, Fry W and Stevenson W. 1998. Phenotypic stability of resistance to late blight in potato clones evaluated at eight sites in the United States American Journal Potato Research 75, pag 211-217.
library(agricolae) data(haynes) str(haynes)library(agricolae) data(haynes) str(haynes)
Incidents and performance of healthy tubers and rotten potato field infested with naturally Ralstonia solanacearum Race 3/Bv 2A, after application of inorganic amendments and a rotation crop in Huanuco Peru, 2006.
data(Hco2006)data(Hco2006)
The format is: List of 2
amendmenta factor
cropa factor
blocka numeric vector, replications
planta numeric vector, number plant
wilt_percenta numeric vector, wilt percentage at 60 days
healtha numeric vector, kg/8m2, 20 plants
rota numeric vector, kg/8m2, 20 plants
Application of inorganic amendment and crop rotation to control bacterial wilt of the potato (MBP).
Experimental field, 2006. Data Kindly provided by Pedro Aley.
International Potato Center. CIP - Lima Peru.
library(agricolae) data(Hco2006) str(Hco2006) wilt<-Hco2006$wilt yield<-Hco2006$yield means <- tapply.stat(wilt[,5],wilt[,1:3],function(x) mean(x,na.rm=TRUE)) names(means)[4]<-"wilt_percent" model <- aov(wilt_percent ~ block + crop, means) anova(model) cv.model(model) yield<-yield[order(paste(yield[,1],yield[,2],yield[,3])),] correlation(means[,4],yield[,4],method="spearman")library(agricolae) data(Hco2006) str(Hco2006) wilt<-Hco2006$wilt yield<-Hco2006$yield means <- tapply.stat(wilt[,5],wilt[,1:3],function(x) mean(x,na.rm=TRUE)) names(means)[4]<-"wilt_percent" model <- aov(wilt_percent ~ block + crop, means) anova(model) cv.model(model) yield<-yield[order(paste(yield[,1],yield[,2],yield[,3])),] correlation(means[,4],yield[,4],method="spearman")
It shows dendrogram of a consensus of a tree generated by hclust.
hcut(consensus,h,group,col.text="blue",cex.text=1,...)hcut(consensus,h,group,col.text="blue",cex.text=1,...)
consensus |
object consensus |
h |
numeric scalar or vector with heights where the tree should be cut. |
group |
an integer scalar with the desired number of group |
col.text |
color of number consensus |
cex.text |
size of number consensus |
... |
Other parameters of the function plot() in cut() |
hcut Returns a data frame with group memberships and consensus tree.
F. de Mendiburu
library(agricolae) data(pamCIP) # only code rownames(pamCIP)<-substr(rownames(pamCIP),1,6) # groups of clusters # output<-consensus(pamCIP,nboot=100) # hcut(output,h=0.4,group=5,main="Group 5") # # hcut(output,h=0.4,group=8,type="t",edgePar=list(lty=1:2,col=2:1),main="group 8" # ,col.text="blue",cex.text=1)library(agricolae) data(pamCIP) # only code rownames(pamCIP)<-substr(rownames(pamCIP),1,6) # groups of clusters # output<-consensus(pamCIP,nboot=100) # hcut(output,h=0.4,group=5,main="Group 5") # # hcut(output,h=0.4,group=8,type="t",edgePar=list(lty=1:2,col=2:1),main="group 8" # ,col.text="blue",cex.text=1)
Determination of heterosis, general combining ability (GCA) and specific combining ability in tuber dry matter, reducing sugars and agronomic characteristics in TPS families.
data(heterosis)data(heterosis)
A data frame with 216 observations on the following 11 variables.
Place1: La Molina, 2=Huancayo
Replicationa numeric vector
Treatmenta numeric vector
Factora factor with levels Control progenie progenitor testigo
Femalea factor with levels Achirana LT-8 MF-I MF-II Serrana TPS-2 TPS-25 TPS-7
Malea factor with levels TPS-13 TPS-67 TS-15
v1Yield (Kg/plant)
v2Reducing sugars (scale):1 low and 5=High
v3Tuber dry matter (percentage)
v4Tuber number/plant
v5Average tuber weight (g)
The study was conducted in 3 environments, La Molina-PERU to 240 masl. during autumn-winter and spring, and in Huancayo-PERU 3180 masl., during summer. The experimental material consisted of 24 families half brother in the form of tubers derived from TPS, obtained crossing between 8 female and 3 male parents. Design used was randomized complete block with three repetitions. The experimental unit was 30 plants in two rows at a distance of 30cm between plants and 90 cm between rows. Variables evaluated were Yield, Tubers number, Dry matter and content and reducing sugars. The analysis was conducted line x tester. The control variety was Desiree.
International Potato Center(CIP). Lima-Peru. Data Kindly provided by of Rolando Cabello.
Tesis "Heterosis, habilidad combinatoria general y especifica para materia seca, azucares reductores y caracteres agronomicos en familias de tuberculos provenientes de semilla sexual de papa. Magister Scientiae Rodolfo Valdivia Lorente. Universidad Nacional Agraria La molina-Lima Peru, Escuela de Post Grado, Mejoramiento genetico de plantas, 2004". Poster: Congreso de la Sociedad Peruana de Genetica - Peru, 2008.
library(agricolae) data(heterosis) str(heterosis) site1<-subset(heterosis,heterosis[,1]==1) site2<-subset(heterosis,heterosis[,1]==2) site3<-subset(heterosis,heterosis[,1]==3) model1<-with(site1,lineXtester(Replication, Female, Male, v1)) DFe <- df.residual(model1) CMe <- deviance(model1)/DFe test1 <- with(site1,HSD.test(v1, Factor,DFe,CMe)) test2 <- with(site1,HSD.test(v1, Treatment,DFe,CMe)) model22<-with(site2,lineXtester(Replication, Female, Male, v3)) model3<-with(site3,lineXtester(Replication, Female, Male, v4))library(agricolae) data(heterosis) str(heterosis) site1<-subset(heterosis,heterosis[,1]==1) site2<-subset(heterosis,heterosis[,1]==2) site3<-subset(heterosis,heterosis[,1]==3) model1<-with(site1,lineXtester(Replication, Female, Male, v1)) DFe <- df.residual(model1) CMe <- deviance(model1)/DFe test1 <- with(site1,HSD.test(v1, Factor,DFe,CMe)) test2 <- with(site1,HSD.test(v1, Treatment,DFe,CMe)) model22<-with(site2,lineXtester(Replication, Female, Male, v3)) model3<-with(site3,lineXtester(Replication, Female, Male, v4))
Returns a vector with group memberships. This function is used by the function consensus of clusters.
hgroups(hmerge)hgroups(hmerge)
hmerge |
The object is components of the hclust |
The merge clusters is printed.
F. de Mendiburu
library(agricolae) data(pamCIP) # only code rownames(pamCIP)<-substr(rownames(pamCIP),1,6) distance <- dist(pamCIP,method="binary") clusters<- hclust( distance, method="complete") # groups of clusters hgroups(clusters$merge)library(agricolae) data(pamCIP) # only code rownames(pamCIP)<-substr(rownames(pamCIP),1,6) distance <- dist(pamCIP,method="binary") clusters<- hclust( distance, method="complete") # groups of clusters hgroups(clusters$merge)
It makes multiple comparisons of treatments by means of Tukey. The level by alpha default is 0.05.
HSD.test(y, trt, DFerror, MSerror, alpha = 0.05, group=TRUE, main = NULL, unbalanced=FALSE,console=FALSE)HSD.test(y, trt, DFerror, MSerror, alpha = 0.05, group=TRUE, main = NULL, unbalanced=FALSE,console=FALSE)
y |
model(aov or lm) or answer of the experimental unit |
trt |
Constant( only y=model) or vector treatment applied to each experimental unit |
DFerror |
Degree free |
MSerror |
Mean Square Error |
alpha |
Significant level |
group |
TRUE or FALSE |
main |
Title |
unbalanced |
TRUE or FALSE. not equal replication |
console |
logical, print output |
It is necessary first makes a analysis of variance.
if y = model, then to apply the instruction:
HSD.test (model, "trt", alpha = 0.05, group = TRUE, main = NULL, unbalanced=FALSE, console=FALSE)
where the model class is aov or lm.
statistics |
Statistics of the model |
parameters |
Design parameters |
means |
Statistical summary of the study variable |
comparison |
Comparison between treatments |
groups |
Formation of treatment groups |
Felipe de Mendiburu
1. Principles and procedures of statistics a biometrical approach
Steel & Torry & Dickey. Third Edition 1997
2. Multiple comparisons theory and methods. Departament of statistics
the Ohio State University. USA, 1996. Jason C. Hsu. Chapman Hall/CRC.
BIB.test, DAU.test, duncan.test,
durbin.test, friedman, kruskal,
LSD.test, Median.test, PBIB.test,
REGW.test, scheffe.test, SNK.test,
waerden.test, waller.test, plot.group
library(agricolae) data(sweetpotato) model<-aov(yield~virus, data=sweetpotato) out <- HSD.test(model,"virus", group=TRUE,console=TRUE, main="Yield of sweetpotato\nDealt with different virus") #stargraph # Variation range: max and min plot(out) #endgraph out<-HSD.test(model,"virus", group=FALSE) print(out$comparison) # Old version HSD.test() df<-df.residual(model) MSerror<-deviance(model)/df with(sweetpotato,HSD.test(yield,virus,df,MSerror, group=TRUE,console=TRUE, main="Yield of sweetpotato. Dealt with different virus"))library(agricolae) data(sweetpotato) model<-aov(yield~virus, data=sweetpotato) out <- HSD.test(model,"virus", group=TRUE,console=TRUE, main="Yield of sweetpotato\nDealt with different virus") #stargraph # Variation range: max and min plot(out) #endgraph out<-HSD.test(model,"virus", group=FALSE) print(out$comparison) # Old version HSD.test() df<-df.residual(model) MSerror<-deviance(model)/df with(sweetpotato,HSD.test(yield,virus,df,MSerror, group=TRUE,console=TRUE, main="Yield of sweetpotato. Dealt with different virus"))
Timing fungicide sprays based on accumulated rainfall thresholds can be a successful component of integrated management packages that include cultivars with moderate or high levels of resistance to late blight. The simplicity of measuring accumulated rainfall means that the technology can potentially be used by resource-poor farmers in developing countries.
data(huasahuasi)data(huasahuasi)
The format is: List of 2 ( AUDPC, YIELD )
blocka factor with levels I II III
trta factor with levels 40mm 7-days Non-application
clona factor with levels C386209.10 C387164.4 Cruza148 Musuq Yungay
y1daa numeric vector, Kgr./plot
y2daa numeric vector, Kgr./plot
y3raa numeric vector, Kgr./plot
d44a numeric vector, 44 days
d51a numeric vector, 51 days
d100a numeric vector, 100 days
The experimental unit was formed by 4 furrows of 1.8 m of length, with distance between furrows from 0.90 m and between plants of 0.30 m. In each furrow was installed 5 plants. The experiment had 3 repetitions. From the beginning of the experiment were fulfilled the following treatments Thresholds 40 mm: Apply the fungicide when 40 precipitation mm accumulates. The minimum interval between applications will be of 7 days. Schedule 7 days: The applications should be carried out every 7 days calendar. Without application: No fungicide application will be made. The evaluation of the severity of the late blight in each treatment started to emergency 80 percentage and then evaluations were made every 7 days until being observed a physiological maturation of the crop.
Experimental field, 2003. Data Kindly provided by Wilmer Perez.
International Potato Center. CIP - Lima Peru.
library(agricolae) data(huasahuasi) names(huasahuasi) str(huasahuasi$AUDPC) str(huasahuasi$YIELD)library(agricolae) data(huasahuasi) names(huasahuasi) str(huasahuasi$AUDPC) str(huasahuasi$YIELD)
calculate AMMI stability value (ASV) and Yield stability index (YSI).
index.AMMI(model)index.AMMI(model)
model |
object AMMI |
AMMI stability value (ASV) was calculated using the following formula, as suggested by Purchase (1997)
ASV = sqrt((SSpc1/SSpc2 * PC1i)^2+(PC2i)^2)
YSI = RASV + RY
RASV = rank(ASV) and RY = rank(Y across by environment)
ASV |
AMMI stability value |
YSI |
Yield stability index |
rASV |
Rank of AMMI stability value |
rYSI |
Rank of yield stability index |
means |
average genotype by environment |
F. de Mendiburu
The use of an AMMI model and its parameters to analyse yield stability in multienvironment trials. N. SABAGHNIA, S.H. SABAGHPOUR AND H. DEHGHANI. Journal of Agricultural Science (2008), 146, 571-581. f 2008 Cambridge University Press 571 doi:10.1017/S0021859608007831 Printed in the United Kingdom
Parametric analysis to describe genotype x environment interaction and yield stability in winter wheat. PURCHASE, J. L. (1997). Ph.D. Thesis, Department of Agronomy, Faculty of Agriculture of the University of the Free State, Bloemfontein, South Africa.
library(agricolae) # Index AMMI data(plrv) model<- with(plrv,AMMI(Locality, Genotype, Rep, Yield, console=FALSE)) Idx<-index.AMMI(model) names(Idx) # Crops with improved stability according AMMI. print(Idx[order(Idx[,3]),]) # Crops with better response and improved stability according AMMI. print(Idx[order(Idx[,4]),])library(agricolae) # Index AMMI data(plrv) model<- with(plrv,AMMI(Locality, Genotype, Rep, Yield, console=FALSE)) Idx<-index.AMMI(model) names(Idx) # Crops with improved stability according AMMI. print(Idx[order(Idx[,3]),]) # Crops with better response and improved stability according AMMI. print(Idx[order(Idx[,4]),])
Scientists use a formula called the biodiversity index to describe the amount of species diversity in a given area.
index.bio(data, method = c("Margalef", "Simpson.Dom", "Simpson.Div", "Berger.Parker", "McIntosh", "Shannon"), level=95, nboot=100, console=TRUE)index.bio(data, method = c("Margalef", "Simpson.Dom", "Simpson.Div", "Berger.Parker", "McIntosh", "Shannon"), level=95, nboot=100, console=TRUE)
data |
number of specimens |
method |
Describe method bio-diversity |
level |
Significant level |
nboot |
size bootstrap |
console |
output console TRUE |
method bio-diversity. "Margalef" "Simpson.Dom" "Simpson.Div" "Berger.Parker" "McIntosh" "Shannon"
Index and confidence intervals.
Felipe de Mendiburu
Magurran, A.E. (1988) Ecological diversity and its measurement. Princeton University Press Efron, B., Tibshirani, R. (1993) An Introduction to the Boostrap. Chapman and Hall/CRC
library(agricolae) data(paracsho) # date 22-06-05 and treatment CON = application with insecticide specimens <- paracsho[1:10,6] output1 <- index.bio(specimens,method="Simpson.Div",level=95,nboot=100) output2 <- index.bio(specimens,method="Shannon",level=95,nboot=100) rbind(output1, output2)library(agricolae) data(paracsho) # date 22-06-05 and treatment CON = application with insecticide specimens <- paracsho[1:10,6] output1 <- index.bio(specimens,method="Simpson.Div",level=95,nboot=100) output2 <- index.bio(specimens,method="Shannon",level=95,nboot=100) rbind(output1, output2)
Smith's index of soil heterogeneity is used primarily to derive optimum plot size. The index gives a single value as a quantitative measure of soil heterogeneity in an area. Graph CV for plot size and shape.
index.smith(data, PLOT=TRUE,...)index.smith(data, PLOT=TRUE,...)
data |
dataframe or matrix |
PLOT |
graphics, TRUE or FALSE |
... |
Parameters of the plot() |
Vx=V(x)/x b
V(x) is the between-plot variance, Vx is the variance per unit area for plot size of x basic unit, and b is the Smith' index of soil heterogeneity.
model |
function pattern of uniformity |
uniformity |
Table of the soil uniformity |
Felipe de Mendiburu
Statistical Procedures for Agriculture Research. Second Edition. Kwanchai A. Gomez and Arturo A. Gomez. 1976. USA
library(agricolae) data(rice) #startgraph table<-index.smith(rice, main="Relationship between CV per unit area and plot size",col="red") #endgraph uniformity <- data.frame(table$uniformity) uniformity # regression variance per unit area an plot size. model <- lm(Vx ~ I(log(Size)),uniformity) coeff <- coef(model) x<-1:max(uniformity$Size) Vx<- coeff[1]+coeff[2]*log(x) #startgraph plot(x,Vx, type="l", col="blue", main="Relationship between variance per unit area and plot size") points(uniformity$Size,uniformity$Vx) #endgraphlibrary(agricolae) data(rice) #startgraph table<-index.smith(rice, main="Relationship between CV per unit area and plot size",col="red") #endgraph uniformity <- data.frame(table$uniformity) uniformity # regression variance per unit area an plot size. model <- lm(Vx ~ I(log(Size)),uniformity) coeff <- coef(model) x<-1:max(uniformity$Size) Vx<- coeff[1]+coeff[2]*log(x) #startgraph plot(x,Vx, type="l", col="blue", main="Relationship between variance per unit area and plot size") points(uniformity$Size,uniformity$Vx) #endgraph
List class intervals.
inter.freq(x)inter.freq(x)
x |
class graph.freq, histogram or numeric |
It show interval classes.
Felipe de Mendiburu
polygon.freq, table.freq, stat.freq,
graph.freq, sturges.freq, join.freq,
ogive.freq, normal.freq
library(agricolae) # example 1 data(growth) h<-hist(growth$height,plot=FALSE) inter.freq(h) # example 2 x<-seq(10,40,5) y<-c(2,6,8,7,3,4) inter.freq(x) histogram <- graph.freq(x,counts=y)library(agricolae) # example 1 data(growth) h<-hist(growth$height,plot=FALSE) inter.freq(h) # example 2 x<-seq(10,40,5) y<-c(2,6,8,7,3,4) inter.freq(x) histogram <- graph.freq(x,counts=y)
In many situations it is required to join classes because of the low .frequency in the intervals. In this process, it is required to join the intervals and ad the .frequencies of them.
join.freq(histogram, join)join.freq(histogram, join)
histogram |
Class graph.freq |
join |
vector |
New histogram with union of classes.
Felipe de Mendiburu
polygon.freq, table.freq, stat.freq,
inter.freq, sturges.freq, graph.freq,
ogive.freq, normal.freq
library(agricolae) data(natives) # histogram h1<-graph.freq(natives$size,plot=FALSE) round(table.freq(h1),4) # Join classes 9, 10,11 and 12 with little frequency. h2<-join.freq(h1,9:12) # new table plot(h2,col="bisque",xlab="Size") round(summary(h2),4)library(agricolae) data(natives) # histogram h1<-graph.freq(natives$size,plot=FALSE) round(table.freq(h1),4) # Join classes 9, 10,11 and 12 with little frequency. h2<-join.freq(h1,9:12) # new table plot(h2,col="bisque",xlab="Size") round(summary(h2),4)
Correlation of Kendall two set. Compute exact p-value with ties.
kendall(data1, data2)kendall(data1, data2)
data1 |
vector |
data2 |
vector |
The correlation of data1, data2 vector with the statistical value and its probability
Felipe de Mendiburu
Numerical Recipes in C. Second Edition. Pag 634
library(agricolae) x <-c(1,1,1,4,2,2,3,1,3,2,1,1,2,3,2,1,1,2,1,2) y <-c(1,1,2,3,4,4,2,1,2,3,1,1,3,4,2,1,1,3,1,2) kendall(x,y)library(agricolae) x <-c(1,1,1,4,2,2,3,1,3,2,1,1,2,3,2,1,1,2,1,2) y <-c(1,1,2,3,4,4,2,1,2,3,1,1,3,4,2,1,1,3,1,2) kendall(x,y)
It makes the multiple comparison with Kruskal-Wallis. The alpha parameter by default is 0.05. Post hoc test is using the criterium Fisher's least significant difference. The adjustment methods include the Bonferroni correction and others.
kruskal(y, trt, alpha = 0.05, p.adj=c("none","holm","hommel", "hochberg", "bonferroni", "BH", "BY", "fdr"), group=TRUE, main = NULL,console=FALSE)kruskal(y, trt, alpha = 0.05, p.adj=c("none","holm","hommel", "hochberg", "bonferroni", "BH", "BY", "fdr"), group=TRUE, main = NULL,console=FALSE)
y |
response |
trt |
treatment |
alpha |
level signification |
p.adj |
Method for adjusting p values (see p.adjust) |
group |
TRUE or FALSE |
main |
Title |
console |
logical, print output |
For equal or different repetition.
For the adjustment methods, see the function p.adjusted.
p-adj = "none" is t-student.
statistics |
Statistics of the model |
parameters |
Design parameters |
means |
Statistical summary of the study variable |
comparison |
Comparison between treatments |
groups |
Formation of treatment groups |
Felipe de Mendiburu
Practical Nonparametrics Statistics. W.J. Conover, 1999
BIB.test, DAU.test, duncan.test,
durbin.test, friedman, HSD.test,
LSD.test, Median.test, PBIB.test,
REGW.test, scheffe.test, SNK.test,
waerden.test, waller.test, plot.group
library(agricolae) data(corn) str(corn) comparison<-with(corn,kruskal(observation,method,group=TRUE, main="corn")) comparison<-with(corn,kruskal(observation,method,p.adj="bon",group=FALSE, main="corn"))library(agricolae) data(corn) str(corn) comparison<-with(corn,kruskal(observation,method,group=TRUE, main="corn")) comparison<-with(corn,kruskal(observation,method,p.adj="bon",group=FALSE, main="corn"))
It obtains the value of the kurtosis for a normally distributed variable. The result is similar to SAS.
kurtosis(x)kurtosis(x)
x |
a numeric vector |
x |
The kurtosis of x |
library(agricolae) x<-c(3,4,5,2,3,4,5,6,4,NA,7) kurtosis(x) # value is -0.1517996library(agricolae) x<-c(3,4,5,2,3,4,5,6,4,NA,7) kurtosis(x) # value is -0.1517996
A special function for the group of treatments in the multiple comparison tests. Use plot.group.
lastC(x)lastC(x)
x |
letters |
x |
Returns the last character of a string |
Felipe de Mendiburu
library(agricolae) x<-c("a","ab","b","c","cd") lastC(x) # "a" "b" "b" "c" "d"library(agricolae) x<-c("a","ab","b","c","cd") lastC(x) # "a" "b" "b" "c" "d"
LATEBLIGHT is a mathematical model that simulates the effect of weather, host growth and resistance, and fungicide use on asexual development and growth of Phytophthora infestans on potato foliage.
lateblight(WS, Cultivar, ApplSys,InocDate, LGR, IniSpor, SR, IE, LP, InMicCol, MatTime=c('EARLYSEASON','MIDSEASON','LATESEASON'),...)lateblight(WS, Cultivar, ApplSys,InocDate, LGR, IniSpor, SR, IE, LP, InMicCol, MatTime=c('EARLYSEASON','MIDSEASON','LATESEASON'),...)
WS |
object weather-severity |
Cultivar |
chr |
ApplSys |
chr |
InocDate |
days |
LGR |
num, see example |
IniSpor |
num |
SR |
num, see example |
IE |
num, Initialization infection |
LP |
num, latent period |
InMicCol |
num |
MatTime |
chr |
... |
plot graphics parameters |
LATEBLIGHT Version LB2004 was created in October 2004 (Andrade-Piedra et al., 2005a, b and c), based on the C-version written by B.E. Ticknor ('BET 21191 modification of cbm8d29.c'), reported by Doster et al. (1990) and described in detail by Fry et al. (1991) (This version is referred as LB1990 by Andrade-Piedra et al. [2005a]). The first version of LATEBLIGHT was developed by Bruhn and Fry (1981) and described in detail by Bruhn et al. (1980).
Ofile |
"Date","nday","MicCol","SimSeverity",... |
Gfile |
"dates","nday","MeanSeverity","StDevSeverity" |
All format data for date is yyyy-mm,dd, for example "2000-04-22". change with function as.Date()
Jorge L. Andrade-Piedra (1) ([email protected]), Gregory A. Forbes (1) ([email protected]), Robert J. Hijmans (2) ([email protected]), William E. Fry (3) ([email protected]) Translation from C language into SAS language: G.A. Forbes Modifications: J.L. Andrade-Piedra and R.J. Hijmans Translation from SAS into R: Felipe de Mendiburu (1) (1) International Potato Center, P.O. Box 1558, Lima 12, Peru (2) University of California, One Shields Avenue, Davis, California 95616, USA (3) Cornell University, 351 Plant Science, Ithaca, NY 14853, USA
Andrade-Piedra, J. L., Hijmans, R. J., Forbes, G. A., Fry, W. E.,
and Nelson, R. J. 2005a. Simulation of potato late blight in the Andes.
I: Modification and parameterization of the LATEBLIGHT model. Phytopathology 95:1191-1199.
Andrade-Piedra, J. L., Hijmans, R. J., Juarez, H. S., Forbes,
G. A., Shtienberg, D., and Fry, W. E. 2005b. Simulation of potato late blight
in the Andes. II: Validation of the LATEBLIGHT model. Phytopathology 95:1200-1208.
Andrade-Piedra, J. L., Forbes, G. A., Shtienberg, D., Grunwald, N. J.,
Chacon, M. G., Taipe, M. V., Hijmans, R. J., and Fry, W. E. 2005c.
Qualification of a plant disease simulation model: Performance of the LATEBLIGHT
model across a broad range of environments. Phytopathology 95:1412-1422.
Bruhn, J.A., Bruck, R.I., Fry, W.E., Arneson, P.A., and Keokosky, E.V. 1980.
User's manual for LATEBLIGHT: a plant disease management game. Cornell University,
Department of Plant Pathology, Ithaca, NY, USA. Mimeo 80-1.
Bruhn, J.A., and Fry, W.E. 1981. Analysis of potato late blight epidemiology
by simulation modeling. Phytopathology 71:612-616.
Doster, M. A., Milgroom, M. G., and Fry, W. E. 1990. Quantification of factors
influencing potato late blight suppression and selection for metalaxyl resistance
in Phytophthora infestans - A simulation approach. Phytopathology 80:1190-1198.
Fry, W.E., Milgroom, M.G., Doster, M.A., Bruhn, J.A., and Bruck, R.I. 1991.
LATEBLIGHT: a plant disease management game - User Manual. Version 3.1.
Microsoft Windows Adaptation by B. E. Ticknor, and P. A. Arneson. Ithaca,
Cornell University, Department of Plant Pathology, Ithaca, NY, USA.
library(agricolae) 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 Lateblight 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" #-------------------------- 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)") # reproduce graph x<- model$Ofile$nday y<- model$Ofile$SimSeverity w<- model$Gfile$nday z<- model$Gfile$MeanSeverity Min<-model$Gfile$MinObs Max<-model$Gfile$MaxObs plot(x,y,type="l",xlim=c(65,95),lwd=1.5,xlab="Time (days after emergence)", ylab="Severity (Percentage)") points(w,z,col="blue",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="blue") } legend("topleft",c("Disease progress curves","Weather-Severity"), title="Description",lty=1,pch=c(3,19),col=c("black","blue"))library(agricolae) 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 Lateblight 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" #-------------------------- 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)") # reproduce graph x<- model$Ofile$nday y<- model$Ofile$SimSeverity w<- model$Gfile$nday z<- model$Gfile$MeanSeverity Min<-model$Gfile$MinObs Max<-model$Gfile$MaxObs plot(x,y,type="l",xlim=c(65,95),lwd=1.5,xlab="Time (days after emergence)", ylab="Severity (Percentage)") points(w,z,col="blue",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="blue") } legend("topleft",c("Disease progress curves","Weather-Severity"), title="Description",lty=1,pch=c(3,19),col=c("black","blue"))
It makes the Line x Tester Genetic Analysis. It also estimates the general and specific combinatory ability effects and the line and tester genetic contribution.
lineXtester(replications, lines, testers, y)lineXtester(replications, lines, testers, y)
replications |
Replications |
lines |
Lines |
testers |
Testers |
y |
Variable, response |
ANOVA with parents and crosses
ANOVA for line X tester analysis
ANOVA for line X tester analysis including parents
GCA Effects: Lines Effects, Testers Effects and SCA Effects.
Standard Errors for Combining Ability Effects.
Genetic Components.
...
Proportional contribution of lines, testers and their interactions to total variance
return anova(formula = Y ~ Replications + Treatments).
where the Treatments contains parents, crosses and crosses vs Parents.
The crosses contains Lines, Testers and its interaction .
Felipe de Mendiburu
Biometrical Methods in Quantitative Genetic Analysis, Singh, Chaudhary. 1979. Hierarchial and factorial mating designs for quantitative genetic analysis in tetrasomic potato. R. Ortis; A.Golmirzaie. Theor Appl Genet (2002) 104:675-679
# see structure line by testers library(agricolae) # example 1 data(heterosis) site1<-subset(heterosis,heterosis[,1]==1) output1<-with(site1,lineXtester(Replication, Female, Male, v2)) # example 2 data(LxT) str(LxT) output2<-with(LxT,lineXtester(replication, line, tester, yield))# see structure line by testers library(agricolae) # example 1 data(heterosis) site1<-subset(heterosis,heterosis[,1]==1) output1<-with(site1,lineXtester(Replication, Female, Male, v2)) # example 2 data(LxT) str(LxT) output2<-with(LxT,lineXtester(replication, line, tester, yield))
Multiple comparisons of treatments by means of LSD and a grouping of treatments. The level by alpha default is 0.05. Returns p-values adjusted using one of several methods
LSD.test(y, trt, DFerror, MSerror, alpha = 0.05, p.adj=c("none","holm","hommel", "hochberg", "bonferroni", "BH", "BY", "fdr"), group=TRUE, main = NULL,console=FALSE)LSD.test(y, trt, DFerror, MSerror, alpha = 0.05, p.adj=c("none","holm","hommel", "hochberg", "bonferroni", "BH", "BY", "fdr"), group=TRUE, main = NULL,console=FALSE)
y |
model(aov or lm) or answer of the experimental unit |
trt |
Constant( only y=model) or vector treatment applied to each experimental unit |
DFerror |
Degrees of freedom of the experimental error |
MSerror |
Means square error of the experimental |
alpha |
Level of risk for the test |
p.adj |
Method for adjusting p values (see p.adjust) |
group |
TRUE or FALSE |
main |
title of the study |
console |
logical, print output |
For equal or different repetition.
For the adjustment methods, see the function p.adjusted.
p-adj ="none" is t-student.
It is necessary first makes a analysis of variance.
if model=y, then to apply the instruction:
LSD.test(model, "trt", alpha = 0.05, p.adj=c("none","holm","hommel",
"hochberg", "bonferroni", "BH", "BY", "fdr"), group=TRUE, main = NULL,console=FALSE)
where the model class is aov or lm.
statistics |
Statistics of the model |
parameters |
Design parameters |
means |
Statistical summary of the study variable |
comparison |
Comparison between treatments |
groups |
Formation of treatment groups |
Felipe de Mendiburu
Steel, R.; Torri,J; Dickey, D.(1997) Principles and Procedures of Statistics A Biometrical Approach. pp178.
BIB.test, DAU.test, duncan.test,
durbin.test, friedman, HSD.test,
kruskal, Median.test, PBIB.test,
REGW.test, scheffe.test, SNK.test,
waerden.test, waller.test, plot.group
library(agricolae) data(sweetpotato) model<-aov(yield~virus, data=sweetpotato) out <- LSD.test(model,"virus", p.adj="bonferroni") #stargraph # Variation range: max and min plot(out) #endgraph # Old version LSD.test() df<-df.residual(model) MSerror<-deviance(model)/df out <- with(sweetpotato,LSD.test(yield,virus,df,MSerror)) #stargraph # Variation interquartil range: Q75 and Q25 plot(out,variation="IQR") #endgraph out<-LSD.test(model,"virus",p.adj="hommel",console=TRUE) plot(out,variation="SD") # variation standard deviationlibrary(agricolae) data(sweetpotato) model<-aov(yield~virus, data=sweetpotato) out <- LSD.test(model,"virus", p.adj="bonferroni") #stargraph # Variation range: max and min plot(out) #endgraph # Old version LSD.test() df<-df.residual(model) MSerror<-deviance(model)/df out <- with(sweetpotato,LSD.test(yield,virus,df,MSerror)) #stargraph # Variation interquartil range: Q75 and Q25 plot(out,variation="IQR") #endgraph out<-LSD.test(model,"virus",p.adj="hommel",console=TRUE) plot(out,variation="SD") # variation standard deviation
Data frame with yield by line x tester.
data(LxT)data(LxT)
A data frame with 92 observations on the following 4 variables.
a numeric vector
a numeric vector
a numeric vector
a numeric vector
Biometrical Methods in Quantitative Genetic Analysis, Singh, Chaudhary. 1979
A partial study on 27 molecular markers.
data(markers)data(markers)
A data frame with 23 observations on the following 27 variables.
marker1a numeric vector
marker2a numeric vector
marker3a numeric vector
marker4a numeric vector
marker5a numeric vector
marker6a numeric vector
marker7a numeric vector
marker8a numeric vector
marker9a numeric vector
marker10a numeric vector
marker11a numeric vector
marker12a numeric vector
marker13a numeric vector
marker14a numeric vector
marker15a numeric vector
marker16a numeric vector
marker17a numeric vector
marker18a numeric vector
marker19a numeric vector
marker20a numeric vector
marker21a numeric vector
marker22a numeric vector
marker23a numeric vector
marker24a numeric vector
marker25a numeric vector
marker26a numeric vector
marker27a numeric vector
International Potato Center Lima-Peru.
International Potato Center Lima-Peru.
library(agricolae) data(markers) str(markers)library(agricolae) data(markers) str(markers)
A nonparametric test for several independent samples. The median test is designed to examine whether several samples came from populations having the same median.
Median.test(y,trt,alpha=0.05,correct=TRUE,simulate.p.value = FALSE, group = TRUE, main = NULL,console=TRUE)Median.test(y,trt,alpha=0.05,correct=TRUE,simulate.p.value = FALSE, group = TRUE, main = NULL,console=TRUE)
y |
Variable response |
trt |
Treatments |
alpha |
error type I |
correct |
a logical indicating whether to apply continuity correction when computing the test statistic for 2 groups. The correction will not be bigger than the differences themselves. No correction is done if simulate.p.value = TRUE. |
simulate.p.value |
a logical indicating whether to compute p-values by Monte Carlo simulation |
group |
TRUE or FALSE |
main |
Title |
console |
logical, print output |
The data consist of k samples of possibly unequal sample size.
Greater: is the number of values that exceed the median of all data and
LessEqual: is the number less than or equal to the median of all data.
statistics |
Statistics of the model |
parameters |
Design parameters |
medians |
Statistical summary of the study variable |
comparison |
Comparison between treatments |
groups |
Formation of treatment groups |
Felipe de Mendiburu
Practical Nonparametrics Statistics. W.J. Conover, 1999
BIB.test, DAU.test, duncan.test,
durbin.test, friedman, HSD.test,
kruskal, LSD.test, PBIB.test,
REGW.test, scheffe.test, SNK.test,
waerden.test, waller.test, plot.group
library(agricolae) # example 1 data(corn) out<-with(corn,Median.test(observation,method,console=FALSE)) z<-bar.err(out$medians,variation = "range",ylim=c(0,120), space=2,border=4,col=3,density=10,angle=45) # example 2 out<-with(corn,Median.test(observation,method,console=FALSE,group=FALSE)) print(out$comparison)library(agricolae) # example 1 data(corn) out<-with(corn,Median.test(observation,method,console=FALSE)) z<-bar.err(out$medians,variation = "range",ylim=c(0,120), space=2,border=4,col=3,density=10,angle=45) # example 2 out<-with(corn,Median.test(observation,method,console=FALSE,group=FALSE)) print(out$comparison)
An irrigation system evaluation by exudation using four varieties of melon, under modality of sowing, SIMPLE ROW. The goal is to analyze the behavior of three hybrid melon varieties and one standard.
data(melon)data(melon)
A data frame with 16 observations on the following 4 variables.
rowa numeric vector
cola numeric vector
varietya factor with levels V1 V2 V3 V4
yielda numeric vector
Varieties: Hibrido Mission (V1); Hibrido Mark (V2); Hibrido Topfligth (V3); Hibrido Hales Best Jumbo (V4).
Tesis. "Evaluacion del sistema de riego por exudacion utilizando cuatro variedades de melon, bajo modalidad de siembra, SIMPLE HILERA". Alberto Angeles L. Universidad Agraria la Molina - Lima Peru.
Universidad Nacional Agraria la molina.
library(agricolae) data(melon) str(melon)library(agricolae) data(melon) str(melon)
Random generation form data, use function density and parameters
montecarlo(data, k, ...)montecarlo(data, k, ...)
data |
vector or object(hist, graph.freq) |
k |
number of simulations |
... |
Other parameters of the function density, only if data is vector |
Generate random numbers with empirical distribution.
Felipe de Mendiburu
library(agricolae) r<-rnorm(50, 10,2) montecarlo(r, k=100, kernel="epanechnikov") # other example h<-hist(r,plot=FALSE) montecarlo(h, k=100) # other example breaks<-c(0, 150, 200, 250, 300) counts<-c(10, 20, 40, 30) op<-par(mfrow=c(1,2),cex=0.8,mar=c(2,3,0,0)) h1<-graph.freq(x=breaks,counts=counts,plot=FALSE) r<-montecarlo(h, k=1000) plot(h1,frequency = 3,ylim=c(0,0.008)) text(90,0.006,"Population\n100 obs.") h2<-graph.freq(r,breaks,frequency = 3,ylim=c(0,0.008)) lines(density(r),col="blue") text(90,0.006,"Montecarlo\n1000 obs.") par(op)library(agricolae) r<-rnorm(50, 10,2) montecarlo(r, k=100, kernel="epanechnikov") # other example h<-hist(r,plot=FALSE) montecarlo(h, k=100) # other example breaks<-c(0, 150, 200, 250, 300) counts<-c(10, 20, 40, 30) op<-par(mfrow=c(1,2),cex=0.8,mar=c(2,3,0,0)) h1<-graph.freq(x=breaks,counts=counts,plot=FALSE) r<-montecarlo(h, k=1000) plot(h1,frequency = 3,ylim=c(0,0.008)) text(90,0.006,"Population\n100 obs.") h2<-graph.freq(r,breaks,frequency = 3,ylim=c(0,0.008)) lines(density(r),col="blue") text(90,0.006,"Montecarlo\n1000 obs.") par(op)
An evaluation of the number, weight and size of 24 native potatoes varieties.
data(natives)data(natives)
A data frame with 876 observations on the following 4 variables.
varietya numeric vector
numbera numeric vector
weighta numeric vector
sizea numeric vector
International Potato Center. CIP - Lima Peru.
library(agricolae) data(natives) str(natives)library(agricolae) data(natives) str(natives)
The resistance for the transformable nonadditivity, due to J. W. Tukey, is based on the detection of a curvilinear relation between y-est(y) and est(y). A freedom degree for the transformable nonadditivity.
nonadditivity(y, factor1, factor2, df, MSerror)nonadditivity(y, factor1, factor2, df, MSerror)
y |
Answer of the experimental unit |
factor1 |
Firts treatment applied to each experimental unit |
factor2 |
Second treatment applied to each experimental unit |
df |
Degrees of freedom of the experimental error |
MSerror |
Means square error of the experimental |
Only two factor: Block and treatment or factor 1 and factor 2.
P, Q and non-additivity analysis of variance
Felipe de Mendiburu
1. Steel, R.; Torri,J; Dickey, D.(1997) Principles and Procedures of Statistics A Biometrical Approach
2. George E.P. Box; J. Stuart Hunter and William G. Hunter. Statistics for experimenters. Wile Series in probability and statistics
library(agricolae) 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))library(agricolae) 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))
A normal distribution graph elaborated from the histogram previously constructed. The average and variance are obtained from the data grouped in the histogram.
normal.freq(histogram, frequency=1, ...)normal.freq(histogram, frequency=1, ...)
histogram |
object constructed by the function hist |
frequency |
1=counts, 2=relative, 3=density |
... |
Other parameters of the function hist |
Felipe de Mendiburu
polygon.freq, table.freq, stat.freq,
inter.freq, sturges.freq, join.freq,
ogive.freq, graph.freq
library(agricolae) data(growth) #startgraph h1<-with(growth,hist(height,col="green",xlim=c(6,14))) normal.freq(h1,col="blue") #endgraph #startgraph h2<-with(growth,graph.freq(height,col="yellow",xlim=c(6,14),frequency=2)) normal.freq(h2,frequency=2) #endgraphlibrary(agricolae) data(growth) #startgraph h1<-with(growth,hist(height,col="green",xlim=c(6,14))) normal.freq(h1,col="blue") #endgraph #startgraph h2<-with(growth,graph.freq(height,col="yellow",xlim=c(6,14),frequency=2)) normal.freq(h2,frequency=2) #endgraph
It plots the cumulative relative .frequencies with the intervals of classes defined in the histogram.
ogive.freq(histogram,type="",xlab="",ylab="",axes="",las=1,...)ogive.freq(histogram,type="",xlab="",ylab="",axes="",las=1,...)
histogram |
object created by the function hist() or graph.freq() |
type |
what type of plot should be drawn. See plot() |
xlab |
x labels |
ylab |
y labels |
axes |
TRUE or FALSE |
las |
values 0,1,2 and 3 are the axis styles. see plot() |
... |
Parameters of the plot() |
Ogive points.
Felipe de Mendiburu
polygon.freq, table.freq, stat.freq,
inter.freq, sturges.freq, join.freq,
graph.freq, normal.freq
library(agricolae) data(growth) h<-graph.freq(growth$height,plot=FALSE) points<-ogive.freq(h,col="red",frame=FALSE, xlab="Height", ylab="Accumulated relative frequency", main="ogive") plot(points,type="b",pch=16,las=1,bty="l")library(agricolae) data(growth) h<-graph.freq(growth$height,plot=FALSE) points<-ogive.freq(h,col="red",frame=FALSE, xlab="Height", ylab="Accumulated relative frequency", main="ogive") plot(points,type="b",pch=16,las=1,bty="l")
This function allows us to compare the treatments averages or the adding of their ranges with the minimal significant difference which can vary from one comparison to another one.
order.group(trt, means, N, MSerror, Tprob, std.err, parameter=1, snk=0, DFerror=NULL,alpha=NULL,sdtdif=NULL,vartau=NULL,console)order.group(trt, means, N, MSerror, Tprob, std.err, parameter=1, snk=0, DFerror=NULL,alpha=NULL,sdtdif=NULL,vartau=NULL,console)
trt |
Treatments |
means |
Means of treatment |
N |
Replications |
MSerror |
Mean square error |
Tprob |
minimum value for the comparison |
std.err |
standard error |
parameter |
Constante 1 (Sd), 0.5 (Sx) |
snk |
Constante = 1 (Student Newman Keuls) |
DFerror |
Degrees of freedom of the experimental error |
alpha |
Level of risk for the test |
sdtdif |
standard deviation of difference in BIB |
vartau |
matrix var-cov in PBIB |
console |
logical, print output |
This function was changed by orderPvalue function that use agricolae. Now the grouping in agricolae is with the probability of the treatments differences and alpha level.
The output is data frame.
trt |
Treatment Levels, Factor |
means |
height, Numeric |
M |
groups levels, Factor |
N |
replications, Numeric |
std.err |
Standard error, Numeric |
It is considered 81 labels as maximum for the formation of groups, greater number will not have label.
Felipe de Mendiburu
library(agricolae) treatments <- c("A","B","C","D","E","F") means<-c(20,40,35,72,49,58) std.err<-c(1.2, 2, 1.5, 2.4, 1, 3.1) replications <- c(4,4,3,4,3,3) MSerror <- 55.8 value.t <- 2.1314 groups<-order.group(treatments,means,replications,MSerror,value.t,std.err,console=FALSE) print(groups)library(agricolae) treatments <- c("A","B","C","D","E","F") means<-c(20,40,35,72,49,58) std.err<-c(1.2, 2, 1.5, 2.4, 1, 3.1) replications <- c(4,4,3,4,3,3) MSerror <- 55.8 value.t <- 2.1314 groups<-order.group(treatments,means,replications,MSerror,value.t,std.err,console=FALSE) print(groups)
When there are treatments and their respective values, these can be compared with a minimal difference of meaning.
orderPvalue(treatment, means, alpha, pvalue, console)orderPvalue(treatment, means, alpha, pvalue, console)
treatment |
treatment |
means |
means of treatment |
alpha |
Alpha value, significante value to comparison |
pvalue |
Matrix of probabilities to comparison |
console |
logical, print output |
The means and groups for treatments.
It is considered 81 labels as maximum for the formation of groups, greater number will not have label.
Felipe de Mendiburu
library(agricolae) treatments <- c("A","B","C") means<-c(2,5,3) alpha <- 0.05 pvalue<-matrix(1,nrow=3,ncol=3) pvalue[1,2]<-pvalue[2,1]<-0.03 pvalue[1,3]<-pvalue[3,1]<-0.10 pvalue[2,3]<-pvalue[3,2]<-0.06 out<-orderPvalue(treatments,means,alpha,pvalue,console=TRUE) barplot(out[,1],names.arg = row.names(out),col=colors()[84:87]) legend("topright",as.character(out$groups),pch=15,col=colors()[84:87],box.col=0)library(agricolae) treatments <- c("A","B","C") means<-c(2,5,3) alpha <- 0.05 pvalue<-matrix(1,nrow=3,ncol=3) pvalue[1,2]<-pvalue[2,1]<-0.03 pvalue[1,3]<-pvalue[3,1]<-0.10 pvalue[2,3]<-pvalue[3,2]<-0.06 out<-orderPvalue(treatments,means,alpha,pvalue,console=TRUE) barplot(out[,1],names.arg = row.names(out),col=colors()[84:87]) legend("topright",as.character(out$groups),pch=15,col=colors()[84:87],box.col=0)
Potato Wild
data(pamCIP)data(pamCIP)
A data frame with 43 observations on the following 107 variables. Rownames: code and genotype's name. column data: molecular markers.
To study the molecular markers in Wild.
Laboratory data.
International Potato Center Lima-Peru (CIP)
library(agricolae) data(pamCIP) str(pamCIP)library(agricolae) data(pamCIP) str(pamCIP)
A locality in Peru. A biodiversity.
data(paracsho)data(paracsho)
A data frame with 110 observations on the following 6 variables.
datea factor with levels 15-12-05 17-11-05 18-10-05 20-09-05 22-06-05 23-08-05 28-07-05
plota factor with levels PARACSHO
Treatmenta factor with levels CON SIN
Ordena factor with levels COLEOPTERA DIPTERA HEMIPTERA HYMENOPTERA LEPIDOPTERA NEUROPTERA NEUROPTERO NOCTUIDAE
Familya factor with levels AGROMYZIDAE ANTHOCORIDAE ANTHOMYIIDAE ANTHOMYLIDAE BLEPHAROCERIDAE BRACONIDAE BROCONIDAE CALUPHORIDAE CECIDOMYIDAE CHENEUMONIDAE CHNEUMONIDAE CHRYOMELIDAE CICADELLIDAE CULICIDAE ERIOCPAMIDAE HEMEROBIIDAE ICHNEUMONIDAE LOUCHAPIDAE MIRIDAE MUSCIDAE MUSICADAE MUSLIDAE MYCETOPHILIDAE MYCETOPHILIIDAE NENPHALIDAE NOCLUIDAE NOCTERIDAE NOCTUIDAE PERALIDAE PIPUNCULIDAE PROCTOTRUPIDAE PSYLLIDAE PYRALIDAE SARCOPHAGIDAE SARCOPILAGIDAE SCATOPHAGIDAE SCATOPHOGIDAE SCIARIDAE SERSIDAE SYRPHIDAE TACHINIDAE TIPULIDAE
Number.of.specimensa numeric vector
Country Peru, Deparment Junin, province Tarma, locality Huasahuasi.
Entomology dataset.
International Potato Center.
library(agricolae) data(paracsho) str(paracsho)library(agricolae) data(paracsho) str(paracsho)
If the cause and effect relationship is well defined, it is possible to represent the whole system of variables in a diagram form known as path-analysis. The function calculates the direct and indirect effects and uses the variables correlation or covariance.
path.analysis(corr.x, corr.y)path.analysis(corr.x, corr.y)
corr.x |
Matrix of correlations of the independent variables |
corr.y |
vector of dependent correlations with each one of the independent ones |
It is necessary first to calculate the correlations.
Direct and indirect effects and residual Effect^2.
Felipe de Mendiburu
Biometrical Methods in Quantitative Genetic Analysis, Singh, Chaudhary. 1979
# Path analysis. Multivarial Analysis. Anderson. Prentice Hall, pag 616 library(agricolae) # Example 1 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") path.analysis(corr.x,corr.y) # Example 2 # data of the progress of the disease related bacterial wilt to the ground # for the component CE Ca K2 Cu data(wilt) data(soil) x<-soil[,c(3,12,14,20)] y<-wilt[,14] cor.y<-correlation(y,x)$correlation cor.x<-correlation(x)$correlation path.analysis(cor.x,cor.y)# Path analysis. Multivarial Analysis. Anderson. Prentice Hall, pag 616 library(agricolae) # Example 1 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") path.analysis(corr.x,corr.y) # Example 2 # data of the progress of the disease related bacterial wilt to the ground # for the component CE Ca K2 Cu data(wilt) data(soil) x<-soil[,c(3,12,14,20)] y<-wilt[,14] cor.y<-correlation(y,x)$correlation cor.x<-correlation(x)$correlation path.analysis(cor.x,cor.y)
Analysis of variance PBIB and comparison mean adjusted. Applied to resoluble designs: Lattices and alpha design.
PBIB.test(block,trt,replication,y,k, method=c("REML","ML","VC"), test = c("lsd","tukey"), alpha=0.05, console=FALSE, group=TRUE)PBIB.test(block,trt,replication,y,k, method=c("REML","ML","VC"), test = c("lsd","tukey"), alpha=0.05, console=FALSE, group=TRUE)
block |
blocks |
trt |
Treatment |
replication |
Replication |
y |
Response |
k |
Block size |
method |
Estimation method: REML, ML and VC |
test |
Comparison treatments |
alpha |
Significant test |
console |
logical, print output |
group |
logical, groups |
Method of comparison treatment.
lsd: least significant difference.
tukey: Honestly significant difference.
Estimate: specifies the estimation method for the covariance parameters.
The REML is the default method. The REML specification performs residual (restricted) maximum likelihood, and
The ML specification performs maximum likelihood, and
the VC specifications apply only to variance component models.
The PBIB.test() function can be called inside a function (improvement by Nelson Nazzicari, Ph.D. Bioinformatician)
ANOVA |
Analysis of variance |
method |
Estimation method: REML, ML and VC |
parameters |
Design parameters |
statistics |
Statistics of the model |
model |
Object: estimation model |
Fstat |
Criterion AIC and BIC |
comparison |
Comparison between treatments |
means |
Statistical summary of the study variable |
groups |
Formation of treatment groups |
vartau |
Variance-Covariance Matrix |
F. de Mendiburu
1. Iterative Analysis of Generalizad Lattice Designs. E.R. Williams (1977) Austral J. Statistics 19(1) 39-42.
2. Experimental design. Cochran and Cox. Second edition. Wiley Classics Library Edition published 1992
BIB.test, DAU.test, duncan.test,
durbin.test, friedman, HSD.test,
kruskal, LSD.test, Median.test,
REGW.test, scheffe.test, SNK.test,
waerden.test, waller.test, plot.group
require(agricolae) # alpha design Genotype<-c(paste("gen0",1:9,sep=""),paste("gen",10:30,sep="")) ntr<-length(Genotype) r<-2 k<-3 s<-10 obs<-ntr*r b <- s*r book<-design.alpha(Genotype,k,r,seed=5) book$book[,3]<- gl(20,3) dbook<-book$book # dataset 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) rm(Genotype) # not run # analysis # require(nlme) # method = REML or LM in PBIB.test and require(MASS) method=VC model <- with(dbook,PBIB.test(block, Genotype, replication, yield, k=3, method="VC")) # model$ANOVA # plot(model,las=2)require(agricolae) # alpha design Genotype<-c(paste("gen0",1:9,sep=""),paste("gen",10:30,sep="")) ntr<-length(Genotype) r<-2 k<-3 s<-10 obs<-ntr*r b <- s*r book<-design.alpha(Genotype,k,r,seed=5) book$book[,3]<- gl(20,3) dbook<-book$book # dataset 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) rm(Genotype) # not run # analysis # require(nlme) # method = REML or LM in PBIB.test and require(MASS) method=VC model <- with(dbook,PBIB.test(block, Genotype, replication, yield, k=3, method="VC")) # model$ANOVA # plot(model,las=2)
Biplot AMMI.
## S3 method for class 'AMMI' plot(x,first=1,second=2,third=3,number=FALSE,gcol=NULL,ecol=NULL, angle=25,lwd=1.8,length=0.1,xlab=NULL,ylab=NULL,xlim=NULL,ylim=NULL,...)## S3 method for class 'AMMI' plot(x,first=1,second=2,third=3,number=FALSE,gcol=NULL,ecol=NULL, angle=25,lwd=1.8,length=0.1,xlab=NULL,ylab=NULL,xlim=NULL,ylim=NULL,...)
x |
object AMMI |
first |
position axis x, 0=Y-dependent, 1=PC1, 2=PC2, 3=PC3 |
second |
position axis y,0=Y-dependent, 1=PC1, 2=PC2, 3=PC3 |
third |
position axis z,0=Y-dependent, 1=PC1, 2=PC2, 3=PC3 |
number |
TRUE or FALSE names or number genotypes |
gcol |
genotype color |
ecol |
environment color |
angle |
angle from the shaft of the arrow to the edge of the arrow head |
lwd |
parameter line width in function arrow |
length |
parameter length in function arrow |
xlab |
x labels |
ylab |
y labels |
xlim |
x limites |
ylim |
y limites |
... |
other parameters of plot |
Produce graphs biplot.
Felipe de Mendiburu
library(agricolae) data(plrv) model<- with(plrv,AMMI(Locality, Genotype, Rep, Yield)) # biplot PC2 vs PC1 plot(model) ## plot PC1 vs Yield plot(model,0,1,gcol="blue",ecol="green")library(agricolae) data(plrv) model<- with(plrv,AMMI(Locality, Genotype, Rep, Yield)) # biplot PC2 vs PC1 plot(model) ## plot PC1 vs Yield plot(model,0,1,gcol="blue",ecol="green")
In many situations it has intervals of class defined with its respective frequencies. By means of this function, the graphic of frequency is obtained and it is possible to superpose the normal distribution, polygon of frequency, Ojiva and to construct the table of complete frequency.
## S3 method for class 'graph.freq' plot(x, breaks=NULL,counts=NULL,frequency=1,plot=TRUE, nclass=NULL,xlab="",ylab="",axes = "",las=1,...)## S3 method for class 'graph.freq' plot(x, breaks=NULL,counts=NULL,frequency=1,plot=TRUE, nclass=NULL,xlab="",ylab="",axes = "",las=1,...)
x |
a vector of values, a object hist(), graphFreq() |
counts |
frequency and x is class intervals |
breaks |
a vector giving the breakpoints between histogram cells |
frequency |
1=counts, 2=relative, 3=density |
plot |
logic |
nclass |
number of classes |
xlab |
x labels |
ylab |
y labels |
axes |
TRUE or FALSE |
las |
values 0,1,2 and 3 are the axis styles. see plot() |
... |
other parameters of plot |
breaks |
a vector giving the breakpoints between histogram cells |
counts |
frequency and x is class intervals |
mids |
center point in class |
relative |
Relative frequency, height |
density |
Density frequency, height |
Felipe de Mendiburu
polygon.freq, table.freq,
stat.freq,inter.freq,sturges.freq,
join.freq,ogive.freq, normal.freq
library(agricolae) data(genxenv) yield <- subset(genxenv$YLD,genxenv$ENV==2) yield <- round(yield,1) h<- graph.freq(yield,axes=FALSE, frequency=1, ylab="frequency",col="yellow") axis(1,h$breaks) axis(2,seq(0,20,0.1)) # To reproduce histogram. h1 <- plot(h, col="blue", frequency=2,border="red", density=8,axes=FALSE, xlab="YIELD",ylab="relative") axis(1,h$breaks) axis(2,seq(0,.4,0.1)) # summary, only frecuency limits <-seq(10,40,5) frequencies <-c(2,6,8,7,3,4) #startgraph h<-graph.freq(limits,counts=frequencies,col="bisque",xlab="Classes") polygon.freq(h,col="red") title( main="Histogram and polygon of frequency", ylab=".frequency") #endgraph # Statistics measures<-stat.freq(h) print(measures) # frequency table full round(table.freq(h),2) #startgraph # ogive ogive.freq(h,col="red",type="b",ylab="Accumulated relative frequency", xlab="Variable") # only frequency polygon h<-graph.freq(limits,counts=frequencies,border=FALSE,col=NULL,xlab=" ",ylab="") title( main="Polygon of frequency", xlab="Variable", ylab="Frecuency") polygon.freq(h,col="blue") grid(col="brown") #endgraph # Draw curve for Histogram h<- graph.freq(yield,axes=FALSE, frequency=3, ylab="f(yield)",col="yellow") axis(1,h$breaks) axis(2,seq(0,0.18,0.03),las=2) lines(density(yield), col = "red", lwd = 2) title("Draw curve for Histogram")library(agricolae) data(genxenv) yield <- subset(genxenv$YLD,genxenv$ENV==2) yield <- round(yield,1) h<- graph.freq(yield,axes=FALSE, frequency=1, ylab="frequency",col="yellow") axis(1,h$breaks) axis(2,seq(0,20,0.1)) # To reproduce histogram. h1 <- plot(h, col="blue", frequency=2,border="red", density=8,axes=FALSE, xlab="YIELD",ylab="relative") axis(1,h$breaks) axis(2,seq(0,.4,0.1)) # summary, only frecuency limits <-seq(10,40,5) frequencies <-c(2,6,8,7,3,4) #startgraph h<-graph.freq(limits,counts=frequencies,col="bisque",xlab="Classes") polygon.freq(h,col="red") title( main="Histogram and polygon of frequency", ylab=".frequency") #endgraph # Statistics measures<-stat.freq(h) print(measures) # frequency table full round(table.freq(h),2) #startgraph # ogive ogive.freq(h,col="red",type="b",ylab="Accumulated relative frequency", xlab="Variable") # only frequency polygon h<-graph.freq(limits,counts=frequencies,border=FALSE,col=NULL,xlab=" ",ylab="") title( main="Polygon of frequency", xlab="Variable", ylab="Frecuency") polygon.freq(h,col="blue") grid(col="brown") #endgraph # Draw curve for Histogram h<- graph.freq(yield,axes=FALSE, frequency=3, ylab="f(yield)",col="yellow") axis(1,h$breaks) axis(2,seq(0,0.18,0.03),las=2) lines(density(yield), col = "red", lwd = 2) title("Draw curve for Histogram")
It plots bars of the averages of treatments to compare. It uses the objects generated by a procedure of comparison like LSD, HSD, Kruskall, Waller-Duncan, Friedman or Durbin. It can also display the 'average' value over each bar in a bar chart.
## S3 method for class 'group' plot(x,variation=c("range","IQR","SE","SD"), decreasing = TRUE, horiz=FALSE,col=NULL,xlim=NULL,ylim=NULL,main=NULL,cex=NULL,hy=0,...)## S3 method for class 'group' plot(x,variation=c("range","IQR","SE","SD"), decreasing = TRUE, horiz=FALSE,col=NULL,xlim=NULL,ylim=NULL,main=NULL,cex=NULL,hy=0,...)
x |
Object created by a test of comparison |
variation |
in lines by range, IQR, standard deviation or error |
decreasing |
Logical, decreasing order of the mean |
horiz |
Horizontal or vertical image |
col |
line colors |
xlim |
optional, axis x limits |
ylim |
optional, axis y limits |
main |
optional, main title |
cex |
optional, group label size |
hy |
optional, default =0, sum group label position |
... |
Parameters of the function barplot() |
The output is a vector that indicates the position of the treatments on the coordinate axes.
Felipe de Mendiburu
BIB.test, DAU.test, duncan.test,
durbin.test, friedman, HSD.test,
kruskal, LSD.test, Median.test,
PBIB.test, REGW.test, scheffe.test,
SNK.test, waerden.test, waller.test
library(agricolae) data(sweetpotato) model<-aov(yield~virus,data=sweetpotato) comparison<- LSD.test(model,"virus",alpha=0.01,group=TRUE) #startgraph op<-par(cex=1.5) plot(comparison,horiz=TRUE,xlim=c(0,50),las=1) title(cex.main=0.8,main="Comparison between\ntreatment means",xlab="Yield",ylab="Virus") #endgraph par(op)library(agricolae) data(sweetpotato) model<-aov(yield~virus,data=sweetpotato) comparison<- LSD.test(model,"virus",alpha=0.01,group=TRUE) #startgraph op<-par(cex=1.5) plot(comparison,horiz=TRUE,xlim=c(0,50),las=1) title(cex.main=0.8,main="Comparison between\ntreatment means",xlab="Yield",ylab="Virus") #endgraph par(op)
Experimental data in blocks, factor A in plots and factor B in sub-plots.
data(plots)data(plots)
A data frame with 18 observations on the following 5 variables.
blocka numeric vector
plota factor with levels p1 p2 p3 p4 p5 p6
Aa factor with levels a1 a2
Ba factor with levels b1 b2 b3
yielda numeric vector
International Potato Center. CIP
library(agricolae) data(plots) str(plots) plots[,1] <-as.factor(plots[,1]) # split-plot analysis model <- aov(yield ~ block + A + Error(plot)+ B + A:B, data=plots) summary(model) b<-nlevels(plots$B) a<-nlevels(plots$A) r<-nlevels(plots$block) dfa <- df.residual(model$plot) Ea <-deviance(model$plot)/dfa dfb <- df.residual(model$Within) Eb <-deviance(model$Within)/dfb Eab <- (Ea +(b-1)*Eb)/(b*r) # Satterthwaite dfab<-(Ea +(b-1)*Eb)^2/(Ea^2/dfa +((b-1)*Eb)^2/dfb) # Comparison A, A(b1), A(b2), A(b3) comparison1 <-with(plots,LSD.test(yield,A,dfa,Ea)) comparison2 <-with(plots,LSD.test(yield[B=="b1"],A[B=="b1"],dfab,Eab)) comparison3 <-with(plots,LSD.test(yield[B=="b2"],A[B=="b2"],dfab,Eab)) comparison4 <-with(plots,LSD.test(yield[B=="b3"],A[B=="b3"],dfab,Eab)) # Comparison B, B(a1), B(a2) comparison5 <-with(plots,LSD.test(yield,B,dfb,Eb)) comparison6 <-with(plots,LSD.test(yield[A=="a1"],B[A=="a1"],dfb,Eb)) comparison7 <-with(plots,LSD.test(yield[A=="a2"],B[A=="a2"],dfb,Eb))library(agricolae) data(plots) str(plots) plots[,1] <-as.factor(plots[,1]) # split-plot analysis model <- aov(yield ~ block + A + Error(plot)+ B + A:B, data=plots) summary(model) b<-nlevels(plots$B) a<-nlevels(plots$A) r<-nlevels(plots$block) dfa <- df.residual(model$plot) Ea <-deviance(model$plot)/dfa dfb <- df.residual(model$Within) Eb <-deviance(model$Within)/dfb Eab <- (Ea +(b-1)*Eb)/(b*r) # Satterthwaite dfab<-(Ea +(b-1)*Eb)^2/(Ea^2/dfa +((b-1)*Eb)^2/dfb) # Comparison A, A(b1), A(b2), A(b3) comparison1 <-with(plots,LSD.test(yield,A,dfa,Ea)) comparison2 <-with(plots,LSD.test(yield[B=="b1"],A[B=="b1"],dfab,Eab)) comparison3 <-with(plots,LSD.test(yield[B=="b2"],A[B=="b2"],dfab,Eab)) comparison4 <-with(plots,LSD.test(yield[B=="b3"],A[B=="b3"],dfab,Eab)) # Comparison B, B(a1), B(a2) comparison5 <-with(plots,LSD.test(yield,B,dfb,Eb)) comparison6 <-with(plots,LSD.test(yield[A=="a1"],B[A=="a1"],dfb,Eb)) comparison7 <-with(plots,LSD.test(yield[A=="a2"],B[A=="a2"],dfb,Eb))
Six environments: Ayacucho, La Molina 02, San Ramon 02, Huancayo, La Molina 03, San Ramon 03.
data(plrv)data(plrv)
A data frame with 504 observations on the following 6 variables.
Genotypea factor with levels 102.18 104.22 121.31 141.28 157.26 163.9 221.19 233.11 235.6 241.2 255.7 314.12 317.6 319.20 320.16 342.15 346.2 351.26 364.21 402.7 405.2 406.12 427.7 450.3 506.2 Canchan Desiree Unica
Localitya factor with levels Ayac Hyo-02 LM-02 LM-03 SR-02 SR-03
Repa numeric vector
WeightPlanta numeric vector
WeightPlota numeric vector
Yielda numeric vector
International Potato Center Lima-Peru
International Potato Center Lima-Peru
library(agricolae) data(plrv) str(plrv)library(agricolae) data(plrv) str(plrv)
The polygon is constructed single or on a histogram. It is necessary to execute the function previously hist.
polygon.freq(histogram, frequency=1, ...)polygon.freq(histogram, frequency=1, ...)
histogram |
Object constructed by the function hist |
frequency |
numeric, counts(1), relative(2) and density(3) |
... |
Other parameters of the function hist |
Felipe de Mendiburu Delgado
polygon.freq, table.freq, stat.freq,
inter.freq, sturges.freq, join.freq,
graph.freq, normal.freq
library(agricolae) data(growth) #startgraph h1<-with(growth,hist(height,border=FALSE,xlim=c(6,14))) polygon.freq(h1,frequency=1,col="red") #endgraph #startgraph h2<-with(growth,graph.freq(height,frequency=2,col="yellow",xlim=c(6,14))) polygon.freq(h2,frequency=2,col="red") #endgraphlibrary(agricolae) data(growth) #startgraph h1<-with(growth,hist(height,border=FALSE,xlim=c(6,14))) polygon.freq(h1,frequency=1,col="red") #endgraph #startgraph h2<-with(growth,graph.freq(height,frequency=2,col="yellow",xlim=c(6,14))) polygon.freq(h2,frequency=2,col="red") #endgraph
A study on the yield of two potato varieties performed at the CIP experimental station.
data(potato)data(potato)
A data frame with 18 observations on the following 4 variables.
datea numeric vector
varietya factor with levels Canchan Unica
harvesta numeric vector
cuttinga numeric vector
Experimental data.
International Potato Center
library(agricolae) data(potato) str(potato)library(agricolae) data(potato) str(potato)
The assessment of the population of R.solanacearum on the floor took place after 48 hours of infestation, during days 15, 29, 43, 58, and 133 days after the infestation soil. More information on soil data(soil).
data(ralstonia)data(ralstonia)
A data frame with 13 observations on the following 8 variables.
placea factor with levels Chmar Chz Cnt1 Cnt2 Cnt3 Hco1 Hco2 Hco3 Hyo1 Hyo2 Namora SR1 SR2
Day2a numeric vector
Day15a numeric vector
Day29a numeric vector
Day43a numeric vector
Day58a numeric vector
Day73a numeric vector
Day133a numeric vector
Logarithm average counts of colonies on plates containing half of M-SMSA 3 repetitions (3 plates by repetition) incubated at 30 degrees centigrade for 48 hours. log(1+UFC/g soil)
Experimental field, 2004. Data Kindly provided by Dr. Sylvie Priou, Liliam Gutarra and Pedro Aley.
International Potato Center. CIP - Lima Peru.
library(agricolae) data(ralstonia) str(ralstonia)library(agricolae) data(ralstonia) str(ralstonia)
It makes the regressions homogeneity test for a group of treatments where each observation presents a linearly dependent reply from another one. There is a linear function in every treatment. The objective is to find out if the linear models of each treatment come from the same population.
reg.homog(trt, x, y)reg.homog(trt, x, y)
trt |
treatment |
x |
independent variable |
y |
dependent variable |
list objects:
Number regressions.
Residual.
Difference of regression.
DF.homgeneity (homogenity degree free).
DF.Residual (degree free error).
F.value. Test statitics.
P.value. P Value (Significant
Criterion. conclusion
Felipe de Mendiburu
Book in Spanish: Metodos estadisticos para la investigacion. Calzada Benza 1960
library(agricolae) data(frijol) evaluation<-with(frijol,reg.homog(technology,index,production)) # Example 2. Applied Regression Analysis a Research tools # 1988. John O.Rawlings. Wadsworth & brooks/cole Advanced Books # & Software. Pacific Grove. Califonia. # Statistics/probability. Series LineNumber<-c(rep("39","30"),rep("52","30")) PlantingDate<-rep(c("16","20","21"),20) HeadWt <- c(2.5,3.0,2.2,2.2,2.8,1.8,3.1,2.8,1.6,4.3,2.7,2.1,2.5,2.6,3.3,4.3, 2.8,3.8,3.8,2.6,3.2,4.3,2.6,3.6,1.7,2.6,4.2,3.1,3.5,1.6,2.0,4.0,1.5,2.4,2.8, 1.4,1.9,3.1,1.7,2.8,4.2,1.3,1.7,3.7,1.7,3.2,3.0,1.6,2.0,2.2,1.4,2.2,2.3,1.0, 2.2,3.8,1.5,2.2,2.0,1.6) Ascorbic <-c(51,65,54,55,52,59,45,41,66,42,51,54,53,41,45,50,45,49,50,51,49, 52,45,55,56,61,49,49,42,68,58,52,78,55,70,75,67,57,70,61,58,84,67,47,71,68, 56,72,58,72,62,63,63,68,56,54,66,72,60,72) trt<-paste(LineNumber,PlantingDate,sep="-") output<-reg.homog(trt,HeadWt,Ascorbic)library(agricolae) data(frijol) evaluation<-with(frijol,reg.homog(technology,index,production)) # Example 2. Applied Regression Analysis a Research tools # 1988. John O.Rawlings. Wadsworth & brooks/cole Advanced Books # & Software. Pacific Grove. Califonia. # Statistics/probability. Series LineNumber<-c(rep("39","30"),rep("52","30")) PlantingDate<-rep(c("16","20","21"),20) HeadWt <- c(2.5,3.0,2.2,2.2,2.8,1.8,3.1,2.8,1.6,4.3,2.7,2.1,2.5,2.6,3.3,4.3, 2.8,3.8,3.8,2.6,3.2,4.3,2.6,3.6,1.7,2.6,4.2,3.1,3.5,1.6,2.0,4.0,1.5,2.4,2.8, 1.4,1.9,3.1,1.7,2.8,4.2,1.3,1.7,3.7,1.7,3.2,3.0,1.6,2.0,2.2,1.4,2.2,2.3,1.0, 2.2,3.8,1.5,2.2,2.0,1.6) Ascorbic <-c(51,65,54,55,52,59,45,41,66,42,51,54,53,41,45,50,45,49,50,51,49, 52,45,55,56,61,49,49,42,68,58,52,78,55,70,75,67,57,70,61,58,84,67,47,71,68, 56,72,58,72,62,63,63,68,56,54,66,72,60,72) trt<-paste(LineNumber,PlantingDate,sep="-") output<-reg.homog(trt,HeadWt,Ascorbic)
Multiple range tests for all pairwise comparisons, to obtain a confident inequalities multiple range tests.
REGW.test(y, trt, DFerror, MSerror, alpha = 0.05, group=TRUE, main = NULL,console=FALSE)REGW.test(y, trt, DFerror, MSerror, alpha = 0.05, group=TRUE, main = NULL,console=FALSE)
y |
model(aov or lm) or answer of the experimental unit |
trt |
Constant( only y=model) or vector treatment applied to each experimental unit |
DFerror |
Degree free |
MSerror |
Mean Square Error |
alpha |
Significant level |
group |
TRUE or FALSE |
main |
Title |
console |
logical, print output |
It is necessary first makes a analysis of variance.
if y = model, then to apply the instruction:
REGW.test (model, "trt", alpha = 0.05, group = TRUE, main = NULL, console = FALSE)
where the model class is aov or lm.
statistics |
Statistics of the model |
parameters |
Design parameters |
regw |
Critical Range Table |
means |
Statistical summary of the study variable |
comparison |
Comparison between treatments |
groups |
Formation of treatment groups |
Felipe de Mendiburu
Multiple comparisons theory and methods. Departament of statistics the Ohio State University. USA, 1996. Jason C. Hsu. Chapman Hall/CRC
BIB.test, DAU.test, duncan.test,
durbin.test, friedman, HSD.test,
kruskal, LSD.test, Median.test,
PBIB.test, scheffe.test, SNK.test,
waerden.test, waller.test, plot.group
library(agricolae) data(sweetpotato) model<-aov(yield~virus,data=sweetpotato) out<- REGW.test(model,"virus", main="Yield of sweetpotato. Dealt with different virus") print(out) REGW.test(model,"virus",alpha=0.05,console=TRUE,group=FALSE)library(agricolae) data(sweetpotato) model<-aov(yield~virus,data=sweetpotato) out<- REGW.test(model,"virus", main="Yield of sweetpotato. Dealt with different virus") print(out) REGW.test(model,"virus",alpha=0.05,console=TRUE,group=FALSE)
This process finds the curve of CV for a different number of markers which allows us to determine the number of optimal markers for a given relative variability. A method of the curvature.
resampling.cv(A, size, npoints)resampling.cv(A, size, npoints)
A |
data frame or matrix of binary data |
size |
number of re-samplings |
npoints |
Number of points to consider the model |
lm(formula = CV ~ I(1/marker))
Table with variation coefficient by number of markers
Felipe de Mendiburu
Efron, B., Tibshirani, R. (1993) An Introduction to the Boostrap. Chapman and Hall/CRC
library(agricolae) #example table of molecular markers data(markers) study<-resampling.cv(markers,size=1,npoints=15) # # Results of the model summary(study$model) coef<-coef(study$model) py<-predict(study$model) Rsq<-summary(study$model)$"r.squared" table.cv <- data.frame(study$table.cv,estimate=py) print(table.cv) # Plot CV #startgraph limy<-max(table.cv[,2])+10 plot(table.cv[,c(1,2)],col="red",frame=FALSE,xlab="number of markers", ylim=c(10,limy), ylab="CV",cex.main=0.8,main="Estimation of the number of markers") ty<-quantile(table.cv[,2],1) tx<-median(table.cv[,1]) tz<-quantile(table.cv[,2],0.95) text(tx,ty, cex=0.8,as.expression(substitute(CV == a + frac(b,markers), list(a=round(coef[1],2),b=round(coef[2],2)))) ) text(tx,tz,cex=0.8,as.expression(substitute(R^2==r,list(r=round(Rsq,3))))) # Plot CV = a + b/n.markers fy<-function(x,a,b) a+b/x x<-seq(2,max(table.cv[,1]),length=50) y <- coef[1] + coef[2]/x lines(x,y,col="blue") #grid(col="brown") rug(table.cv[,1]) #endgraphlibrary(agricolae) #example table of molecular markers data(markers) study<-resampling.cv(markers,size=1,npoints=15) # # Results of the model summary(study$model) coef<-coef(study$model) py<-predict(study$model) Rsq<-summary(study$model)$"r.squared" table.cv <- data.frame(study$table.cv,estimate=py) print(table.cv) # Plot CV #startgraph limy<-max(table.cv[,2])+10 plot(table.cv[,c(1,2)],col="red",frame=FALSE,xlab="number of markers", ylim=c(10,limy), ylab="CV",cex.main=0.8,main="Estimation of the number of markers") ty<-quantile(table.cv[,2],1) tx<-median(table.cv[,1]) tz<-quantile(table.cv[,2],0.95) text(tx,ty, cex=0.8,as.expression(substitute(CV == a + frac(b,markers), list(a=round(coef[1],2),b=round(coef[2],2)))) ) text(tx,tz,cex=0.8,as.expression(substitute(R^2==r,list(r=round(Rsq,3))))) # Plot CV = a + b/n.markers fy<-function(x,a,b) a+b/x x<-seq(2,max(table.cv[,1]),length=50) y <- coef[1] + coef[2]/x lines(x,y,col="blue") #grid(col="brown") rug(table.cv[,1]) #endgraph
This process consists of finding the values of P-value by means of a re-sampling (permutation) process along with the values obtained by variance analysis.
resampling.model(model,data,k,console=FALSE)resampling.model(model,data,k,console=FALSE)
model |
model in R |
data |
data for the study of the model |
k |
number of re-samplings |
console |
logical, print output |
Model solution with resampling.
Felipe de Mendiburu
Efron, B., Tibshirani, R. (1993) An Introduction to the Boostrap. Chapman and Hall/CRC Phillip I. Good, (2001) Resampling Methods. Birkhauser. Boston . Basel . Berlin
#example 1 Simple linear regression library(agricolae) data(clay) model<-"ralstonia ~ days" analysis<-resampling.model(model,clay,k=2,console=TRUE) #example 2 Analysis of variance: RCD data(sweetpotato) model<-"yield~virus" analysis<-resampling.model(model,sweetpotato,k=2,console=TRUE) #example 3 Simple linear regression data(Glycoalkaloids) model<-"HPLC ~ spectrophotometer" analysis<-resampling.model(model,Glycoalkaloids,k=2,console=TRUE) #example 4 Factorial in RCD 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=2,console=TRUE)#example 1 Simple linear regression library(agricolae) data(clay) model<-"ralstonia ~ days" analysis<-resampling.model(model,clay,k=2,console=TRUE) #example 2 Analysis of variance: RCD data(sweetpotato) model<-"yield~virus" analysis<-resampling.model(model,sweetpotato,k=2,console=TRUE) #example 3 Simple linear regression data(Glycoalkaloids) model<-"HPLC ~ spectrophotometer" analysis<-resampling.model(model,Glycoalkaloids,k=2,console=TRUE) #example 4 Factorial in RCD 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=2,console=TRUE)
The data correspond to the yield of rice variety IR8 (g/m2) for land uniformity studies. The growing area is 18x36 meters.
data(rice)data(rice)
A data frame with 36 observations on the following 18 variables.
V1a numeric vector
V2a numeric vector
V3a numeric vector
V4a numeric vector
V5a numeric vector
V6a numeric vector
V7a numeric vector
V8a numeric vector
V9a numeric vector
V10a numeric vector
V11a numeric vector
V12a numeric vector
V13a numeric vector
V14a numeric vector
V15a numeric vector
V16a numeric vector
V17a numeric vector
V18a numeric vector
Table 12.1 Measuring Soil Heterogeneity
Statistical Procedures for Agriculture Research. Second Edition. Kwanchai A. Gomez and Arturo A. Gomez. 1976. USA Pag. 481.
Statistical Procedures for Agriculture Research. Second Edition. Kwanchai A. Gomez and Arturo A. Gomez. 1976. USA
library(agricolae) data(rice) str(rice)library(agricolae) data(rice) str(rice)
Mother/Baby Trials allow farmers and researchers to test best-bet technologies or new cultivars. Evaluation of advanced Clones of potato in the Valley of Rio Chillon - PERU (2004)
data(RioChillon)data(RioChillon)
The format is list of 2:
1. mother: data.frame: 30 obs. of 3 variables:
- block (3 levels)
- clon (10 levels)
- yield (kg.)
2. babies: data.frame: 90 obs. of 3 variables:
- farmer (9 levels)
- clon (10 levels)
- yield (kg.)
1. Replicated researcher-managed "mother trials" with typically 10 treatments
evaluated in small plots.
2. Unreplicated "baby trials" with 10 treatments evaluated in large plots.
3. The "baby trials" with a subset of the treatments in the mother trial.
Experimental field.
International Potato Center. CIP - Lima Peru.
# Analisys the Mother/Baby Trial Design library(agricolae) data(RioChillon) # First analysis the Mother Trial Design. model<-aov(yield ~ block + clon, RioChillon$mother) anova(model) cv.model(model) comparison<-with(RioChillon$mother,LSD.test(yield,clon, 18, 4.922, group=TRUE)) # Second analysis the babies Trial. comparison<-with(RioChillon$babies,friedman(farmer,clon, yield, group=TRUE)) # Third # The researcher makes use of data from both mother and baby trials and thereby obtains # information on suitability of new technologies or cultivars # for different agro-ecologies and acceptability to farmers.# Analisys the Mother/Baby Trial Design library(agricolae) data(RioChillon) # First analysis the Mother Trial Design. model<-aov(yield ~ block + clon, RioChillon$mother) anova(model) cv.model(model) comparison<-with(RioChillon$mother,LSD.test(yield,clon, 18, 4.922, group=TRUE)) # Second analysis the babies Trial. comparison<-with(RioChillon$babies,friedman(farmer,clon, yield, group=TRUE)) # Third # The researcher makes use of data from both mother and baby trials and thereby obtains # information on suitability of new technologies or cultivars # for different agro-ecologies and acceptability to farmers.
Scheffe 1959, method is very general in that all possible contrasts can be tested for significance and confidence intervals can be constructed for the corresponding linear. The test is conservative.
scheffe.test(y, trt, DFerror, MSerror, Fc, alpha = 0.05, group=TRUE, main = NULL, console=FALSE )scheffe.test(y, trt, DFerror, MSerror, Fc, alpha = 0.05, group=TRUE, main = NULL, console=FALSE )
y |
model(aov or lm) or answer of the experimental unit |
trt |
Constant( only y=model) or vector treatment applied to each experimental unit |
DFerror |
Degrees of freedom |
MSerror |
Mean Square Error |
Fc |
F Value |
alpha |
Significant level |
group |
TRUE or FALSE |
main |
Title |
console |
logical, print output |
It is necessary first makes a analysis of variance.
if y = model, then to apply the instruction:
scheffe.test (model, "trt", alpha = 0.05, group = TRUE, main = NULL, console = FALSE)
where the model class is aov or lm.
statistics |
Statistics of the model |
parameters |
Design parameters |
means |
Statistical summary of the study variable |
comparison |
Comparison between treatments |
groups |
Formation of treatment groups |
Felipe de Mendiburu
Robert O. Kuehl. 2nd ed. Design of experiments. Duxbury, copyright 2000.
Steel, R.; Torri,J; Dickey, D.(1997) Principles and Procedures of Statistics
A Biometrical Approach. pp189
BIB.test, DAU.test, duncan.test,
durbin.test, friedman, HSD.test,
kruskal, LSD.test, Median.test,
PBIB.test, REGW.test, SNK.test,
waerden.test, waller.test, plot.group
library(agricolae) data(sweetpotato) model<-aov(yield~virus, data=sweetpotato) comparison <- scheffe.test(model,"virus", group=TRUE,console=TRUE, main="Yield of sweetpotato\nDealt with different virus") # Old version scheffe.test() df<-df.residual(model) MSerror<-deviance(model)/df Fc<-anova(model)["virus",4] out <- with(sweetpotato,scheffe.test(yield, virus, df, MSerror, Fc)) print(out)library(agricolae) data(sweetpotato) model<-aov(yield~virus, data=sweetpotato) comparison <- scheffe.test(model,"virus", group=TRUE,console=TRUE, main="Yield of sweetpotato\nDealt with different virus") # Old version scheffe.test() df<-df.residual(model) MSerror<-deviance(model)/df Fc<-anova(model)["virus",4] out <- with(sweetpotato,scheffe.test(yield, virus, df, MSerror, Fc)) print(out)
It finds the similarity matrix of binary tables (1 and 0).
similarity(A)similarity(A)
A |
Matrix, data binary |
Distance matrix. Class = dist.
Felipe de Mendiburu
#example table of molecular markers library(agricolae) data(markers) distance<-similarity(markers) #startgraph tree<-hclust(distance,method="mcquitty") plot(tree,col="blue") #endgraph#example table of molecular markers library(agricolae) data(markers) distance<-similarity(markers) #startgraph tree<-hclust(distance,method="mcquitty") plot(tree,col="blue") #endgraph
This process consists of validating the variance analysis results using a simulation process of the experiment. The validation consists of comparing the calculated values of each source of variation of the simulated data with respect to the calculated values of the original data. If in more than 50 percent of the cases they are higher than the real one, then it is considered favorable and the probability reported by the ANOVA is accepted, since the P-Value is the probability of (F > F.value).
simulation.model(model,file, categorical = NULL,k,console=FALSE)simulation.model(model,file, categorical = NULL,k,console=FALSE)
model |
Model in R |
file |
Data for the study of the model |
categorical |
position of the columns of the data that correspond to categorical variables |
k |
Number of simulations |
console |
logical, print output |
model |
ouput linear model, lm |
simulation |
anova simulation |
Felipe de Mendiburu
library(agricolae) #example 1 data(clay) model<-"ralstonia ~ days" simulation.model(model,clay,k=15,console=TRUE) #example 2 data(sweetpotato) model<-"yield~virus" simulation.model(model,sweetpotato,categorical=1,k=15,console=TRUE) #example 3 data(Glycoalkaloids) model<-"HPLC ~ spectrophotometer" simulation.model(model,Glycoalkaloids,k=15,console=TRUE) #example 4 data(potato) model<-"cutting~date+variety" simulation.model(model,potato,categorical=c(1,2,3),k=15,console=TRUE)library(agricolae) #example 1 data(clay) model<-"ralstonia ~ days" simulation.model(model,clay,k=15,console=TRUE) #example 2 data(sweetpotato) model<-"yield~virus" simulation.model(model,sweetpotato,categorical=1,k=15,console=TRUE) #example 3 data(Glycoalkaloids) model<-"HPLC ~ spectrophotometer" simulation.model(model,Glycoalkaloids,k=15,console=TRUE) #example 4 data(potato) model<-"cutting~date+variety" simulation.model(model,potato,categorical=c(1,2,3),k=15,console=TRUE)
Data frame for AMMI analysis with 50 genotypes in 5 environments.
data(sinRepAmmi)data(sinRepAmmi)
A data frame with 250 observations on the following 3 variables.
ENVa factor with levels A1 A2 A3 A4 A5
GENa numeric vector
YLDa numeric vector
Experimental data.
International Potato Center - Lima Peru.
library(agricolae) data(sinRepAmmi) str(sinRepAmmi)library(agricolae) data(sinRepAmmi) str(sinRepAmmi)
It returns the skewness of a distribution. It is similar to SAS.
skewness(x)skewness(x)
x |
a numeric vector |
The skewness of x.
library(agricolae) x<-c(3,4,5,2,3,4,NA,5,6,4,7) skewness(x) # value is 0,3595431, is slightly asimetrica (positive) to the rightlibrary(agricolae) x<-c(3,4,5,2,3,4,NA,5,6,4,7) skewness(x) # value is 0,3595431, is slightly asimetrica (positive) to the right
SNK is derived from Tukey, but it is less conservative (finds more differences). Tukey controls the error for all comparisons, where SNK only controls for comparisons under consideration. The level by alpha default is 0.05.
SNK.test(y, trt, DFerror, MSerror, alpha = 0.05, group=TRUE, main = NULL,console=FALSE)SNK.test(y, trt, DFerror, MSerror, alpha = 0.05, group=TRUE, main = NULL,console=FALSE)
y |
model(aov or lm) or answer of the experimental unit |
trt |
Constant( only y=model) or vector treatment applied to each experimental unit |
DFerror |
Degree free |
MSerror |
Mean Square Error |
alpha |
Significant level |
group |
TRUE or FALSE |
main |
Title |
console |
logical, print output |
It is necessary first makes a analysis of variance.
if y = model, then to apply the instruction:
SNK.test (model, "trt", alpha = 0.05, group = TRUE, main = NULL, console = FALSE)
where the model class is aov or lm.
statistics |
Statistics of the model |
parameters |
Design parameters |
snk |
Critical Range Table |
means |
Statistical summary of the study variable |
comparison |
Comparison between treatments |
groups |
Formation of treatment groups |
Felipe de Mendiburu
1. Principles and procedures of statistics a biometrical approach
Steel & Torry & Dickey. Third Edition 1997
2. Multiple comparisons theory and methods. Departament of statistics
the Ohio State University. USA, 1996. Jason C. Hsu. Chapman Hall/CRC.
BIB.test, DAU.test, duncan.test,
durbin.test, friedman, HSD.test,
kruskal, LSD.test, Median.test,
PBIB.test, REGW.test, scheffe.test,
waerden.test, waller.test, plot.group
library(agricolae) data(sweetpotato) model<-aov(yield~virus,data=sweetpotato) out <- SNK.test(model,"virus", console=TRUE, main="Yield of sweetpotato. Dealt with different virus") print(SNK.test(model,"virus", group=FALSE)) # version old SNK.test() df<-df.residual(model) MSerror<-deviance(model)/df out <- with(sweetpotato,SNK.test(yield,virus,df,MSerror, group=TRUE)) print(out$groups)library(agricolae) data(sweetpotato) model<-aov(yield~virus,data=sweetpotato) out <- SNK.test(model,"virus", console=TRUE, main="Yield of sweetpotato. Dealt with different virus") print(SNK.test(model,"virus", group=FALSE)) # version old SNK.test() df<-df.residual(model) MSerror<-deviance(model)/df out <- with(sweetpotato,SNK.test(yield,virus,df,MSerror, group=TRUE)) print(out$groups)
We analyzed the physical and chemical properties of different soils, as full characterization of soil and special analysis of micro-elements. These analyses were conducted in the laboratory analysis of soils, plants, water and fertilizers in the La Molina National Agrarian University (UNALM). To which the different soil samples were dried to the environment, screened (mesh 0.5xo, 5 mm) and sterilized by steam 4 to 5 hours with a Lindinger Steam aerator SA150 and SA700, with the possible aim of eliminating bacteria saprophytic or antagonists to prevent the growth of bacteria (R.solanacearum).
data(soil)data(soil)
A data frame with 13 observations on the following 23 variables.
placea factor with levels Chmar Chz Cnt1 Cnt2 Cnt3 Hco1 Hco2 Hco3 Hyo1 Hyo2 Namora SR1 SR2
pHa numeric vector
ECa numeric vector, electrical conductivity
CaCO3a numeric vector
MOa numeric vector
CICa numeric vector
Pa numeric vector
Ka numeric vector
sanda numeric vector
slimea numeric vector
claya numeric vector
Caa numeric vector
Mga numeric vector
K2a numeric vector
Naa numeric vector
Al_Ha numeric vector
K_Mga numeric vector
Ca_Mga numeric vector
Ba numeric vector
Cua numeric vector
Fea numeric vector
Mna numeric vector
Zna numeric vector
Cnt1= Canete, Cnt2=Valle Dulce(Canete), Cnt3=Valle Grande(Canete), Chz=Obraje-Carhuaz(Ancash), Chmar=Chucmar-Chota(Huanuco, Hco1= Mayobamba-Chinchao(Huanuco), Hco2=Nueva Independencia-Chinchao(Huanuco), Hco3=San Marcos-Umari(Huanuco), Hyo1=La Victoria-Huancayo(Junin), Hyo1=El Tambo-Huancayo(Junin), Namora=Namora(Cajamarca), SR1= El Milagro-San Ramon(Junin), Sr2=La Chinchana-San Ramon(Junin).
Experimental field, 2004. Data Kindly provided by Dr. Sylvie Priou, Liliam Gutarra and Pedro Aley.
International Potato Center - Lima, PERU.
library(agricolae) data(soil) str(soil)library(agricolae) data(soil) str(soil)
The variance analysis of a split plot design is divided into two parts: the plot-factor analysis and the sub-plot factor analysis.
sp.plot(block, pplot, splot, Y)sp.plot(block, pplot, splot, Y)
block |
replications |
pplot |
main-plot Factor |
splot |
sub-plot Factor |
Y |
Variable, response |
The split-plot design is specifically suited for a two-factor experiment on of the factors is assigned to main plot (main-plot factor), the second factor, called the subplot factor, is assigned into subplots. The model is mixed, the blocks are random and the study factors are fixed applied according to the design.
ANOVA: Splip plot analysis
Felipe de Mendiburu
Statistical procedures for agricultural research. Kwanchai A. Gomez, Arturo A. Gomez. Second Edition. 1984.
ssp.plot, strip.plot, design.split,
design.strip
library(agricolae) data(plots) model<-with(plots,sp.plot(block,A,B,yield)) # with aov plots[,1]<-as.factor(plots[,1]) AOV <- aov(yield ~ block + A*B + Error(block/A),data=plots) summary(AOV)library(agricolae) data(plots) model<-with(plots,sp.plot(block,A,B,yield)) # with aov plots[,1]<-as.factor(plots[,1]) AOV <- aov(yield ~ block + A*B + Error(block/A),data=plots) summary(AOV)
The variance analysis of a split-split plot design is divided into three parts: the main-plot, subplot and sub-subplot analysis.
ssp.plot(block, pplot, splot, ssplot, Y)ssp.plot(block, pplot, splot, ssplot, Y)
block |
replications |
pplot |
Factor main plot |
splot |
Factor subplot |
ssplot |
Factor sub-subplot |
Y |
Variable, response |
The split-split-plot design is an extension of the split-plot design to accommodate a third factor: one factor in main-plot, other in subplot and the third factor in sub-subplot. The model is mixed, the blocks are random and the study factors are fixed applied according to the design.
ANOVA: Splip Split plot analysis
Felipe de Mendiburu
Statistical procedures for agricultural research. Kwanchai A. Gomez, Arturo A. Gomez. Second Edition. 1984.
sp.plot, strip.plot, design.split,
design.strip
# Statistical procedures for agricultural research, pag 143 # Grain Yields of Three Rice Varieties Grown under #Three Management practices and Five Nitrogen levels; in a #split-split-plot design with nitrogen as main-plot, #management practice as subplot, and variety as sub-subplot #factores, with three replications. library(agricolae) f <- system.file("external/ssp.csv", package="agricolae") ssp<-read.csv(f) model<-with(ssp,ssp.plot(block,nitrogen,management,variety,yield)) gla<-model$gl.a; glb<-model$gl.b; glc<-model$gl.c Ea<-model$Ea; Eb<-model$Eb; Ec<-model$Ec op<-par(mfrow=c(1,3),cex=0.6) out1<-with(ssp,LSD.test(yield,nitrogen,gla,Ea,console=TRUE)) out2<-with(ssp,LSD.test(yield,management,glb,Eb,console=TRUE)) out3<-with(ssp,LSD.test(yield,variety,glc,Ec,console=TRUE)) plot(out1,xlab="Nitrogen",las=1,variation="IQR") plot(out2,xlab="Management",variation="IQR") plot(out3,xlab="Variety",variation="IQR") # with aov ssp$block<-factor(ssp$block) ssp$nitrogen<-factor(ssp$nitrogen) ssp$management<-factor(ssp$management) ssp$variety<-factor(ssp$variety) AOV<-aov(yield ~ block + nitrogen*management*variety + Error(block/nitrogen/management),data=ssp) summary(AOV) par(op)# Statistical procedures for agricultural research, pag 143 # Grain Yields of Three Rice Varieties Grown under #Three Management practices and Five Nitrogen levels; in a #split-split-plot design with nitrogen as main-plot, #management practice as subplot, and variety as sub-subplot #factores, with three replications. library(agricolae) f <- system.file("external/ssp.csv", package="agricolae") ssp<-read.csv(f) model<-with(ssp,ssp.plot(block,nitrogen,management,variety,yield)) gla<-model$gl.a; glb<-model$gl.b; glc<-model$gl.c Ea<-model$Ea; Eb<-model$Eb; Ec<-model$Ec op<-par(mfrow=c(1,3),cex=0.6) out1<-with(ssp,LSD.test(yield,nitrogen,gla,Ea,console=TRUE)) out2<-with(ssp,LSD.test(yield,management,glb,Eb,console=TRUE)) out3<-with(ssp,LSD.test(yield,variety,glc,Ec,console=TRUE)) plot(out1,xlab="Nitrogen",las=1,variation="IQR") plot(out2,xlab="Management",variation="IQR") plot(out3,xlab="Variety",variation="IQR") # with aov ssp$block<-factor(ssp$block) ssp$nitrogen<-factor(ssp$nitrogen) ssp$management<-factor(ssp$management) ssp$variety<-factor(ssp$variety) AOV<-aov(yield ~ block + nitrogen*management*variety + Error(block/nitrogen/management),data=ssp) summary(AOV) par(op)
A method based on the statistical ranges of the study variable per environment for the stability analysis.
stability.nonpar(data, variable = NULL, ranking = FALSE, console=FALSE)stability.nonpar(data, variable = NULL, ranking = FALSE, console=FALSE)
data |
First column the genotypes following environment |
variable |
Name of variable |
ranking |
logical, print ranking |
console |
logical, print output |
ranking |
data frame |
statistics |
Statistical analysis chi square test |
Felipe de Mendiburu
Haynes K G, Lambert D H, Christ B J, Weingartner D P, Douches D S, Backlund J E, Fry W and Stevenson W. 1998. Phenotypic stability of resistance to late blight in potato clones evaluated at eight sites in the United States American Journal Potato Research 75, pag 211-217.
library(agricolae) data(haynes) stability.nonpar(haynes,"AUDPC",ranking=TRUE,console=TRUE) # Example 2 data(CIC) data1<-CIC$comas[,c(1,6,7,17,18)] data2<-CIC$oxapampa[,c(1,6,7,19,20)] cic <- rbind(data1,data2) means <- by(cic[,5], cic[,c(2,1)], function(x) mean(x,na.rm=TRUE)) means <-as.data.frame(means[,]) cic.mean<-data.frame(genotype=row.names(means),means) cic.mean<-delete.na(cic.mean,"greater") out<-stability.nonpar(cic.mean) out$ranking out$statisticslibrary(agricolae) data(haynes) stability.nonpar(haynes,"AUDPC",ranking=TRUE,console=TRUE) # Example 2 data(CIC) data1<-CIC$comas[,c(1,6,7,17,18)] data2<-CIC$oxapampa[,c(1,6,7,19,20)] cic <- rbind(data1,data2) means <- by(cic[,5], cic[,c(2,1)], function(x) mean(x,na.rm=TRUE)) means <-as.data.frame(means[,]) cic.mean<-data.frame(genotype=row.names(means),means) cic.mean<-delete.na(cic.mean,"greater") out<-stability.nonpar(cic.mean) out$ranking out$statistics
This procedure calculates the stability variations as well as the statistics of selection for the yield and the stability. The averages of the genotype through the different environment repetitions are required for the calculations. The mean square error must be calculated from the joint variance analysis.
stability.par(data,rep,MSerror,alpha=0.1,main=NULL,cova = FALSE,name.cov=NULL, file.cov=0,console=FALSE)stability.par(data,rep,MSerror,alpha=0.1,main=NULL,cova = FALSE,name.cov=NULL, file.cov=0,console=FALSE)
data |
matrix of averages, by rows the genotypes and columns the environment |
rep |
Number of repetitions |
MSerror |
Mean Square Error |
alpha |
Label significant |
main |
Title |
cova |
Covariable |
name.cov |
Name covariable |
file.cov |
Data covariable |
console |
logical, print output |
Stable (i) determines the contribution of each genotype to GE interaction by calculating var(i); (ii) assigns ranks to genotypes from highest to lowest yield receiving the rank of 1; (iii) calculates protected LSD for mean yield comparisons; (iv) adjusts yield rank according to LSD (the adjusted rank labeled Y); (v) determines significance of var(i) usign an aproximate F-test; (vi) assigns stability rating (S) as follows: -8, -4 and -2 for var(i) significant at the 0.01, 0.05 and 0.10 probability levels, and 0 for nonsignificant var(i) ( the higher the var(i), the less stable the genotype); (vii) sums adjusted yield rank, Y, and stability rating, S, for each genotype to determine YS(i) statistic; and (viii) calculates mean YS(i) and identifies genotypes (selection) with YS(i) > mean YS(i).
analysis |
Analysis of variance |
statistics |
Statistics of the model |
stability |
summary stability analysis |
Felipe de Mendiburu
Kang, M. S. 1993. Simultaneous selection for yield and stability: Consequences for growers. Agron. J. 85:754-757. Manjit S. Kang and Robert Mangari. 1995. Stable: A basic program for calculating stability and yield-stability statistics. Agron. J. 87:276-277
library(agricolae) # example 1 # Experimental data, # replication rep= 4 # Mean square error, MSerror = 1.8 # 12 environment # 17 genotype = 1,2,3,.., 17 # yield averages of 13 genotypes in localities 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) #example 2 covariable. precipitation precipitation<- c(1000,1100,1200,1300,1400,1500,1600,1700,1800,1900,2000,2100) stability.par(dataStb, rep=4, MSerror=1.8, alpha=0.1, main="Genotype", cova=TRUE, name.cov="Precipitation", file.cov=precipitation,console=TRUE)library(agricolae) # example 1 # Experimental data, # replication rep= 4 # Mean square error, MSerror = 1.8 # 12 environment # 17 genotype = 1,2,3,.., 17 # yield averages of 13 genotypes in localities 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) #example 2 covariable. precipitation precipitation<- c(1000,1100,1200,1300,1400,1500,1600,1700,1800,1900,2000,2100) stability.par(dataStb, rep=4, MSerror=1.8, alpha=0.1, main="Genotype", cova=TRUE, name.cov="Precipitation", file.cov=precipitation,console=TRUE)
By this process the variance and central measures ar found: average, medium and mode of grouped data.
stat.freq(histogram)stat.freq(histogram)
histogram |
Object create by function hist() |
Statistics of grouped data.
Felipe de mendiburu
polygon.freq, table.freq, graph.freq,
inter.freq, sturges.freq, join.freq,
ogive.freq, normal.freq
library(agricolae) data(growth) grouped<-with(growth,hist(height,plot=FALSE)) measures<-stat.freq(grouped) print(measures)library(agricolae) data(growth) grouped<-with(growth,hist(height,plot=FALSE)) measures<-stat.freq(grouped) print(measures)
The variance analysis of a strip-plot design is divided into three parts: the horizontal-factor analysis, the vertical-factor analysis, and the interaction analysis.
strip.plot(BLOCK, COL, ROW, Y)strip.plot(BLOCK, COL, ROW, Y)
BLOCK |
replications |
COL |
Factor column |
ROW |
Factor row |
Y |
Variable, response |
The strip-plot design is specifically suited for a two-factor experiment in which the desired precision for measuring the interaction effects between the two factors is higher than that for measuring the main efect two factors
Data and analysis of the variance of the strip plot design.
Felipe de Mendiburu
Statistical procedures for agricultural research. Kwanchai A. Gomez, Arturo A. Gomez. Second Edition. 1984.
ssp.plot, sp.plot, design.split,
design.strip
# Yield library(agricolae) data(huasahuasi) YIELD<-huasahuasi$YIELD market <- YIELD$y1da + YIELD$y2da non_market <- YIELD$y3da yield <- market + non_market model<-with(YIELD,strip.plot(block, clon, trt, yield)) out1<-with(YIELD,LSD.test(yield,clon,model$gl.a,model$Ea)) oldpar<-par(mar=c(3,8,1,1),cex=0.8) plot(out1,xlim=c(0,80),horiz=TRUE,las=1) out2<-with(YIELD,LSD.test(yield,trt,model$gl.b,model$Eb)) plot(out2,xlim=c(0,80),horiz=TRUE,las=1) par(oldpar)# Yield library(agricolae) data(huasahuasi) YIELD<-huasahuasi$YIELD market <- YIELD$y1da + YIELD$y2da non_market <- YIELD$y3da yield <- market + non_market model<-with(YIELD,strip.plot(block, clon, trt, yield)) out1<-with(YIELD,LSD.test(yield,clon,model$gl.a,model$Ea)) oldpar<-par(mar=c(3,8,1,1),cex=0.8) plot(out1,xlim=c(0,80),horiz=TRUE,las=1) out2<-with(YIELD,LSD.test(yield,trt,model$gl.b,model$Eb)) plot(out2,xlim=c(0,80),horiz=TRUE,las=1) par(oldpar)
if k=0 then classes: k = 1 + log(n,2). if k > 0, fixed nclass.
sturges.freq(x,k=0)sturges.freq(x,k=0)
x |
vector |
k |
constant |
Statistics of sturges for a histogram.
Felipe de mendiburu
Reza A. Hoshmand. 1988. Statistical Methods for Agricultural Sciences, Timber Press, Incorporated, pag 18-21.
polygon.freq, table.freq, stat.freq,
inter.freq, graph.freq, join.freq,
ogive.freq, normal.freq
library(agricolae) data(natives) classes<-with(natives,sturges.freq(size)) # information of the classes breaks <- classes$breaks breaks #startgraph # Histogram with the established classes h<-with(natives,graph.freq(size,breaks,frequency=1, col="yellow",axes=FALSE, xlim=c(0,0.12),main="",xlab="",ylab="")) axis(1,breaks,las=2) axis(2,seq(0,400,50),las=2) title(main="Histogram of frequency\nSize of the tubercule of the Oca", xlab="Size of the oca", ylab="Frequency") #endgraphlibrary(agricolae) data(natives) classes<-with(natives,sturges.freq(size)) # information of the classes breaks <- classes$breaks breaks #startgraph # Histogram with the established classes h<-with(natives,graph.freq(size,breaks,frequency=1, col="yellow",axes=FALSE, xlim=c(0,0.12),main="",xlab="",ylab="")) axis(1,breaks,las=2) axis(2,seq(0,400,50),las=2) title(main="Histogram of frequency\nSize of the tubercule of the Oca", xlab="Size of the oca", ylab="Frequency") #endgraph
It finds the absolute, relative and accumulated frequencies with the class intervals defined from a previously calculated histogram by the "hist" of R function.
## S3 method for class 'graph.freq' summary(object,...)## S3 method for class 'graph.freq' summary(object,...)
object |
Object by function graph.freq() |
... |
other parameters of graphic |
Frequency table.
Lower |
Lower limit class |
Upper |
Upper limit class |
Main |
class point |
Frequency |
Frequency |
Percentage |
Percentage frequency |
CF |
Cumulative frequency |
CPF |
Cumulative Percentage frequency |
Felipe de Mendiburu
polygon.freq, stat.freq, graph.freq,
inter.freq, sturges.freq, join.freq,
ogive.freq, normal.freq
library(agricolae) data(growth) h2<-with(growth,graph.freq(height,plot=FALSE)) print(summary(h2),row.names=FALSE)library(agricolae) data(growth) h2<-with(growth,graph.freq(height,plot=FALSE)) print(summary(h2),row.names=FALSE)
The data correspond to an experiment with costanero sweetpotato made at the locality of the Tacna department, southern Peru. The effect of two viruses (Spfmv and Spcsv) was studied. The treatments were the following: CC (Spcsv) = Sweetpotato chlorotic dwarf, FF (Spfmv) = Feathery mottle, FC (Spfmv y Spcsv) = Viral complex and OO (witness) healthy plants. In each plot, 50 sweetpotato plants were sown and 12 plots were employed. Each treatment was made with 3 repetitions and at the end of the experiment the total weight in kilograms was evaluated. The virus transmission was made in the cuttings and these were sown in the field.
data(sweetpotato)data(sweetpotato)
A data frame with 12 observations on the following 2 variables.
virusa factor with levels cc fc ff oo
yielda numeric vector
Experimental field.
International Potato Center. CIP - Lima Peru
library(agricolae) data(sweetpotato) str(sweetpotato)library(agricolae) data(sweetpotato) str(sweetpotato)
It finds the absolute, relative and accumulated frequencies with the class intervals defined from a previously calculated histogram by the "hist" of R function.
table.freq(object)table.freq(object)
object |
Object by function graph.freq() |
Frequency table.
Lower |
Lower limit class |
Upper |
Upper limit class |
Main |
class point |
Frequency |
Frequency |
Percentage |
Percentage frequency |
CF |
Cumulative frequency |
CPF |
Cumulative Percentage frequency |
Felipe de Mendiburu
polygon.freq, stat.freq, graph.freq,
inter.freq, sturges.freq, join.freq,
ogive.freq, normal.freq
library(agricolae) data(growth) h2<-with(growth,graph.freq(height,plot=FALSE)) print(table.freq(h2),row.names=FALSE)library(agricolae) data(growth) h2<-with(growth,graph.freq(height,plot=FALSE)) print(table.freq(h2),row.names=FALSE)
This process lies in finding statistics which consist of more than one variable, grouped or crossed by factors. The table must be organized by columns between variables and factors.
tapply.stat(y, x, stat = "mean")tapply.stat(y, x, stat = "mean")
y |
data.frame variables |
x |
data.frame factors |
stat |
Method |
Statistics of quantitative variables by categorical variables.
Felipe de Mendiburu
library(agricolae) # case of 1 single factor data(sweetpotato) tapply.stat(sweetpotato[,2],sweetpotato[,1],mean) with(sweetpotato,tapply.stat(yield,virus,sd)) with(sweetpotato,tapply.stat(yield,virus,function(x) max(x)-min(x))) with(sweetpotato,tapply.stat(yield,virus, function(x) quantile(x,0.75,6)-quantile(x,0.25,6))) # other case data(cotton) with(cotton,tapply.stat(yield,cotton[,c(1,3,4)],mean)) with(cotton,tapply.stat(yield,cotton[,c(1,4)],max)) # Height of pijuayo data(growth) with(growth,tapply.stat(height, growth[,2:1], function(x) mean(x,na.rm=TRUE)))library(agricolae) # case of 1 single factor data(sweetpotato) tapply.stat(sweetpotato[,2],sweetpotato[,1],mean) with(sweetpotato,tapply.stat(yield,virus,sd)) with(sweetpotato,tapply.stat(yield,virus,function(x) max(x)-min(x))) with(sweetpotato,tapply.stat(yield,virus, function(x) quantile(x,0.75,6)-quantile(x,0.25,6))) # other case data(cotton) with(cotton,tapply.stat(yield,cotton[,c(1,3,4)],mean)) with(cotton,tapply.stat(yield,cotton[,c(1,4)],max)) # Height of pijuayo data(growth) with(growth,tapply.stat(height, growth[,2:1], function(x) mean(x,na.rm=TRUE)))
The Kendall method in order to find the K variance.
vark(x, y)vark(x, y)
x |
Vector |
y |
vector |
Script in C to R.
variance of K for Kendall's tau
Felipe de Mendiburu
Numerical Recipes in C. Second Edition.
cor.matrix, cor.vector, cor.mv
library(agricolae) x <-c(1,1,1,4,2,2,3,1,3,2,1,1,2,3,2,1,1,2,1,2) y <-c(1,1,2,3,4,4,2,1,2,3,1,1,3,4,2,1,1,3,1,2) vark(x,y)library(agricolae) x <-c(1,1,1,4,2,2,3,1,3,2,1,1,2,3,2,1,1,2,1,2) y <-c(1,1,2,3,4,4,2,1,2,3,1,1,3,4,2,1,1,3,1,2) vark(x,y)
A nonparametric test for several independent samples.
waerden.test(y, trt, alpha=0.05, group=TRUE, main=NULL,console=FALSE)waerden.test(y, trt, alpha=0.05, group=TRUE, main=NULL,console=FALSE)
y |
Variable response |
trt |
Treatments |
alpha |
Significant level |
group |
TRUE or FALSE |
main |
Title |
console |
logical, print output |
The data consist of k samples of possibly unequal sample size.
The post hoc test is using the criterium Fisher's least
significant difference (LSD).
statistics |
Statistics of the model |
parameters |
Design parameters |
means |
Statistical summary of the study variable |
comparison |
Comparison between treatments |
groups |
Formation of treatment groups |
Felipe de Mendiburu
Practical Nonparametrics Statistics. W.J. Conover, 1999
BIB.test, DAU.test, duncan.test,
durbin.test, friedman, HSD.test,
kruskal, LSD.test, Median.test,
PBIB.test, REGW.test, scheffe.test,
SNK.test, waller.test, plot.group
library(agricolae) # example 1 data(corn) out1<-with(corn,waerden.test(observation,method,group=TRUE)) print(out1$groups) plot(out1) out2<-with(corn,waerden.test(observation,method,group=FALSE)) print(out2$comparison) # example 2 data(sweetpotato) out<-with(sweetpotato,waerden.test(yield,virus,alpha=0.01,group=TRUE)) print(out)library(agricolae) # example 1 data(corn) out1<-with(corn,waerden.test(observation,method,group=TRUE)) print(out1$groups) plot(out1) out2<-with(corn,waerden.test(observation,method,group=FALSE)) print(out2$comparison) # example 2 data(sweetpotato) out<-with(sweetpotato,waerden.test(yield,virus,alpha=0.01,group=TRUE)) print(out)
A Bayes rule for the symmetric multiple comparisons problem.
waller(K, q, f, Fc)waller(K, q, f, Fc)
K |
Is the loss ratio between type I and type II error |
q |
Numerator Degrees of freedom |
f |
Denominator Degrees of freedom |
Fc |
F ratio from an analysis of variance |
K-RATIO (K): value specifies the Type 1/Type 2 error seriousness ratio for the Waller-Duncan test. Reasonable values for KRATIO are 50, 100, and 500, which roughly correspond for the two-level case to ALPHA levels of 0.1, 0.05, and 0.01. By default, the procedure uses the default value of 100.
Waller value for the Waller and Duncan test.
Felipe de Mendiburu
Waller, R. A. and Duncan, D. B. (1969). A Bayes Rule for the Symmetric Multiple Comparison Problem, Journal of the American Statistical Association 64, pages 1484-1504.
Waller, R. A. and Kemp, K. E. (1976) Computations of Bayesian t-Values for Multiple Comparisons, Journal of Statistical Computation and Simulation, 75, pages 169-172.
Principles and procedures of statistics a biometrical approach Steel & Torry & Dickey. Third Edition 1997.
# Table Duncan-Waller K=100, F=1.2 pag 649 Steel & Torry library(agricolae) K<-100 Fc<-1.2 q<-c(8,10,12,14,16,20,40,100) f<-c(seq(4,20,2),24,30,40,60,120) 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) print(W.D)# Table Duncan-Waller K=100, F=1.2 pag 649 Steel & Torry library(agricolae) K<-100 Fc<-1.2 q<-c(8,10,12,14,16,20,40,100) f<-c(seq(4,20,2),24,30,40,60,120) 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) print(W.D)
The Waller-Duncan k-ratio t test is performed on all main effect means in the MEANS statement. See the K-RATIO option for information on controlling details of the test.
waller.test(y, trt, DFerror, MSerror, Fc, K = 100, group=TRUE, main = NULL, console=FALSE)waller.test(y, trt, DFerror, MSerror, Fc, K = 100, group=TRUE, main = NULL, console=FALSE)
y |
model(aov or lm) or answer of the experimental unit |
trt |
Constant( only y=model) or vector treatment applied to each unit |
DFerror |
Degrees of freedom |
MSerror |
Mean Square Error |
Fc |
F Value |
K |
K-RATIO |
group |
TRUE or FALSE |
main |
Title |
console |
logical, print output |
It is necessary first makes a analysis of variance.
K-RATIO (K): value specifies the Type 1/Type 2 error seriousness ratio for
the Waller-Duncan test. Reasonable values for KRATIO are 50, 100, and 500,
which roughly correspond for the two-level case to ALPHA levels of 0.1, 0.05,
and 0.01. By default, the procedure uses the default value of 100.
if y = model, then to apply the instruction:
waller.test (model, "trt", alpha = 0.05, group = TRUE, main = NULL, console = FALSE)
where the model class is aov or lm.
statistics |
Statistics of the model |
parameters |
Design parameters |
means |
Statistical summary of the study variable |
comparison |
Comparison between treatments |
groups |
Formation of treatment groups |
Felipe de Mendiburu
Waller, R. A. and Duncan, D. B. (1969).
A Bayes Rule for the Symmetric Multiple Comparison Problem,
Journal of the American Statistical Association 64, pages 1484-1504.
Waller, R. A. and Kemp, K. E. (1976)
Computations of Bayesian t-Values for Multiple Comparisons,
Journal of Statistical Computation and Simulation, 75, pages 169-172.
Steel & Torry & Dickey. Third Edition 1997 Principles and procedures of statistics a biometrical approach
BIB.test, DAU.test, duncan.test,
durbin.test, friedman, HSD.test,
kruskal, LSD.test, Median.test,
PBIB.test, REGW.test, scheffe.test,
SNK.test, waerden.test, plot.group
library(agricolae) data(sweetpotato) model<-aov(yield~virus, data=sweetpotato) out <- waller.test(model,"virus", group=TRUE) #startgraph oldpar<-par(mfrow=c(2,2)) # variation: SE is error standard # variation: range is Max - Min bar.err(out$means,variation="SD",horiz=TRUE,xlim=c(0,45),bar=FALSE, col=colors()[25],space=2, main="Standard deviation",las=1) bar.err(out$means,variation="SE",horiz=FALSE,ylim=c(0,45),bar=FALSE, col=colors()[15],space=2,main="SE",las=1) bar.err(out$means,variation="range",ylim=c(0,45),bar=FALSE,col="green", space=3,main="Range = Max - Min",las=1) bar.group(out$groups,horiz=FALSE,ylim=c(0,45),density=8,col="red", main="Groups",las=1) #endgraph # Old version HSD.test() df<-df.residual(model) MSerror<-deviance(model)/df Fc<-anova(model)["virus",4] out <- with(sweetpotato,waller.test(yield, virus, df, MSerror, Fc, group=TRUE)) print(out) par(oldpar)library(agricolae) data(sweetpotato) model<-aov(yield~virus, data=sweetpotato) out <- waller.test(model,"virus", group=TRUE) #startgraph oldpar<-par(mfrow=c(2,2)) # variation: SE is error standard # variation: range is Max - Min bar.err(out$means,variation="SD",horiz=TRUE,xlim=c(0,45),bar=FALSE, col=colors()[25],space=2, main="Standard deviation",las=1) bar.err(out$means,variation="SE",horiz=FALSE,ylim=c(0,45),bar=FALSE, col=colors()[15],space=2,main="SE",las=1) bar.err(out$means,variation="range",ylim=c(0,45),bar=FALSE,col="green", space=3,main="Range = Max - Min",las=1) bar.group(out$groups,horiz=FALSE,ylim=c(0,45),density=8,col="red", main="Groups",las=1) #endgraph # Old version HSD.test() df<-df.residual(model) MSerror<-deviance(model)/df Fc<-anova(model)["virus",4] out <- with(sweetpotato,waller.test(yield, virus, df, MSerror, Fc, group=TRUE)) print(out) par(oldpar)
Weather and Severity
weatherSeverity(weather,severity,dates,EmergDate,EndEpidDate,NoReadingsH, RHthreshold)weatherSeverity(weather,severity,dates,EmergDate,EndEpidDate,NoReadingsH, RHthreshold)
weather |
object, see example |
severity |
object, see example |
dates |
vector dates |
EmergDate |
date |
EndEpidDate |
date |
NoReadingsH |
num, 1 |
RHthreshold |
num, percentage |
Weather and severity
Wfile |
"Date","Rainfall","Tmp","HumidHrs","humidtmp" |
Sfile |
"Cultivar","ApplSys","dates","nday","MeanSeverity","StDevSeverity" |
EmergDate |
date |
EndEpidDate |
date |
All format data for date is yyyy-mm,dd, for example "2000-04-22". change with function as.Date()
library(agricolae) 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 and threshold 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)library(agricolae) 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 and threshold 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)
Percentage of bacterial wilt and area under the curve of disease progression (AUDPC) relative tomato plants transplanted in different soil types artificially infested with R.solanacearum 133 days before.
data(wilt)data(wilt)
A data frame with 13 observations on the following 15 variables.
placea factor with levels Chmar Chz Cnt1 Cnt2 Cnt3 Hco1 Hco2 Hco3 Hyo1 Hyo2 Namora SR1 SR2
Day7a numeric vector
Day11a numeric vector
Day15a numeric vector
Day19a numeric vector
Day23a numeric vector
Day27a numeric vector
Day31a numeric vector
Day35a numeric vector
Day39a numeric vector
Day43a numeric vector
Day47a numeric vector
Day51a numeric vector
AUDPCa numeric vector
relativea numeric vector
Percentajes bacterial wilt. Day7 = evaluated to 7 days, Days11 = evaluated to 11 days. see data(soil) and data(ralstonia)
Experimental field, 2004. Data Kindly provided by Dr. Sylvie Priou, Liliam Gutarra and Pedro Aley.
International Potato Center. CIP - Lima Peru.
library(agricolae) data(wilt) days<-c(7,11,15,19,23,27,31,35,39,43,47,51) AUDPC<-audpc(wilt[,-1],days) relative<-audpc(wilt[,-1],days,type="relative")library(agricolae) data(wilt) days<-c(7,11,15,19,23,27,31,35,39,43,47,51) AUDPC<-audpc(wilt[,-1],days) relative<-audpc(wilt[,-1],days,type="relative")
The yacon (Smallanthus sonchifolius) is a plant native to the Andes, considered a traditional crop in Peru and natural source of FOS, which is a type of carbohydrate that can not be digested by the and the human body that have joined several beneficial properties in health, such as improve the absorption of calcium, reducing the level of triglycerides and cholesterol and stimulate better gastrointestinal function.
data(yacon)data(yacon)
A data frame with 432 observations on the following 19 variables.
localitya factor with levels, Cajamarca, Lima, Oxapampa in PERU
sitea numeric vector
dosea factor with levels F0 F150 F80
entrya factor with levels AKW5075 AMM5136 AMM5150 AMM5163 ARB5125 CLLUNC118 P1385 SAL136
replicationa numeric vector, replications
heighta numeric vector, plant height, centimeters
stalksa numeric vector, number of stalks
wfra numeric vector, weight of fresh roots, grams
wffa numeric vector, weight of fresh foliage, grams
wfka numeric vector, weight fresh kroner, grams
rootsa numeric vector, matter of dried roots, grams
FOSa numeric vector, fructo-oligosaccharides, percentaje
glucosea numeric vector, percentaje
fructosea numeric vector, percentaje
sucrosea numeric vector, percentaje
brixa numeric vector, degrees Brix
foliagea numeric vector, matter dry foliage, grams
drya numeric vector, dry matter kroner, grams
IHa numeric vector, Index harvest, 0 to 1
Proportion or fraction of the plant that is used (seeds, fruit, root) on dry basis. Part usable in a proportion of total mass dissected. Plant of frijol, weight = 100g and frijol = 50g then, IH = 50/100 = 0.5 or 50 percentaje. Degrees Brix is a measurement of the mass ratio of dissolved sugar to water in a liquid.
CIP. Experimental field, 2003, Data Kindly provided by Ivan Manrique and Carolina Tasso.
International Potato Center. CIP - Lima Peru.
library(agricolae) data(yacon) str(yacon)library(agricolae) data(yacon) str(yacon)
applied to designs: complete block, latin square, graeco, split plot, strip plot, lattice, alpha lattice, Augmented block, cyclic, Balanced Incomplete Block and factorial.
zigzag(outdesign)zigzag(outdesign)
outdesign |
output design |
fieldbook |
Remuneration of serpentine plots. |
Felipe de Mendiburu
design.ab, design.alpha,design.bib,
design.split , design.cyclic , design.dau ,
design.graeco, design.lattice, design.lsd,
design.rcbd, design.strip
library(agricolae) trt<-letters[1:5] r<-4 outdesign <- design.rcbd(trt,r,seed=9) fieldbook <- zigzag(outdesign)library(agricolae) trt<-letters[1:5] r<-4 outdesign <- design.rcbd(trt,r,seed=9) fieldbook <- zigzag(outdesign)