Type:	Package
Title:	Agricultural Datasets
Version:	1.24
Description:	Datasets from books, papers, and websites related to agriculture. Example graphics and analyses are included. Data come from small-plot trials, multi-environment trials, uniformity trials, yield monitors, and more.
License:	MIT + file LICENSE
URL:	https://kwstat.github.io/agridat/
BugReports:	https://github.com/kwstat/agridat/issues
Suggests:	AER, agricolae, betareg, broom, car, coin, corrgram, desplot, dplyr, effects, emmeans, equivalence, FrF2, gam, gge, ggplot2, gnm, gstat, HH, knitr, lattice, latticeExtra, lme4, lucid, mapproj, maps, MASS, MCMCglmm, metafor, mgcv, NADA, nlme, nullabor, ordinal, pbkrtest, pls, pscl, qicharts, qtl, reshape2, rmarkdown, sp, SpATS, survival, testthat, vcd
VignetteBuilder:	knitr
Encoding:	UTF-8
Language:	en-US
LazyData:	true
RoxygenNote:	7.3.2
NeedsCompilation:	no
Packaged:	2024-10-27 16:04:47 UTC; wrightkevi
Author:	Kevin Wright [aut, cre, cph]
Maintainer:	Kevin Wright <kw.stat@gmail.com>
Repository:	CRAN
Date/Publication:	2024-10-27 17:00:02 UTC

Datasets from agricultural experiments

Description

This package contains datasets from publications relating to agriculture, including field crops, tree crops, animal studies, and a few others.

Details

If you use these data, please cite both the agridat package and the original source of the data.

Abbreviations in the 'other' column include: xy = coordinates, pls = partial least squares, rsm = response surface methodology, row-col = row-column design, ts = time series,

Uniformity trials with a single genotype

Yield monitor

Animals

Trees

Field and horticulture crops

Time series

Other

Summaries:

Diallel experiments:

Factorial experiments:

Multi-environment trials with multi-genotype,loc,rep,year:

Data with markers: hadasch.lettuce.markers, steptoe.morex.geno

Data with pedigree: butron.maize

Author(s)

Kevin Wright, with support from many people who granted permission to include their data in this package.

References

J. White and Frits van Evert. (2008). Publishing Agronomic Data. Agron J. 100, 1396-1400. https://doi.org/10.2134/agronj2008.0080F

Barley heights and environmental covariates in Norway

Description

Average height for 15 genotypes of barley in each of 9 years. Also 19 covariates in each of the 9 years.

Usage

data("aastveit.barley.covs")
data("aastveit.barley.height")

Format

The 'aastveit.barley.covs' dataframe has 9 observations on the following 20 variables.

year

year

R1

avg rainfall (mm/day) in period 1

R2

avg rainfall (mm/day) in period 2

R3

avg rainfall (mm/day) in period 3

R4

avg rainfall (mm/day) in period 4

R5

avg rainfall (mm/day) in period 5

R6

avg rainfall (mm/day) in period 6

S1

daily solar radiation (ca/cm^2) in period 1

S2

daily solar radiation (ca/cm^2) in period 2

S3

daily solar radiation (ca/cm^2) in period 3

S4

daily solar radiation (ca/cm^2) in period 4

S5

daily solar radiation (ca/cm^2) in period 5

S6

daily solar radiation (ca/cm^2) in period 6

ST

sowing time, measured in days after April 1

T1

avg temp (deg Celsius) in period 1

T2

avg temp (deg Celsius) in period 2

T3

avg temp (deg Celsius) in period 3

T4

avg temp (deg Celsius) in period 4

T5

avg temp (deg Celsius) in period 5

T6

avg temp (deg Celsius) in period 6

The 'aastveit.barley.height' dataframe has 135 observations on the following 3 variables.

year

year, 9 years spanning from 1974 to 1982

gen

genotype, 15 levels

height

height (cm)

Details

Experiments were conducted at As, Norway.

The height dataframe contains average plant height (cm) of 15 varieties of barley in each of 9 years.

The growth season of each year was divided into eight periods from sowing to harvest. Because the plant stop growing about 20 days after ear emergence, only the first 6 periods are included here.

Used with permission of Harald Martens.

Source

Aastveit, A. H. and Martens, H. (1986). ANOVA interactions interpreted by partial least squares regression. Biometrics, 42, 829–844. https://doi.org/10.2307/2530697

References

J. Chadoeuf and J. B. Denis (1991). Asymptotic variances for the multiplicative interaction model. J. App. Stat., 18, 331-353. https://doi.org/10.1080/02664769100000032

Examples

## Not run: 

library(agridat)
data("aastveit.barley.covs")
data("aastveit.barley.height")

libs(reshape2, pls)
  
  # First, PCA of each matrix separately

  Z <- acast(aastveit.barley.height, year ~ gen, value.var="height")
  Z <- sweep(Z, 1, rowMeans(Z))
  Z <- sweep(Z, 2, colMeans(Z)) # Double-centered
  sum(Z^2)*4 # Total SS = 10165
  sv <- svd(Z)$d
  round(100 * sv^2/sum(sv^2),1) # Prop of variance each axis
  # Aastveit Figure 1.  PCA of height
  biplot(prcomp(Z),
         main="aastveit.barley - height", cex=0.5)
  
  U <- aastveit.barley.covs
  rownames(U) <- U$year
  U$year <- NULL
  U <- scale(U) # Standardized covariates
  sv <- svd(U)$d
  # Proportion of variance on each axis
  round(100 * sv^2/sum(sv^2),1)

  # Now, PLS relating the two matrices
  m1 <- plsr(Z~U)
  loadings(m1)
  # Aastveit Fig 2a (genotypes), but rotated differently
  biplot(m1, which="y", var.axes=TRUE)
  # Fig 2b, 2c (not rotated)
  biplot(m1, which="x", var.axes=TRUE)

  # Adapted from section 7.4 of Turner & Firth,
  # "Generalized nonlinear models in R: An overview of the gnm package"
  # who in turn reproduce the analysis of Chadoeuf & Denis (1991),
  # "Asymptotic variances for the multiplicative interaction model"

  libs(gnm)
  dath <- aastveit.barley.height
  dath$year = factor(dath$year)

  set.seed(42)
  m2 <- gnm(height ~ year + gen + Mult(year, gen), data = dath)
  # Turner: "To obtain parameterization of equation 1, in which sig_k is the
  # singular value for component k, the row and column scores must be constrained
  # so that the scores sum to zero and the squared scores sum to one.
  # These contrasts can be obtained using getContrasts"
  gamma <- getContrasts(m2, pickCoef(m2, "[.]y"),
                        ref = "mean", scaleWeights = "unit")
  delta <- getContrasts(m2, pickCoef(m2, "[.]g"),
                        ref = "mean", scaleWeights = "unit")
  # estimate & std err
  gamma <- gamma$qvframe
  delta <- delta$qvframe
  # change sign of estimate
  gamma[,1] <- -1 * gamma[,1]
  delta[,1] <- -1 * delta[,1]
  # conf limits based on asymptotic normality, Chadoeuf table 8, p. 350, 
  round(cbind(gamma[,1], gamma[, 1] +
                           outer(gamma[, 2], c(-1.96, 1.96))) ,3)  
  round(cbind(delta[,1], delta[, 1] +
                           outer(delta[, 2], c(-1.96, 1.96))) ,3)

## End(Not run)

Multi-environment trial evaluating 36 maize genotypes in 9 locations

Description

Multi-environment trial evaluating 36 maize genotypes in 9 locations

Usage

data("acorsi.grayleafspot")

Format

A data frame with 324 observations on the following 3 variables.

gen

genotype, 36 levels

env

environment, 9 levels

rep

replicate, 2 levels

y

grey leaf spot severity

Details

Experiments conducted in 9 environments in Brazil in 2010-11. Each location had an RCB with 2 reps.

The response variable is the percentage of leaf area affected by gray leaf spot within each experimental unit (plot).

Acorsi et al. use this data to illustrate the fitting of a generalized AMMI model with non-normal data.

Source

C. R. L. Acorsi, T. A. Guedes, M. M. D. Coan, R. J. B. Pinto, C. A. Scapim, C. A. P. Pacheco, P. E. O. Guimaraes, C. R. Casela. (2016). Applying the generalized additive main effects and multiplicative interaction model to analysis of maize genotypes resistant to grey leaf spot. Journal of Agricultural Science. https://doi.org/10.1017/S0021859616001015

Electronic data and R code kindly provided by Marlon Coan.

References

None

Examples

## Not run: 

  library(agridat)
  data(acorsi.grayleafspot)
  dat <- acorsi.grayleafspot
  
  # Acorsi figure 2. Note: Acorsi used cell means
  op <- par(mfrow=c(2,1), mar=c(5,4,3,2))
  libs(lattice)
  boxplot(y ~ env, dat, las=2,
          xlab="environment", ylab="GLS severity")
  title("acorsi.grayleafspot")
  boxplot(y ~ gen, dat, las=2,
          xlab="genotype", ylab="GLS severity")
  par(op)
  
  # GLM models
  
  # glm main-effects model with logit u(1-u) and wedderburn u^2(1-u)^2
  # variance functions
  # glm1 <- glm(y~ env/rep + gen + env, data=dat, family=quasibinomial)
  # glm2 <- glm(y~ env/rep + gen + env, data=dat, family=wedderburn)
  # plot(glm2, which=1); plot(glm2, which=2)
  
  # GAMMI models of Acorsi. See also section 7.4 of Turner
  # "Generalized nonlinear models in R: An overview of the gnm package"
  
  # full gnm model with wedderburn, seems to work
  libs(gnm)
  set.seed(1)
  gnm1 <- gnm(y ~  env/rep + env + gen + instances(Mult(env,gen),2),
              data=dat,
              family=wedderburn, iterMax =800)
  deviance(gnm1) # 433.8548
  # summary(gnm1)
  # anova(gnm1, test ="F")  # anodev, Acorsi table 4
  ##                          Df Deviance Resid. Df Resid. Dev       F    Pr(>F)    
  ## NULL                                       647     3355.5                      
  ## env                       8  1045.09       639     2310.4 68.4696 < 2.2e-16 ***
  ## env:rep                   9    12.33       630     2298.1  0.7183    0.6923    
  ## gen                      35  1176.23       595     1121.9 17.6142 < 2.2e-16 ***
  ## Mult(env, gen, inst = 1) 42   375.94       553      745.9  4.6915 < 2.2e-16 ***
  ## Mult(env, gen, inst = 2) 40   312.06       513      433.9  4.0889 3.712e-14 ***


  # maybe better, start simple and build up the model
  gnm2a <- gnm(y ~  env/rep + env + gen,
               data=dat,
               family=wedderburn, iterMax =800)

  # add first interaction term
  res2a <- residSVD(gnm2a, env, gen, 2)
  gnm2b <- update(gnm2a, . ~ . + Mult(env,gen,inst=1),
                  start = c(coef(gnm2a), res2a[, 1]))
  deviance(gnm2b) # 692.19

  # add second interaction term
  res2b <- residSVD(gnm2b, env, gen, 2)
  gnm2c <- update(gnm2b, . ~ . + Mult(env,gen,inst=1) + Mult(env,gen,inst=2),
                  start = c(coef(gnm2a), res2a[, 1], res2b[,1]))
  deviance(gnm2c) # 433.8548
  # anova(gnm2c) # weird error message

  # note, to build the ammi biplot, use the first column of res2a to get
  # axis 1, and the FIRST column of res2b to get axis 2. Slightly confusing
  emat <- cbind(res2a[1:9, 1], res2b[1:9, 1])
  rownames(emat) <- gsub("fac1.", "", rownames(emat))
  
  gmat <- cbind(res2a[10:45, 1], res2b[10:45, 1])
  rownames(gmat) <- gsub("fac2.", "", rownames(gmat))

  # match Acorsi figure 4
  biplot(gmat, emat, xlim=c(-2.2, 2.2), ylim=c(-2.2, 2.2), expand=2, cex=0.5,
         xlab="Axis 1", ylab="Axis 2",
         main="acorsi.grayleafspot - GAMMI biplot")

## End(Not run)

Multi-environment trial of sorghum at 3 locations across 5 years

Description

Multi-environment trial of sorghum at 3 locations across 5 years

Format

A data frame with 289 observations on the following 6 variables.

gen

genotype, 28 levels

trial

trial, 2 levels

env

environment, 13 levels

yield

yield kg/ha

year

year, 2001-2005

loc

location, 3 levels

Details

Sorghum yields at 3 locations across 5 years. The trials were carried out at three locations in dry, hot lowlands of Ethiopia:

Melkassa (39 deg 21 min E, 8 deg 24 min N)

Mieso (39 deg 22 min E, 8 deg 41 min N)

Kobo (39 deg 37 min E, 12 deg 09 min N)

Trial 1 was 14 hybrids and one open-pollinated variety.

Trial 2 was 12 experimental lines.

Used with permission of Asfaw Adugna.

Source

Asfaw Adugna (2008). Assessment of yield stability in sorghum using univariate and multivariate statistical approaches. Hereditas, 145, 28–37. https://doi.org/10.1111/j.0018-0661.2008.2023.x

Examples

## Not run: 

library(agridat)
data(adugna.sorghum)
dat <- adugna.sorghum

libs(lattice)
redblue <- colorRampPalette(c("firebrick", "lightgray", "#375997"))
levelplot(yield ~ env*gen, data=dat, main="adugna.sorghum gxe heatmap",
          col.regions=redblue)

# Genotype means match Adugna
tapply(dat$yield, dat$gen, mean)

# CV for each genotype.  G1..G15 match, except for G2.
# The table in Adugna scrambles the means for G16..G28
libs(reshape2)
mat <- acast(dat, gen~env,  value.var='yield')
round(sqrt(apply(mat, 1, var, na.rm=TRUE)) / apply(mat, 1, mean, na.rm=TRUE) * 100,2)

# Shukla stability.  G1..G15 match Adugna.  Can't match G16..G28.
dat1 <- droplevels(subset(dat, trial=="T1"))
mat1 <- acast(dat1, gen~env,  value.var='yield')
w <- mat1; k=15; n=8  # k=p gen, n=q env
w <- sweep(w, 1, rowMeans(mat1, na.rm=TRUE))
w <- sweep(w, 2, colMeans(mat1, na.rm=TRUE))
w <- w + mean(mat1, na.rm=TRUE)
w <- rowSums(w^2, na.rm=TRUE)
sig2 <- k*w/((k-2)*(n-1)) - sum(w)/((k-1)*(k-2)*(n-1))
round(sig2/10000,1) # Genotypes in T1 are divided by 10000

## End(Not run)

Multi-environment trial of cereal with lodging data

Description

Percent lodging is given for 32 genotypes at 7 environments.

Format

A data frame with 224 observations on the following 3 variables.

env

environment, 1-7

gen

genotype, 1-32

y

percent lodged

Details

This data is for the first year of a three-year study.

Used with permission of Chris Glasbey.

Source

D. J. Allcroft and C. A. Glasbey, 2003. Analysis of crop lodging using a latent variable model. Journal of Agricultural Science, 140, 383–393. https://doi.org/10.1017/S0021859603003332

Examples

## Not run: 

library(agridat)
data(allcroft.lodging)
dat <- allcroft.lodging

# Transformation
dat$sy <- sqrt(dat$y)
# Variety 4 has no lodging anywhere, so add a small amount
dat[dat$env=='E5' & dat$gen=='G04',]$sy <- .01

libs(lattice)
dotplot(env~y|gen, dat, as.table=TRUE,
        xlab="Percent lodged (by genotype)", ylab="Variety",
        main="allcroft.lodging")

# Tobit model
libs(AER)
m3 <- tobit(sy ~ 1 + gen + env, left=0, right=100, data=dat)

# Table 2 trial/variety means
preds <- expand.grid(gen=levels(dat$gen), env=levels(dat$env))
preds$pred <- predict(m3, newdata=preds)
round(tapply(preds$pred, preds$gen, mean),2)
round(tapply(preds$pred, preds$env, mean),2)


## End(Not run)

For the 34 sheep sires, the number of lambs in each of 5 foot shape classes.

Description

For the 34 sheep sires, the number of lambs in each of 5 foot shape classes.

Usage

data("alwan.lamb")

Format

A data frame with 340 observations on the following 11 variables.

year

numeric 1980/1981

breed

breed PP, BRP, BR

sex

sex of lamb M/F

sire0

sire ID according to Alwan

shape

sire ID according to Gilmour

count

number of lambs

sire

shape of foot

yr

numeric contrast for year

b1

numeric contrast for breeds

b2

numeric contrast for breeds

b3

numeric contrast for breeds

Details

There were 2513 lambs classified on the presence of deformities in their feet. The lambs represent the offspring of 34 sires, 5 strains, 2 years.

The variables yr, b1, b2, b3 are numeric contrasts for the fixed effects as defined in the paper by Gilmour (1987) and used in the SAS example. Gilmour does not explain the reason for the particular contrasts. The counts for classes LF1, LF2, LF3 were combined.

Source

Mohammed Alwan (1983). Studies of the flock mating performance of Booroola merino crossbred ram lambs, and the foot conditions in Booroola merino crossbreds and Perendale sheep grazed on hill country. Thesis, Massey University. https://hdl.handle.net/10179/5900 Appendix I, II.

References

Gilmour, Anderson, and Rae (1987). Variance components on an underlying scale for ordered multiple threshold categorical data using a generalized linear mixed model. Journal of Animal Breeding and Genetics, 104, 149-155. https://doi.org/10.1111/j.1439-0388.1987.tb00117.x

SAS/STAT(R) 9.2 Users Guide, Second Edition Example 38.11 Maximum Likelihood in Proportional Odds Model with Random Effects https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm

Examples

## Not run: 

  library(agridat)
  data(alwan.lamb)
  dat <- alwan.lamb

  # merge LF1 LF2 LF3 class counts, and combine M/F
  dat$shape <- as.character(dat$shape)
  dat$shape <- ifelse(dat$shape=="LF2", "LF3", dat$shape)
  dat$shape <- ifelse(dat$shape=="LF1", "LF3", dat$shape)
  dat <- aggregate(count ~ year+breed+sire0+sire+shape+yr+b1+b2+b3,
                   dat, FUN=sum)

  dat <- transform(dat,
                   year=factor(year), breed=factor(breed),
                   sire0=factor(sire0), sire=factor(sire))
  # LF5 or LF3 first is a bit arbitary...affects the sign of the coefficients
  dat <- transform(dat, shape=ordered(shape, levels=c("LF5","LF4","LF3")))
  
  # View counts by year and breed
  libs(latticeExtra)
  dat2 <- aggregate(count ~ year+breed+shape, dat, FUN=sum)
  useOuterStrips(barchart(count ~ shape|year*breed, data=dat2,
                          main="alwan.lamb"))

  # Model used by Gilmour and SAS
  dat <- subset(dat, count > 0) 
  libs(ordinal)
  m1 <- clmm(shape ~ yr + b1 + b2 + b3 + (1|sire), data=dat,
             weights=count, link="probit", Hess=TRUE)
  summary(m1) # Very similar to Gilmour results
  ordinal::ranef(m1) # sign is opposite of SAS

  ## SAS var of sires .04849
  ## Effect 	Shape 	Estimate 	Standard Error 	DF 	t Value 	Pr > |t|
  ## Intercept 	1 	0.3781 	0.04907 	29 	7.71 	<.0001
  ## Intercept 	2 	1.6435 	0.05930 	29 	27.72 	<.0001
  ## yr 	  	0.1422 	0.04834 	2478 	2.94 	0.0033
  ## b1 	  	0.3781 	0.07154 	2478 	5.28 	<.0001
  ## b2 	  	0.3157 	0.09709 	2478 	3.25 	0.0012
  ## b3 	  	-0.09887 	0.06508 	2478 	-1.52 	0.1289
  
  ## Gilmour results for probit analysis
  ## Int1   .370 +/- .052
  ## Int2  1.603 +/- .061
  ## Year  -.139 +/- .052
  ## B1    -.370 +/- .076
  ## B2    -.304 +/- .103
  ## B3     .098 +/- .070

  # Plot random sire effects with intervals, similar to SAS example
  plot.random <- function(model, random.effect, ylim=NULL, xlab="", main="") {
    tab <- ordinal::ranef(model)[[random.effect]]
    tab <- data.frame(lab=rownames(tab), est=tab$"(Intercept)")
    tab <- transform(tab,
                     lo = est - 1.96 * sqrt(model$condVar),
                     hi = est + 1.96 * sqrt(model$condVar))
    # sort by est, and return index
    ix <- order(tab$est)
    tab <- tab[ix,]
    
    if(is.null(ylim)) ylim <- range(c(tab$lo, tab$hi))
    n <- nrow(tab)
    plot(1:n, tab$est, axes=FALSE, ylim=ylim, xlab=xlab,
         ylab="effect", main=main, type="n")
    text(1:n, tab$est, labels=substring(tab$lab,2) , cex=.75)
    axis(1)
    axis(2)
    segments(1:n, tab$lo, 1:n, tab$hi, col="gray30")
    abline(h=c(-.5, -.25, 0, .25, .5), col="gray")
    return(ix)  
  }
  ix <- plot.random(m1, "sire")

  # foot-shape proportions for each sire, sorted by estimated sire effects
  # positive sire effects tend to have lower proportion of lambs in LF4 and LF5
  tab <- prop.table(xtabs(count ~ sire+shape, dat), margin=1)
  tab <- tab[ix,]
  tab <- tab[nrow(tab):1,] # reverse the order
  lattice::barchart(tab,
                    horizontal=FALSE, auto.key=TRUE,
                    main="alwan.lamb", xlab="Sire", ylab="Proportion of lambs",
                    scales=list(x=list(rot=70)),
                    par.settings = simpleTheme(col=c("yellow","orange","red")) )
  
  detach("package:ordinal") # to avoid VarCorr clash with lme4
  

## End(Not run)

Uniformity trial of wheat

Description

Uniformity trial of wheat in India in 1940.

Usage

data("ansari.wheat.uniformity")

Format

A data frame with 768 observations on the following 3 variables.

row

row

col

column

yield

yield of grain per plot, in half-ounces

Details

An experiment was conducted at the Government Research Farm, Raya (Muttra District), during the rainy season of 1939-40.

"Wheat was sown over an area of 180 ft. x 243 ft. with 324 rows on a field of average fertility. It had wheat during 1938-39 rabi and was fallow during 1939-40 kharif. The seed was sown behind desi plough in rows 9 inches apart, the length of each row being 180 feet".

"At the time of harvest, 18 rows on both sides and 10 feet at the end of the field were discarded to eliminate border effects and an area of 160 feet x 216 feet with 288 rows was harvested in small units, each being 2 feet 3 inches broad with three rows 20 feet long. There were 96 units across the rows and eight units along the rows. The total number of unit plots thus obtained was 768. The yield of grain for each unit plot was weighed and recorded separately and is given in the appendix."

Field width: 96 plots * 2.25 feet = 216 feet.

Field length: 8 plots * 20 feet = 160 feet.

Comment: There seems to be a strong cyclical patern to the fertility gradient. "History of the field reveals no explanation for this phenomenon, as an average field usually found on the farm was selected for the trial."

Source

Ansari, M. A. A., and G. K. Sant (1943). A Study of Soil Heterogeneity in Relation to Size and Shape of Plots in a Wheat Field at Raya (Muhra District). Ind. J. Agr. Sci, 13, 652-658. https://archive.org/details/in.ernet.dli.2015.271748

References

None

Examples

## Not run: 

  library(agridat)
  data(ansari.wheat.uniformity)
  dat <- ansari.wheat.uniformity

  # match Ansari figure 3
  libs(desplot)
  desplot(dat, yield ~ col*row,
          flip=TRUE, aspect=216/160, # true aspect
          main="ansari.wheat.uniformity") 


## End(Not run)

Uniformity trial of groundnut

Description

Uniformity trial of groundnut

Usage

data("arankacami.groundnut.uniformity")

Format

A data frame with 96 observations on the following 3 variables.

row

row

col

column

yield

yield, kg/plot

Details

The year of the experiment is unknown, but before 1995.

Basic plot size is 0.75 m (rows) x 4 m (columns).

Source

Ira Arankacami, R. Rangaswamy. (1995). A Text Book of Agricultural Statistics. New Age International Publishers. Table 19.1. Page 237. https://www.google.com/books/edition/A_Text_Book_of_Agricultural_Statistics/QDLWE4oakSgC

References

None

Examples

## Not run: 
library(agridat)
data(arankacami.groundnut.uniformity)
dat <- arankacami.groundnut.uniformity

require(desplot)
desplot(dat, yield ~ col*row,
        flip=TRUE, aspect=(12*.75)/(8*4),
        main="arankacami.groundnut.uniformity")


## End(Not run)

Split-split plot experiment of apple trees

Description

Split-split plot experiment of apple trees with different spacing, root stock, and cultivars.

Format

A data frame with 120 observations on the following 10 variables.

rep

block, 5 levels

row

row

pos

position within each row

spacing

spacing between trees, 6,10,14 feet

stock

rootstock, 4 levels

gen

genotype, 2 levels

yield

yield total, kg/tree from 1975-1979

trt

treatment code

Details

In rep 1, the 10-foot-spacing main plot was split into two non-contiguous pieces. This also happened in rep 4. In the analysis of Cornelius and Archbold, they consider each row x within-row-spacing to be a distinct main plot. (Also true for the 14-foot row-spacing, even though the 14-foot spacing plots were contiguous.)

The treatment code is defined as 100 * spacing + 10 * stock + gen, where stock=0,1,6,7 for Seedling,MM111,MM106,M0007 and gen=1,2 for Redspur,Golden, respectively.

Source

D Archbold and G. R. Brown and P. L. Cornelius. (1987). Rootstock and in-row spacing effects on growth and yield of spur-type delicious and Golden delicious apple. Journal of the American Society for Horticultural Science, 112, 219-222.

References

Cornelius, PL and Archbold, DD, 1989. Analysis of a split-split plot experiment with missing data using mixed model equations. Applications of Mixed Models in Agriculture and Related Disciplines. Pages 55-79.

Examples

## Not run: 

library(agridat)
data(archbold.apple)
dat <- archbold.apple

# Define main plot and subplot
dat <- transform(dat, rep=factor(rep), spacing=factor(spacing), trt=factor(trt),
                 mp = factor(paste(row,spacing,sep="")),
                 sp = factor(paste(row,spacing,stock,sep="")))

# Due to 'spacing', the plots are different sizes, but the following layout
# shows the relative position of the plots and treatments. Note that the
# 'spacing' treatments are not contiguous in some reps.
libs(desplot)
desplot(dat, spacing~row*pos,
        col=stock, cex=1, num=gen, # aspect unknown
        main="archbold.apple")


libs(lme4, lucid)  
m1 <- lmer(yield ~ -1 + trt + (1|rep/mp/sp), dat)

vc(m1)  # Variances/means on Cornelius, page 59
##         grp        var1 var2   vcov sdcor
## sp:(mp:rep) (Intercept) <NA>  193.3 13.9
##      mp:rep (Intercept) <NA>  203.8 14.28
##         rep (Intercept) <NA>  197.3 14.05
##    Residual        <NA> <NA> 1015   31.86
  

## End(Not run)

Multi-environment trial of early white food corn

Description

Multi-environment trial of early white food corn for 60 white hybrids.

Format

A data frame with 540 observations on the following 9 variables.

loc

location, 9 levels

gen

gen, 60 levels

yield

yield, bu/ac

stand

stand, percent

rootlodge

root lodging, percent

stalklodge

stalk lodging, percent

earht

ear height, inches

flower

days to flower

moisture

moisture, percent

Details

Data are the average of 3 replications.

Yields were measured for each plot and converted to bushels / acre and adjusted to 15.5 percent moisture.

Stand is expressed as a percentage of the optimum plant stand.

Lodging is expressed as a percentage of the total plants for each hybrid.

Ear height was measured from soil level to the top ear leaf collar. Heights are expressed in inches.

Days to flowering is the number of days from planting to mid-tassel or mid-silk.

Moisture of the grain was measured at harvest.

Source

L. Darrah, R. Lundquist, D. West, C. Poneleit, B. Barry, B. Zehr, A. Bockholt, L. Maddux, K. Ziegler, and P. Martin. (1996). White Food Corn 1996 Performance Tests. Agricultural Research Service Special Report 502.

Examples

## Not run: 
  
  library(agridat)
  
  data(ars.earlywhitecorn96)
  dat <- ars.earlywhitecorn96

  libs(lattice)
  # These views emphasize differences between locations
  dotplot(gen~yield, dat, group=loc, auto.key=list(columns=3),
          main="ars.earlywhitecorn96")
  ## dotplot(gen~stalklodge, dat, group=loc, auto.key=list(columns=3),
  ##         main="ars.earlywhitecorn96")
  splom(~dat[,3:9], group=dat$loc, auto.key=list(columns=3),
        main="ars.earlywhitecorn96")
  
  # MANOVA
  m1 <- manova(cbind(yield,earht,moisture) ~ gen + loc, dat)
  m1
  summary(m1)
  

## End(Not run)

Multi-environment trial of soybean in Australia

Description

Yield and other traits of 58 varieties of soybeans, grown in four locations across two years in Australia. This is four-way data of Year x Loc x Gen x Trait.

Format

A data frame with 464 observations on the following 10 variables.

env

environment, 8 levels, first character of location and last two characters of year

loc

location

year

year

gen

genotype of soybeans, 1-58

yield

yield, metric tons / hectare

height

height (meters)

lodging

lodging

size

seed size, (millimeters)

protein

protein (percentage)

oil

oil (percentage)

Details

Measurement are available from four locations in Queensland, Australia in two consecutive years 1970, 1971.

The 58 different genotypes of soybeans consisted of 43 lines (40 local Australian selections from a cross, their two parents, and one other which was used a parent in earlier trials) and 15 other lines of which 12 were from the US.

Lines 1-40 were local Australian selections from Mamloxi (CPI 172) and Avoyelles (CPI 15939).

Note on the data in Basford and Tukey book. The values for line 58 for Nambour 1970 and Redland Bay 1971 are incorrectly listed on page 477 as 20.490 and 15.070. They should be 17.350 and 13.000, respectively. In the data set made available here, these values have been corrected.

Used with permission of Kaye Basford, Pieter Kroonenberg.

Source

Basford, K. E., and Tukey, J. W. (1999). Graphical analysis of multiresponse data illustrated with a plant breeding trial. Chapman and Hall/CRC.

Retrieved from: https://three-mode.leidenuniv.nl/data/soybeaninf.htm

References

K E Basford (1982). The Use of Multidimensional Scaling in Analysing Multi-Attribute Genotype Response Across Environments, Aust J Agric Res, 33, 473–480.

Kroonenberg, P. M., & Basford, K. E. B. (1989). An investigation of multi-attribute genotype response across environments using three-mode principal component analysis. Euphytica, 44, 109–123.

Marcin Kozak (2010). Use of parallel coordinate plots in multi-response selection of interesting genotypes. Communications in Biometry and Crop Science, 5, 83-95.

Examples

