# ---------------------------------------------------------------
 # Exercise 3.3
 # ---------------------------------------------------------------

data <- scan(what =list(treat="", moisture=0), multi.line=F)
          NaCl               80.5
          NaCl               79.3
          NaCl               79.0
    formicAcid               89.1
    formicAcid               75.7
    formicAcid               81.2
      beetPulp               77.8
      beetPulp               79.5
      beetPulp               77.0
       control               76.7
       control               77.2
       control               78.6

data

Ex3.3 <- data.frame(data)
rm(data)

Ex3.3$treat <- as.factor(Ex3.3$treat)
Ex3.3$treat <- relevel(Ex3.3$treat, ref="control")     # setting the reference level for Dunnett,
                                                       # must be defined before calling the lm function

Ex3.3.lm <- lm(moisture ~ treat, data=Ex3.3)

anova(Ex3.3.lm) 
summary(Ex3.3.lm)

 # check model assumptions

par(mfrow=c(2,2))  
plot(Ex3.3.lm)
par(mfrow=c(1,1))  

 # ---------------------------------------------------------------
 # Exercise 5.5 
 # ---------------------------------------------------------------


 # Dunnett's multiple comparison method

library(multcomp)

Ex3.3.Dunnett <- glht(Ex3.3.lm,linfct=mcp(treat="Dunnett"))

summary(Ex3.3.Dunnett)

confint(Ex3.3.Dunnett)

plot(Ex3.3.Dunnett,sub="Dunnetts method")


 # ---------------------------------------------------------------
 # Exercise 6.1 
 # ---------------------------------------------------------------


data <- scan(what =list(treat="", pctBluestem=0), multi.line=F)
         1-non                    97
         1-irr                    83
         2-non                    85
         3-non                    64
         4-non                    52
         4-irr                    48
         1-non                    96
         1-irr                    87
         2-non                    84
         3-non                    72
         4-non                    56
         4-irr                    58
         1-non                    92
         1-irr                    78
         2-non                    78
         3-non                    63
         4-non                    44
         4-irr                    49
         1-non                    95
         1-irr                    81
         2-non                    79
         3-non                    74
         4-non                    50
         4-irr                    53

data

Ex6.1 <- data.frame(data)
rm(data)

Ex6.1$treat <- as.factor(Ex6.1$treat)
Ex6.1$treat <- relevel(Ex6.1$treat, ref="1-non")     # setting the reference level for Dunnett
                                                     # must be defined before calling the lm function

Ex6.1.lm <- lm(pctBluestem ~ treat, data=Ex6.1)

anova(Ex6.1.lm) 
summary(Ex6.1.lm)

 # check model assumptions, part a

par(mfrow=c(2,2))  
plot(Ex6.1.lm)
par(mfrow=c(1,1))  

 # contrasts and multiple comparisons, parts c and d

library(lsmeans)

levels(Ex6.1$treat)

Ex6.1.lsm <- lsmeans(Ex6.1.lm, "treat")

Contrasts <- list(Quad = c(1, 0, -1, -1, 0, 1),IrrEffect=c(1,-1,0,0,-1,1))

contrast(Ex6.1.lsm,Contrasts,adjust="none")

 # part e

library(multcomp)

Ex6.1.CompBest <- glht(Ex6.1.lm,linfct=mcp(treat="Dunnett"),alternative="less")  # a one-sided Dunnett method is used
                                                                                # to compare to the best method

summary(Ex6.1.CompBest)

confint(Ex6.1.CompBest)

plot(Ex6.1.CompBest,sub="Comparison with the best method")


 # ---------------------------------------------------------------
 # Exercise 7.2 
 # ---------------------------------------------------------------

powercr <- function(ngrp,rat,alpha,n1,nlast,ninc)
     {
        n       <- seq(n1,nlast,by=ninc)
        power   <- numeric(length(n))
        ndf     <-  ngrp -1
        ddf         <-  ngrp*(n-1)
        fcr         <-  qf(1-alpha,ndf,ddf)
        lambda      <-  n*rat
       # print(c(n,ddf,fcr,lambda))
        power    <-  1 - pf(fcr,ndf,ddf,lambda) 
        print(cbind(n,ddf,lambda,fcr,power))
        plot(n,power,type="b",xlab="sample size per group",ylab="Power")
      }

 # for the homework problem using the powercr function

powercr(3,.1667,.05,2,100,1)   


 # or using the power.anova.test function

var(c(10, 11, 11))
 # = .3333, this is the between group variance

sampsz  <- seq(2,100,1)
powvals <- sapply(sampsz, function (x) power.anova.test(groups=3,within.var=4.0,
                                      between.var=.3333,sig.level=.05,n=x)$power)
plot(sampsz,powvals,type="b")

power.anova.test(groups=3,within.var=4.0,between.var=.3333,sig.level=.05,n=77)