Tutorial 6.3: Linear regression on the FEV dataset
As an exercise on linear regression, we will analyse the FEV dataset.
1 The FEV dataset
The FEV, which is an acronym for forced expiratory volume, is a measure of how much air a person can exhale (in liters) during a forced breath. In this dataset, the FEV of 606 children, between the ages of 6 and 17, were measured. The dataset also provides additional information on these children: their age, their height, their gender and, most importantly, whether the child is a smoker or a non-smoker.
The goal of this experiment was to find out whether or not smoking has an effect on the FEV of children.
Note: to analyse this dataset properly, we will need some relatively advanced modeling techniques. At the end of this week, you will have seen all three required steps to analyse such a dataset! For now, we will limit ourselves to exploring the data.
Load the required libraries
library(tidyverse)2 Import data
fev <- read_tsv("https://raw.githubusercontent.com/GTPB/PSLS20/master/data/fev.txt")## Parsed with column specification:
## cols(
## age = col_double(),
## fev = col_double(),
## height = col_double(),
## gender = col_character(),
## smoking = col_double()
## )head(fev)## # A tibble: 6 x 5
## age fev height gender smoking
## <dbl> <dbl> <dbl> <chr> <dbl>
## 1 9 1.71 57 f 0
## 2 8 1.72 67.5 f 0
## 3 7 1.72 54.5 f 0
## 4 9 1.56 53 m 0
## 5 9 1.90 57 m 0
## 6 8 2.34 61 f 0There are a few things in the formatting of the data that can be improved upon:
- Both the
genderandsmokingcan be transformed to factors. - The
heightvariable is written in inches. Assuming that this audience is mainly Portuguese/Belgian, inches are hard to interpret. Let’s add a new column,height_cm, with the values converted to centimeter using themutatefunction.
fev <- fev %>%
mutate(gender = as.factor(gender)) %>%
mutate(smoking = as.factor(smoking)) %>%
mutate(height_cm = height*2.54)
head(fev)## # A tibble: 6 x 6
## age fev height gender smoking height_cm
## <dbl> <dbl> <dbl> <fct> <fct> <dbl>
## 1 9 1.71 57 f 0 145.
## 2 8 1.72 67.5 f 0 171.
## 3 7 1.72 54.5 f 0 138.
## 4 9 1.56 53 m 0 135.
## 5 9 1.90 57 m 0 145.
## 6 8 2.34 61 f 0 155.3 Data Exploration
Now, let’s make a first explorative boxplot, showing only the FEV for both smoking categories.
fev %>%
ggplot(aes(x=smoking,y=fev,fill=smoking)) +
scale_fill_manual(values=c("dimgrey","firebrick")) +
theme_bw() +
geom_boxplot(outlier.shape=NA) +
geom_jitter(width = 0.2, size=0.1) +
ggtitle("Boxplot of FEV versus smoking") +
ylab("fev (l)") +
xlab("smoking status")
Did you expect these results?
It appears that children that smoke have a higher median FEV than children that do not smoke. Should we change legislations worldwide and make smoking obligatory for children?
Maybe there is something else going on in the data. Now, we will generate a similar plot, but we will stratify the data based on age (age as factor).
fev %>%
ggplot(aes(x=as.factor(age),y=fev,fill=smoking)) +
geom_boxplot(outlier.shape=NA) +
geom_point(width = 0.2, size = 0.1, position = position_jitterdodge()) +
theme_bw() +
scale_fill_manual(values=c("dimgrey","firebrick")) +
ggtitle("Boxplot of FEV versus smoking, stratified on age") +
ylab("fev (l)") +
xlab("smoking status")## Warning: Ignoring unknown parameters: width
This plot seems to already give us a more plausible picture. First, it seems that we do not have any smoking children of ages 6, 7 or 8. Second, when looking at the results per age “category”, it seems no longer the case that smokers have a much higher FEV than non-smokers; for the higher ages, the contrary seems true.
This shows that taking into account confounders (in this case) is crucial! If we simply analyse the dataset based on the smoking status and FEV values only, our inference might be incorrect.
Can we provide an even better visualization of the data, taking into account more useful explanatory variables with respect to the FEV?
fev %>%
ggplot(aes(x=as.factor(age),y=fev,fill=smoking)) +
geom_boxplot(outlier.shape=NA) +
geom_point(width = 0.2, size = 0.1, position = position_jitterdodge()) +
theme_bw() +
scale_fill_manual(values=c("dimgrey","firebrick")) +
ggtitle("Boxplot of FEV versus smoking, stratified on age and gender") +
ylab("fev (l)") +
xlab("smoking status") +
facet_grid(rows = vars(gender))## Warning: Ignoring unknown parameters: width
This plot holds one extra level of information, the gender of the child. Especially for higher ages, the median FEV is higher for males as compared to females.
The only source of information that is lacking is height. To look at the effect of height, we could simply make a scatterplot displaying the FEV in function of a child’s height (in cm). Additionally, we could color the dots based on gender.
fev %>%
ggplot(aes(x=height_cm,y=fev,color=gender)) +
geom_point() +
scale_color_manual(values=c("darkorchid","olivedrab4")) +
theme_bw() +
ggtitle("Boxplot of FEV versus height") +
ylab("fev (l)") +
xlab("height (cm)")
There is a clear relationship between height and FEV. In addition, we see that for the large height values (>175cm), we mainly find male subjects.
4 Analysis
4.1 Assumptions of linear regression
List the assumptions:
- The observations are independent of each other
- Linearity between the response and predictor variable
- The residues of the model must be normally distributed
- Homoscedasticity of the data
4.1.1 Independence
The first assumption is not met. There are patterns of relatedness in the data. Indeed, we expect that fish of the same species have a more similar response to the poison as compared to fish of another species. Equivalently, we expect that fish of the same weight have a more similar response to the poison as compared to fish of different weight.
For now, we will ignore this problem. In a later tutorial, we will see how we can account for this dependence by including the species data and weigth data in the linear model (class on multiple regression). As such, you may assume the assumption is met.
We will check the other assumptions by first fitting the linear model and plotting the output. As such, we will get all the required diagnostic plots.
model <- lm(fev~smoking, data=fev)
## display the diagnostic plots of the model
plot(model)
We have four diagnostic plots:
4.1.2 Linearity with the Residuals vs fitted plot
- predictor of predictions \(\hat\beta_0+\hat\beta_1 x\) on \(X\)-axis
- residuals on \(Y\)-as \[e_i=y_i-\hat{g}(x_i)=y_i-\hat\beta_0-\hat\beta_1\times x_i,\]
If there would be a linear relationship in the data, the residuals are expected to lie on the y=0 line for the entire range of predicted values. Based on the first diagnostic plot, the linearity assumption is met.
4.1.3 Normal Q-Q
- QQ-plot of the residuals \(e_i\).
The residuals of the linear regression model should be normally distributed. Based on the second diagnostic plot, the normality assumption is not met.
4.1.4 Homoscedasticity
- Square-root of the absolute value of standardized residuals in function of the fitted values
To meet the third assumption of linear regression, the variance on the Square-root of the absolute value of standardized residuals must be similar over the entire range of fitted values. The smoother in the plot helps us with looking at this; it should be nicely horizontal over the entire range of fitted values. This is clearly not the case here; based on the third diagnostic plot, the homosccedasticity assumption is met.
4.2 Log transformation
The normality assumption was not met, so we cannot continue with performing a linear regression analysis.
One possibility is to log-transform the data. Indeed, data with heavy right tails can often be normalized by applying this strategy
Generate a new column in the poison dataset, containing the log-transformed survival times.
fev <- fev %>%
mutate(log.fev = log(fev))Refit the (log-)linear model. Generate diagnostic plots. Assess the assumptions.
log.model <- lm(log.fev~smoking,data=fev)
plot(log.model)
Upon log-transformation, all the required assumption (except the one for independence, see above) are met.
Look at the output of the log-linear model:
summary(log.model)##
## Call:
## lm(formula = log.fev ~ smoking, data = fev)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.7701 -0.1958 0.0048 0.1953 0.8338
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.92286 0.01270 72.649 < 2e-16 ***
## smoking1 0.22591 0.04004 5.642 2.58e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2966 on 604 degrees of freedom
## Multiple R-squared: 0.05007, Adjusted R-squared: 0.0485
## F-statistic: 31.84 on 1 and 604 DF, p-value: 2.582e-08Compute the 95% confidence interval on the model parameters:
confint(log.model)## 2.5 % 97.5 %
## (Intercept) 0.8979164 0.9478117
## smoking1 0.1472777 0.30454235 Conclusion
5.1 Interpretation on the log scale
Currently, all the outcomes should be interpreted on the log-scale. Indeed, since we are now modelling the log FEV, we don’t have direct inference on the FEV as such
We may interpret the output as follows:
- The slope:
If the smoking status changes from 0 (non-smoking) to 1 (smoking), the natural logarithm of FEV will increase, on average, by 0.22591 units (95% CI [0.1472777,0.3045423]). This increase is significant on the 5% significance level (p=2.58e-08).
- The intercept:
The average natural logarithm of FEV for children that have smoking status 0 is 0.92286 (95% CI [0.1472777; 0.3045423]). This value is significantly different from 0 on the 5% significance level (< 2e-16).
5.2 Interpretation on the original scale
The interpretation the log-scale is quite difficult: the natural logarithm survival time is not exactly a casual measure. To ease the interpetation, we may backtransform the results to the original scale (time in minutes). This we can do by taking the exponent of the outcomes:
exp(summary(log.model)$coefficients[,"Estimate"])## (Intercept) smoking1
## 2.516487 1.253463and of their confidence intervals:
exp(confint(log.model))## 2.5 % 97.5 %
## (Intercept) 2.454484 2.580058
## smoking1 1.158676 1.356004Now, we can interpret the results in terms of the geometric mean:
- The slope:
If the smoking status changes from 0 (non-smoking) to 1 (smoking), the geometric mean of the FEV will increase, on average, by 1.253463 liters (95% CI [1.158676; 1.356004]). This increase is significant on the 5% significance level (p=2.58e-08).
- The intercept:
The geometric mean of the FEV for children that have smoking status 0 is 2.516487 liters (95% CI [2.454484; 2.580058]). This value is significantly different from 0 on the 5% significance level (< 2e-16).
Again, note that in this dataset the assumption of independence is actually not met. In the tutorial of multiple regression, we will revisit this exercise. More specifically, we will there study the association between smoking and FEV while accounting for differences in the age, height and gender of the children.