Assignment.Rmd

---
title: "Assignment 3 - Thames Valle 2020-2022"
author: "u5538165"
date: "`r Sys.Date()`"
output:
  pdf_document:
    extra_dependencies: ["float"]
geometry: margin = 25mm
fontsize: 11pt 
---
Word Count: 2986
```{r setup, include=FALSE, }
knitr::opts_chunk$set(echo = TRUE,
                      fig.pos = "H",
                      out.extra = "",
                      fig.width = 5,
                      fig.height = 3,
                      out.width = "70%")
library(dplyr)
library(lubridate)
library(tidyr)
library(knitr)
library(glmnetUtils)
library(gt)
library(MASS)
library(car)
library(pROC)
library(glmnet)
library(ggplot2)
load("C:/Users/ajbev/OneDrive/Desktop/Year3/ST404/assignment 3/StopSearchThamesV.Rdata")
data <- StopSearchThamesV
```

```{r cleaning, include = FALSE}
# make a season column and time column
data <- data %>%
  mutate(
    season = case_when(
      month(Date) %in% c(12,1,2) ~ "Winter",
      month(Date) %in% c(3,4,5) ~ "Spring",
      month(Date) %in% c(6,7,8) ~ "Summer",
      month(Date) %in% c(9,10,11) ~ "Autumn"
    ),
    time = format(Date, "%H:%M:%S")
  )
#exclude date as now not needed
data$Date <- NULL
#looking at the categories
unique(data$Outcome)
unique(data$Legislat)
unique(data$SearchPurpose)

#convert all blank data entries as NA
data$AgeRange[data$AgeRange == ""] <- NA
data$Outcome[data$Outcome== ""] <- NA

colSums(is.na(data))
#2249 missing in age range, 2806 missing in outcome, 881 missing in sex
# we omit as the total data set has 110087 rows so omitting is not impactful on the data
data <- na.omit(data)

nrow(data)
# make the column entries cleaner and create new time of day column
data$male <- ifelse(data$male == "Male",1,0)
data$female <- 1 - data$male
data <- data[data$AgeRange != "",]
data$AgeRange <- droplevels(data$AgeRange)
data$hour <- as.integer(substr(data$time, 1,2))
data$time <- cut(data$hour,
                 breaks = c(-1,5,10,16,20,23),
                 labels=c("Night", "Morning", "Afternoon", "Evening", "Late Evening"))
data$hour <- NULL

data <- data[data$Year %in% c(2020,2021,2022),]

data2 <- data
                
```
# Abstract

In this report I aim to explore the relationship between factors that are recorded in police stop and search data. This also includes building predictive models to give probabilities of an individual being Female (for example) given that you know all other variables. These variables include: Type of Search, Age, Outcome of Search, Sex, Ethnicity, Legislation, Search Purpose, Season of Year, Time of Day. I used various model selection and variable selection techniques in order to ensure all parameters in the model were justified to be in there. Cross validation - splitting data into subsets to train and test - was also implemented which ensured that I was not over fitting my model (including too many unnecessary parameters). Throughout the entire report it was clear that the biggest factor in both questions was ethnicity. In female stop and searches, we see that the odds of an individual being White is much higher than any other ethnicity given all else constant. In arrests, we see searches of Black individuals have a 22% higher odds of arrest than White individuals given all else constant. Another important thing to mention is the prominence of Drugs in this data set and hence in the UK. Despite not being a standout in the arrest rate, we see Drugs as a major reason for the rate of stop and searches in women as seen in Table 2 and Figure 1 and also in the UK in general due to the large count. Interestingly, we also see a positive association between age and arrest rate - over 34's have a 51% higher odds of being arrested than 18 - 24 year old's when everything else is constant. This could be down to laws/police being more relaxed with younger people or that older individuals tend to be involved in more severe cases.

# Question 1: Sex differences

### Summary

In the question we are told that the probability that an individual being in the data set was a woman is around 8% (national average). In this current subset of data that this project is focusing on, the percentage of women is 10.5 % so slightly higher. In the next table, it will help us understand better the proportions of women being stopped and searched for each variable and its levels.  

```{r 1, echo = FALSE}

# a function that gets the percentage of females for that variable. Inputs the variable and its name.
get_pct_female <- function(var, varname) {
  #creates a cross tab table
  tab_counts <- table(var, data$male)
  # then the row proportions
  tab_prop <- prop.table(tab_counts, margin = 1)
  
  n_female <- tab_counts[, 1]  #0=female
  n_male   <- tab_counts[, 2]  #1=male
  # then build the output data frame
  data.frame(
    Variable  = varname,
    Category  = rownames(tab_counts),
    N         = as.integer(n_female) + as.integer(n_male),
    Pct_Female = round(tab_prop[, 1] * 100, 1),
    row.names = NULL
  )
}

#then using the function above we see the percentage of women for each category of each variable rounded
table1 <- bind_rows(
  get_pct_female(data$Type,             "Type of Search"),
  get_pct_female(data$AgeRange,         "Age Range"),
  get_pct_female(data$OfficerDefEthnic, "Officer-Defined Ethnicity"),
  get_pct_female(data$Legislat,         "Legislation"),
  get_pct_female(data$time, "Time of Day"),
  get_pct_female(data$SearchPurpose,    "Search Purpose"),
  get_pct_female(data$season,           "Season")
)

colnames(table1) <- c("Variable", "Category", "Count", "% Female")

# for presenting purposes we use kable.
kable(table1, 
      caption = "Female Percentage by variable",
      col.names = c("Variable", "Category", "Count", "% Female"),
      digits = 1, 
      booktabs = TRUE,
      row.names = FALSE)

```