## Not run: 

  library(agridat)
  data(australia.soybean)
  dat <- australia.soybean

  libs(reshape2)
  dm <- melt(dat, id.var=c('env', 'year','loc','gen'))

  # Joint plot of genotypes & traits. Similar to Figure 1 of Kroonenberg 1989
  dmat <- acast(dm, gen~variable, fun=mean)
  dmat <- scale(dmat)
  biplot(princomp(dmat), main="australia.soybean trait x gen biplot", cex=.75)


  # Figure 1 of Kozak 2010, lines 44-58
  libs(reshape2, lattice, latticeExtra)
  data(australia.soybean)
  dat <- australia.soybean
  dat <- melt(dat, id.var=c('env', 'year','loc','gen'))
  dat <- acast(dat, gen~variable, fun=mean)
  dat <- scale(dat)
  dat <- as.data.frame(dat)[,c(2:6,1)]
  dat$gen <- rownames(dat)
  # data for the graphic by Kozak
  dat2 <- dat[44:58,]
  dat3 <- subset(dat2, is.element(gen, c("G48","G49","G50","G51")))

  parallelplot( ~ dat3[,1:6]|dat3$gen, main="australia.soybean",
               as.table=TRUE, horiz=FALSE) +
    parallelplot( ~ dat2[,1:6], horiz=FALSE, col="gray80") +
    parallelplot( ~ dat3[,1:6]|dat3$gen,
                 as.table=TRUE, horiz=FALSE, lwd=2)


## End(Not run)

Trial of wheat with nitrogen fertilizer in two fertility zones

Description

Trial of wheat with nitrogen fertilizer in two fertility zones

Usage

data("bachmaier.nitrogen")

Format

A data frame with 88 observations on the following 3 variables.

nitro

nitrogen fertilizer, kg/ha

yield

wheat yield, Mg/ha

zone

fertility zone

Details

Data from a wheat fertilizer experiment in Germany in two yield zones. In each zone, the design was an RCB with 4 blocks and 11 nitrogen levels. The yield of each plot was measured.

Electronic data originally downloaded from http://www.tec.wzw.tum.de/bachmaier/vino.zip (no longer available).

Source

Bachmaier, Martin. 2009. A Confidence Set for That X-Coordinate Where a Quadratic Regression Model Has a Given Gradient. Statistical Papers 50: 649–60. https://doi.org/10.1007/s00362-007-0104-1.

References

Bachmaier, Martin. Test and confidence set for the difference of the x-coordinates of the vertices of two quadratic regression models. Stat Papers (2010) 51:285–296, https://doi.org/10.1007/s00362-008-0159-7

Examples

library(agridat)
data(bachmaier.nitrogen)
dat <- bachmaier.nitrogen

# Fit a quadratic model for the low-fertility zone
dlow <- subset(dat, zone=="low")
m1 <- lm(yield ~ nitro + I(nitro^2), dlow)

# Slope of tangent line for economic optimum
m <- .005454 # = (N 0.60 euro/kg) / (wheat 110 euro/Mg)
# x-value of tangent point
b1 <- coef(m1)[2]
b2 <- coef(m1)[3]
opt.bach <- (m-b1)/(2*b2)
round(opt.bach, 0)

# conf int for x value of tangent point
round(vcovs <- vcov(m1), 7)
b1b1 <- vcovs[2,2] # estimated var of b1
b1b2 <- vcovs[2,3] # estimated cov of b1,b2
b2b2 <- vcovs[3,3]
tval <- qt(1 - 0.05/2, nrow(dlow)-3)
A <- b2^2 - b2b2 * tval^2
B <- (b1-m)*b2 - b1b2 * tval^2
C <- ((b1-m)^2 - b1b1 * tval^2)/4
D <- B^2 - 4*A*C
x.lo <- -2*C / (B-sqrt(B^2-4*A*C))
x.hi <- (-B + sqrt(B^2-4*A*C))/(2*A)
ci.bach <- c(x.lo, x.hi)
round(ci.bach,0) # 95% CI 173,260 Matches Bachmaier

# Plot raw data, fitted quadratic, optimum, conf int
plot(yield~nitro, dlow)
p1 <- data.frame(nitro=seq(0,260, by=1))
p1$pred <- predict(m1, new=p1)
lines(pred~nitro, p1)
abline(v=opt.bach, col="blue")
abline(v=ci.bach, col="skyblue")
title("Economic optimum with 95 pct confidence interval")

Uniformity trial of cotton in Egypt

Description

Uniformity trial of cotton in Egypt 1921-1923.

Usage

data("bailey.cotton.uniformity")

Format

A data frame with 794 observations on the following 5 variables.

row

row ordinate

col

column ordinate

yield

yield, in rotls

year

year

loc

location

Details

Two pickings were taken. The weights of seeds cotton for first and second pickings were totaled. Yields were measured in "rotl", which "are on the order of a pound".

Layout at Sakha and Gemmeiza (page 9): Total area 4.86 feddans. Each bed was 20 ridges of 7 m each, total dimension 15 m x 7 m. Add 1.5m for irrigation channel. Center-to-center distances 15m x 8.5m.

Charts 3 & 5 show yield of "Selected Average Plants". These data are not used here.

Chart 1: Sakha 1921, 8 x 20. Bed yield in rotls. Length 20 ridges * .75 m = 15m. Width = 7m.

Chart 2: Gemmeiza 1921, 8 x 20.

Chart 3: Total S.A.P. yield in grams. (not used here)

Chart 4: Gemmeiza 1922, 8 x 20.

Chart 5: Total S.A.P. yield in grams. (not used here)

Layout at Giza (page 10)

Beds were 8 ridges of 7 m each, total dimension 6m x 7m. Add 1.5m for irrigation channel. Center-to-center distance 6m x 8.5m

Chart 6 - Giza 1921, 14 x 11 = 154 plots

Chart 7 - Giza 1923, 20 x 8 = 160 plots

Bailey said the results at Giza 1921 were not suitable for reliability experiments.

Data were typed and proofread by KW 2023.01.11

Source

Bailey, M. A., and Trought, T. (1926). An account of experiments carried out to determine the experimental error of field trials with cotton in Egypt. Egypt Ministry of Agriculture, Technical and Science Service Bulletin 63, Min. Agriculture Egypt Technical and Science Bulletin 63. https://www.google.com/books/edition/Bulletin/xBQlAQAAIAAJ?pg=PA46-IA205

References

None

Examples

## Not run: 
  library(agridat)
  data(bailey.cotton.uniformity)
  dat <- bailey.cotton.uniformity
  dat <- transform(dat, env=paste(year,loc))

  # Data check. Matches Bailey 1926 Table 1. 28.13, , 46.02, 31.74, 13.52
  libs(dplyr)
  # dat 

  libs(desplot)
  desplot(dat, yield ~ col*row|env, main="bailey.cotton.uniformity")

  # The yield scales are quite different at each loc, and the dimensions
  # are different, so plot each location separately.
  # Note: Bailey does not say if plots are 7x15 meters, or 15x7 meters.
  # The choices here seem most likely in our opinion.
  desplot(dat, yield ~ col*row, subset= env=="1921 Sakha",
    main="1921 Sakha", aspect=(20*8.5)/(8*15))
  desplot(dat, yield ~ col*row, subset= env=="1921 Gemmeiza",
    main="1921 Gemmeiza", aspect=(20*8.5)/(8*15))
  desplot(dat, yield ~ col*row, subset= env=="1922 Gemmeiza",
    main="1922 Gemmeiza", aspect=(20*8.5)/(8*15))
  desplot(dat, yield ~ col*row, subset= env=="1921 Giza",
    main="1921 Giza", aspect=(11*6)/(14*8.5))
  # 1923 Giza has alternately hi/lo yield rows. Not noticed by Bailey.
  desplot(dat, yield ~ col*row, subset= env=="1923 Giza",
    main="1923 Giza", aspect=(20*6)/(8*8.5))
  

## End(Not run)  

Uniformity trials of barley, 10 years on same ground

Description

Uniformity trials of barley at Davis, California, 1925-1935, 10 years on same ground.

Format

A data frame with 570 observations on the following 4 variables.

row

row

col

column

year

year

yield

yield, pounds/acre

Details

Ten years of uniformity trials were sown on the same ground. Baker (1952) shows a map of the field, in which gravel subsoil extended from the upper right corner diagonally lower-center. This part of the field had lower yields on the 10-year average map.

Plot 41 in 1928 is missing.

Field width: 19 plots = 827 ft

Field length: 3 plots * 161 ft + 2 alleys * 15 feet = 513 ft

Source

Baker, GA and Huberty, MR and Veihmeyer, FJ. (1952) A uniformity trial on unirrigated barley of ten years' duration. Agronomy Journal, 44, 267-270. https://doi.org/10.2134/agronj1952.00021962004400050011x

Examples

## Not run: 

library(agridat)

data(baker.barley.uniformity)
dat <- baker.barley.uniformity

# Ten-year average
dat2 <- aggregate(yield ~ row*col, data=dat, FUN=mean, na.rm=TRUE)

libs(desplot)
desplot(dat, yield~col*row|year,
        aspect = 513/827, # true aspect
        main="baker.barley.uniformity - heatmaps by year")

desplot(dat2, yield~col*row,
        aspect = 513/827, # true aspect
        main="baker.barley.uniformity - heatmap of 10-year average")
# Note low yield in upper right, slanting to left a bit due to sandy soil
# as shown in Baker figure 1.


# Baker fig 2, stdev vs mean
dat3 <- aggregate(yield ~ row*col, data=dat, FUN=sd, na.rm=TRUE)
plot(dat2$yield, dat3$yield, xlab="Mean yield", ylab="Std Dev yield",
     main="baker.barley.uniformity")

# Baker table 4, correlation of plots across years
# libs(reshape2)
# mat <- acast(dat, row+col~year)
# round(cor(mat, use='pair'),2)


## End(Not run)

Uniformity trial of strawberry

Description

Uniformity trial of strawberry

Usage

data("baker.strawberry.uniformity")

Format

A data frame with 700 observations on the following 4 variables.

trial

Factor for trial

row

row ordinate

col

column ordinate

yield

yield per plant/plot in grams

Details

Trial T1:

200 plants were grown in two double-row beds at Davis, California, in 1946. The rows were 1 foot apart. The beds were 42 inches apart. The plants were 10 inches apart within a row, each row consisting of 50 plants.

Field length: 50 plants * 10 inches = 500 inches.

Field width: 12 in + 42 in + 12 in = 66 inches.

The layout of the experiment in Table 1 shows 4 columns. There is 12 inches between column 1 and column 2, then 42 inches, then 12 inches between column 3 and column 4. For the data in this R package, we added 3 to the right two columns index values to indicate this layout. (Should be 3.5, but we want to have an integer).

Trial T2:

500 plants were grown in single beds. The beds were 30 inches apart. Each bed was 50 plants long with 10 inches between plants.

Field length: 50 plants * 10 in = 500 in.

Field width: 10 beds * 30 in = 300 in.

Source

G. A. Baker and R. E. Baker (1953). Strawberry Uniformity Yield Trials. Biometrics, 9, 412-421. https://doi.org/10.2307/3001713

References

None

Examples

## Not run: 

library(agridat)

data(baker.strawberry.uniformity)
dat <- baker.strawberry.uniformity

# Match mean and cv of Baker p 414.
libs(dplyr)
dat <- group_by(dat, trial)
summarize(dat, mn=mean(yield), cv=sd(yield)/mean(yield))

libs(desplot)
desplot(dat, yield ~ col*row, subset=trial=="T1",
        flip=TRUE, aspect=500/66, tick=TRUE,
        main="baker.strawberry.uniformity - trial T1")
desplot(dat, yield ~ col*row, subset=trial=="T2",
        flip=TRUE, aspect=500/300, tick=TRUE,
        main="baker.strawberry.uniformity - trial T2")


## End(Not run)

Uniformity trial of wheat

Description

Uniformity trial of wheat

Usage

data("baker.wheat.uniformity")

Format

A data frame with 225 observations on the following 3 variables.

row

row

col

col

yield

yield (grams)

Details

Data was collected in 1939-1940. The trial consists of sixteen 40 ft. x 40 ft. blocks subdivided into nine plots each. The data were secured in 1939-1940 from White Federation wheat. The design of the experiment was square with alleys 20 feet wide between blocks. The plots were 10 feet long with two guard rows on each side.

Morning glories infested the middle two columns of blocks, uniformly over the blocks affected.

The data here include missing values for the alleys so that the field map is approximately the correct shape and size.

Field width: 4 blocks of 40 feet + 3 alleys of 20 feet = 220 feet.

Field length: 4 blocks of 40 feet + 3 alleys of 20 feet = 220 feet.

Source

G. A. Baker, E. B. Roessler (1957). Implications of a uniformity trial with small plots of wheat. Hilgardia, 27, 183-188. https://hilgardia.ucanr.edu/Abstract/?a=hilg.v27n05p183 https://doi.org/10.3733/hilg.v27n05p183

References

None

Examples

## Not run: 
  
  library(agridat)
  data(baker.wheat.uniformity)
  dat <- baker.wheat.uniformity

  libs(desplot)
  desplot(dat, yield ~ col*row,
          flip=TRUE, aspect=1,
          main="baker.wheat.uniformity")


## End(Not run)

Uniformity trial of peanuts

Description

Uniformity trial of peanuts in Alabama, 1946.

Usage

data("bancroft.peanut.uniformity")

Format

A data frame with 216 observations on the following 5 variables.

row

row

col

column

yield

yield, pounds per plot

block

block

Details

The data are obtained from two parts of the same field, located at Wiregrass Substation, Headland, Alabama, USA. Each part had 18 rows, 3 feet wide, 100 feet long. Plots were harvested in 1946. Green weights in pounds were recorded.

Each plot was 16.66 linear feet of row and 3 feet in width, 50 sq feet.

Field width: 6 plots * 16.66 feet = 100 feet

Field length: 18 plots * 3 feet = 54 feet

Conclusions: Based on the relative efficiencies, increasing the size of the plot along the row is better than across the row. Narrow, rectangular plots are more efficient.

Source

Bancroft, T. A. et a1., (1948). Size and Shape of Plots and Distribution of Plot Yield for Field Experiments with Peanuts. Alabama Agricultural Experiment Station Progress Report, sec. 39. Table 4, page 6. https://aurora.auburn.edu/bitstream/handle/11200/1345/0477PROG.pdf;sequence=1

References

None

Examples

## Not run: 

library(agridat)
data(bancroft.peanut.uniformity)
dat <- bancroft.peanut.uniformity
  
# match means Bancroft page 3
## dat 
## # A tibble: 2 x 2
##   block    mn
##   <chr> <dbl>
## 1 B1     2.46
## 2 B2     2.05
  
libs(desplot)
desplot(dat, yield ~ col*row|block,
        flip=TRUE, aspect=(18*3)/(6*16.66), # true aspect
        main="bancroft.peanut.uniformity")


## End(Not run)  

Multi-environment trial of maize in Texas.

Description

Multi-environment trial of maize in Texas.

Usage

data("barrero.maize")

Format

A data frame with 14568 observations on the following 15 variables.

year

year of testing, 2000-2010

yor

year of release, 2000-2010

loc

location, 16 places in Texas

env

environment (year+loc), 107 levels

rep

replicate, 1-4

gen

genotype, 847 levels

daystoflower

numeric

plantheight

plant height, cm

earheight

ear height, cm

population

plants per hectare

lodged

percent of plants lodged

moisture

moisture percent

testweight

test weight kg/ha

yield

yield, Mt/ha

Details

This is a large (14500 records), multi-year, multi-location, 10-trait dataset from the Texas AgriLife Corn Performance Trials.

These data are from 2-row plots approximately 36in wide by 25 feet long.

Barrero et al. used this data to estimate the genetic gain in maize hybrids over a 10-year period of time.

Used with permission of Seth Murray.

Source

Barrero, Ivan D. et al. (2013). A multi-environment trial analysis shows slight grain yield improvement in Texas commercial maize. Field Crops Research, 149, Pages 167-176. https://doi.org/10.1016/j.fcr.2013.04.017

References

None.

Examples

## Not run: 
  library(agridat)
  data(barrero.maize)
  dat <- barrero.maize

  library(lattice)
  bwplot(yield ~ factor(year)|loc, dat,
         main="barrero.maize - Yield trends by loc",
         scales=list(x=list(rot=90)))
  
  # Table 6 of Barrero. Model equation 1.
  if(require("asreml", quietly=TRUE)){
    libs(dplyr,lucid)
    dat <- arrange(dat, env)
    dat <- mutate(dat,
                  yearf=factor(year), env=factor(env),
                  loc=factor(loc), gen=factor(gen), rep=factor(rep))
  
    m1 <- asreml(yield ~ loc + yearf + loc:yearf, data=dat,
                 random = ~ gen + rep:loc:yearf +
                   gen:yearf + gen:loc +
                   gen:loc:yearf,
                 residual = ~ dsum( ~ units|env),
                 workspace="500mb")
  
    # Variance components for yield match Barrero table 6.
    lucid::vc(m1)[1:5,]
    ##        effect component std.error z.ratio bound 
    ## rep:loc:yearf   0.111     0.01092    10       P 0  
    ##           gen   0.505     0.03988    13       P 0  
    ##     gen:yearf   0.05157   0.01472     3.5     P 0  
    ##       gen:loc   0.02283   0.0152      1.5     P 0.2
    ## gen:loc:yearf   0.2068    0.01806    11       P 0  
    
    summary(vc(m1)[6:112,"component"]) # Means match last row of table 6
    ##   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    ## 0.1286  0.3577  0.5571  0.8330  1.0322  2.9867 
  }

## End(Not run)

Uniformity trials of apples, lemons, oranges, and walnuts

Description

Uniformity trials of apples, lemons, oranges, and walnuts, in California & Utah, 1915-1918.

Format

Each dataset has the following format

row

row

col

column

yield

yield per tree in pounds

Details

A few of the trees affected by disease were eliminated and the yield was replaced by the average of the eight surrounding trees.

The following details are from Batchelor (1918).

Jonathan Apples

"The apple records were obtained from a 10-year old Jonathan apple orchard located at Providence, Utah. The surface soil of this orchard is very uniform to all appearances except on the extreme eastern edge, where the percentage of gravel increases slightly. The trees are planted 16 feet apart, east and west, and 30 feet apart north and south."

Note: The orientation of the field is not given in the paper, but all other fields in the paper have north at the top, so that is assumed to be true for this field as well. Yields may be from 1916.

Field width: 8 trees * 16 feet = 128 feet

Field length: 28 rows * 30 feet = 840 feet

Eureka Lemon

The lemon (Citrus limonia) tree yields were obtained from a grove of 364 23-year-old trees, located at Upland, California. The records extend from October 1, 1915, to October 1, 1916. The grove consists of 14 rows of 23-year-old trees, extending north and south, with 26 trees in a row, planted 24 by 24 feet apart. This grove presents the most uniform appearance of any under consideration [in this paper]. The land is practically level, and the soil is apparently uniform in texture. The records show a grouping of several low-yielding trees; yet a field observation gives one the impression that the grove as a whole is remarkably uniform.

Field width: 14 trees * 24 feet = 336 feet

Field length: 26 trees * 24 feet = 624 feet

Navel 1 at Arlington

These records were of the 1915-16 yields of one thousand 24-year-old navel-orange trees near Arlington station, Riverside, California. The grove consists of 20 rows of trees from north to south, with 50 trees in a row, planted 22 by 22 feet. A study of the records shows certain distinct high- and low-yielding areas. The northeast corner and the south end contain notably high-yielding trees. The north two-thirds of the west side contains a large number of low-yielding trees. These areas are apparently correlated with soil variation. Variations from tree to tree also occur, the cause of which is not evident. These variations, which are present in every orchard, bring uncertainty into the results offield experiments.

Field width: 20 trees * 22 feet = 440 feet

Field length: 50 trees * 22 feet = 1100 feet

Navel 2 at Antelope

The navel-orange grove later referred to as the Antelope Heights navels is a plantation of 480 ten-yearold trees planted 22 by 22 feet, located at Naranjo, California. The yields are from 1916. The general appearance of the trees gives a visual impression of uniformity greater than a comparison of the individual tree production substantiates.

Field width: 15 trees * 22 feet = 330 feet

Field length: 33 trees * 22 feet = 726 feet

Valencia Orange

The Valencia orange grove is composed of 240 15-year-old trees, planted 21 feet 6 inches by 22 feet 6 inches, located at Villa Park, California. The yields were obtained in 1916.

Field width: 12 rows * 22 feet = 264 feet

Field length: 20 rows * 22 feet = 440 feet

Walnut

The walnut (Juglans regia) yields were obtained during the seasons of 1915 and 1916 from a 24-year-old Santa Barbara softshell seedling grove, located at Whittier, California. [Note, The yields here appear to be the 1915 yields.] The planting is laid out 10 trees wide and 32 trees long, entirely surrounded by additional walnut plantings, except on a part of one side which is adjacent to an orange grove. The trees are planted on the square system, 50 feet apart.

Field width: 10 trees * 50 feet = 500 feet

Field length: 32 trees * 50 feet = 1600 feet

Source

L. D. Batchelor and H. S. Reed. (1918). Relation of the variability of yields of fruit trees to the accuracy of field trials. J. Agric. Res, 12, 245–283. https://books.google.com/books?id=Lil6AAAAMAAJ&lr&pg=PA245

References

McCullagh, P. and Clifford, D., (2006). Evidence for conformal invariance of crop yields, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science, 462, 2119–2143. https://doi.org/10.1098/rspa.2006.1667

Examples

## Not run: 

library(agridat)
libs(desplot)

  # Apple
  data(batchelor.apple.uniformity)
  desplot(batchelor.apple.uniformity, yield~col*row,
          aspect=840/128, tick=TRUE, # true aspect
          main="batchelor.apple.uniformity")

  # Lemon
  data(batchelor.lemon.uniformity)
  desplot(batchelor.lemon.uniformity, yield~col*row,
          aspect=624/336, # true aspect
          main="batchelor.lemon.uniformity")

  # Navel1 (Arlington)
  data(batchelor.navel1.uniformity)
  desplot(batchelor.navel1.uniformity, yield~col*row,
          aspect=1100/440, # true aspect
          main="batchelor.navel1.uniformity - Arlington")

  # Navel2 (Antelope)
  data(batchelor.navel2.uniformity)
  desplot(batchelor.navel2.uniformity, yield~col*row,
          aspect=726/330, # true aspect
          main="batchelor.navel2.uniformity - Antelope")

  # Valencia
  data(batchelor.valencia.uniformity)
  desplot(batchelor.valencia.uniformity, yield~col*row,
          aspect=440/264, # true aspect
          main="batchelor.valencia.uniformity")

  # Walnut
  data(batchelor.walnut.uniformity)
  desplot(batchelor.walnut.uniformity, yield~col*row,
          aspect=1600/500, # true aspect
          main="batchelor.walnut.uniformity")


## End(Not run)

Survey and satellite data for corn and soy areas in Iowa

Description

Survey and satellite data for corn and soy areas in Iowa

Usage

data("battese.survey")

Format

A data frame with 37 observations on the following 9 variables.

county

county name

segment

sample segment number (within county)

countysegs

number of segments in county

cornhect

hectares of corn in segment

soyhect

hectares of soy

cornpix

pixels of corn in segment

soypix

pixels of soy

cornmean

county mean of corn pixels per segment

soymean

county mean of soy pixels per segment

Details

The data are for 12 counties in north-central Iowa in 1978.

The USDA determined the area of soybeans in 37 area sampling units (called 'segments'). Each segment is about one square mile (about 259 hectares). The number of pixels of that were classified as corn and soybeans came from Landsat images obtained in Aug/Sep 1978. Each pixel represents approximately 0.45 hectares.

Data originally compiled by USDA.

This data is also available in R packages: 'rsae::landsat' and 'JoSAE::landsat'.

Source

Battese, George E and Harter, Rachel M and Fuller, Wayne A. (1988). An error-components model for prediction of county crop areas using survey and satellite data. Journal of the American Statistical Association, 83, 28-36. https://doi.org/10.2307/2288915

Battese (1982) preprint version. https://www.une.edu.au/__data/assets/pdf_file/0017/15542/emetwp15.pdf

References

Pushpal K Mukhopadhyay and Allen McDowell. (2011). Small Area Estimation for Survey Data Analysis Using SAS Software SAS Global Forum 2011.

Examples

## Not run: 

library(agridat)
data(battese.survey)
dat <- battese.survey

# Battese fig 1 & 2.  Corn plot shows outlier in Hardin county
libs(lattice)
dat <- dat[order(dat$cornpix),]
xyplot(cornhect ~ cornpix, data=dat, group=county, type=c('p','l'),
       main="battese.survey", xlab="Pixels of corn", ylab="Hectares of corn",
       auto.key=list(columns=3))

dat <- dat[order(dat$soypix),]
xyplot(soyhect ~ soypix, data=dat, group=county, type=c('p','l'),
       main="battese.survey", xlab="Pixels of soy", ylab="Hectares of soy",
       auto.key=list(columns=3))

libs(lme4, lucid)
  
# Fit the models of Battese 1982, p.18.  Results match
m1 <- lmer(cornhect ~ 1 + cornpix + (1|county), data=dat)
fixef(m1)
## (Intercept)     cornpix 
##   5.4661899   0.3878358 
vc(m1)
##      grp        var1 var2   vcov  sdcor
##   county (Intercept) <NA>  62.83  7.926
## Residual        <NA> <NA> 290.4  17.04 
m2 <- lmer(soyhect ~ 1 + soypix + (1|county), data=dat)
fixef(m2)
## (Intercept)      soypix 
##  -3.8223566   0.4756781 
vc(m2)
##      grp        var1 var2  vcov sdcor
##   county (Intercept) <NA> 239.2 15.47
## Residual        <NA> <NA> 180   13.42
  
# Predict for Humboldt county as in Battese 1982 table 2
5.4662+.3878*290.74
# 118.2152 # mu_i^0
5.4662+.3878*290.74+ -2.8744
# 115.3408 # mu_i^gamma
(185.35+116.43)/2
# 150.89 # y_i bar
  
# Survey regression estimator of Battese 1988
  
# Delete the outlier
dat2 <- subset(dat, !(county=="Hardin" & soyhect < 30))
  
# Results match top-right of Battese 1988, p. 33
m3 <- lmer(cornhect ~ cornpix + soypix + (1|county), data=dat2)
fixef(m3)
## (Intercept)     cornpix      soypix 
##  51.0703979   0.3287217  -0.1345684 
vc(m3)
##      grp        var1 var2  vcov sdcor
##   county (Intercept) <NA> 140   11.83
## Residual        <NA> <NA> 147.3 12.14
m4 <- lmer(soyhect ~ cornpix + soypix + (1|county), data=dat2)
fixef(m4)
##  (Intercept)      cornpix       soypix 
## -15.59027098   0.02717639   0.49439320 
vc(m4)
##      grp        var1 var2  vcov sdcor
##   county (Intercept) <NA> 247.5 15.73
## Residual        <NA> <NA> 190.5 13.8 


## End(Not run)

Counts of webworms in a beet field, with insecticide treatments.

Description

Counts of webworms in a beet field, with insecticide treatments.

Usage

data("beall.webworms")

Format

A data frame with 1300 observations on the following 7 variables.

row

row

col

column

y

count of webworms

block

block

trt

treatment

spray

spray treatment yes/no

lead

lead treatment yes/no

Details

The beet webworm lays egg masses as small as 1 egg, seldom exceeding 5 eggs. The larvae can move freely, but usually mature on the plant on which they hatch.

Each plot contained 25 unit areas, each 1 row by 3 feet long. The row width is 22 inches. The arrangement of plots within the blocks seems certain, but the arrangement of the blocks/treatments is not certain, since the authors say "since the plots were 5 units long and 5 wide it is only practicable to combine them into groups of 5 in one direction or the other".

Treatment 1 = None. Treatment 2 = Contact spray. Treatment 3 = Lead arsenate. Treatment 4 = Both spray, lead arsenate.

Source

Beall, Geoffrey (1940). The fit and significance of contagious distributions when applied to observations on larval insects. Ecology, 21, 460-474. Table 6. https://doi.org/10.2307/1930285

References

Michal Kosma et al. (2019). Over-dispersed count data in crop and agronomy research. Journal of Agronomy and Crop Science. https://doi.org/10.1111/jac.12333

Examples

## Not run: 

library(agridat)
data(beall.webworms)
dat <- beall.webworms

# Match Beall table 1
# with(dat, table(y,trt))

libs(lattice)
histogram(~y|trt, data=dat, layout=c(1,4), as.table=TRUE,
          main="beall.webworms")

# Visualize Beall table 6.  Block effects may exist, but barely.
libs(desplot)
grays <- colorRampPalette(c("white","#252525"))
desplot(dat, y ~ col*row,
        col.regions=grays(10),
        at=0:10-0.5,
        out1=block, out2=trt, num=trt, flip=TRUE, # aspect unknown
        main="beall.webworms (count of worms)")

# Following plot suggests interaction is needed
# with(dat, interaction.plot(spray, lead, y))

# Try the models of Kosma et al, Table 1.

# Poisson model
m1 <- glm(y ~ block + spray*lead, data=dat, family="poisson")
logLik(m1) # -1497.719 (df=16)

# Negative binomial model
# libs(MASS)
# m2 <- glm.nb(y ~ block + spray*lead, data=dat)
# logLik(m2) # -1478.341 (df=17)

# # Conway=Maxwell-Poisson model (takes several minutes)
# libs(spaMM)
# # estimate nu parameter
# m3 <- fitme(y ~ block + spray*lead, data=dat, family = COMPoisson())
# logLik(m3) # -1475.999 
# # Kosma logLik(m3)=-1717 seems too big. Typo? Different model?


## End(Not run)

Yields of 8 barley varieties in 1913 as used by Student.

Description

Yields of 8 barley varieties in 1913.

Usage

data("beaven.barley")

Format

A data frame with 160 observations on the following 4 variables.

row

row

col

column

gen

genotype

yield

yield (grams)

Details

Eight races of barley were grown on a regular pattern of plots.

These data were prepared from Richey (1926) because the text was cleaner.

Each plot was planted 40 inches on a side, but only the middle square 36 inches on a side was harvested.

Field width: 32 plots * 3 feet = 96 feet

Field length: 5 plots * 3 feet = 15 feet

Source

Student. (1923). On testing varieties of cereals. Biometrika, 271-293.

https://doi.org/10.1093/biomet/15.3-4.271

References

Frederick D. Richey (1926). The moving average as a basis for measuring correlated variation in agronomic experiments. Jour. Agr. Research, 32, 1161-1175.

Examples

## Not run: 

library(agridat)

data(beaven.barley)
dat <- beaven.barley

# Match the means shown in Richey table IV
tapply(dat$yield, dat$gen, mean)
##       a       b       c       d       e       f       g       h
## 298.080 300.710 318.685 295.260 306.410 276.475 304.605 271.820

# Compare to Student 1923, diagram I,II
libs(desplot)
desplot(dat, yield ~ col*row,
        aspect=15/96, # true aspect
        main="beaven.barley - variety trial", text=gen)

## End(Not run)

Mating crosses of chickens

Description

Mating crosses of chickens

Usage

data("becker.chicken")

Format

A data frame with 45 observations on the following 3 variables.

male

male parent

female

female parent

weight

weight (g) at 8 weeks

Details

From a large flock White Rock chickens, five male sires were chosen and mated to each of three female dams, producing 3 female progeny. The data are body weights at eight weeks of age.

Becker (1984) used these data to demonstrate the calculation of heritability.

Source

Walter A. Becker (1984). Manual of Quantitative Genetics, 4th ed. Page 83.

References

None

Examples

## Not run: 

  library(agridat)
  data(becker.chicken)
  dat <- becker.chicken
  
  libs(lattice)
  dotplot(weight ~ female, data=dat, group=male,
          main="becker.chicken - progeny weight by M*F",
          xlab="female parent",ylab="progeny weight",
          auto.key=list(columns=5))

  # Sums match Becker
  # sum(dat$weight)
  # aggregate(weight ~  male + female, dat, FUN=sum)

  # Variance components
  libs(lme4,lucid)
  m1 <- lmer(weight ~  (1|male) + (1|female), data=dat)
  # vc(m1)
  ## grp        var1 var2      vcov    sdcor
  ## 1   female (Intercept) <NA> 1096   33.1
  ## 2     male (Intercept) <NA>  776.8 27.87
  ## 3 Residual        <NA> <NA> 5524   74.32

  # Calculate heritabilities
  # s2m <- 776  # variability for males
  # s2f <- 1095 # variability for females
  # s2w <- 5524 # variability within crosses
  # vp <- s2m + s2f + s2w # 7395
  # 4*s2m/vp # .42 male heritability
  #4*s2f/vp # .59 female heritability


