# ---------------------------------------------------------- 
# Exercise 2.8                                               
# ---------------------------------------------------------- 

data <- scan(what = list(signal="",delay=0), multi.line=F)
 pretimed 36.6 pretimed 39.2 pretimed 30.4 pretimed 37.1 pretimed 34.1  
 semiactuated 17.5 semiactuated 20.6 semiactuated 18.7 semiactuated 25.7 semiactuated 22.0
 fullyactuated 15.0 fullyactuated 10.4 fullyactuated 18.9 fullyactuated 10.5 fullyactuated 15.2 

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

traffic$signal <- as.factor(traffic$signal)
anova(lm(delay ~ signal, data=traffic))

# check model assumptions
par(mfrow=c(2,2))  
plot(lm(delay ~ signal, data=traffic))
par(mfrow=c(1,1))  

# calculate quantities as shown in the SAS code, then call:

powercr <- function(ngrp,rat,alpha,r1,rlast,rinc)
     {
        rnumber <- round( (rlast -r1)/rinc)
        r       <- numeric(rnumber)
        power   <- numeric(rnumber)
        for (i in seq(r1,rlast,by=rinc))
           {
              r[i]        <-  i
              ndf         <-  ngrp -1
              ddf         <-  (ngrp)*(i-1)
              fcr         <-  qf(1-alpha,ndf,ddf)
              lambda      <-  i*rat
              power[i]    <-  1 - pf(fcr,ndf,ddf,lambda) 
            }
        print(cbind(r,power))
        plot(r,power,type="b")
      }


powercr(3,.6963,.01,2,35,1)     


# ---------------------------------------------------------- 
# Exercise 4.1                                               
# ---------------------------------------------------------- 

data <- matrix(scan(),ncol=2,byrow=TRUE)
1520 1953 1520 2135 1520 2471 1520 4727 1520 6134 1520 6314
1620 1190 1620 1286 1620 1550 1620 2125 1620 2557 1620 2845
1660 651 1660 837 1660 848 1660 1038 1660 1361 1660 1543
1708 511 1708 651 1708 651 1708 652 1708 688 1708 729

data

Heaters <- data.frame(data)
colnames(Heaters) <- c("temp","hours")
rm(data)

Heaters$temp <- as.factor(Heaters$temp)  # to make sure that site is a factor, not numeric

boxplot(hours ~ temp, data=Heaters,main="Hermit Crab data")

anova(lm(hours ~ temp, data=Heaters))

 # part a

windows()
par(mfrow=c(2,2))  
plot(lm(hours ~ temp, data=Heaters))
par(mfrow=c(1,1))  


  # part b: the first approach from the text, using regression of log(mean) on log(std)

Heatmeans <- by(Heaters$hours,Heaters$temp,mean)

Heatstds <- by(Heaters$hours,Heaters$temp,sd)

lm(log(Heatstds) ~ log(Heatmeans))


  # the Box Cox method is available using the MASS library

library(MASS)

 # plotting the log-likelihood function for the Box-Cox transformation
 # start with a broader range
boxcox(hours ~ temp, data=Heaters,lambda = seq(-2.00, 2.00, length = 50))

 # zoom in for a closer look
boxcox(hours ~ temp, data=Heaters,lambda = seq(-1.5, -.5, length = 50))


   # part c: try a transformed value

Heaters$invhours <- 1/(Heaters$hours)         # inverse transformation

boxplot(invhours ~ temp, data=Heaters,main="Transformed Heater data")

anova(lm(invhours ~ temp, data=Heaters))

windows()
par(mfrow=c(2,2))  
plot(lm(invhours ~ temp, data=Heaters))
par(mfrow=c(1,1))  



# ---------------------------------------------------------- 
# Exercise 5.2                                               
# ---------------------------------------------------------- 

data <- matrix(scan(),ncol=2,byrow=TRUE)
 1 167.3  1 166.7   2 186.7  2 184.2  3 100.0  3 107.9 
 4 214.5  4 215.3   5 148.5  5 149.5  6 171.5  6 167.3 
 7 161.5  7 159.4   8 243.6  8 245.5 
   
data
Ex3.5 <- data.frame(data)
colnames(Ex3.5) <- c("patient","chol")
rm(data)

Ex3.5$patient <- as.factor(Ex3.5$patient)

library(nlme)

 # note: this gives estimates of the standard deviations, not variances
Ex3.5.lme  <- lme(chol ~ 1, data = Ex3.5, random = ~ 1 | patient )
summary(Ex3.5.lme)
intervals(Ex3.5.lme,level=.9)

 # get residual by predicted plot
plot(Ex3.5.lme)

 # ICC
 42.33^2/(42.33^2 + 2.45^2)



# ---------------------------------------------------------- 
# Exercise 5.4                                               
# ---------------------------------------------------------- 

data <- scan(what = list(lot="",conc=0), multi.line=F) 
  3469-72 39 3469-72 57 3469-72 63 3469-72 66 
  3849-52 56 3849-52 13 3849-52 25 3849-52 31 
  3721-24 64 3721-24 83 3721-24 88 3721-24 71 
  3477-80 29 3477-80 55 3477-80 21 3477-80 51 
  3669-72 38 3669-72 66 3669-72 53 3669-72 81 
  3873-76 11 3873-76 49 3873-76 34 3873-76 10 
  3777-80 23 3777-80 0 3777-80 5 3777-80 20 
  3461-64 10 3461-64 11 3461-64 23 3461-64 37 

data
Ex5.4 <- data.frame(data)
rm(data)

Ex5.4$lot <- as.factor(Ex5.4$lot)

library(nlme)

 # note: this gives estimates of the standard deviations, not variances
Ex5.4.lme  <- lme(conc ~ 1, data = Ex5.4, random = ~ 1 | lot )
summary(Ex5.4.lme)

 # get residual by predicted plot
plot(Ex5.4.lme)

 # total variance estimate
 20.77^2 + 15.2^2

 # ICC
  20.77^2/(20.77^2 + 15.2^2)

 # total standard deviation estimate
 sqrt(20.77^2 + 15.2^2)