ANCOVA in R

### Analysis of covariance, cricket example

Input = ("
         Species   Temp   Pulse
         ex       20.8   67.9
         ex       20.8   65.1
         ex       24     77.3
         ex       24     78.7
         ex       24     79.4
         ex       24     80.4
         ex       26.2   85.8
         ex       26.2   86.6
         ex       26.2   87.5
         ex       26.2   89.1
         ex       28.4   98.6
         ex       29    100.8
         ex       30.4   99.3
         ex       30.4  101.7
         niv      17.2   44.3
         niv      18.3   47.2
         niv      18.3   47.6
         niv      18.3   49.6
         niv      18.9   50.3
         niv      18.9   51.8
         niv      20.4   60
         niv      21     58.5
         niv      21     58.9
         niv      22.1   60.7
         niv      23.5   69.8
         niv      24.2   70.9
         niv      25.9   76.2
         niv      26.5   76.1
         niv      26.5   77
         niv      26.5   77.7
         niv      28.6   84.7
         ")

Data = read.table(textConnection(Input),header=TRUE)

#generalized model
model.1 = lm (Pulse ~ Temp + Species + Temp:Species, data = Data)

#ancova
model.2 = lm (Pulse ~ Temp + Species, data = Data)

#general linear F test
anova(model.1,model.2)

## Analysis of Variance Table
## 
## Model 1: Pulse ~ Temp + Species + Temp:Species
## Model 2: Pulse ~ Temp + Species
##   Res.Df    RSS Df Sum of Sq     F Pr(>F)
## 1     27 85.074                          
## 2     28 89.350 -1   -4.2758 1.357 0.2542

#reduced regression model
model.3 = lm (Pulse ~ Temp, data = Data)

#general linear F test
anova(model.2,model.3)

## Analysis of Variance Table
## 
## Model 1: Pulse ~ Temp + Species
## Model 2: Pulse ~ Temp
##   Res.Df    RSS Df Sum of Sq     F    Pr(>F)    
## 1     28  89.35                                 
## 2     29 687.35 -1      -598 187.4 6.272e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Final model is ANCOVA
summary(model.2)

## 
## Call:
## lm(formula = Pulse ~ Temp + Species, data = Data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0128 -1.1296 -0.3912  0.9650  3.7800 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.21091    2.55094  -2.827  0.00858 ** 
## Temp          3.60275    0.09729  37.032  < 2e-16 ***
## Speciesniv  -10.06529    0.73526 -13.689 6.27e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.786 on 28 degrees of freedom
## Multiple R-squared:  0.9896, Adjusted R-squared:  0.9888 
## F-statistic:  1331 on 2 and 28 DF,  p-value: < 2.2e-16

# Plot
I.nought = -7.21091
I1 = I.nought + 0
I2 = I.nought + -10.06529
B  = 3.60275

plot(x   = Data$Temp,
     y   = Data$Pulse,
     col = as.factor(Data$Species),
     pch = 16,
     xlab = "Temperature",
     ylab = "Pulse")

legend('bottomright',
       legend = levels(as.factor(Data$Species)),
       col = 1:2,
       cex = 1,   
       pch = 16)

abline(I1, B, lty=1, lwd=2, col = 1)

abline(I2, B, lty=1, lwd=2, col = 2)