In Table 1, it shows all variables mentioned in the question alongside their individual levels. Some key areas to point out: The percentage of women being stopped and searched does seem to increase slightly as the age range increases. We can also see that in the "Asian" level, women are stopped and searched (2.9%) much less than in the "Black" level (6.2%) along with "White" women (14.0%). In the "White" ethnicity, the percentage of women being stopped is very high at 14.0%. We can identify categories that are very rare in women being stopped and searched such as "CriminalJPOAct1994" (0.0%), "SearchPurpose:Other" and "Weapons"(8.6%) - so anything other than drugs and theft. And also as mentioned earlier Asian and Black women. 

We can also identify that there is not a strong correlation between any of the categories in Season with Summer (10.9% highest) and Spring (9.9% lowest). In terms of time of day we see a rise in the proportion of women being stopped and searched in the morning and a dip in the afternoon. The other grouped times are all very similar. 

Initially the data set had just a date column in the format "Year - Month - Date HH:MM:SS". I created two new columns extracting the hour of the day (24 hour clock) and the month of the year. I then grouped the months into the 4 traditional seasons in order to reduce the number of parameters in the model and also to be able to identify any association between daylight/temperature with the percentage of women stopped and searched. To deal with the hour column, I grouped them into Morning, Afternoon, Evening, Late Evening and Night (as seen in table 1). With the hope to identify any patterns with the time of day and proportion of women stopped and searched along with reducing the number of parameters compared to including all 24 hours separately. Grouping the time of day this way was intuitive to me because you would expect there to be more stop and searches in the afternoon/late evenings and less in the mornings (due to the amount of people in public during these times). 
```{r , echo = FALSE, warning = FALSE, fig.cap="Stacked bar chart of percentage of stop and searches being female."}

stacked_var_data <- table1 %>%
  mutate(
    Female_Count = Count * (`% Female` / 100)
  ) %>%
  group_by(Variable) %>%
  mutate(
    Category_pct_females = Female_Count / sum(Female_Count) * 100
  ) %>% 
  ungroup()

ggplot(stacked_var_data, aes(x = Variable,
                             y = Category_pct_females,
                             fill = Category)) + 
  geom_col(width = 0.75)+
  geom_text(
    aes(label = ifelse(Category_pct_females >= 5, Category, "")),position = position_stack(vjust = 0.5), size = 3, colour = "black" )+
  labs(title = "Distribution of Female cases within each variable", x = "Variable", y = "Percentage of Female cases", fill = "Category") +
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none")

```

Figure 1 shows the distribution of female counts across the categories of each stop and search variable. Within age range, for example, we can see there to be approximately an even spread throughout however in ethnicity there is a large proportion of this taken up by counts recorded as white. This plot enables us to see where female cases are concentrated in the data set. 


```{r ,fig.pos="H", echo = FALSE, fig.cap="Line graph of the percentage of females being stopped by time of day and year", fig.height=4, fig.width=5}

data %>%
  group_by(time,season) %>%
  summarise(pct_female = mean(female, na.rm = TRUE) * 100,
            n = n(), .groups = "drop") %>%
  ggplot(aes(x = time, y = pct_female, group = season, colour = season)) +
  geom_line() +
  labs(title = "% female by time of day and year(season)",
       x = "Time of day",
       y = "% female", 
       colour = "Season") + 
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))


```

From the plot on time of day and year we can see that in the spring there is a clear dip in the proportion of stop and searches being female. It also allows us to see that in the Autumn and Winter there is a big relative spike in the morning for the stop and search percentage in women. As mentioned earlier, this plot allows us to identify the dip in the afternoon, alongisde allowing us to see how the rate changes with time and season together.  


 

### Modelling 


Moving onto model selection, as we want a model that gives accurate predictions and we only have 7 variables in the model, we will start with AIC step wise regression. Under these circumstances, step wise is useful because we have a lot of observations and a relatively small number of possible models to consider. After performing both directions of AIC step wise regression we see that this has kept all variables apart from search type. 

When using only step wise, selected coefficients are larger in absolute value than the true population values as the selection chooses variables with a strong association by chance. Also, the standard errors do not account for the uncertainty from the selection process itself which leads to falsely narrow confidence intervals.  

***Note: The reference categories are the most common categories. This is important to remember for model interpretation. "18-24", "White","Drugs", "Legislat:Other","Summer","Night".***

