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(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
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

1. Loading the Dataset

I am continuing to use the Texas Department of Family and Protective Services (DFPS) dataset from Homework 6. I want to examine whether fiscal year and victim confirmation status are related to the number of victims.

dfps_data <- read.csv("dfps_data.csv", check.names = FALSE)

head(dfps_data)
##   Fiscal Year Program     County    Region Confirmed Victims # of Victims
## 1        2025     DCI Hutchinson 1-Lubbock     Not Confirmed            1
## 2        2025     DCI    Lubbock 1-Lubbock     Not Confirmed           36
## 3        2025     DCI    Lubbock 1-Lubbock  Confirmed Victim            4
## 4        2025     DCI  Ochiltree 1-Lubbock     Not Confirmed            1
## 5        2025     DCI     Potter 1-Lubbock     Not Confirmed           13
## 6        2025     DCI     Potter 1-Lubbock  Confirmed Victim            1
str(dfps_data)
## 'data.frame':    3121 obs. of  6 variables:
##  $ Fiscal Year      : int  2025 2025 2025 2025 2025 2025 2025 2025 2025 2025 ...
##  $ Program          : chr  "DCI" "DCI" "DCI" "DCI" ...
##  $ County           : chr  "Hutchinson" "Lubbock" "Lubbock" "Ochiltree" ...
##  $ Region           : chr  "1-Lubbock" "1-Lubbock" "1-Lubbock" "1-Lubbock" ...
##  $ Confirmed Victims: chr  "Not Confirmed" "Not Confirmed" "Confirmed Victim" "Not Confirmed" ...
##  $ # of Victims     : chr  "1" "36" "4" "1" ...

2. Creating the Linear Model

My dependent variable is the number of victims. My independent variables are fiscal year and confirmation status.

regression_data <- data.frame( victims = as.numeric(dfps_data[["# of Victims"]]), fiscal_year = as.numeric(dfps_data[["Fiscal Year"]]), confirmation_status = as.factor(dfps_data[["Confirmed Victims"]]) )
## Warning in data.frame(victims = as.numeric(dfps_data[["# of Victims"]]), : NAs
## introduced by coercion
regression_data <- na.omit(regression_data)
victim_model <- lm(victims ~ fiscal_year + confirmation_status, data = regression_data)

summary(victim_model)
## 
## Call:
## lm(formula = victims ~ fiscal_year + confirmation_status, data = regression_data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -28.93 -22.13  -7.62  -0.34 895.21 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      -1670.3258   648.8012  -2.574  0.01008 *  
## fiscal_year                          0.8295     0.3211   2.583  0.00984 ** 
## confirmation_statusNot Confirmed    20.6207     1.9741  10.446  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 51.4 on 3113 degrees of freedom
## Multiple R-squared:  0.03571,    Adjusted R-squared:  0.03509 
## F-statistic: 57.65 on 2 and 3113 DF,  p-value: < 2.2e-16

In Homework 6, my model had an R-squared of approximately 0.0357. This means fiscal year and confirmation status explained about 3.57% of the variation in victim counts. Although the predictors were statistically significant, the model explained relatively little of the overall variation.

One limitation is that victim counts are count data rather than strictly continuous measurements. I am retaining this model to continue my previous analysis, but this should be considered when interpreting the results.

3. Testing the Assumptions of Linear Regression

A. Linearity

I use the residuals-versus-fitted plot and Rainbow test to examine linearity.

plot(victim_model, which = 1)

The residuals are generally scattered around zero, and the red line provides some visual support for linearity. However, there are some unusually large positive residuals.

victim_linearity_model <- lm(victims ~ fiscal_year, data = regression_data)

raintest(victim_linearity_model)
## 
##  Rainbow test
## 
## data:  victim_linearity_model
## Rain = 0.79169, df1 = 1558, df2 = 1556, p-value = 1

The Rainbow test p-value is 1, which is greater than 0.05. This means the test does not identify a significant linearity problem. The relationship between fiscal year and victim counts appears linear.

B. Independence of Errors

I use the Durbin-Watson test to examine independence of errors.

durbinWatsonTest(victim_model)
##  lag Autocorrelation D-W Statistic p-value
##    1      0.05954429       1.88081    0.01
##  Alternative hypothesis: rho != 0

The p-value is 0.012, which is below 0.05. This means the errors are not independent. The D-W statistic is 1.88081, which is relatively close to 2, but the significant p-value indicates autocorrelation.

This assumption is violated.

C. Homoscedasticity

I use the Scale-Location plot and Breusch-Pagan test to examine whether the variance of the residuals is consistent.

plot(victim_model, which = 3)

The red line and the spread of the points help determine whether the model is homoscedastic.

bptest(victim_model)
## 
##  studentized Breusch-Pagan test
## 
## data:  victim_model
## BP = 19.912, df = 2, p-value = 4.745e-05

The p-value is 4.745e-05, which is below 0.05. This means the model is heteroscedastic, and the assumption of homoscedasticity is violated.

D. Normality of Residuals

I use a Q-Q plot and Shapiro-Wilk test to examine normality.

plot(victim_model, which = 2)

shapiro.test(victim_model$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  victim_model$residuals
## W = 0.4269, p-value < 2.2e-16

The Shapiro-Wilk p-value is below 0.05, indicating that the residuals are not normally distributed.

The unusually large positive residuals also suggest departures from normality, particularly in the upper tail.

This assumption is violated.

E. No Multicollinearity

I use VIF to examine whether the independent variables are strongly correlated.

vif(victim_model)
##         fiscal_year confirmation_status 
##            1.000087            1.000087

Both VIF values are approximately 1.000087, which is well below 5. This means there is no significant multicollinearity problem between my predictors.

I also examine the correlation between my numeric variables.

numeric_variables <- regression_data[c("victims", "fiscal_year")]

cor(numeric_variables)
##               victims fiscal_year
## victims     1.0000000   0.0437523
## fiscal_year 0.0437523   1.0000000

The correlation shows the direction and strength of the relationship between fiscal year and victim counts. Confirmation status is categorical, so I do not include it in this Pearson correlation calculation.

4. Does My Model Meet the Assumptions?

My model meets the linearity assumption based on the Rainbow test, and the VIF results show no multicollinearity problem.

However, the Durbin-Watson test indicates that the errors are not independent. The Breusch-Pagan test indicates heteroscedasticity, and the Shapiro-Wilk test indicates that the residuals are not normally distributed.

Overall, my model does not meet all five assumptions.

5. Which Assumptions Are Violated?

The assumptions that are violated are:

  1. Independence of errors: The Durbin-Watson p-value is 0.012.
  2. Homoscedasticity: The Breusch-Pagan p-value is 4.745e-05.
  3. Normality: The Shapiro-Wilk p-value is below 0.05.

Linearity is supported by the Rainbow test, and multicollinearity does not appear to be a problem because both VIF values are approximately 1.

6. What Would I Do to Mitigate These Assumptions?

I will try a log transformation and a robust linear model to see whether they improve the results.

A. Log Transformation

I use a log transformation to try to improve linearity, homoscedasticity, and normality. I add 1 to the victim count to account for possible zero values.

victim_model_log <- lm(log(victims + 1) ~ fiscal_year + confirmation_status, data = regression_data)

summary(victim_model_log)
## 
## Call:
## lm(formula = log(victims + 1) ~ fiscal_year + confirmation_status, 
##     data = regression_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6051 -0.8268 -0.2835  0.7201  4.6429 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      -44.915472  15.190445  -2.957  0.00313 ** 
## fiscal_year                        0.022931   0.007519   3.050  0.00231 ** 
## confirmation_statusNot Confirmed   0.778282   0.046219  16.839  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.203 on 3113 degrees of freedom
## Multiple R-squared:  0.08573,    Adjusted R-squared:  0.08515 
## F-statistic:   146 on 2 and 3113 DF,  p-value: < 2.2e-16

The transformed model has an R-squared of 0.08573, compared with approximately 0.0357 in my original model. This means the transformed model explains approximately 8.57% of the variation in the transformed outcome.

Checking Linearity

plot(victim_model_log, which = 1)

victim_linearity_model_log <- lm(log(victims + 1) ~ fiscal_year, data = regression_data)

raintest(victim_linearity_model_log)
## 
##  Rainbow test
## 
## data:  victim_linearity_model_log
## Rain = 1.0116, df1 = 1558, df2 = 1556, p-value = 0.4099

The Rainbow test p-value is 0.4099, which is greater than 0.05. This means the transformed model does not show a significant linearity problem.

The log transformation appears to maintain the linearity assumption.

Checking Independence of Errors

durbinWatsonTest(victim_model_log)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.2011539      1.597035       0
##  Alternative hypothesis: rho != 0

The p-value is reported as 0, which is below 0.05. This means the errors are still not independent.

The log transformation did not resolve this problem.

Checking Homoscedasticity

plot(victim_model_log, which = 3)

bptest(victim_model_log)
## 
##  studentized Breusch-Pagan test
## 
## data:  victim_model_log
## BP = 198.91, df = 2, p-value < 2.2e-16

The Breusch-Pagan p-value is below 2.2e-16, which is below 0.05. This means the transformed model is still heteroscedastic.

The log transformation did not resolve this assumption.

Checking Normality

plot(victim_model_log, which = 2)

shapiro.test(victim_model_log$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  victim_model_log$residuals
## W = 0.9331, p-value < 2.2e-16

The Shapiro-Wilk p-value is below 0.05, meaning the residuals are still not normally distributed.

However, the W statistic increased from 0.4269 in the original model to 0.9331 in the transformed model. This suggests that the transformation improved normality, although the assumption is still violated.

Checking Multicollinearity

vif(victim_model_log)
##         fiscal_year confirmation_status 
##            1.000087            1.000087

The VIF values remain approximately 1.000087. Multicollinearity is still not a problem.

Overall, the log transformation improved the R-squared value and the normality statistic, but it did not resolve independence or homoscedasticity.

B. Robust Linear Model (RLM)

I also try a robust linear model because my original model contained unusually large residuals.

victim_model_rlm <- rlm(victims ~ fiscal_year + confirmation_status,data = regression_data)

summary(victim_model)
## 
## Call:
## lm(formula = victims ~ fiscal_year + confirmation_status, data = regression_data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -28.93 -22.13  -7.62  -0.34 895.21 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      -1670.3258   648.8012  -2.574  0.01008 *  
## fiscal_year                          0.8295     0.3211   2.583  0.00984 ** 
## confirmation_statusNot Confirmed    20.6207     1.9741  10.446  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 51.4 on 3113 degrees of freedom
## Multiple R-squared:  0.03571,    Adjusted R-squared:  0.03509 
## F-statistic: 57.65 on 2 and 3113 DF,  p-value: < 2.2e-16
summary(victim_model_rlm)
## 
## Call: rlm(formula = victims ~ fiscal_year + confirmation_status, data = regression_data)
## Residuals:
##     Min      1Q  Median      3Q     Max 
##  -9.153  -6.334  -2.606   5.903 911.666 
## 
## Coefficients:
##                                  Value     Std. Error t value  
## (Intercept)                      -326.4907  116.3886    -2.8052
## fiscal_year                         0.1637    0.0576     2.8424
## confirmation_statusNot Confirmed    5.0736    0.3541    14.3269
## 
## Residual standard error: 9.155 on 3113 degrees of freedom

I compare the regular and robust regression coefficients to see whether unusually large observations are influencing the results.

The robust model provides another way to examine the data, but it does not automatically resolve the independence or homoscedasticity problems.

Unlike regular linear regression, the robust model does not provide the same p-values or R-squared values.

Conclusion

My original model met the linearity and multicollinearity assumptions, but violated independence, homoscedasticity, and normality.

I tried a log transformation, which improved the R-squared from approximately 0.0357 to 0.08573 and improved the normality statistic. However, the transformed model still violated independence, homoscedasticity, and normality.

I also fitted a robust linear model to examine the influence of unusually large residuals.

Overall, the log transformation improved some aspects of the model, but neither approach establishes that all regression assumptions have been satisfied.