# This file is posted as wheat8.ssc

#  Enter the data into a data frame and 
#  define factors

wheat <- read.table("c:\\courses\\st511\\snew\\wheatw.dat", 
   col.names=c("tray","moist","fert","y"))

# create factors while retaining the original
# columns in the data frame as continuous variables
# (non factors)

wheat$trayf <- as.factor(wheat$tray)
wheat$moistf <- as.factor(wheat$moist)
wheat$fertf <- as.factor(wheat$fert)

# Multistratum models may be fitted using aov() if
# the design is balanced and all factors are random. 
# Models are specified by a model formula
#
#  response ~ mean.formula + Error(strata.formular)

summary(wheat.aov <- aov(y ~ moistf*fertf + 
          Error(trayf/moistf),data=wheat))


# Use the lme function to fit the split-plot
# model with random block effects.  Here
# each fertilizer level is used once in each
# tray. 

options(contrasts=c(factor="contr.treatment",
                    ordered="contr.poly"))
wheat.lme <- lme( y ~ moistf*fertf , data=wheat,
    random = ~ 1 | trayf )

# Print F-tests for fixed effects
anova(wheat.lme)

# Print estimates of variance components
VarCorr(wheat.lme)

# Print estimates of fixed effects
fixef(wheat.lme)

# Print predictions of random effects
ranef(wheat.lme)

# Print approximate confidence intervals
intervals(wheat.lme)

# Plot residuals versus fitted values
par(fin=c(6,6),cex=1.5,mkh=1.2,mex=1.5)
plot(wheat.lme, resid(.,type="p") ~ fitted(.),
      adj=-0.3)

# Plot observations versus fitted values
plot(wheat.lme, y ~ fitted(.), adj=-0.3)

# Normal probability plot of within whole
# plot residuals
qqnorm( wheat.lme, ~resid(.))

# Normal probability plot of random
# effects
qqnorm( wheat.lme, ~ranef(.,level=1))

# Print covariance matrix for  estimates of fixed effects
wheat.lme$varFix

# Construct an approximate 95% confidence intervals
# for estimable functions of parameters.  Use 
# conservative degrees of freedom when the variance
# estimate is based on more than one MS.  You could
# add code for Satterthwaite degrees of freedom  


c1 <- matrix(c(1, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 
               0, 0, 0, 0, 0),nrow=1)
c1.est <- c1%*%wheat.lme$coef$fixed
c1.std <- sqrt(c1%*%wheat.lme$varFix%*%t(c1)) 
df.est <- 8
c1.est.lower <- c1.est - qt(.975,df.est)*c1.std
c1.est.upper <- c1.est + qt(.975,df.est)*c1.std

data.frame(c1.est.lower, c1.est, c1.est.upper)


c1 <- matrix(c(0, -1, 0, 0, 0, 0, 0, -.25, 0 , 0, -.25, 
               0, 0, -.25, 0, 0),nrow=1)
c1.est <- c1%*%wheat.lme$coef$fixed
c1.std <- sqrt(c1%*%wheat.lme$varFix%*%t(c1)) 
df.est <- 8
c1.est.lower <- c1.est - qt(.975,df.est)*c1.std
c1.est.upper <- c1.est + qt(.975,df.est)*c1.std

data.frame(c1.est.lower, c1.est, c1.est.upper)

# A within whole plot contrast. Change the 
# degrees of freedom


c1 <- matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 1 , 0, 0, 
               0, 0, 0, 0, 0),nrow=1)
c1.est <- c1%*%wheat.lme$coef$fixed
c1.std <- sqrt(c1%*%wheat.lme$varFix%*%t(c1)) 
df.est <- 36
c1.est.lower <- c1.est - qt(.975,df.est)*c1.std
c1.est.upper <- c1.est + qt(.975,df.est)*c1.std

data.frame(c1.est.lower, c1.est, c1.est.upper)


# Construct a profile plot


means <- tapply(wheat$y,list(wheat$moistf,wheat$fertf),mean)

#  Set up the axes and title of the profile plot

par(fin=c(7,7),cex=1.2,lwd=3,mex=1.5)
x.axis <- unique(wheat$moist)
matplot(c(10,40), c(0,16), type="n", 
     xlab="Moisture",  ylab="Mean Dry Weight", 
     main= "Wheat Plants") 

#  Add a profile for each soil type

matlines(x.axis,means,type='l',lty=c(1,3,5,7),lwd=3)     

# Plot points for the individual observations

matpoints(x.axis,means, pch=c(15,16,18,20))    

# Add a legend to the plot

legend(-5,-4,
   legend=c('Fert(2 mg)','Fert(4 mg)','Fert(6 mg)','Fert(8 mg)'),
          lty=c(1,3,5,7),bty='n')


means <- tapply(wheat$y,list(wheat$fertf,wheat$moistf),mean)

#  Set up the axes and title of the profile plot

par(fin=c(7,7),cex=1.2,lwd=3,mex=1.5)
x.axis <- unique(wheat$fert)
matplot(c(2,8), c(0,16), type="n", 
     xlab="Fertilizer",  ylab="Mean Dry Weight", 
     main= "Wheat Plants") 

#  Add a profile for each moisture level 

matlines(x.axis,means,type='l',lty=c(1,3,5,7),lwd=3)     

# Plot points for the individual observations

matpoints(x.axis,means, pch=c(15,16,18,20))    

# Add a legend to the plot

legend(-5,-4,
   legend=c('Moisture(10)','Moisture(20)','Moisture(30)','Moisture(40)'),
          lty=c(1,3,5,7),bty='n')


# Fit a two-dimensional quadratic model
# First center the explantory variables

wheat$df <- wheat$fert-mean(wheat$fert)
wheat$dm <- wheat$moist-mean(wheat$moist)
 

#  Fit the quadratic model

wheat.q <- lme( y ~ df+dm+df^2+dm^2+df*dm , data=wheat,
    random = ~ 1 | trayf )

# Print F-tests for fixed effects
anova(wheat.q)

# Print estimates of variance components
VarCorr(wheat.q)

# Print estimates of fixed effects
fixef(wheat.q)

# Print predictions of random effects
ranef(wheat.q)

# Print approximate confidence intervals
intervals(wheat.q)

# You should make some residual plots

# Estimate the mean dry weight when fert=5mg
# and moisture=15/ml/pot/day

meanf<-mean(wheat$fert)
meanm<-mean(wheat$moist)
c1 <- matrix(c(1, 5-meanf, 15-meanm, (5-meanf)^2, (15-meanm)^2,
               (5-meanf)*(15-meanm)),nrow=1)
c1.est <- c1%*%wheat.q$coef$fixed
c1.std <- sqrt(c1%*%wheat.q$varFix%*%t(c1)) 
df.est <- 9
c1.est.lower <- c1.est - qt(.975,df.est)*c1.std
c1.est.upper <- c1.est + qt(.975,df.est)*c1.std

data.frame(c1.est.lower, c1.est, c1.est.upper)


# Compare the fit of the quadratic model to the 
# model from part b