### Lecture 37 - Code for running ANCOVA in R


### import data
datum=read.csv(file.choose())
head(datum)

### Demonstration of 'swamping' 

# Analyze effect of sex by itself
results=lm(Size~Sex,data=datum)
summary(results)
# Note that effect of sex is not significant, ability to detect effect of sex is
# 'swamped' out by effect of age

### Visual depiction of effect of swamping
plot(datum$Sex,datum$Size) # Plot of relationship between sex and size
plot(datum$Age,datum$Size) # plot of relationship between age and size

# color code points to see differences between sexes
points(datum$Age[1:31],datum$Size[1:31],col="blue")
points(datum$Age[32:61],datum$Size[32:61],col="red")


### Run ANCOVA using lm to get effect of age and sex on size

results=lm(Size~Age+Sex,data=datum)
summary(results)
confint(results)



### example 2

### import data
datum=read.csv(file.choose())
head(datum)

### Run ANCOVA using lm
results=lm(Biomass~Rainfall+Habitat,data=datum)
summary(results)
confint(results)

###Is 'Habitat' significant as a whole (ANOVA) ?

results2=lm(Biomass~Rainfall,data=datum) #reduced model without habitat
anova(results,results2) #test for whether habitat significantly improves the fit of the model



####Predictions

#When making predictions, every x variable in model must be in new data
NewData=data.frame(Rainfall=10, Habitat="Pine")
predict(results,newdata=NewData,interval="prediction")



####Change of reference so that you can compare 'pine' to 'hardwood'

datum$Habitat=as.factor(datum$Habitat)
results3=lm(Biomass~Rainfall+relevel(Habitat,ref="Pine"),data=datum)
#Note that rainfall is still in model
#"Habitat" is replaced with the relevel function
summary(results3)