---
title: "Using Papadakis covariates for nearest neighbor analysis"
author: "Kevin Wright"
date: "`r Sys.Date()`"
output:
  pdf_document: default
  html_document:
    highlight: tango
bibliography: agridat.bib
vignette: |
  % \VignetteIndexEntry{Using Papadakis covariates for nearest neighbor analysis} 
  % \VignetteEncoding{UTF-8} 
  % \VignetteEngine{knitr::rmarkdown}
---
<style type="text/css">
body, td { font-size: 14pt; color: #000; padding:1em; }
code, pre, #ws, #message { font-size: 12pt;
  color: #000; background-color: #f5f5f5; }
code {border: 1px solid #e0e0e0; padding: 0 5px;}
pre > code { border: 0; }
</style>

# R setup

```{r setup, message=FALSE, echo=FALSE}
knitr::opts_chunk$set(fig.width=6, fig.height=6)
options(width=90)
```
@papadakis1937method believed that traditional blocking in field experiments
did not adequately represent the patchiness of soil fertility patterns and he
instead proposed adjusting the yield of each plot by the performance of the
neighboring plots.
  
If there is heterogeneity in the field that is of a scale smaller than the
block (but larger than the individual plots) then adjacent plots will be
positively correlated and this information about the neighboring plots can
be used to reduce the effect of spatial heterogeneity and increase the
accuracy of the treatment effects.
  
The Papadakis method is a nearest neighbor method that uses a residual
covariate in the analysis. In essence, the method follows the following steps.

1. Fit a treatment model and calculate the residuals from the model.

2. Calculate covariates that are the average of the neighboring residuals.

3. Fit a model with additional covariate terms for the residuals.

The left-right (LR) covariate for the (i,j)th plot is the
average of the residuals for the plots immediately to the left and
right of the (i,j)th plot.  If one of these neighbors is missing, then
the covariate is constructed from the single remaining neighboring
residual.  Border plots use only one neighboring residual.  The up-down (UD)
covariate is similarly constructed from residuals for plots
immediately up or down from the (i,j)th plot.

```{r func}
papcov <- function(resid,x,y){

  # Make sure x and y are numeric
  if(is.factor(x)) x <- as.numeric(as.character(x))
  if(is.factor(y)) y <- as.numeric(as.character(y))
  xy <- paste(x,y,sep=":")

  # Average neighboring residuals in up/down direction
  xym1 <- paste(x,y-1,sep=":")
  xyp1 <- paste(x,y+1,sep=":")
  rm1 <- resid[match(xym1,xy,NA)]
  rp1 <- resid[match(xyp1,xy,NA)]
  ud <- (rm1+rp1)/2
  # If only one neighboring residual is available, then just use it
  ud <- ifelse(is.na(ud) & !is.na(rm1),rm1,ud)
  ud <- ifelse(is.na(ud) & !is.na(rp1),rp1,ud)

  # Average neighboring residuals in left/right direction
  xm1y <- paste(x-1,y,sep=":")
  xp1y <- paste(x+1,y,sep=":")
  cm1 <- resid[match(xm1y,xy,NA)]
  cp1 <- resid[match(xp1y,xy,NA)]
  lr <- (cm1+cp1)/2
  # If only one neighboring residual is available, then just use it
  lr <- ifelse(is.na(lr) & !is.na(cm1),cm1,lr)
  lr <- ifelse(is.na(lr) & !is.na(cp1),cp1,lr)

  return(list(LR=lr, UD=ud))
}
```

# Reproduce Hinz 1987 case 2

@hinz1987nearest used the Papadakis technique to analzye a field experiment of
tobacco.

```{r hinz}
library("agridat")
data(federer.tobacco)
dat <- federer.tobacco
dat <- transform(dat, height=height-600) # For simplicity

# Model 1 - RCB
m1 <- aov(height ~ factor(block) + factor(dose), dat)
anova(m1)

# Model 2 - Row/Col as class variables
m2 <- aov(height ~ factor(block) + factor(dose) + factor(row), dat)
anova(m2)

# Model 3 - Two-step Papadakis
m3 <- aov(height ~ factor(dose), dat)
dat <- cbind(dat, papcov(m3$resid, dat$block, dat$row))
m4 <- aov(height ~ factor(dose) + LR + UD, data=dat)
anova(m4)
# Resid MS uses 1 df less to account for covariates.  Matches Hinz.
582073 / 46 
```

# Iterated example as given in Stroup et al, Table 2

@stroup1994removing used the Papadakis tecnique in an iterative manner.

```{r stroup}
library("agridat")
data(stroup.nin)
dat2 <- stroup.nin
dat2 <- subset(dat2,!is.na(yield))
n.gen <- nlevels(dat2$gen)

# RCB model, ranks match Stroup Table 2, RCB Alliance
m5 <-  lm(yield ~ gen -1 + rep, data=dat2)
pred.rcb <- coef(m5)[1:n.gen] # RCB adj means
rev(57-sort(rank(pred.rcb)))

# Initial genotype model (no blocks)
m6 <-  lm(yield ~ gen -1, data=dat2)
# Calculate Papadakis covariates
pp <- papcov(resid(m6), dat2$col, dat2$row)
dat2$LR <- pp$LR
dat2$UD <- pp$UD
# Single iteration of Papadakis model
m7 <- lm(yield ~ gen - 1 + LR + UD, data=dat2)
# Papadakis adjusted means
adjmn <- coef(m7)[1:n.gen]
# Residual = observed - adjusted mean
resid <- dat2$yield - adjmn[match(paste0("gen",dat2$gen),names(adjmn))]

# Now iterate Papadakis method to convergence
iter <- 0
notConv <- TRUE
while(notConv){
  iter <- iter + 1
  # Covariates based on residuals
  pp <- papcov(resid, dat2$col, dat2$row)
  dat2$LR <- pp$LR
  dat2$UD <- pp$UD
  m8 <- lm(yield ~ gen - 1 + LR + UD, data=dat2)
  # Check convergence of adjusted means
  prevmn <- adjmn
  adjmn <- coef(m8)[1:n.gen]
  tol <- sum((adjmn - prevmn)^2)
  cat("Iteration: ",iter," tol: ",tol,"\n")
  notConv <- tol > .001
  resid <- dat2$yield - adjmn[match(paste0("gen",dat2$gen),names(adjmn))]
}
pred.pap <- adjmn

# Ranks almost match Stroup et al, Table 2, Alliance, RCB+NNA-PAP
all <- data.frame(rcb=57-rank(pred.rcb), nna=57-rank(pred.pap))
all[order(all$rcb),]

# Visually compare the coefficients from the two methods
lims=range(c(pred.rcb,pred.pap))
plot(pred.rcb, pred.pap, xlim=lims, ylim=lims,
     xlab="RCB",ylab="Papadakis",type='n')
text(pred.rcb,pred.pap, substring(names(pred.rcb),4),cex=0.5)
title("Iterated Papadakis vs. RCB")
abline(0,1)
```
The variety 'Buckskin' has a large adjustment when using the Papadakis
method.  This makes sense, because Buckskin ended up in the high-yielding part
of each of the four reps.

# References