## End(Not run)

Multi-environment trial of wheat with Augmented design

Description

Multi-environment trial of wheat in Nebraska with Augmented design

Usage

data("belamkar.augmented")

Format

A data frame with 2700 observations on the following 9 variables.

loc

location

rep

replicate

iblock

incomplete block

gen_new

new genotype (1=yes, 0=no)

gen_check

check genotype (0=no)

gen

genotype name

col

column ordinate

row

row ordinate

yield

yield, bu/ac

Details

The experiment had 8 locations with 270 new, experimental lines (genotypes) and 3 check lines. There were 10 incomplete blocks at each location. There were 2 replicate blocks at Alliance and 1 block at all other locations. Each plot was 3 m long by 1.2 m wide.

The electronic data were found in supplement S4 downloaded from https://doi.org/10.25387/g3.6249410 The license for the data is CC-BY 4.0.

Source

Vikas Belamkar, Mary J. Guttieri, Waseem Hussain, Diego Jarquín, Ibrahim El-basyoni, Jesse Poland, Aaron J. Lorenz, P. Stephen Baenziger (2018). Genomic Selection in Preliminary Yield Trials in a Winter Wheat Breeding Program. G3 Genes|Genomes|Genetics, 8, Pages 2735–2747. https://doi.org/10.1534/g3.118.200415

References

Same data appear in ASRtriala package: https://vsni.co.uk/free-software/asrtriala

Examples

## Not run: 
  library(agridat)
  data(belamkar.augmented)
  dat <- belamkar.augmented

  libs(desplot)
  desplot(dat, yield ~ col*row|loc, out1=rep, out2=iblock)
  # Experiment design showing check placement
  dat$gen_check <- factor(dat$gen_check)
  desplot(dat, gen_check ~ col*row|loc, out1=rep, out2=iblock,
          main="belamkar.augmented")

  # Belamkar supplement S3 has R code for analysis
  if(require("asreml", quietly=TRUE)){
    library(asreml)

    # AR1xAR1 model to calculate BLUEs for a single loc
    d1 <- droplevels(subset(dat, loc=="Lincoln"))
    d1$colf <- factor(d1$col)
    d1$rowf <- factor(d1$row)
    d1$gen <- factor(d1$gen)
    d1$gen_check <- factor(d1$gen_check)
    d1 <- d1[order(d1$col),]
    d1 <- as.data.frame(d1)
    m1 <- asreml(fixed=yield ~ gen_check, data=d1,
                 random = ~ gen_new:gen,
                 residual = ~ar1(colf):ar1v(rowf) )
    p1 <- predict(m1, classify="gen")
    head(p1$pvals)
  }

## End(Not run)

RCB experiment of spring barley in United Kingdom

Description

RCB experiment of spring barley in United Kingdom

Format

A data frame with 225 observations on the following 4 variables.

col

column (also blocking factor)

row

row

yield

yield

gen

variety/genotype

Details

RCB design, each column is one rep.

Used with permission of David Higdon.

Source

Besag, J. E., Green, P. J., Higdon, D. and Mengersen, K. (1995). Bayesian computation and stochastic systems. Statistical Science, 10, 3-66. https://www.jstor.org/stable/2246224

References

Davison, A. C. 2003. Statistical Models. Cambridge University Press. Pages 534-535.

Examples

## Not run: 
  
  library(agridat)
  data(besag.bayesian)
  dat <- besag.bayesian

  # Yield values were scaled to unit variance
  # var(dat$yield, na.rm=TRUE)
  # .999

  # Besag Fig 2. Reverse row numbers to match Besag, Davison
  dat$rrow <- 76 - dat$row
  libs(lattice)
  xyplot(yield ~ rrow|col, dat, layout=c(1,3), type='s',
         xlab="row", ylab="yield", main="besag.bayesian")

  if(require("asreml", quietly=TRUE)) {
    libs(asreml, lucid)

    # Use asreml to fit a model with AR1 gradient in rows  
    dat <- transform(dat, cf=factor(col), rf=factor(rrow))
    m1 <- asreml(yield ~ -1 + gen, data=dat, random= ~ rf)
    m1 <- update(m1, random= ~ ar1v(rf))
    m1 <- update(m1)
    m1 <- update(m1)
    m1 <- update(m1)
    lucid::vc(m1)
  
    # Visualize trends, similar to Besag figure 2.
    # Need 'as.vector' because asreml uses a named vector
    dat$res <- unname(m1$resid)
    dat$geneff <- coef(m1)$fixed[as.numeric(dat$gen)]
    dat <- transform(dat, fert=yield-geneff-res)
    libs(lattice)
    xyplot(geneff ~ rrow|col, dat, layout=c(1,3), type='s',
           main="besag.bayesian - Variety effects", ylim=c(5,15 ))
    xyplot(fert ~ rrow|col, dat, layout=c(1,3), type='s',
           main="besag.bayesian - Fertility", ylim=c(-2,2))
    xyplot(res ~ rrow|col, dat, layout=c(1,3), type='s',
           main="besag.bayesian - Residuals", ylim=c(-4,4))
  }

## End(Not run) 

Competition experiment in beans with height measurements

Description

Competition experiment in beans with height measurements

Usage

data("besag.beans")

Format

A data frame with 152 observations on the following 6 variables.

gen

genotype / variety

height

plot height, cm

yield

plot yield, g

row

row / block

rep

replicate factor

col

column

Details

Field beans of regular height were grown beside shorter varieties. In each block, each variety occurred once as a left-side neighbor and once as a right-side neighbor of every variety (including itself). Border plots were placed at the ends of each block. Each block with 38 adjacent plots. Each plot was one row, 3 meters long with 50 cm spacing between rows. No gaps between plots. Spacing between plants was 6.7 cm. Four blocks (rows) were used, each with six replicates.

Plot yield and height was recorded.

Kempton and Lockwood used models that adjusted yield according to the difference in height of neighboring plots.

Field length: 4 plots * 3m = 12m

Field width: 38 plots * 0.5 m = 19m

Source

Julian Besag and Rob Kempton (1986). Statistical Analysis of Field Experiments Using Neighbouring Plots. Biometrics, 42, 231-251. Table 6. https://doi.org/10.2307/2531047

References

Kempton, RA and Lockwood, G. (1984). Inter-plot competition in variety trials of field beans (Vicia faba L.). The Journal of Agricultural Science, 103, 293–302.

Examples

## Not run: 

library(agridat)

data(besag.beans)
dat = besag.beans

libs(desplot)
desplot(dat, yield ~ col*row,
        aspect=12/19, out1=row, out2=rep, num=gen, cex=1, # true aspect
        main="besag.beans")


libs(reshape2)
# Add a covariate = excess height of neighbors
mat <- acast(dat, row~col, value.var='height')
mat2 <- matrix(NA, nrow=4, ncol=38)
mat2[,2:37] <- (mat[,1:36] + mat[,3:38] - 2*mat[,2:37])
dat2 <- melt(mat2)
colnames(dat2) <- c('row','col','cov')
dat <- merge(dat, dat2)
  
# Drop border plots
dat <- subset(dat, rep != 'R0')
  
libs(lattice)
# Plot yield vs neighbors height advantage
xyplot(yield~cov, data=dat, group=gen,
       main="besag.beans",
       xlab="Mean excess heights of neighbor plots",
       auto.key=list(columns=3))
  
# Trial mean.
mean(dat$yield) # 391 matches Kempton table 3
  
# Mean excess height of neighbors for each genotype
# tapply(dat$cov, dat$gen, mean)/2 # Matches Kempton table 4

# Variety means, matches Kempton table 4 mean yield
m1 <- lm(yield ~ -1 + gen, dat)
coef(m1)

# Full model used by Kempton, eqn 5.  Not perfectly clear.
# Appears to include rep term, perhaps within block
dat$blk <- factor(dat$row)
dat$blkrep <- factor(paste(dat$blk, dat$rep))
m2 <- lm(yield ~ -1 + gen + blkrep + cov, data=dat)
coef(m2) # slope 'cov' = -.72, while Kempton says -.79


## End(Not run)

Check variety yields in winter wheat.

Description

Check variety yields in winter wheat.

Usage

data("besag.checks")

Format

A data frame with 364 observations on the following 4 variables.

yield

yield, units of 10g

row

row

col

column

gen

genotype/variety

Details

This data was used by Besag to show the spatial variation in a field experiment, but Besag did not use the data for any analysis.

Yields of winter wheat varieties (Bounty and Huntsman) at the Plant Breeding Institute, Cambridge, in 1980. These data are the 'checks' genotypes in a larger variety trial.

There is a column of checks, then five columns of new varieties. Repeat.

Plot dimensions approx 1.5 by 4.5 metres

Field length: 52 rows * 4.5 m = 234 m

Field width: 31 columns * 1.5 m = 46.5

Electronic version of data supplied by David Clifford.

Source

Besag, J.E. & Kempton R.A. (1986). Statistical analysis of field experiments using neighbouring plots. Biometrics, 42, 231-251. https://doi.org/10.2307/2531047

References

Kempton, Statistical Methods for Plant Variety Evaluation, page 91–92

Examples

## Not run: 
library(agridat)
data(besag.checks)
dat <- besag.checks
  libs(desplot)
  desplot(dat, yield~col*row,
          num=gen, aspect=234/46.5, # true aspect
          main="besag.checks")

## End(Not run)

RCB experiment of wheat, 50 varieties in 3 blocks with strong spatial trend.

Description

RCB experiment of wheat, 50 varieties in 3 blocks with strong spatial trend.

Format

A data frame with 150 observations on the following 4 variables.

yield

yield of wheat

gen

genotype, factor with 50 levels

col

column/block

row

row

Details

RCB experiment on wheat at El Batan, Mexico. There are three single-column replicates with 50 varieties in each replicate.

Plot dimensions are not given by Besag.

Data retrieved from https://web.archive.org/web/19991008143232/www.stat.duke.edu/~higdon/trials/elbatan.dat

Used with permission of David Higdon.

Source

Julian Besag and D Higdon, 1999. Bayesian Analysis of Agricultural Field Experiments, Journal of the Royal Statistical Society: Series B,61, 691–746. Table 1. https://doi.org/10.1111/1467-9868.00201

References

Wilkinson 1984.

Besag & Seheult 1989.

Examples

## Not run: 

library(agridat)
data(besag.elbatan)
dat <- besag.elbatan

libs(desplot)
desplot(dat, yield~col*row,
        num=gen, # aspect unknown
        main="besag.elbatan - wheat yields")


# Besag figure 1
library(lattice)
xyplot(yield~row|col, dat, type=c('l'),
       layout=c(1,3), main="besag.elbatan wheat yields")


# RCB
m1 <- lm(yield ~ 0 + gen + factor(col), dat)
p1 <- coef(m1)[1:50]

# Formerly used gam package, but as of R 3.1, Rcmd check --as-cran
# is complaining
# Calls: plot.gam ... model.matrix.gam -> predict -> predict.gam -> array
# but it works perfectly in interactive mode !!!
# Remove the FALSE to run the code below
if(is.element("gam", search())) detach(package:gam)
libs(mgcv)
m2 <- mgcv::gam(yield ~ -1 + gen + factor(col) + s(row), data=dat)
plot(m2, residuals=TRUE, main="besag.elbatan")
pred <- cbind(dat, predict(m2, dat, type="terms"))
# Need to correct for the average loess effect, which is like
# an overall intercept term.
adjlo <-  mean(pred$"s(row)")
p2 <- coef(m2)[1:50] + adjlo

# Compare estimates
lims <- range(c(p1,p2))
plot(p1, p2, xlab="RCB prediction",
     ylab="RCB with smooth trend (predicted)",
     type='n', xlim=lims, ylim=lims,
     main="besag.elbatan")
text(p1, p2, 1:50, cex=.5)
abline(0,1,col="gray")


## End(Not run)

Presence of footroot disease in an endive field

Description

Presence of footroot disease in an endive field

Format

A data frame with 2506 observations on the following 3 variables.

col

column

row

row

disease

plant is diseased, Y=yes,N=no

Details

In a field of endives, does each plant have footrot, or not? Data are binary on a lattice of 14 x 179 plants.

Modeled as an autologistic distribution.

We assume the endives are a single genotype.

Besag (1978) may have had data taken at 4 time points. This data was extracted from Friel and Pettitt. It is not clear what, if any, time point was used.

Friel does not give the dimensions. Besag is not available.

Source

J Besag (1978). Some Methods of Statistical Analysis for Spatial Data. Bulletin of the International Statistical Institute, 47, 77-92.

References

N Friel & A. N Pettitt (2004). Likelihood Estimation and Inference for the Autologistic Model. Journal of Computational and Graphical Statistics, 13:1, 232-246. https://doi.org/10.1198/1061860043029

Examples

## Not run: 
  
  library(agridat)
  data(besag.endive)
  dat <- besag.endive

  # Incidence map.  Figure 2 of Friel and Pettitt
  libs(desplot)
  grays <- colorRampPalette(c("#d9d9d9","#252525"))
  desplot(dat, disease~col*row,
          col.regions=grays(2),
          aspect = 0.5, # aspect unknown
          main="besag.endive - Disease incidence")
  
  
  # Besag (2000) "An Introduction to Markov Chain Monte Carlo" suggested
  # that the autologistic model is not a very good fit for this data.
  # We try it anyway.  No idea if this is correct or how to interpret...
  
  libs(ngspatial)
  A = adjacency.matrix(179,14)
  X = cbind(x=dat$col, y=dat$row)
  Z = as.numeric(dat$disease=="Y")
  m1 <- autologistic(Z ~ 0+X, A=A, control=list(confint="none"))
  
  summary(m1)
  ## Coefficients:
  ##      Estimate Lower Upper MCSE
  ## Xx  -0.007824    NA    NA   NA
  ## Xy  -0.144800    NA    NA   NA
  ## eta  0.806200    NA    NA   NA

  
  if(require("asreml", quietly=TRUE)) {
    libs(asreml,lucid)
    
    # Now try an AR1xAR1 model.
    dat2 <- transform(dat, xf=factor(col), yf=factor(row),
                      pres=as.numeric(disease=="Y"))
    
    m2 <- asreml(pres ~ 1, data=dat2,
                 resid = ~ar1(xf):ar1(yf))
    # The 0/1 response is arbitrary, but there is some suggestion
    # of auto-correlation in the x (.17) and y (.10) directions,
    # suggesting the pattern is more 'patchy' than just random noise,
    # but is it meaningful?
    
    lucid::vc(m2)
    ##       effect component std.error z.ratio bound 
    ##     xf:yf(R)   0.1301   0.003798    34       P   0
    ## xf:yf!xf!cor   0.1699   0.01942      8.7     U   0
    ## xf:yf!yf!cor   0.09842  0.02038      4.8     U   0
  }


## End(Not run)

Multi-environment trial of corn, incomplete-block design

Description

Multi-environment trial of corn, incomplete-block designlocation.

Format

A data frame with 1152 observations on the following 7 variables.

county

county

row

row

col

column

rep

rep

block

incomplete block

yield

yield

gen

genotype, 1-64

Details

Multi-environment trial of 64 corn hybrids in six counties in North Carolina. Each location had 3 replicates in in incomplete-block design with an 18x11 lattice of plots whose length-to-width ratio was about 2:1.

Note: In the original data, each county had 6 missing plots. This data has rows for each missing plot that uses the same county/block/rep to fill-out the row, sets the genotype to G01, and sets the yield to missing. These missing values were added to the data so that asreml could more easily do AR1xAR1 analysis using rectangular regions.

Each location/panel is:

Field length: 18 rows * 2 units = 36 units.

Field width: 11 plots * 1 unit = 11 units.

Retrieved from https://web.archive.org/web/19990505223413/www.stat.duke.edu/~higdon/trials/nc.dat

Used with permission of David Higdon.

Source

Julian Besag and D Higdon, 1999. Bayesian Analysis of Agricultural Field Experiments, Journal of the Royal Statistical Society: Series B, 61, 691–746. Table 1. https://doi.org/10.1111/1467-9868.00201

Examples

## Not run: 

  library(agridat)
  data(besag.met)
  dat <- besag.met

  libs(desplot)
  desplot(dat, yield ~ col*row|county,
          aspect=36/11, # true aspect
          out1=rep, out2=block,
          main="besag.met")


  # Average reps
  datm <- aggregate(yield ~ county + gen, data=dat, FUN=mean)
  
  # Sections below fit heteroskedastic variance models (variance for each variety)
  # asreml takes 1 second, lme 73 seconds, SAS PROC MIXED 30 minutes



  # lme
  # libs(nlme)
  # m1l <- lme(yield ~ -1 + gen, data=datm, random=~1|county,
  #            weights = varIdent(form=~ 1|gen))
  # m1l$sigma^2 * c(1, coef(m1l$modelStruct$varStruct, unc = FALSE))^2
  ##           G02    G03    G04    G05    G06    G07    G08
  ##  91.90 210.75  63.03 112.05  28.39 237.36  72.72  42.97
  ## ... etc ...
  
  if(require("asreml", quietly=TRUE)) {
   libs(asreml, lucid)

   # Average reps
   datm <- aggregate(yield ~ county + gen, data=dat, FUN=mean)
   #  asreml Using 'rcov' ALWAYS requires sorting the data
   datm <- datm[order(datm$gen),]
   
   m1 <- asreml(yield ~ gen, data=datm,
                random = ~ county,
                residual = ~ dsum( ~ units|gen))
   vc(m1)[1:7,]
   ##      effect component std.error z.ratio bound 
   ##    county   1324       836.1      1.6     P 0.2
   ## gen_G01!R     91.98     58.91     1.6     P 0.1
   ## gen_G02!R    210.6     133.6      1.6     P 0.1
   ## gen_G03!R     63.06     40.58     1.6     P 0.1
   ## gen_G04!R    112.1      71.59     1.6     P 0.1
   ## gen_G05!R     28.35     18.57     1.5     P 0.2
   ## gen_G06!R    237.4     150.8      1.6     P 0  
  
   # We get the same results from asreml & lme
   # plot(m1$vparameters[-1],
   #      m1l$sigma^2 * c(1, coef(m1l$modelStruct$varStruct, unc = FALSE))^2)
   
   # The following example shows how to construct a GxE biplot
   # from the FA2 model.
   
   
   dat <- besag.met
   dat <- transform(dat, xf=factor(col), yf=factor(row))
   dat <- dat[order(dat$county, dat$xf, dat$yf), ]
   
   # First, AR1xAR1
   m1 <- asreml(yield ~ county, data=dat,
                random = ~ gen:county,
                residual = ~ dsum( ~ ar1(xf):ar1(yf)|county))
   # Add FA1
   m2 <- update(m1, random=~gen:fa(county,1)) # rotate.FA=FALSE
   # FA2
   m3 <- update(m2, random=~gen:fa(county,2))
   asreml.options(extra=50)
   m3 <- update(m3, maxit=50)
   asreml.options(extra=0)
   
   # Use the loadings to make a biplot
   vars <- vc(m3)
   psi <- vars[grepl("!var$", vars$effect), "component"]
   la1 <- vars[grepl("!fa1$", vars$effect), "component"]
   la2 <- vars[grepl("!fa2$", vars$effect), "component"]
   mat <- as.matrix(data.frame(psi, la1, la2))
   # I tried using rotate.fa=FALSE, but it did not seem to
   # give orthogonal vectors.  Rotate by hand.
   rot <- svd(mat[,-1])$v # rotation matrix
   lam <- mat[,-1] 
   colnames(lam) <- c("load1", "load2")
   
   co3 <- coef(m3)$random # Scores are the GxE coefficients
   ix1 <- grepl("_Comp1$", rownames(co3))
   ix2 <- grepl("_Comp2$", rownames(co3))
   sco <- matrix(c(co3[ix1], co3[ix2]), ncol=2, byrow=FALSE)
   sco <- sco 
   dimnames(sco) <- list(levels(dat$gen) , c('load1','load2'))
   rownames(lam) <- levels(dat$county)
   sco[,1:2] <- -1 * sco[,1:2]
   lam[,1:2] <- -1 * lam[,1:2]
   biplot(sco, lam, cex=.5, main="FA2 coefficient biplot (asreml)")
   # G variance matrix
   gvar <- lam 
  
   # Now get predictions and make an ordinary biplot
   p3 <- predict(m3, data=dat, classify="county:gen")
   p3 <- p3$pvals
   libs("gge")  
   bi3 <- gge(p3, predicted.value ~ gen*county, scale=FALSE)
   if(interactive()) dev.new()
   # Very similar to the coefficient biplot
   biplot(bi3, stand=FALSE, main="SVD biplot of FA2 predictions")
  }
  

## End(Not run)

Four-way factorial agronomic experiment in triticale

Description

Four-way factorial agronomic experiment in triticale

Usage

data("besag.triticale")

Format

A data frame with 54 observations on the following 7 variables.

yield

yield, g/m^2

row

row

col

column

gen

genotype / variety, 3 levels

rate

seeding rate, kg/ha

nitro

nitrogen rate, kw/ha

regulator

growth regulator, 3 levels

Details

Experiment conducted as a factorial on the yields of triticale. Fully randomized. Plots were 1.5m x 5.5m, but the orientation is not clear.

Besag and Kempton show how accounting for neighbors changes non-significant genotype differences into significant differences.

Source

Julian Besag and Rob Kempton (1986). Statistical Analysis of Field Experiments Using Neighbouring Plots. Biometrics, 42, 231-251. Table 2. https://doi.org/10.2307/2531047

References

None.

Examples

## Not run: 

  library(agridat)
  data(besag.triticale)
  dat <- besag.triticale
  dat <- transform(dat, rate=factor(rate), nitro=factor(nitro))
  dat <- transform(dat, xf=factor(col), yf=factor(row))

  libs(desplot)
  desplot(dat, yield ~ col*row,
          # aspect unknown
          main="besag.triticale")

  # Besag & Kempton are not perfectly clear on the model, but
  # indicate that there was no evidence of any two-way interactions.
  # A reduced, main-effect model had genotype effects that were
  # "close to significant" at the five percent level.
  # The model below has p-value of gen at .04, so must be slightly
  # different than their model.
  m2 <- lm(yield ~ gen + rate + nitro + regulator + yf, data=dat)
  anova(m2)

  # Similar, but not exact, to Besag figure 5
  dat$res <- resid(m2)
  libs(lattice)
  xyplot(res ~ col|as.character(row), data=dat,
         as.table=TRUE, type="s", layout=c(1,3),
         main="besag.triticale")
  
  if(require("asreml", quietly=TRUE)) {
    libs(asreml)

    # Besag uses an adjustment based on neighboring plots.
    # This analysis fits the standard AR1xAR1 residual model
    
    dat <- dat[order(dat$xf, dat$yf), ]
    m3 <- asreml(yield ~ gen + rate + nitro + regulator +
                   gen:rate + gen:nitro + gen:regulator +
                   rate:nitro + rate:regulator +
                   nitro:regulator + yf, data=dat,
                 resid = ~ ar1(xf):ar1(yf))
    wald(m3) # Strongly significant gen, rate, regulator
    ##                 Df Sum of Sq Wald statistic Pr(Chisq)    
    ## (Intercept)      1   1288255        103.971 < 2.2e-16 ***
    ## gen              2    903262         72.899 < 2.2e-16 ***
    ## rate             1    104774          8.456  0.003638 ** 
    ## nitro            1       282          0.023  0.880139    
    ## regulator        2    231403         18.676 8.802e-05 ***
    ## yf               2      3788          0.306  0.858263    
    ## gen:rate         2      1364          0.110  0.946461    
    ## gen:nitro        2     30822          2.488  0.288289    
    ## gen:regulator    4     37269          3.008  0.556507    
    ## rate:nitro       1      1488          0.120  0.728954    
    ## rate:regulator   2     49296          3.979  0.136795    
    ## nitro:regulator  2     41019          3.311  0.191042    
    ## residual (MS)          12391                             
  }
  

## End(Not run)

Multi-environment trial of wheat, conventional and semi-dwarf varieties

Description

Multi-environment trial of wheat, conventional and semi-dwarf varieties, 7 locs with low/high fertilizer levels.

Format

A data frame with 168 observations on the following 5 variables.

gen

genotype

loc

location

nitro

nitrogen fertilizer, low/high

yield

yield (g/m^2)

type

type factor, conventional/semi-dwarf

Details

Conducted in U.K. in 1975. Each loc had three reps, two nitrogen treatments.

Locations were Begbroke, Boxworth, Crafts Hill, Earith, Edinburgh, Fowlmere, Trumpington.

At the two highest-yielding locations, Earith and Edinburgh, yield was _lower_ for the high-nitrogen treatment. Blackman et al. say "it seems probable that effects on development and structure of the crop were responsible for the reductions in yield at high nitrogen".

Source

Blackman, JA and Bingham, J. and Davidson, JL (1978). Response of semi-dwarf and conventional winter wheat varieties to the application of nitrogen fertilizer. The Journal of Agricultural Science, 90, 543–550. https://doi.org/10.1017/S0021859600056070

References

Gower, J. and Lubbe, S.G. and Gardner, S. and Le Roux, N. (2011). Understanding Biplots, Wiley.

Examples

## Not run: 

library(agridat)
data(blackman.wheat)
dat <- blackman.wheat

libs(lattice)

# Semi-dwarf generally higher yielding than conventional
# bwplot(yield~type|loc,dat, main="blackman.wheat")

# Peculiar interaction--Ear/Edn locs have reverse nitro response
dotplot(gen~yield|loc, dat, group=nitro, auto.key=TRUE,
        main="blackman.wheat: yield for low/high nitrogen")

# Height data from table 6 of Blackman.  Height at Trumpington loc.
# Shorter varieties have higher yields, greater response to nitro.
heights <- data.frame(gen=c("Cap", "Dur", "Fun", "Hob", "Hun", "Kin",
                            "Ran", "Spo", "T64", "T68","T95", "Tem"),
                      ht=c(101,76,76,80,98,88,98,81,86,73,78,93))
dat$height <- heights$ht[match(dat$gen, heights$gen)]
xyplot(yield~height|loc,dat,group=nitro,type=c('p','r'),
       main="blackman.wheat",
       subset=loc=="Tru", auto.key=TRUE)


libs(reshape2)
# AMMI-style biplot Fig 6.4 of Gower 2011
dat$env <- factor(paste(dat$loc,dat$nitro,sep="-"))
datm <- acast(dat, gen~env, value.var='yield')
datm <- sweep(datm, 1, rowMeans(datm))
datm <- sweep(datm, 2, colMeans(datm))
biplot(prcomp(datm), main="blackman.wheat AMMI-style biplot")



## End(Not run)

Corn borer infestation under four treatments

Description

Corn borer infestation under four treatments

Format

A data frame with 48 observations on the following 3 variables.

borers

number of borers per hill

treat

treatment factor

freq

frequency of the borer count

Details

Four treatments to control corn borers. Treatment 1 is the control.

In 15 blocks, for each treatment, 8 hills of plants were examined, and the number of corn borers present was recorded. The data here are aggregated across blocks.

Bliss mentions that the level of infestation varied significantly between the blocks.

Source

C. Bliss and R. A. Fisher. (1953). Fitting the Negative Binomial Distribution to Biological Data. Biometrics, 9, 176–200. Table 3. https://doi.org/10.2307/3001850

Geoffrey Beall. 1940. The Fit and Significance of Contagious Distributions when Applied to Observations on Larval Insects. Ecology, 21, 460-474. Page 463. https://doi.org/10.2307/1930285

Examples

## Not run: 

library(agridat)
data(bliss.borers)
dat <- bliss.borers

# Add 0 frequencies
dat0 <- expand.grid(borers=0:26, treat=c('T1','T2','T3','T4'))
dat0 <- merge(dat0,dat, all=TRUE)
dat0$freq[is.na(dat0$freq)] <- 0

# Expand to individual (non-aggregated) counts for each hill
dd <- data.frame(borers = rep(dat0$borers, times=dat0$freq),
                 treat = rep(dat0$treat, times=dat0$freq))

libs(lattice)
histogram(~borers|treat, dd, type='count', breaks=0:27-.5,
          layout=c(1,4), main="bliss.borers", xlab="Borers per hill")


libs(MASS)
  m1 <- glm.nb(borers~0+treat, data=dd)
  # Bliss, table 3, presents treatment means, which are matched by:
  exp(coef(m1)) # 4.033333 3.166667 1.483333 1.508333
  # Bliss gives treatment values k = c(1.532,1.764,1.333,1.190).
  # The mean of these is 1.45, similar to this across-treatment estimate
  m1$theta # 1.47


# Plot observed and expected distributions for treatment 2
libs(latticeExtra)
  xx <- 0:26
  yy <- dnbinom(0:26, mu=3.17, size=1.47)*120 # estimates are from glm.nb
  histogram(~borers, dd, type='count', subset=treat=='T2',
            main="bliss.borers - trt T2 observed and expected",
            breaks=0:27-.5) +
              xyplot(yy~xx, col='navy', type='b')


# "Poissonness"-type plot
libs(vcd)
  dat2 <- droplevels(subset(dat, treat=='T2'))
  vcd::distplot(dat2$borers, type = "nbinomial",
           main="bliss.borers neg binomialness plot")
  # Better way is a rootogram
  g1 <- vcd::goodfit(dat2$borers, "nbinomial")
  plot(g1, main="bliss.borers - Treatment 2")


## End(Not run)

Diallel cross of winter beans

Description

Diallel cross of winter beans

Format

A data frame with 36 observations on the following 3 variables.

female

female parent

male

male parent

yield

yield, grams/plot

stems

stems per plot

nodes

podded nodes per stem

pods

pods per podded node

seeds

seeds per pod

weight

weight (g) per 100 seeds

height

height (cm) in April

width

width (cm) in April

flower

mean flowering date in May

Details

Yield in grams/plot for full diallel cross between 6 inbred lines of winter beans. Values are means over two years.

Source

D. A. Bond (1966). Yield and components of yield in diallel crosses between inbred lines of winter beans (Viciafaba). The Journal of Agricultural Science, 67, 325–336. https://doi.org/10.1017/S0021859600017329

References

Peter John, Statistical Design and Analysis of Experiments, p. 85.

Examples

## Not run: 
  
  library(agridat)
  data(bond.diallel)
  dat <- bond.diallel
  
  # Because these data are means, we will not be able to reproduce
  # the anova table in Bond. More useful as a multivariate example.

  libs(corrgram)
  corrgram(dat[ , 3:11], main="bond.diallel",
           lower=panel.pts)

  # Multivariate example from sommer package
  corrgram(dat[,c("stems","pods","seeds")],
           lower=panel.pts, upper=panel.conf, main="bond.diallel")
  
  libs(sommer)           
  m1 <- mmer(cbind(stems,pods,seeds) ~ 1,
             random= ~ vs(female)+vs(male),
             rcov= ~ vs(units),
             dat)

  #### genetic variance covariance
  cov2cor(m1$sigma$`u:female`)
  cov2cor(m1$sigma$`u:male`)
  cov2cor(m1$sigma$`u:units`)


## End(Not run)

Uniformity trials of barley, wheat, lentils

Description

Uniformity trials of barley, wheat, lentils in India 1930-1932.

Usage

data("bose.multi.uniformity")

Format

A data frame with 1170 observations on the following 5 variables.

year

year

crop

crop

row

row ordinate

col

column ordinate

yield

yield per plot in grams

Details

