# NOT RUN {
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
if(require(reshape2) & require(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")
# }
# NOT RUN {
# 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.
if(require(desplot)){
desplot(yield ~ xord*plot|harv, d2,
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(yield ~ xord + plot, agg,
aspect=536/639, # true aspect
main="harris.multi.uniformity fertility")
}
# }
# NOT RUN {
# dontrun
# }
Run the code above in your browser using DataLab