```{r 3, include = FALSE, echo = FALSE}
# relevel for the most common categories
data$AgeRange <- relevel(factor(data$AgeRange), ref = "18-24")
data$OfficerDefEthnic <- relevel(factor(data$OfficerDefEthnic), ref = "White")
data$SearchPurpose <- relevel(factor(data$SearchPurpose), ref = "Drugs")
data$Legislat <- relevel(factor(data$Legislat), ref = "Other")
data$Type <- relevel(factor(data$Type), ref = "PersonVehicle")
data$season <- relevel(factor(data$season), ref = "Summer")
data$time <- relevel(factor(data$time), ref = "Night")


# null
model_null <- glm(female ~ 1, data = data, family = binomial)


# full
model_full <- glm(female ~ Type + time + season + Legislat + SearchPurpose + AgeRange + OfficerDefEthnic, data = data, family = binomial)

model_step <- stepAIC(model_full, direction = "both", trace = TRUE)
```

After performing stepwise AIC with both directions we see that this method removed Type of search as a variable reducing the AIC to 25830. Hence the final model as of now includes time of day, season, legislation, search purpose, age, ethnicity. To make my work consistent I will perform ridge and elastic net regression with type removed.

Looking at the model summary for the AIC optimised model, it shows that in the age range 25 - 34, this does not differ from the baseline of 18 - 24 and because these age groups are of similar lifestyles I will combine them. The main point of concern and limitation in the model is the Legislation variable. Looking further into the excel document about the "other" column it contains a wide range of severity including public stores act or sporting events act to aviation security act and customs and excise management act. Because of this I will collapse the Drugs,FireArms and CriminalJPO Acts into Other as they are not significant in our model. Despite there being no women stop and searches under the criminal JPO act, there were only 51 cases of this so to improve the model I decided it best to collapse this. We also see in the time variable that the counts for Morning, Evening and Late Evening are low compared to other and Afternoon (hence the error is high). Because of this and their non significance I will collapse these into a new Other category for time. Finally due to the small variation in the season variable, I will combine Winter and Autumn into one, also because these two seasons are similar in terms of temperature and daylight hours. 

To improve accuracy as stated in the question and to stop potential overfitting, I applied regularised logistic regression techniques using Elastic net (alpha from zero to 0.5). I looped over these alpha values because I want the most accurate model and LASSO would remove variables entirely.


```{r elastic net, echo = FALSE, fig.cap="Cross validation curve for Elastic net, alpha = 0.5"}

data$AgeRange2 <- as.character(data$AgeRange)
data$AgeRange2[data$AgeRange2 %in% c("18-24", "25-34")] <- "18-34"
data$AgeRange2 <- as.factor(data$AgeRange2)


data$Legislat2 <- as.character(data$Legislat)
data$Legislat2[data$Legislat2 %in% c("CriminalJPOAct1994", "FirearmsAct1968","DrugsAct1971","Other")] <- "Other"
data$Legislat2 <- as.factor(data$Legislat2)

data$time2 <- as.character(data$time)
data$time2[data$time2 %in% 
  c("Morning", "Evening", "Late Evening")] <- "Other"
data$time2 <- as.factor(data$time2)

data$season2 <- as.character(data$season)
data$season2[data$season2 %in% 
  c("Winter", "Autumn")] <- "Winter/Autumn"
data$season2 <- as.factor(data$season2)


levels(data$Legislat2)[levels(data$Legislat2) == "EvidenceAct1984"] <- "Evidence"

X <- model.matrix(female ~ time2 + season2 + Legislat2 + SearchPurpose + AgeRange2 + OfficerDefEthnic, data = data)[, -1]

y <- data$female

lambdastotry <- 10^seq(-5,2,length.out = 50)
alphastotry <- seq(0,0.5,by=0.1)

elasticNetcv <- cva.glmnet(y =y, x = X, alpha = alphastotry, lambda = lambdastotry, nfolds = 10, family = "binomial")

plot(elasticNetcv, main = "Alpha Tuning for Cross Validation")



```

From the figure we can see there is not much difference in minimal binomial deviance between differing alpha values. As we want the most accurate model and also because in stepwise all but the time variable was kept, we will use ridge regression. 

After choosing ridge, we now choose a lambda. Often, I would choose lambda.1se as it is simpler and still gives accurate predictions however as the question states most accurate then I will choose lambda.min which is $0.00052$. So this method has confirmed to us that the predictors in our GLM are stable and not overfitted as our coefficients are very similar (small lambda) . 

```{r ridge,message=FALSE, echo = FALSE, fig.cap = "ROC Curve for our model after ridge regression."}

fmodel <- elasticNetcv$modlist[[1]]

opt_lambda <- fmodel$lambda.min

pred_probs <- predict(fmodel, newx = X, s = opt_lambda, type = "response")

roc <- roc(y,as.vector(pred_probs))

auc(roc)

plot(roc, main = paste("ROC Curve - Model AUC =",round(auc(roc),3)))

fmodel_glm2 <- glm(female ~ time2 + season2 + Legislat2 + SearchPurpose + AgeRange2 + OfficerDefEthnic, data = data, family = binomial)

```
This AUC value suggests the model correctly distinguishes between male and female stop and searches 64.8% of the time. This is an acceptable value however reasons for it not being good could be down to the small percentage of females in the data or the lack of variables that account for other contexts.