A field about 1/4 acre was sown in three consecutive years (beginning in 1929-1930) with barley, wheat, and lentil.

At harvest, borders 3 feet on east and west and 6 feet on north and south were removed. The field was divided into plots four feet square, which were harvested separately, measured in grams.

Fertility contours of the field were somewhat similar across years, with correlation values across years 0.45, 0.48, 0.21.

Field width: 15 plots * 4 feet = 60 feet.

Field length: 26 plots * 4 feet = 104 feet.

Conclusions:

"An experimental field which may be sensibly uniform for one crop or for one season may not be so for another crop or in a different season" p. 592.

Source

Bose, R. D. (1935). Some soil heterogeneity trials at Pusa and the size and shape of experimental plots. Ind. J. Agric. Sci., 5, 579-608. Table 1 (p. 585), Table 4 (p. 589), Table 5 (p. 590). https://archive.org/details/in.ernet.dli.2015.271739

References

Shaw (1935). Handbook of Statistics for Use in Plant-Breeding and Agricultural Problems, p. 149-170. https://krishikosh.egranth.ac.in/handle/1/21153

Examples

## Not run: 

  library(agridat)
  data(bose.multi.uniformity)
  dat <- bose.multi.uniformity

  # match sum at bottom of Bose tables 1, 4, 5
  # library(dplyr)
  # dat 

  libs(desplot, dplyr)
  # Calculate percent of mean yield for each year
  dat <- group_by(dat, year)
  dat <- mutate(dat, pctyld = (yield-mean(yield))/mean(yield))

  dat <- ungroup(dat)
  dat <- mutate(dat, year=as.character(year))
  # Bose smoothed the data by averaging 2x3 plots together before drawing
  # contour maps.  Heatmaps of raw data have similar structure to Bose Fig 1.
  desplot(dat, pctyld ~ col*row|year,
          tick=TRUE, flip=TRUE, aspect=(26)/(15),
          main="bose.multi.* - Percent of mean yield")
  
  # contourplot() results need to be mentally flipped upside down
  # contourplot(pctyld ~ col*row|year, dat,
  #   region=TRUE, as.table=TRUE, aspect=26/15)


## End(Not run)

Weight of cork samples on four sides of trees

Description

The cork data gives the weights of cork borings of the trunk for 28 trees on the north (N), east (E), south (S) and west (W) directions.

Format

Data frame with 28 observations on the following 5 variables.

tree

tree number

dir

direction N,E,S,W

y

weight of cork deposit (centigrams), north direction

Source

C.R. Rao (1948). Tests of significance in multivariate analysis. Biometrika, 35, 58-79. https://doi.org/10.2307/2332629

References

K.V. Mardia, J.T. Kent and J.M. Bibby (1979) Multivariate Analysis, Academic Press.

Russell D Wolfinger, (1996). Heterogeneous Variance: Covariance Structures for Repeated Measures. Journal of Agricultural, Biological, and Environmental Statistics, 1, 205-230.

Examples

## Not run: 

  library(agridat)
  data(box.cork)
  dat <- box.cork

  libs(reshape2, lattice)
  dat2 <- acast(dat, tree ~ dir, value.var='y')
  splom(dat2, pscales=3,
        prepanel.limits = function(x) c(25,100),
        main="box.cork", xlab="Cork yield on side of tree",
        panel=function(x,y,...){
          panel.splom(x,y,...)
          panel.abline(0,1,col="gray80")
        })


  ## Radial star plot, each tree is one line
  libs(plotrix)
  libs(reshape2)
  dat2 <- acast(dat, tree ~ dir, value.var='y')
  radial.plot(dat2, start=pi/2, rp.type='p', clockwise=TRUE,
              radial.lim=c(0,100), main="box.cork",
              lwd=2, labels=c('North','East','South','West'),
              line.col=rep(c("royalblue","red","#009900","dark orange",
                             "#999999","#a6761d","deep pink"),
                           length=nrow(dat2)))

  if(require("asreml", quietly=TRUE)) {  
    libs(asreml, lucid)
    
    # Unstructured covariance
    dat$dir <- factor(dat$dir)
    dat$tree <- factor(dat$tree)  
    dat <- dat[order(dat$tree, dat$dir), ]
    
    # Unstructured covariance matrix
    m1 <- asreml(y~dir, data=dat, residual = ~ tree:us(dir))
    
    lucid::vc(m1)
    
    # Note: 'rcor' is a personal function to extract the correlations
    # into a matrix format
    # round(kw::rcor(m1)$dir, 2)
    #        E      N      S      W
    # E 219.93 223.75 229.06 171.37
    # N 223.75 290.41 288.44 226.27
    # S 229.06 288.44 350.00 259.54
    # W 171.37 226.27 259.54 226.00
    
    # Note: Wolfinger used a common diagonal variance
    
    # Factor Analytic with different specific variances
    # fixme: does not work with asreml4
    # m2 <- update(m1, residual = ~tree:facv(dir,1))
    # round(kw::rcor(m2)$dir, 2)
    #       E       N      S      W
    # E 219.94 209.46 232.85 182.27
    # N 209.46 290.41 291.82 228.43
    # S 232.85 291.82 349.99 253.94
    # W 182.27 228.43 253.94 225.99
  }
  

## End(Not run)

Uniformity trial of 4 crops on the same land

Description

Uniformity trial of 4 crops on the same land in Trinidad.

Usage

data("bradley.multi.uniformity")

Format

A data frame with 440 observations on the following 5 variables.

row

row

col

column

yield

yield, pounds per plot

season

season

crop

crop

Details

Experiments conducted in Trinidad.

Plots were marked in May 1939 in Fields 1, 2, and 3. Prior to 1939 it was difficult to obtain significant results on this land.

Plots were 1/40 acre each, 33 feet square. Discard between blocks (the rows) was 7 feet and between plots (the columns) was 4 feet. For roadways, a gap of 14 feet is between blocks 10 and 11 and a gap of 10 feet between plots E/F (which we call columns 5/6).

Data was collected for 4 crops. Two other crops had poor germination and were omitted.

Field width: 10 plots * 33 feet + 8 gaps * 4 feet + 1 gap * 10 = 372 feet

Field length: 11 blocks (plots) * 33 feet + 9 gaps * 7 feet + 1 gap * 14 feet = 440 feet

Crop 1. Woolly Pyrol. Crop cut at flowering and weighed in pounds. Note, woolly pyrol appears to be a bean also called black gram, phaseolus mungo.

Crop 2. Woolly Pyrol. Crop cut at flowering and weighed in pounds.

Crop 3. Maize. Net weight of cobs in pounds. Source document also has number of cobs.

Crop 4. Yams. Weights in pounds. Source document has weight to 1/4 pound, which has here been rounded to the nearest pound. (Half pounds were rounded to nearest even pound.) Source document also has number of yams.

Notes by Bradley.

The edges of the field tended to be slightly higher yielding. Thought to be due to the heavier cultivation which the edges recieve (p. 18).

The plot in row 9, col 7 (9G in Bradley) is higher yielding than its neighbors, thought to be the site of a saman tree dug up and burned when the field was plotted. Bits of charcoal were still in the soil.

Bradley also examined soil samples on selected plots and looked at nutrients, moisture, texture, etc. The selected plots were 4 high-yielding plots and 4 low-yielding plots. Little difference was observed. Unexpectedly, yams gave higher yield on plots with more compaction.

Source

P. L. Bradley (1941). A study of the variation in productivity over a number of fixed plots in field 2. Dissertation: The University of the West Indies. Appendix 1a, 1b, 1c, 1d. https://uwispace.sta.uwi.edu/items/e874561d-52e5-4e39-8416-ff8c1756049c https://hdl.handle.net/2139/41259

The data are repeated in: C. E. Wilson. Study of the plots laid out on field II with a view to obtaining plot-fertility data for use in future experiments on these plots, season 1940-41. Dissertation: The University of the West Indies. Page 36-39. https://uwispace.sta.uwi.edu/dspace/handle/2139/43658

References

None

Examples

## Not run: 

library(agridat)
data(bradley.multi.uniformity)
dat <- bradley.multi.uniformity

# figures similar to Bradley, pages 11-15
libs(desplot)
desplot(dat, yield ~ col*row, subset=season==1,
        flip=TRUE, aspect=433/366, # true aspect (omits roadways)
        main="bradley.multi.uniformity - season 1, woolly pyrol")

desplot(dat, yield ~ col*row, subset=season==2,
        flip=TRUE, aspect=433/366, # true aspect (omits roadways)
        main="bradley.multi.uniformity - season 2, woolly pyrol")

desplot(dat, yield ~ col*row, subset=season==3,
        flip=TRUE, aspect=433/366, # true aspect (omits roadways)
        main="bradley.multi.uniformity - season 3, maize")


desplot(dat, yield ~ col*row, subset=season==4,
        flip=TRUE, aspect=433/366, # true aspect (omits roadways)
        main="bradley.multi.uniformity - season 4, yams")

dat1 <- subset(bradley.multi.uniformity, season==1)
dat2 <- subset(bradley.multi.uniformity, season==2)
dat3 <- subset(bradley.multi.uniformity, season==3)
dat4 <- subset(bradley.multi.uniformity, season==4)
  # to combine plots across seasons, each yield value was converted to percent
  # of maximum yield in that season. Same as Bradley, page 17.
  dat1$percent <- dat1$yield / max(dat1$yield) * 100
  dat2$percent <- dat2$yield / max(dat2$yield) * 100
  dat3$percent <- dat3$yield / max(dat3$yield) * 100
  dat4$percent <- dat4$yield / max(dat4$yield) * 100
  # make sure data is in same order, then combine
  dat1 <- dat1[order(dat1$col, dat1$row),]
  dat2 <- dat2[order(dat2$col, dat2$row),]
  dat3 <- dat3[order(dat3$col, dat3$row),]
  dat4 <- dat4[order(dat4$col, dat4$row),]
  dat14 <- dat1[,c('row','col')]
  dat14$fertility <- dat1$percent + dat2$percent + dat3$percent + dat4$percent

  libs(desplot)
  desplot(dat14, fertility ~ col*row,
          tick=TRUE, flip=TRUE, aspect=433/366, # true aspect (omits roadways)
          main="bradley.multi.uniformity - fertility")


## End(Not run)

Multi-environment trial of rape in Manitoba

Description

Rape seed yields for 5 genotypes, 3 years, 9 locations.

Format

A data frame with 135 observations on the following 4 variables.

gen

genotype

year

year, numeric

loc

location, 9 levels

yield

yield, kg/ha

Details

The yields are the mean of 4 reps.

Note, in table 2 of Brandle, the value of Triton in 1985 at Bagot is shown as 2355, but should be 2555 to match the means reported in the paper.

Used with permission of P. McVetty.

Source

Brandle, JE and McVetty, PBE. (1988). Genotype x environment interaction and stability analysis of seed yield of oilseed rape grown in Manitoba. Canadian Journal of Plant Science, 68, 381–388.

Examples

## Not run: 

library(agridat)
data(brandle.rape)
dat <- brandle.rape

libs(lattice)
dotplot(gen~yield|loc, dat, group=year, auto.key=list(columns=3),
        main="brandle.rape, yields per location", ylab="Genotype")

# Matches table 4 of Brandle
# round(tapply(dat$yield, dat$gen, mean),0)

# Brandle reports variance components:
# sigma^2_gl: 9369  gy: 14027 g: 72632 resid: 150000
# Brandle analyzed rep-level data, so the residual variance is different.
# The other components are matched by the following analysis.

libs(lme4)
libs(lucid)
dat$year <- factor(dat$year)
m1 <- lmer(yield ~ year + loc + year:loc + (1|gen) +
             (1|gen:loc) + (1|gen:year), data=dat)
vc(m1)
##      grp        var1 var2  vcov  sdcor
##  gen:loc (Intercept) <NA>  9363  96.76
## gen:year (Intercept) <NA> 14030 118.4
##      gen (Intercept) <NA> 72630 269.5
## Residual        <NA> <NA> 75010 273.9
  

## End(Not run)

Switchback experiment on dairy cattle, milk yield for two treatments

Description

Switchback experiment on dairy cattle, milk yield for two treatments

Usage

data("brandt.switchback")

Format

A data frame with 30 observations on the following 5 variables.

group

group: A,B

cow

cow, 10 levels

trt

treatment, 2 levels

period

period, 3 levels

yield

milk yield, pounds

Details

In this experiment, 10 cows were selected from the Iowa State College Holstein-Friesian herd and divided into two equal groups. Care was taken to have the groups as nearly equal as possible with regard to milk production, stage of gestation, body weight, condition and age. These cows were each given 10 pounds of timothy hay and 30 pounds of corn silage daily but were fed different grain mixtures. Treatment T1, then, consisted of feeding a grain mixture of 1 part of corn and cob meal to 1 part of ground oats, while treatment T2 consisted of feeding a grain mixture of 4 parts corn and cob meal, 4 parts of ground oats and 3 parts of gluten feed. The three treatment periods covered 105 days – three periods of 35 days each. The yields for the first 7 days of each period were not considered because of the possible effect of the transition from one treatment to the other. The data, together with sums and differences which aid in the calculations incidental to testing, are given in table 2.

It seems safe to conclude that the inclusion of gluten feed in the grain mixture fed in a timothy hay ration to Holstein-Friesian cows increased the production of milk. The average increase was 21.7 pounds per cow for a 28-day period.

Source

A.E. Brandt (1938). Tests of Significance in Reversal or Switchback Trials Iowa State College, Agricultural Research Bulletins. Bulletin 234. Book 22. https://lib.dr.iastate.edu/ag_researchbulletins/22/

Examples

## Not run: 
  
library(agridat)

data(brandt.switchback)
dat <- brandt.switchback

# In each period, treatment 2 is slightly higher
# bwplot(yield~trt|period,dat, layout=c(3,1), main="brandt.switchback",
#     xlab="Treatment", ylab="Milk yield")

# Yield at period 2 (trt T2) is above the trend in group A,
# below the trend (trt T1) in group B.
# Equivalently, treatment T2 is above the trend line
libs(lattice)
xyplot(yield~period|group, data=dat, group=cow, type=c('l','r'),
    auto.key=list(columns=5), main="brandt.switchback",
    xlab="Period.  Group A: T1,T2,T1.  Group B: T2,T1,T2",
    ylab="Milk yield (observed and trend) per cow")

# Similar to Brandt Table 10
m1 <- aov(yield~period+group+cow:group+period:group, data=dat)
anova(m1)


## End(Not run)

Multi-environment trial of cucumbers in a latin square design

Description

Cucumber yields in latin square design at two locs.

Format

A data frame with 32 observations on the following 5 variables.

loc

location

gen

genotype/cultivar

row

row

col

column

yield

weight of marketable fruit per plot

Details

Conducted at Clemson University in 1985. four cucumber cultivars were grown in a latin square design at Clemson, SC, and Tifton, GA.

Separate variances are modeled each location.

Plot dimensions are not given.

Bridges (1989) used this data to illustrate fitting a heterogeneous mixed model.

Used with permission of William Bridges.

Source

William Bridges (1989). Analysis of a plant breeding experiment with heterogeneous variances using mixed model equations. Applications of mixed models in agriculture and related disciplines, S. Coop. Ser. Bull, 45–51.

Examples

## Not run: 

  library(agridat)
  data(bridges.cucumber)
  dat <- bridges.cucumber
  dat <- transform(dat, rowf=factor(row), colf=factor(col))

  libs(desplot)
  desplot(dat, yield~col*row|loc,
          # aspect unknown
          text=gen, cex=1,
          main="bridges.cucumber")

  # Graphical inference test for heterogenous variances
  libs(nullabor)
  # Create a lineup of datasets
  fun <- null_permute("loc")
  dat20 <- lineup(fun, dat, n=20, pos=9)

  # Now plot
  libs(lattice)
  bwplot(yield ~ loc|factor(.sample), dat20,
         main="bridges.cucumber - graphical inference")

  if(require("asreml", quietly=TRUE)) {
    libs(asreml,lucid)
    
    ## Random row/col/resid. Same as Bridges 1989, p. 147
    m1 <- asreml(yield ~ 1 + gen + loc + loc:gen,
                 random = ~ rowf:loc + colf:loc, data=dat)
  
    lucid::vc(m1)
    ##   effect component std.error z.ratio bound 
    ## rowf:loc     31.62     23.02     1.4     P   0
    ## colf:loc     18.08     15.32     1.2     P   0
    ## units(R)     31.48     12.85     2.4     P   0
    
    ## Random row/col/resid at each loc. Matches p. 147
    m2 <- asreml(yield ~ 1 + gen + loc + loc:gen,
                 random = ~ at(loc):rowf + at(loc):colf, data=dat,
                 resid = ~ dsum( ~ units|loc))
    lucid::vc(m2)
    ##                effect component std.error z.ratio bound 
    ## at(loc, Clemson):rowf     32.32    36.58     0.88     P   0
    ##  at(loc, Tifton):rowf     30.92    28.63     1.1      P   0
    ## at(loc, Clemson):colf     22.55    28.78     0.78     P   0
    ##  at(loc, Tifton):colf     13.62    14.59     0.93     P   0
    ##        loc_Clemson(R)     46.85    27.05     1.7      P   0
    ##         loc_Tifton(R)     16.11     9.299    1.7      P   0
    
    predict(m2, data=dat, classify='loc:gen')$pvals
    ##       loc      gen predicted.value std.error    status
    ## 1 Clemson   Dasher            45.6      5.04 Estimable
    ## 2 Clemson Guardian            31.6      5.04 Estimable
    ## 3 Clemson Poinsett            21.4      5.04 Estimable
    ## 4 Clemson   Sprint            26        5.04 Estimable
    ## 5  Tifton   Dasher            50.5      3.89 Estimable
    ## 6  Tifton Guardian            38.7      3.89 Estimable
    ## 7  Tifton Poinsett            33        3.89 Estimable
    ## 8  Tifton   Sprint            39.2      3.89 Estimable
    
    # Is a heterogeneous model justified? Maybe not.
    # m1$loglik
    ## -67.35585
    # m2$loglik
    ## -66.35621
  }
  

## End(Not run)

Long term wheat yields on Broadbalk fields at Rothamsted.

Description

Long term wheat yields on Broadbalk fields at Rothamsted.

Format

A data frame with 1258 observations on the following 4 variables.

year

year

plot

plot

grain

grain yield, tonnes

straw

straw yield, tonnes

Details

Note: This data is only 1852-1925. You can find recent data for these experiments at the Electronic Rothamsted Archive: https://www.era.rothamsted.ac.uk/

Rothamsted Experiment station conducted wheat experiments on the Broadbalk Fields beginning in 1844 with data for yields of grain and straw collected from 1852 to 1925. Ronald Fisher was hired to analyze data from the agricultural trials. Organic manures and inorganic fertilizer treatments were applied in various combinations to the plots.

N1 is 48kg, N1.5 is 72kg, N2 is 96kg, N4 is 192kg nitrogen.

Electronic version of the data was retrieved from http://lib.stat.cmu.edu/datasets/Andrews/

Source

D.F. Andrews and A.M. Herzberg. 1985. Data: A Collection of Problems from Many Fields for the Student and Research Worker. Springer.

References

Broadbalk Winter Wheat Experiment. https://www.era.rothamsted.ac.uk/index.php?area=home&page=index&dataset=4

Examples

## Not run: 
  
library(agridat)
data(broadbalk.wheat)
dat <- broadbalk.wheat

libs(lattice)
## xyplot(grain~straw|plot, dat, type=c('p','smooth'), as.table=TRUE,
##        main="broadbalk.wheat")
xyplot(grain~year|plot, dat, type=c('p','smooth'), as.table=TRUE,
       main="broadbalk.wheat") # yields are decreasing

# See the treatment descriptions to understand the patterns
redblue <- colorRampPalette(c("firebrick", "lightgray", "#375997"))
levelplot(grain~year*plot, dat, main="broadbalk.wheat: Grain", col.regions=redblue)


## End(Not run)

Uniformity trial of corn at 3 locations in Iowa.

Description

Uniformity trial of corn at 3 locations in Iowa.

Usage

data("bryan.corn.uniformity")

Format

A data frame with 1728 observations on the following 4 variables.

expt

experiment (variety/orientation)

row

row

col

column

yield

yield, pounds per plot

Details

Three varieties of corn were planted. Each experiment was 48 rows, each row 48 hills long, .65 acres. A "hill" is a single hole with possibly multiple seeds. Spacing between the hills would be sqrt(43560 sq ft * .64) / 48 = 3.5 feet.

In the experiment code, K=Krug, I=Iodent, M=McCulloch (varieties of corn), 23=1923, 25=1925, E=East/West, N=North/South.

Each experiment was aggregated into experimental units by combining 8 hills, both in East/West direction and also in North/South direction. Thus, each field is represented twice in the data, once with "E" in the field name and once with "N".

Source

Arthur Bryan (1933). Factors Affecting Experimental Error in Field Plot Tests With Corn. Agricultural Experiment Station, Iowa State College. Tables 22-27. https://hdl.handle.net/2027/uiug.30112019568168

References

None

Examples

## Not run: 
  library(agridat)
  data(bryan.corn.uniformity)
  dat <- bryan.corn.uniformity
  
  libs(desplot)
  desplot(dat, yield ~ col*row|expt,
          main="bryan.corn.uniformity",
          aspect=(48*3.5/(6*8*3.5)), # true aspect
          flip=TRUE, tick=TRUE)

  # CVs in Table 5, column 8 hills
  # libs(dplyr)
  # dat 
  #   summarize(cv=sd(yield)/mean(yield)*100)
  ##   expt    cv
  ## 1 K23E  10.9
  ## 2 K23N  10.9
  ## 3 I25E  16.3
  ## 4 I25N  17.0
  ## 5 M25E  16.2
  ## 6 M25N  17.2


## End(Not run)

Multi-environment trial of wheat in Sweden in 2016.

Description

Multi-environment trial of wheat in Sweden in 2016.

Usage

data("buntaran.wheat")

Format

A data frame with 1069 observations on the following 7 variables.

zone

Geographic zone: south, middle, north

loc

Location

rep

Block replicate (up to 4)

alpha

Incomplete-block in the alpha design

gen

Genotype (cultivar)

yield

Dry matter yield, kg/ha

Details

Dry matter yield from wheat trials in Sweden in 2016. The experiments in each location were multi-rep with incomplete blocks in an alpha design.

Electronic data are from the online supplement of Buntaran (2020) and also from the "init" package at https://github.com/Flavjack/inti.

Source

Buntaran, Harimurti et al. (2020). Cross-validation of stagewise mixed-model analysis of Swedish variety trials with winter wheat and spring barley. Crop Science, 60, 2221-2240. http://doi.org/10.1002/csc2.20177

References

None.

Examples

## Not run: 
data(buntaran.wheat)
library(agridat)
dat <- buntaran.wheat
library(lattice)
bwplot(yield~loc|zone, dat, layout=c(1,3),
       scales=list(x=list(rot=90)),
       main="buntaran.wheat")

## End(Not run)

Incomplete block alpha design

Description

Incomplete block alpha design

Usage

data("burgueno.alpha")

Format

A data frame with 48 observations on the following 6 variables.

rep

rep, 3 levels

block

block, 12 levels

row

row

col

column

gen

genotype, 16 levels

yield

yield

Details

A field experiment with 3 reps, 4 blocks per rep, laid out as an alpha design.

The plot size is not given.

Electronic version of the data obtained from CropStat software.

Used with permission of Juan Burgueno.

Source

J Burgueno, A Cadena, J Crossa, M Banziger, A Gilmour, B Cullis. 2000. User's guide for spatial analysis of field variety trials using ASREML. CIMMYT. https://books.google.com/books?id=PR_tYCFyLCYC&pg=PA1

Examples

## Not run: 

  library(agridat)
  data(burgueno.alpha)
  dat <- burgueno.alpha

  libs(desplot)
  desplot(dat, yield~col*row,
          out1=rep, out2=block, # aspect unknown
          text=gen, cex=1,shorten="none",
          main='burgueno.alpha')


  libs(lme4,lucid)
  # Inc block model
  m0 <- lmer(yield ~ gen + (1|rep/block), data=dat)
  vc(m0) # Matches Burgueno p. 26
  ##        grp        var1 var2   vcov sdcor
  ##  block:rep (Intercept) <NA>  86900 294.8
  ##        rep (Intercept) <NA> 200900 448.2
  ##   Residual        <NA> <NA> 133200 365  


  if(require("asreml", quietly=TRUE)) {
    libs(asreml)
    
    dat <- transform(dat, xf=factor(col), yf=factor(row))
    dat <- dat[order(dat$xf, dat$yf),]                 
    
    # Sequence of models on page 36 of Burgueno
    
    m1 <- asreml(yield ~  gen, data=dat)
    m1$loglik # -232.13
    
    m2 <- asreml(yield ~  gen, data=dat,
                 random = ~ rep)
    m2$loglik # -223.48
    
    # Inc Block model
    m3 <- asreml(yield ~  gen, data=dat,
                 random = ~ rep/block)
    m3$loglik # -221.42
    m3$coef$fixed # Matches solution on p. 27
    
    # AR1xAR1 model
    m4 <- asreml(yield ~ 1 + gen, data=dat,
                 resid = ~ar1(xf):ar1(yf))
    m4$loglik # -221.47
    plot(varioGram(m4), main="burgueno.alpha") # Figure 1
    
    m5 <- asreml(yield ~ 1 + gen, data=dat,
                 random= ~ yf, resid = ~ar1(xf):ar1(yf))
    m5$loglik # -220.07
    
    m6 <- asreml(yield ~ 1 + gen + pol(yf,-2), data=dat,
                 resid = ~ar1(xf):ar1(yf))
    m6$loglik # -204.64
    
    m7 <- asreml(yield ~ 1 + gen + lin(yf), data=dat,
                 random= ~ spl(yf), resid = ~ar1(xf):ar1(yf))
    m7$loglik # -212.51
    
    m8 <- asreml(yield ~ 1 + gen + lin(yf), data=dat,
                 random= ~ spl(yf))
    m8$loglik # -213.91
    
    # Polynomial model with predictions
    m9 <- asreml(yield ~ 1 + gen + pol(yf,-2) + pol(xf,-2), data=dat,
                 random= ~ spl(yf),
                 resid = ~ar1(xf):ar1(yf))
    m9 <- update(m9)
    m9$loglik # -191.44 vs -189.61
  
    m10 <- asreml(yield ~ 1 + gen + lin(yf)+lin(xf), data=dat,
                  resid = ~ar1(xf):ar1(yf))
    m10$loglik # -211.56
    
    m11 <- asreml(yield ~ 1 + gen + lin(yf)+lin(xf), data=dat,
                  random= ~ spl(yf),
                  resid = ~ar1(xf):ar1(yf))
    m11$loglik # -208.90
    
    m12 <- asreml(yield ~ 1 + gen + lin(yf)+lin(xf), data=dat,
                  random= ~ spl(yf)+spl(xf),
                  resid = ~ar1(xf):ar1(yf))
    m12$loglik # -206.82
    
    m13 <- asreml(yield ~ 1 + gen + lin(yf)+lin(xf), data=dat,
                  random= ~ spl(yf)+spl(xf))
    m13$loglik # -207.52
  }
  

## End(Not run)

Row-column design

Description

Row-column design

Usage

data("burgueno.rowcol")

Format

A data frame with 128 observations on the following 5 variables.

rep

rep, 2 levels

row

row

col

column

gen

genotype, 64 levels

yield

yield, tons/ha

Details

A field experiment with two contiguous replicates in 8 rows, 16 columns.

The plot size is not given.

Electronic version of the data obtained from CropStat software.

Used with permission of Juan Burgueno.

Source

J Burgueno, A Cadena, J Crossa, M Banziger, A Gilmour, B Cullis (2000). User's guide for spatial analysis of field variety trials using ASREML. CIMMYT.

Examples

## Not run: 

  library(agridat)
  data(burgueno.rowcol)
  dat <- burgueno.rowcol

  # Two contiguous reps in 8 rows, 16 columns
  libs(desplot)
  desplot(dat, yield ~ col*row,
          out1=rep, # aspect unknown
          text=gen, shorten="none", cex=.75,
          main="burgueno.rowcol")

  libs(lme4,lucid)
  
  # Random rep, row and col within rep
  # m1 <- lmer(yield ~ gen + (1|rep) + (1|rep:row) + (1|rep:col), data=dat)
  # vc(m1) # Match components of Burgueno p. 40
  ##      grp        var1 var2   vcov  sdcor
  ##  rep:col (Intercept) <NA> 0.2189 0.4679
  ##  rep:row (Intercept) <NA> 0.1646 0.4057
  ##      rep (Intercept) <NA> 0.1916 0.4378
  ## Residual        <NA> <NA> 0.1796 0.4238
  
  if(require("asreml", quietly=TRUE)) {
    libs(asreml,lucid)
    
    # AR1 x AR1 with linear row/col effects, random spline row/col
    dat <- transform(dat, xf=factor(col), yf=factor(row))
    dat <- dat[order(dat$xf,dat$yf),]
    m2 <- asreml(yield ~ gen + lin(yf) + lin(xf), data=dat,
                 random = ~ spl(yf) + spl(xf),
                 resid = ~ ar1(xf):ar1(yf))
    m2 <- update(m2) # More iterations
    
    # Scaling of spl components has changed in asreml from old versions
    lucid::vc(m2) # Match Burgueno p. 42
    ##       effect component std.error z.ratio bound 
    ##      spl(yf)  0.09077    0.08252   1.1       P 0
    ##      spl(xf)  0.08107    0.08209   0.99      P 0
    ##     xf:yf(R)  0.1482     0.03119   4.8       P 0
    ## xf:yf!xf!cor  0.1152     0.2269    0.51      U 0.1
    ## xf:yf!yf!cor  0.009467   0.2414    0.039     U 0.9
    
    plot(varioGram(m2), main="burgueno.rowcol")
  }
  

## End(Not run)

Field experiment with unreplicated genotypes plus one repeated check.

Description

Field experiment with unreplicated genotypes plus one repeated check.

Usage

data("burgueno.unreplicated")

Format

A data frame with 434 observations on the following 4 variables.

gen

genotype, 281 levels

col

column

row

row

yield

yield, tons/ha

Details

A field experiment with 280 new genotypes. A check genotype is planted in every 4th column.

The plot size is not given.

Electronic version of the data obtained from CropStat software.

Used with permission of Juan Burgueno.

Source

J Burgueno, A Cadena, J Crossa, M Banziger, A Gilmour, B Cullis (2000). User's guide for spatial analysis of field variety trials using ASREML. CIMMYT.

Examples

## Not run: 

  library(agridat)
  data(burgueno.unreplicated)
  dat <- burgueno.unreplicated

  # Define a 'check' variable for colors
  dat$check <- ifelse(dat$gen=="G000", 2, 1)
  # Every fourth column is the 'check' genotype
  libs(desplot)
  desplot(dat, yield ~ col*row,
          col=check, num=gen, #text=gen, cex=.3, # aspect unknown
          main="burgueno.unreplicated")

  if(require("asreml", quietly=TRUE)) {
    libs(asreml,lucid)

    # AR1 x AR1 with random genotypes
    dat <- transform(dat, xf=factor(col), yf=factor(row))
    dat <- dat[order(dat$xf,dat$yf),]
    m2 <- asreml(yield ~ 1, data=dat, random = ~ gen,
                 resid = ~ ar1(xf):ar1(yf))
    lucid::vc(m2)
    ##       effect component std.error z.ratio bound 
    ##          gen    0.9122   0.127       7.2     P 0  
    ##     xf:yf(R)    0.4993   0.05601     8.9     P 0  
    ## xf:yf!xf!cor   -0.2431   0.09156    -2.7     U 0  
    ## xf:yf!yf!cor    0.1255   0.07057     1.8     U 0.1
    
    # Note the strong saw-tooth pattern in the variogram.  Seems to
    # be column effects.
    plot(varioGram(m2), xlim=c(0,15), ylim=c(0,9), zlim=c(0,0.5),
         main="burgueno.unreplicated - AR1xAR1")
    # libs(lattice) # Show how odd columns are high
    # bwplot(resid(m2) ~ col, data=dat, horizontal=FALSE)
    
    # Define an even/odd column factor as fixed effect
    # dat$oddcol <- factor(dat$col 
    # The modulus operator throws a bug, so do it the hard way.
    dat$oddcol <- factor(dat$col - floor(dat$col / 2) *2 )
  
    m3 <- update(m2, yield ~ 1 + oddcol)
    m3$loglik # Matches Burgueno table 3, line 3
    
    plot(varioGram(m3), xlim=c(0,15), ylim=c(0,9), zlim=c(0,0.5),
         main="burgueno.unreplicated - AR1xAR1 + Even/Odd")
    # Much better-looking variogram
  }
  

