library(readr)
library(readxl)
library(skimr)
library(flextable)
library(dplyr)
library(stringr)
library(tidyr)
library(lubridate)
library(stringr)
library(ggplot2)
library(forcats)
library(datadictionary)
library(tidyverse)

MassShootings <- read_csv("Mother Jones - Mass Shootings Database, 1982 - 2022 - Sheet1.csv")

State_level_data <- read_excel("Sate level data & factor score.xlsx")

skim(MassShootings)

skim(State_level_data)

# Removing the two values as they were not provided and given as character. So replacing them with 0

MassShootings[7,6] = "0"
MassShootings[7,7] = "0"

# Changing the variable type character to numeric to perform group level summary statistics

MassShootings$injured <- as.numeric(as.character(MassShootings$injured))
MassShootings$total_victims <- as.numeric(as.character(MassShootings$total_victims))


# separating the location variable into two variables city and state and removing the space before state values using str_trim

# Using lubridate package and the function mdy to modify the date variable
MassShootings<- separate(MassShootings,location,c("City","STATE"),sep=",") %>% mutate(STATE=str_trim(STATE),date=mdy(date), Year=year(date))



# Making adjustments to the race variable that are capitalized differently and has no values using mutate and stringr function str_replace


MassShootings <- MassShootings %>% mutate(across('race',str_replace,"black","Black"))
MassShootings <- MassShootings %>% mutate(across('race',str_replace,"unclear","Unknown"))
MassShootings <- MassShootings %>% mutate(across('race',str_replace,"-","Unknown"))
MassShootings <- MassShootings %>% mutate(across('race',str_replace,"white","White"))

# Modifying the place variable with different capitalizations and wordings

MassShootings <- MassShootings %>% mutate(across('Place',str_replace,"workplace","Workplace"))
MassShootings <- MassShootings %>% mutate(across('Place',str_replace,"\nWorkplace","Workplace"))
MassShootings <- MassShootings %>% mutate(across('Place',str_replace,"religious","Religious"))
MassShootings <- MassShootings %>% mutate(across('Place',str_replace,"Other\n","Other"))

# creating a data dictionary using the function datadictionary and keeping only required data using slice function

Dictionary1<-datadictionary::create_dictionary(MassShootings)

# Slicing only required rows of the data dictionary

Dictionary1<- Dictionary1 %>% slice(3,5,7,9,11,13,18,23,28,30,32,34,36,38,40,42,44,46,48,50,52,54,56,58,65)

# updating the description of variables to the dictionary

Dictionary1$label<-c("case name assigned to the mass shooting ",
                     "City", 
                     "State",
                     "Date of incident",
                    "Summary of the Mass shooting",
                    "The number of people dead",
                    "The number of people injured",
                    "Total number of victims",
                     "Type of the place",
                    "Age of the shooter",
                    "Any prior signs for the offender's mental issues",
                    "Mental health details of the shooter",
                     "Legality of the weapons obtained",
                    "where was the weapons obtained from",
                    "Type of the weapon used",
                    "Details of the weapon",
                    "Race of the shooter",
                     "Gender of the shooter",
                    "Source of the data",
                    "Sources of the mental health",
                    "Additional sources",
                    "Latitude of the place of incident",
                     "Longitude of the place of incident",
                    "Type of Mass shooting",
                    "Year of mass shooting")

Dictionary1<- select(Dictionary1,1:3)

# renaming the variables

Dictionary1<- Dictionary1 %>% rename(Variable = item,Description=label,Type=class)
Dictionary1[Dictionary1 == "character"] <- "Categorical"
Dictionary1[Dictionary1 == "numeric"] <- "Quantitative"

flextable(Dictionary1,cwidth = 4,cheight = 1)

# Creating data dictionary for 2nd dataset of Household firearm estimate

Dictionary2<-datadictionary::create_dictionary(State_level_data)

#Slicing the data dictionary to keep the required rows 
Dictionary2<- Dictionary2 %>% slice(3,8,13,15,20,25,30,35,40,45,50,55,60,65,70,75,80,85,90,95)

