library(ggplot2)
library(tibble)
library(broom)
library(readr)
library(skimr)
library(dplyr)
library(RwR)3 Birth weight—a case study
The question that we address in this case study is how birth weight of children is associated with a number of other observable variables. The dataset comes from a sub-study of The Danish National Birth Cohort Study. The Danish National Birth Cohort Study was a nationwide study of pregnant women and their offspring (Olsen et al. 2001). The included pregnant women completed a computer assisted telephone interview, scheduled to take place in pregnancy weeks 12–16 (for some, the interview took place later). We consider data on women interviewed before pregnancy week 17, who were still pregnant in week 17. One of the original purposes of the sub-study was to investigate if fever episodes during pregnancy were associated with fetal death.
Olsen, Jørn, Mads Melbye, Sjurdur F. Olsen, et al. 2001. “The Danish National Birth Cohort - Its Background, Structure and Aim.” Scandinavian Journal of Public Health 29 (4): 300–307.
We focus on birth weight as the outcome variable of interest. If \(Y\) denotes the birth weight of a child, the objective is to find a good model of \(Y\) given a relevant set of predictor variables \(X\). What we believe to be relevant can depend upon many things, for instance, that the variables used as predictors should be observable when we want to make a prediction. If we want to answer causal questions, e.g., does coffee drinking during pregnancy affect the birth weight of the child, we need to pay attention to assumptions required to answer such a question. Importantly, in a regression modeling context we should then also condition on confounders (often called adjusting in the literature).
In this chapter, we focus particularly on the relation between birth weight and the gestational age (age of fetus) at birth, and, given gestational age, how much the additional predictors add to the model.
3.1 Outline and prerequisites
The chapter primarily relies on the following R packages.
The chapter covers an exploratory analysis of the birth weight data included in the RwR package and then a number of univariate and multivariate linear regression models.
3.2 Exploratory data analysis
First, we obtain the dataset by reading it from the csv file included in the RwR package (Requires RwR version >= 0.3.1). The data is also available in the package as a dataset, but it is instructive to see how it is imported from a file as well.
Code
birth_weight_file <- system.file("extdata/birth_weight.txt", package = "RwR")
birth_weight <- read_table(
birth_weight_file,
col_types = cols(
interviewWeek = col_character(),
fetalDeath = col_factor(levels = c("0", "1")),
abortions = col_factor(levels = c("0", "1", "2", "3")),
children = col_factor(levels = c("0", "1")),
smoking = col_factor(levels = c("1", "2", "3")),
coffee = col_factor(levels = c("1", "2", "3")),
feverEpisodes = col_integer())
)Mistakes are easily made if the classes of the columns in the data frame are not appropriate.
The read_table() function infers the column types automatically if they are not specified. For this dataset, several of the categorical variables are encoded using integers, and the column classes therefore require explicit specification above.
It is always a good idea to check that the data was read correctly, that the columns of the resulting data frame have the correct names and are of the correct classes, and to check the dimensions of the resulting data frame using dim(), nrow() or ncol(). This dataset has 12 variables and 11817 cases (rows).
| interviewWeek | fetalDeath | age | abortions | children | gestationalAge | smoking | alcohol | coffee | length | weight | feverEpisodes |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 14 | 0 | 36.73374 | 0 | 1 | 40 | 1 | 0.0 | 1 | NA | NA | 0 |
| 12 | 0 | 34.98700 | 0 | 1 | 41 | 3 | 2.0 | 2 | 53 | 3900 | 2 |
| 13 | 1 | 33.70021 | 0 | 0 | 35 | 1 | 0.0 | 1 | NA | NA | 0 |
| 16 | 0 | 33.06229 | 0 | 1 | 38 | 1 | 4.0 | 2 | 48 | 2800 | 0 |
| 15 | 0 | 31.83299 | 0 | 0 | 40 | 1 | 1.0 | 1 | 54 | 4040 | 0 |
| 13 | 0 | 33.02943 | 0 | 1 | 42 | 1 | 0.5 | 1 | 53 | 3600 | 1 |
NA, partly due to fetal death.
The 12 variables and their encoding in the birth weight dataset:
interviewWeek: Pregnancy week at interview.fetalDeath: Indicator of fetal death (1 = death).age: Mother’s age at conception in years.abortions: Number of previous spontaneous abortions (0, 1, 2, 3+).children: Indicator of previous children (1 = previous children).gestationalAge: Gestational age in weeks at end of pregnancy.smoking: Smoking status; 0, 1–10 or 11+ cigs/day encoded as 1, 2 or 3.alcohol: Number of weekly drinks during pregnancy.coffee: Coffee consumption; 0, 1–7 or 8+ cups/day encoded as 1, 2 or 3.length: Birth length in cm.weight: Birth weight in gram.feverEpisodes: Number of mother’s fever episodes before interview.
3.2.1 Descriptive summaries
The first step is to summarize the variables in the dataset using simple descriptive statistics. This is to get an idea about the data and the variable ranges, but also to discover potential issues that we need to take into account in the further analysis. The list of issues we should be aware of includes, but is not limited to,
- extreme observations and potential outliers
- missing values
- skewness or asymmetry of marginal distributions
Anything worth noticing should be noticed. It should not necessarily be written down in a final report, but figures and tables should be prepared to reveal and not conceal.
A quick summary of the variables in a data frame can be obtained with the summary() function. It prints quantile information for numeric variables and frequencies for factor variables. This is the first example where the class of the columns matter for the result that R produces. Information on the number of missing observations for each variable is also given. The skim() function from the skimr package provides a more detailed view on the data with a better formatting, and it clearly breaks down the summary according to the variable type.
Code
skim(birth_weight)Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| interviewWeek | 0 | 1 | 1 | 2 | 0 | 10 | 0 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| fetalDeath | 39 | 1 | FALSE | 2 | 0: 11659, 1: 119 |
| abortions | 0 | 1 | FALSE | 4 | 0: 9598, 1: 1709, 2: 395, 3: 115 |
| children | 0 | 1 | FALSE | 2 | 1: 6513, 0: 5304 |
| smoking | 0 | 1 | FALSE | 3 | 1: 8673, 2: 1767, 3: 1377 |
| coffee | 4 | 1 | FALSE | 3 | 1: 7821, 2: 3624, 3: 368 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| age | 0 | 1.00 | 29.63 | 4.20 | 16.29 | 26.64 | 29.48 | 32.48 | 44.91 | ▁▆▇▃▁ |
| gestationalAge | 0 | 1.00 | 39.44 | 2.82 | 17.00 | 39.00 | 40.00 | 41.00 | 47.00 | ▁▁▁▇▁ |
| alcohol | 1 | 1.00 | 0.51 | 0.92 | 0.00 | 0.00 | 0.00 | 1.00 | 15.00 | ▇▁▁▁▁ |
| length | 538 | 0.95 | 51.78 | 5.89 | 0.00 | 51.00 | 52.00 | 54.00 | 99.00 | ▁▁▇▁▁ |
| weight | 538 | 0.95 | 3572.45 | 593.45 | 0.00 | 3250.00 | 3600.00 | 3950.00 | 6140.00 | ▁▁▇▆▁ |
| feverEpisodes | 0 | 1.00 | 0.20 | 0.48 | 0.00 | 0.00 | 0.00 | 0.00 | 10.00 | ▇▁▁▁▁ |
skim() function from the skimr package.
Further investigations of the marginal distributions of the variables in the dataset can be obtained by using histograms, density estimates, tabulations and barplots. Barplots are preferable over histograms for numerical variables that take only a small number of different values, e.g., counts. This is the case for the feverEpisodes variable. Before such figures and tables are produced—or perhaps after they have been produced once, but before they enter a final report—we may prefer to clean the data a little.
We can observe from the data summary in Table 3.1 that for some cases weight or length is missing or registered as \(0\)—and in some other cases weight or length is unrealistically small. These are most likely registration mistakes. Likewise, some lengths are registered as 99, and further scrutiny, e.g., by a scatter plot, reveals an observation with weight 3550 gram but gestationalAge registered as 18. We clean the data by excluding the cases with unrealistic values for some variables. This will leave four cases with missing values of coffee and one with a missing value of alcohol. We exclude those five cases as well, and we also remove the columns interviewWeek and fetalDeath as they will not be used subsequently.
Code
birth_weight <- filter(
birth_weight,
weight > 32,
length > 10 & length < 99,
gestationalAge > 18
) |>
na.omit() |>
select(- c(interviewWeek, fetalDeath))We present density or histogram estimates of five numerical variables, see Figure 3.1. The density estimates, as the majority of the figures presented in this book, were produced using the ggplot2 package. Readers familiar with ordinary R graphics can easily produce histograms with the hist() function or density estimates with the density() function. For this simple task ggplot() function does not offer much of an advantage—besides the fact that figures have the same style as other ggplot figures. However, the well-thought-out design and entire functionality of ggplot has resulted in plotting methods that are powerful and expressive. The benefit is that with ggplot it is possible to produce quite complicated figures with clear and logical R expressions—and without the need to mess with low-level technical plotting details.
What is most noteworthy in Figure 3.1 is that the distribution of alcohol is extremely skewed, with more than half of the cases not drinking alcohol at all. This is noteworthy since little variation in a predictor makes it more difficult to detect whether it is associated with the outcome.
Code
tmp <- lapply(
names(birth_weight),
\(x) ggplot(data = birth_weight[, x, drop = FALSE], aes(.data[[x]])) +
xlab(x) + ylab("")
)
gd <- geom_density(adjust = 2, fill = gray(0.7))
gb <- geom_bar(fill = gray(0.7))
gridExtra::grid.arrange(
tmp[[1]] + gd,
tmp[[4]] + gb,
tmp[[6]] + gb,
tmp[[8]] + gd,
tmp[[9]] + gd,
ncol = 5
)
For the categorical or discrete count variables we produce barplots instead of density estimates. Figure 3.2 shows that all these variables except children have quite skewed distributions.
In summary, the typical pregnant woman does not smoke or drink alcohol or coffee, nor has she had any previous spontaneous abortions or any fever episodes. About one-third drinks coffee or alcohol or smokes. These observations may not be surprising—they reflect what is to be expected for a random sample of cases. Little variation of a predictor can, however, make estimation and detection of associations between the outcome and the predictors more difficult. In this case the dataset is quite large, and cases with the least frequently occurring predictor values are present in reasonable numbers.
Code
gridExtra::grid.arrange(
tmp[[2]] + gb,
tmp[[3]] + gb,
tmp[[5]] + gb,
tmp[[7]] + gb,
tmp[[10]] + aes(x = factor(feverEpisodes)) + gb + xlab("feverEpisodes"),
ncol = 5
)
# The coercion of 'feverEpisodes' to a factor fixes the x-axes.
3.2.2 Pairwise associations
The next step is to investigate associations of the variables. We are still not attempting to build a predictive model and the outcome does not yet play a special role. One purpose is again to get better acquainted with the data—this time by focusing on covariation—but there is also one particular issue that we should be aware of.
- Near collinearity of predictors.
Add this bullet point to the previous list of issues. If two or more predictors in the dataset are strongly correlated, they contain, from a predictive point of view, more or less the same information, but perhaps encoded in slightly different ways. Strongly correlated predictors result in the same problem as predictors with little variation. It can become difficult to estimate and detect joint association of the predictors with the outcome.
A technical consequence of near collinearity is that statistical tests of whether one of the predictors could be excluded become non-significant if the other is included, whereas a test of joint exclusion of the predictors can be highly significant. Thus it will become difficult to determine, on statistical grounds, if one predictor should be included over the other. It is best to know about such potential issues upfront. Perhaps it is, by subject matter arguments, possible to choose one of the predictors over the other as the most relevant to include in the model.
A scatter plot matrix is a useful graphical summary of the pairwise association of numerical variables. It can be supplemented with computations of Pearson correlations.
Code
ggally_hexbin <- function(data, mapping, ...)
ggplot(data, mapping) + stat_binhex(...)
GGally::ggpairs(
birth_weight,
columns = c("age", "alcohol", "gestationalAge", "length", "weight"),
lower = list(continuous = GGally::wrap("hexbin", bins = 20)),
diag = list(continuous = "blankDiag"),
)
The scatter plots, Figure 12.1, show that length and weight are (not surprisingly) very correlated, and that both of these variables are also highly correlated with gestationalAge. The alcohol and age variables are mildly correlated, but they are virtually uncorrelated with the other three variables.
The scatter plot matrix was produced using the ggpairs() function from the GGally package. The dataset is quite large and a naive plot of a scatter plot matrix results in a lot of overplotting and large pdf graphics files. Figures can be saved as high-resolution png files instead of pdf files to remedy problems with file size. The actual plotting may, however, still be slow, and the information content in the plot may be limited due to the overplotting. A good way to deal with overplotting is to use hexagonal binning of data points. This was done using the custom ggally_hexbin() function based on the stat_hexbin().
If two categorical variables are strongly dependent the corresponding vectors of dummy variable encoding of the categorical levels will be collinear. In extreme cases where only certain pairwise combinations of the categorical variables are observed, the resulting dummy variable vectors will be perfectly collinear. Dependence between categorical variables may be investigated by cross-tabulation and formal tests of independence can, for instance, be computed. However, neither test statistics of independence nor corresponding \(p\)-values are measures of the degree of dependence—they scale with the size of the dataset and become more and more extreme for larger datasets. There is no single suitable substitute for the Pearson correlation that applies to categorical variables in general. In this particular example all the categorical variables are, in fact, ordinal. In this case we can use the Spearman correlation. The Spearman correlation is simply the Pearson correlation between the ranks of the observations. Since we only need to be able to sort observations to compute ranks, the Spearman correlation is well defined for ordinal as well as numerical variables.
Code
cp <- cor(
data.matrix(birth_weight),
use = "complete.obs",
method = "spearman"
)
corrplot::corrplot(
cp,
diag = FALSE,
order = "hclust",
addrect = 4,
tl.srt = 45,
tl.col = "black",
tl.cex = 0.8
)
Figure 3.4 shows Spearman correlations of all variables—categorical as well as numerical. For continuous variables the Spearman correlation is invariant to monotone transformations and thus less sensitive to outliers than the Pearson correlation. This invariance and its broader applicability make the Spearman correlation attractive as a means for exploratory investigations of pairwise association.
For the production of the plot of the correlation matrix, Figure 3.4, we used a hierarchical clustering of the variables. The purpose was to sort the variables so that the large correlations are concentrated around the diagonal. Since there is no natural order of the variables, the correlation matrix could be plotted using any order. We want to choose an order that brings highly correlated variables close together to make the figure easier to read. Hierarchical clustering can be useful for this purpose.
What we see most clearly from Figure 3.4 are three groupings of positively correlated variables. The weight, length and gestationalAge group, a group consisting of age, children and abortions (not surprising), and a grouping of alcohol, smoking and coffee with mainly coffee being correlated with the two others.
An alternative way to study the relation between a continuous and a categorical variable is to look at the distribution of the continuous variable stratified according to the values of the categorical variable. This can be done using boxplots.
Figure 3.5 shows boxplots for the continuous variable weight. From this figure we see that weight seems to be larger if the mother has had children before and to be negatively related to coffee drinking and smoking.
Code
birth_weight_long <- tidyr::pivot_longer(
birth_weight[, c("weight", "abortions", "children", "smoking", "coffee", "feverEpisodes")],
cols = !weight,
values_transform = as.factor
)
ggplot(birth_weight_long, aes(x = value, y = weight)) +
geom_boxplot(fill = gray(0.8), outlier.alpha = 0.1) + xlab("") +
facet_wrap(~ name, scale = "free_x", ncol = 5)
3.2.3 Birth weight and gestational age
As mentioned in the introduction, our primary interest is to build a model of the relation between gestational age at birth and the birth weight. The scatter plot of weight against gestational age, which is part of the scatter plot matrix in Figure 3.7, shows the strong marginal association. Figure 3.6 gives an alternative graphical illustration of this association.
Code
p_bar <- ggplot(birth_weight, aes(as.factor(gestationalAge))) +
geom_bar(fill = gray(0.7)) + theme(
axis.title.x = element_blank(),
axis.text.x = element_blank(),
)
p_box <- ggplot(birth_weight, aes(x = as.factor(gestationalAge), y = weight)) +
geom_boxplot(fill = gray(0.8), outlier.alpha = 0.1) + xlab("Gestational age")
gridExtra::grid.arrange(p_bar, p_box, heights = c(1, 4))
Figure 3.6 corroborates what we have already observed. There is a strong relation between birth weight and gestational age. The overall trend is also clear: birth weight increases with fetus age. However, the figure also indicates that this relation is not linear in the gestational age variable, and it even suggests that past week 44, the weight might decrease.
3.3 Linear models of birth weight
This section will build a multivariate linear regression model with weight as the outcome variable. Before building a joint regression model we will explore the marginal associations between the other variables and weight in greater details via univariate regression models with one model for each predictor.
3.3.1 Univariate regression models
We analyze the data here via a number of separate univariate linear regression models. The purpose of this analysis is not to do variable selection or draw any strong conclusions—it serves instead as a supplement to the exploratory analysis in the previous sections. In contrast to correlation considerations, the univariate regression models allow us to study single predictor associations with the outcome for categorical as well as numerical predictors, and the analysis can be generalized to models where the outcome is categorical.
The analysis uses the add1() function to add single predictors from a template formula (form) to avoid tedious repetitions of fitting the models using lm().
Code
form <- weight ~ gestationalAge +
length +
age +
children +
coffee +
alcohol +
smoking +
abortions +
feverEpisodes
null_model <- lm(weight ~ 1, data = birth_weight)
one_term_models <- add1(null_model, form, test = "F")Table 3.3 shows the result of testing if inclusion of each of the predictors by themselves is significant. That is, we test the model with only an intercept against the alternative where a single predictor is included. The test used is the \(F\)-test—see Section 2.6.2 for details on the theory. For each of the categorical predictor variables the encoding requires df (degrees of freedom) dummy variables in addition to the intercept to encode the inclusion of a variable with df + 1 levels.
Code
one_term_models |>
tidy() |>
arrange(p.value, desc(statistic)) |>
filter(term != "<none>") |>
select(-AIC) |>
knitr::kable()| term | df | sumsq | rss | statistic | p.value |
|---|---|---|---|---|---|
| length | 1 | 2351722932.3 | 1262184175 | 2.077487e+04 | 0.0000000 |
| gestationalAge | 1 | 932314502.1 | 2681592605 | 3.876542e+03 | 0.0000000 |
| children | 1 | 97619557.6 | 3516287550 | 3.095475e+02 | 0.0000000 |
| smoking | 2 | 53321111.0 | 3560585996 | 8.348023e+01 | 0.0000000 |
| coffee | 2 | 21990295.9 | 3591916811 | 3.412799e+01 | 0.0000000 |
| abortions | 3 | 6272999.9 | 3607634107 | 6.461428e+00 | 0.0002294 |
| age | 1 | 3954121.7 | 3609952986 | 1.221303e+01 | 0.0004764 |
| feverEpisodes | 1 | 1085988.3 | 3612821119 | 3.351611e+00 | 0.0671660 |
| alcohol | 1 | 169962.2 | 3613737145 | 5.244098e-01 | 0.4689818 |
Figure 3.7 shows the scatter plots of weight against the 4 continuous predictors. This is actually just the last row in the scatter plot matrix in Figure 12.1, but this time we have added the regression line. For the continuous variables the tests reported in Table 3.3 are tests of whether the regression line has slope \(0\).
Code
birth_weight_long <- tidyr::pivot_longer(
select(birth_weight, c(age, alcohol, gestationalAge, length, weight)),
cols = !weight,
names_to = "variable"
)
binScale <- scale_fill_continuous(
breaks = c(1, 10, 100, 1000),
low = "gray80",
high = "black",
trans = "log",
guide = "none"
)
ggplot(birth_weight_long, aes(value, weight)) +
stat_binhex(bins = 20) + binScale + xlab("") +
facet_wrap(~ variable, scales = "free_x", ncol = 4) +
geom_smooth(linewidth = 1, method = "lm")
Figure 3.7 illustrates the univariate linear regression models of weight for each of the four continuous predictor variables together with the scatter plot. It is clear from these plots that birth weight is only weakly, if at all, marginally associated with age of the mother or alcohol consumption, while birth weight is strongly associated with gestational age and birth length.
3.3.2 A multivariate linear regression model
To build a multivariate regression model of the outcome variable weight, we need to decide which predictors we want to include. Our initial analyses showed that length is (unsurprisingly) a very good predictor of weight, but it is also close to being an equivalent “body size” measurement, and it will be affected in similar ways as weight by variables that affect fetus growth. From a predictive modeling point of view it is in most cases useless, as it will not be observable unless weight is also observable.
We are particularly interested in the association of weight with gestationalAge, but one might argue that gestationalAge is likewise of little interest if we want to predict weight early in pregnancy. We are, however, interested in what the other predictors add on top of gestationalAge. The exploratory analysis showed that gestationalAge is virtually unrelated to the other predictors, and the remaining predictors are not strongly correlated, and we have not found reasons to exclude any of them. We will thus fit a main effects linear model with eight predictors. We include all the predictors as they are, that is, we do not transform any of them.
The main effects model.
Code
form <- update(form, . ~ . - length)
birth_weight_lm <- lm(form, data = birth_weight)
tidy(birth_weight_lm) |>
knitr::kable()| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -2169.444992 | 98.598974 | -22.0027137 | 0.0000000 |
| gestationalAge | 145.160059 | 2.303741 | 63.0105736 | 0.0000000 |
| age | -2.002013 | 1.204674 | -1.6618710 | 0.0965668 |
| children1 | 185.954556 | 9.896906 | 18.7891615 | 0.0000000 |
| coffee2 | -65.539288 | 10.390137 | -6.3078367 | 0.0000000 |
| coffee3 | -141.784319 | 27.241898 | -5.2046418 | 0.0000002 |
| alcohol | -2.749242 | 5.088446 | -0.5402910 | 0.5890072 |
| smoking2 | -101.952291 | 13.049484 | -7.8127448 | 0.0000000 |
| smoking3 | -131.188927 | 14.910537 | -8.7984041 | 0.0000000 |
| abortions1 | 27.842979 | 13.088798 | 2.1272373 | 0.0334223 |
| abortions2 | 48.763060 | 25.445584 | 1.9163663 | 0.0553440 |
| abortions3 | -50.029260 | 45.803151 | -1.0922668 | 0.2747395 |
| feverEpisodes | 6.363578 | 9.393500 | 0.6774448 | 0.4981378 |
Table 3.4 shows the estimated \(\beta\)-parameters among other things. Note that all categorical variables (specifically, those that are encoded as factors in the data frame) are included via a dummy variable representation. The precise encoding is determined by a linear constraint, known as a contrast. By default, the first factor level is constrained to have parameter \(0\), in which case the remaining parameters represent differences to this base level. In this case it is only occasionally of interest to look at the \(t\)-tests for testing if a single parameter is \(0\).
Table 3.5 shows instead \(F\)-tests of excluding any one of the predictors. It shows that the predictors basically fall into two groups; the strong predictors gestationalAge, children, smoking and coffee, and the weak predictors abortions, age, feverEpisodes and alcohol. The table was obtained using the drop1() function. We resist the temptation to use the tests for model reduction or model selection, as this will generally invalidate subsequent statistical analyses unless properly accounted for.
Code
drop1(birth_weight_lm, test = "F") |>
tidy() |>
arrange(p.value, desc(statistic)) |>
filter(term != "<none>") |>
select(-AIC) |>
knitr::kable()| term | df | sumsq | rss | statistic | p.value |
|---|---|---|---|---|---|
| gestationalAge | 1 | 904499510.74 | 3442125857 | 3970.3323795 | 0.0000000 |
| children | 1 | 80425963.18 | 2618052310 | 353.0325909 | 0.0000000 |
| smoking | 2 | 26628256.86 | 2564254603 | 58.4428345 | 0.0000000 |
| coffee | 2 | 12916211.03 | 2550542557 | 28.3480811 | 0.0000000 |
| abortions | 3 | 2072100.55 | 2539698447 | 3.0318527 | 0.0280891 |
| age | 1 | 629181.69 | 2538255528 | 2.7618151 | 0.0965668 |
| feverEpisodes | 1 | 104551.28 | 2537730898 | 0.4589315 | 0.4981378 |
| alcohol | 1 | 66502.33 | 2537692849 | 0.2919143 | 0.5890072 |
To investigate the model assumptions, we turn to model diagnostics. Several aspects of the statistical analysis presented so far rely on these assumptions, though the theory is postponed to the subsequent sections. Most notably, the distribution of the test statistics, and thus the \(p\)-values, depend on the strong set of assumptions, A3 + A5. We cannot hope to prove that the assumptions are fulfilled, but we can check—mostly using graphical methods—that they are either not obviously wrong, and if they appear to be wrong, the methods should reveal how we can improve the model.
There exists a plot method for lm objects that can produce a range of standard diagnostic plots. The method is useful for quick interactive usage, but the resulting plots are difficult to customize for publication quality. We develop alternatives based on ggplot2.
Model diagnostics for the linear model are mostly based on the residuals, which are estimates of the unobserved errors \(\varepsilon_i\), or the standardized residuals, which are estimates of \(\varepsilon_i / \sigma\). Plots of the standardized residuals against the fitted values, or against any one of the predictors, are useful to detect deviations from A1 or A2. For A3 we consider a qq-plot against the standard normal distribution. The assumptions A4 or A5 are more difficult to investigate. If we do not have a specific idea about how the errors, and thus the observations, might be correlated, it is very difficult to do anything.
Code
birth_weight_diag <- augment(birth_weight_lm)
p1 <- ggplot(birth_weight_diag, aes(.fitted, .std.resid)) +
stat_binhex(bins = 20) + binScale + geom_smooth(linewidth = 1, fill = "blue") +
xlab("fitted values") + ylab("standardized residuals")
p2 <- ggplot(birth_weight_diag, aes(gestationalAge, .std.resid)) +
stat_binhex(bins = 20) + binScale + geom_smooth(linewidth = 1, fill = "blue") +
xlab("gestationalAge") + ylab("")
p3 <- ggplot(birth_weight_diag, aes(sample = .std.resid)) +
geom_abline(intercept = 0, slope = 1, color = "blue", linewidth = 1) +
geom_qq() +
xlab("theoretical quantiles") + ylab("")
gridExtra::grid.arrange(p1, p2, p3, ncol = 3)
gestationalAge, and a qq-plot against the normal distribution.
The residual plot in Figure 3.8 shows that the model is not spot on. The plot of the residuals against gestationalAge shows that there is a nonlinear effect that the linear model does not catch. Thus A1 is not fulfilled. We address this specific issue in Section 7.3.2, where we solve the problem using splines. The qq-plot shows that the tails of the residuals are heavier than the normal distribution and right skewed. However, given the problems with A1, this issue is of secondary interest.
The diagnostics considered above address if the data set as a whole does not comply to the model assumptions. Single observations can also be extreme and, for instance, have a large influence on how the model is fitted. For this reason we should also be aware of single extreme observations in the residual plots and the qq-plot.
3.3.3 Including interactions
One standard way to try to improve on the multivariate linear regression model from the previous section is to include interactions between the predictors. The inclusion of interactions results in a substantial increase in the complexity of the models, even if we have only a few predictors. Moreover, it becomes possible to construct an overwhelming number of comparisons of models. Searching haphazardly through thousands of models with various combinations of interactions is not recommended. It will result in spurious discoveries that will be difficult to reproduce in other studies. The strategy we follow here is to include interactions of the strongest predictors from the main effects model and investigate if such interactions add any value to the model.
To comprehend an interaction model it is helpful to visualize the model to the extent it is possible. This is a point where the ggplot2 package is really strong. It supports a number of ways to stratify a plot according to different variables.
Code
ggplot(birth_weight, aes(gestationalAge, weight)) +
facet_grid(coffee ~ children + smoking, label = label_both) +
binScale + stat_binhex(bins = 20) +
geom_smooth(method = "lm", linewidth = 1, se = FALSE, fullrange = TRUE)
weight against gestationalAge stratified according to the values of smoking, children and coffee
Figure 3.9 shows a total of 18 scatter plots where the stratification is according to children, smoking and coffee. A regression line was fitted separately for each plot. This corresponds to a model with a third order interaction between the four strong predictors (and with the weak predictors left out). Variations between the regression lines are seen across the different plots, which is an indication of interaction effects. For better comparison of the regression lines it can be beneficial to plot them differently. Figure 3.10 shows an example where the stratification according to coffee is visualized by color coding the levels of coffee.
Code
ggplot(birth_weight, aes(gestationalAge, weight, color = coffee)) +
facet_grid(. ~ children + smoking, label = label_both) +
geom_smooth(method = "lm", linewidth = 1, se = FALSE, fullrange = TRUE)
gestationalAge stratified according to the values of smoking, coffee and children
Below we specify a linear regression model via the formula interface that includes a third order interaction among the strong predictors. The parameters are estimated via lm(), but we will not print out all of the 42 parameters. Instead we compare the model to the main effects model with no interactions via an \(F\)-test, which is done with the anova() function.
Code
form <- weight ~ smoking * coffee * children * gestationalAge +
age +
alcohol +
abortions +
feverEpisodes
birth_weight_lm_int <- lm(form, data = birth_weight)
anova(
birth_weight_lm,
birth_weight_lm_int
) |> knitr::kable()| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 11139 | 2537626346 | NA | NA | NA | NA |
| 11110 | 2514285831 | 29 | 23340516 | 3.55641 | 0 |
Table 3.6 shows that the \(F\)-test of the third order interaction model against the main effects model is clearly significant, that is, we can reject the null hypothesis of no interactions against the interaction model.
We reconsider model diagnostics for the extended model, where we have included the interactions. Figure 3.11 shows the residual plot. The inclusion of the interactions did not solve the earlier observed problems with the model fit. This is hardly surprising as the problem with the model appears to be related to a nonlinear relation between weight and gestationalAge. An apparent nonlinearity could be explained by interaction effects, but this would require a strong correlation between some predictors, e.g., that heavy coffee drinkers (coffee = 3) have large values of gestationalAge. We already established that this was not the case.
Since there is some lack of model fit, we should be skeptical about the conclusions from formal hypothesis tests. However, deviations from A1 result in an increased residual variance, which will generally result in more conservative tests. That is, it will become harder to reject a null hypothesis, and thus, in this case, conclude that inclusion of the interactions is significant.
Before we conclude the analysis, we test the hypothesis where we exclude the four weak predictors together. Table 3.7 shows that the test results in a borderline \(p\)-value of around 5%. On the basis of this we might exclude the four weak predictors even though Table 3.5 suggested that the number of abortions is related to weight. The highly skewed distribution of abortions resulted in large standard errors and low power despite the size of the data set. In combination with the different signs on the estimated parameters in Table 3.4, depending upon whether the woman had had 1, 2 or 3+ spontaneous abortions, the study is inconclusive on how abortions is related to weight.
Code
form <- weight ~ smoking * coffee * children * gestationalAge
birth_weight_lm_int_sub <- lm(form, data = birth_weight)
anova(
birth_weight_lm_int_sub,
birth_weight_lm_int
) |> knitr::kable()| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 11116 | 2517158539 | NA | NA | NA | NA |
| 11110 | 2514285831 | 6 | 2872708 | 2.11563 | 0.0482525 |
In conclusion, we have developed a predictive model of weight given in terms of a third order interaction of the four predictors gestationalAge, smoking, coffee and children and additive linear effects of the four weak predictors age, alcohol, abortions and feverEpisodes. The model is not a perfect fit, as it does not catch a nonlinear relation between weight and gestationalAge. The fitted model can be visualized as in the Figure 3.9 or Figure 3.10.
We note that the formal \(F\)-test of the interaction model against the main effects model justifies the need for the increased model complexity. It is, however, clear from the figures that the actual differences in slopes are small, and the significance of the test reflects that we have a large data set, which makes us able to detect even small differences.
There is no clear-cut interpretation of the interactions either. The regression lines in the figures should, preferably, be equipped with confidence bands. This can be achieved by removing the se = FALSE argument to the geom_smooth() function, but this will result in a separate variance estimate for each combination of smoking, coffee and children. If we want to use the pooled variance estimate obtained by our model, we have to do something else. How this is achieved is shown in Section 7.3.2, where we also consider how to deal with the nonlinearity using spline basis expansions.
Exercises
Exercise 3.1 This exercise is on the modeling of annual wage. The data set wage from the R package RwR is from a US survey with 1125 respondents. It contains five variables. The outcome variable is wage, which is the respondent’s annual wage in USD. In addition to wage the dataset includes the variables: hours (annual hours worked), educ (education, 1 = no high school diploma, 2 = high school graduate, 3 = college graduate), age (age of respondent in years) and children (number of children living with the respondent). The respondents with hours and wage equal to \(0\) are unemployed.
The objective of this exercise is to fit a predictive regression model of wage using some or all of the other variables in the dataset.
There are few respondents with many children, and it might therefore be a good idea to collapse respondents with three or more children, say, into one group.
Fit and report a linear, additive model of the mean wage given some of the other variables using data for the employed respondents only.
Investigate how well the linear, additive model fits the data.
Investigate the effect of nonlinear variable transformations and inclusion of interaction effects. Choose a model to report.
Use the model chosen to predict the wage distribution for the unemployed respondents should they get a job. Compare with the wage distribution for the employed respondents. Discuss reasons that would invalidate such wage predictions for the unemployed.