library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.3     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(readxl)
library(pastecs)
## 
## Attaching package: 'pastecs'
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(MASS)
## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select
load("NSDUH_2023.Rdata")
nsduh <- data
nsduh$drinking_days <- ifelse(
  nsduh$IRALCFM %in% c(91, 93), 0,
  ifelse(
    nsduh$IRALCFM %in% 1:30,
    nsduh$IRALCFM,
    NA_real_
  )
)
nsduh$age_group <- factor(nsduh$AGE3)
nsduh$sex_group <- factor(nsduh$IRSEX)
model2 <- lm(drinking_days ~ age_group + sex_group, data = nsduh)
summary(model2)
## 
## Call:
## lm(formula = drinking_days ~ age_group + sex_group, data = nsduh)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.384 -4.185 -1.178  0.478 30.056 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.52913    0.10301   5.137  2.8e-07 ***
## age_group2   0.15880    0.13797   1.151  0.24976    
## age_group3   0.42195    0.13975   3.019  0.00253 ** 
## age_group4   1.33330    0.13335   9.998  < 2e-16 ***
## age_group5   3.65563    0.13065  27.979  < 2e-16 ***
## age_group6   3.68228    0.14178  25.972  < 2e-16 ***
## age_group7   4.16164    0.13959  29.814  < 2e-16 ***
## age_group8   4.40893    0.13033  33.828  < 2e-16 ***
## age_group9   4.70760    0.11412  41.252  < 2e-16 ***
## age_group10  4.80511    0.13204  36.392  < 2e-16 ***
## age_group11  4.85448    0.13148  36.922  < 2e-16 ***
## sex_group2  -1.00712    0.05109 -19.712  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.048 on 56618 degrees of freedom
##   (75 observations deleted due to missingness)
## Multiple R-squared:  0.08655,    Adjusted R-squared:  0.08637 
## F-statistic: 487.7 on 11 and 56618 DF,  p-value: < 2.2e-16
plot(model2, which = 1)

raintest(model2)
## 
##  Rainbow test
## 
## data:  model2
## Rain = 0.97315, df1 = 28315, df2 = 28303, p-value = 0.989
dwtest(model2)
## 
##  Durbin-Watson test
## 
## data:  model2
## DW = 1.9959, p-value = 0.3139
## alternative hypothesis: true autocorrelation is greater than 0
plot(model2, which = 3)

library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
bptest(model2)
## 
##  studentized Breusch-Pagan test
## 
## data:  model2
## BP = 3100.6, df = 11, p-value < 2.2e-16
plot(model2, which = 2)

set.seed(123)

r_sample <- sample(residuals(model2), 5000)

shapiro.test(r_sample)
## 
##  Shapiro-Wilk normality test
## 
## data:  r_sample
## W = 0.72398, p-value < 2.2e-16
plot(model2, which = 2)

vif(model2)
##               GVIF Df GVIF^(1/(2*Df))
## age_group 1.004114 10        1.000205
## sex_group 1.004114  1        1.002055

  1. Based on my results, linear regression model meets some assumption but violates others. For linearity, the Rainbow test had a p value of 0.989, which is greater than 0.05. This suggests that there is no significant evidence of a nonlinear relationship, so the assumption of the linearity appears to be met. For independence of errors, the Durbin Watson statistic was 1.996 with a p value of 0.624. Since the statistic is close to 2 and the p value us greater than 0.05, there is no evidence of significant correlation. Therefore, the assumption appears to be met. For homoscedasticity, the Breush Pagan test has a p value less than 0.001. This means the residuals do not have constant variance, so my model vioates the assumption. For normality of residuals, the Shapiro Wilk test had a W value of 0.724 and a p value less than 0.001. This suggests that the residuals are not normally distributed, meaning the model also violates the assumption of normality. Lastly, the multicollinearity, the adjusted GVIF values for age group and sex were approximately 1.00. Since these values are well below 5, there is no evidence of problematic multicollinearity, so this assumption is met.

  2. Based on my results, my model violates two assumptions: homoscedasticity and normality of residuals. The Breusch-Pagan test had a p-value less than 0.001, indicating that the residuals do not have constant variance. The Shapiro-Wilk test also had a p-value less than 0.001, suggesting that the residuals are not normally distributed. These violations may be related to the number of days participants reported drinking alcohol, especially since some participants reported zero drinking days.

model2_log <- lm(log1p(drinking_days) ~ age_group + sex_group,
                 data = nsduh)

summary(model2_log)
## 
## Call:
## lm(formula = log1p(drinking_days) ~ age_group + sex_group, data = nsduh)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1408 -0.9597 -0.1004  0.6034  3.3641 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.081954   0.016439   4.985 6.20e-07 ***
## age_group2   0.056020   0.022019   2.544    0.011 *  
## age_group3   0.153281   0.022303   6.873 6.37e-12 ***
## age_group4   0.409981   0.021282  19.264  < 2e-16 ***
## age_group5   0.976772   0.020851  46.844  < 2e-16 ***
## age_group6   0.977446   0.022627  43.198  < 2e-16 ***
## age_group7   1.022254   0.022277  45.889  < 2e-16 ***
## age_group8   1.045175   0.020800  50.249  < 2e-16 ***
## age_group9   1.058868   0.018213  58.139  < 2e-16 ***
## age_group10  1.012558   0.021072  48.052  < 2e-16 ***
## age_group11  0.921844   0.020983  43.933  < 2e-16 ***
## sex_group2  -0.134827   0.008154 -16.536  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9652 on 56618 degrees of freedom
##   (75 observations deleted due to missingness)
## Multiple R-squared:  0.1436, Adjusted R-squared:  0.1434 
## F-statistic: 862.7 on 11 and 56618 DF,  p-value: < 2.2e-16
plot(model2_log, which = 3)

bptest(model2_log)
## 
##  studentized Breusch-Pagan test
## 
## data:  model2_log
## BP = 7821.4, df = 11, p-value < 2.2e-16
plot(model2_log, which = 2)

r_sample_log <- sample(residuals(model2_log), 5000)

shapiro.test(r_sample)
## 
##  Shapiro-Wilk normality test
## 
## data:  r_sample
## W = 0.72398, p-value < 2.2e-16

To address the assumptions my original model violated, I tried a log transformation of my dependent variable, drinking days. Since some participants reported zero drinking days, I tested the assumptions again using the Breusch-Pagan and Shapiro-Wilk tests. The Shapiro-Wilk W value increased from 0.724 to 0.909, suggesting that the residuals became closer to a normal distribution. However, the p-value was still less than 0.001, meaning the model continued to violate the assumption of normality. The Breusch-Pagan test also remained significant, with a p-value less than 0.001, showing that the model still had unequal variance.Overall, the log transformation helped improve the normality of the residuals, but it did not completely fix the violations. This suggests that a different approach may be needed to better account for the distribution of drinking days.