## End(Not run)

Multi-environment trial of maize with pedigrees

Description

Maize yields in a multi-environment trial. Pedigree included.

Format

A data frame with 245 observations on the following 5 variables.

gen

genotype

male

male parent

female

female parent

env

environment

yield

yield, Mg/ha

Details

Ten inbreds were crossed to produce a diallel without reciprocals. The 45 F1 crosses were evaluated along with 4 checks in a triple-lattice 7x7 design. Pink stem borer infestation was natural.

Experiments were performed in 1995 and 1996 at three sites in northwestern Spain: Pontevedra (42 deg 24 min N, 8 deg 38 min W, 20 m over sea), Pontecaldelas (42 deg 23 N, 8 min 32 W, 300 m above sea), Ribadumia (42 deg 30 N, 8 min 46 W, 50 m above sea).

A two-letter location code and the year are concatenated to define the environment.

The average number of larvae per plant in each environment:

Used with permission of Ana Butron.

Source

Butron, A and Velasco, P and Ordas, A and Malvar, RA (2004). Yield evaluation of maize cultivars across environments with different levels of pink stem borer infestation. Crop Science, 44, 741-747. https://doi.org/10.2135/cropsci2004.7410

Examples

## Not run: 

  library(agridat)
  data(butron.maize)
  dat <- butron.maize

  libs(reshape2)
  mat <- acast(dat, gen~env, value.var='yield')
  mat <- sweep(mat, 2, colMeans(mat))
  mat.svd <- svd(mat)
  # Calculate PC1 and PC2 scores as in Table 4 of Butron
  # Comment out to keep Rcmd check from choking on '
  # round(mat.svd$u[,1:2] 

  biplot(princomp(mat), main="butron.maize", cex=.7) # Figure 1 of Butron


  if(require("asreml", quietly=TRUE)) {

    # Here we see if including pedigree information is helpful for a
    # multi-environment model
    # Including the pedigree provided little benefit
    
    # Create the pedigree
    ped <- dat[, c('gen','male','female')]
    ped <- ped[!duplicated(ped),] # remove duplicates
    unip <- unique(c(ped$male, ped$female)) # Unique parents
    unip <- unip[!is.na(unip)]
    # We have to define parents at the TOP of the pedigree
    ped <- rbind(data.frame(gen=c("Dent","Flint"), # genetic groups
                            male=c(0,0),
                            female=c(0,0)),
                 data.frame(gen=c("A509","A637","A661","CM105","EP28",
                                  "EP31","EP42","F7","PB60","Z77016"),
                            male=rep(c('Dent','Flint'),each=5),
                            female=rep(c('Dent','Flint'),each=5)),
                 ped)
    ped[is.na(ped$male),'male'] <- 0
    ped[is.na(ped$female),'female'] <- 0

    libs(asreml)
    ped.ainv <- ainverse(ped)
      
    m0 <- asreml(yield ~ 1+env, data=dat, random = ~ gen)
    m1 <- asreml(yield ~ 1+env, random = ~ vm(gen, ped.ainv), data=dat)
    m2 <- update(m1, random = ~ idv(env):vm(gen, ped.ainv))
    m3 <- update(m2, random = ~ diag(env):vm(gen, ped.ainv))
    m4 <- update(m3, random = ~ fa(env,1):vm(gen, ped.ainv))
    #summary(m0)$aic
    #summary(m4)$aic
    ##    df      AIC
    ## m0  2 229.4037
    ## m1  2 213.2487
    ## m2  2 290.6156
    ## m3  6 296.8061
    ## m4 11 218.1568
    
    p0 <- predict(m0, data=dat, classify="gen")$pvals
    p1 <- predict(m1, data=dat, classify="gen")$pvals
    p1par <- p1[1:12,]   # parents
    p1 <- p1[-c(1:12),]  # remove parents
    # Careful!  Need to manually sort the predictions
    p0 <- p0[order(as.character(p0$gen)),]
    p1 <- p1[order(as.character(p1$gen)),]
    
    # lims <- range(c(p0$pred, p1$pred)) * c(.95,1.05)
    lims <- c(6,8.25) # zoom in on the higher-yielding hybrids
    plot(p0$predicted.value, p1$predicted.value,
         pch="", xlim=lims, ylim=lims, main="butron.maize",
         xlab="BLUP w/o pedigree", ylab="BLUP with pedigree")
    abline(0,1,col="lightgray")
    text(x=p0$predicted.value, y=p1$predicted.value,
         p0$gen, cex=.5, srt=-45)
    text(x=min(lims), y=p1par$predicted.value, p1par$gen, cex=.5, col="red")
    round( cor(p0$predicted.value, p1$predicted.value), 3)
    # 0.994
    # Including the pedigree provided very little change
  }
  

## End(Not run)

Diameters of apples

Description

Measurements of the diameters of apples

Format

A data frame with 480 observations on the following 6 variables.

tree

tree, 10 levels

apple

apple, 24 levels

size

size of apple

appleid

unique id number for each apple

time

time period, 1-6 = (week/2)

diameter

diameter, inches

Details

Experiment conducted at the Winchester Agricultural Experiment Station of Virginia Polytechnic Institute and State University. Twentyfive apples were chosen from each of ten apple trees.

Of these, there were 80 apples in the largest size class, 2.75 inches in diameter or greater.

The diameters of the apples were recorded every two weeks over a 12-week period.

Source

Schabenberger, Oliver and Francis J. Pierce. 2002. Contemporary Statistical Models for the Plant and Soil Sciences. CRC Press, Boca Raton, FL.

Examples

## Not run: 
  
  library(agridat)
  data(byers.apple)
  dat <- byers.apple

  libs(lattice)
  xyplot(diameter ~ time | factor(appleid), data=dat, type=c('p','l'),
         strip=strip.custom(par.strip.text=list(cex=.7)),
         main="byers.apple")

  # Overall fixed linear trend, plus random intercept/slope deviations
  # for each apple.  Observations within each apple are correlated.
  libs(nlme)
  libs(lucid)
  m1 <- lme(diameter ~ 1 + time, data=dat,
            random = ~ time|appleid, method='ML',
            cor = corAR1(0, form=~ time|appleid),
            na.action=na.omit)
  vc(m1)
  ##       effect   variance   stddev corr
  ##  (Intercept) 0.007354   0.08575    NA
  ##         time 0.00003632 0.006027 0.83
  ##     Residual 0.0004555  0.02134    NA

## End(Not run)

Multi-environment trial of maize with fertilization

Description

Maize fertilization trial on Antigua and St. Vincent.

Format

A data frame with 612 observations on the following 7 variables.

isle

island, 2 levels

site

site

block

block

plot

plot, numeric

trt

treatment factor combining N,P,K

ears

number of ears harvested

yield

yield in kilograms

N

nitrogen fertilizer level

P

phosphorous fertilizer level

K

potassium fertilizer level

Details

Antigua is a coral island in the Caribbean with sufficient level land for experiments and a semi-arid climate, while St. Vincent is volcanic and level areas are uncommon, but the rainfall can be seasonally heavy.

There are 8-9 sites on each island.

Plots were 16 feet by 18 feet. A central area 12 feet by 12 feet was harvested and recorded.

The number of ears harvested was only recorded on the isle of Antigua.

The actual amounts of N, P, and K are not given. Only 0, 1, 2, 3.

The digits of the treatment represent the levels of nitrogen, phosphorus, and potassium fertilizer, respectively.

The TEAN site suffered damage from goats on plot 27, 35 and 36.

The LFAN site suffered damage from cattle on one boundary–plots 9, 18, 27, 36.

Electronic version of the data was retrieved from http://lib.stat.cmu.edu/datasets/Andrews/ https://www2.stat.duke.edu/courses/Spring01/sta114/data/andrews.html

Source

D.F. Andrews and A.M. Herzberg. 1985. Data: A Collection of Problems from Many Fields for the Student and Research Worker. Springer. Table 58.1 and 58.2.

References

Also in the DAAG package as data sets antigua and stVincent.

Examples

library(agridat)
data(caribbean.maize)
dat <- caribbean.maize

# Yield and ears are correlated
libs(lattice)
xyplot(yield~ears|site, dat, ylim=c(0,10), subset=isle=="Antigua",
       main="caribbean.maize - Antiqua")

# Some locs show large response to nitrogen (as expected), e.g. UISV, OOSV
dotplot(trt~yield|site, data=dat, main="caribbean.maize treatment response")

# Show the strong N*site interaction with little benefit on Antiqua, but
# a strong response on St.Vincent.
dat <- transform(dat, env=paste(substring(isle,1,1),site,sep="-"))
bwplot(yield~N|env, dat,
       main="caribbean.maize", xlab="nitrogen")

Germination of alfalfa seeds at various salt concentrations

Description

Germination of alfalfa seeds at various salt concentrations

Usage

data("carlson.germination")

Format

A data frame with 120 observations on the following 3 variables.

gen

genotype factor, 15 levels

germ

germination percent, 0-100

nacl

salt concentration percent, 0-2

Details

Data are means averaged over 5, 10, 15, and 20 day counts. Germination is expressed as a percent of the no-salt control to account for differences in germination among the cultivars.

Source

Carlson, JR and Ditterline, RL and Martin, JM and Sands, DC and Lund, RE. (1983). Alfalfa Seed Germination in Antibiotic Agar Containing NaCl. Crop science, 23, 882-885. https://doi.org/10.2135/cropsci1983.0011183X002300050016x

Examples

## Not run: 

library(agridat)
data(carlson.germination)
dat <- carlson.germination
dat$germ <- dat$germ/100 # Convert to percent

# Separate response curve for each genotype.
# Really, we should use a glmm with random int/slope for each genotype
m1 <- glm(germ~ 0 + gen*nacl, data=dat, family=quasibinomial)

# Plot data and fitted model
libs(latticeExtra)
newd <- data.frame(expand.grid(gen=levels(dat$gen), nacl=seq(0,2,length=100)))
newd$pred <- predict(m1, newd, type="response")
xyplot(germ~nacl|gen, dat, as.table=TRUE, main="carlson.germination",
       xlab="Percent NaCl", ylab="Fraction germinated") +
xyplot(pred~nacl|gen, newd, type='l', grid=list(h=1,v=0))


# Calculate LD50 values.  Note, Carlson et al used quadratics, not glm.
# MASS::dose.p cannot handle multiple slopes, so do a separate fit for
# each genotype.  Results are vaguely similar to Carlson table 5.
## libs(MASS)
## for(ii in unique(dat$gen)){
##   cat("\n", ii, "\n")
##   mm <- glm(germ ~ 1 + nacl, data=dat, subset=gen==ii, family=quasibinomial(link="probit"))
##   print(dose.p(mm))
## }
##              Dose         SE
## Anchor    1.445728  0.05750418
## Apollo    1.305804  0.04951644
## Baker     1.444153  0.07653989
## Drylander 1.351201  0.03111795
## Grimm     1.395735  0.04206377


## End(Not run)

Nonlinear maize yield-density model

Description

Nonlinear maize yield-density model.

Format

A data frame with 32 observations on the following 3 variables.

gen

genotype/hybrid, 8 levels

pop

population (plants)

yield

yield, pounds per hill

Details

Eight single-cross hybrids were in the experiment–Hy2xOh7 and WF9xC103 were included because it was believed they had optimum yields at relatively high and low populations. Planted in 1963. Plots were thinned to 2, 4, 6, 8 plants per hill, giving densities 8, 16, 24, 32 thousand plants per acre. Hills were in rows 40 inches apart. One hill = 1/4000 acre. Split-plot design with 5 reps, density is main plot and subplot was hybrid.

Source

S G Carmer and J A Jackobs (1965). An Exponential Model for Predicting Optimum Plant Density and Maximum Corn Yield. Agronomy Journal, 57, 241–244. https://doi.org/10.2134/agronj1965.00021962005700030003x

Examples

library(agridat)
data(carmer.density)
dat <- carmer.density
dat$gen <- factor(dat$gen, levels=c('Hy2x0h7','WF9xC103','R61x187-2',
                             'WF9x38-11','WF9xB14','C103xB14',
                             '0h43xB37','WF9xH60'))

# Separate analysis for each hybrid
# Model: y = x * a * k^x.  Table 1 of Carmer and Jackobs.
out <- data.frame(a=rep(NA,8), k=NA)
preds <- NULL
rownames(out) <- levels(dat$gen)
newdat <- data.frame(pop=seq(2,8,by=.1))
for(i in levels(dat$gen)){
  print(i)
  dati <- subset(dat, gen==i)
  mi <- nls(yield ~ pop * a * k^pop, data=dati, start=list(a=10,k=1))
  out[i, ] <- mi$m$getPars()
  # Predicted values
  pi <- cbind(gen=i, newdat, pred= predict(mi, newdat=newdat))
  preds <- rbind(preds, pi)
}
# Optimum plant density is -1/log(k)
out$pop.opt <- -1/log(out$k)
round(out, 3)
##               a     k pop.opt
## Hy2x0h7   0.782 0.865   6.875
## WF9xC103  1.039 0.825   5.192
## R61x187-2 0.998 0.798   4.441
## WF9x38-11 1.042 0.825   5.203
## WF9xB14   1.067 0.806   4.647
## C103xB14  0.813 0.860   6.653
## 0h43xB37  0.673 0.862   6.740
## WF9xH60   0.858 0.854   6.358


# Fit an overall fixed-effect with random deviations for each hybrid.
libs(nlme)
m1 <- nlme(yield ~ pop * a * k^pop,
           fixed = a + k ~ 1,
           random = a + k ~ 1|gen,
           data=dat, start=c(a=10,k=1))
# summary(m1) # Random effect for 'a' probably not needed


libs(latticeExtra)
# Plot Data, fixed-effect prediction, random-effect prediction.
pdat <- expand.grid(gen=levels(dat$gen), pop=seq(2,8,length=50))
pdat$pred <- predict(m1, pdat)
pdat$predf <- predict(m1, pdat, level=0)

xyplot(yield~pop|gen, dat, pch=16, as.table=TRUE,
       main="carmer.density models",
       key=simpleKey(text=c("Data", "Fixed effect","Random effect"),
         col=c("blue", "red","darkgreen"), columns=3, points=FALSE)) +
  xyplot(predf~pop|gen, pdat, type='l', as.table=TRUE, col="red") +
  xyplot(pred~pop|gen, pdat, type='l', col="darkgreen", lwd=2)

Relative cotton yield for different soil potassium concentrations

Description

Relative cotton yield for different soil potassium concentrations

Format

A data frame with 24 observations on the following 2 variables.

yield

Relative yield

potassium

Soil potassium, ppm

Details

Cate & Nelson used this data to determine the minimum optimal amount of soil potassium to achieve maximum yield.

Note, Fig 1 of Cate & Nelson does not match the data from Table 2. It sort of appears that points with high-concentrations of potassium were shifted left to a truncation point. Also, the calculations below do not quite match the results in Table 1. Perhaps the published data were rounded?

Source

Cate, R.B. and Nelson, L.A. (1971). A simple statistical procedure for partitioning soil test correlation data into two classes. Soil Science Society of America Journal, 35, 658–660. https://doi.org/10.2136/sssaj1971.03615995003500040048x

Examples

## Not run: 

library(agridat)
data(cate.potassium)
dat <- cate.potassium
names(dat) <- c('y','x')

CateNelson <- function(dat){
  dat <- dat[order(dat$x),] # Sort the data by x
  x <- dat$x
  y <- dat$y

  # Create a data.frame to store the results
  out <- data.frame(x=NA, mean1=NA, css1=NA, mean2=NA, css2=NA, r2=NA)

  css <- function(x) { var(x) * (length(x)-1) }
  tcss <- css(y) # Total corrected sum of squares

  for(i in 2:(length(y)-2)){
    y1 <- y[1:i]
    y2 <- y[-(1:i)]

    out[i, 'x'] <- x[i]
    out[i, 'mean1'] <- mean(y1)
    out[i, 'mean2'] <- mean(y2)
    out[i, 'css1'] <- css1 <- css(y1)
    out[i, 'css2'] <- css2 <- css(y2)
    out[i, 'r2'] <-  ( tcss - (css1+css2)) / tcss
  }
  return(out)
}

cn <- CateNelson(dat)
ix <- which.max(cn$r2)
with(dat, plot(y~x, ylim=c(0,110), xlab="Potassium", ylab="Yield"))
title("cate.potassium - Cate-Nelson analysis")
abline(v=dat$x[ix], col="skyblue")
abline(h=(dat$y[ix] + dat$y[ix+1])/2, col="skyblue")

  # another approach with similar results
  # https://joe.org/joe/2013october/tt1.php
  libs("rcompanion")
  cateNelson(dat$x, dat$y, plotit=0)

## End(Not run)

Factorial experiment of rice, 3x5x3x3

Description

Factorial experiment of rice, 3x5x3x3.

Usage

data("chakravertti.factorial")

Format

A data frame with 405 observations on the following 7 variables.

block

block/field

yield

yield

date

planting date, 5 levels

gen

genotype/variety, 3 levels

treat

treatment combination, 135 levels

seeds

number of seeds per hole, 3 levels

spacing

spacing, inches, 3 levels

Details

There were 4 treatment factors:

3 Genotypes (varieties): Nehara, Bhasamanik, Bhasakalma

5 Planting dates: Jul 16, Aug 1, Aug 16, Sep 1, Sep 16

3 Spacings: 6 in, 9 in, 12 inches

3 Seedlings per hole: 1, 2, local method

There were 3x5x3x3=135 treatment combinations. The experiment was divided in 3 blocks (fields). Total 405 plots.

"The plots of the same sowing date within each block were grouped together, and the position occupied by the sowing date groups within Within the blocks were determined at random. This grouping together of plots of the same sewing date was adopted to facilitate cultural operations. For the same reason, the three varieties were also laid out in compact rows. The nine combinations of spacings and seedling numbers were then thrown at random within each combination of date of planting and variety as shown in the diagram."

Note: The diagram appears to show the treatment combinations, NOT the physical layout.

Basically, date is a whole-plot effect, genotype is a sub-plot effect, and the 9 treatments (spacings * seedlings) are completely randomized withing the sub-plot effect.

Source

Chakravertti, S. C. and S. S. Bose and P. C. Mahalanobis (1936). A complex experiment on rice at the Chinsurah farm, Bengal, 1933-34. The Indian Journal of Agricultural Science, 6, 34-51. https://archive.org/details/in.ernet.dli.2015.271737/page/n83/mode/2up

References

None

Examples

## Not run: 

  libs(agridat)
  data(chakravertti.factorial)
  dat <- chakravertti.factorial
  
  # Simple means for each factor. Same as Chakravertti Table 3
  group_by(dat, gen) 
  group_by(dat, date) 
  group_by(dat, spacing) 
  group_by(dat, seeds) 

  libs(HH)
  interaction2wt(yield ~ gen + date + spacing + seeds, data=dat, main="chakravertti.factorial")

  # ANOVA matches Chakravertti table 2
  # This has a very interesting error structure.
  # block:date is error term for date
  # block:date:gen is error term for gen and date:gen
  # Residual is error term for all other tests (not needed inside Error())
  dat <- transform(dat,spacing=factor(spacing))
  m2 <- aov(yield ~ block + date + 
              gen + date:gen + 
              spacing + seeds +
              seeds:spacing + date:seeds + date:spacing + gen:seeds + gen:spacing +
              date:gen:seeds + date:gen:spacing + date:seeds:spacing + gen:seeds:spacing +
              date:gen:seeds:spacing + Error(block/(date + date:gen)),
            data=dat)
  summary(m2)
  

## End(Not run)

Fractional factorial of sugarcane, 1/3 3^5 = 3x3x3x3x3

Description

Fractional factorial of sugarcane, 1/3 3^5 = 3x3x3x3x3.

Usage

data("chinloy.fractionalfactorial")

Format

A data frame with 81 observations on the following 10 variables.

yield

yield

block

block

row

row position

col

column position

trt

treatment code

n

nitrogen treatment, 3 levels 0, 1, 2

p

phosphorous treatment, 3 levels 0, 1, 2

k

potassium treatment, 3 levels 0, 1, 2

b

bagasse treatment, 3 levels 0, 1, 2

m

filter press mud treatment, 3 levels 0, 1, 2

Details

An experiment grown in 1949 at the Worthy Park Estate in Jamaica.

There were 5 treatment factors:

3 Nitrogen levels: 0, 3, 6 hundred-weight per acre.

3 Phosphorous levels: 0, 4, 8 hundred-weight per acre.

3 Potassium (muriate of potash) levels: 0, 1, 2 hundred-weight per acre.

3 Bagasse (applied pre-plant) levels: 0, 20, 40 tons per acre.

3 Filter press mud (applied pre-plant) levels: 0, 10, 20 tons per acre.

Each plot was 18 yards long by 6 yards (3 rows) wide. Plots were arranged in nine columns of nine, a 2-yard space separating plots along the rows and two guard rows separating plots across the rows.

Field width: 6 yards * 9 plots + 4 yards * 8 gaps = 86 yards

Field length: 18 yards * 9 plots + 2 yards * 8 gaps = 178 yards

Source

T. Chinloy, R. F. Innes and D. J. Finney. (1953). An example of fractional replication in an experiment on sugar cane manuring. Journ Agricultural Science, 43, 1-11. https://doi.org/10.1017/S0021859600044567

References

None

Examples

## Not run: 

library(agridat)
data(chinloy.fractionalfactorial)
dat <- chinloy.fractionalfactorial

# Treatments are coded with levels 0,1,2. Make sure they are factors
dat <- transform(dat,
                 n=factor(n), p=factor(p), k=factor(k), b=factor(b), m=factor(m))

# Experiment layout
libs(desplot)
desplot(dat, yield ~ col*row,
        out1=block, text=trt, shorten="no", cex=0.6,
        aspect=178/86,
        main="chinloy.fractionalfactorial")

# Main effect and some two-way interactions. These match Chinloy table 6.
# Not sure how to code terms like p^2k=b^2m
m1 <- aov(yield ~ block + n + p + k + b + m + n:p + n:k + n:b + n:m, dat)
anova(m1)


## End(Not run)

Competition between varieties in cotton

Description

Competition between varieties in cotton, measurements taken for each row.

Usage

data("christidis.competition")

Format

A data frame with 270 observations on the following 8 variables.

plot

plot

plotrow

row within plot

block

block

row

row, only 1 row

col

column

gen

genotype

yield

yield, kg

height

height, cm

Details

Nine genotypes/varieties of cotton were used in a variety test. The plots were 100 meters long and 2.40 meters wide, each plot having 3 rows 0.80 meters apart.

The layout was an RCB of 5 blocks, each block having 2 replicates of every variety (with the original intention of trying 2 seed treatments). Each row was harvested/weighed separately. After the leaves of the plants had dried up and fallen, the mean height of each row was measured.

Christidis found significant competition between varieties, but not due to height differences. Crude analysis.

TODO: Find a better analysis of this data which incorporates field trends AND competition effects, maybe including a random effect for border rows of all genotype pairs (as neighbors)?

Source

Christidis, Basil G (1935). Intervarietal competition in yield trials with cotton. The Journal of Agricultural Science, 25, 231-237. Table 1. https://doi.org/10.1017/S0021859600009710

References

None

Examples

## Not run: 

library(agridat)
data(christidis.competition)
dat <- christidis.competition

# Match Christidis Table 2 means
# aggregate(yield ~ gen, aggregate(yield ~ gen+plot, dat, sum), mean)

# Each RCB block has 2 replicates of each genotype
# with(dat, table(block,gen))

libs(lattice)

# Tall plants yield more
# xyplot(yield ~ height|gen, data=dat)

# Huge yield variation across field. Also heterogeneous variance.
xyplot(yield ~ col, dat, group=gen, auto.key=list(columns=5),
       main="christidis.competition")


libs(mgcv)
if(is.element("package:gam", search())) detach("package:gam")
# Simple non-competition model to remove main effects
m1 <- gam(yield ~ gen + s(col), data=dat)
p1 <- as.data.frame(predict(m1, type="terms"))
names(p1) <- c('geneff','coleff')
dat2 <- cbind(dat, p1)
dat2 <- transform(dat2, res=yield-geneff-coleff)
libs(lattice)
xyplot(res ~  col, data=dat2, group=gen,
       main="christidis.competition - residuals")


## End(Not run)

Uniformity trial of cotton

Description

Uniformity trial of cotton in Greece, 1938

Usage

data("christidis.cotton.uniformity")

Format

A data frame with 1024 observations on the following 4 variables.

col

column

row

row

yield

yield, kg/unit

block

block factor

Details

The experiment was conducted in 1938 at Sindos by the Greek Cotton Research Institute.

Each block consisted of 20 rows, 1 meter apart and 66 meters long. Two rows on each side and 1 meter on each end were removed for borders. Each row was divided into 4 meter-lengths and harvested separately. There were 4 blocks, oriented at 0, 30, 60, 90 degrees.

Each block contained 16 rows, each 64 meters long.

Field width: 16 units * 4 m = 64 m

Field depth: 16 rows * 1 m = 16 m

Source

Christidis, B. G. (1939). Variability of Plots of Various Shapes as Affected by Plot Orientation. Empire Journal of Experimental Agriculture 7: 330-342. Table 1.

References

None

Examples

## Not run: 

library(agridat)
data(christidis.cotton.uniformity)
dat <- christidis.cotton.uniformity
  
# Match the mean yields in table 2. Not sure why '16' is needed
# sapply(split(dat$yield, dat$block), mean)*16
  
libs(desplot)
dat$yld <- dat$yield/4*1000 # re-scale to match Christidis fig 1
desplot(dat, yld ~ col*row|block,
        flip=TRUE, aspect=(16)/(64),
        main="christidis.cotton.uniformity")


## End(Not run)  

Uniformity trial of wheat

Description

Uniformity trial of wheat at Cambridge, UK in 1931.

Usage

data("christidis.wheat.uniformity")

Format

A data frame with 288 observations on the following 3 variables.

row

row

col

column

yield

yield

Details

Two blocks, 24 rows each. In block A, each 90-foot row was divided into 12 units, each unit 7.5 feet long. Rows were 8 inches wide.

Field width: 12 units * 7.5 feet = 90 feet

Field length: 24 rows * 8 inches = 16 feet

Source

Christidis, Basil G (1931). The importance of the shape of plots in field experimentation. The Journal of Agricultural Science, 21, 14-37. Table VI, p. 28. https://dx.doi.org/10.1017/S0021859600007942

References

None

Examples

## Not run: 

library(agridat)
data(christidis.wheat.uniformity)
dat <- christidis.wheat.uniformity
  
# sum(dat$yield) # Matches Christidis
  
 libs(desplot)
 desplot(dat, yield ~  col*row,
         flip=TRUE, aspect=16/90, # true aspect
         main="christidis.wheat.uniformity")


## End(Not run)

Soil resistivity in a field

Description

Soil resistivity in a field

Format

A data frame with 8641 observations on the following 5 variables.

northing

y ordinate

easting

x ordinate

resistivity

Soil resistivity, ohms

is.ns

Indicator of north/south track

track

Track number

Details

Resistivity is related to soil salinity.

Electronic version of the data was retrieved from http://lib.stat.cmu.edu/datasets/Andrews/

Cleaned version from Luke Tierney https://homepage.stat.uiowa.edu/~luke/classes/248/examples/soil

Source

William Cleveland, (1993). Visualizing Data.

Examples

## Not run: 

  library(agridat)
  data(cleveland.soil)
  dat <- cleveland.soil

  # Similar to Cleveland fig 4.64
  ## libs(latticeExtra)
  ## redblue <- colorRampPalette(c("firebrick", "lightgray", "#375997"))
  ## levelplot(resistivity ~ easting + northing, data = dat,
  ##           col.regions=redblue,
  ##           panel=panel.levelplot.points,
  ##           aspect=2.4, xlab= "Easting (km)", ylab= "Northing (km)",
  ##           main="cleveland")
  
  # 2D loess plot. Cleveland fig 4.68
  sg1 <- expand.grid(easting = seq(.15, 1.410, by = .02),
                     northing = seq(.150, 3.645, by = .02))
  lo1 <- loess(resistivity~easting*northing, data=dat, span = 0.1, degree = 2)
  fit1 <- predict(lo1, sg1)
  libs(lattice)
  redblue <- colorRampPalette(c("firebrick", "lightgray", "#375997"))
  levelplot(fit1 ~ sg1$easting * sg1$northing,
            col.regions=redblue,
            cuts = 9,
            aspect=2.4, xlab = "Easting (km)", ylab = "Northing (km)",
            main="cleveland.soil - 2D smooth of Resistivity")
  
  # 3D loess plot with data overlaid
  libs(rgl)
  bg3d(color = "white")
  clear3d()
  points3d(dat$easting, dat$northing, dat$resistivity / 100,
           col = rep("gray50", nrow(dat)))
  rgl::surface3d(seq(.15, 1.410, by = .02),
                 seq(.150, 3.645, by = .02),
                 fit1/100, alpha=0.9, col=rep("wheat", length(fit1)),
                 front="fill", back="fill")
  close3d()

## End(Not run)

Yield and number of plants in a sugarbeet fertilizer experiment

Description

Yield and number of plants in a sugarbeet fertilizer experiment.

Usage

data("cochran.beets")

Format

A data frame with 42 observations on the following 4 variables.

fert

fertilizer treatment

block

block

yield

yield, tons/acres

plants

number of plants per plot

Details

Yield (tons/acre) and number of beets per plot. Fertilizer treatments combine superphosphate (P), muriate of potash (K), and sodium nitrate (N).

Source

George Snedecor (1946). Statisitcal Methods, 4th ed. Table 12.13, p. 332.

References

H. Fairfield Smith (1957). Interpretation of Adjusted Treatment Means and Regressions in Analysis of Covariance. Biometrics, 13, 282-308. https://doi.org/10.2307/2527917

Examples

## Not run: 

library(agridat)
data(cochran.beets)
dat = cochran.beets

# P has strong effect
libs(lattice)
xyplot(yield ~ plants|fert, dat, main="cochran.beets") 


## End(Not run)

Multi-environment trial of corn, balanced incomplete block design

Description

Balanced incomplete block design in corn

Format

A data frame with 52 observations on the following 3 variables.

loc

location/block, 13 levels

gen

genotype/line, 13 levels

yield

yield, pounds/plot

Details

Incomplete block design. Each loc/block has 4 genotypes/lines. The blocks are planted at different locations.

Conducted in 1943 in North Carolina.

Source

North Carolina Agricultural Experiment Station, United States Department of Agriculture.

References

Cochran, W.G. and Cox, G.M. (1957), Experimental Designs, 2nd ed., Wiley and Sons, New York, p. 448.

Examples

## Not run: 

library(agridat)

data(cochran.bib)
dat <- cochran.bib

