# 

setwd('U:/public_html/588/Class17')    
X <- read.csv('bats.csv', header = T)

attach(X)
lmass <- log(MASS)
lener <- log(ENERGY)
 n <- length(MASS)

 plot(MASS, ENERGY, pch = as.numeric(TYPE), cex = 2)


 plot(lmass, lener, pch = as.numeric(TYPE), cex = 2)

 lm1 <- lm(lener ~ TYPE)     # this is, of course, the usual ANOVA analysis
 anova(lm1) 

 (Types <- levels(TYPE))
 bird <- as.numeric(TYPE == Types[3])   # this sets up indicator/"dummy" variables
 ebat <- as.numeric(TYPE == Types[1])   

 lm0 <- lm(lener ~ 1)  # this is the "empty" model containing only beta_0
 lm2 <- lm(lener ~ bird + ebat)
  summary(lm2) 
 
 anova(lm0 ,lm2)     # this is the "Extra SS" F-test which of course gives the same answer 

## now we add the 'lmass' variable

 lm3a <- lm(lener ~ lmass + bird + ebat)  # this is how it's done in Chapter 10
  summary(lm3a)
 lm3b <- lm(lener ~ lmass + TYPE)          # this is a more convenient ANCOVA interface
  summary(lm3b)
  anova(lm3b)      # in the ANOVA table, they pulled two dummy coefficients into one 'TYPE' group


  lm4 <- lm(lener ~ lmass + TYPE + lmass:TYPE)
  anova(lm4)


 SST <- sum((lener - mean(lener))^2)   # this is the SS for the lm0 model
 SSE <- 0.5533                         # this is the REsidual SS for the lm3a or lm3b
 MSE <- SSE/(n-4)
 MST <- SST/(n-1)

 (  R2 <- (SST - SSE)/SST   )
 (  R2adj <- (MST - MSE)/MST   )


 # predictive intervals for, say, birds

 xrange <- seq(10,800,10)
 lrange <- log(xrange)
 gridbirds <- data.frame( lmass = lrange, bird = rep(1,n), ebat = rep(0,n) )
           # i.e. put the bird indicator at 1
 
 pred <- predict(lm3a, gridbirds)
 pred.plim <- predict(lm3a, gridbirds, interval="prediction")
 pred.clim <- predict(lm3a, gridbirds, interval="confidence")
 
matplot(xrange, exp(cbind(pred.clim, pred.plim[,-1])), type="l")
points(MASS, ENERGY, pch = as.numeric(TYPE) )