### Model summary

As we are using logistic regression and our outcome is binary we are working on the log scale. In order to present odds ratios we need to exponentiate them. The confidence intervals were calculated through the "confint" function in R directly applied to the model. 

```{r table, echo= FALSE, message = FALSE}

OR_confint <- exp(cbind(OR = coef(fmodel_glm2), confint(fmodel_glm2)))

glm2_summary <- summary(fmodel_glm2)$coefficients
pvals <- glm2_summary[, 4]

summary_table <- data.frame(
  OR = round(OR_confint[,1],3),
  Lower = round(OR_confint[,2],3),
  Upper = round(OR_confint[,3], 3),
  p_value = ifelse(pvals < 0.001, "<0.001", round(pvals,3))
) %>%
  mutate(CI = paste0("[", Lower, ",", Upper,"]"))

summary_table <- summary_table[-1,]

kable(summary_table, caption = "Model Summary for Female Stop and Search ")

```
### Ethnicity 
Here we can see that with the reference being White, the odds that a stop and search of a female being Black,Asian or Other is around 61%,82% and 32% less odds compared to people that are White respectively. 


### Search Purpose 
Again, with the reference level being Drugs, we can see that the odds that a female's search purpose is Other,Weapons or Theft is 37%,27% and 21% less odds compared to those who have a search purpose of drugs. 

### Table of probabiltiies

```{r probtable, echo = FALSE}

pred_table <- expand.grid(
  time2 = c("Night","Other"),
  season2 = "Summer",
  Legislat2 = "Other",
  SearchPurpose = c("Drugs"),
  AgeRange2 = c("10-17","over 34"),
  OfficerDefEthnic = c("White", "Asian","Black","Other")
)

predictions <- predict(fmodel_glm2, newdata = pred_table, type = "link", se.fit = TRUE)

pred_table$FemaleProbability <- round(plogis(predictions$fit), 4)

pred_table$SE <- round(plogis(predictions$fit) * (1 - plogis(predictions$fit)) * predictions$se.fit, 4)

pred_table$MaleProbability <- round(1 - plogis(predictions$fit), 4)

kable(pred_table,caption = "Prediction Table for Sex", col.names = c("Time","Season", "Legislation","Search","Age","Ethnicity","F_Prob","SE","M_Prob"))

```

In the table of probabilities shown above, I decided to contrast Time of day to the most different time ( Night and Other which includes Morning, Evening and Late evening). I also decided to vary the Age to the most contrasting categories whilst keeping Summer,Other and Drugs fixed for Season, Legislation and Search Purpose (most common). In the table we see given that the variables are Night - Summer - Other - Drugs - over 34 - White, the probability that a stop and search is female is 0.198 whereas if we just change the ethnicity to Asian the probability drops to 0.043. 

 



# Question 2: Ethnic differences

### Summary

```{r q2, echo = FALSE}
# filter to just person searches and in the specific years

data2$arrest <- ifelse(data2$Outcome %in% c("Suspect arrested", "Arrest"), 1, 0)

# object of search, time of day, season, legislation, sex, age, ethnicity      

```


We can calculate that the percentage of cases where there was an arrest is 12.5% (5022 out of 40084). This is below the national average of 14% and I found this by creating a new "Arrest" indicator variable from the "Outcome" variables that had value "suspect arrested" and "arrest". 

```{r summary, echo = FALSE}

get_pct_arrest <- function(v, name){
  tab_counts <- table(v, data2$arrest)
  tab_prop <- prop.table(tab_counts, margin = 1)
  
  N_noarrests <- tab_counts[,1]
  N_arrests <- tab_counts[,2]
  
  data.frame(
    Variable = name,
    Category = rownames(tab_counts),
    Total_Count = as.integer(N_noarrests) + as.integer(N_arrests),
    `Pct_Arrests` = round(tab_prop[,2] *100, 2), row.names = NULL
  )
}

table2 <- bind_rows(
  get_pct_arrest(data2$OfficerDefEthnic, "Officer defined Ethnicity"),
  get_pct_arrest(data2$AgeRange,         "Age range"),
  get_pct_arrest(data2$SearchPurpose,    "Search purpose"),
  get_pct_arrest(data2$Legislat,         "Legislation"),
  get_pct_arrest(data2$male,             "Male"),
  get_pct_arrest(data2$season,           "Season"),
  get_pct_arrest(data2$time,             "Time")
)

kable(table2,
      caption =  "Arrest percentage within each category of each variable",
      col.names = c("Variable","Category","Total Count","Pct Arrested"))

```