# Show the incomplete-block structure
libs(lattice)
redblue <- colorRampPalette(c("firebrick", "lightgray", "#375997"))
levelplot(yield~loc*gen, dat,
          col.regions=redblue,
          xlab="loc (block)", main="cochran.bib - incomplete blocks")

with(dat, table(gen,loc))
rowSums(as.matrix(with(dat, table(gen,loc))))
colSums(as.matrix(with(dat, table(gen,loc))))

m1 = aov(yield ~ gen + Error(loc), data=dat)
summary(m1)

libs(nlme)
m2 = lme(yield ~ -1 + gen, data=dat, random=~1|loc)


## End(Not run)

Potato scab infection with sulfur treatments

Description

Potato scab infection with sulfur treatments

Format

A data frame with 32 observations on the following 5 variables.

inf

infection percent

trt

treatment factor

row

row

col

column

Details

The experiment was conducted to investigate the effect of sulfur on controlling scab disease in potatoes. There were seven treatments. Control, plus spring and fall application of 300, 600, 1200 pounds/acre of sulfur. The response variable was infection as a percent of the surface area covered with scab. A completely randomized design was used with 8 replications of the control and 4 replications of the other treatments.

Although the original analysis did not show significant differences in the sulfur treatments, including a polynomial trend in the model uncovered significant differences (Tamura, 1988).

Source

W.G. Cochran and G. Cox, 1957. Experimental Designs, 2nd ed. John Wiley, New York.

References

Tamura, R.N. and Nelson, L.A. and Naderman, G.C., (1988). An investigation of the validity and usefulness of trend analysis for field plot data. Agronomy Journal, 80, 712-718.

https://doi.org/10.2134/agronj1988.00021962008000050003x

Examples

## Not run: 

library(agridat)
data(cochran.crd)
dat <- cochran.crd

# Field plan
libs(desplot)
desplot(dat, inf~col*row,
        text=trt, cex=1, # aspect unknown
        main="cochran.crd")

# CRD anova.  Table 6 of Tamura 1988
contrasts(dat$trt) <- cbind(c1=c(1,1,1,-6,1,1,1),   # Control vs Sulf
                            c2=c(-1,-1,-1,0,1,1,1)) # Fall vs Sp
m1 <- aov(inf ~ trt, data=dat)
anova(m1)
summary(m1, split=list(trt=list("Control vs Sulf"=1, "Fall vs Spring"=2)))

# Quadratic polynomial for columns...slightly different than Tamura 1988
m2 <- aov(inf ~ trt + poly(col,2), data=dat)
anova(m2)
summary(m2, split=list(trt=list("Control vs Sulf"=1, "Fall vs Spring"=2)))


## End(Not run)

Counts of eelworms before and after fumigant treatments

Description

Counts of eelworms before and after fumigant treatments

Format

A data frame with 48 observations on the following 7 variables.

block

block factor, 4 levels

row

row

col

column

fumigant

fumigant factor

dose

dose, Numeric 0,1,2. Maybe should be a factor?

initial

count of eelworms pre-treatment

final

count of eelworms post-treatment

grain

grain yield in pounds

straw

straw yield in pounds

weeds

ratio of weeds to total oats

Details

A soil fumigation experiment on Spring Oats, conducted in 1935.

Each plot is 30 links x 41.7 links, but it is not clear which side of the plot has a specific length.

Treatment codes: Con = Control, Chl = Chlorodinitrobenzen, Cym = Cymag, Car = Carbon Disulphide jelly, See = Seekay.

Experiment was conducted in 1935 at Rothamsted Experiment Station. In early March 400 grams of soil (4 x 100g) were sampled and the number of eelworm cysts were counted. Fumigants were added to the soil, oats were sown and later harvested. In October, the plots were again sampled and the final count of cysts recorded.

The Rothamsted report concludes that "Car" and "Cym" produced higher yields, due partly to the nitrogen in the fumigant, while "Chl" decreased the yield. All fumigants reduced weeds. The crop was 'unusually weedy'. "Car" and "See" decreased the number of eelworm cysts in the soil.

The original data can be found in the Rothamsted Report. The report notes the position of the blocks in the field were slightly different than shown.

The experiment plan shown in Bailey (2008, p. 73), shows columns 9-11 shifted slightly upward. It is not clear why.

Thanks to U.Genschel for identifying a typo.

Source

Cochran and Cox, 1950. Experimental Designs. Table 3.1.

References

R. A. Bailey (2008). Design of Comparative Experiments. Cambridge.

Other Experiments at Rothamsted (1936). Report For 1935, Rothamsted Research. pp 174 - 193. https://doi.org/10.23637/ERADOC-1-67

Examples

## Not run: 

  library(agridat)
  data(cochran.eelworms)
  dat <- cochran.eelworms

  libs(lattice)
  splom(dat[ , 5:10],
        group=dat$fumigant, auto.key=TRUE,
        main="cochran.eelworms")
  
  libs(desplot)
  desplot(dat, fumigant~col*row, text=dose, flip=TRUE, cex=2)
  
  # Very strong spatial trends
  desplot(dat, initial ~ col*row,
          flip=TRUE, # aspect unknown
          main="cochran.eelworms")


  # final counts are strongly related to initial counts
  libs(lattice)
  xyplot(final~initial|factor(dose), data=dat, group=fumigant,
         main="cochran.eelworms - by dose (panel) & fumigant",
         xlab="Initial worm count",
         ylab="Final worm count", auto.key=list(columns=5))
  
  # One approach...log transform, use 'initial' as covariate, create 9 treatments
  dat <- transform(dat, trt=factor(paste0(fumigant, dose)))
  m1 <- aov(log(final) ~ block + trt + log(initial), data=dat)
  anova(m1)


## End(Not run)

Factorial experiment of beans, 2x2x2x2

Description

Factorial experiment of beans, 2x2x2x2.

Usage

data("cochran.factorial")

Format

A data frame with 32 observations on the following 4 variables.

rep

rep factor

block

block factor

trt

treatment factor, 16 levels

yield

yield (pounds)

d

dung treatment, 2 levels

n

nitrogen treatment, 2 levels

p

phosphorous treatment, 2 levels

k

potassium treatment, 2 levels

Details

Conducted by Rothamsted Experiment Station in 1936.

There were 4 treatment factors:

2 d dung levels: None, 10 tons/acre.

2 n nitrochalk levels: None, 0.4 hundredweight nitrogen per acre.

2 p superphosphate levels: None, 0.6 hundredweight per acre

2 k muriate of potash levels: None, 1 hundredweight K20 per acres.

The response variable is the yield of beans.

Source

Cochran, W.G. and Cox, G.M. (1957), Experimental Designs, 2nd ed., Wiley and Sons, New York, p. 160.

Examples

## Not run: 

library(agridat)
data(cochran.factorial)
dat <- cochran.factorial

# Ensure factors
dat <- transform(dat, d=factor(d), n=factor(n), p=factor(p), k=factor(k))

# Cochran table 6.5.
m1 <- lm(yield ~ rep * block + (d+n+p+k)^3, data=dat)
anova(m1)


libs(FrF2)
aliases(m1)
MEPlot(m1, select=3:6,
       main="cochran.factorial - main effects plot")


## End(Not run)

Latin square design in wheat

Description

Six wheat plots were sampled by six operators and shoot heights measured. The operators sampled plots in six ordered sequences. The dependent variate was the difference between measured height and true height of the plot.

Format

A data frame with 36 observations on the following 4 variables.

row

row

col

column

operator

operator factor

diff

difference between measured height and true height

Source

Cochran, W.G. and Cox, G.M. (1957), Experimental Designs, 2nd ed., Wiley and Sons, New York.

Examples

## Not run: 

library(agridat)
data(cochran.latin)
dat <- cochran.latin

libs(desplot)
desplot(dat, diff~col*row,
        text=operator, cex=1, # aspect unknown
        main="cochran.latin")


dat <- transform(dat, rf=factor(row), cf=factor(col))
aov.dat <- aov(diff ~ operator + Error(rf*cf), dat)
summary(aov.dat)
model.tables(aov.dat, type="means")

## End(Not run)

Balanced lattice experiment in cotton

Description

Balanced lattice experiment in cotton

Usage

data("cochran.lattice")

Format

A data frame with 80 observations on the following 5 variables.

y

percent of affected flower buds

rep

replicate

row

row

col

column

trt

treatment factor

Details

The experiment is a balanced lattice square with 16 treatments in a 4x4 layout in each of 5 replicates. The treatments were applied to cotton plants. Each plot was ten rows wide by 70 feet long (about 1/18 of an acre). (Estimated plot width is 34.5 feet.) Data were collected from the middle 4 rows. The data are the percentages of squares showing attack by boll weevils. A 'square' is the name given to a young flower bud.

The plot orientation is not clear.

Source

William G. Cochran, Gertrude M. Cox. Experimental Designs, 2nd Edition. Page 490.

Originally from: F. M. Wadley (1946). Incomplete block designs in insect population problems. J. Economic Entomology, 38, 651–654.

References

Walter Federer. Combining Standard Block Analyses With Spatial Analyses Under a Random Effects Model. Cornell Univ Tech Report BU-1373-MA. https://hdl.handle.net/1813/31971

Examples

## Not run: 

library(agridat)
data(cochran.lattice)
dat <- cochran.lattice

libs(desplot)
desplot(dat, y~row*col|rep,
        text=trt, # aspect unknown, should be 2 or .5
         main="cochran.lattice")


# Random rep,row,column model often used by Federer
libs(lme4)
dat <- transform(dat, rowf=factor(row), colf=factor(col))
m1 <-  lmer(y ~ trt + (1|rep) + (1|rep:row) + (1|rep:col), data=dat)
summary(m1)


## End(Not run)

Wireworms controlled by fumigants in a latin square

Description

Wireworms controlled by fumigants in a latin square

Format

A data frame with 25 observations on the following 4 variables.

row

row

col

column

trt

fumigant treatment, 5 levels

worms

count of wireworms per plot

Details

Plots were approximately 22 cm by 13 cm. Layout of the experiment was a latin square. The number of wireworms in each plot was counted, following soil fumigation the previous year.

Source

W. G. Cochran (1938). Some difficulties in the statistical analysis of replicated experiments. Empire Journal of Experimental Agriculture, 6, 157–175.

References

Ron Snee (1980). Graphical Display of Means. The American Statistician, 34, 195-199. https://www.jstor.org/stable/2684060 https://doi.org/10.1080/00031305.1980.10483028

W. Cochran (1940). The analysis of variance when experimental errors follow the Poisson or binomial laws. The Annals of Mathematical Statistics, 11, 335-347. https://www.jstor.org/stable/2235680

G W Snedecor and W G Cochran, 1980. Statistical Methods, Iowa State University Press. Page 288.

Examples

## Not run: 

library(agridat)
data(cochran.wireworms)
dat <- cochran.wireworms

libs(desplot)
desplot(dat, worms ~ col*row,
        text=trt, cex=1, # aspect unknown
        main="cochran.wireworms")

# Trt K is effective, but not the others.  Really, this says it all.
libs(lattice)
bwplot(worms ~ trt, dat, main="cochran.wireworms", xlab="Treatment")

# Snedecor and Cochran do ANOVA on sqrt(x+1).
dat <- transform(dat, rowf=factor(row), colf=factor(col))
m1 <- aov(sqrt(worms+1) ~ rowf + colf + trt, data=dat)
anova(m1)

# Instead of transforming, use glm
m2 <- glm(worms ~ trt + rowf + colf, data=dat, family="poisson")
anova(m2)

# GLM with random blocking.
libs(lme4)
m3 <- glmer(worms ~ -1 +trt +(1|rowf) +(1|colf), data=dat, family="poisson")
summary(m3)
## Fixed effects:
##      Estimate Std. Error z value Pr(>|z|)    
## trtK   0.1393     0.4275   0.326    0.745    
## trtM   1.7814     0.2226   8.002 1.22e-15 ***
## trtN   1.9028     0.2142   8.881  < 2e-16 ***
## trtO   1.7147     0.2275   7.537 4.80e-14 ***


## End(Not run)

Potato yields in single-drill plots

Description

Potato yields in single-drill plots

Usage

data("connolly.potato")

Format

A data frame with 80 observations on the following 6 variables.

rep

block

gen

variety

row

row

col

column

yield

yield, kg/ha

matur

maturity group

Details

Connolly et el use this data to illustrate how yield can be affected by competition from neighboring plots.

This data uses M1, M2, M3 for maturity, while Connolly et al use FE (first early), SE (second early) and M (maincrop).

The trial was 20 sections, each of which was an independent row of 20 drills. The data here are four reps of single-drill plots from sections 1, 6, 11, and 16.

The neighbor covariate for a plot is defined as the average of the plots to the left and right. For drills at the edge of the trial, the covariate was the average of the one neighboring plot yield and the section (i.e. rep) mean.

It would be interesting to fit a model that uses differences in maturity between a plot and its neighbor as the actual covariate.

https://doi.org/10.1111/j.1744-7348.1993.tb04099.x

Used with permission of Iain Currie.

Source

Connolly, T and Currie, ID and Bradshaw, JE and McNicol, JW. (1993). Inter-plot competition in yield trials of potatoes Solanum tuberosum L. with single-drill plots. Annals of Applied Biology, 123, 367-377.

Examples

library(agridat)
data(connolly.potato)
dat <- connolly.potato

# Field plan
libs(desplot)
desplot(dat, yield~col*row,
        out1=rep, # aspect unknown
        main="connolly.potato yields (reps not contiguous)")


# Later maturities are higher yielding
libs(lattice)
bwplot(yield~matur, dat, main="connolly.potato yield by maturity")

# Observed raw means. Matches Connolly table 2.
mn <- aggregate(yield~gen, data=dat, FUN=mean)
mn[rev(order(mn$yield)),]

# Create a covariate which is the average of neighboring plot yields
libs(reshape2)
mat <- acast(dat, row~col, value.var='yield')
mat2 <- matrix(NA, nrow=4, ncol=20)
mat2[,2:19] <- (mat[ , 1:18] + mat[ , 3:20])/2
mat2[ , 1] <- (mat[ , 1] + apply(mat, 1, mean))/2
mat2[ , 20] <- (mat[ , 20] + apply(mat, 1, mean))/2
dat2 <- melt(mat2)
colnames(dat2) <- c('row','col','cov')
dat <- merge(dat, dat2)
# xyplot(yield ~ cov, data=dat, type=c('p','r'))

# Connolly et al fit a model with avg neighbor yield as a covariate
m1 <- lm(yield ~ 0 + gen + rep + cov, data=dat)
coef(m1)['cov'] # = -.303  (Connolly obtained -.31)

# Block names and effects
bnm <- c("R1","R2","R3","R4")
beff <- c(0, coef(m1)[c('repR2','repR3','repR4')])
# Variety names and effects
vnm <- paste0("V", formatC(1:20, width=2, flag='0'))
veff <- coef(m1)[1:20]

# Adjust yield for variety and block effects
dat <- transform(dat, yadj = yield - beff[match(rep,bnm)]
                - veff[match(gen,vnm)])

# Similar to Connolly Fig 1.  Point pattern doesn't quite match
xyplot(yadj~cov, data=dat, type=c('p','r'),
       main="connolly.potato",
       xlab="Avg yield of nearest neighbors",
       ylab="Yield, adjusted for variety and block effects")

Uniformity trial of rice in Malaysia

Description

Uniformity trial of rice in Malaysia

Usage

data("coombs.rice.uniformity")

Format

A data frame with 54 observations on the following 3 variables.

row

row

col

column

yield

yield in gantangs per plot

Details

Estimated harvest date is 1915 or earlier.

Field length, 18 plots * 1/2 chain.

Field width, 3 plots * 1/2 chain.

Source

Coombs, G. E. and J. Grantham (1916). Field Experiments and the Interpretation of their results. The Agriculture Bulletin of the Federated Malay States, No 7. https://www.google.com/books/edition/The_Agricultural_Bulletin_of_the_Federat/M2E4AQAAMAAJ

References

None

Examples

## Not run: 
  library(agridat)
  data(coombs.rice.uniformity)
  dat <- coombs.rice.uniformity

  # Data check. Matches Coombs 709.4
  # sum(dat$yield)

  # There are an excess number of 12s and 14s in the yield
  libs(lattice)
  qqmath( ~ yield, dat) # weird
  
  libs(desplot)
  desplot(dat, yield ~ col*row,
          main="coombs.rice.uniformity",
          flip=TRUE, aspect=(18 / 3))

## End(Not run)

Multi-environment trial of maize for 9 cultivars at 20 locations.

Description

Maize yields for 9 cultivars at 20 locations.

Usage

data("cornelius.maize")

Format

A data frame with 180 observations on the following 3 variables.

env

environment factor, 20 levels

gen

genotype/cultivar, 9 levels

yield

yield, kg/ha

Details

Cell means (kg/hectare) for the CIMMYT EVT16B maize yield trial.

Source

P L Cornelius and J Crossa and M S Seyedsadr. (1996). Statistical Tests and Estimators of Multiplicative Models for Genotype-by-Environment Interaction. Book: Genotype-by-Environment Interaction. Pages 199-234.

References

Forkman, Johannes and Piepho, Hans-Peter. (2014). Parametric bootstrap methods for testing multiplicative terms in GGE and AMMI models. Biometrics, 70(3), 639-647. https://doi.org/10.1111/biom.12162

Examples

## Not run: 

library(agridat)
data(cornelius.maize)
dat <- cornelius.maize

# dotplot(gen~yield|env,dat) # We cannot compare genotype yields easily
# Subtract environment mean from each observation
libs(reshape2)
mat <- acast(dat, gen~env)
mat <- scale(mat, scale=FALSE)
dat2 <- melt(mat)
names(dat2) <- c('gen','env','yield')
libs(lattice)
bwplot(yield ~ gen, dat2,
       main="cornelius.maize - environment centered yields")

if(0){
# This reproduces the analysis of Forkman and Piepho.

test.pc <- function(Y0, type="AMMI", n.boot=10000, maxpc=6) {

  # Test the significance of Principal Components in GGE/AMMI

  # Singular value decomposition of centered/double-centered Y
  Y <- sweep(Y0, 1, rowMeans(Y0)) # subtract environment means
  if(type=="AMMI") {
    Y <- sweep(Y, 2, colMeans(Y0)) # subtract genotype means
    Y <- Y + mean(Y0)
  }
  lam <- svd(Y)$d

  # Observed value of test statistic.
  # t.obs[k] is the proportion of variance explained by the kth term out of
  # the k...M terms, e.g. t.obs[2] is lam[2]^2 / sum(lam[2:M]^2)
  t.obs <- { lam^2/rev(cumsum(rev(lam^2))) } [1:(M-1)]
  t.boot <- matrix(NA, nrow=n.boot, ncol=M-1)

  # Centering rows/columns reduces the rank by 1 in each direction.
  I <- if(type=="AMMI") nrow(Y0)-1 else nrow(Y0)
  J <- ncol(Y0)-1
  M <- min(I, J) # rank of Y, maximum number of components
  M <- min(M, maxpc) # Optional step: No more than 5 components

  for(K in 0:(M-2)){ # 'K' multiplicative components in the svd

    for(bb in 1:n.boot){
      E.b <- matrix(rnorm((I-K) * (J-K)), nrow = I-K, ncol = J-K)
      lam.b <- svd(E.b)$d
      t.boot[bb, K+1] <- lam.b[1]^2 / sum(lam.b^2)
    }

  }

  # P-value for each additional multiplicative term in the SVD.
  # P-value is the proportion of time bootstrap values exceed t.obs
  colMeans(t.boot > matrix(rep(t.obs, n.boot), nrow=n.boot, byrow=TRUE))
}

dat <- cornelius.maize

# Convert to matrix format
libs(reshape2)
dat <- acast(dat, env~gen, value.var='yield')

## R> test.pc(dat,"AMMI")
## [1] 0.0000 0.1505 0.2659 0.0456 0.1086 # Forkman: .00 .156 .272 .046 .111

## R> test.pc(dat,"GGE")
## [1] 0.0000 0.2934 0.1513 0.0461 0.2817 # Forkman: .00 .296 .148 .047 .285

}

## End(Not run)

Multi-environment trial of corn

Description

The data is the yield (kg/acre) of 20 genotypes of corn at 7 locations.

Format

A data frame with 140 observations on the following 3 variables.

gen

genotype, 20 levels

loc

location, 7 levels

yield

yield, kg/acre

Details

The data is used by Corsten & Denis (1990) to illustrate two-way clustering by minimizing the interaction sum of squares.

In their paper, the labels on the location dendrogram have a slight typo. The order of the loc labels shown is 1 2 3 4 5 6 7. The correct order of the loc labels is 1 2 4 5 6 7 3.

Used with permission of Jean-Baptiste Denis.

Source

L C A Corsten and J B Denis, (1990). Structuring Interaction in Two-Way Tables By Clustering. Biometrics, 46, 207–215. Table 1. https://doi.org/10.2307/2531644

Examples

## Not run: 

library(agridat)
data(corsten.interaction)
dat <- corsten.interaction

libs(reshape2)
m1 <- melt(dat, measure.var='yield')
dmat <- acast(m1, loc~gen)

# Corsten (1990) uses this data to illustrate simultaneous row and
# column clustering based on interaction sums-of-squares.
# There is no (known) function in R to reproduce this analysis
# (please contact the package maintainer if this is not true).
# For comparison, the 'heatmap' function clusters the rows and
# columns _independently_ of each other.
heatmap(dmat, main="corsten.interaction")


## End(Not run)

Strip-split-plot of barley with fertilizer, calcium, and soil factors.

Description

Strip-split-plot of barley with fertilizer, calcium, and soil factors.

Format

A data frame with 96 observations on the following 5 variables.

rep

replicate, 4 levels

soil

soil, 3 levels

fert

fertilizer, 4 levels

calcium

calcium, 2 levels

yield

yield of winter barley

Details

Four different fertilizer treatments are laid out in vertical strips, which are then split into subplots with different levels of calcium. Soil type is stripped across the split-plot experiment, and the entire experiment is then replicated three times.

Sometimes called a split-block design.

Source

Comes from the notes of Gertrude Cox and A. Rotti.

References

SAS/STAT(R) 9.2 User's Guide, Second Edition. Example 23.5 Strip-Split Plot. https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_anova_sect030.htm

Examples

## Not run: 

library(agridat)
data(cox.stripsplit)
dat <- cox.stripsplit

# Raw means
# aggregate(yield ~ calcium, data=dat, mean)
# aggregate(yield ~ soil, data=dat, mean)
# aggregate(yield ~ calcium, data=dat, mean)

libs(HH)
interaction2wt(yield ~ rep + soil + fert + calcium, dat,
               x.between=0, y.between=0,
               main="cox.stripsplit")

# Traditional AOV
m1 <- aov(yield~ fert*calcium*soil +
          Error(rep/(fert+soil+calcium:fert+soil:fert)),
          data=dat)
summary(m1)

# With balanced data, the following are all basically identical

libs(lme4)
# The 'rep:soil:fert' term causes problems...so we drop it.
m2 <- lmer(yield ~ fert*soil*calcium + (1|rep) + (1|rep:fert) +
             (1|rep:soil) + (1|rep:fert:calcium), data=dat)


if(0){
  # afex uses Kenword-Rogers approach for denominator d.f.
  libs(afex)
  mixed(yield ~ fert*soil*calcium + (1|rep) + (1|rep:fert) +
          (1|rep:soil) + (1|rep:fert:calcium) + (1|rep:soil:fert), data=dat,
        control=lmerControl(check.nobs.vs.rankZ="ignore"))
  ##              Effect      stat ndf     ddf F.scaling p.value
  ## 1       (Intercept) 1350.8113   1  3.0009         1  0.0000
  ## 2              fert    3.5619   3  9.0000         1  0.0604
  ## 3              soil    3.4659   2  6.0000         1  0.0999
  ## 4           calcium    1.8835   1 12.0000         1  0.1950
  ## 5         fert:soil    1.2735   6 18.0000         1  0.3179
  ## 6      fert:calcium    4.4457   3 12.0000         1  0.0255
  ## 7      soil:calcium    0.2494   2 24.0000         1  0.7813
  ## 8 fert:soil:calcium    0.3504   6 24.0000         1  0.9027
}


## End(Not run)

Cucumber yields and quantitative traits

Description

Cucumber yields and quantitative traits

Usage

data("cramer.cucumber")

Format

A data frame with 24 observations on the following 9 variables.

cycle

cycle

rep

replicate

plants

plants per plot

flowers

number of pistillate flowers

branches

number of branches

leaves

number of leaves

totalfruit

total fruit number

culledfruit

culled fruit number

earlyfruit

early fruit number

Details

The data are used to illustrate path analysis of the correlations between phenotypic traits.

Used with permission of Christopher Cramer.

Source

Christopher S. Cramer, Todd C. Wehner, and Sandra B. Donaghy. 1999. Path Coefficient Analysis of Quantitative Traits. In: Handbook of Formulas and Software for Plant Geneticists and Breeders, page 89.

References

Cramer, C. S., T. C. Wehner, and S. B. Donaghy. 1999. PATHSAS: a SAS computer program for path coefficient analysis of quantitative data. J. Hered, 90, 260-262 https://doi.org/10.1093/jhered/90.1.260

Examples

## Not run: 

  library(agridat)
  data(cramer.cucumber)
  dat <- cramer.cucumber

  libs(lattice)
  splom(dat[3:9], group=dat$cycle,
        main="cramer.cucumber - traits by cycle",
        auto.key=list(columns=3))


  # derived traits
  dat <- transform(dat,
                   marketable = totalfruit-culledfruit,
                   branchesperplant = branches/plants,
                   nodesperbranch = leaves/(branches+plants),
                   femalenodes = flowers+totalfruit)
  dat <- transform(dat,
                   perfenod = (femalenodes/leaves),
                   fruitset = totalfruit/flowers,
                   fruitperplant = totalfruit / plants,
                   marketableperplant = marketable/plants,
                   earlyperplant=earlyfruit/plants)
  # just use cycle 1
  dat1 <- subset(dat, cycle==1)

  # define independent and dependent variables
  indep <- c("branchesperplant", "nodesperbranch", "perfenod", "fruitset")
  dep0 <- "fruitperplant"
  dep <- c("marketable","earlyperplant")

  # standardize trait data for cycle 1
  sdat <- data.frame(scale(dat1[1:8, c(indep,dep0,dep)]))

  # slopes for dep0 ~ indep
  X <- as.matrix(sdat[,indep])
  Y <- as.matrix(sdat[,c(dep0)])
  # estdep <- solve(t(X) 
  estdep <- solve(crossprod(X), crossprod(X,Y))
  estdep
  ## branchesperplant 0.7160269
  ## nodesperbranch   0.3415537
  ## perfenod         0.2316693
  ## fruitset         0.2985557

  # slopes for dep ~ dep0
  X <- as.matrix(sdat[,dep0])
  Y <- as.matrix(sdat[,c(dep)])
  # estind2 <- solve(t(X) 
  estind2 <- solve(crossprod(X), crossprod(X,Y))
  estind2
  ##  marketable earlyperplant
  ##     0.97196     0.8828393

  # correlation coefficients for indep variables
  corrind=cor(sdat[,indep])
  round(corrind,2)
  ##                  branchesperplant nodesperbranch perfenod fruitset
  ## branchesperplant             1.00           0.52    -0.24     0.09
  ## nodesperbranch               0.52           1.00    -0.44     0.14
  ## perfenod                    -0.24          -0.44     1.00     0.04
  ## fruitset                     0.09           0.14     0.04     1.00

  # Correlation coefficients for dependent variables
  corrdep=cor(sdat[,c(dep0, dep)])
  round(corrdep,2)
  ##               fruitperplant marketable earlyperplant
  ## fruitperplant          1.00       0.97          0.88
  ## marketable             0.97       1.00          0.96
  ## earlyperplant          0.88       0.96          1.00

  result = corrind
  result = result*matrix(estdep,ncol=4,nrow=4,byrow=TRUE)
  round(result,2) # match SAS output columns 1-4
  ##                  branchesperplant nodesperbranch perfenod fruitset
  ## branchesperplant             0.72           0.18    -0.06     0.03
  ## nodesperbranch               0.37           0.34    -0.10     0.04
  ## perfenod                    -0.17          -0.15     0.23     0.01
  ## fruitset                     0.07           0.05     0.01     0.30

  resdep0 = rowSums(result)
  resdep <- cbind(resdep0,resdep0)*matrix(estind2, nrow=4,ncol=2,byrow=TRUE)
  colnames(resdep) <- dep
  # slightly different from SAS output last 2 columns
  round(cbind(fruitperplant=resdep0, round(resdep,2)),2)
  ##                  fruitperplant marketable earlyperplant
  ## branchesperplant          0.87       0.84          0.76
  ## nodesperbranch            0.65       0.63          0.58
  ## perfenod                 -0.08      -0.08         -0.07
  ## fruitset                  0.42       0.41          0.37

## End(Not run)

Weight gain in pigs for different treatments

Description

Weight gain in pigs for different treatments, with initial weight and feed eaten as covariates.

Usage

data("crampton.pig")

Format

A data frame with 50 observations on the following 5 variables.

treatment

feed treatment

rep

replicate

weight1

initial weight

feed

feed eaten

weight2

final weight

Details

A study of the effect of initial weight and feed eaten on the weight gaining ability of pigs with different feed treatments.

The data are extracted from Ostle. It is not clear that 'replicate' is actually a blocking replicate as opposed to a repeated measurement. The original source document needs to be consulted.

Source

Crampton, EW and Hopkins, JW. (1934). The Use of the Method of Partial Regression in the Analysis of Comparative Feeding Trial Data, Part II. The Journal of Nutrition, 8, 113-123. https://doi.org/10.1093/jn/8.3.329

References

Bernard Ostle. Statistics in Research, Page 458. https://archive.org/details/secondeditionsta001000mbp

Goulden (1939). Methods of Statistical Analysis, 1st ed. Page 256-259. https://archive.org/details/methodsofstatist031744mbp

Examples

## Not run: 
  
  library(agridat)

  data(crampton.pig)
  dat <- crampton.pig

  dat <- transform(dat, gain=weight2-weight1)
  libs(lattice)
  # Trt 4 looks best
  xyplot(gain ~ feed, dat, group=treatment, type=c('p','r'),
         auto.key=list(columns=5),
         xlab="Feed eaten", ylab="Weight gain", main="crampton.pig")
  
  # Basic Anova without covariates
  m1 <- lm(weight2 ~ treatment + rep, data=dat)
  anova(m1)
  # Add covariates
  m2 <- lm(weight2 ~ treatment + rep + weight1 + feed, data=dat)
  anova(m2)
  # Remove treatment, test this nested model for significant treatments
  m3 <- lm(weight2 ~ rep + weight1 + feed, data=dat)
  anova(m2,m3) # p-value .07. F=2.34 matches Ostle

## End(Not run)

Multi-environment trial of wheat for 18 genotypes at 25 locations

Description

Wheat yields for 18 genotypes at 25 locations

Format

A data frame with 450 observations on the following 3 variables.

loc

location

locgroup

location group: Grp1-Grp2

gen

genotype

gengroup

genotype group: W1, W2, W3

yield

grain yield, tons/ha

Details

Grain yield from the 8th Elite Selection Wheat Yield Trial to evaluate 18 bread wheat genotypes at 25 locations in 15 countries.

Cross et al. used this data to cluster loctions into 2 mega-environments and clustered genotypes into 3 wheat clusters.

Locations

Used with permission of Jose' Crossa.

Source

Crossa, J and Fox, PN and Pfeiffer, WH and Rajaram, S and Gauch Jr, HG. (1991). AMMI adjustment for statistical analysis of an international wheat yield trial. Theoretical and Applied Genetics, 81, 27–37. https://doi.org/10.1007/BF00226108

