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(janitor)
## 
## Attaching package: 'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test
k12 <- read_csv("Violence Project - K-12 SCHOOLS DATABASE 2026 - Event.csv")
## Rows: 398 Columns: 58
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (46): Full Date, Day of Week, Time of Day, During Class?, School Name, ...
## dbl  (10): ID #, Month, Day, Year, Latitude, Longitude, Number of
## Victims Ki...
## num   (1): Enrollment
## time  (1): Time Homicide
## Occurred
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
higher_ed <- read_csv("Violence Project - HIGHER ED DATABASE 2026 - Event.csv")
## Rows: 280 Columns: 60
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (48): Full Date, Day of Week, Time of Day, During Class?, School Name, ...
## dbl  (10): ID #, Month, Day, Year, Latitude, Longitude, Number of
## Victims Ki...
## num   (1): Enrollment
## time  (1): Time Homicide
## Occurred
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
k12 <- k12 %>%
  mutate(School_Type_Group = "K-12")

higher_ed <- higher_ed %>%
  mutate(School_Type_Group = "Higher Ed")

all_shootings <- bind_rows(k12, higher_ed)

all_shootings <- all_shootings %>%
  clean_names()

View(all_shootings)
# Create a linear model "lm()" from the variables, with a continuous dependent variable as the outcome
# Research Question: Have school shootings become more deadly over time? (Independent variable: Year, Dependent variable: num fatalities per shooting per year)

#average per year
yearly_data <- all_shootings %>%
  group_by(year) %>%
  summarise(avg_fatalities = mean(number_of_victims_killed, na.rm = TRUE))

View(yearly_data)
#make a linear model
model1 <- lm(avg_fatalities ~ year, data = yearly_data)

summary(model1)
## 
## Call:
## lm(formula = avg_fatalities ~ year, data = yearly_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.39267 -0.19896 -0.14046 -0.03545  2.04786 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  9.381985  27.268602   0.344    0.734
## year        -0.003987   0.013550  -0.294    0.771
## 
## Residual standard error: 0.5182 on 24 degrees of freedom
## Multiple R-squared:  0.003594,   Adjusted R-squared:  -0.03792 
## F-statistic: 0.08657 on 1 and 24 DF,  p-value: 0.7711

My p-vale is 0.771 which is above .05 which means the difference by year is not statistically significant. I do not have enough evidence to conclude that average fatalities per shooting have changed over time.

The coefficient, which is -0.003987, is the number that tells me if the relationship is negative or positive. Since it’s negative, it actually tells me that there has been a very, very small decrease over time. This is interesting because this is different than what one of my scholarly sources was telling me.

My r-squared is 0.003594 whih tells me how much this model can explain the variation in my dependent varaible- the change in the number of folks killed on average by year. My model only explains .36% which is astonishingly little.

All in all, none of these numbers show that school shooting fatalities have increased over time.

#For question 3, I need to check the following assumptions: Linearity (plot and raintest)

#plot
plot(model1, which = 1)

# now the raintest. First, I load the package. Then I perform test
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
raintest(model1)
## 
##  Rainbow test
## 
## data:  model1
## Rain = 0.24789, df1 = 13, df2 = 11, p-value = 0.9901

My plot test showed that most of my dots are on the red line, but a few plot points are not on the line. This mostly supports the assumption of linearity.

Then, I ran a raintest. The test shows a p-value of 0.9901 which is greater than .05. This means that the test did not show significant evidence of linearity violation. In the end, this shows that I essentially could use a straight line to represent my data. The answer to beginning question of whether shootings have gotten more deadly over time is shown by this data to be “no”. There are some outliers but not enough to prove that my model violates the assumption of linearity.

# Now I have to test: Independence of errors (durbin-watson)

dwtest(model1)
## 
##  Durbin-Watson test
## 
## data:  model1
## DW = 2.1268, p-value = 0.5438
## alternative hypothesis: true autocorrelation is greater than 0

Is my model able to relatively accurately predict in the future the number of fatalities in school shootings each year. Knowing how off my predictions are in one year can help me understand if they will be off the next year too. I ran the Durbin-Watson test to test the assumption of independence. I got a statistic of 2.1268 (2= little evidence of autocorrelation, 0= errors move together, 4= errors move in opposite directions). I got a p-value of 0.5438 which is greater than 0.05. This shows me that there isn’t significant evidence of positive autocorrelation among residuals Therefore, the assumption of independence of errors= reasonable.

# Now I need to test how spread out my errors are. The spread should be consistent. Test Homoscedasticity (plot, bptest)

plot(model1, which = 3)

bptest(model1)
## 
##  studentized Breusch-Pagan test
## 
## data:  model1
## BP = 0.62201, df = 1, p-value = 0.4303