Looking at the table , we can clearly see an increasing trend in arrests with age range. The earlier age ranges have a lower arrest rate than the older age ranges. Intuitively, this is what you expect as older individuals may be stopped with stronger grounds for suspicion and hence arrest whereas younger people may be stopped more speculatively. We can see Stolen/Theft search purpose has a higher arrest rate, this could be because theft is easily spotted and searches are based on stronger evidence whereas drug searches for example are more speculative.In Sex there is a slight difference between male and females, which is slightly higher in men, as expected. The seasons seem to have an evenly weighted arrest rate, with Summer being slightly lower. In the evening and late evening we see a lower arrest rate than any other time (peaking in the morning). We see a higher arrest rate in the evidence act which is similar to that of Firearms however the firearms count is very low in proportion to our dataset. Same as before, CriminalJPOAct has a very low count which we will take into account further on.


```{r diag, echo = FALSE, fig.cap="Bar chart of percentage of arrests per ethnic category."}

ethical_data <- as.data.frame(prop.table(table(data2$OfficerDefEthnic, data2$arrest), margin = 1))

names(ethical_data) <- c("Ethnicity", "Arrest", "proportion")

ethical_data <- ethical_data[ethical_data$Arrest==1,]

ggplot(ethical_data, aes(x = reorder(Ethnicity,-proportion), y = proportion * 100, fill = Ethnicity)) + geom_bar(stat = "identity") + labs(title = "Arrest rate by ethnicity", subtitle = "Thames Valley person searches in 2020-2022 period", x = "Ethnicity", y = "% arrested", fill = "") + geom_text(aes(label = paste0(round(proportion *100, 1), "%")))



```

Figure 4 shows the arrest rate by officer defined ethnicity, where the overall rate is 12%. The Black ethnicity has the highest arrest rate at 14.3% and Asian has the lowest arrest rate at 11.2%.  Comparing this to question 1, the differences between ethnicities are a lot less significant. However the figure allows us to see clearly the slight pattern between arrest rate and ethnicity, despite not including details about any other variable. 

### Modelling

As question 2 does not explicitly state that the prediction accuracy is the most important here we will focus on understanding the ethnic differences in arrest rates. Hence, I will begin by using a stepwise BIC in order to reduce model complexity. The same problems with stepwise as in question 1 will also apply here as well (I wont repeat them). 

```{r, echo = FALSE, include = FALSE}

full_model <- glm(arrest ~ SearchPurpose + time + season + Legislat + male + AgeRange + OfficerDefEthnic, data = data2, family = "binomial")

null_model <- glm(arrest ~ 1, data = data2, family = "binomial")


data2$Legislat <- relevel(factor(data2$Legislat), ref = "Other")
data2$OfficerDefEthnic <- relevel(factor(data2$OfficerDefEthnic), ref = "White")
data2$AgeRange <- relevel(factor(data2$AgeRange), ref = "18-24")
data2$SearchPurpose <- relevel(factor(data2$SearchPurpose), ref = "Drugs")
data2$male <- relevel(factor(data2$male), ref = "1")
data2$time <- relevel(factor(data2$time), ref = "Night")



model_bic_step <- stepAIC(full_model, direction = "both", k = log(nrow(data2)))

```
After performing this method it has dropped the season variable which makes sense given how evenly weighted the percentages were for each season. In the model output the non significant categories are morning, FirearmsAct and Ethnicity:Other. I will not produce the current model summary due to spacing and tidiness. I will now combine Night and Morning due to their similar rates and I will collapse CriminalJPO and FirearmsAct (due to low counts) into Other for the same reasons as the previous question.

I decided to use Lasso regression after the step wise BIC because my goal is looking into ethnic difference between my outcome variable. The lasso penalty performs a variable selection by shrinking irrelevant coefficients to zero producing a more parsimonious model. This will help us as a simpler model that still has good predictive power is a lot easier to interpret for readers. So we do the same as before however this time we do not tune for alpha. 

```{r lasso, echo = FALSE, fig.cap="Lasso CV Plot"}
data2$time2 <- as.character(data2$time)
data2$time2[data2$time2 %in% c("Morning", "Night")] <- "Night/Morning"
data2$time2 <- as.factor(data2$time2)

data2$Legislat2 <- as.character(data2$Legislat)
data2$Legislat2[data2$Legislat2 %in% 
  c("FirearmsAct1968", "Other","CriminalJPOAct1994")] <- "Other"
data2$Legislat2 <- as.factor(data2$Legislat2)

levels(data2$SearchPurpose)[levels(data2$SearchPurpose) == "StolenTheft"] <- "Theft"

levels(data2$Legislat2)[levels(data2$Legislat2) == "EvidenceAct1984"] <- "Evidence"

glm_model2 <- glm(arrest ~ SearchPurpose + time2 + Legislat2 + male + AgeRange + OfficerDefEthnic, data = data2, family = "binomial")


X2 <- model.matrix(arrest~SearchPurpose + time2 + Legislat2 + male + AgeRange + OfficerDefEthnic, data = data2)[,-1]

y2 <- data2$arrest

lasso_cv <- cv.glmnet(X2, y2, alpha = 1, lambda = lambdastotry, nfolds = 10, family = "binomial")

par(mar = c(5,4,6,2))
plot(lasso_cv, main = "")
title("Lasso CV Plot", line = 4)

lambda <- lasso_cv$lambda.1se

lasso_table <- coef(lasso_cv, s = lambda)

```
This figure shows how the binomial deviance changes as lambda is decreased. I will choose the lambda.1se (0.0072 the left dotted line) as the difference between the two lambda values is very small. Additionally, we also want a simpler model.