References

Jean-Louis Laffont, Kevin Wright and Mohamed Hanafi (2013). Genotype + Genotype x Block of Environments (GGB) Biplots. Crop Science, 53, 2332-2341. https://doi.org/10.2135/cropsci2013.03.0178

Examples

## Not run: 

  library(agridat)
  data(crossa.wheat)
  dat <- crossa.wheat
  
  # AMMI biplot.  Fig 3 of Crossa et al.
  libs(agricolae)
  m1 <- with(dat, AMMI(E=loc, G=gen, R=1, Y=yield))
  b1 <- m1$biplot[,1:4]
  b1$PC1 <- -1 * b1$PC1 # Flip vertical
  plot(b1$yield, b1$PC1, cex=0.0,
       text(b1$yield, b1$PC1, cex=.5, labels=row.names(b1),col="brown"),
       main="crossa.wheat AMMI biplot",
       xlab="Average yield", ylab="PC1", frame=TRUE)
  mn <- mean(b1$yield)
  abline(h=0, v=mn, col='wheat')

  g1 <- subset(b1,type=="GEN")
  text(g1$yield, g1$PC1, rownames(g1), col="darkgreen", cex=.5)
  
  e1 <- subset(b1,type=="ENV")
  arrows(mn, 0,
         0.95*(e1$yield - mn) + mn, 0.95*e1$PC1,
         col= "brown", lwd=1.8,length=0.1)
  
  # GGB example
  library(agridat)
  data(crossa.wheat)
  dat2 <- crossa.wheat
  libs(gge)
  # Specify env.group as column in data frame
  m2 <- gge(dat2, yield~gen*loc,
            env.group=locgroup, gen.group=gengroup,
            scale=FALSE)
  biplot(m2, main="crossa.wheat - GGB biplot")
  

## End(Not run)

Germination of Orobanche seeds for two genotypes and two treatments.

Description

Number of Orobanche seeds tested/germinated for two genotypes and two treatments.

Format

plate

Factor for replication

gen

Factor for genotype with levels O73, O75

extract

Factor for extract from bean, cucumber

germ

Number of seeds that germinated

n

Total number of seeds tested

Details

Egyptian broomrape, orobanche aegyptiaca is a parasitic plant family. The plants have no chlorophyll and grow on the roots of other plants. The seeds remain dormant in soil until certain compounds from living plants stimulate germination.

Two genotypes were studied in the experiment, O. aegyptiaca 73 and O. aegyptiaca 75. The seeds were brushed with one of two extracts prepared from either a bean plant or cucmber plant.

The experimental design was a 2x2 factorial, each with 5 or 6 reps of plates.

Source

Crowder, M.J., 1978. Beta-binomial anova for proportions. Appl. Statist., 27, 34-37. https://doi.org/10.2307/2346223

References

N. E. Breslow and D. G. Clayton. 1993. Approximate inference in generalized linear mixed models. Journal of the American Statistical Association, 88:9-25. https://doi.org/10.2307/2290687

Y. Lee and J. A. Nelder. 1996. Hierarchical generalized linear models with discussion. J. R. Statist. Soc. B, 58:619-678.

Examples

## Not run: 

  library(agridat)
  data(crowder.seeds)
  dat <- crowder.seeds
  m1.glm <- m1.glmm <- m1.glmmtmb <- m1.hglm <- NA


  # ----- Graphic
  libs(lattice)
  dotplot(germ/n~gen|extract, dat, main="crowder.seeds")


  # --- GLMM.  Assumes Gaussian random effects
  libs(MASS)
  m1.glmm <- glmmPQL(cbind(germ, n-germ) ~ gen*extract, random= ~1|plate,
                     family=binomial(), data=dat)
  summary(m1.glmm)
  ## round(summary(m1.glmm)$tTable,2)
  ##                        Value Std.Error DF t-value p-value
  ## (Intercept)            -0.44      0.25 17   -1.80    0.09
  ## genO75                 -0.10      0.31 17   -0.34    0.74
  ## extractcucumber         0.52      0.34 17    1.56    0.14
  ## genO75:extractcucumber  0.80      0.42 17    1.88    0.08


  # ----- glmmTMB
  libs(glmmTMB)
  m1.glmmtmb <- glmmTMB(cbind(germ, n-germ) ~ gen*extract + (1|plate),
                        data=dat,
                        family=binomial)
  summary(m1.glmmtmb)
  ## round(summary(m1.glmmtmb)$coefficients$cond , 2)
  ##                        Estimate Std. Error z value Pr(>|z|)
  ## (Intercept)               -0.45       0.22   -2.03     0.04
  ## genO75                    -0.10       0.28   -0.35     0.73
  ## extractcucumber            0.53       0.30    1.74     0.08
  ## genO75:extractcucumber     0.81       0.38    2.11     0.04
 

## End(Not run)

Early generation variety trial in wheat

Description

Early generation variety trial in wheat

Format

A data frame with 670 observations on the following 5 variables.

gen

genotype factor

row

row

col

column

entry

entry (genotype) number

yield

yield of each plot, kg/ha

weed

weed score

Details

The data are from a field experiment conducted at Tullibigeal, New South Wales, Australia in 1987-88. The aim of these trials is to identify and retain the top (10-20 percent) lines for further testing.

Most genotypes are unreplicated, with some augmented genotypes. In each row, every 6th plot was variety 526 = 'Kite'. Six other varieties 527-532 were randomly placed in the trial, with 3 to 5 plots of each. Each plot was 15m x 1.8m, "oriented with the longest side with rows".

The 'weed' variable is a visual score on a 0 to 10 scale, 0 = no weeds, 10 = 100 percent weeds.

Cullis et al. (1989) presented an analysis of early generation variety trials that included a one-dimensional spatial analysis. Below, a two-dimensional spatial analysis is presented.

Note: The 'row' and 'col' variables are as in the VSN link below (switched compared to the paper by Cullis et al.)

Field width: 10 rows * 15 m = 150 m

Field length: 67 plots * 1.8 m = 121 m

The orientation is not certain, but the alternative orientation would have a field roughly 20m x 1000m, which seems unlikely.

Source

Brian R. Cullis, Warwick J. Lill, John A. Fisher, Barbara J. Read and Alan C. Gleeson (1989). A New Procedure for the Analysis of Early Generation Variety Trials. Journal of the Royal Statistical Society. Series C (Applied Statistics), 38, 361-375. https://doi.org/10.2307/2348066

References

Unreplicated early generation variety trial in Wheat. https://www.vsni.co.uk/software/asreml/htmlhelp/asreml/xwheat.htm

Examples

## Not run: 

  library(agridat)
  data(cullis.earlygen)
  dat <- cullis.earlygen

  # Show field layout of checks.  Cullis Table 1.
  dat$check <- ifelse(dat$entry < 8, dat$entry, NA)
  libs(desplot)
  desplot(dat, check ~ col*row,
          num=entry, cex=0.5, flip=TRUE, aspect=121/150, # true aspect
          main="cullis.earlygen (yield)")

  desplot(dat, yield ~ col*row,
          num="check", cex=0.5, flip=TRUE, aspect=121/150, # true aspect
          main="cullis.earlygen (yield)")

  grays <- colorRampPalette(c("white","#252525"))
  desplot(dat, weed ~ col*row,
          at=0:6-0.5, col.regions=grays(7)[-1],
          flip=TRUE, aspect=121/150, # true aspect
          main="cullis.earlygen (weed)")

  libs(lattice)
  bwplot(yield ~ as.character(weed), dat,
         horizontal=FALSE,
         xlab="Weed score", main="cullis.earlygen")

  # Moving Grid
  libs(mvngGrAd)
  shape <- list(c(1),
                c(1),
                c(1:4),
                c(1:4))
  # sketchGrid(10,10,20,20,shapeCross=shape, layers=1, excludeCenter=TRUE)
  m0 <- movingGrid(rows=dat$row, columns=dat$col, obs=dat$yield,
                   shapeCross=shape, layers=NULL)
  dat$mov.avg <- fitted(m0)

  if(require("asreml", quietly=TRUE)) {
    libs(asreml,lucid)

    # Start with the standard AR1xAR1 analysis
    dat <- transform(dat, xf=factor(col), yf=factor(row))
    dat <- dat[order(dat$xf, dat$yf),]
    m2 <- asreml(yield ~ weed, data=dat, random= ~gen,
                 resid = ~ ar1(xf):ar1(yf))
    
    # Variogram suggests a polynomial trend
    m3 <- update(m2, fixed= yield~weed+pol(col,-1))
    
    # Now add a nugget variance
    m4 <- update(m3, random= ~ gen + units)
    
    lucid::vc(m4)
    ##       effect component std.error z.ratio bound 
    ##          gen  73780    10420         7.1     P 0  
    ##        units  30440     8073         3.8     P 0.1
    ##     xf:yf(R)  54730    10630         5.1     P 0  
    ## xf:yf!xf!cor      0.38     0.115     3.3     U 0  
    ## xf:yf!yf!cor      0.84     0.045    19       U 0  
    
    ## # Predictions from models m3 and m4 are non-estimable.  Why?
    ## # Use model m2 for predictions
    ## predict(m2, classify="gen")$pvals
    ## ##         gen predicted.value std.error    status
    ## ## 1     Banks        2723.534  93.14719 Estimable
    ## ## 2    Eno008        2981.056 162.85241 Estimable
    ## ## 3    Eno009        2978.008 161.57129 Estimable
    ## ## 4    Eno010        2821.399 153.96943 Estimable
    ## ## 5    Eno011        2991.612 161.53507 Estimable
    
    
    ## # Compare AR1 with Moving Grid
    ## dat$ar1 <- fitted(m2)
    ## head(dat[ , c('yield','ar1','mov.avg')])
    ## ##    yield      ar1       mg
    ## ## 1   2652 2467.980 2531.998
    ## ## 11  3394 3071.681 3052.160
    ## ## 21  3148 2826.188 2807.031
    ## ## 31  3426 3026.985 3183.649
    ## ## 41  3555 3070.102 3195.910
    ## ## 51  3453 3006.352 3510.511
    ## pairs(dat[ , c('yield','ar1','mg')])
  }
  

## End(Not run)

Incomplete-block experiment of maize in Ethiopia.

Description

Incomplete-block experiment of maize in Ethiopia.

Usage

data("damesa.maize")

Format

A data frame with 264 observations on the following 8 variables.

site

site, 4 levels

rep

replicate, 3 levels

block

incomplete block

plot

plot number

gen

genotype, 22 levels

row

row ordinate

col

column ordinate

yield

yield, t/ha

Details

An experiment harvested in 2012, evaluating drought-tolerant maize hybrids at 4 sites in Ethiopia. At each site, an incomplete-block design was used.

Damesa et al use this data to compare single-stage and two-stage analyses.

Source

Tigist Mideksa Damesa, Jens Möhring, Mosisa Worku, Hans-Peter Piepho (2017). One Step at a Time: Stage-Wise Analysis of a Series of Experiments. Agronomy J, 109, 845-857. https://doi.org/10.2134/agronj2016.07.0395

References

None

Examples

## Not run: 
  library(agridat)
  data(damesa.maize)
  libs(desplot)
  desplot(damesa.maize,
          yield ~ col*row|site,
          main="damesa.maize",
          out1=rep, out2=block, num=gen, cex=1)

  if(require("asreml", quietly=TRUE)) {
    # Fit the single-stage model in Damesa
    libs(asreml,lucid)
    m0 <- asreml(data=damesa.maize,
                 fixed = yield ~ gen,
                 random = ~ site + gen:site + at(site):rep/block,
                 residual = ~ dsum( ~ units|site) )
    lucid::vc(m0) # match Damesa table 1 column 3
    ##                 effect component std.error z.ratio bound 
    ##       at(site, S1):rep   0.08819   0.1814     0.49     P 0  
    ##       at(site, S2):rep   1.383     1.426      0.97     P 0  
    ##       at(site, S3):rep   0              NA      NA     B 0  
    ##       at(site, S4):rep   0.01442   0.02602    0.55     P 0  
    ##                   site  10.45      8.604      1.2      P 0.1
    ##               gen:site   0.1054    0.05905    1.8      P 0.1
    ## at(site, S1):rep:block   0.3312    0.3341     0.99     P 0  
    ## at(site, S2):rep:block   0.4747    0.1633     2.9      P 0  
    ## at(site, S3):rep:block   0              NA      NA     B 0  
    ## at(site, S4):rep:block   0.06954   0.04264    1.6      P 0  
    ##              site_S1!R   1.346     0.3768     3.6      P 0  
    ##              site_S2!R   0.1936    0.06628    2.9      P 0  
    ##              site_S3!R   1.153     0.2349     4.9      P 0  
    ##              site_S4!R   0.1112    0.03665    3        P 0  
  }
  

## End(Not run)

Darwin's maize data of crossed/inbred plant heights

Description

Darwin's maize data of crossed/inbred plant heights.

Format

A data frame with 30 observations on the following 4 variables.

pot

Pot factor, 4 levels

pair

Pair factor, 12 levels

type

Type factor, self-pollinated, cross-pollinated

height

Height, in inches (measured to 1/8 inch)

Details

Charles Darwin, in 1876, reported data from an experiment that he had conducted on the heights of corn plants. The seeds came from the same parents, but some seeds were produced from self-fertilized parents and some seeds were produced from cross-fertilized parents. Pairs of seeds were planted in pots. Darwin hypothesized that cross-fertilization produced produced more robust and vigorous offspring.

Darwin wrote, "I long doubted whether it was worth while to give the measurements of each separate plant, but have decided to do so, in order that it may be seen that the superiority of the crossed plants over the self-fertilised, does not commonly depend on the presence of two or three extra fine plants on the one side, or of a few very poor plants on the other side. Although several observers have insisted in general terms on the offspring from intercrossed varieties being superior to either parent-form, no precise measurements have been given;* and I have met with no observations on the effects of crossing and self-fertilising the individuals of the same variety. Moreover, experiments of this kind require so much time–mine having been continued during eleven years–that they are not likely soon to be repeated."

Darwin asked his cousin Francis Galton for help in understanding the data. Galton did not have modern statistical methods to approach the problem and said, "I doubt, after making many tests, whether it is possible to derive useful conclusions from these few observations. We ought to have at least 50 plants in each case, in order to be in a position to deduce fair results".

Later, R. A. Fisher used Darwin's data in a book about design of experiments and showed that a t-test exhibits a significant difference between the two groups.

Source

Darwin, C. R. 1876. The effects of cross and self fertilisation in the vegetable kingdom. London: John Murray. Page 16. https://darwin-online.org.uk/converted/published/1881_Worms_F1357/1876_CrossandSelfFertilisation_F1249/1876_CrossandSelfFertilisation_F1249.html

References

R. A. Fisher, (1935) The Design of Experiments, Oliver and Boyd. Page 30.

Examples

## Not run: 

library(agridat)

data(darwin.maize)
dat <- darwin.maize

# Compare self-pollination with cross-pollination
libs(lattice)
bwplot(height~type, dat, main="darwin.maize")

  libs(reshape2)
  dm <- melt(dat)
  d2 <- dcast(dm, pot+pair~type)
  d2$diff <- d2$cross-d2$self
  t.test(d2$diff)
  ## 	One Sample t-test
  ## t = 2.148, df = 14, p-value = 0.0497
  ## alternative hypothesis: true mean is not equal to 0
  ## 95 percent confidence interval:
  ##  0.003899165 5.229434169


## End(Not run)

Multi-environment trial of maize

Description

Multi-environment trial of maize with 3 reps.

Usage

data("dasilva.maize")

Format

A data frame with 1485 observations on the following 4 variables.

env

environment

rep

replicate block, 3 per env

gen

genotype

yield

yield (tons/hectare)

Details

Each location had 3 blocks. Block numbers are unique across environments.

NOTE! The environment codes in the supplemental data file of da Silva 2015 do not quite match the environment codes of the paper, but are mostly off by 1.

DaSilva Table 1 has a footnote "Machado et al 2007". This reference appears to be:

Machado et al. Estabilidade de producao de hibridos simples e duplos de milhooriundos de um mesmo conjunto genico. Bragantia, 67, no 3. www.scielo.br/pdf/brag/v67n3/a10v67n3.pdf

In DaSilva Table 1, the mean of E1 is 10.803. This appears to be a copy of the mean from row 1 of Table 1 in Machado. Using the supplemental data from this paper, the correct mean is 8.685448.

Source

A Bayesian Shrinkage Approach for AMMI Models. Carlos Pereira da Silva, Luciano Antonio de Oliveira, Joel Jorge Nuvunga, Andrezza Kellen Alves Pamplona, Marcio Balestre. Plos One. Supplemental material. https://doi.org/10.1371/journal.pone.0131414

Used via license: Creative Commons BY-SA.

References

J.J. Nuvunga, L.A. Oliveira, A.K.A. Pamplona, C.P. Silva, R.R. Lima and M. Balestre. Factor analysis using mixed models of multi-environment trials with different levels of unbalancing. Genet. Mol. Res. 14.

Examples

library(agridat)
data(dasilva.maize)
dat <- dasilva.maize

# Try to match Table 1 of da Silva 2015.
# aggregate(yield ~ env, data=dat, FUN=mean)
##   env     yield
## 1  E1  6.211817  # match E2 in Table 1
## 2  E2  4.549104  # E3
## 3  E3  5.152254  # E4
## 4  E4  6.245904  # E5
## 5  E5  8.084609  # E6
## 6  E6 13.191890  # E7
## 7  E7  8.895721  # E8
## 8  E8  8.685448  
## 9  E9  8.737089  # E9

# Unable to match CVs in Table 2, but who knows what they used
# for residual variance.
# aggregate(yield ~ env, data=dat, FUN=function(x) 100*sd(x)/mean(x))

# Match DaSilva supplement 2, ANOVA
# m1 <- aov(yield ~ env + gen + rep:env + gen:env, dat)
# anova(m1)
## Response: yield
##            Df Sum Sq Mean Sq  F value    Pr(>F)    
## env         8 8994.2 1124.28 964.1083 < 2.2e-16 ***
## gen        54  593.5   10.99   9.4247 < 2.2e-16 ***
## env:rep    18   57.5    3.19   2.7390 0.0001274 ***
## env:gen   432  938.1    2.17   1.8622 1.825e-15 ***
## Residuals 972 1133.5    1.17                       

Uniformity trial of soybean

Description

Uniformity trial of soybean in Brazil, 1970.

Usage

data("dasilva.soybean.uniformity")

Format

A data frame with 1152 observations on the following 3 variables.

row

row

col

column

yield

yield, grams/plot

Details

Field length: 48 rows * .6 m = 28.8 m

Field width: 24 columns * .6 m = 14.4 m

Source

Enedino Correa da Silva. (1974). Estudo do tamanho e forma de parcelas para experimentos de soja (Plot size and shape for soybean yield trials). Pesquisa Agropecuaria Brasileira, Serie Agronomia, 9, 49-59. Table 3, page 52-53. https://seer.sct.embrapa.br/index.php/pab/article/view/17250

References

Humada-Gonzalez, G.G. (2013). Estimação do tamanho otimo de parcela experimental em experimento com soja. Dissertation, Universidade Federal de Lavras. http://repositorio.ufla.br/jspui/handle/1/744

Examples

## Not run: 

library(agridat)
data(dasilva.soybean.uniformity)
dat <- dasilva.soybean.uniformity

libs(desplot)
desplot(dat, yield ~ col*row,
        flip=TRUE, aspect=28.8/14.4, 
        main="dasilva.soybean.uniformity")
  

## End(Not run)

Growth of soybean varieties in 3 years

Description

Growth of soybean varieties in 3 years

Usage

data("davidian.soybean")

Format

A data frame with 412 observations on the following 5 variables.

plot

plot code

variety

variety, F or P

year

1988-1990

day

days after planting

weight

weight of soybean leaves

Details

This experiment compared the growth patterns of two genotypes of soybean varieties: F=Forrest (commercial variety) and P=Plant Introduction number 416937 (experimental variety).

Data were collected in 3 consecutive years.

At the start of each growing season, 16 plots were seeded (8 for each variety). Data were collected approximately weekly. At each timepoint, six plants were randomly selected from each plot. The leaves from these 6 plants were weighed, and average leaf weight per plant was reported. (We assume that the data collection is destructive and different plants are sampled at each date).

Note: this data is the same as the "nlme::Soybean" data.

Source

Marie Davidian and D. M. Giltinan, (1995). Nonlinear Models for Repeated Measurement Data. Chapman and Hall, London.

Electronic version retrieved from https://www4.stat.ncsu.edu/~davidian/data/soybean.dat

References

Pinheiro, J. C. and Bates, D. M. (2000). Mixed-Effects Models in S and S-PLUS. Springer, New York.

Examples

## Not run: 

  library(agridat)
  data(davidian.soybean)
  dat <- davidian.soybean
  dat$year <- factor(dat$year)

  libs(lattice)
  xyplot(weight ~ day|variety*year, dat,
         group=plot, type='l',
         main="davidian.soybean")


  # The only way to keep your sanity with nlme is to use groupedData objects
  # Well, maybe not.  When I use "devtools::run_examples",
  # the "groupedData" function creates a dataframe with/within(?) an
  # environment, and then "nlsList" cannot find datg, even though
  # ls() shows datg is visible and head(datg) is fine.
  # Also works fine in interactive mode. It is driving me insane.
  # reid.grasses has the same problem
  # Use if(0){} to block this code from running.
  if(0){
    libs(nlme)
    datg <- groupedData(weight ~ day|plot, dat)
    # separate fixed-effect model for each plot
    # 1988P6 gives unusual estimates
    m1 <- nlsList(SSlogis, data=datg,
                  subset = plot != "1988P6")
    # plot(m1) # seems heterogeneous
    plot(intervals(m1), layout=c(3,1)) # clear year,variety effects in Asym

    # A = maximum, B = time of half A = steepness of curve
    # C = sharpness of curve (smaller = sharper curve)

    # switch to mixed effects
    m2 <- nlme(weight ~ A / (1+exp(-(day-B)/C)),
               data=datg,
               fixed=list(A ~ 1, B ~ 1, C ~ 1),
               random = A +B +C ~ 1,
               start=list(fixed = c(17,52,7.5))) # no list!

    # add covariates for A,B,C effects, correlation, weights
    # not necessarily best model, but it shows the syntax
    m3 <- nlme(weight ~ A / (1+exp(-(day-B)/C)),
               data=datg,
               fixed=list(A ~ variety + year,
                          B ~ year,
                        C ~ year),
               random = A +B +C ~ 1,
               start=list(fixed= c(19,0,0,0,
                                   55,0,0,
                                   8,0,0)),
               correlation = corAR1(form = ~ 1|plot),
               weights=varPower(), # really helps
               control=list(mxMaxIter=200))

    plot(augPred(m3), layout=c(8,6),
      main="davidian.soybean - model 3")
  } # end if(0)

## End(Not run)

Uniformity trial of pasture.

Description

Uniformity trial of pasture in Australia.

Usage

data("davies.pasture.uniformity")

Format

A data frame with 760 observations on the following 3 variables.

row

row

col

column

yield

yield per plot, grams

Details

Conducted at the Waite Agricultural Research Institute in 1928. A rectangle 250 x 200 links was selected, divided into 1000 plots measuring 10 x 5 links, that is 1/2000th acre. Plots were hand harvested for herbage and air-dried. Cutting began Tue, 25 Sep and ended Sat, 29 Sep, by which time 760 plots had been harvested. Rain fell, harvesting ceased.

The minimum recommended plot size is 150 square links. The optimum recommended plot size is 450 square links, 5 x 90 links in size.

Note, there were 4 digits that were hard to read in the original document. Best estimates of these digits were used for the yields of the affects plots. The yields were digitally watermarked with an extra .01 added to the yield value.

The botanical composition of species clearly influenced the total herbage.

Field length: 40 plots * 5 links = 200 links

Field width: 19 plots * 10 links = 190 links

Source

J. Griffiths Davies (1931). The Experimental Error of the Yield from Small Plots of Natural Pasture. Council for Scientific and Industrial Research (Aust.) Bulletin 48. Table 1.

References

None

Examples

## Not run: 

library(agridat)

  data(davies.pasture.uniformity)
  dat <- davies.pasture.uniformity
  
  # range(dat$yield) # match Davies
  # mean(dat$yield) # 227.77, Davies has 221.7
  # sd(dat$yield)/mean(dat$yield) # 33.9, Davies has 32.5

  # libs(lattice)
  # qqmath( ~ yield, dat) # clearly non-normal, skewed right

  libs(desplot)
  desplot(dat, yield ~ col*row,
          flip=TRUE, aspect=(40*5)/(19*10), # true aspect
          main="davies.pasture.uniformity") 


## End(Not run)

Uniformity trial of wheat

Description

Uniformity trial of wheat in 1903 in Missouri.

Usage

data("day.wheat.uniformity")

Format

A data frame with 3090 observations on the following 4 variables.

row

row

col

col

grain

grain weight, grams per plot

straw

straw weight, grams per plot

Details

These data are from the Shelbina field of the Missouri Agricultural Experiment Station. The field (plat) was about 1/4 acre in area and apparently uniform throughout. In the fall of 1912, wheat was drilled in rows 8 inches apart, each row 155 feet long. The wheat was harvested in June, in 5-foot segments. The gross weight and the grain weight was measured, the straw weight was calculated by subtraction.

Field width: 31 series * 5 feet = 155 feet

Field length: 100 rows, 8 inches apart = 66.66 feet

Source

James Westbay Day (1916). The relation of size, shape, and number of replications of plats to probable error in field experimentation. Dissertation, University of Missouri. Table 1, page 22. https://hdl.handle.net/10355/56391

References

James W. Day (1920). The relation of size, shape, and number of replications of plats to probable error in field experimentation. Agronomy Journal, 12, 100-105. https://doi.org/10.2134/agronj1920.00021962001200030002x

Examples

## Not run: 

library(agridat)
data(day.wheat.uniformity)
dat <- day.wheat.uniformity

libs(desplot)
desplot(dat, grain~col*row,
        flip=TRUE, aspect=(100*8)/(155*12), # true aspect
        main="day.wheat.uniformity - grain yield")
  
# similar to Day table IV
libs(lattice)
xyplot(grain~straw, data=dat, main="day.wheat.uniformity", type=c('p','r'))
# cor(dat$grain, dat$straw) # .9498 # Day calculated 0.9416
  
libs(desplot)
desplot(dat, straw~col*row,
        flip=TRUE, aspect=(100*8)/(155*12), # true aspect
        main="day.wheat.uniformity - straw yield")
  
# Day fig 2
coldat <- aggregate(grain~col, dat, sum) 
xyplot(grain ~ col, coldat, type='l', ylim=c(2500,6500))
dat$rowgroup <- round((dat$row +1)/3,0)
rowdat <- aggregate(grain~rowgroup, dat, sum) 
xyplot(grain ~ rowgroup, rowdat, type='l', ylim=c(2500,6500))


## End(Not run)  

Multi-environment trial with structured missing values

Description

Grain yield was measured on 5 genotypes in 26 environments. Missing values were non-random, but structured.

Format

env

environment, 26 levels

gen

genotype factor, 5 levels

yield

yield

Used with permission of Jean-Baptists Denis.

Source

Denis, J. B. and C P Baril, 1992, Sophisticated models with numerous missing values: The multiplicative interaction model as an example. Biul. Oceny Odmian, 24–25, 7–31.

References

H P Piepho, (1999) Stability analysis using the SAS system, Agron Journal, 91, 154–160. https://doi.og/10.2134/agronj1999.00021962009100010024x

Examples

## Not run: 

library(agridat)
data(denis.missing)
dat <- denis.missing

# view missingness structure
libs(reshape2)
acast(dat, env~gen, value.var='yield')


libs(lattice)
redblue <- colorRampPalette(c("firebrick", "lightgray", "#375997"))
levelplot(yield ~ gen*env, data=dat,
          col.regions=redblue,
          main="denis.missing - incidence heatmap")

# stability variance (Table 3 in Piepho)
libs(nlme)
m1 <- lme(yield ~ -1 + gen, data=dat, random= ~ 1|env,
          weights = varIdent(form= ~ 1|gen),
          na.action=na.omit)
svar <- m1$sigma^2 * c(1, coef(m1$modelStruct$varStruct, unc = FALSE))^2
round(svar, 2)
##          G5    G3    G1    G2
## 39.25 22.95 54.36 12.17 23.77


## End(Not run)

Multi-environment trial of perennial ryegrass in France

Description

Plant strength of perennial ryegrass in France for 21 genotypes at 7 locations.

Format

A data frame with 147 observations on the following 3 variables.

gen

genotype, 21 levels

loc

location, 7 levels

strength

average plant strength * 100

Details

INRA conducted a breeding trial in western France with 21 genotypes at 7 locations. The observed data is 'strength' averaged over 7-10 plants per plot and three plots per location (after adjusting for blocking effects). Each plant was scored on a scale 0-9.

The original data had a value of 86.0 for genotype G1 at location L4–this was replaced by an additive estimated value of 361.2 as in Gower and Hand (1996).

Source

Jean-Baptiste Denis and John C. Gower, 1996. Asymptotic confidence regions for biadditive models: interpreting genotype-environment interaction, Applied Statistics, 45, 479-493. https://doi.org/10.2307/2986069

References

Gower, J.C. and Hand, D.J., 1996. Biplots. Chapman and Hall.

Examples

library(agridat)
data(denis.ryegrass)
dat <- denis.ryegrass

# biplots (without ellipses) similar to Denis figure 1
libs(gge)
m1 <- gge(dat, strength ~ gen*loc, scale=FALSE)
biplot(m1, main="denis.ryegrass biplot")

Latin square of four breeds of sheep with four diets

Description

Latin square of four breeds of sheep with four diets

Usage

data("depalluel.sheep")

Format

A data frame with 32 observations on the following 5 variables.

food

diet

animal

animal number

breed

sheep breed

weight

weight, pounds

date

months after start

Details

This may be the earliest known Latin Square experiment.

Four sheep from each of four breeds were randomized to four feeds and four slaughter dates.

Sheep that eat roots will eat more than sheep eating corn, but each acre of land produces more roots than corn.

de Palleuel said: In short, by adopting the use of roots, instead of corn, for the fattening of all sorts of cattle, the farmers in the neighborhood of the capital will not only gain great profit themselves, but will also very much benefit the public by supplying this great city with resources, and preventing the sudden rise of meat in her markets, which is often considerable.

Source

M. Crette de Palluel (1788). On the advantage and economy of feeding sheep in the house with roots. Annals of Agriculture, 14, 133-139. https://books.google.com/books?id=LXIqAAAAYAAJ&pg=PA133

References

None

Examples

## Not run: 

library(agridat)

data(depalluel.sheep)
dat <- depalluel.sheep

# Not the best view...weight gain is large in the first month, then slows down
# and the linear line hides this fact
libs(lattice)
xyplot(weight ~ date|food, dat, group=animal, type='l', auto.key=list(columns=4),
       xlab="Months since start",
       main="depalluel.sheep")

## End(Not run)

Graeco-Latin Square experiment in pine

Description

Graeco-Latin Square experiment in pine

Usage

data("devries.pine")

Format

A data frame with 36 observations on the following 6 variables.

block

block

row

row

col

column

spacing

spacing treatment

thinning

thinning treatment

volume

stem volume in m^3/ha

growth

annual stem volume increment m^3/ha at age 11

Details

Experiment conducted on Caribbean Pine at Coebiti in Surinam (Long 55 28 30 W, Lat 5 18 5 N). Land was cleared in Jan 1965 and planted May 1965. Each experimental plot was 60m x 60m. Roads 10 m wide run between the rows. Each block is thus 180m wide and 200m deep. Data were collected only on 40m x 40m plots in the center of each experimental unit. Plots were thinned in 1972 and 1975. The two treatment factors (spacing, thinning) were assigned in a Graeco-Latin Square design.

