Motivating example: Evalaute the association between Age and Glucose level
Let’s suppose our main research question is to determine whether age was associated with glucose level. We will set the glucose level as our dependent variable and age as our independent variable (or predictor of interest). Then, we update our linear regression model’s structural form:
We call this the “expected” value because we are predicting this using a model. All regression models take existing data and attempt to make predictions. However, if your assumptions are violated, then these predictions are erroneous. As George Box once wrote, “All models are wrong, but some are useful.”
where \(E[Glucose_{i} | Age_{i}]\) denotes the expected Glucose level for subject \(i\) given the Age of subject \(i\).
We can also represent this relationship in a directed acyclic graph (DAG) diagram.
I used DAGitty to generate the DAG diagram. DAGitty can be found here.
DAG diagram illustrating the causal relationship between Age and Glucose level.
Visualize the association between Age and Glucose level
Let’s look at how Age is related to Glucose level by plotting their relationship. We’ll do this using ggplot. As Age increases, the Glucose level also increases. There appears to be a positive relationship between Age and Glucose level.
### Plot the association between the subject's age and glucose level
ggplot(diabetes.data, aes(x = Age, y = Glucose)) +
geom_point() +
stat_smooth()
In general, when describing a regression model, \(Y\) is denoted as the outcome and \(X\) is denoted as the predictor of interest.
Constructing the linear regression model
Now, we can construct a linear regression model with Glucose level as the dependent variable and Age as the independent variable (or predictor of interest). We will use the lm() function with Glucose level as the \(Y\) variable and Age as the \(X\) variable. The formula for a regression model in R uses the ~ symbol. For example, if was want to regress Age on Glucose level, we use the notation Glucose ~ Age.
By using the lm() function, we can construct the linear regression model: lm(Glucose ~ Age, data = diabetes.data). We can create an object that will contain this linear model; I called this object linear.model1.
We generate the 95% confidence interval (CI) by using the confint() function.
\(\beta_{1}\) coefficient is in the
linear regression output as Age, which is
0.71642.
Here is how we put all of this together in R:
### Linear regression model (Y = Glucose, X = Age)
linear.model1 <- lm(Glucose ~ Age, data = diabetes.data)
summary(linear.model1)##
## Call:
## lm(formula = Glucose ~ Age, data = diabetes.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -126.453 -20.849 -3.058 18.304 86.159
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 97.08016 3.34095 29.06 < 2e-16 ***
## Age 0.71642 0.09476 7.56 1.15e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30.86 on 766 degrees of freedom
## Multiple R-squared: 0.06944, Adjusted R-squared: 0.06822
## F-statistic: 57.16 on 1 and 766 DF, p-value: 1.15e-13The lm() function does not generate 95% CI, so you will
need to use the confint() function.
### Generate the 95% CI
confint(linear.model1)## 2.5 % 97.5 %
## (Intercept) 90.5216601 103.6386585
## Age 0.5304001 0.9024361Interpret the linear regression output
Make sure to have the library gtsummary loaded to create
tables from the regression output.
We are interested in the coefficients. To make interpreting the output easier, we can create a table to visualize the critical elements.
#### Present the output in a table
model1 <- tbl_regression(linear.model1, intercept = TRUE)
as_gt(model1) %>%
gt::tab_header("Table 2. Linear regression model output (Glucose ~ Age)") %>%
gt::tab_options(table.align='left')| Table 2. Linear regression model output (Glucose ~ Age) | |||
| Characteristic | Beta | 95% CI | p-value |
|---|---|---|---|
| (Intercept) | 97 | 91, 104 | <0.001 |
| Age | 0.72 | 0.53, 0.90 | <0.001 |
| Abbreviation: CI = Confidence Interval | |||
The Intercept denotes the \(Y\) intercept when \(X\) is equal to zero. In this case, it would be where Glucose level would be on the linear plot when Age is equal to zero, which is 97.08-units of Glucose.
When possible, present the coefficient’s point estimate and the 95%
CI. For example, a 1-year increase in Age is associated with a 0.72-unit
increase in Glucose level (95% CI: 0.53, 0.90).
The Age coefficient denotes the change in Glucose level for a one-unit increase in Age. In other words, a 1-year increase in Age is associated with a 0.72-unit increase in Glucose level. Since the 95% CI is between 0.53-units and 0.90-units of Glucose, it does not include zero, so this association is statistically significant. We can also look a the p-value of the Age coefficient to determine whether this is statistically significant (<0.0001). However, it is preferable to present the 95% CI when describing the association between \(X\) and \(Y\).
The Adjusted R squared denotes the amount of data that are explained by the linear regression model. In other words, the current linear regression model explains 6.8% of the data.
Visualize predicted model
We can plot the linear form of the model against the actual data using the ggPredict() function.
ggPredict(linear.model1, digits = 1, show.point = TRUE, se = TRUE, xpos = 0.5)
Here is a version without the scatterplot.
ggPredict(linear.model1, digits = 1, show.point = FALSE, se = TRUE, xpos = 0.5)
Adding a confounder
We generate a new variable called pregnancy.history that
is a dichotomous variable (0 = no history of pregnancy and 1 = history
of pregnancy).
Let’s add to our current model by including a confounder. We have a variable called Pregnancies but this provides the number of past pregnancies. We want to create a new variable that has a dichotomous outcome: History of Pregnancy or No history of pregnancy. To do this, we need to to subset the data and create rules. Anyone with 1 or more pregnancies will be coded as 1. Anyone who does not have a history of pregnancy will be coded as 0.
#### Generate groups based on pregnancy history (Group 0 = 0, Group 1 = 1 or more pregnancies)
diabetes.data$pregnancy.history[diabetes.data$Pregnancies == 0] = 0
diabetes.data$pregnancy.history[diabetes.data$Pregnancies > 0] = 1
table(diabetes.data$pregnancy.history)##
## 0 1
## 111 657We see that there are 657 women who a history of pregnancy and 111 women with no history of pregnancy.
A DAG diagram illustrating the relationship between Pregnancy History as a confounder on the Age to Glucose relationship can be drawn.
DAG diagram illustrating the causal relationship between Age and Glucose level and Pregnancy History as a confounder.
In our linear regression model, we have the following structural form:
where \(E[Glucose_{i} | Age_{i}]\) denotes the expected Glucose level for subject \(i\) given the Age of subject \(i\).
But we can include a confounder pregnancy.history:
where \(E[Glucose_{i} | Age_{i}, PregnancyHistory_{i}]\) denotes the expected Glucose level for subject \(i\) given the Age of subject \(i\) controlling for Pregnancy History of subject \(i\).
All you need to do to add another varialbe in the regression model is
to use the + symbol. For example
lm(Glucose ~ Age + pregnancy.history, data = diabetes.data).
Using the lm() function, we can add pregnancy.history to the linear regression model:
Here is how we put all of this together in R:
### Linear regression model (Y = Glucose, X1 = Age, X2 = Pregnancy History)
linear.model2 <- lm(Glucose ~ Age + pregnancy.history, data = diabetes.data)
summary(linear.model2)##
## Call:
## lm(formula = Glucose ~ Age + pregnancy.history, data = diabetes.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -125.715 -20.546 -2.991 17.316 87.734
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 102.00752 3.95100 25.818 < 2e-16 ***
## Age 0.76050 0.09638 7.891 1.04e-14 ***
## pregnancy.history -7.47264 3.22137 -2.320 0.0206 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30.77 on 765 degrees of freedom
## Multiple R-squared: 0.07594, Adjusted R-squared: 0.07352
## F-statistic: 31.43 on 2 and 765 DF, p-value: 7.592e-14confint(linear.model2)## 2.5 % 97.5 %
## (Intercept) 94.2514334 109.7636012
## Age 0.5712954 0.9497004
## pregnancy.history -13.7964201 -1.1488593We can present the model output into a table.
#### Present the output in a table
model2 <- tbl_regression(linear.model2, intercept = TRUE)
as_gt(model2) %>%
gt::tab_header("Table 3. Linear regression model output with confounder (Glucose ~ Age + Pregnancy History)") %>%
gt::tab_options(table.align='left')| Table 3. Linear regression model output with confounder (Glucose ~ Age + Pregnancy History) | |||
| Characteristic | Beta | 95% CI | p-value |
|---|---|---|---|
| (Intercept) | 102 | 94, 110 | <0.001 |
| Age | 0.76 | 0.57, 0.95 | <0.001 |
| pregnancy.history | -7.5 | -14, -1.1 | 0.021 |
| Abbreviation: CI = Confidence Interval | |||
You can see that the Age coefficient is slightly different from our first model. It is 0.76 with a 95% CI of 0.57, 0.95. Compare this to the previous model’s result, which was 0.72; 95% CI: 0.53, 0.90.
#### Merge the two linear regression model's outputs
model1 <- tbl_regression(linear.model1, intercept = TRUE)
model2 <- tbl_regression(linear.model2, intercept = TRUE)
table1 <- tbl_merge(tbls = list(model1, model2),
tab_spanner = c("**Model 1**", "**Model 2**"))
as_gt(table1) %>%
gt::tab_header("Table 4. Comparison between linear regression models [Model 1 (crude) v. Model 2 (adjusted)]") %>%
gt::tab_options(table.align='left')| Table 4. Comparison between linear regression models [Model 1 (crude) v. Model 2 (adjusted)] | ||||||
| Characteristic |
Model 1
|
Model 2
|
||||
|---|---|---|---|---|---|---|
| Beta | 95% CI | p-value | Beta | 95% CI | p-value | |
| (Intercept) | 97 | 91, 104 | <0.001 | 102 | 94, 110 | <0.001 |
| Age | 0.72 | 0.53, 0.90 | <0.001 | 0.76 | 0.57, 0.95 | <0.001 |
| pregnancy.history | -7.5 | -14, -1.1 | 0.021 | |||
| Abbreviation: CI = Confidence Interval | ||||||
Model 1 is considered the crude model or the unadjusted model. Model 2 is the adjusted model because it is adjusting based on the Pregnancy History confounder. Notice that the \(\beta_{1, unadjusted}\) is 0.72 which is lower than the \(\beta_{1, adjusted}\) result which is 0.76.
Additionally, the Adjusted R squared is higher in model 2 (7.35%) compared to model 1, which was 6.82%. This means that Model 2 does a better job of explaining the data than Model 1.
Let’s plot Model 2’s results.
ggPredict(linear.model2, digits = 1, show.point = FALSE, se = TRUE, xpos = 0.5)
You can derive the difference between the groups with a history of
pregnancy and without a history of pregnancy by substracting 102 - 94.5,
which is 7.5, the \(\beta_{2}\) for
pregnancy.history.
We can see that the group that had a history of pregnancy is lower than the group that did not have a history of pregnancy. This makes sense when you look at the pregnancy.history coefficient. It is -7.5, which means that a subject with a history of pregnancy is associated with a 7.5 decrease in Glucose level (95% CI: -13.80, -1.15) compared to a subject without a history of pregnancy controlling for age. Therefore, for all ranges of Age, the group with a history of pregnancy will have Glucose levels that are 7.5 units lower than a group without a history of pregnancy. You can visualize this on the plot; the linear lines do not cross and remain constant across all ranges of Age. But there is a positive correlation between Age and Glucose level.
Homoscedasticity v. Heteroscedasticity.
