analysis_of_covariance_ANCOVA
Chapter 6 ANCOVA Linear Regression Practice Script This script is for practicing using a analysis of covariance (ANCOVA) from chapter 6 of the Beckerman textbook.
See https://rpubs.com/jorgelizarazo/Rdemo5_wekk5_chapter6 for lab demo
ANCOVA combines a categorical explanatory variable with a continuous explanatory variable. We are essentially combining regression and one-way ANOVA
Clear R’s brain:
Import libraries:
##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Warning: package 'ggfortify' was built under R version 4.6.1Import the dataset and call it limp then view it:
limp <- read.csv("/Users/jroweibekwe/Library/Mobile Documents/com~apple~CloudDocs/MSc/Bio 9915/Bio_9915_R_Practice/textbook_data_9915/limpet.csv")
glimpse(limp)## Rows: 24
## Columns: 3
## $ DENSITY <int> 8, 8, 8, 8, 8, 8, 15, 15, 15, 15, 15, 15, 30, 30, 30, 30, 30, …
## $ SEASON <chr> "spring", "spring", "spring", "summer", "summer", "summer", "s…
## $ EGGS <dbl> 2.875, 2.625, 1.750, 2.125, 1.500, 1.875, 2.600, 1.866, 2.066,…The data relates egg production by limpets to four density conditions in two seasons. The response variable (y) is egg production (EGGS) and the independent variables are DENSITY (continuous) and SEASON (categorical).
This analysis is trying to solve if the density dependence of egg production differs between spring and summer?
#Creating a plot: The response (y) variable is egg production, DENSITY is the continuous independent variable (x), and the values are also classified by the category SEASON. This suggests that a plot of egg production against density, distinguishing the two seasonal categories, is what we need.
ggplot(limp, aes(x = DENSITY, y = EGGS, colour = SEASON)) +
geom_point() +
scale_color_manual(values = c(spring="green", summer="red")) + #lets you choose the exact colours assigned to categories in your plot
theme_bw()
On our graph, if we were to describe the relationship between egg
production and density by a straight line (ignore the season for the
moment), then y is egg production, x is density, b is where the line
crosses the y-axis (i.e. the egg production at 0 density—admittedly an
odd concept!), and m is the slope of the egg production density
relationship. This slope represents the change in egg production per
unit change in density.
Overall, egg numbers decline with increasing density—that is, there is a negative slope. Second, there may be a seasonal difference(the intercept) the value of egg production at zero density, is different in each season.
Possible slopes: 1. The number of eggs produced by limpets does not vary with density or season. *Common intercept, Zero slope(s)
The number of eggs produced by limpets does not vary with density but is reduced in the summer season (dashed). *Different intercepts, Zero slope(s), Parallel horizontal line
The number of eggs produced declines with density, but the maximum number of eggs (intercept) and the rate (slope) do not vary with season *Same intercept, Same slope, Same lines
The number of eggs produced declines with density and the number of eggs at density zero (intercept) differs between seasons but the rate (slope) does not vary with season.
- Different intercepts, Same (negative) slope, Parallel lines
- The number of eggs at density zero (intercept) and the rate (slope) vary with season *Different intercepts, Different slopes - ‘Crossing’ lines
#Constructing the ANCOVA: Since we want to include an effect of DENSITY (main effect), an effect of SEASON (main effect), and the potential for the effect of DENSITY to depend on SEASON (interaction). The model that captures this best is:
limp.mod <- lm(EGGS ~ DENSITY * SEASON, data = limp)
#lm creates linear models
#the * symbol between the two explanatory variables## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "contrasts" "xlevels" "call" "terms"
## [13] "model"#R has stored for us the coeffcients (intercept and slope estimates), the residuals, the fitted values, and a number of other values. #Check the Model assumptions:
## Warning: `fortify(<lm>)` was deprecated in ggplot2 4.0.0.
## ℹ Please use `broom::augment(<lm>)` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
## Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
## Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).
Interpreting the Model Output: anova() produces an ANOVA table listing sums of squares, mean squares, F-values, and p-values. It does not perform an ANOVA. summary() returns a table of coefficients (slopes and intercepts), standard errors, and t-values.
## Analysis of Variance Table
##
## Response: EGGS
## Df Sum Sq Mean Sq F value Pr(>F)
## DENSITY 1 5.0241 5.0241 30.1971 2.226e-05 ***
## SEASON 1 3.2502 3.2502 19.5350 0.0002637 ***
## DENSITY:SEASON 1 0.0118 0.0118 0.0711 0.7925333
## Residuals 20 3.3275 0.1664
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#This table reveals that there is no additional, significant variation in egg number explained by allowing different slopes for each season (since the p-value for the DENSITY:SEASON row is 0.79). We fail to reject the null hypothesis that the slopes are the same, i.e. we have failed to find evidence for the alternative hypothesis.
#Our result could actually be written:
#We tested the hypothesis that the effect of density on egg production in limpets depends on the season in which they are reproducing. We found no evidence for an interaction between density and season (F = 0.0711; df = 1, 20; p = 0.79), indicating that the effects of density and season are additive.##
## Call:
## lm(formula = EGGS ~ DENSITY * SEASON, data = limp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.65468 -0.25021 -0.03318 0.28335 0.57532
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.664166 0.234118 11.380 3.45e-10 ***
## DENSITY -0.033650 0.008259 -4.074 0.000591 ***
## SEASONsummer -0.812282 0.331092 -2.453 0.023450 *
## DENSITY:SEASONsummer 0.003114 0.011680 0.267 0.792533
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4079 on 20 degrees of freedom
## Multiple R-squared: 0.7135, Adjusted R-squared: 0.6705
## F-statistic: 16.6 on 3 and 20 DF, p-value: 1.186e-05#Provides a statement of your model by repeating your usage of lm()—the call.
#Next it provides the range, interquartile range, and median of the residuals as a quick, additional check of the residuals.
#The coefficient table provides insight into the estimates for the lines that you have fitted—these are the coeffcients that quantify the important biological relationships.#PUTTING THE LINES ONTO THE FIGURE: To start this process, we invoke two functions: expand.grid() and predict().
expand.grid() is a function that generates a grid of numbers — essentially it builds a factorial representation of any variables you provide to it, and returns a data frame
predict() is a function that generates fitted values from a model. It requires at least one argument, but works most effectively for us with three: a model, a set of new ‘x’ (explanatory) values at whichwewant to know ‘y’ values, generated according to our model, and another argument, interval, which provides a rapid way to get the 95% confidence interval, doing some hard work for us.
## FIRST SECOND
## 1 A 1
## 2 B 1
## 3 A 2
## 4 B 2## 1 2 3 4 5 6 7 8
## 2.3949692 2.3949692 2.3949692 1.6075953 1.6075953 1.6075953 2.1594217 2.1594217
## 9 10 11 12 13 14 15 16
## 2.1594217 1.3938428 1.3938428 1.3938428 1.6546769 1.6546769 1.6546769 0.9358016
## 17 18 19 20 21 22 23 24
## 0.9358016 0.9358016 1.1499321 1.1499321 1.1499321 0.4777604 0.4777604 0.4777604First, we create the ‘new x’s’ (DENSITY) at which we want values of y (EGGS). Here we use expand.grid() to generate a set of ‘new x’ values, and we take special care to label them as the column name in the original dataset, DENSITY.
We know the density varies from 8 to 45, so let’s use the function seq() to get ten numbers in between 8 and 45:
to check that it worked:
## DENSITY
## 1 8.00000
## 2 12.11111
## 3 16.22222
## 4 20.33333
## 5 24.44444
## 6 28.55556## [1] "factor"Let’s add SEASON to this grid now. SEASON has two levels, spring and summer. We can add these levels to the ‘new x’s’ as follows. Note how we use the function levels(), against the dataset, to extract the information we desire:
new.x <- expand.grid(DENSITY = seq(from = 8, to = 45, length.out = 10), SEASON = levels(limp$SEASON))To check if it worked:
## DENSITY SEASON
## 1 8.00000 spring
## 2 12.11111 spring
## 3 16.22222 spring
## 4 20.33333 spring
## 5 24.44444 spring
## 6 28.55556 springuse predict() with three arguments: a model, a value for newdata, and confidence intervals. We assign what predict returns to an object called new.y:
To check if it worked:
## fit lwr upr
## 1 2.394969 2.019285 2.770654
## 2 2.256632 1.931230 2.582034
## 3 2.118294 1.834274 2.402315
## 4 1.979957 1.724062 2.235852
## 5 1.841619 1.595998 2.087241
## 6 1.703282 1.447918 1.958646Housekeeping step: combine the new y’s with the new x’s, so that we have a clear picture of what it is you’ve made, using the function data.frame(). Call the data frame addThese, because we are gonna add these to the plot. Rename the fit produced by predict to EGGS to match the original data using rename() - dplyr function.
To check if it worked:
## DENSITY SEASON EGGS lwr upr
## 1 8.00000 spring 2.394969 2.019285 2.770654
## 2 12.11111 spring 2.256632 1.931230 2.582034
## 3 16.22222 spring 2.118294 1.834274 2.402315
## 4 20.33333 spring 1.979957 1.724062 2.235852
## 5 24.44444 spring 1.841619 1.595998 2.087241
## 6 28.55556 spring 1.703282 1.447918 1.958646This creates a data frame that contains the grid of seasons and densities, as well as the predicted values and 95% CI at each of these combinations.
Before we finish with adding these data to the figure, we want to show how flexible this method is. Let’s modify the code so that we get the predicted number of eggs for each season at the mean limpet density:
#Making the plot:
ggplot(limp, aes(x = DENSITY, y = EGGS, colour = SEASON)) + # first add the points
geom_point(size = 5) + geom_smooth(data = addThese, #geom_smooth() adds the fits and CIs
aes(ymin = lwr, ymax = upr, # note we don't need to specify DENSITY AND EGG
fill = SEASON), stat = 'identity') + #stat = ‘identity’ tells geom_smooth() not to try and calculate anything fancy-shmancy for us
scale_colour_manual(values = c(spring="green", summer="red")) + #now adjust the colours
scale_fill_manual(values = c(spring="green", summer="red")) +
# theme it
theme_bw()
#In Summary: 1. Enter your data in a spreadsheet program, check it, and save it as a comma-separated values file. 2. Start R, wipe its brain of all previous data — rm(list=ls())—and get your libraries sorted. 3. Import your data into R and check this has worked, for example using glimpse(). 4. Make a picture that reflects your questions and that will help you specify your model to answer your question. 5. Specify a statistical model to capture the hypothesis you are testing. 6. Assess whatever assumptions are important to the type of modelling you have decided to implement (use autoplot()). 7. Interpret the model with anova() and summary(). 8. Add fitted lines and CIs, when appropriate, to the plot; make the figure that is worth 1000 words. 9. Keep your data file and script file very safe (and backed up in a separate location). With these, you can recreate all your graphs