Spacing: A=2.5, B=3, C=3.5. Thinning: Z=low, M=medium, S=heavy.

Field width: 4 blocks x 180 m = 720 m

Field length: 1 block x 200 m = 200 m.

Source

P.G. De Vries, J.W. Hildebrand, N.R. De Graaf. (1978). Analysis of 11 years growth of carribbean pine in a replicated Graeco-Latin square spacing-thinning experiment in Surinam. Page 46, 51. https://edepot.wur.nl/287590

References

None

Examples

## Not run: 

  library(agridat)
  data(devries.pine)
  dat <- devries.pine

  libs(desplot)
  desplot(dat, volume ~ col*row,
          main="devries.pine - expt design and tree volume",
          col=spacing, num=thinning, cex=1, out1=block, aspect=200/720)
  
  libs(HH)
  HH::interaction2wt(volume ~ spacing+thinning, dat,
                     main="devries.pine")

  # ANOVA matches appendix 5 of DeVries
  m1 <- aov(volume ~ block + spacing + thinning + block:factor(row) +
              block:factor(col), data=dat)
  anova(m1)


## End(Not run)

Multi-environment trial of wheat

Description

Yield of 10 spring wheat varieties for 17 locations in 1976.

Format

A data frame with 134 observations on the following 3 variables.

gen

genotype, 10 levels

env

environment, 17 levels

yield

yield (t/ha)

Details

Yield of 10 spring wheat varieties for 17 locations in 1976.

Used to illustrate modified joint regression.

Source

Digby, P.G.N. (1979). Modified joint regression analysis for incomplete variety x environment data. Journal of Agricultural Science, 93, 81-86. https://doi.org/10.1017/S0021859600086159

References

Hans-Pieter Piepho, 1997. Analyzing Genotype-Environment Data by Mixed-Models with Multiplicative Terms. Biometrics, 53, 761-766. https://doi.org/10.2307/2533976

RJOINT procedure in GenStat. https://www.vsni.co.uk/software/genstat/htmlhelp/server/RJOINT.htm

Examples

## Not run: 
  
  library(agridat)
  data(digby.jointregression)
  dat <- digby.jointregression
  
  # Simple gen means, ignoring unbalanced data.
  # Matches Digby table 2, Unadjusted Mean
  round(tapply(dat$yield, dat$gen, mean),3)
  
  # Two-way model. Matches Digby table 2, Fitting Constants
  m00 <- lm(yield ~ 0 + gen + env, dat)
  round(coef(m00)[1:10]-2.756078+3.272,3) # Adjust intercept
  # genG01 genG02 genG03 genG04 genG05 genG06 genG07 genG08 genG09 genG10 
  #  3.272  3.268  4.051  3.724  3.641  3.195  3.232  3.268  3.749  3.179 
  
  n.gen <- nlevels(dat$gen)
  n.env <- nlevels(dat$env)
  
  # Estimate theta (env eff)
  m0 <- lm(yield ~ -1 + env + gen, dat)
  thetas <- coef(m0)[1:n.env]
  thetas <- thetas-mean(thetas) # center env effects
  # Add env effects to the data
  dat$theta <- thetas[match(paste("env",dat$env,sep=""), names(thetas))]
  
  # Initialize beta (gen slopes) at 1
  betas <- rep(1, n.gen)
  
  done <- FALSE
  while(!done){
    
    betas0 <- betas
    
    # M1: Fix thetas (env effects), estimate beta (gen slope)
    m1 <- lm(yield ~ -1 + gen + gen:theta, data=dat)
    betas <- coef(m1)[-c(1:n.gen)]
    dat$beta <- betas[match(paste("gen",dat$gen,":theta",sep=""), names(betas))]
    # print(betas)

    # M2: Fix betas (gen slopes), estimate theta (env slope)
    m2 <- lm(yield ~ env:beta + gen -1, data=dat)
    thetas <- coef(m2)[-c(1:n.gen)]
    thetas[is.na(thetas)] <- 0  # Change last coefficient from NA to 0
    dat$theta <- thetas[match(paste("env",dat$env,":beta",sep=""), names(thetas))]
    # print(thetas)

    # Check convergence
    chg <- sum(((betas-betas0)/betas0)^2)
    cat("Relative change in betas",chg,"\n")
    if(chg < .0001) done <- TRUE
    
  }

  libs(lattice)
  xyplot(yield ~ theta|gen, data=dat, xlab="theta (environment effect)",
         main="digby.jointregression - stability plot")

  # Dibgy Table 2, modified joint regression
  
  # Genotype sensitivities (slopes)
  round(betas,3) # Match Digby table 2, Modified joint regression sensitivity
  # genG01 genG02 genG03 genG04 genG05 genG06 genG07 genG08 genG09 genG10
  #  0.953  0.739  1.082  1.024  1.142  0.877 1.089  0.914  1.196  0.947

  # Env effects. Match Digby table 3, Modified joint reg
  round(thetas,3)+1.164-.515 # Adjust intercept to match
  # envE01 envE02 envE03 envE04 envE05 envE06 envE07 envE08 envE09 envE10
  # -0.515 -0.578 -0.990 -1.186  1.811  1.696 -1.096  0.046  0.057  0.825
  # envE11 envE12 envE13 envE14 envE15 envE16 envE17
  # -0.576  1.568 -0.779 -0.692  0.836 -1.080  0.649

  # Using 'gnm' gives similar results.
  # libs(gnm)
  # m3 <- gnm(yield ~ gen + Mult(gen,env), data=dat) # slopes negated
  # round(coef(m3)[11:20],3)

  # Using 'mumm' gives similar results, though gen is random and the
  # coeffecients are shrunk toward 0 a bit.
  if(require("mumm", quietly=TRUE)) {
    libs(mumm)
    m1 <- mumm(yield ~ -1 + env + mp(gen, env), dat)
    round(1 + ranef(m1)$`mp gen:env`,2)
  }
  

## End(Not run)

Bodyweight of cows in a 2-by-2 factorial experiment

Description

Bodyweight of cows in a 2-by-2 factorial experiment.

Format

A data frame with 598 observations on the following 5 variables.

animal

Animal factor, 26 levels

iron

Factor with levels Iron, NoIron

infect

Factor levels Infected, NonInfected

weight

Weight in (rounded to nearest 5) kilograms

day

Days after birth

Details

Diggle et al., 1994, pp. 100-101, consider an experiment that studied how iron dosing (none/standard) and micro-organism (infected or non-infected) influence the weight of cows.

Twenty-eight cows were allocated in a 2-by-2 factorial design with these factors. Some calves were inoculated with tuberculosis at six weeks of age. At six months, some calves were maintained on supplemental iron diet for a further 27 months.

The weight of each animal was measured at 23 times, unequally spaced. One cow died during the study and data for another cow was removed.

Source

Diggle, P. J., Liang, K.-Y., & Zeger, S. L. (1994). Analysis of Longitudinal Data. Page 100-101.

Retrieved Oct 2011 from https://www.maths.lancs.ac.uk/~diggle/lda/Datasets/

References

Lepper, AWD and Lewis, VM, 1989. Effects of altered dietary iron intake in Mycobacterium paratuberculosis-infected dairy cattle: sequential observations on growth, iron and copper metabolism and development of paratuberculosis. Research in veterinary science, 46, 289–296.

Arunas P. Verbyla and Brian R. Cullis and Michael G. Kenward and Sue J. Welham, (1999), The analysis of designed experiments and longitudinal data by using smoothing splines. Appl. Statist., 48, 269–311.

SAS/STAT(R) 9.2 User's Guide, Second Edition. https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_glimmix_sect018.htm

Examples

## Not run: 

  library(agridat)
  data(diggle.cow)
  dat <- diggle.cow
  
  # Figure 1 of Verbyla 1999
  libs(latticeExtra)
  useOuterStrips(xyplot(weight ~ day|iron*infect, dat, group=animal,
                        type='b', cex=.5, 
                        main="diggle.cow"))
  
  # Scaling
  dat <- transform(dat, time = (day-122)/10)


  if(require("asreml", quietly=TRUE)) {
    libs(asreml, latticeExtra)
    
    ## # Smooth for each animal.  No treatment effects. Similar to SAS Output 38.6.9
    
    m1 <- asreml(weight ~ 1 + lin(time) + animal + animal:lin(time), data=dat,
                 random = ~ animal:spl(time))
    p1 <- predict(m1, data=dat, classify="animal:time",
                  design.points=list(time=seq(0,65.9, length=50)))
    p1 <- p1$pvals
    p1 <- merge(dat, p1, all=TRUE) # to get iron/infect merged in
    foo1 <- xyplot(weight ~ day|iron*infect, dat, group=animal,
                   main="diggle.cow")
    foo2 <- xyplot(predicted.value ~ day|iron*infect, p1, type='l', group=animal)
    print(foo1+foo2)
  }
  

## End(Not run)

Uniformity trial of safflower

Description

Uniformity trial of safflower in Arizona in 1958.

Usage

data("draper.safflower.uniformity")

Format

A data frame with 640 observations on the following 4 variables.

expt

experiment

row

row

col

column

yield

yield per plot (grams)

Details

Experiments were conducted at the Agricultural Experiment Station Farm at Eloy, Arizona. The crop was harvested in July 1958.

The crop was planted in two rows 12 inches apart on vegetable beds 40 inches center to center.

In each test, the end ranges and one row of plots on one side were next to alleys, and those plots gave estimates of border effects.

Experiment E4 (four foot test)

Sandy streaks were present in the field. Average yield was 1487 lb/ac. A diagonal fertility gradient was in this field. Widening the plot was equally effective as lengthening the plot to reduce variability. The optimum plot size was 1 bed wide, 24 feet long. Considering economic costs, the optimum size was 1 bed, 12 feet long.

Field width: 16 beds * 3.33 feet = 53 feet

Field length: 18 ranges * 4 feet = 72 feet

Experiment E5 (five foot test)

Average yield 2517 lb/ac, typical for this crop. Combining plots lengthwise was more effective than widening the plots, in order to reduce variability. The optimum plot size was 1 bed wide, 25 feet long. Considering economic costs, the optimum size was 1 bed, 18 feet long.

Field width: 14 beds * 3.33 feet = 46.6 feet.

Field length: 18 ranges * 5 feet = 90 feet.

Data are from Table A & B of Draper, p. 53-56. Typed by K.Wright.

Source

Arlen D. Draper. (1959). Optimum plot size and shape for safflower yield tests. Dissertation. University of Arizona. https://hdl.handle.net/10150/319371 Page 53-56.

References

None

Examples

## Not run: 

  library(agridat)
  data(draper.safflower.uniformity)
  dat4 <- subset(draper.safflower.uniformity, expt=="E4")
  dat5 <- subset(draper.safflower.uniformity, expt=="E5")
  
  libs(desplot)
  desplot(dat4, yield~col*row,
          flip=TRUE, tick=TRUE, aspect=72/53, # true aspect
          main="draper.safflower.uniformity (four foot)")
  
  desplot(dat5, yield~col*row,
          flip=TRUE, tick=TRUE, aspect=90/46, # true aspect
          main="draper.safflower.uniformity (five foot)")

  # Draper appears to removed the border plots, but it is difficult to
  # match his results exactly
  dat4 <- subset(dat4, row>1 & row<20)
  dat4 <- subset(dat4, col>1 & col<17)
  dat5 <- subset(dat5, row>1 & row<20)
  dat5 <- subset(dat5, col<15)
  # Convert gm/plot to pounds/acre. Draper (p. 20) says 1487 pounds/acre
  mean(dat4$yield) / 453.592 / (3.33*4) * 43560 # 1472 lb/ac
  
  libs(agricolae)
  libs(reshape2)
  
  s4 <- index.smith(acast(dat4, row~col, value.var='yield'),
                    main="draper.safflower.uniformity (four foot)",
                    col="red")$uni
  s4 # match Draper table 2, p 22
  
  ## s5 <- index.smith(acast(dat5, row~col, value.var='yield'),
  ##                   main="draper.safflower.uniformity (five foot)",
  ##                   col="red")$uni
  ## s5 # match Draper table 1, p 21
  

## End(Not run)

Uniformity trial of groundnut

Description

Uniformity trial of groundnut.

Usage

data("ducker.groundnut.uniformity")

Format

A data frame with 215 observations on the following 3 variables.

row

row ordinate

col

column ordinate

yield

yield, pounds per plot

Details

The experiment was grown in Nyasaland, Cotton Experiment Station, Domira Bay, 1942-43. There were 44x5 identical plots, each 1/220 acre in area. Single ridge plots each one chain in length, and one yard apart. Two rows of groundnuts are planted per ridge, staggered at 1 foot between holes. Holes are spaced 18 inches x 12 inches. Two seeds are planted per hole.

The yield values are pounds of nuts in shell.

Field length: 5 plots, 22 yards each = 110 yards.

Field width: 44 plots, 1 yard each = 44 yards.

This data was made available with special help from the staff at Rothamsted Research Library.

Data typed by K.Wright and checked by hand.

Source

Rothamsted Research Library, Box STATS17 WG Cochran, Folder 2.

References

None

Examples

## Not run: 
  library(agridat)
  data(ducker.groundnut.uniformity)
  dat <- ducker.groundnut.uniformity
  libs(desplot)
  desplot(dat, yield ~ col*row,
          flip=TRUE, aspect=110/44,
          main="ducker.groundnut.uniformity")

## End(Not run)

Sugar beet yields with competition effects

Description

Sugar beet yields with competition effects

Format

A data frame with 114 observations on the following 5 variables.

gen

Genotype factor, 36 levels plus Border

col

Column

block

Row/Block

wheel

Position relative to wheel tracks

yield

Root yields, kg/plot

Details

This sugar-beet trial was conducted in 1979.

Single-row plots, 12 m long, 0.5 m between rows. Each block is made up of all 36 genotypes laid out side by side. Guard/border plots are at each end. Root yields were collected.

Wheel tracks are located between columns 1 and 2, and between columns 5 and 6, for each set of six plots. Each genotype was randomly allocated once to each pair of plots (1,6), (2,5), (3,4) across the three reps. Wheel effect were not significant in _this_ trial.

Field width: 18m + 1m guard rows = 19m

Field length: 3 blocks * 12m + 2*0.5m spacing = 37m Retrieved from https://www.ma.hw.ac.uk/~iain/research/JAgSciData/data/Trial1.dat

Used with permission of Iain Currie.

Source

Durban, M., Currie, I. and R. Kempton, 2001. Adjusting for fertility and competition in variety trials. J. of Agricultural Science, 136, 129–140.

Examples

## Not run: 

library(agridat)

data(durban.competition)
dat <- durban.competition

# Check that genotypes were balanced across wheel tracks.
with(dat, table(gen,wheel))

libs(desplot)
desplot(dat, yield ~ col*block,
        out1=block, text=gen, col=wheel, aspect=37/19, # true aspect
        main="durban.competition")


# Calculate residual after removing block/genotype effects
m1 <- lm(yield ~ gen + block, data=dat)
dat$res <- resid(m1)

## desplot(dat, res ~ col*block, out1=block, text=gen, col=wheel,
##         main="durban.competition - residuals")

# Calculate mean of neighboring plots
dat$comp <- NA
dat$comp[3:36] <- ( dat$yield[2:35] + dat$yield[4:37] ) / 2
dat$comp[41:74] <- ( dat$yield[40:73] + dat$yield[42:75] ) / 2
dat$comp[79:112] <- ( dat$yield[78:111] + dat$yield[80:113] ) / 2

# Demonstrate the competition effect
# Competitor plots have low/high yield -> residuals are negative/positive
libs(lattice)
xyplot(res~comp, dat, type=c('p','r'), main="durban.competition",
       xlab="Average yield of neighboring plots", ylab="Residual")


## End(Not run)

Row-column experiment of spring barley, many varieties

Description

Row-column experiment of spring barley, many varieties

Format

A data frame with 544 observations on the following 5 variables.

row

row

bed

bed (column)

rep

rep, 2 levels

gen

genotype, 272 levels

yield

yield, tonnes/ha

Details

Spring barley variety trial of 272 entries (260 new varieties, 12 control). Grown at the Scottish Crop Research Institute in 1998. Row-column design with 2 reps, 16 rows (north/south) by 34 beds (east/west). The land sloped downward from row 16 to row 1. Plot yields were converted to tonnes per hectare.

Plot dimensions are not given.

Used with permission of Maria Durban.

Source

Durban, Maria and Hackett, Christine and McNicol, James and Newton, Adrian and Thomas, William and Currie, Iain. 2003. The practical use of semiparametric models in field trials, Journal of Agric Biological and Envir Stats, 8, 48-66. https://doi.org/10.1198/1085711031265

References

Edmondson, Rodney (2020). Multi-level Block Designs for Comparative Experiments. J of Agric, Biol, and Env Stats. https://doi.org/10.1007/s13253-020-00416-0

Examples

## Not run: 

  library(agridat)
  data(durban.rowcol)
  dat <- durban.rowcol
  
  libs(desplot)
  desplot(dat, yield~bed*row,
          out1=rep, num=gen, # aspect unknown
          main="durban.rowcol")
  

  # Durban 2003 Figure 1
  m10 <- lm(yield~gen, data=dat)
  dat$resid <- m10$resid
  ## libs(lattice)
  ## xyplot(resid~row, dat, type=c('p','smooth'), main="durban.rowcol")
  ## xyplot(resid~bed, dat, type=c('p','smooth'), main="durban.rowcol")
  
  # Figure 3
  libs(lattice)
  xyplot(resid ~ bed|factor(row), data=dat,
         main="durban.rowcol",
         type=c('p','smooth'))
  
  

  # Figure 5 - field trend
  # note, Durban used gam package like this
  # m1lo <- gam(yield ~ gen + lo(row, span=10/16) + lo(bed, span=9/34), data=dat)
  libs(mgcv)
  m1lo <- gam(yield ~ gen + s(row) + s(bed, k=5), data=dat)
  new1 <- expand.grid(row=unique(dat$row),bed=unique(dat$bed))
  new1 <- cbind(new1, gen="G001")
  p1lo <- predict(m1lo, newdata=new1)
  libs(lattice)
  wireframe(p1lo~row+bed, new1, aspect=c(1,.5), main="Field trend")


  if(require("asreml", quietly=TRUE)) {
    libs(asreml)

    dat <- transform(dat, rowf=factor(row), bedf=factor(bed))
    dat <- dat[order(dat$rowf, dat$bedf),]

    m1a1 <- asreml(yield~gen + lin(rowf) + lin(bedf), data=dat,
                   random=~spl(rowf) + spl(bedf) + units,
                   family=asr_gaussian(dispersion=1))
    m1a2 <- asreml(yield~gen + lin(rowf) + lin(bedf), data=dat,
                   random=~spl(rowf) + spl(bedf) + units,
                   resid = ~ar1(rowf):ar1(bedf))
    m1a2 <- update(m1a2)
    m1a3 <- asreml(yield~gen, data=dat, random=~units,
                   resid = ~ar1(rowf):ar1(bedf))
    
    # Figure 7
    libs(lattice)
    v7a <- asr_varioGram(x=dat$bedf, y=dat$rowf, z=m1a3$residuals)
    wireframe(gamma ~ x*y, v7a, aspect=c(1,.5)) # Fig 7a
    
    v7b <- asr_varioGram(x=dat$bedf, y=dat$rowf, z=m1a2$residuals)
    wireframe(gamma ~ x*y, v7b, aspect=c(1,.5)) # Fig 7b
    
    v7c <- asr_varioGram(x=dat$bedf, y=dat$rowf, z=m1lo$residuals)
    wireframe(gamma ~ x*y, v7c, aspect=c(1,.5)) # Fig 7c
  }
  

## End(Not run)

Split-plot experiment of barley with fungicide treatments

Description

Split-plot experiment of barley with fungicide treatments

Format

A data frame with 560 observations on the following 6 variables.

yield

yield, tonnes/ha

block

block, 4 levels

gen

genotype, 70 levels

fung

fungicide, 2 levels

row

row

bed

bed (column)

Details

Grown in 1995-1996 at the Scottish Crop Research Institute. Split-plot design with 4 blocks, 2 whole-plot fungicide treatments, and 70 barley varieties or variety mixes. Total area was 10 rows (north/south) by 56 beds (east/west).

Used with permission of Maria Durban.

Source

Durban, Maria and Hackett, Christine and McNicol, James and Newton, Adrian and Thomas, William and Currie, Iain. 2003. The practical use of semiparametric models in field trials, Journal of Agric Biological and Envir Stats, 8, 48-66. https://doi.org/10.1198/1085711031265.

Examples

## Not run: 

  library(agridat)
  data(durban.splitplot)
  dat <- durban.splitplot

  libs(desplot)
  desplot(dat, yield~bed*row,
          out1=block, out2=fung, num=gen, # aspect unknown
          main="durban.splitplot")


  # Durban 2003, Figure 2
  m20 <- lm(yield~gen + fung + gen:fung, data=dat)
  dat$resid <- m20$resid
  ## libs(lattice)
  ## xyplot(resid~row, dat, type=c('p','smooth'), main="durban.splitplot")
  ## xyplot(resid~bed, dat, type=c('p','smooth'), main="durban.splitplot")

  # Figure 4 doesn't quite match due to different break points
  libs(lattice)
  xyplot(resid ~ bed|factor(row), data=dat,
         main="durban.splitplot",
         type=c('p','smooth'))


  # Figure 6 - field trend
  # note, Durban used gam package like this
  # m2lo <- gam(yield ~ gen*fung + lo(row, bed, span=.082), data=dat)
  libs(mgcv)
  m2lo <- gam(yield ~ gen*fung + s(row, bed,k=45), data=dat)
  new2 <- expand.grid(row=unique(dat$row), bed=unique(dat$bed))
  new2 <- cbind(new2, gen="G01", fung="F1")
  p2lo <- predict(m2lo, newdata=new2)
  libs(lattice)
  wireframe(p2lo~row+bed, new2, aspect=c(1,.5),
            main="durban.splitplot - Field trend")

  if(require("asreml", quietly=TRUE)) {
    libs(asreml,lucid)
    
    # Table 5, variance components.  Table 6, F tests
    dat <- transform(dat, rowf=factor(row), bedf=factor(bed))
    dat <- dat[order(dat$rowf, dat$bedf),]
    m2a2 <- asreml(yield ~ gen*fung, random=~block/fung+units, data=dat,
                   resid =~ar1v(rowf):ar1(bedf))
    m2a2 <- update(m2a2)
    
    lucid::vc(m2a2)
    ##             effect component std.error z.ratio bound 
    ##              block   0              NA      NA     B  NA
    ##         block:fung   0.01206  0.01512      0.8     P   0
    ##              units   0.02463  0.002465    10       P   0
    ##       rowf:bedf(R)   1              NA      NA     F   0
    ## rowf:bedf!rowf!cor   0.8836   0.03646     24       U   0
    ## rowf:bedf!rowf!var   0.1261   0.04434      2.8     P   0
    ## rowf:bedf!bedf!cor   0.9202   0.02846     32       U   0
    
    wald(m2a2)
  }
  

## End(Not run)

Height of barley plants in a study of non-normal data

Description

Height of barley plants in a study of non-normal data.

Usage

data("eden.nonnormal")

Format

A data frame with 256 observations on the following 3 variables.

pos

position within block

block

block (numeric)

height

height of wheat plant

Details

This data was used in a very early example of a permutation test.

Eden & Yates used this data to consider the impact of non-normal data on the validity of a hypothesis test that assumes normality. They concluded that the skew data did not negatively affect the analysis of variance.

Grown at Rothamsted. Eight blocks of Yeoman II wheat. Sampling of the blocks was quarter-meter rows, four times in each row. Rows were selected at random. Position within the rows was partly controlled to make use of the whole length of the block. Plants at both ends of the sub-unit were measured. Shoot height is measured from ground level to the auricle of the last expanded leaf.

Source

T. Eden, F. Yates (1933). On the validity of Fisher's z test when applied to an actual example of non-normal data. Journal of Agric Science, 23, 6-17. https://doi.org/10.1017/S0021859600052862

References

Kenneth J. Berry, Paul W. Mielke, Jr., Janis E. Johnston Permutation Statistical Methods: An Integrated Approach.

Examples

## Not run: 

library(agridat)
data(eden.nonnormal)
dat <- eden.nonnormal
mean(dat$height) # 55.23 matches Eden table 1

# Eden figure 2
libs(dplyr, lattice)
# Blocks had different means, so substract block mean from each datum
dat <- group_by(dat, block)
dat <- mutate(dat, blkmn=mean(height))
dat <- transform(dat, dev=height-blkmn)

histogram( ~ dev, data=dat, breaks=seq(from=-40, to=30, by=2.5),
          xlab="Deviations from block means",
          main="eden.nonnormal - heights skewed left")

  # calculate skewness, permutation 
  
  libs(dplyr, lattice, latticeExtra)
  
  # Eden table 1
  # anova(aov(height ~ factor(block), data=dat))
  
  # Eden table 2,3. Note, this may be a different definition of skewness
  # than is commonly used today (e.g. e1071::skewness).
  skew <- function(x){
    n <- length(x)
    x <- x - mean(x)
    s1 = sum(x)
    s2 = sum(x^2)
    s3 = sum(x^3)
    k3=n/((n-1)*(n-2)) * s3 -3/n*s2*s1 + 2/n^2 * s1^3
    return(k3)
  }
  # Negative values indicate data are skewed left
  dat <- group_by(dat, block)
  summarize(dat, s1=sum(height),s2=sum(height^2), mean2=var(height), k3=skew(height))
  ##   block     s1       s2     mean2         k3
  ##   <int>  <dbl>    <dbl>     <dbl>      <dbl>
  ## 1     1 1682.0  95929.5 242.56048 -1268.5210
  ## 2     2 1858.0 111661.5 121.97984 -1751.9919
  ## 3     3 1809.5 108966.8 214.36064 -3172.5284
  ## 4     4 1912.0 121748.5 242.14516 -2548.2194
  ## 5     5 1722.0  99026.5 205.20565  -559.0629
  ## 6     6 1339.0  63077.0 227.36190  -801.2740
  ## 7     7 1963.0 123052.5  84.99093  -713.2595
  ## 8     8 1854.0 112366.0 159.67339 -1061.9919

  # Another way to view skewness with qq plot. Panel 3 most skewed.
  qqmath( ~ dev|factor(block), data=dat,
         as.table=TRUE,
         ylab="Deviations from block means",
         panel = function(x, ...) {
           panel.qqmathline(x, ...)
           panel.qqmath(x, ...)
         })

  # Now, permutation test.
  # Eden: "By a process of amalgamation the eight sets of 32 observations were
  # reduced to eight sets of four and the data treated as a potential
  # layout for a 32-plot trial".
  dat2 <- transform(dat, grp = rep(1:4, each=8))
  dat2 <- aggregate(height ~ grp+block, dat2, sum)
  dat2$trt <- rep(letters[1:4], 8)
  dat2$block <- factor(dat2$block)

  # Treatments were assigned at random 1000 times
  set.seed(54323)
  fobs <- rep(NA, 1000)
  for(i in 1:1000){
    # randomize treatments within each block
    # trick from https://stackoverflow.com/questions/25085537
    dat2$trt <- with(dat2, ave(trt, block, FUN = sample))
    fobs[i] <- anova(aov(height ~ block + trt, dat2))["trt","F value"]
  }

  # F distribution with 3,21 deg freedom
  # Similar to Eden's figure 4, but on a different horizontal scale
  xval <- seq(from=0,to=max(fobs), length=50)
  yval <- df(xval, df1 = 3, df2 = 21)
  # Re-scale, 10 = max of historgram, 0.7 = max of density
  histogram( ~ fobs, breaks=xval,
            xlab="F value",
            main="Observed (histogram) & theoretical (line) F values") +
    xyplot((10/.7)* yval ~ xval, type="l", lwd=2)


## End(Not run)

Potato yields in response to potash and nitrogen fertilizer

Description

Potato yields in response to potash and nitrogen fertilizer. Data from Fisher's 1929 paper Studies in Crop Variation 6. A different design was used each year.

Format

A data frame with 225 observations on the following 9 variables.

year

year/type factor

yield

yield, pounds per plot

block

block

row

row

col

column

trt

treatment factor

nitro

nitrogen fertilizer, cwt/acre

potash

potash fertilizer, cwt/acre

ptype

potash type

Details

The data is of interest to show the gradual development of experimental designs in agriculture.

In 1925/1926 the potato variety was Kerr's Pink. In 1927 Arran Comrade.

In the 1925a/1926a qualitative experiments, the treatments are O=None, S=Sulfate, M=Muriate, P=Potash manure salts. The design was a Latin Square.

The 1925/1926b/1927 experiments were RCB designs with treatment codes defining the amount and type of fertilizer used. Note: the 't' treatment was not defined in the original paper.

Source

T Eden and R A Fisher, 1929. Studies in Crop Variation. VI. Experiments on the response of the potato to potash and nitrogen. Journal of Agricultural Science, 19: 201-213.

References

McCullagh, P. and Clifford, D., (2006). Evidence for conformal invariance of crop yields, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science, 462, 2119–2143. https://doi.org/10.1098/rspa.2006.1667

Examples

## Not run: 

library(agridat)
data(eden.potato)
dat <- eden.potato

# 1925 qualitative
d5a <- subset(dat, year=='1925a')
libs(desplot)
desplot(d5a, trt~col*row,
        text=yield, cex=1, shorten='no', # aspect unknown
        main="eden.potato: 1925 qualitative")
anova(m5a <- aov(yield~trt+factor(row)+factor(col), d5a)) # table 2

# 1926 qualitative
d6a <- subset(dat, year=='1926a')
libs(desplot)
desplot(d6a, trt~col*row,
        text=yield, cex=1, shorten='no', # aspect unknown
        main="eden.potato: 1926 qualitative")
anova(m6a <- aov(yield~trt+factor(row)+factor(col), d6a)) # table 4

# 1925 quantitative
d5 <- subset(dat, year=='1925b')
libs(desplot)
desplot(d5, yield ~ col*row,
        out1=block, text=trt, cex=1, # aspect unknown
        main="eden.potato: 1925 quantitative")

# Trt 't' not defined, seems to be the same as 'a'
libs(lattice)
dotplot(trt~yield|block, d5,
        # aspect unknown
        main="eden.potato: 1925 quantitative")
anova(m5 <- aov(yield~trt+block, d5)) # table 6

# 1926 quantitative
d6 <- subset(dat, year=='1926b')
libs(desplot)
desplot(d6, yield ~ col*row,
        out1=block, text=trt, cex=1, # aspect unknown
        main="eden.potato: 1926 quantitative")
anova(m6 <- aov(yield~trt+block, d6)) # table 7

# 1927 qualitative + quantitative
d7 <- droplevels(subset(dat, year==1927))
libs(desplot)
desplot(d7, yield ~ col*row,
        out1=block, text=trt, cex=1, col=ptype, # aspect unknown
      main="eden.potato: 1927 qualitative + quantitative")

# Table 8.  Anova, mean yield tons / acre
anova(m7 <- aov(yield~trt+block+ptype + ptype:potash, d7))
libs(reshape2)
me7 <- melt(d7, measure.vars='yield')
acast(me7, potash~nitro, fun=mean) * 40/2240 # English ton = 2240 pounds
acast(me7, potash~ptype, fun=mean) * 40/2240


## End(Not run)

Uniformity trial of tea

Description

Uniformity trial of tea in Ceylon.

Usage

data("eden.tea.uniformity")

Format

A data frame with 144 observations on the following 4 variables.

entry

entry number

yield

yield

row

row

col

column