This test returned a p-value of .4303 which is greater than .05. This means that the test did not find significant evidence that the spear of my prediction errors is inconsistent. Still, there are couple of outliers. This suggest that the assumption of homoscedasticity is reasonable.

# Now I will test the normality of residuals (QQ plot, shapiro test)

# First, I ceate the QQ plot
plot(model1, which = 2)

# Shapiro test
shapiro.test(model1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  model1$residuals
## W = 0.60049, p-value = 2.995e-07

It looks like the impact of those outliers if finally showing up in the tests. In the QQ plot, the dots should stick close to the line, but some observations are far from it. Here, my p-value is 2.995e-7. This is less than .05 which means there is evidence of statistical significance. This test shows that my model violates the assumption of normality of residuals.

# Now I need to run the kitchen sink model. So far, I have only been examining one independent variable, so now I will bring in other factors I am interested in looking at to create a model for this question. 

# Kitchen sink model
kitchen_sink <- lm(number_of_victims_killed ~ school_type_group + armed_person_on_scene + accessibility_to_public + during_class, data = all_shootings)

summary(kitchen_sink)
## 
## Call:
## lm(formula = number_of_victims_killed ~ school_type_group + armed_person_on_scene + 
##     accessibility_to_public + during_class, data = all_shootings)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6960 -0.6497 -0.0232 -0.0048 30.1862 
## 
## Coefficients:
##                                                                Estimate
## (Intercept)                                                     1.03698
## school_type_groupK-12                                          -0.03220
## armed_person_on_sceneYes                                        0.64494
## accessibility_to_publicRestricted Space and Didn't Have Access  1.86841
## accessibility_to_publicRestricted Space but Had Access         -0.01381
## accessibility_to_publicUnknown                                 -0.03698
## during_classYes                                                 0.79062
##                                                                Std. Error
## (Intercept)                                                       0.13095
## school_type_groupK-12                                             0.15473
## armed_person_on_sceneYes                                          0.28061
## accessibility_to_publicRestricted Space and Didn't Have Access    0.37557
## accessibility_to_publicRestricted Space but Had Access            0.19161
## accessibility_to_publicUnknown                                    0.94992
## during_classYes                                                   0.18792
##                                                                t value Pr(>|t|)
## (Intercept)                                                      7.919 9.98e-15
## school_type_groupK-12                                           -0.208   0.8352
## armed_person_on_sceneYes                                         2.298   0.0218
## accessibility_to_publicRestricted Space and Didn't Have Access   4.975 8.30e-07
## accessibility_to_publicRestricted Space but Had Access          -0.072   0.9426
## accessibility_to_publicUnknown                                  -0.039   0.9690
## during_classYes                                                  4.207 2.94e-05
##                                                                   
## (Intercept)                                                    ***
## school_type_groupK-12                                             
## armed_person_on_sceneYes                                       *  
## accessibility_to_publicRestricted Space and Didn't Have Access ***
## accessibility_to_publicRestricted Space but Had Access            
## accessibility_to_publicUnknown                                    
## during_classYes                                                ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.882 on 671 degrees of freedom
## Multiple R-squared:  0.08513,    Adjusted R-squared:  0.07695 
## F-statistic: 10.41 on 6 and 671 DF,  p-value: 4.764e-11
# Now I have to check multicollinearity
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
vif(kitchen_sink)
##                             GVIF Df GVIF^(1/(2*Df))
## school_type_group       1.111409  1        1.054234
## armed_person_on_scene   1.029911  1        1.014845
## accessibility_to_public 1.277475  3        1.041659
## during_class            1.219111  1        1.104134

Using the kitchen sink method, I examined school type, armed personnel on scene, public accessibility, and whether the shooting occured during class. The GVIF values are all close to one (not 5 or 10 which would be a serious concern). This means that my independent variables are not stongly correlated to one another (one variable doesn’t perfectly predict another variable which could indicate an issue with the data because nothing is ever perfect). These tests do not show an issue with multicollinearity.

I could not use cor() on this assignment because none of my numbers are continuous numerical, they are all whole numbers or yes/no type of answers. The best I could do was calculate yearly averages so at least I was working with decimals.

Questions 4-6: How does your model meet those assumptions? If your model violates an assumption, which one? What would you do to mitigate this assumption? Show your work.

The only test from above in which my model violated the assumption of residuals which I had to use the Shapiro-Wilk test to find.

First and foremost, I have to understand the concept behind what these outliers are.During this HW assignment, I googled “Are school shootings usually small” and found that- yes! they are!. The dataset I have been working with is cleary showing me that. Most school shootings kill one person, and mass school shootings are an outlier. I want to explore this more as I add onto my Lit Review for my paper.

On a mathematical level, one thing I can do is perform a logarithmic code to reel the impact of the outliers in.