dat <- harris.multi.uniformity
# Combine year/crop into 'harvest'
dat <- transform(dat, harv = factor(paste0(year,".",crop)))
# Convert 1911 from tons to pounds
dat$yield[dat$year==1911] <- 340 * dat$yield[dat$year==1911]
# Average yields. Harris 1928, table 2.
aggregate(yield~harv, dat, mean)
# Corrgram
require(reshape2)
mat <- acast(dat, series+plot~harv, value.var='yield')
if(require(corrgram)) corrgram(mat)
# Compare to Harris 1928, table 4. More positive than negative correlations.
densityplot(as.vector(cor(mat)))
# Standardize yields for each year
mats <- scale(mat)
# Melt and re-name columns so we can make field maps. Obvious spatial
# patterns that persist over years
d2 <- melt(mats)
names(d2) <- c('ord','harv','yield')
d2$series <- as.numeric(substring(d2$ord,1,1))
d2$plot <- as.numeric(substring(d2$ord,3))
# True field shape: aspect=.84 (536 ft North-South by 639 ft E-W)
# Series 2 is on the east side, so switch 2 and 3 for correct plotting
d2$xord <- 5 - dat$series
# Note that for alfalfa, higher-yielding plots in 1912-1914 were
# lower-yielding in 1922-1923.
desplot(yield ~ xord*plot|harv, data=d2, aspect=.84, flip=TRUE,
main="harris.multi.uniformity")
# Crude fertility map by averaging across years shows probable
# sub-surface water effects
agg <- aggregate(yield ~ xord + plot, data=d2, mean)
desplot(yield ~ xord + plot, agg, main="harris.multi.uniformity fertility")Run the code above in your browser using DataLab