```{r comparison table, echo = FALSE}

glm_sum <- summary(glm_model2)$coefficients
glm_sf <- ifelse(glm_sum[-1,4] < 0.05, "Significant","Not")
glm_direct <- ifelse(glm_sum[-1,1] > 0 ,"+","-")


lasso_coef <- as.matrix(coef(lasso_cv, s = lambda))
lasso_stat <- ifelse(lasso_coef[-1] ==0, "Zeroed","Kept")
lasso_direct <- ifelse(lasso_coef[-1]>0,"+","-")

consis <- ifelse(
  (glm_sf == "Significant"& lasso_stat == "Kept" & glm_direct == lasso_direct) | (glm_sf == "Not" & lasso_stat == "Zeroed"),"Yes","No"
)

compare_table <- data.frame(
  GLM_SF = glm_sf,
  GLM_Direction = glm_direct,
  Lasso_Decision = lasso_stat,
  Lasso_Direction = ifelse(lasso_stat == "Zeroed","-",lasso_direct),
  Consistent = consis,
  row.names = rownames(glm_sum[-1,])
)

kable(compare_table, caption = "Comparison table for the Stepwise BIC and Lasso techniques.")

```
In the table, its clear that most of the predictors were consistent across our GLM and Lasso method (for variable selection). However we did have 3 variables that were zeroed by Lasso but seen as significant in the GLM. 

As the Night/Morning level has nearly half of the data set and it is significant in our GLM I will not remove it from the model and the same goes for the EvidenceAct, however it is important to note this as a limitation. 

The biggest limitation here is the Search Purpose. The other category has 246 observations which is very small, however I will not remove or collapse it into another level as this does not make sense due to the range of search purposes "other" could contain, so it is important to remember this. We essentially have three main search purposes: Drugs, Weapons and Theft. 

### Model summary


```{r summary2, echo = FALSE, message = FALSE}

# Get GLM summary for p-values
glm_sum2 <- summary(glm_model2)$coefficients

# Get ORs and CIs
odds_conf <- exp(cbind(OR = coef(glm_model2), confint(glm_model2)))

# Combine
results_table <- data.frame(
  OR      = round(odds_conf[-1, 1], 3),
  Lower  = round(odds_conf[-1, 2], 3),
  Upper = round(odds_conf[-1, 3], 3),
  P_value = ifelse(glm_sum[-1, 4] < 0.001, "<0.001", 
                   round(glm_sum[-1, 4], 3))
) %>%
  mutate(CI = paste0("[", Lower, ", ", Upper, "]")) 

kable(results_table,
             col.names = c("OR", "Lower","Upper", "P-value", "CI"),
             caption   = "Summary Results for Arrests")


```

***Note: Again the reference levels are Search Purpose Drugs, Afternoon, Drugs Act, male, 18-24, White. Our interpretations are for when every other variable is fixed.***

### Ethnicity
We see that the odds of a stop and search leading to an arrest in Black, Asian, and other ethnicities are 22% more likely, 16% less likely and 6.4% more likely compared to that of White ethnicity. 

### Age
The odds of a stop and search leading to an arrest in the age ranges 10-17, 25-34 and over 34 is 32% lower odds, 46% higher odds and 51% higher odds compared to those in the age range 18-24. 

### Table of probabilities

```{r predictions, echo = FALSE}


prediction_table <- expand.grid(
  time2 = c("Evening", "Afternoon"),
  Legislat2 = "Other",
  SearchPurpose = "Drugs",
  AgeRange = c("18-24","over 34"),
  OfficerDefEthnic = c("White","Other","Asian","Black"),
  male = c(0,1))

prediction_table$male <- factor(prediction_table$male, levels = c("0","1"))

preds <- predict(glm_model2, newdata = prediction_table, type = "link", se.fit = TRUE)

prediction_table$Arrest_Prob <- round(plogis(preds$fit), 4)

prediction_table$SE <- round(plogis(preds$fit) *(1-plogis(preds$fit)) * preds$se.fit, 4)

kable(prediction_table, caption = "Prediction table for Arrests", col.names = c("Time", "Legislat", "Search","Age","Ethnicity", "Male", "Arrest_Prob", "SE"))


```

Table 7 is rather long however it was important to vary the time of day and also the age range alongside the difference between male and female arrest probabilities. The tables shown so far have always focused on Drugs and 18-24 and over 34 years old. This is because Drugs is the most common reason for people being stopped and searched and also arrested in this data set. Focusing on Ethnicity, it is interesting to see the difference between Asian and Black. In Asian women we see an extremely low arrest rate compared to that of Black women and men in particular. With the highest arrest probability being Black men, over 34, in the afternoon with a search purpose of drugs.  


# Conclusion and Discussion