# Selecting only required columns of the data ditionary
Dictionary2<- select(Dictionary2,1:3)

# renaming the variables

Dictionary2<- Dictionary2 %>% rename(Variable = item,Description=label,Type=class)
Dictionary2[Dictionary2 == "character"] <- "Categorical"
Dictionary2[Dictionary2 == "numeric"] <- "Quantitative"

# updating the description of variables to the dictionary

Dictionary2$Description<-c("State FIP value ",
                           "Year",
                          "Name of the State",
                          "Factor scores for household firearm ownership latent factor",
                          "Standard errors of factor scores for household firearm ownership latent factor",
                          "State has universal background checks law (1=yes; 0=no)",
                          "State has permit to purchase law (1=yes; 0=no)",
                          "Female Firearm Suicide/Total Female Suicide *100",
                          "Male Firearm Suicide/Total Male Suicide *100",
                          "BRFSS state-level survey estimate",
                          "Gallup state-level survey estimate",
                          "General Social Survey state-level survey estimate",
                          "Pew state-level survey estimate",
                          "Square root of Hunting Licenses/Population",
                          "Within-year standardization of the square root of Guns and Ammo Subscriptions/Population",
                          "Within-year standardization of Background Checks (without checks for permits)/Population",
                          "Indicator for Pew firearm ownership question wording change",
                          "Blended linear spline 1 represents roughly 1980–1992",
                          "Blended linear spline 1 represents roughly 1993–2004",
                          "Blended linear spline 1 represents roughly 2005–2016")

flextable(Dictionary2,cwidth = 4,cheight = 1)

# plot 1: Line chart showing in which year the most shootings took place

# creating a new df Yearcount that has the Count of mass shootings in the specific year

YearCount <- MassShootings %>% count(year)

massshootings_year<-YearCount %>% ggplot(aes(x = year, y = n)) +
  geom_line(size=1,color="Red") +
  scale_y_continuous(breaks = c(2,4,6,8,10,12))+
  scale_x_continuous(breaks = c(1980,1985,1990,1995,2000,2005,2010,2015,2020,2025))+
  labs(title = "Mass shootings in the United States from 1980-2022",
       x = "Year",
       y = "Number of incidents",
       caption = "Data source: www.motherjones.com") +
  theme(plot.title = element_text(hjust = 0.5))+
  theme_get()

massshootings_year

# creating a new df locationcount that has the values of number of mass shootings occured in each state
locationCount <- MassShootings %>% count(STATE)

massshootings_state<-locationCount %>% ggplot(aes(x = fct_reorder(STATE, n), y = n)) +
  geom_col(color="black",fill="#87CEEB") +
  coord_flip()+
  scale_y_continuous(breaks = c(0,4,8,12,16,20,24,28))+
  labs(title = "Number of shootings per state in the United States",
       x = "States",
       y = "Number of incidents",
       caption = "Data source: www.motherjones.com") +
  theme_bw()+
  theme(axis.text.y  = element_text(size = 7))

massshootings_state

# Removing the Las vegas mass shooting data as it was only big outlier that alters the results of analysis and for better visualization 
massshootings2 <- MassShootings[-45,]

massshootings_race<-massshootings2 %>% mutate(race=fct_relevel(race,"White","Black","Unknown","Latino","Asian","Other","Native American"))%>%
  ggplot(aes(x=race,y=total_victims,fill=race)) +
  geom_boxplot(color="black") +
  scale_y_continuous(breaks = c(0,10,20,30,40,50,60,70,80,90,100,110,120))+
  scale_fill_viridis_d()+
  labs(title = "Victims in mass shootings by suspects of different ethnicities ",
       x = "Race of the criminal",
       y = "Number of victims",
       caption = "Data source: www.motherjones.com")+
  theme_light()+
  theme( legend.position = "bottom", panel.border = element_rect(color = "black", fill = NA),
         axis.line = element_line(color = "black"),
      axis.ticks = element_line(color = "black"),
      axis.text = element_text(color = "black",size = 7))

massshootings_race

# Creating a new table from the data set with data required for the plot from the state of mississippi

