Covariation Analysis
library(dplyr)
library(tidyr)
library(ggplot2)
Learning Objectives
- examine how covariates reveal patterns in traits
- plot relationships adjusting for covariates
- determine how covariates change interpretation of relationships
Read Phenotypes
Here we read in data and reduce to species with counts at least 10.
surveys <- read.csv("../data/portal_data_joined.csv")
## Remove missing data.
surveys_complete <- surveys %>%
filter(species_id != "", !is.na(weight)) %>%
filter(!is.na(hindfoot_length), sex != "")
# count records per species
species_counts <- surveys_complete %>%
group_by(species_id) %>%
tally
# get names of the species with counts >= 10
frequent_species <- species_counts %>%
filter(n >= 10) %>%
select(species_id)
# filter out the less-frequent species
surveys_complete <- surveys_complete %>%
filter(species_id %in% frequent_species$species_id)
Now focus on one genus, Dipodomys.
surveys_dip <- surveys_complete %>%
filter(genus == "Dipodomys")
Correlation across species in genus
Here we expand from species to genus and look at correlation, first by showing a scatterplot.
surveys_plot <- ggplot(surveys_dip,
aes(x = weight,
y = hindfoot_length))
surveys_plot +
geom_point(alpha = 1, shape = 1)
Add a straight line does not seem to capture the relationship.
surveys_plot +
geom_point(alpha = 1, shape = 1, col = "gray") +
geom_smooth(method = "lm")
Adding a smooth line does a better job, but is this what we really want?
surveys_plot +
geom_point(alpha = 1, shape = 1, col = "gray") +
geom_smooth()
Relate Species to a Phenotype
The above plot has a disturbing pattern that the data seem to be from more than one group, which might be found with species in the genus. Here we examine one phenotype at a time with species, showing boxplots and densities.
ggplot(surveys_dip,
aes(x = species, y = weight)) +
geom_jitter(col="lightgrey") +
geom_boxplot() +
coord_flip()
Given that variability seems to increase with weight, it often makes sense to use a log transform to put weight on a relative scale. Here we will try log10 in two ways.
ggplot(surveys_dip,
aes(x = species, y = weight)) +
geom_jitter(col="lightgrey") +
geom_boxplot() +
coord_flip() +
scale_y_log10()
ggplot(surveys_dip,
aes(x = species, y = log10(weight))) +
geom_jitter(col="lightgrey") +
geom_boxplot() +
coord_flip()
The plots visually look the same except for the scale on the horizontal axis.
ggplot(surveys_dip,
aes(x = weight, col = species, )) +
geom_density() +
scale_x_log10()
Notice how the shapes of the densities are similar, which is supportive of the use of log10 transformation. Sometimes it is easier to see patterns by using facets.
ggplot(surveys_dip,
aes(x = weight, col = species, )) +
geom_density() +
facet_grid(species~.) +
scale_x_log10()
We can formally compare weight among species using analysis of variance. We first fit using the linear model verb lm, then use anova to get an analysis of variance table. [There is another verb, aov, but most people just use the more flexible lm.]
fit <- lm(log10(weight) ~ species, surveys_dip)
anova(fit)
## Analysis of Variance Table
##
## Response: log10(weight)
## Df Sum Sq Mean Sq F value Pr(>F)
## species 2 327.77 163.885 24617 < 2.2e-16 ***
## Residuals 14537 96.78 0.007
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Paradoxically, more information comes with the summary of fit.
summary(fit)
##
## Call:
## lm(formula = log10(weight) ~ species, data = surveys_dip)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.99095 -0.02642 0.01497 0.05276 0.20862
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.6284793 0.0008273 1968.43 <2e-16 ***
## speciesordii 0.0538433 0.0017523 30.73 <2e-16 ***
## speciesspectabilis 0.4416509 0.0019938 221.51 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08159 on 14537 degrees of freedom
## Multiple R-squared: 0.772, Adjusted R-squared: 0.772
## F-statistic: 2.462e+04 on 2 and 14537 DF, p-value: < 2.2e-16
The broom package has many nice ways to turn output of a test into a data frame for later reuse.
broom::tidy(fit)
## # A tibble: 3 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 1.63 0.000827 1968. 0
## 2 speciesordii 0.0538 0.00175 30.7 6.15e-201
## 3 speciesspectabilis 0.442 0.00199 222. 0
If we are really interested in the means, one way to get them is to remove the intercept, which is done with a 0 or -1 in the formula.
fit0 <- lm(log10(weight) ~ 0 + species, surveys_dip)
broom::tidy(fit0)
## # A tibble: 3 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 speciesmerriami 1.63 0.000827 1968. 0
## 2 speciesordii 1.68 0.00154 1089. 0
## 3 speciesspectabilis 2.07 0.00181 1141. 0
However, it is much easier to use the lsmeans package:
lsmeans::lsmeans(fit, ~ species)
## Error in `loadNamespace()`:
## ! there is no package called 'lsmeans'
Recall that these are means of the log10 of weight. Best practice if you want in original units with indication of precision is to find the confidence interval and back-transform (10^value) those. It makes no sense to back-transform the SE.
There is another package called car (Companion to Applied Regression) that goes along with the Applied Regression book by Fox and Weisberg. This is a powerful package, but takes some getting used to. Here is their analysis of variance function.
car::Anova(fit)
## Error in `loadNamespace()`:
## ! there is no package called 'car'
Challenge
Make boxplots and density plots for hindfoot_length. How do these compare with the weight plots?
Plot two phenotypes adjusting for species
surveys_plot <- ggplot(surveys_dip,
aes(x = weight,
y = hindfoot_length,
col = species))
surveys_plot +
geom_point(alpha = 1, shape = 1) +
scale_color_brewer(type="qualitative", palette = "Dark2")
Add a smooth and straight line
surveys_plot +
geom_point(alpha = 1, shape = 1) +
scale_color_brewer(type="qualitative", palette = "Dark2") +
geom_smooth(se=FALSE, col="black") +
geom_smooth(se=FALSE, col="red", method="lm")
This is not very satisfying. The straight line (which we saw in the previous lesson) misses key aspects of the relationship, and the curved line, while going through the data, does this without explicit regard to species. Now we allow separate lines by species. Since our surveys_plot object already identifies color with species, the change is very simple.
surveys_plot +
geom_point(alpha = 0.2, shape = 1) +
scale_color_brewer(type="qualitative", palette = "Dark2") +
geom_smooth(se=FALSE, method="lm")
It might be easier to see this in separate facets.
surveys_plot +
geom_point(alpha = 1, shape = 1) +
scale_color_brewer(type="qualitative", palette = "Dark2") +
geom_smooth(se=FALSE, method="lm", col = "black") +
facet_wrap(~species) +
theme(legend.position="none")
Adjusting correlation by covariate
The correlation is differs across species:
(cor_species <- surveys_dip %>%
group_by(species) %>%
summarize(cor = cor(weight, hindfoot_length)))
## # A tibble: 3 × 2
## species cor
## <chr> <dbl>
## 1 merriami 0.357
## 2 ordii 0.238
## 3 spectabilis 0.475
Compare these with the overall correlation
with(surveys_dip, cor(sqrt(weight), hindfoot_length))
## [1] 0.9051034
Challenge
Why do the species-specific and overall correlations differ so much?
Adjusting group mean by other phenotype
We can turn this around and ask what is the mean for log2(WBC) by Sex.
(surveys_mean <- surveys_dip %>%
group_by(species) %>%
summarize_at(vars(weight, hindfoot_length), mean))
## # A tibble: 3 × 3
## species weight hindfoot_length
## <chr> <dbl> <dbl>
## 1 merriami 43.1 36.0
## 2 ordii 48.9 35.6
## 3 spectabilis 120. 50.0
We can add these to a plot with straight lines to verify that the regression lines go through the means.
ggplot(surveys_dip, aes(x=weight, y=hindfoot_length, col=species)) +
geom_point(alpha = 1, shape = 1) +
scale_color_brewer(type="qualitative", palette = "Dark2") +
geom_smooth(method="lm", col="black") +
geom_point(data = surveys_mean,
aes(x=weight, y=hindfoot_length),
size = 4, col = "black", shape = 1) +
facet_wrap(~species) +
theme(legend.position="none")
Linear Model Fits with additive covariate
Here we repeat what we had before, but we provide an analysis of variance (ANOVA) summary wtih anova() rather than the usual summary(). We are using linear models and ANOVA without explaining many details. That would be a different lesson.
fit <- lm(hindfoot_length ~ weight, surveys_dip)
anova(fit)
## Analysis of Variance Table
##
## Response: hindfoot_length
## Df Sum Sq Mean Sq F value Pr(>F)
## weight 1 320822 320822 74767 < 2.2e-16 ***
## Residuals 14538 62382 4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Now we add species to the model. Note that this model assumes the correlations within species are the same. Put another way, we assume the lines in plots are parallel.
fit_species <- lm(hindfoot_length ~ species + weight, surveys_dip)
anova(fit_species)
## Analysis of Variance Table
##
## Response: hindfoot_length
## Df Sum Sq Mean Sq F value Pr(>F)
## species 2 346184 173092 77724.6 < 2.2e-16 ***
## weight 1 4648 4648 2087.3 < 2.2e-16 ***
## Residuals 14536 32372 2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We can formally compare these fits to see if we improved the fit.
anova(fit,fit_species, test="F")
## Analysis of Variance Table
##
## Model 1: hindfoot_length ~ weight
## Model 2: hindfoot_length ~ species + weight
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 14538 62382
## 2 14536 32372 2 30010 6737.8 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The p value PR(>F) shows that the fit is dramatically improved.
Linear Model with interacting covariate
Are the slopes the same? Sometimes the relationship between two phenotypes shifts with level of a covarate. For instance, males and females have different correlations, and we want to account for that in our linear model. We do that by replacing the plus (+) with an asterisk (*).
fit_species_int <- lm(hindfoot_length ~ species * weight, surveys_dip)
anova(fit_species_int)
## Analysis of Variance Table
##
## Response: hindfoot_length
## Df Sum Sq Mean Sq F value Pr(>F)
## species 2 346184 173092 78610.485 < 2.2e-16 ***
## weight 1 4648 4648 2111.131 < 2.2e-16 ***
## species:weight 2 369 185 83.838 < 2.2e-16 ***
## Residuals 14534 32002 2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We can compare this to the parallel line model.
anova(fit_species, fit_species_int, test="F")
## Analysis of Variance Table
##
## Model 1: hindfoot_length ~ species + weight
## Model 2: hindfoot_length ~ species * weight
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 14536 32372
## 2 14534 32002 2 369.21 83.838 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Notice that the F statistic is the same here as for the species:weight term in the earlier table. This interaction term, species:weight, is highly significant. When there are significant interactions, it is difficult to interpret the main effects. It is usually best in that situation to report the significant interaction and then evaluate each level of the covariate (each sex in this case) separately. That is, interaction means the relationship is different for different levels of the covariate, and you cannot interpret the relationship averaged over those levels.
Looking back at the plot, we see that for merriami species has a steeper slope.
Relationship of Correlation to Slope
The plots above are a bit confusing for species spectabilis, when you notice that this species has the largest correlation of weight to hindfoot_length, 0.4754936. The reason the slope seems off is that the relationship of slope to correlation is as follows:
slope(hf ~ wt) = cor(wt, hf) * sd(hf) / sd(wt)
Here is the calculation done separately by species.
surveys_dip %>%
group_by(species) %>%
summarize(weight_sd = sd(weight),
hf_length_sd = sd(hindfoot_length),
cor = cor(weight, hindfoot_length)) %>%
ungroup %>%
mutate(slope = cor * hf_length_sd / weight_sd) %>%
select(species,cor,slope,weight_sd,hf_length_sd)
## # A tibble: 3 × 5
## species cor slope weight_sd hf_length_sd
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 merriami 0.357 0.0764 6.82 1.46
## 2 ordii 0.238 0.0501 7.93 1.67
## 3 spectabilis 0.475 0.0424 23.1 2.06
Here are the intercept and slope for species spectabilis.
coef(lm(hindfoot_length ~ weight,
surveys_dip %>%
filter(species == "spectabilis")))
## (Intercept) weight
## 44.89605282 0.04239113
The following is a much more advanced approach to get the slope and intercept using the lm() function with another function called do():
(coef_species <- surveys_dip %>%
group_by(species) %>%
do(data.frame(coef = coef(lm(hindfoot_length ~ weight, data=.)))) %>%
mutate(term = rep(c("intercept", "slope"), length.out = n())) %>%
spread(term, coef))
## # A tibble: 3 × 3
## # Groups: species [3]
## species intercept slope
## <chr> <dbl> <dbl>
## 1 merriami 32.7 0.0764
## 2 ordii 33.1 0.0501
## 3 spectabilis 44.9 0.0424
Challenge
Now that you have the numbers, explain the paradox between correlations and slopes, particularly for spectabilis.
Recap
Correlation is useful to look for relationships, but it is important to consider other factors that might affect relationships. Sometimes a single value is adequate, but in this case the correlation between weight and hindfoot_length was confounded by species, and appeared stronger than it really is. That is, we were seeing correlation across species rather than within species.