if (FALSE) {
library(agridat)
data(harris.multi.uniformity)
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
libs(reshape2,corrgram)
mat <- acast(dat, series+plot~harv, value.var='yield')
corrgram(mat, main="harris.multi.uniformity - correlation of crop yields")
# Compare to Harris 1928, table 4. More positive than negative correlations.
# densityplot(as.vector(cor(mat)), xlab="correlations",
# main="harris.multi.uniformity")
# 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))
# 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.
libs(desplot)
desplot(d2, yield ~ xord*plot|harv,
aspect=536/639, flip=TRUE, # true aspect
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(agg, yield ~ xord + plot,
aspect=536/639, # true aspect
main="harris.multi.uniformity fertility")
}
Run the code above in your browser using DataLab