Mississippi_firearm_suicide<-filter(State_level_data,STATE=="Mississippi")%>%select("STATE","Year","Male_FS_S","Fem_FS_S")

# Multiplying the variables with 100 as they are very small in number and hard to visualize if plotted with same values
Mississippi_firearm_suicide <- Mississippi_firearm_suicide %>% mutate(Male_FS_S=Male_FS_S*100,Fem_FS_S=Fem_FS_S*100)

# renaming column names to male and female instead of Male_FS_S, Fem_FS_S

Mississippi_firearm_suicide <-  rename(Mississippi_firearm_suicide, Male=Male_FS_S,Female=Fem_FS_S)

# Pivoting the variables to make it easier to visualize both genders in the same plot
Mississippi_firearm_suicide <- Mississippi_firearm_suicide %>% pivot_longer(cols = c(Male,Female),names_to = "Gender",values_to ="Firearm_sucide_rate")

# Plot 4
mississippi_suicide_gender<-Mississippi_firearm_suicide %>% ggplot(aes(x=Year,y=Firearm_sucide_rate,color=Gender))+
  geom_line(size=1)+
  geom_point(shape=16, color="black", fill="black", size=1)+
  scale_x_continuous(breaks = c(1980,1985,1990,1995,2000,2005,2010,2015,2020))+
  scale_color_discrete(labels = c("Female", "Male"))+
  labs(title = "Firearm suicide rate in mississippi based on gender",
       x="Year",
       Y="Rate of firearm suicide",
       caption = "RAND : State-Level Estimates of Household Firearm Ownership")+
  theme_light()+
  theme(legend.position = "bottom",panel.border = element_rect(color = "black", fill = NA),
         axis.line = element_line(color = "black"),
      axis.ticks = element_line(color = "black"),
      axis.text = element_text(color = "black",size = 7))

mississippi_suicide_gender

# Joining the two data sets by Year and STATE variables to keep the Variables Year,STATE,Fem_FS_S,Male_FS_S to find out male suicides from 1980-2020

Newtable <- left_join(MassShootings,select(State_level_data,c(Year,STATE,Fem_FS_S,Male_FS_S)),by=c("Year","STATE"))

# Removing the NA values in the newtable as there are  53 missing values after left joining the data sets

Newtable <- Newtable[complete.cases(Newtable),]

# Multiplying the variable by 100 to make it visually better in the plot

Newtable$Fem_FS_S <- Newtable$Fem_FS_S*100
Newtable$Male_FS_S <- Newtable$Male_FS_S*100


# The  Male_FS_S    variable is - Male Firearm Suicide/Total Male Suicide *100

# creating a plot to show the male fire arm suicide percentage

male_suicide_percentage<-Newtable %>% ggplot(aes(x=Year,y=Male_FS_S))+
  geom_col(color="red",fill="red")+
  scale_x_continuous(breaks = c(1980,1985,1990,1995,2000,2005,2010,2015,2020))+
  labs(title = " Percentage of Male firearm suicide per year ",
       subtitle = "The male firearm suicide variable is multiplied by 100 to make it visually better",
       x = "Year",
       y = "percentage of Male firearm suicides",
       caption = "Data source: www.motherjones.com & RAND : State-Level Estimates of Household Firearm Ownership") +
  theme_bw()
  

male_suicide_percentage

# 1- This data frame is for the summary statistics of firearm suicide rate in mississipi based on gender

firearm_suicide_gender<- Mississippi_firearm_suicide %>% group_by(Gender) %>% 
  summarise(mean_f=round(mean(Firearm_sucide_rate),2),
            median_f=round(median(Firearm_sucide_rate),2),
            sd_f=round(sd(Firearm_sucide_rate),2),
            max_f=round(max(Firearm_sucide_rate),2),
            min_f=round(min(Firearm_sucide_rate),2),.groups = 'drop') %>%
  as.data.frame()

 firearm_suicide_gender

