Loading the Data

setwd("~/Documents/R Files/Homework 6")
districts <- read_csv("district_info.csv", skip=1)

2) Create a linear model “lm()” from the variables, with a continuous dependent variable as the outcome

We will be using the same model we used in Homework 6

model <- lm(DPETRSKP ~ DPETECOP + DPETLEPP + DPETSPEP + DPETALLC,
            data = districts)

summary(model)

Call:
lm(formula = DPETRSKP ~ DPETECOP + DPETLEPP + DPETSPEP + DPETALLC, 
    data = districts)

Residuals:
    Min      1Q  Median      3Q     Max 
-41.294  -6.432  -0.065   6.209  77.657 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 6.189e+00  1.562e+00   3.963 7.85e-05 ***
DPETECOP    3.503e-01  2.035e-02  17.214  < 2e-16 ***
DPETLEPP    5.570e-01  3.060e-02  18.202  < 2e-16 ***
DPETSPEP    6.803e-01  7.247e-02   9.387  < 2e-16 ***
DPETALLC    5.799e-05  3.200e-05   1.812   0.0703 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.22 on 1203 degrees of freedom
Multiple R-squared:  0.5218,    Adjusted R-squared:  0.5203 
F-statistic: 328.2 on 4 and 1203 DF,  p-value: < 2.2e-16

3) Check the following assumptions:

3a) Linearity (plot and raintest)

plot(model, which = 1)

raintest(model)

    Rainbow test

data:  model
Rain = 1.0406, df1 = 604, df2 = 599, p-value = 0.3128

3b) Independence of errors (durbin-watson)

durbinWatsonTest(model)
 lag Autocorrelation D-W Statistic p-value
   1       0.0800619        1.8396   0.008
 Alternative hypothesis: rho != 0

3c) Homoscedasticity (plot, bptest)

plot(model, which = 1)

bptest(model)

    studentized Breusch-Pagan test

data:  model
BP = 72.162, df = 4, p-value = 7.931e-15

3d) Normality of residuals (QQ plot, shapiro test)

plot(model, which = 2)

shapiro.test(residuals(model))

    Shapiro-Wilk normality test

data:  residuals(model)
W = 0.92362, p-value < 2.2e-16

3e) No multicolinarity (VIF, cor)

vif(model)
DPETECOP DPETLEPP DPETSPEP DPETALLC 
1.322277 1.503743 1.119100 1.078609 
cor(districts[, c("DPETECOP",
                  "DPETLEPP",
                  "DPETSPEP",
                  "DPETALLC")],
    use = "complete.obs")
           DPETECOP   DPETLEPP    DPETSPEP    DPETALLC
DPETECOP 1.00000000  0.4506561  0.05043960  0.01633756
DPETLEPP 0.45065607  1.0000000 -0.26397594  0.24665550
DPETSPEP 0.05043960 -0.2639759  1.00000000 -0.05921799
DPETALLC 0.01633756  0.2466555 -0.05921799  1.00000000

4) does your model meet those assumptions? You don’t have to be perfectly right, just make a good case. The model meets the multicollinearity assumption, with all VIF values below 2 and relatively low correlations among the independent variables. However, the model shows some assumption violations. The Durbin-Watson test was significant (p = 0.004), suggesting some dependence in the errors. The Breusch-Pagan test was also significant, indicating heteroscedasticity. Finally, the Q-Q plot and Shapiro-Wilk test show that the residuals are not normally distributed.

5) If your model violates an assumption, which one? The clearest violations are homoscedasticity and normality of residuals. The residual plot shows uneven spread, and the Breusch-Pagan test is significant. The Q-Q plot also shows substantial departures from the reference line, especially in the tails, and the Shapiro-Wilk test is significant. The Durbin-Watson test also suggests some autocorrelation.

6) What would you do to mitigate this assumption? Show your work.

I would focus on heteroscedasticity and non-normal residuals because both tests were statistically significant. To try to improve these assumptions, I square-root transformed the dependent variable and reran the model. I then repeated the residual, Breusch-Pagan, Q-Q, and Shapiro-Wilk checks to see whether the model improved.

districts$sqrt_atrisk <- sqrt(districts$DPETRSKP)

model2 <- lm(sqrt_atrisk ~ DPETECOP + DPETLEPP + DPETSPEP + DPETALLC,
             data = districts)

plot(model2, which = 1)

bptest(model2)

    studentized Breusch-Pagan test

data:  model2
BP = 72.94, df = 4, p-value = 5.432e-15
plot(model2, which = 2)

shapiro.test(residuals(model2))

    Shapiro-Wilk normality test

data:  residuals(model2)
W = 0.92831, p-value < 2.2e-16

After log-transforming district enrollment and rerunning the model, the Breusch-Pagan test remained statistically significant (p < 2.2e-16). This means the heteroscedasticity was not fully resolved, even after the transformation.