Case 1.2 - Statistics workflow in R continued
In this case study, we will again be looking into the bodyweight of mice, but this time we will include the sex variable. Our hypothesis is that both of our explanatory variables, diet and sex, have an effect on our response, weight, so it will be interesting to analyse these variables jointly.
library(readr) # For reading data.
library(tibble) # The "tibble" datastructure.
library(dplyr) # Working with the tibble datastructure.
library(ggplot2) # Plotting.
library(magrittr) # Need this for the %>% (pipe) operator.
library(devtools) # For printing the session info at the end of the notebook.Read and pre-process data
Assuming we have already downloaded the data in the previous case study, we simply read it.
We simplify the names, as last time, and convert character variables to factors.
# Change column names.
mice <- mice %>% rename('sex' = 'Sex', 'diet' = 'Diet', 'weight' = 'Bodyweight')
# Convert character variables to factors.
mice <- mice %>% mutate(sex=as.factor(sex), diet=as.factor(diet))Exploratory data analysis
First, let’s calculate the mean weight of the mice stratefied by diet and sex. We stratify the observations using the group_by() function, providing the variables we want to group by. The summarize() function makes calculations based on the groups, and the results are stored in a tibble together with the group variables diet and sex variables. We also count the number of samples in each group with the n() function.
mice %>% group_by(diet, sex) %>%
summarize(mean=mean(weight, na.rm=TRUE), samples=n())Let’s plot the distribution of the data, as in the previous case studies, with box plots, histograms, and QQ-plots. This time, we’ve used facet_grid() to group the observations with respect to both diet and sex, and arranged the plots in rows and columns.
ggplot(mice, aes(weight)) +
geom_histogram(na.rm=TRUE, bins=15) +
labs(title='Histogram of weight grouped by diet and sex') +
facet_grid(rows=vars(diet), cols=vars(sex))
We make a box plot as in the previous case study, but use facet_grid() to group in terms of sex as well.
ggplot(mice, aes(x=diet, y=weight)) +
geom_boxplot(na.rm=TRUE) +
labs(x='Diet', y='Weight', title='Box plot of weight grouped by diet and sex') +
facet_grid(~ sex)
Similarly, we make QQ-plots.
ggplot(mice, aes(sample=weight)) +
stat_qq(na.rm=TRUE) + stat_qq_line(na.rm=TRUE) +
labs(title='QQ-plot of weight grouped by diet and sex') +
facet_grid(rows=vars(diet), cols=vars(sex))
On thing that’s very interesting here it that it seems more reasonably now to treat the groups as normally distributed, as the lower tails are not as heavy as before. The box plots in a very illustrative way indicate that both the diet and sex have an effect on weight.
Modelling
We already know, from the previous case study, that diet has an effect on weight, and have strong indications that sex does as well.
We want to fit a linear model between the response \(y\) (weight) and the explanatory variables \(X\) (diet and sex), with intercept \(\alpha\), slope \(\beta\), and assuming normally distributed noise \(\epsilon\).
\[ y = \beta X + \alpha + \epsilon \]
We do this using the lm() function, passing it the data, and a formula relating our variables of interest, and in this case, using the ~ character, we relate to R that we want a linear model of the type we have written above. R will then build the model on the data we pass it.
First, we try a model with both variables but no interactions, using the formula weight ~ diet + sex. We use the summary() function to get some useful information about the model fit.
model1 <- lm(weight ~ diet + sex, data=mice)
summary(model1)
Call:
lm(formula = weight ~ diet + sex, data = mice)
Residuals:
Min 1Q Median 3Q Max
-11.9589 -3.1278 -0.2036 2.5975 23.8264
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.5425 0.2725 86.401 <2e-16 ***
diethf 3.1211 0.3225 9.679 <2e-16 ***
sexM 7.7753 0.3218 24.163 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.666 on 838 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.4459, Adjusted R-squared: 0.4446
F-statistic: 337.2 on 2 and 838 DF, p-value: < 2.2e-16At a glance, this model looks good. Let’s try to fit an interaction between the two variables. Note that the expression weight ~ diet * sex is equivalent to weight ~ diet + sex + diet:sex, where diet:sex is the interaction term.
model2 <- lm(weight ~ diet * sex, data=mice)
summary(model2)
Call:
lm(formula = weight ~ diet * sex, data = mice)
Residuals:
Min 1Q Median 3Q Max
-12.3679 -3.1089 -0.2834 2.5766 24.2211
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.8934 0.3102 77.024 < 2e-16 ***
diethf 2.3755 0.4522 5.253 1.9e-07 ***
sexM 7.0704 0.4397 16.081 < 2e-16 ***
diethf:sexM 1.5086 0.6432 2.345 0.0192 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.653 on 837 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.4495, Adjusted R-squared: 0.4476
F-statistic: 227.9 on 3 and 837 DF, p-value: < 2.2e-16The diethf and sexM estimates are significant at 0.1%, and although the diethf:sexM term has a much lower p-value than the others, it is still significant at 5%.
Note that lm() ignored observations with missing data.
The estimates for the coefficients should be understood the following way. The diethf estimate indicates that the high-fat group is 2.3755 units of measurement heavier than the chow group on average. The sexM estimate says that the males are 7.0504 units of measurement heavier than the females. The interaction term, diethf:sexM indicates that the mice that are both male and fed a high-fat diet are an additional 1.5086 units of measurement heavier.
Let’s check the 95% confidence interval for these measurements. We do this using the confint() function, providing the model, a list of estimates, and the desired percentile. It seems it is uncertain whether the interaction term has any meaningful contribution to the mice weights, as it has a wide confidence interval; on the other hand this confidence interval does not include 0, indicating that the effect is likely non-zero.
confint(model2, c('diethf', 'sexM', 'diethf:sexM'), 0.95) 2.5 % 97.5 %
diethf 1.4879360 3.263098
sexM 6.2074251 7.933443
diethf:sexM 0.2460406 2.771157The last thing we want to do before we wrap this case study up is do a little bit of model inspection by studying the residuals. We assume that the error is normally distributed with zero mean, \(\epsilon \in N(0, \sigma_{\epsilon})\). If this assumption is violated, then the model is a poor fit on the data, regardless of how impressive our p-values are.
Let’s use the summary() function to get some quick summary stats. Turns out the residuals’ mean is very close to zero, but has a high variance.
residuals <- model2$residuals
summary(residuals, digits=22) Min. 1st Qu. Median Mean 3rd Qu. Max.
-12.3679 -3.1089 -0.2834 0.0000 2.5766 24.2211 Because we can’t see the significant digits in the summary() output above, let’s calculate the mean again.
mean(residuals)[1] -1.133758e-16From the box plot and QQ-plot below, it seems that the distribution of the residuals has a long upper tail, but seems fairly normally distributed.
temp <- tibble(residuals=residuals)
ggplot(temp, aes(x='', residuals)) +
geom_boxplot(na.rm=TRUE) +
labs(x='', y='Residuals', title='Box plot of residuals')
temp <- tibble(residuals=residuals)
ggplot(temp, aes(sample=residuals)) +
stat_qq(na.rm=TRUE) + stat_qq_line(na.rm=TRUE) +
labs(x='theoretical', y='sample', title='QQ-plot')
We have fit a linear model to the mice weights with the diet and sex variables, including an interaction term. We have shown that assumptions about normal distributed data holds. We have inspected the resulting models through p-values and confidence intervals of coefficient estimates, as well as the distribution of the residuals. We can safely conclude that both diet and gender has an effect on the weight of this mouse population. We also know that the interaction between these variable has an effect, but are not sure about the true effect size.
devtools::session_info()─ Session info ──────────────────────────────────────────────────────────────────────────────────────────────────────
setting value
version R version 3.4.4 (2018-03-15)
os Ubuntu 16.04.6 LTS
system x86_64, linux-gnu
ui RStudio
language en_US
collate en_US.UTF-8
ctype en_US.UTF-8
tz Atlantic/Faroe
date 2019-08-23
─ Packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────
package * version date lib source
assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.4.4)
backports 1.1.4 2019-04-10 [1] CRAN (R 3.4.4)
base64enc 0.1-3 2015-07-28 [1] CRAN (R 3.4.4)
callr 3.2.0 2019-03-15 [1] CRAN (R 3.4.4)
cli 1.1.0 2019-03-19 [1] CRAN (R 3.4.4)
codetools 0.2-15 2016-10-05 [4] CRAN (R 3.3.1)
colorspace 1.4-1 2019-03-18 [1] CRAN (R 3.4.4)
crayon 1.3.4 2017-09-16 [1] CRAN (R 3.4.4)
desc 1.2.0 2018-05-01 [1] CRAN (R 3.4.4)
devtools * 2.0.2 2019-04-08 [1] CRAN (R 3.4.4)
digest 0.6.19 2019-05-20 [1] CRAN (R 3.4.4)
dplyr * 0.8.1 2019-05-14 [1] CRAN (R 3.4.4)
evaluate 0.14 2019-05-28 [1] CRAN (R 3.4.4)
fansi 0.4.0 2018-10-05 [1] CRAN (R 3.4.4)
foreach * 1.4.4 2017-12-12 [1] CRAN (R 3.4.4)
fs 1.3.1 2019-05-06 [1] CRAN (R 3.4.4)
ggplot2 * 3.1.1 2019-04-07 [1] CRAN (R 3.4.4)
glue 1.3.1 2019-03-12 [1] CRAN (R 3.4.4)
gtable 0.3.0 2019-03-25 [1] CRAN (R 3.4.4)
hms 0.4.2 2018-03-10 [1] CRAN (R 3.4.4)
htmltools 0.3.6 2017-04-28 [1] CRAN (R 3.4.4)
iterators * 1.0.10 2018-07-13 [1] CRAN (R 3.4.4)
itertools * 0.1-3 2014-03-12 [1] CRAN (R 3.4.4)
jsonlite 1.6 2018-12-07 [1] CRAN (R 3.4.4)
knitr 1.23 2019-05-18 [1] CRAN (R 3.4.4)
labeling 0.3 2014-08-23 [1] CRAN (R 3.4.4)
lazyeval 0.2.2 2019-03-15 [1] CRAN (R 3.4.4)
lubridate * 1.7.4 2018-04-11 [1] CRAN (R 3.4.4)
magrittr * 1.5 2014-11-22 [1] CRAN (R 3.4.4)
memoise 1.1.0 2017-04-21 [1] CRAN (R 3.4.4)
missForest * 1.4 2013-12-31 [1] CRAN (R 3.4.4)
munsell 0.5.0 2018-06-12 [1] CRAN (R 3.4.4)
pillar 1.4.1 2019-05-28 [1] CRAN (R 3.4.4)
pkgbuild 1.0.3 2019-03-20 [1] CRAN (R 3.4.4)
pkgconfig 2.0.2 2018-08-16 [1] CRAN (R 3.4.4)
pkgload 1.0.2 2018-10-29 [1] CRAN (R 3.4.4)
plyr 1.8.4 2016-06-08 [1] CRAN (R 3.4.4)
prettyunits 1.0.2 2015-07-13 [1] CRAN (R 3.4.4)
processx 3.3.1 2019-05-08 [1] CRAN (R 3.4.4)
ps 1.3.0 2018-12-21 [1] CRAN (R 3.4.4)
purrr 0.3.2 2019-03-15 [1] CRAN (R 3.4.4)
R6 2.4.0 2019-02-14 [1] CRAN (R 3.4.4)
randomForest * 4.6-14 2018-03-25 [1] CRAN (R 3.4.4)
Rcpp 1.0.1 2019-03-17 [1] CRAN (R 3.4.4)
readr * 1.3.1 2018-12-21 [1] CRAN (R 3.4.4)
remotes 2.0.4 2019-04-10 [1] CRAN (R 3.4.4)
reshape2 1.4.3 2017-12-11 [1] CRAN (R 3.4.4)
rlang 0.3.4 2019-04-07 [1] CRAN (R 3.4.4)
rmarkdown 1.13 2019-05-22 [1] CRAN (R 3.4.4)
rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.4.4)
rstudioapi 0.10 2019-03-19 [1] CRAN (R 3.4.4)
scales 1.0.0 2018-08-09 [1] CRAN (R 3.4.4)
sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.4.4)
stringi 1.4.3 2019-03-12 [1] CRAN (R 3.4.4)
stringr 1.4.0 2019-02-10 [1] CRAN (R 3.4.4)
testthat 2.1.1 2019-04-23 [1] CRAN (R 3.4.4)
tibble * 2.1.2 2019-05-29 [1] CRAN (R 3.4.4)
tidyselect 0.2.5 2018-10-11 [1] CRAN (R 3.4.4)
usethis * 1.5.0 2019-04-07 [1] CRAN (R 3.4.4)
utf8 1.1.4 2018-05-24 [1] CRAN (R 3.4.4)
vctrs 0.1.0 2018-11-29 [1] CRAN (R 3.4.4)
withr 2.1.2 2018-03-15 [1] CRAN (R 3.4.4)
xfun 0.7 2019-05-14 [1] CRAN (R 3.4.4)
yaml 2.2.0 2018-07-25 [1] CRAN (R 3.4.4)
zeallot 0.1.0 2018-01-28 [1] CRAN (R 3.4.4)
[1] /home/olavur/R/x86_64-pc-linux-gnu-library/3.4
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/library