# 2- Summary statistics data frame of number of victims in mass shootings by ethnicity of suspect
total_victims_race<- massshootings2 %>% group_by(race) %>% 
  summarise(mean_v=round(mean(total_victims),2),
            median_v=round(median(total_victims),2),
            sd_v=round(sd(total_victims),2),
            max_v=round(max(total_victims),2),
            min_v=round(min(total_victims),2),.groups = 'drop') %>%
  as.data.frame()

total_victims_race

freq_table<- table(MassShootings$Place, MassShootings$race)
freq_table 

as.data.frame(freq_table)

# 2b) Implement at least one permutation test based on a traditional hypothesis test, such as a two-sample t-test or a chi-squared test of independence, to test a hypothesis of interest for your data.
# creating a new variable shooter which determines if the shooter is minor or adult

MassShootings <- mutate(MassShootings, Shooter = ifelse(age_of_shooter >= 18, "Adult", "Minor"))

# creating a new dataframe with fatalities and shooter variables

TableA <- data.frame(MassShootings$fatalities,MassShootings$Shooter)

TableA<-TableA %>% rename(fatalities=MassShootings.fatalities,shooter=MassShootings.Shooter)

table(TableA$Shooter)

# creating a box plot of fatalities by those shooters

boxplot(TableA$fatalities~TableA$shooter, las = 1, 
        ylab = "fatalities ", 
        xlab = "shooter",
        main = "fatalities by age of shooter")

# difference in sample MEANS
mean(TableA$fatalities[TableA$shooter == "Minor"]) # mean for minor

mean(TableA$fatalities[TableA$shooter == "Adult"]) # mean for Adult

# absolute difference in means
test.stat1 <- abs(mean(TableA$fatalities[TableA$shooter == "Adult"]) - 
                  mean(TableA$fatalities[TableA$shooter == "Minor"]))
test.stat1

# difference in sample MEDIANS
median(TableA$fatalities[TableA$shooter == "Adult"]) # median for Adult

median(TableA$fatalities[TableA$shooter == "Minor"]) # median for minor

# absolute difference in medians
test.stat2 <- abs(median(TableA$fatalities[TableA$shooter == "Adult"]) - 
                  median(TableA$fatalities[TableA$shooter == "Minor"]))  
test.stat2

# Permutation Test
  
# for getting the same results by using the seed
set.seed(1979)  
  
# number of observations to sample
n <- length(TableA$shooter)  
  
# number of permutation samples to take
P <- 100000
  
# the variable we will resample from 
variable <- TableA$fatalities  
  
# initialize a matrix to store the permutation data
PermSamples <- matrix(0, nrow = n, ncol = P)

  
# each column is a permutation sample of data. now, get those permutation samples, using a loop

for(i in 1:P)
  {
    PermSamples[, i] <- sample(variable, 
                               size = n, 
                               replace = FALSE)
  }
  
# first 5 of the samples of first 5 observations
PermSamples[1:5, 1:5]

# initialize vectors to store all of the Test-stats
Perm.test.stat1 <- Perm.test.stat2 <- rep(0, P)
  
# loop thru, and calculate the test-stats
for (i in 1:P)
  {
    # calculate the perm-test-stat1 and save it
    Perm.test.stat1[i] <- abs(mean(PermSamples[TableA$shooter == "Adult",i]) - 
                              mean(PermSamples[TableA$shooter == "Minor",i]))
      
    # calculate the perm-test-stat2 and save it
    Perm.test.stat2[i] <- abs(median(PermSamples[TableA$shooter == "Adult",i]) - 
                              median(PermSamples[TableA$shooter == "Minor",i]))
  }
  

# The test stats

test.stat1; test.stat2

# First 15 observations of both the test stats
round(Perm.test.stat1[1:15], 1)

round(Perm.test.stat2[1:15], 1)

# Caluculating p value
(Perm.test.stat1 >= test.stat1)[1:15]

mean((Perm.test.stat1 >= test.stat1)[1:15])

mean(Perm.test.stat1 >= test.stat1)

mean(Perm.test.stat2 >= test.stat2)

# Simulating data from distribution
set.seed(1989)
n <- 135
observed <- MassShootings$total_victims


# Create a raincloud plot to visualize the simulated data. Describe the data in terms of its number of peaks (e.g., unimodal, bimodal, etc.) and its symmetry.

