library(tidyverse)
## Warning: package 'purrr' was built under R version 4.4.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── 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)
## Warning: package 'lmtest' was built under R version 4.4.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.4.3
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(sandwich)
## Warning: package 'sandwich' was built under R version 4.4.3
traffic_data <- read.csv("traffic_data.csv")
data_clean <- traffic_data %>%
drop_na()
unique(data_clean$Insurance)
## [1] "No" "Yes"
model <- lm(ViolFine ~ InsuranceOffense, data = data_clean)
bp_test <- bptest(model)
set.seed(123) # for reproducibility
resids_sample <- sample(residuals(model), size = 5000)
shapiro.test(resids_sample)
##
## Shapiro-Wilk normality test
##
## data: resids_sample
## W = 0.88915, p-value < 2.2e-16
shapiro_test <- shapiro.test(resids_sample)
data_clean$LogFine <- log(data_clean$ViolFine)
model1<- lm(ViolFine ~ InsuranceOffense, data = data_clean)
robust_se <- vcovHC(model1, type = "HC3")
summary(model1)
##
## Call:
## lm(formula = ViolFine ~ InsuranceOffense, data = data_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -163.42 -84.11 13.89 46.89 435.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 84.1053 0.8404 100.08 <2e-16 ***
## InsuranceOffenseYes 79.3098 2.2366 35.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 76.1 on 9546 degrees of freedom
## Multiple R-squared: 0.1164, Adjusted R-squared: 0.1163
## F-statistic: 1257 on 1 and 9546 DF, p-value: < 2.2e-16
coeftest(model1, vcov = vcovHC(model1, type = "HC3"))
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 84.10530 0.81671 102.981 < 2.2e-16 ***
## InsuranceOffenseYes 79.30975 2.53463 31.291 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
bp_test_log <- bptest(model1)
shapiro_test_log <- shapiro.test(resids_sample)
summary(model1)
##
## Call:
## lm(formula = ViolFine ~ InsuranceOffense, data = data_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -163.42 -84.11 13.89 46.89 435.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 84.1053 0.8404 100.08 <2e-16 ***
## InsuranceOffenseYes 79.3098 2.2366 35.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 76.1 on 9546 degrees of freedom
## Multiple R-squared: 0.1164, Adjusted R-squared: 0.1163
## F-statistic: 1257 on 1 and 9546 DF, p-value: < 2.2e-16
bp_test_log
##
## studentized Breusch-Pagan test
##
## data: model1
## BP = 52.477, df = 1, p-value = 4.354e-13
shapiro_test_log
##
## Shapiro-Wilk normality test
##
## data: resids_sample
## W = 0.88915, p-value < 2.2e-16
coeftest(model1, vcov = robust_se)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 84.10530 0.81671 102.981 < 2.2e-16 ***
## InsuranceOffenseYes 79.30975 2.53463 31.291 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1