Exploratory data analysis

Exercise 1

Create a table1 and use GGally::ggpairs to explore your data.

# create the table 1
label(nsduh$demog_sex)       <- "Sex"
label(nsduh$alc_agefirst)    <- "Alcohol (age first used)   "
label(nsduh$demog_age_cat6)  <- "Age"
label(nsduh$demog_income)    <- "Income"

table1(~ demog_sex + alc_agefirst + demog_age_cat6 + demog_income | mj_lifetime, 
       data=nsduh, overall = F)
No
(N=491)		 Yes
(N=509)
Sex
Female 285 (58.0%) 249 (48.9%)
Male 206 (42.0%) 260 (51.1%)
Alcohol (age first used)
Mean (SD) 19.7 (5.09) 16.0 (3.08)
Median [Min, Max] 19.0 [3.00, 45.0] 16.0 [3.00, 29.0]
Missing 145 (29.5%) 12 (2.4%)
Age
18-25 53 (10.8%) 80 (15.7%)
26-34 52 (10.6%) 87 (17.1%)
35-49 129 (26.3%) 136 (26.7%)
50-64 109 (22.2%) 132 (25.9%)
65+ 148 (30.1%) 74 (14.5%)
Income
Less than $20,000 76 (15.5%) 85 (16.7%)
$20,000 - $49,999 166 (33.8%) 127 (25.0%)
$50,000 - $74,999 64 (13.0%) 81 (15.9%)
$75,000 or more 185 (37.7%) 216 (42.4%) 
# use GGally::ggpairs to explore your data
ggpairs(
  na.omit(nsduh)[,-2],
  upper = list(continuous = "density", discrete = wrap("colbar", size = 2, alignPercent = 1), 
               combo = "facethist"),
  title = "The National Survey on Drug Use and Health (NSDUH)",
  columnLabels = c("Marijuana use", "Alcohol (age first used)", "Age", "Sex", "Income"),
  aes(colour=mj_lifetime)) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1,size = 8, color = "black"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning in geom_text(mapping = mapping_text, stat = "prop", position = position_fill(0.5, : Ignoring unknown parameters: `alignPercent`
## Ignoring unknown parameters: `alignPercent`
## Ignoring unknown parameters: `alignPercent`
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning in geom_text(mapping = mapping_text, stat = "prop", position = position_fill(0.5, : Ignoring unknown parameters: `alignPercent`
## Ignoring unknown parameters: `alignPercent`
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning in geom_text(mapping = mapping_text, stat = "prop", position =
## position_fill(0.5, : Ignoring unknown parameters: `alignPercent`
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Calculating odds ratio

Exercise 2

Use the contingency table below to calculate (a) probability of lifetime marijuana among sample subjects (independent of sex), among sample females and among sample males, (b) Odds of lifetime marijuana use among sample subjects, among sample females and among sample males. (c) What is the risk difference, the risk ratio and the odds ratio when comparing lifetime marijuana use among males and females in our sample?

# showing absolute numbers
nsduh |> 
  tabyl(demog_sex, mj_lifetime) |> 
  adorn_title()
##            mj_lifetime    
##  demog_sex          No Yes
##     Female         285 249
##       Male         206 260
# showing percentages
nsduh |> 
  tabyl(demog_sex, mj_lifetime) |> 
  adorn_percentages() |> 
  adorn_pct_formatting() |> 
  adorn_title()
##            mj_lifetime      
##  demog_sex          No   Yes
##     Female       53.4% 46.6%
##       Male       44.2% 55.8%
# Replace the blanks below with appropriate values
tribble(
  ~Prediction, ~General, ~Female, ~Male, 
  "Probability", 0.509, 0.4662921, 0.5579399,  # 0.509 = 509/(509+491); similar methods for each sex
  "Odds", 1.03666, 0.8736841, 1.262136, # 1.03666 = 0.509/(1-0.509); similar methods for each sex
  "Log Odds", 0.03600401, -0.1350364, 0.2328055 # 0.03600401 = log(1.03666); similar methods for each sex
) |> gt()
Prediction General Female Male
Probability 0.50900000 0.4662921 0.5579399
Odds 1.03666000 0.8736841 1.2621360
Log Odds 0.03600401 -0.1350364 0.2328055
 
# Replace the blanks below with appropriate values
tribble(
  ~Contrast, ~`Males vs. Females`, 
  "Risk difference", 0.0916478,# 0.5579399-0.4662921
  "Risk Ratio",      1.196546, # 0.5579399/0.4662921
  "Odds ratio",      1.444614, # 1.262136 / 0.8736841
  "Log odds ratio",  0.3678422 # log(1.444614)
) |> gt()
Contrast Males vs. Females
Risk difference 0.0916478
Risk Ratio 1.1965460
Odds ratio 1.4446140
Log odds ratio 0.3678422

Exercise 3

You will now estimate the same measures from the previous question, but this time using two types of regressions: a logistic regression and a regular least squares linear regression. Show your code and present your results in a table adding the confidence intervals (see below). Explain your findings: what are the differences between the three approaches (using a contingency table as in the question above, using logistic regression and using linear regression)? What explains these differences (or lack thereof)? What are the advantages/drawbacks of using one approach over the other?

####################
# The null model 
####################

# Using logistic regression
mdl.glm.null <- glm(
  mj_lifetime ~ 1, 
  data = nsduh, 
  family = binomial
  )

summary(mdl.glm.null)
## 
## Call:
## glm(formula = mj_lifetime ~ 1, family = binomial, data = nsduh)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.03600    0.06326   0.569    0.569
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1386  on 999  degrees of freedom
## Residual deviance: 1386  on 999  degrees of freedom
## AIC: 1388
## 
## Number of Fisher Scoring iterations: 3
# Log - odds for lifetime marijuana (everyone)
tidy(
  mdl.glm.null, 
  conf.int = TRUE
  )
## # A tibble: 1 × 7
##   term        estimate std.error statistic p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 (Intercept)   0.0360    0.0633     0.569   0.569  -0.0880     0.160
# Odds of lifetime marijuana (everyone)
tidy(
  mdl.glm.null, 
  exponentiate = TRUE,
  conf.int = TRUE
  )
## # A tibble: 1 × 7
##   term        estimate std.error statistic p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 (Intercept)     1.04    0.0633     0.569   0.569    0.916      1.17
# Probability of lifetime marijuana (everyone)
tidy(
  mdl.glm.null, 
  conf.int = TRUE
  ) |> 
  select(estimate, conf.low, conf.high) |> 
  mutate_all(plogis)
## # A tibble: 1 × 3
##   estimate conf.low conf.high
##      <dbl>    <dbl>     <dbl>
## 1    0.509    0.478     0.540
####################
# The unadjusted model
####################

# The unadjusted logistic regression model
mdl.glm <- glm(mj_lifetime.numeric ~ demog_sex, data = nsduh, family = binomial)
summary(mdl.glm)
## 
## Call:
## glm(formula = mj_lifetime.numeric ~ demog_sex, family = binomial, 
##     data = nsduh)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   -0.13504    0.08675  -1.557  0.11954   
## demog_sexMale  0.36784    0.12738   2.888  0.00388 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1386.0  on 999  degrees of freedom
## Residual deviance: 1377.6  on 998  degrees of freedom
## AIC: 1381.6
## 
## Number of Fisher Scoring iterations: 3
# Log - odds for lifetime marijuana
tidy(
  mdl.glm, 
  conf.int = TRUE
  )
## # A tibble: 2 × 7
##   term          estimate std.error statistic p.value conf.low conf.high
##   <chr>            <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 (Intercept)     -0.135    0.0867     -1.56 0.120     -0.305    0.0348
## 2 demog_sexMale    0.368    0.127       2.89 0.00388    0.119    0.618
# Odds of lifetime marijuana
tidy(
  mdl.glm, 
  exponentiate = TRUE,
  conf.int = TRUE
  )
## # A tibble: 2 × 7
##   term          estimate std.error statistic p.value conf.low conf.high
##   <chr>            <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 (Intercept)      0.874    0.0867     -1.56 0.120      0.737      1.04
## 2 demog_sexMale    1.44     0.127       2.89 0.00388    1.13       1.86
# Probability of lifetime marijuana (everyone)
tidy(
  mdl.glm, 
  conf.int = TRUE
  ) |> 
  select(estimate, conf.low, conf.high) |> 
  mutate_all(plogis)
## # A tibble: 2 × 3
##   estimate conf.low conf.high
##      <dbl>    <dbl>     <dbl>
## 1    0.466    0.424     0.509
## 2    0.591    0.530     0.650
# Calculating the confidence intervals
avg_predictions(mdl.glm, variables = "demog_sex")
## 
##  demog_sex Estimate Pr(>|z|)   S 2.5 % 97.5 %
##     Female    0.466   0.1195 3.1 0.424  0.509
##     Male      0.558   0.0126 6.3 0.512  0.602
## 
## Columns: demog_sex, estimate, p.value, s.value, conf.low, conf.high 
## Type:  invlink(link)

Now building the table using the results from the logistic regression model. To calculate the results for female and male, we use the equations derived from the model: log odds female = -0.13504 (demog_sexMale = 0) log odds male = -0.13504 + 0.36784 (demog_sexMale = 1)

tribble(
  ~Prediction, ~General, ~Female, ~Male, 
  "Probability", "0.509 CI [0.478, 0.540]", " CI 0.466 [0.424, 0.509]", " CI 0.558 [0.513, 0.603]",
  "Odds", "1.04 CI [0.916, 1.17]", "0.874 CI [0.737, 1.04]", "1.2623 CI [1.053, 1.518]", 
  "Log Odds", "0.036 CI [0.0880, 0.160]", "-0.135 CI [-0.305, 0.0348]", "0.233 CI [0.051, 0.41]"
) |> gt()
Prediction General Female Male
Probability 0.509 CI [0.478, 0.540] CI 0.466 [0.424, 0.509] CI 0.558 [0.513, 0.603]
Odds 1.04 CI [0.916, 1.17] 0.874 CI [0.737, 1.04] 1.2623 CI [1.053, 1.518]
Log Odds 0.036 CI [0.0880, 0.160] -0.135 CI [-0.305, 0.0348] 0.233 CI [0.051, 0.41]

Now calculating the results using linear regression. Notice: we will be violating assumptions. Oh no! 😱😱😱😱😱

####################
# The null model
####################

# Linear regression
mdl.lm.null <- lm(mj_lifetime.numeric ~ 1, data = nsduh)
summary(mdl.lm.null)
## 
## Call:
## lm(formula = mj_lifetime.numeric ~ 1, data = nsduh)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -0.509 -0.509  0.491  0.491  0.491 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.50900    0.01582   32.18   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5002 on 999 degrees of freedom
# Probability of lifetime marijuana (everyone) using linear regression null
tidy(
  mdl.lm.null, 
  conf.int = TRUE
  )
## # A tibble: 1 × 7
##   term        estimate std.error statistic   p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 (Intercept)    0.509    0.0158      32.2 1.75e-156    0.478     0.540
####################
# The unadjusted model
####################

# unadjusted linear model
mdl.lm <- lm(mj_lifetime.numeric ~ demog_sex, data = nsduh)
summary(mdl.lm)
## 
## Call:
## lm(formula = mj_lifetime.numeric ~ demog_sex, data = nsduh)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5579 -0.4663  0.4421  0.4421  0.5337 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.46629    0.02156  21.623   <2e-16 ***
## demog_sexMale  0.09165    0.03159   2.901   0.0038 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4983 on 998 degrees of freedom
## Multiple R-squared:  0.008363,   Adjusted R-squared:  0.00737 
## F-statistic: 8.417 on 1 and 998 DF,  p-value: 0.003799
# Probability of lifetime marijuana (everyone) using unadjusted linear regression 
tidy(
  mdl.lm, 
  conf.int = TRUE
  )
## # A tibble: 2 × 7
##   term          estimate std.error statistic  p.value conf.low conf.high
##   <chr>            <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)     0.466     0.0216     21.6  2.40e-85   0.424      0.509
## 2 demog_sexMale   0.0916    0.0316      2.90 3.80e- 3   0.0297     0.154
# Calculating the confidence intervals
avg_predictions(mdl.lm, variables = "demog_sex")
## 
##  demog_sex Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
##     Female    0.466     0.0216 21.6   <0.001 342.0 0.424  0.509
##     Male      0.558     0.0231 24.2   <0.001 426.3 0.513  0.603
## 
## Columns: demog_sex, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
## Type:  response

Now building the table using the results from the linear regression model. To calculate the results for female and male, we use the equations derived from the model: prob female = 0.46629 (demog_sexMale = 0) prob male = 0.46629 + 0.09165 (demog_sexMale = 1

tribble(
  ~Prediction, ~General, ~Female, ~Male, 
  "Probability", "0.509 CI [0.478, 0.540]", " CI 0.466 [0.424, 0.509]", " CI 0.558 [0.513, 0.603]",
  "Odds",        "1.04 CI [0.916, 1.17]", "0.874 CI [0.737, 1.04]", "1.2623 CI [1.053, 1.518]", 
  "Log Odds",    "0.036 CI [0.0880, 0.160]", "-0.135 CI [-0.305, 0.0348]", "0.233 CI [0.051, 0.41]"
) |> gt()
Prediction General Female Male
Probability 0.509 CI [0.478, 0.540] CI 0.466 [0.424, 0.509] CI 0.558 [0.513, 0.603]
Odds 1.04 CI [0.916, 1.17] 0.874 CI [0.737, 1.04] 1.2623 CI [1.053, 1.518]
Log Odds 0.036 CI [0.0880, 0.160] -0.135 CI [-0.305, 0.0348] 0.233 CI [0.051, 0.41]

Both models provide the same results. Now calculating the contrasts table using the logistic regression model. The log odds ratio male vs female in the logistic regression model is exp(demog_sexMale), from that we can calculate the rest.

# Calculating odds ratio (exp(demox_sexMale))
mdl.glm |>
  tidy(conf.int = TRUE, exponentiate = TRUE)
## # A tibble: 2 × 7
##   term          estimate std.error statistic p.value conf.low conf.high
##   <chr>            <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 (Intercept)      0.874    0.0867     -1.56 0.120      0.737      1.04
## 2 demog_sexMale    1.44     0.127       2.89 0.00388    1.13       1.86
# Calculating risk ratio and CI
avg_comparisons(
  model = mdl.glm, 
  variables = "demog_sex", # This is the exposure variable
  comparison = c("ratio")
)
## 
##       Term      Contrast Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
##  demog_sex Male / Female      1.2     0.0742 16.1   <0.001 192.0  1.05   1.34
## 
## Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
## Type:  response
# Calculating risk difference and CI
avg_comparisons(
  model = mdl.glm, 
  variables = "demog_sex", # This is the exposure variable
  comparison = "difference"
)
## 
##       Term      Contrast Estimate Std. Error   z Pr(>|z|)   S  2.5 % 97.5 %
##  demog_sex Male - Female   0.0916     0.0315 2.9  0.00367 8.1 0.0298  0.153
## 
## Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
## Type:  response

Contrasts using the logistic regression.

tribble(
  ~Contrast, ~`Males vs. Females`, 
  "Risk difference", "0.0916 CI [0.0298, 0.153]", 
  "Risk Ratio",      "1.2 CI [1.05, 1.34]",
  "Odds ratio",      "1.44 CI [1.13, 1.86]", 
  "Log odds ratio",  "0.368 CI [0.119, 0.618]"
) |> gt()
Contrast Males vs. Females
Risk difference 0.0916 CI [0.0298, 0.153]
Risk Ratio 1.2 CI [1.05, 1.34]
Odds ratio 1.44 CI [1.13, 1.86]
Log odds ratio 0.368 CI [0.119, 0.618]

Contrasts table using linear regression.

avg_comparisons(
  model = mdl.lm, 
  variables = "demog_sex",
  comparison = c("lnor")
)
## 
##       Term                      Contrast Estimate Std. Error    z Pr(>|z|)   S
##  demog_sex ln(odds(Male) / odds(Female))    0.368      0.128 2.88  0.00393 8.0
##  2.5 % 97.5 %
##  0.118  0.618
## 
## Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
## Type:  response
avg_comparisons(
  model = mdl.lm, 
  variables = "demog_sex",
  comparison = "ratio"
)
## 
##       Term      Contrast Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
##  demog_sex Male / Female      1.2     0.0742 16.1   <0.001 191.7  1.05   1.34
## 
## Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
## Type:  response
avg_comparisons(
  model = mdl.lm, 
  variables = "demog_sex", # This is the exposure variable
  comparison = "difference"
)
## 
##       Term      Contrast Estimate Std. Error   z Pr(>|z|)   S  2.5 % 97.5 %
##  demog_sex Male - Female   0.0916     0.0316 2.9  0.00372 8.1 0.0297  0.154
## 
## Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
## Type:  response
tribble(
  ~Contrast, ~`Males vs. Females`, 
  "Risk difference", "0.0916 CI [0.0297, 0.154]", 
  "Risk Ratio",      "1.2 CI [1.05, 1.34]",
  "Odds ratio",      "1.44 CI [1.13, 1.86]", 
  "Log odds ratio",  "0.368 CI [0.118, 0.618]"
) |> gt()
Contrast Males vs. Females
Risk difference 0.0916 CI [0.0297, 0.154]
Risk Ratio 1.2 CI [1.05, 1.34]
Odds ratio 1.44 CI [1.13, 1.86]
Log odds ratio 0.368 CI [0.118, 0.618]

Answer: We can see that all three methods gives us the same (or very similar) results. The main difference in using the contingency table and the models is that with the models we can calculate some statistics (error, p-value, CI) very easily with the models and check if our results are significant. In this case we see no difference in using the linear or logistic model, but using linear regression in this case should not be preferred, due to the nature of the data.

Exercise 4

Write up your results using the template shown below.

Conclusion: Males have significantly higher odds of lifetime marijuana use than females (OR = 1.44 ; 95% CI = 1.13 , 1.86 p = 0.0038). Males have 0.09% higher odds of lifetime marijuana use than females.

For this simple case with a binary outcome and a binary predictor, the odds ratio you calculated should be the same as the odds ratio obtained via the cross-product from a 2 x 2 contingency table. But what if your predictor variable is not binary, but continuous?

Exercise 5

Whereas in the previous questions we explored the relationship between lifetime marijuana use and sex, we will now explore the association between lifetime marijuana use (mj_lifetime) and age at first use of alcohol (alc_agefirst)? Fit a logistic regression and interpret the results, using the template shown below.

# fit a unjusted model
mdl.glm.unadjusted <- glm(
  mj_lifetime ~ alc_agefirst,
  data = nsduh,
  family = binomial()
)

summary(mdl.glm.unadjusted)
## 
## Call:
## glm(formula = mj_lifetime ~ alc_agefirst, family = binomial(), 
##     data = nsduh)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    5.3407     0.4747   11.25   <2e-16 ***
## alc_agefirst  -0.2835     0.0267  -10.62   <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: 1141.45  on 842  degrees of freedom
## Residual deviance:  968.44  on 841  degrees of freedom
##   (157 observations deleted due to missingness)
## AIC: 972.44
## 
## Number of Fisher Scoring iterations: 5
tidy(
  mdl.glm.unadjusted, 
  exponentiate = TRUE,
  conf.int = TRUE
  )
## # A tibble: 2 × 7
##   term         estimate std.error statistic  p.value conf.low conf.high
##   <chr>           <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)   209.       0.475       11.3 2.29e-29   84.7     545.   
## 2 alc_agefirst    0.753    0.0267     -10.6 2.46e-26    0.713     0.792

Answer: Age at first alcohol use is significantly negatively associated with lifetime marijuana use (OR = 0.753 ; 95% CI = 0.713 , 0.792 ; p < 2.46e-26). Individuals who first used alcohol at age 19 years have 24.7 % lower odds of having ever used marijuana than those who first used alcohol at age 18 years. In contrast, the reduction in odds associated with starting drinking alcohol 3 years later is 75.3 % lower odds.

Adjusted logistic regression

The main reason for combining multiple predictors is to address the problem of confounding bias. Confounding bias occurs when a predictor affects both outcome and exposure. Adding a confounder to as a predictor to the model removes this bias. With multiple predictors in a logistic regression model, the resulting odds ratios are called Adjusted Odds Ratio (AOR).

Exercise 6

What is the association between lifetime marijuana use (mj_lifetime) and age at first use of alcohol (alc_agefirst), adjusted for age (demog_age_cat6), sex (demog_sex), and income (demog_income)? Fit a model with and without an interaction term (alc_agefirst:demog_sex) and compute the AOR, its 95% CI, and the p-value that tests the significance of age at first use of alcohol. Report also their AORs of the other predictors, their 95% CIs and p-values. Since there are some categorical predictors with more than two levels, use car::Anova(model, type = 3) function to compute p-values for multiple degrees of freedom. In your answer, fill in the template below.

# fit logistic model
mdl.glm.adjusted <- glm(
  mj_lifetime ~ alc_agefirst + demog_age_cat6 + demog_sex + demog_income,
  data = nsduh,
  family = binomial()
)

summary(mdl.glm.adjusted)
## 
## Call:
## glm(formula = mj_lifetime ~ alc_agefirst + demog_age_cat6 + demog_sex + 
##     demog_income, family = binomial(), data = nsduh)
## 
## Coefficients:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    6.25417    0.59136  10.576  < 2e-16 ***
## alc_agefirst                  -0.27541    0.02756  -9.992  < 2e-16 ***
## demog_age_cat626-34           -0.29619    0.32865  -0.901  0.36746    
## demog_age_cat635-49           -0.80427    0.29656  -2.712  0.00669 ** 
## demog_age_cat650-64           -0.68991    0.29854  -2.311  0.02084 *  
## demog_age_cat665+             -1.27475    0.30429  -4.189  2.8e-05 ***
## demog_sexMale                 -0.06089    0.16181  -0.376  0.70669    
## demog_income$20,000 - $49,999 -0.53088    0.26641  -1.993  0.04629 *  
## demog_income$50,000 - $74,999 -0.07931    0.30488  -0.260  0.79476    
## demog_income$75,000 or more   -0.36120    0.25323  -1.426  0.15376    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1141.45  on 842  degrees of freedom
## Residual deviance:  935.55  on 833  degrees of freedom
##   (157 observations deleted due to missingness)
## AIC: 955.55
## 
## Number of Fisher Scoring iterations: 5
tidy(
  mdl.glm.adjusted, 
  exponentiate = TRUE,
  conf.int = TRUE
  )
## # A tibble: 10 × 7
##    term                 estimate std.error statistic  p.value conf.low conf.high
##    <chr>                   <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
##  1 (Intercept)           520.       0.591     10.6   3.85e-26  169.     1722.   
##  2 alc_agefirst            0.759    0.0276    -9.99  1.65e-23    0.718     0.800
##  3 demog_age_cat626-34     0.744    0.329     -0.901 3.67e- 1    0.387     1.41 
##  4 demog_age_cat635-49     0.447    0.297     -2.71  6.69e- 3    0.247     0.791
##  5 demog_age_cat650-64     0.502    0.299     -2.31  2.08e- 2    0.275     0.891
##  6 demog_age_cat665+       0.279    0.304     -4.19  2.80e- 5    0.152     0.502
##  7 demog_sexMale           0.941    0.162     -0.376 7.07e- 1    0.684     1.29 
##  8 demog_income$20,000…    0.588    0.266     -1.99  4.63e- 2    0.347     0.987
##  9 demog_income$50,000…    0.924    0.305     -0.260 7.95e- 1    0.507     1.68 
## 10 demog_income$75,000…    0.697    0.253     -1.43  1.54e- 1    0.421     1.14
car::Anova(mdl.glm.adjusted, type = 3)
## Analysis of Deviance Table (Type III tests)
## 
## Response: mj_lifetime
##                LR Chisq Df Pr(>Chisq)    
## alc_agefirst    149.610  1  < 2.2e-16 ***
## demog_age_cat6   23.964  4   8.12e-05 ***
## demog_sex         0.142  1     0.7066    
## demog_income      5.510  3     0.1380    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# fit a model with interaction
mdl.glm.adjusted.interaction <- glm(
  mj_lifetime ~ alc_agefirst + demog_age_cat6 + demog_sex + demog_income + 
    alc_agefirst:demog_sex,
  data = nsduh,
  family = binomial()
)

summary(mdl.glm.adjusted.interaction)
## 
## Call:
## glm(formula = mj_lifetime ~ alc_agefirst + demog_age_cat6 + demog_sex + 
##     demog_income + alc_agefirst:demog_sex, family = binomial(), 
##     data = nsduh)
## 
## Coefficients:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    7.37354    0.82921   8.892  < 2e-16 ***
## alc_agefirst                  -0.33782    0.04229  -7.988 1.38e-15 ***
## demog_age_cat626-34           -0.29506    0.33145  -0.890   0.3733    
## demog_age_cat635-49           -0.81603    0.29870  -2.732   0.0063 ** 
## demog_age_cat650-64           -0.68657    0.30071  -2.283   0.0224 *  
## demog_age_cat665+             -1.26004    0.30636  -4.113 3.91e-05 ***
## demog_sexMale                 -2.13813    0.98473  -2.171   0.0299 *  
## demog_income$20,000 - $49,999 -0.53035    0.26822  -1.977   0.0480 *  
## demog_income$50,000 - $74,999 -0.09306    0.30727  -0.303   0.7620    
## demog_income$75,000 or more   -0.36288    0.25466  -1.425   0.1542    
## alc_agefirst:demog_sexMale     0.11893    0.05552   2.142   0.0322 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1141.5  on 842  degrees of freedom
## Residual deviance:  930.9  on 832  degrees of freedom
##   (157 observations deleted due to missingness)
## AIC: 952.9
## 
## Number of Fisher Scoring iterations: 5
tidy(
  mdl.glm.adjusted.interaction, 
  exponentiate = TRUE,
  conf.int = TRUE
  )
## # A tibble: 11 × 7
##    term                 estimate std.error statistic  p.value conf.low conf.high
##    <chr>                   <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
##  1 (Intercept)          1593.       0.829      8.89  5.99e-19 337.      8717.   
##  2 alc_agefirst            0.713    0.0423    -7.99  1.38e-15   0.654      0.772
##  3 demog_age_cat626-34     0.744    0.331     -0.890 3.73e- 1   0.385      1.42 
##  4 demog_age_cat635-49     0.442    0.299     -2.73  6.30e- 3   0.243      0.785
##  5 demog_age_cat650-64     0.503    0.301     -2.28  2.24e- 2   0.275      0.897
##  6 demog_age_cat665+       0.284    0.306     -4.11  3.91e- 5   0.153      0.511
##  7 demog_sexMale           0.118    0.985     -2.17  2.99e- 2   0.0167     0.800
##  8 demog_income$20,000…    0.588    0.268     -1.98  4.80e- 2   0.346      0.991
##  9 demog_income$50,000…    0.911    0.307     -0.303 7.62e- 1   0.498      1.66 
## 10 demog_income$75,000…    0.696    0.255     -1.42  1.54e- 1   0.419      1.14 
## 11 alc_agefirst:demog_…    1.13     0.0555     2.14  3.22e- 2   1.01       1.26
lrtest(mdl.glm.adjusted, mdl.glm.adjusted.interaction)
## Likelihood ratio test
## 
## Model 1: mj_lifetime ~ alc_agefirst + demog_age_cat6 + demog_sex + demog_income
## Model 2: mj_lifetime ~ alc_agefirst + demog_age_cat6 + demog_sex + demog_income + 
##     alc_agefirst:demog_sex
##   #Df  LogLik Df  Chisq Pr(>Chisq)  
## 1  10 -467.78                       
## 2  11 -465.45  1 4.6556    0.03095 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpreting the output:

The AOR for our primary predictor alc_agefirst is 0.759. This represents the OR for lifetime marijuana use comparing those with a one-year difference in age at first use of alcohol, adjusted for age, sex, and income.

The remaining AORs compare levels of categorical predictors to their reference level, adjusted for the other predictors in the model. For example, comparing individuals with the same age of first alcohol use, sex, and income, 35-49 year-olds have 55.3 % lower odds of lifetime marijuana use than 18-25 year-olds (OR = 0.447 ; 95% CI = 0.247 , 0.791 ; p = 6.69e-3 ). An overall 4 df p-value for age can be read from the Type III Test table. It is significant at p = 8.12e-05. An overall, 3 df p-value for income is not significant at (p = 0.1380).

Conclusion: After adjusting for age, sex, and income, age at first alcohol use is significantly negatively associated with lifetime marijuana use (AOR = 0.759; 95% CI = 0.718 , 0.800 ; p < 3.85e-26). Individuals who first used alcohol at a given age have 24.1 % lower odds of having ever used marijuana than those who first used alcohol one year earlier. The association between age of first alcohol use and lifetime marijuana use differs significantly between males and females (p = 0.0322). A likelihood ratio test shows that the adjusted model with the interaction is superior to the one without the interaction (p = 0.03095).

Exercise 7

Replicate the model summary shown below, and discuss its main features. Use the likelihood ratio test to determine which model is best, and discuss the association between age of first alcohol and lifetime marijuana use, when comparing males to females.

# summarize model
modelsummary(
  list(
    "unadjusted model" = mdl.glm.unadjusted,
    "adjusted model" = mdl.glm.adjusted,
    "interaction model" = mdl.glm.adjusted.interaction
  ),
  stars = T
)
unadjusted model adjusted model interaction model
(Intercept) 5.341*** 6.254*** 7.374***
(0.475) (0.591) (0.829)
alc_agefirst −0.284*** −0.275*** −0.338***
(0.027) (0.028) (0.042)
demog_age_cat626-34 −0.296 −0.295
(0.329) (0.331)
demog_age_cat635-49 −0.804** −0.816**
(0.297) (0.299)
demog_age_cat650-64 −0.690* −0.687*
(0.299) (0.301)
demog_age_cat665+ −1.275*** −1.260***
(0.304) (0.306)
demog_sexMale −0.061 −2.138*
(0.162) (0.985)
demog_income$20,000 - $49,999 −0.531* −0.530*
(0.266) (0.268)
demog_income$50,000 - $74,999 −0.079 −0.093
(0.305) (0.307)
demog_income$75,000 or more −0.361 −0.363
(0.253) (0.255)
alc_agefirst × demog_sexMale 0.119*
(0.056)
Num.Obs. 843 843 843
AIC 972.4 955.6 952.9
BIC 981.9 1002.9 1005.0
Log.Lik. −484.219 −467.776 −465.448
F 112.743 14.804 13.139
RMSE 0.44 0.43 0.43
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 
# LR test
anova(mdl.glm.unadjusted, mdl.glm.adjusted, mdl.glm.adjusted.interaction, test = "LRT")
## Analysis of Deviance Table
## 
## Model 1: mj_lifetime ~ alc_agefirst
## Model 2: mj_lifetime ~ alc_agefirst + demog_age_cat6 + demog_sex + demog_income
## Model 3: mj_lifetime ~ alc_agefirst + demog_age_cat6 + demog_sex + demog_income + 
##     alc_agefirst:demog_sex
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       841     968.44                          
## 2       833     935.55  8   32.886 6.457e-05 ***
## 3       832     930.90  1    4.656   0.03095 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Answer: The LRT indicated that the interaction model outperformed the other two models. The significant difference observed in the interaction term suggests that younger individuals who first consume alcohol are at a higher risk of using marijuana in the future. Furthermore, females face even greater risk than males when they are young.

Exercise 8

Replicate the forest plot shown below, using centralized variables scaled such that standard deviations are 0.5 (you can use the function datawizard::standardize(two_sd = TRUE). What is the advantage of this plot compared to the model summary tables shown above?

nsduh %>% 
  datawizard::standardize(two_sd = TRUE) -> nsduh.std

mdl.glm.adjusted.interaction.std <- glm(
  mj_lifetime ~ alc_agefirst + demog_age_cat6 + demog_sex + demog_income + 
    alc_agefirst:demog_sex,
  data = nsduh.std,
  family = binomial()
)

mdl.glm.adjusted.interaction.std %>% 
  model_parameters() %>% 
  plot() + 
  labs(x = "Log odds ratio (standardized)")

Answer: The plot clearly shows that the effect of the interaction term is different from other effects.

Exercise 9

We want to continue exploring the association between age of first alcohol use and lifetime marijuana use, and how this association varies between males and females. Try to replicate the following graphs, and interpret what you see. Note also that something interesting happens at the age of 18. What is special about that age, according to the two graphs? Is the difference between the predicted outcomes of males and females significant at any age bracket? If so, what is the age bracket?

# probability plot
plot_predictions(
  mdl.glm.adjusted.interaction, 
  condition = c("alc_agefirst", "demog_sex")
  ) + 
  labs(x = "Alcohol (age first used)", 
       y = "Predicted probability of lifetime marijuana use",
       color = "Sex",
       fill = "Sex",
       title = "Predicted probability") + 
  theme_bw()

# RD
plot_comparisons(
  mdl.glm.adjusted.interaction, 
  variables = list(
    demog_sex = c("Male", "Female")
    ),
  comparison = "difference",
  condition = "alc_agefirst") + 
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(x = "Alcohol (age first used)", 
       y = "Pr(y|female) - Pr(y|male)",
       title = "Probability difference") + 
  theme_bw()

# OR
plot_comparisons(
  mdl.glm.adjusted.interaction, 
  variables = list(
    demog_sex = c("Male", "Female")
    ), 
  comparison = "lnor",
  transform = "exp",
  condition = "alc_agefirst") + 
  geom_hline(yintercept = 1, linetype = "dashed") + 
  labs(x = "Alcohol (age first used)", 
       y = "odds(y|female)/odds(y|male)",
       title = "Odds ratio") + 
  theme_bw()

# log OR
plot_comparisons(
  mdl.glm.adjusted.interaction, 
  variables = list(
    demog_sex = c("Male", "Female")
    ), 
  comparison = "lnor",
  condition = "alc_agefirst") + 
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(x = "Alcohol (age first used)", 
       y = "log(odds(y|female)/odds(y|male))",
       title = "Log Odds ratio") + 
  theme_bw()

# RR
plot_comparisons(
  mdl.glm.adjusted.interaction, 
  variables = list(
    demog_sex = c("Male", "Female")
    ), 
  comparison = "ratio",
  condition = "alc_agefirst") + 
  geom_hline(yintercept = 1, linetype = "dashed") + 
  labs(x = "Alcohol (age first used)", 
       y = "Pr(y|female)/Pr(y|male)",
       title = "Probability ratio") + 
  theme_bw()

Answer: The figure clearly demonstrates that the risk of using marijuana decreases as age decreases, and this difference is significant between males and females below and above the age of 18. Males show a smoother increase in risk as they get younger, while females experience a sharp increase in risk. The risk difference plot indicates a significant difference between males and females aged 10 to 13. Additionally, the probability ratio plot shows a significant difference across all ages.