tibble(value=observed) %>% ggplot(aes(y = value)) +
  ggdist::stat_halfeye(adjust = .5, width = 2*.3, .width = c(0.5, 1)) + 
    geom_boxplot(width = .3, outlier.shape = NA) +
  ggdist::stat_dots(side = "left", dotsize = 6, justification = 1.05, binwidth = .1,
                    color = "black") +
    coord_flip() +
    labs(y = "Value",
         title = "simulated values from non-normal distribution") + 
  theme_bw() +
  theme(legend.position = "none")

# What is the sample median for this data?

median(observed)

# Number of bootstrap samples
B <- 10000

# Instantiating matrix for bootstrap samples
boots <- matrix(NA, nrow = n, ncol = B)

# Sampling with replacement B times
for(b in 1:B) {
boots[, b] <- observed[sample(1:n, size = n, replace = TRUE)]
}


#  Using the generated bootstrap samples, create a bootstrap distribution of sample medians, and visualize this distribution using a histogram.

# Installing vector for bootstrap medians
bootmedians <- vector(length=B)

#sampling with replacement B times

for (b in 1:B){
  bootmedians[b]<-median(boots[,b])
}

# creating a histogram

tibble(Median=bootmedians) %>% ggplot(aes(x = Median)) +
  geom_histogram(color="white")+
  labs(y = "frequency",
       title = "Distribution of bootstrap medians") + 
  theme_bw()

# Use the bootstrap samples to obtain a nonparametric estimate of the standard error of the sample median.

SEestimate<-sd(bootmedians)
SEestimate

# Use the bootstrap samples to obtain a nonparametric 95% confidence interval for the population median.

lowerBoundMed <- quantile(bootmedians,probs=0.025)
upperBoundMed <- quantile(bootmedians,probs=0.975)

# creating a histogram

tibble(Median=bootmedians) %>% ggplot(aes(x = Median)) +
  geom_histogram(color="white")+
  geom_vline(xintercept= lowerBoundMed,size=1,
             color="dodgerblue",linetype="solid")+
  geom_vline(xintercept= upperBoundMed,size=1,
             color="dodgerblue",linetype="solid")+
  labs(y = "frequency",
       title = "Distribution of bootstrap medians") + 
  theme_bw()

# Generate B=10,000 samples each of size 30 from a normal distribution using the estimated mean and standard deviation from the observed data set, and visualize this distribution using a histogram.

B <- 10000

# Instantiating matrix for bootstrap samples
paramBoots <- matrix(NA, nrow = n, ncol = B)
xBar<-mean(observed)
s<-sd(observed)

# Sampling a normal set of n values, B times
for(b in 1:B) {
paramBoots[, b] <- rnorm(n=n,mean=xBar,sd=s)
}


# Installing vector for bootstrap medians
bootParmedians <- vector(length=B)

# caluculating median for each simulated data set
for (b in 1:B){
  bootParmedians[b]<-median(paramBoots[,b])
}


# creating a histogram

tibble(Median=bootParmedians) %>% ggplot(aes(x = Median)) +
  geom_histogram(color="white")+
  labs(y = "frequency",
       title = "Distribution of Parametric bootstrap medians") + 
  theme_bw()

# Obtain a parametric bootstrap estimate of the standard error of the sample median.

SEparestimate<-sd(bootmedians)
SEparestimate

# Obtain a parametric bootstrap 95% confidence interval for the sample median.

lowerBoundParMed <- quantile(bootParmedians,probs=0.025)
upperBoundParMed <- quantile(bootParmedians,probs=0.975)


# creating a histogram

tibble(Median=bootParmedians) %>% ggplot(aes(x = Median)) +
  geom_histogram(color="white")+
  geom_vline(xintercept= lowerBoundParMed,size=1,
             color="dodgerblue",linetype="solid")+
   geom_vline(xintercept= upperBoundParMed,size=1,
             color="dodgerblue",linetype="solid")+
  labs(y = "frequency",
       title = "Distribution of parametric bootstrap medians") + 
  theme_bw()