In conclusion to go back over the two questions asked for the report we first were asked to make a model to provide accurate predictions for the probability an observation was female given we know the rest of the variables. Second, we were asked to create a new indicator variable for arrests made as a result of the stop and searches and to investigate how the probability of arrest vary as certain variables are changed - most importantly for ethnicity. 

In the first question, ethnicity was the most significant important variable, as seen before, individuals being recorded as Asian, Black and Other had a much lower odds of being female compared to White individuals. Which tells us that female stop and searches are concentrated on White individuals. In the second question, individuals recorded as Black has a 22% higher odds of arrest compared to White individuals after keeping everything the same, whilst Asian people had 16% lower odds. To strengthen these findings they were consistent across our GLM and regularised models. 

Drugs is the most dominant search purpose across the data set in terms of pure count however drug searches are the lowest to lead to arrest and also they are most associated with females in the first question. This could suggest that drug searches are most speculative meaning officers will search somebody with little actual evidence to do so rather than for a search purpose such as theft which is more obvious to see. 

As mentioned, there is a clear positive association between age and arrest probability with people over 34 having 51% higher odds of arrest than those in 18-24, keeping everything else the same. This finding would need further investigation however one reason could be that younger people are more likely to receive cautions not arrests. Also the odds of arrest in males is slightly higher than that of females which is consistent with what you would expect however the magnitude of the effect suggests sex is not a key variable when predicting arrest rates. 

A further analysis that would interest me is to cover different more zoomed in areas of the UK where I can have access to variables of the city/town such that deprivation levels, quality of life, demographic composition. Then compare these results with different cities to look for similarities and differences. This would also explain the results better as for example in extremely deprived areas you would expect the crime rate to be high which would impact the stop and search count/rate and the arrest count/rate in certain ethnicities that live there. This would then explain why certain ethnicities have a higher arrest rate/stop and search rate. 

### Limitations

Staring with the limitations of our data, it only includes one police force - Thames Valley - which may not correspond similar to other police forces (different policing styles or local environments). Officer Defined Ethnicity has room for human error, we are relying on officers to correctly identify someones ethnicity which gives our data uncertainty. As mentioned and explained in the conclusion, the data does not include area characteristics. Stepwise is known to produce unstable results as the method selects variables greedily without considering all combinations. Despite checking our model variable selection further through penalised regression, this is still a limitation in our models. Again with step wise we also have the risk of overfitting, producing narrow confidence intervals. The decision to collapse certain levels of variables was up to me under the guidance of significance testing and ridge/lasso regression however this is still subjective and different decisions will lead to different models. Despite using cross validation to select our lambda, the final model coefficients were based on the data set. To reduce bias, future models should be validated on a held out sample. 


# Appendix

In terms of resources, I mainly used the lecture notes from ST404 and the computer practical solutions for ideas or ways to select a model/select variables. This was useful however it was difficult to keep going back through notes. If I needed questions answered that were not Statistics related such as what certain legislations mean I would use my search engine google chrome. If I had trouble when knitting a the rmd file despite the code running and could not fix it easily, I would ask AI to fix the error - usually it would be to do with R packages or the YAML section. I would also ask for help when I had overlapping lables in figures.

References:

Prabhakaran, S. (2016) *Model Selection in R - AIC, BIC, and Stepwise Methods*. 
r-statistics.co. Available at: https://r-statistics.co/Model-Selection-in-R.html

Kassambara, A. (2018) *Penalized Logistic Regression Essentials in R: Ridge, Lasso 
and Elastic Net*. STHDA. Available at: https://www.sthda.com/english/articles/36-classification-methods-essentials/149-penalized-logistic-regression-essentials-in-r-ridge-lasso-and-elastic-net/

