# seedling analysis

seedlingdata_read.table("seedlingdata.asc",header=T)
seedlingdata$hdiff_seedlingdata$h8-seedlingdata$h1
seedlingdata$pdiff_(seedlingdata$h8-seedlingdata$h1)/seedlingdata$h1
attach(seedlingdata)
seed_seedlingdata

# a plot of the data...
yl_range(seed[,4:11],na.rm=T)
par(mfcol=c(4,4),oma=c(1,5,5,1),mar=c(5,5,2,2))
for(k in 1:4){
  ok_seed$species=="jp" & seed$blk==k & seed$trt=="D"
  matplot(1:8,t(seed[ok,4:11]),type="l",col=1,lty=1,ylim=yl,
  xlab="week",ylab="height in cm")
}
for(k in 1:4){
  ok_seed$species=="jp" & seed$blk==k & seed$trt=="W"
  matplot(1:8,t(seed[ok,4:11]),type="l",col=1,lty=1,ylim=yl,
  xlab="week",ylab="height in cm")
}
for(k in 1:4){
  ok_seed$species=="as" & seed$blk==k & seed$trt=="D"
  matplot(1:8,t(seed[ok,4:11]),type="l",col=1,lty=1,ylim=yl,
  xlab="week",ylab="height in cm")
}
for(k in 1:4){
  ok_seed$species=="as" & seed$blk==k & seed$trt=="W"
  matplot(1:8,t(seed[ok,4:11]),type="l",col=1,lty=1,ylim=yl,
  xlab="week",ylab="height in cm")
}
mtext(c("block 1","block 2","block 3","block 4"), side=2, 
 at=c(282.63547,199.55112,117.02814, 34.50518), line=1,cex=1.2,
 outer=T,srt=90)
mtext(c("stressed","control","stressed","control"), side=3, 
 at=c(-25.726543,-15.718289, -5.656226,  4.728684), line=0,cex=0.9,
 outer=T,srt=0)
mtext(c("Jackpine","Aspen"), side=3, 
 at=c(-20.7, -.65), line=2,cex=1.2, outer=T,srt=0)

# fitted means
hdmat_matrix(seed$hdiff,nrow=5)
mfit_apply(hdmat,2,mean,na.rm=T)

# make a plot of the diffs by block and treatment...
par(mfrow=c(2,2),mar=c(5,5,4,2))
plot(as.numeric(seed$blk)-.2,as.numeric(seed$hdiff),type="n",xlab="Block",ylab="change in height",
  xlim=c(.8,4.2))
ok_seed$trt=="D" & seed$species=="jp"
points(seed$blk[ok]-.2,seed$hdiff[ok],pch="D")
ok_seed$trt=="W" & seed$species=="jp"
points(seed$blk[ok]+.2,seed$hdiff[ok],pch="W")
mtext("Jackpines",side=3,line=1)
matlines(rbind(1:4-.2,1:4+.2),
  matrix(c(mfit[1:2],mfit[5:6],mfit[9:10],mfit[13:14]),nrow=2),
  lty=1,col=1,lwd=2)
matpoints(rbind(1:4-.2,1:4+.2),
  matrix(c(mfit[1:2],mfit[5:6],mfit[9:10],mfit[13:14]),nrow=2),
  pch=9)

plot(as.numeric(seed$blk)-.2,as.numeric(seed$hdiff),type="n",xlab="Block",
  ylab="change in height",xlim=c(.8,4.2))
ok_seed$trt=="D" & seed$species=="as"
points(seed$blk[ok]-.2,seed$hdiff[ok],pch="D")
ok_seed$trt=="W" & seed$species=="as"
points(seed$blk[ok]+.2,seed$hdiff[ok],pch="W")
mtext("Aspen",side=3,line=1)
matlines(rbind(1:4-.2,1:4+.2),
  matrix(c(mfit[3:4],mfit[7:8],mfit[11:12],mfit[15:16]),nrow=2),
  lty=1,col=1,lwd=2)
matpoints(rbind(1:4-.2,1:4+.2),
  matrix(c(mfit[3:4],mfit[7:8],mfit[11:12],mfit[15:16]),nrow=2),
  pch=9)

# fit a simple AOV to the Jackpine data...

seedlingdata$blk_as.factor(seedlingdata$blk)
 aov(formula = hdiff ~ blk + trt + trt:blk, 
     data = seedlingdata, subset = species ==
     "jp", na.action = na.omit)

Terms:
                     blk      trt  trt:blk Residuals 
 Sum of Squares  19.5608  43.4723  24.7647  257.4320
Deg. of Freedom        3        1        3        32

Residual standard error: 2.836327 
Estimated effects are balanced

          Df Sum of Sq  Mean Sq  F Value     Pr(F) 
      blk  3   19.5608  6.52025 0.810498 0.4974680
      trt  1   43.4723 43.47225 5.403804 0.0265970
  trt:blk  3   24.7647  8.25492 1.026125 0.3940754
Residuals 32  257.4320  8.04475                   

# shows a treatment effect!...


# fit a simple AOV to the Aspen data...

 aov(formula = hdiff ~ blk + trt + trt:blk, 
     data = seedlingdata, subset = species == "as", 
     na.action = na.omit)

Terms:
                     blk      trt  trt:blk Residuals 
 Sum of Squares  342.868    8.526  203.267  1798.918
Deg. of Freedom        3        1        3        30

Residual standard error: 7.743638 
Estimated effects may be unbalanced

          Df Sum of Sq  Mean Sq  F Value     Pr(F) 
      blk  3   342.868 114.2893 1.905968 0.1499123
      trt  1     8.526   8.5263 0.142191 0.7087671
  trt:blk  3   203.267  67.7556 1.129939 0.3526952
Residuals 30  1798.918  59.9639                   

# no sign of any treatment effect...

Perhaps the simplest approach is to reduce the data to a single
measure of growth over this 8 week period.  Using h8-h1 works
fine, as does something like (h8-h1)/h1.

Both the ANOVA models and the plot show that the Jackpines
appear to have a treatment effect.  The Aspens show no treatment
effect.  
Est and CI for the treatment effect: 
Jackpines:  1.04cm  CI:(0.15,1.94)
Aspens:     0.5 cm  CI:(-2,3)   

It's important not to fit the data with a single model with a
common error variance since the variability appears to be rather
different for the two species.  Fitting such a model will obscure
the treatment effect in the Jackpines.

Since the data are balanced (or nearly so), the ANOVA inference
will be robust to non-normality.  So the inference is fairly
trustworthy here.

Also, there is no obvious pattern in the growth trajectories when
comparing the D to the W groups.