library(haven)
library(broom)Data analysis for epidemiology
PUST Campus
Introduction
The aim of this analysis is to investigate whether smoking is associated with reduced lung function (FEV) in children using linear regression models. Smoking status is considered the primary exposure, while age, height and sex are evaluated as potential confounders or effect modifiers.
Load library
Load dataset
lung <- read_dta("lung_function_data.dta")
lung# A tibble: 655 × 5
age fev height sex smoke
<dbl> <dbl> <dbl> <dbl+lbl> <dbl+lbl>
1 9 1.71 57 1 [Female] 0 [No]
2 8 1.72 67.5 1 [Female] 0 [No]
3 7 1.72 54.5 1 [Female] 0 [No]
4 9 1.56 53 0 [Male] 0 [No]
5 9 1.89 57 0 [Male] 0 [No]
6 8 2.34 61 1 [Female] 0 [No]
7 6 1.92 58 1 [Female] 0 [No]
8 6 1.41 56 1 [Female] 0 [No]
9 8 1.99 58.5 1 [Female] 0 [No]
10 9 1.94 60 1 [Female] 0 [No]
# ℹ 645 more rowsConvert categorical variables
lung$smoke <- factor(lung$smoke)
lung$sex <- factor(lung$sex)Question 1.1: Uni-variate Linear Regression
The crude association between each explanatory variable and FEV is assessed using separate linear regression models.
Smoking Status (Crude Association)
summary(lm(fev ~ smoke, data = lung))
Call:
lm(formula = fev ~ smoke, data = lung)
Residuals:
Min 1Q Median 3Q Max
-1.7751 -0.6339 -0.1021 0.4804 3.2269
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.56614 0.03466 74.037 < 2e-16 ***
smoke1 0.71072 0.10994 6.464 1.99e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8412 on 652 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.06023, Adjusted R-squared: 0.05879
F-statistic: 41.79 on 1 and 652 DF, p-value: 1.993e-10With Age
summary(lm(fev ~ age, data = lung))
Call:
lm(formula = fev ~ age, data = lung)
Residuals:
Min 1Q Median 3Q Max
-1.57539 -0.34567 -0.04989 0.32124 2.12786
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.431648 0.077895 5.541 4.36e-08 ***
age 0.222041 0.007518 29.533 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5675 on 652 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.5722, Adjusted R-squared: 0.5716
F-statistic: 872.2 on 1 and 652 DF, p-value: < 2.2e-16With Height
summary(lm(fev ~ height, data = lung))
Call:
lm(formula = fev ~ height, data = lung)
Residuals:
Min 1Q Median 3Q Max
-1.75167 -0.26619 -0.00401 0.24474 2.11936
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.432679 0.181460 -29.94 <2e-16 ***
height 0.131976 0.002955 44.66 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4307 on 652 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.7537, Adjusted R-squared: 0.7533
F-statistic: 1995 on 1 and 652 DF, p-value: < 2.2e-16With Sex
summary(lm(fev ~ sex, data = lung))
Call:
lm(formula = fev ~ sex, data = lung)
Residuals:
Min 1Q Median 3Q Max
-2.01615 -0.69366 -0.06417 0.58083 2.98085
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.8121 0.0462 60.876 < 2e-16 ***
sex1 -0.3610 0.0663 -5.445 7.35e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.848 on 653 degrees of freedom
Multiple R-squared: 0.04343, Adjusted R-squared: 0.04196
F-statistic: 29.65 on 1 and 653 DF, p-value: 7.35e-08Question 1.2 An exploration of confounding and effect modification
Step 1. Crude Model
model1 <- lm(fev ~ smoke, data = lung) # Fit the crude model with smoking only.
summary(model1)
Call:
lm(formula = fev ~ smoke, data = lung)
Residuals:
Min 1Q Median 3Q Max
-1.7751 -0.6339 -0.1021 0.4804 3.2269
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.56614 0.03466 74.037 < 2e-16 ***
smoke1 0.71072 0.10994 6.464 1.99e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8412 on 652 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.06023, Adjusted R-squared: 0.05879
F-statistic: 41.79 on 1 and 652 DF, p-value: 1.993e-10Step 2. Assess Confounding
Determine Confounding
A variable was considered a potential confounder if its inclusion in the regression model changed the estimated regression coefficient for smoking by 10% or more compared with the crude model (the change-in-estimate criterion). Variables meeting this criterion were retained in the multivariable model as confounders.
Smoking + Age
model2 <- lm(fev ~ smoke + age, data = lung)
summary(model2)
Call:
lm(formula = fev ~ smoke + age, data = lung)
Residuals:
Min 1Q Median 3Q Max
-1.6653 -0.3564 -0.0508 0.3494 2.0894
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.367373 0.081436 4.511 7.65e-06 ***
smoke1 -0.208995 0.080745 -2.588 0.00986 **
age 0.230605 0.008184 28.176 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5651 on 651 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.5766, Adjusted R-squared: 0.5753
F-statistic: 443.3 on 2 and 651 DF, p-value: < 2.2e-16Smoking + Height
model3 <- lm(fev ~ smoke + height, data = lung)
summary(model3)
Call:
lm(formula = fev ~ smoke + height, data = lung)
Residuals:
Min 1Q Median 3Q Max
-1.7505 -0.2660 -0.0041 0.2447 2.1207
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.427620 0.187577 -28.935 <2e-16 ***
smoke1 0.006319 0.058686 0.108 0.914
height 0.131883 0.003081 42.808 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.431 on 651 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.7537, Adjusted R-squared: 0.7529
F-statistic: 995.9 on 2 and 651 DF, p-value: < 2.2e-16Smoking + Sex
model4 <- lm(fev ~ smoke + sex, data = lung)
summary(model4)
Call:
lm(formula = fev ~ smoke + sex, data = lung)
Residuals:
Min 1Q Median 3Q Max
-1.95758 -0.63880 -0.03123 0.52159 3.03942
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.75358 0.04541 60.643 < 2e-16 ***
smoke1 0.76070 0.10726 7.092 3.45e-12 ***
sex1 -0.39571 0.06420 -6.163 1.25e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8183 on 651 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.112, Adjusted R-squared: 0.1093
F-statistic: 41.07 on 2 and 651 DF, p-value: < 2.2e-16Calculate the percentage change
beta_crude <- coef(model1)["smoke1"]
beta_age <- coef(model2)["smoke1"]
abs(beta_age - beta_crude) *100/ beta_crude smoke1
129.4061 Since the regression coefficient for smoking changed by more than 10% after adjusting for age, age was identified as a potential confounder.
beta_Height <- coef(model3)["smoke1"]
abs(beta_Height - beta_crude) *100/ beta_crude smoke1
99.11085 Adjustment for height changed the estimated regression coefficient for smoking by 99.1% compared with the crude model. Since this change exceeded the 10% threshold, height was identified as a potential confounder.
beta_sex <- coef(model4)["smoke1"]
abs(beta_sex - beta_crude) *100/ beta_crude smoke1
7.032875 Since the change in the smoking regression coefficient was less than 10% after adjusting for sex, sex was not considered a confounder according to the change-in-estimate criterion.
The next step is to examine whether these variables were also associated with the primary exposure (smoking status), which is an important characteristic of a confounding variable.
model5 <- glm(smoke ~ age, data = lung, family = binomial)
summary(model5)
Call:
glm(formula = smoke ~ age, family = binomial, data = lung)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -7.74391 0.70890 -10.924 <2e-16 ***
age 0.48364 0.05513 8.773 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 423.45 on 653 degrees of freedom
Residual deviance: 318.56 on 652 degrees of freedom
(1 observation deleted due to missingness)
AIC: 322.56
Number of Fisher Scoring iterations: 6model6 <- glm(smoke ~ height, data = lung, family = binomial)
summary(model6)
Call:
glm(formula = smoke ~ height, family = binomial, data = lung)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -15.3262 2.0171 -7.598 3.01e-14 ***
height 0.2067 0.0308 6.712 1.92e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 423.45 on 653 degrees of freedom
Residual deviance: 365.17 on 652 degrees of freedom
(1 observation deleted due to missingness)
AIC: 369.17
Number of Fisher Scoring iterations: 6The results showed that both age and height were significantly associated with the primary exposure (smoking status), supporting their potential role as confounding variables. Therefor we need to retained age and height in the subsequent multivariable analyses.
Are age, height, or sex effect modifiers of the association between smoking and FEV?
Look at the interaction term:
smoke:age
smoke:height
smoke:sex
If the p-value < 0.05, there is evidence of effect modification (interaction). If the p-value is ≥ 0.05, there is no evidence that the effect of smoking differs across levels of that variable.
Smoking × Age
model_age_int <- lm(fev ~ smoke * age, data = lung)
summary(model_age_int)
Call:
lm(formula = fev ~ smoke * age, data = lung)
Residuals:
Min 1Q Median 3Q Max
-1.76645 -0.34947 -0.03364 0.33679 2.05990
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.253396 0.082651 3.066 0.00226 **
smoke1 1.943571 0.414285 4.691 3.31e-06 ***
age 0.242558 0.008332 29.113 < 2e-16 ***
smoke1:age -0.162703 0.030738 -5.293 1.65e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5537 on 650 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.5941, Adjusted R-squared: 0.5922
F-statistic: 317.1 on 3 and 650 DF, p-value: < 2.2e-16The interaction term between smoking and age was statistically significant (p < 0.001), indicating that age is a potential effect modifier of the association between smoking and FEV.
Smoking × Height
model_height_int <- lm(fev ~ smoke * height, data = lung)
summary(model_height_int)
Call:
lm(formula = fev ~ smoke * height, data = lung)
Residuals:
Min 1Q Median 3Q Max
-1.74332 -0.26960 -0.00462 0.23909 2.12941
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.364866 0.190421 -28.174 <2e-16 ***
smoke1 -2.036782 1.128464 -1.805 0.0716 .
height 0.130847 0.003128 41.832 <2e-16 ***
smoke1:height 0.031062 0.017133 1.813 0.0703 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4302 on 650 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.7549, Adjusted R-squared: 0.7538
F-statistic: 667.3 on 3 and 650 DF, p-value: < 2.2e-16The interaction term between smoking and height was not statistically significant (p > 0.05), indicating no evidence of effect modification by age.
Smoking × Sex
model_sex_int <- lm(fev ~ smoke * sex, data = lung)
summary(model_sex_int)
Call:
lm(formula = fev ~ smoke * sex, data = lung)
Residuals:
Min 1Q Median 3Q Max
-2.04923 -0.64038 -0.02221 0.50700 3.05862
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.73438 0.04638 58.959 < 2e-16 ***
smoke1 1.00885 0.16672 6.051 2.43e-09 ***
sex1 -0.35517 0.06739 -5.271 1.85e-07 ***
smoke1:sex1 -0.42211 0.21745 -1.941 0.0527 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8166 on 650 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.1172, Adjusted R-squared: 0.1131
F-statistic: 28.75 on 3 and 650 DF, p-value: < 2.2e-16The interaction term between smoking and sex was not statistically significant (p = 0.053). Therefore, there was no statistically significant evidence that sex modified the association between smoking and FEV in the unadjusted interaction model.
Question: 1.3 Build an appropriate final multivariable linear regression model, estimate the effect of being a smoker on FEV level, controlling for confounders and effect modifiers and answer the question:
is there any evidence that smoking leads to worse lung function in children?
is there any evidence that the effect of smoking on lung function is modified by sex of the child?
Recommended multivariable model building approach: forward step-wise approach to adjust for confounders, where confounders are identified in the second part of the analysis. # For question ii), you will need to fit and test an interaction between sex and smoking status in the regression model, then extract the appropriate stratum-specific coefficients
model_final1 <- lm(fev ~ smoke+age+smoke*age, data = lung)
summary(model_final1)
Call:
lm(formula = fev ~ smoke + age + smoke * age, data = lung)
Residuals:
Min 1Q Median 3Q Max
-1.76645 -0.34947 -0.03364 0.33679 2.05990
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.253396 0.082651 3.066 0.00226 **
smoke1 1.943571 0.414285 4.691 3.31e-06 ***
age 0.242558 0.008332 29.113 < 2e-16 ***
smoke1:age -0.162703 0.030738 -5.293 1.65e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5537 on 650 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.5941, Adjusted R-squared: 0.5922
F-statistic: 317.1 on 3 and 650 DF, p-value: < 2.2e-16The interaction term between smoking and age remained statistically significant (p < 0.001) after adjusting for age, indicating that the association between smoking and FEV varied according to age.
Adding height in the model_final1
model_final2 <- lm(fev ~ smoke+age+smoke*age+height, data = lung)
summary(model_final2)
Call:
lm(formula = fev ~ smoke + age + smoke * age + height, data = lung)
Residuals:
Min 1Q Median 3Q Max
-1.49605 -0.26445 -0.01816 0.25104 1.98173
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.598716 0.229069 -20.076 < 2e-16 ***
smoke1 0.005672 0.325785 0.017 0.986
age 0.061191 0.010375 5.898 5.92e-09 ***
height 0.108581 0.004931 22.019 < 2e-16 ***
smoke1:age -0.008796 0.024300 -0.362 0.717
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4192 on 649 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.7677, Adjusted R-squared: 0.7662
F-statistic: 536.1 on 4 and 649 DF, p-value: < 2.2e-16After adjusting for height, which had been identified as a potential confounder, the interaction between smoking and age was no longer statistically significant (p = 0.717). This indicates that the apparent effect modification by age observed in the previous model was largely explained by confounding due to height. Therefore, there was no evidence that age modified the association between smoking and FEV after controlling for height, and the interaction term was not retained in the final model.
Final model
model_final <- lm(fev ~ smoke+age+height, data = lung)
summary(model_final)
Call:
lm(formula = fev ~ smoke + age + height, data = lung)
Residuals:
Min 1Q Median 3Q Max
-1.50182 -0.26305 -0.01882 0.24989 1.98535
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.616007 0.223883 -20.618 < 2e-16 ***
smoke1 -0.110232 0.060017 -1.837 0.0667 .
age 0.059741 0.009563 6.247 7.57e-10 ***
height 0.109095 0.004720 23.115 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4189 on 650 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.7676, Adjusted R-squared: 0.7665
F-statistic: 715.7 on 3 and 650 DF, p-value: < 2.2e-16Since the interaction between smoking and age was no longer statistically significant after adjusting for height, the interaction term was removed from the model. Following the principle of model parsimony, the final multivariable linear regression model included smoking, age, and height as the main effects.
The final multivariable linear regression model
Age and height were retained because they were identified as important confounders using the change-in-estimate criterion (>10%) and were significantly associated with both smoking status and FEV. The smoking-by-age interaction was not retained because it was no longer statistically significant after adjustment for height (p = 0.717). Therefore, this model was selected as the final model to estimate the adjusted association between smoking and lung function.
Conclusion
The crude linear regression model indicated a positive association between smoking and FEV (β = 0.711, p < 0.001). After adjustment for age, the smoking coefficient changed to β = −0.209, representing a 129.4% change from the crude estimate. Adjustment for height reduced the smoking coefficient to β = 0.006, representing a 99.1% change. These substantial changes indicate that both age and height are important confounders of the association between smoking and FEV. In contrast, adjustment for sex resulted in only a 7.0% change in the smoking coefficient, indicating that sex is not an important confounder according to the 10% change-in-estimate criterion. Effect modification was assessed by including interaction terms between smoking and each potential modifier. There was strong evidence of an interaction between smoking and age (β = −0.163, p < 0.001), indicating that the effect of smoking on FEV varies with age. There was no statistically significant evidence of interaction between smoking and height (p = 0.070) or between smoking and sex (p = 0.053).