***I will not repeat any code that was used twice but with columns and values changed such as for the prediction tables.***
```{r appendixcode, eval = FALSE, echo=TRUE}
# relevel for the most common categories so this # is now our baseline levels
data$AgeRange <- relevel(factor(data$AgeRange), ref = "18-24")
data$OfficerDefEthnic <- relevel(factor(data$OfficerDefEthnic), ref = "White")
data$SearchPurpose <- relevel(factor(data$SearchPurpose), ref = "Drugs")
data$Legislat <- relevel(factor(data$Legislat), ref = "Other")
data$Type <- relevel(factor(data$Type), ref = "PersonVehicle")
data$season <- relevel(factor(data$season), ref = "Summer")
data$time <- relevel(factor(data$time), ref = "Night")

# null
model_null <- glm(female ~ 1, data = data, family = binomial)

# full
model_full <- glm(female ~ Type + time + season + 
Legislat + SearchPurpose + AgeRange + OfficerDefEthnic, data = data, family = binomial)
#stepwise aic
model_step <- stepAIC(model_full, direction = "both", trace = TRUE)


# create design matrix and prediction column to

#implement cross validation as seen in the lecture notes
X <- model.matrix(female ~ time2 + season2 + Legislat2 + 
SearchPurpose + AgeRange2 + OfficerDefEthnic, data = data)[, -1]
y <- data$female
# give a range of lambdas
lambdastotry <- 10^seq(-5,2,length.out = 50)
# tune for the right alpha
alphastotry <- seq(0,0.5,by=0.1)
# Seen in the notes, elastic net , alpha tuning.
elasticNetcv <- cva.glmnet(y =y, x = X, alpha = alphastotry, 
                           
lambda = lambdastotry, nfolds = 10, family = "binomial")
# finally plot the different 
# alpha curves against binomial deviance.
plot(elasticNetcv, main = "Alpha Tuning for Cross Validation")


# calculate the odds ratios and confidence intervals through 
#exponentiating the coefficients and then using the confint function.
OR_confint <- exp(cbind(OR = coef(fmodel_glm2), confint(fmodel_glm2)))
# extract the coefficients from the model summary and also the p values.
glm2_summary <- summary(fmodel_glm2)$coefficients
pvals <- glm2_summary[, 4]
# Create the summary table as a data frame and create the columns 
#and fill them in from the values we just calculated.
summary_table <- data.frame(
  OR = round(OR_confint[,1],3),
  Lower = round(OR_confint[,2],3),
  Upper = round(OR_confint[,3], 3),
  p_value = ifelse(pvals < 0.001, "<0.001", round(pvals,3))
) %>%
  mutate(CI = paste0("[", Lower, ",", Upper,"]"))
# delete the intercept for tidiness and also not very interpretable
summary_table <- summary_table[-1,]
#use kable to give a readable output table.
kable(summary_table, caption = "Model Summary for Female Stop and Search ")


# similar to before create a grid that we can fill values into 
#with the chosen contrasting variables.
pred_table <- expand.grid(
  time2 = c("Night","Other"),
  season2 = "Summer",
  Legislat2 = "Other",
  SearchPurpose = c("Drugs"),
  AgeRange2 = c("10-17","over 34"),
  OfficerDefEthnic = c("White", "Asian","Black","Other")
)
# use the predict function from our glm to get predictions and the SE
predictions <- predict(fmodel_glm2, newdata = pred_table, type = "link", se.fit = TRUE)

#As our predictions are on the log odds scale we convert using the logistic function.
pred_table$FemaleProbability <- round(plogis(predictions$fit), 4)

#We then use the delta method to convert to standard error on the correct scale, not log.
pred_table$SE <- round(plogis(predictions$fit) *
(1 - plogis(predictions$fit)) * predictions$se.fit, 4)

#simple to calculate the male probability
pred_table$MaleProbability <- round(1 - plogis(predictions$fit), 4)

# seen before use kable and caption and rename column for presentation.
kable(pred_table,
      caption = "Prediction Table for Sex",
      col.names = c("Time","Season", "Legislation",

              "Search","Age","Ethnicity","F_Prob","SE","M_Prob"))

```

```{r appendixb, eval = FALSE, echo=TRUE}

# I created the arrest binary variable by looking through all levels in outcome and filtering the two levels seen below.
data2$arrest <- ifelse(data2$Outcome %in% c("Suspect arrested", "Arrest"), 1, 0)

# very similar to AIC however you need to specify the k value being equal to log(n).
model_bic_step <- stepAIC(full_model, direction = "both", k = log(nrow(data2)))

#without repeating very similar code to before, we now change alpha to be 1 to use LASSO and 
lasso_cv <- cv.glmnet(X2, y2, alpha = 1, lambda = lambdastotry, 
nfolds = 10, family = "binomial")

#Then plot the figure to show how the deviance changes with lambda.
plot(lasso_cv, main = "Lasso CV Plot")


#extract the glm coefficients
glm_sum <- summary(glm_model2)$coefficients
# Create a vector saying if each level is significant or not.
glm_sf <- ifelse(glm_sum[-1,4] < 0.05, "Significant","Not")
# give the direction of the coefficient - negative or positive to check whether they agree
glm_direct <- ifelse(glm_sum[-1,1] > 0 ,"+","-")

# same with lasso obtain vectors with the coefficients and
# whether or not LASSO zeroed out the level.
lasso_coef <- as.matrix(coef(lasso_cv, s = lambda))
lasso_stat <- ifelse(lasso_coef[-1] ==0, "Zeroed","Kept")
#coefficient direction negative or positive
lasso_direct <- ifelse(lasso_coef[-1]>0,"+","-")

#create a rule for whether or not LASSO and GLM agree or not. 
consis <- ifelse(
  (glm_sf == "Significant"& lasso_stat == "Kept" & glm_direct == lasso_direct) 
| (glm_sf == "Not" & lasso_stat == "Zeroed"),"Yes","No"
)

#very similar to past table constructions however catch when LASSO zeros out
#a level this just gives it a dash for direction (it has no direction).
compare_table <- data.frame(
  GLM_SF = glm_sf,
  GLM_Direction = glm_direct,
  Lasso_Decision = lasso_stat,
  Lasso_Direction = ifelse(lasso_stat == "Zeroed","-",lasso_direct),
  Consistent = consis,
  row.names = rownames(glm_sum[-1,])
)

#kable used again for presentation.
kable(compare_table, caption = "Comparison table for the Stepwise BIC and Lasso techniques.")

#Rest of the code is very similar to part A.


```