data <- scan(what = list(row=0,col=0,treatment="",yield=0), multi.line=F)
  1 1 E 59.45  1 2 A 47.28  1 3 C 54.44  1 4 B 50.14  1 5 D 59.45 
  2 1 C 55.16  2 2 D 60.89  2 3 B 56.59  2 4 E 60.17  2 5 A 48.71 
  3 1 B 44.41  3 2 C 53.72  3 3 D 55.87  3 4 A 47.99  3 5 E 59.45 
  4 1 A 42.26  4 2 B 50.14  4 3 E 55.87  4 4 D 58.74  4 5 C 55.87 
  5 1 D 60.89  5 2 E 59.45  5 3 A 49.43  5 4 C 59.45  5 5 B 57.31 

data


table87   <- data.frame(data)
rm(data)

table87$row         <- as.factor(table87$row)
table87$col         <- as.factor(table87$col)
table87$treatment   <- as.factor(table87$treatment)


table87ls  <- lm(yield ~ row +col +treatment, data = table87)

anova(table87ls)

library(phia)

 # polynomial contrasts

lineartrt.contrast <- list(treatment = c(-2,-1,0,1,2))
quadtrt.contrast <- list(treatment = c(2,-1,-2,-1,2))

testInteractions(table87ls ,custom=lineartrt.contrast)
testInteractions(table87ls ,custom=quadtrt.contrast)

 # means and standard errors

interactionMeans(table87ls,factors="treatment")

 # check model assumptions

par(mfrow=c(2,2))  
plot(table87ls)
par(mfrow=c(1,1))  


 # observed means and standard errors of differences of observed means 
 # (rarely of interest for unbalanced designs)

model.tables(aov(table87ls),"means",se=T)