---
title: Supplement to Demonstrative Evidence and the Use of Algorithms in Jury Trials
short-title: Algorithms in Jury Trials # running title
author:
  - first: Rachel                 # first name
    first-init: R.              # initial of first name
    last: Rogers                  # last name
    marker: 1               # identifier for the affiliation
    email: rachel.rogers@huskers.unl.edu          # optional email address
    footnote-id: 1
  - first: Susan
    first-init: S.
    last: VanderPlas
    marker: 1
footnote:
  - content: Corresponding author
    type: corresp               # set type for the corresponding author
    id: 1
affiliation:
  - institution: University of Nebraska-Lincoln
    prefix: Department of Statistics        # optional
    country: United States of America
    marker: 1
classoption:
  - letterpaper    # or a4paper, etc.
  - inpress        # or review
# - linenumber     # add line number for submission
preamble: >
  \usepackage{amsfonts,amsmath,amssymb,amsthm}
  \usepackage{booktabs}
  \usepackage{lipsum}
  \usepackage{xr}
  \externaldocument{algs_jury_rogers}
abstract: >
  We investigate how the use of bullet comparison algorithms and demonstrative evidence may affect juror perceptions of reliability, credibility, and understanding of expert witnesses and presented evidence. The use of statistical methods in forensic science is motivated by a lack of scientific validity and error rate issues present in many forensic analysis methods. We explore what our study says about how this type of forensic evidence is perceived in the courtroom – where individuals unfamiliar with advanced statistical methods are asked to evaluate results in order to assess guilt. In the course of our initial study, we found that individuals overwhelmingly provided high Likert scale ratings in reliability, credibility, and scientificity regardless of experimental condition. This discovery of scale compression - where responses are limited to a few values on a larger scale, despite experimental manipulations - limits statistical modeling but provides opportunities for new experimental manipulations which may improve future studies in this area.
keywords:
# alphabetical; excluding anything in the title already; lower case
  - explainable machine learning
  - jury perception
bibliography: bibliography.bib
output: jds.rmd::pdf_article
---

```{r setup,echo=FALSE, message=FALSE, warning=FALSE}

# for reading in the data
library(readr) 

# for data cleaning
library(tidyr) 
library(plyr)
library(dplyr)

# for plotting/tables
library(ggplot2)
#library(ggblend)
library(patchwork)
library(gt)

# for statistical analysis
library(MASS)
library(VGAM)
library(emmeans)
library(mgcv)

knitr::opts_chunk$set(dpi = 300)

# Reading in data
merged_results_allresponses<-read_csv("data/merged_results_clean.csv")

# Removing person who took notes for on the first page that did not correspond to their scenario number
merged_results_allresponses <-merged_results_allresponses %>% subset(clean_prints!=231)

# Removing those who did not pass the attention checks
merged_results <- merged_results_allresponses %>% filter(check1 =="9mm" & check2 == "Moderately reliable") %>% dplyr::rename("Algorithm"="algorithm", "Conclusion"="conclusion", "Picture"="picture")

# Ordering the Conclusion Variable
merged_results$Conclusion <-
  factor(merged_results$Conclusion,
         levels = c("NoMatch", "Inconclusive", "Match"))

# Reading in the demographic dataset
demographics <-
  read_csv("data/demographic_clean.csv")

```


```{r data-cleaning, warning= FALSE, message= FALSE, echo= FALSE}


# Ordering the Likert scale
merged_results$firetestcred = factor(
  merged_results$firetestcred,
  levels = c(
    "Extremely noncredible",
    "Moderately noncredible",
    "Weakly noncredible",
    "Neither credible nor noncredible",
    "Weakly credible",
    "Moderately credible",
    "Extremely credible"
  )
)

# Creating line breaks in the Likert scale
levels(merged_results$firetestcred) <-
  gsub(" ", "\n", levels(merged_results$firetestcred))


# Ordering Likert scale
merged_results$firerel = factor(
  merged_results$firerel,
  levels = c(
    "Extremely unreliable",
    "Moderately unreliable",
    "Weakly unreliable",
    "Neither reliable nor unreliable",
    "Weakly reliable",
    "Moderately reliable",
    "Extremely reliable"
  )
)

# Creating line breaks
levels(merged_results$firerel) <-
  gsub(" ", "\n", levels(merged_results$firerel))


# Ordering Likert scale
merged_results$examsci = factor(
  merged_results$examsci,
  levels = c(
    "Extremely unscientific",
    "Moderately unscientific",
    "Weakly unscientific",
    "Neither scientific nor unscientific",
    "Weakly scientific",
    "Moderately scientific",
    "Extremely scientific"
  )
)

# Creating line breaks
levels(merged_results$examsci) <-
  gsub(" ", "\n", levels(merged_results$examsci))


# Creating table of algorithm understanding and conclusion
algunder_table<-table(merged_results$algunder, merged_results$Conclusion)
algunder_df<-data.frame(algunder_table, Question="Algorithm")
colnames(algunder_df)[2] <- "Conclusion"

# Creating table of examiner understanding (when the algorithm is present) and conclusion
expunder_table<-table(merged_results[merged_results$Algorithm=="Yes",]$expunder,merged_results[merged_results$Algorithm=="Yes",]$Conclusion)
expunder_df<- data.frame(expunder_table, Question="Examiner")
colnames(expunder_df)[2] <- "Conclusion"

under_df<-rbind(algunder_df, expunder_df)

# Cleaning up the understanding label
levels(under_df$Var1) <-
  gsub("<br/>", "", levels(under_df$Var1))

levels(under_df$Var1) <-
  gsub("  ", " ", levels(under_df$Var1))


levels(under_df$Var1) <- gsub(" ", "\n", levels(under_df$Var1))

# Putting gun probability and defendant probability in a single column
merged_probcombine<- merged_results %>% 
  dplyr::select(Algorithm, gunprob, Conclusion, probability) %>%
  pivot_longer(cols=c(gunprob, probability), names_to = "Type", values_to="Probability")

merged_probcombine$Type <- factor(merged_probcombine$Type, levels = c("probability", "gunprob"), 
                  labels = c("Cole", "Gun"))

```

# Ordinal Logistic Regression (Likert Scales)

Because there are not enough observations in all categories, only categories with enough observations are considered.
Most analyses are also limited to main effects.
Due to the scale compression mentioned throughout the article, this analytical approach is not recommended, as it ignores key aspects of the data collection process (such as the inclusion of the complete scale).
If there are only two categories for consideration, the response is considered as binomial.
If there are more than two categories, ordered logistic regression is first implemented using the 'VGAM' package (the 'polr' package implementation failed to find starting values in several cases), and assumptions of proportional odds are tested by comparing the likelihood to the model without the parallel odds assumption.
Unless otherwise noted, the parallel odds assumption holds.
In a few cases, there were not enough observations for the model to be computed without the parallel odds assumption.

## Credibility

### How credible did you find the testimony of Terry Smith (the firearm examiner)?

Figure \ref{fig:credible} indicates that only the top two categories of the Likert scale have enough data for a formal analysis.
Thus, only the top two categories were considered using a binomial generalized linear model.
This model does not fully consider the responses or response options given to participants, and is not recommended. 
There were not significant differences between conditions.

```{r credible, fig.cap= "Histogram of Firearms Examiner Credibility", fig.width= 7, fig.height= 4, fig.align= "center", warning= FALSE, message= FALSE, echo= FALSE}

# Credibility plot
ggplot(subset(merged_results,!is.na(firetestcred)),
         aes(x = firetestcred, fill = Conclusion)) +
  geom_bar(
    mapping = aes(y =after_stat(prop), group = Conclusion),
    position = position_dodge(preserve = "single"), color="black"
  ) +
  facet_grid(Algorithm ~ ., labeller =  purrr::partial(label_both, sep = ":\n")) +
  ggtitle("How credible did you find the testimony of Terry Smith (the firearm examiner)?") +
  scale_fill_manual(values = c("#FF8E00", "grey20", "#037AC7")) +
  theme_bw()+
  ylab("Proportion")+
  theme(axis.title.x = element_blank())


```

```{r firetestcred, echo=FALSE, warning=FALSE}

merged_results$Conclusion <- factor(merged_results$Conclusion)
merged_results$Algorithm <- factor(merged_results$Algorithm)
merged_results$Picture <- factor(merged_results$Picture)

firetestcred_subset <- filter(merged_results, 
                              firetestcred %in% c("Moderately\ncredible", "Extremely\ncredible"))

olr_firetestcred<- glm(firetestcred~Algorithm+Conclusion+Picture, 
                       family="binomial", data=firetestcred_subset)

anova(olr_firetestcred, test="Chisq")

```

### How credible did you find the testimony of Adrian Jones (the algorithm expert)?

Similar to Figure \ref{fig:credible}, Figure \ref{fig:algcred} shows that most individuals only selected the top two categories of the Likert scale.
Thus, as before, only the top two categories will be considered in statistical analysis (although this does not reflect how the data was collected).

```{r algcred, fig.cap= "Histogram of Algorithm Expert Credibility", fig.width= 7, fig.height= 4, fig.align= "center", warning= FALSE, message= FALSE, echo= FALSE}


algorithm_results <- merged_results %>% filter(Algorithm == "Yes")
algorithm_results$algtestcred = factor(
  algorithm_results$algtestcred,
  levels = c(
    "Extremely noncredible",
    "Moderately noncredible",
    "Weakly noncredible",
    "Neither credible nor noncredible",
    "Weakly credible",
    "Moderately credible",
    "Extremely credible"
  )
)

levels(algorithm_results$algtestcred) <-
  gsub(" ", "\n", levels(algorithm_results$algtestcred))


ggplot(subset(algorithm_results,!is.na(algtestcred)),
         aes(x = algtestcred, fill = Conclusion)) +
  geom_bar(mapping = aes(y = after_stat(prop), group = Conclusion),
           position = position_dodge(preserve = "single"), color="black") +
  facet_grid(Picture ~ ., labeller = label_both) +
  ggtitle("How credible did you find the testimony of Adrian Jones\n(the algorithm expert)?") +
  scale_fill_manual(values = c("#FF8E00", "grey20", "#037AC7")) +
  theme_bw()+
  ylab("Proportion")+
  theme(axis.title.x = element_blank())

```

```{r algtestcred, echo=FALSE, warning=FALSE}

algorithm_results$Conclusion <- factor(algorithm_results$Conclusion)
algorithm_results$Picture <- factor(algorithm_results$Picture)

algtestcred_subset <- filter(algorithm_results, 
                             algtestcred %in% c("Moderately\ncredible", "Extremely\ncredible"))

olr_algtestcred<- glm(algtestcred~Conclusion+Picture, family="binomial", data=algtestcred_subset)

anova(olr_algtestcred, test="Chisq")
```

## Reliability

### In general, how reliable do you think firearm evidence is?

Figure \ref{fig:reliable} has observations from each condition combination in the top three categories of the Likert scale (weakly reliable, moderately reliable, and extremely reliable), so an ordered logistic regression using the three top categories is used.

```{r reliable, fig.cap= "Histogram of perceived firearm reliability as a field", fig.width= 7, fig.height= 3, fig.align= "center", echo= FALSE, message= FALSE, warning= FALSE}

#Plotting
ggplot(subset(merged_results,!is.na(firerel)),
       aes(x = firerel, fill = Conclusion)) +
  geom_bar(
    mapping = aes(y = after_stat(prop), group = Conclusion),
    position = position_dodge(preserve = "single"), color="black"
  ) +
  facet_grid(Algorithm ~ ., labeller = label_both) +
  ggtitle("In general, how reliable do you think firearm evidence is?") +
  scale_fill_manual(values = c("#FF8E00", "grey20", "#037AC7")) +
  theme_bw()+
  ylab("Proportion")+
  theme(axis.title.x = element_blank())

```

```{r firerel, echo=FALSE, warning=FALSE}

firerel_subset <- filter(merged_results, 
                         firerel %in% c("Weakly\nreliable", "Moderately\nreliable",
                                        "Extremely\nreliable"))
 
vglm_firerel<- vglm(firerel~Conclusion+Picture+Algorithm, firerel_subset,
                                       family = cumulative(link = "logit", 
                                                           parallel = TRUE, reverse = TRUE))

coef(summary(vglm_firerel))

vglm_firerel2<- vglm(firerel~Conclusion+Picture+Algorithm, firerel_subset,
                                       family = cumulative(link = "logit", 
                                                           parallel = FALSE, reverse = TRUE))

# lrtest(vglm_firerel2,vglm_firerel)

```

### How reliable do you think the firearm evidence in this case is?

Based on Figure \ref{fig:caserel}, the top three categories contain results and are used in analysis (using ordered logistic regression).

```{r caserel, fig.cap= "Histogram of overall case reliability", fig.width= 7, fig.height= 3, fig.align= "center", echo= FALSE, message= FALSE, warning= FALSE}


df.overrel <- ddply(
  algorithm_results,
  .(Conclusion),
  summarise,
  prop = prop.table(table(overrel)),
  overrel = names(table(overrel))
)

df.overrel$overrel = factor(
  df.overrel$overrel,
  levels = c(
    "Extremely unreliable",
    "Moderately unreliable",
    "Weakly unreliable",
    "Neither reliable nor unreliable",
    "Weakly reliable",
    "Moderately reliable",
    "Extremely reliable"
  )
)

#Reorder factors (match/inconclusive/nomatch)
#gsub on axis (replace spaces with new lines)
levels(df.overrel$overrel) <-
  gsub(" ", "\n", levels(df.overrel$overrel))
ggplot(df.overrel, aes(overrel, prop, fill = Conclusion)) +
  geom_bar(stat = "identity", position = position_dodge(preserve = "single"), color="black") +
  ggtitle("How reliable do you think the firearm evidence in this case is?") +
  scale_fill_manual(values = c("#FF8E00", "grey20", "#037AC7")) +
  theme_bw()+
  ylab("Proportion")+
  theme(axis.title.x = element_blank())

```

```{r overrel, echo=FALSE, warning=FALSE}

overrel_subset <- filter(merged_results, 
                         overrel %in% c("Weakly reliable", 
                                        "Moderately reliable", "Extremely reliable"))

overrel_subset$overrel <- factor(overrel_subset$overrel, 
                                 levels=c("Weakly reliable", "Moderately reliable", 
                                          "Extremely reliable"))

vglm_overrel<- vglm(overrel~Conclusion+Picture, overrel_subset, 
                    family = cumulative(link = "logit", parallel = TRUE, reverse = TRUE))

coef(summary(vglm_overrel))

vglm_overrel2<- vglm(overrel~Conclusion+Picture, overrel_subset,
                     family = cumulative(link = "logit", parallel = FALSE, reverse = TRUE))

# lrtest(vglm_overrel2,vglm_overrel)

```

### How reliable do you think the firearms examiner's subjective opinion of the bullet comparison is, in this case?

```{r examrel, fig.cap= "Histogram of perceived firearm exam reliability", fig.width= 7, fig.height= 4, fig.align= "center", echo= FALSE, message= FALSE, warning= FALSE}

merged_results$examrel = factor(
  merged_results$examrel,
  levels = c(
    "Extremely unreliable",
    "Moderately unreliable",
    "Weakly unreliable",
    "Neither reliable nor unreliable",
    "Weakly reliable",
    "Moderately reliable",
    "Extremely reliable"
  )
)

levels(merged_results$examrel) <-
  gsub(" ", "\n", levels(merged_results$examrel))

ggplot(subset(merged_results,!is.na(examrel)),
       aes(x = examrel, fill = Algorithm)) +
  geom_bar(mapping = aes(y = ..prop.., group = Algorithm),
           position = position_dodge(preserve = "single"), color="black") +
  facet_grid(Conclusion ~ ., labeller =  purrr::partial(label_both, sep = ":\n")) +
  ggtitle("How reliable do you think the firearms examiner's subjective opinion of 
          the bullet comparison is, in this case?") +
  scale_fill_manual(values = c("grey20", "tomato")) +
  theme_bw()+
  ylab("Proportion")+
  theme(axis.title.x = element_blank())

```

In this case, there may not be enough observations in the "weakly reliable" category for analysis - so only the highest two categories of the Likert scale are analyzed.

```{r examreltest, echo=FALSE, warning=FALSE}

examrel_subset <- filter(merged_results, 
                         examrel %in% c("Weakly\nreliable", "Moderately\nreliable", 
                                        "Extremely\nreliable"))

vglm_examrel<- vglm(examrel~Conclusion+Picture+Algorithm, examrel_subset,
                                       family = cumulative(link = "logit", 
                                                           parallel = TRUE, reverse = TRUE))

coef(summary(vglm_examrel))

vglm_examrel2<- vglm(examrel~Conclusion+Picture+Algorithm, examrel_subset,
                                       family = cumulative(link = "logit", 
                                                           parallel = FALSE, reverse = TRUE))

# lrtest(vglm_examrel2,vglm_examrel)

```

### How reliable do you think the firearm algorithm evidence is, in this case?

In Figure \ref{fig:algrel}, there are enough observations in the three highest categories for ordered logistic regression.

```{r algrel, fig.cap= "Histogram of perceived algorithm reliability", fig.width= 7, fig.height= 3, fig.align= "center", echo= FALSE, message= FALSE, warning= FALSE}

df.algrel <- ddply(
  algorithm_results,
  .(Conclusion),
  summarise,
  prop = prop.table(table(algrel)),
  algrel = names(table(algrel))
)

df.algrel$algrel = factor(
  df.algrel$algrel,
  levels = c(
    "Extremely unreliable",
    "Moderately unreliable",
    "Weakly unreliable",
    "Neither reliable nor unreliable",
    "Weakly reliable",
    "Moderately reliable",
    "Extremely reliable"
  )
)

levels(df.algrel$algrel) <-
  gsub(" ", "\n", levels(df.algrel$algrel))

ggplot(df.algrel, aes(algrel, prop, fill = Conclusion)) +
  geom_bar(stat = "identity", position = position_dodge(preserve = "single"), color="black") +
  ggtitle("How reliable do you think the firearm algorithm evidence is,\n in this case?") +
  scale_fill_manual(values = c("#FF8E00", "grey20", "#037AC7")) +
  theme_bw()+
  ylab("Proportion")+
  theme(axis.title.x = element_blank())

```

```{r algreltest, echo=FALSE, warning=FALSE}

algrel_subset <- filter(algorithm_results, 
                        overrel %in% c("Weakly reliable", "Moderately reliable", 
                                       "Extremely reliable"))

algrel_subset$algrel <- factor(algrel_subset$algrel, 
                               levels=c("Weakly reliable", "Moderately reliable", 
                                        "Extremely reliable"))

vglm_algrel<- vglm(algrel~Conclusion+Picture, algrel_subset, 
                    family = cumulative(link = "logit", parallel = TRUE, reverse = TRUE))

coef(summary(vglm_algrel))

vglm_algrel2<- vglm(algrel~Conclusion+Picture, algrel_subset,
                     family = cumulative(link = "logit", parallel = FALSE, reverse = TRUE))

# lrtest(vglm_algrel2,vglm_algrel)

```

## Scientificity

### In general, how scientific do you think firearm evidence is?

The top three categories are used for analysis (Figure \ref{fig:gensci}).

```{r gensci, fig.cap= "Histogram of perceived firearm scientificity as a field", fig.width= 7, fig.height= 4, fig.align= "center", echo= FALSE, message= FALSE, warning= FALSE}

merged_results$firesci = factor(
  merged_results$firesci,
  levels = c(
    "Extremely unscientific",
    "Moderately unscientific",
    "Weakly unscientific",
    "Neither scientific nor unscientific",
    "Weakly scientific",
    "Moderately scientific",
    "Extremely scientific"
  )
)

levels(merged_results$firesci) <-
  gsub(" ", "\n", levels(merged_results$firesci))

ggplot(subset(merged_results,!is.na(firesci)),
       aes(x = firesci, fill = Algorithm)) +
  geom_bar(mapping = aes(y = ..prop.., group = Algorithm),
           position = position_dodge(preserve = "single"), color="black") +
  facet_grid(Conclusion ~ ., labeller =  purrr::partial(label_both, sep = ":\n")) +
  ggtitle("In general, how scientific do you think the firearm evidence is?") +
  scale_fill_manual(values = c("grey20", "tomato")) +
  theme_bw()+
  ylab("Proportion")+
  theme(axis.title.x = element_blank())


```

```{r firesci, echo=FALSE, warning=FALSE}

firesci_subset <- filter(merged_results, 
                         firesci %in% c("Weakly\nscientific", "Moderately\nscientific",
                                        "Extremely\nscientific"))

vglm_firesci<- vglm(firesci~Conclusion+Picture+Algorithm, firesci_subset,
                                       family = cumulative(link = "logit", 
                                                           parallel = TRUE, reverse = TRUE))

coef(summary(vglm_firesci))

vglm_firesci2<- vglm(firesci~Conclusion+Picture+Algorithm, firesci_subset,
                     family = cumulative(link = "logit", parallel = FALSE, reverse = TRUE))

# lrtest(vglm_firesci2,vglm_firesci)

```

### How scientific do you think the firearm evidence in this case is?

Shown in Figure \ref{fig:oversci}, the top three categories are considered for analysis.

```{r oversci, fig.cap= "Histogram of perceived overall scientificity in this case", fig.width= 7, fig.height= 3, fig.align= "center", echo= FALSE, message= FALSE, warning= FALSE}


algorithm_results$oversci = factor(
  algorithm_results$oversci,
  levels = c(
    "Extremely unscientific",
    "Moderately unscientific",
    "Weakly unscientific",
    "Neither scientific nor unscientific",
    "Weakly scientific",
    "Moderately scientific",
    "Extremely scientific"
  )
)

levels(algorithm_results$oversci) <-
  gsub(" ", "\n", levels(algorithm_results$oversci))

ggplot(subset(algorithm_results,!is.na(oversci)),
       aes(x = oversci, fill = Conclusion)) +
  geom_bar(mapping = aes(y = ..prop.., group = Conclusion),
           position = position_dodge(preserve = "single"), color="black") +
  facet_grid(Picture ~ ., labeller = label_both) +
  ggtitle("How scientific do you think the firearm evidence in this case is?") +
  scale_fill_manual(values = c("#FF8E00", "grey20", "#037AC7")) +
  theme_bw()+
  ylab("Proportion")+
  theme(axis.title.x = element_blank())


```

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

oversci_subset <- filter(algorithm_results, 
                         oversci %in% c("Weakly\nscientific", "Moderately \nscientific",
                                        "Extremely\nscientific"))

oversci_subset$oversci <- factor(oversci_subset$oversci, 
                                 levels=c("Weakly\nscientific", "Moderately \nscientific",
                                          "Extremely\nscientific"))

vglm_oversci<- vglm(oversci~Conclusion+Picture, oversci_subset, 
                    family = cumulative(link = "logit", parallel = TRUE, reverse = TRUE))

coef(summary(vglm_oversci))

vglm_oversci2<- vglm(oversci~Conclusion+Picture, oversci_subset,
                     family = cumulative(link = "logit", parallel = FALSE, reverse = TRUE))

# lrtest(vglm_oversci2,vglm_oversci)

```


### How scientific do you think the firearms examiner's subjective opinion of the bullet comparison is, in this case?

Figure \ref{fig:science} shows this graph of scientificity. 
There are observations in all three of the top categories, so an ordered logistic regression was used in this case as well.

```{r science, fig.cap= "Histogram of perceived scientificity of the bullet comparison of the firearm examiner", fig.width= 7, fig.height= 4, fig.align= "center", echo= FALSE, message= FALSE, warning= FALSE}

# Plotting
ggplot(subset(merged_results,!is.na(examsci)),
       aes(x = examsci, fill = Algorithm)) +
  geom_bar(
    mapping = aes(y = ..prop.., group = Algorithm),
    position = position_dodge(preserve = "single"), color="black"
  ) +
  facet_grid(Conclusion ~ ., labeller =  purrr::partial(label_both, sep = ":\n")) +
  ggtitle(
    "How scientific do you think the firearms examiner\'s subjective opinion \n of the bullet comparison evidence is, in this case?"
  ) +
  scale_fill_manual(values = c("grey20", "tomato")) +
  theme_bw()+
  ylab("Proportion")+
  theme(axis.title.x = element_blank())


```

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

examsci_subset <- filter(merged_results, 
                         examsci %in% c("Weakly\nscientific", "Moderately \nscientific", 
                                        "Extremely\nscientific"))

examsci_subset$examsci <- factor(examsci_subset$examsci, 
                                 levels=c("Weakly\nscientific", "Moderately \nscientific", 
                                          "Extremely\nscientific"))

vglm_examsci<- vglm(examsci~Conclusion+Picture+Algorithm, examsci_subset, 
                    family = cumulative(link = "logit", parallel = TRUE, reverse = TRUE))

coef(summary(vglm_examsci))

vglm_examsci2<- vglm(examsci~Conclusion+Picture+Algorithm, examsci_subset,
                     family = cumulative(link = "logit", parallel = FALSE, reverse = TRUE))

# lrtest(vglm_examsci2,vglm_examsci)

```

### How scientific do you think the firearm algorithm evidence is, in this case?

Figure \ref{fig:algsci} indicates that the top three categories may be used.

```{r algsci, fig.cap= "Histogram of perceived algorithm scientificity in this case", fig.width= 7, fig.height= 3, fig.align= "center", echo= FALSE, message= FALSE, warning= FALSE}

algorithm_results$algsci = factor(
  algorithm_results$algsci,
  levels = c(
    "Extremely unscientific",
    "Moderately unscientific",
    "Weakly unscientific",
    "Neither scientific nor unscientific",
    "Weakly scientific",
    "Moderately scientific",
    "Extremely scientific"
  )
)

levels(algorithm_results$algsci) <-
  gsub(" ", "\n", levels(algorithm_results$algsci))

ggplot(subset(algorithm_results,!is.na(algsci)),
       aes(x = algsci, fill = Conclusion)) +
  geom_bar(mapping = aes(y = ..prop.., group = Conclusion),
           position = position_dodge(preserve = "single"), color="black") +
  facet_grid(Picture ~ ., labeller = label_both) +
  ggtitle("How scientific do you think the firearm algorithm evidence is, \n in this case?") +
  scale_fill_manual(values = c("#FF8E00", "grey20", "#037AC7")) +
  theme_bw()+
  ylab("Proportion")+
  theme(axis.title.x = element_blank())


```

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

algsci_subset <- filter(algorithm_results, 
                        algsci %in% c("Weakly\nscientific", "Moderately \nscientific",
                                      "Extremely\nscientific"))

algsci_subset$algsci <- factor(algsci_subset$algsci, 
                               levels=c("Weakly\nscientific", "Moderately \nscientific",
                                        "Extremely\nscientific"))

vglm_algsci<- vglm(algsci~Conclusion+Picture, algsci_subset, 
                    family = cumulative(link = "logit", parallel = TRUE, reverse = TRUE))

coef(summary(vglm_algsci))

vglm_algsci2<- vglm(algsci~Conclusion+Picture, algsci_subset,
                     family = cumulative(link = "logit", parallel = FALSE, reverse = TRUE))

# lrtest(vglm_algsci2,vglm_algsci)

```

## Understanding

### Based on this testimony, how would you rate your understanding of the method described for the examiner's personal bullet comparison?

The top three categories are used, based on Figure \ref{fig:expunder}.

```{r expunder, fig.cap= "Histogram of understanding for the explanation of the firearms examiner", fig.width= 7, fig.height= 4, fig.align= "center", echo= FALSE, message= FALSE, warning= FALSE}

merged_results$expunder = factor(
  merged_results$expunder,
  levels = c(
    "1 <br/> I understood nothing",
    "2.0",
    "3 <br/> I understood about half of the method",
    "4.0",
    "5 <br/> I understood everything"
  )
)

levels(merged_results$expunder) <-
  gsub("<br/>", "", levels(merged_results$expunder))

levels(merged_results$expunder) <-
  gsub("  ", " ", levels(merged_results$expunder))

levels(merged_results$expunder) <-
  gsub(" ", "\n", levels(merged_results$expunder))

ggplot(subset(merged_results,!is.na(expunder)),
       aes(x = expunder, fill = Algorithm)) +
  geom_bar(mapping = aes(y = ..prop.., group = Algorithm),
           position = position_dodge(preserve = "single"), color="black") +
  facet_grid(Conclusion ~ ., labeller =  purrr::partial(label_both, sep = ":\n")) +
  ggtitle("Based on this testimony, how would you rate your understanding of the 
          method described for the examiner's personal bullet comparison?") +
  scale_fill_manual(values = c("grey20", "tomato")) +
  theme_bw()+
  theme(axis.title.x = element_blank())+ 
  ylab("Proportion")+
  scale_x_discrete(labels=c('1\nI\nunderstood\nnothing', '2', '3', '4',
                            '5\nI\nunderstood\neverything'))


```

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

expunder_subset <- filter(merged_results, expunder %in% c('3\nI\nunderstood\nabout\nhalf\nof\nthe\nmethod','4.0', '5\nI\nunderstood\neverything'))

vglm_expunder<- vglm(expunder~Conclusion+Picture+Algorithm, expunder_subset,
                                       family = cumulative(link = "logit", 
                                                           parallel = TRUE, reverse = TRUE))

coef(summary(vglm_expunder))

vglm_expunder2<- vglm(expunder~Conclusion+Picture+Algorithm, expunder_subset,
                      family = cumulative(link = "logit", parallel = FALSE, reverse = TRUE))

# lrtest(vglm_expunder2,vglm_expunder)

```


### Based on this testimony, how would you rate your understanding of the method described for the bullet matching algorithm?

Here, all but the lowest category of understanding have observations, shown in Figure \ref{fig:algunder}.

```{r algunder, fig.cap= "Histogram of understanding for the algorithm explanation", fig.width= 7, fig.height= 3, fig.align= "center", echo= FALSE, message= FALSE, warning= FALSE}


algorithm_results$algunder = factor(
  algorithm_results$algunder,
  levels = c(
    "1 <br/> I understood nothing",
    "2.0",
    "3 <br/> I understood about half of the method",
    "4.0",
    "5 <br/> I understood everything"
  )
)

levels(algorithm_results$algunder) <-
  gsub("<br/>", "", levels(algorithm_results$algunder))

levels(algorithm_results$algunder) <-
  gsub("  ", " ", levels(algorithm_results$algunder))

levels(algorithm_results$algunder) <-
  gsub(" ", "\n", levels(algorithm_results$algunder))

ggplot(subset(algorithm_results,!is.na(algunder)),
       aes(x = algunder, fill = Conclusion)) +
  geom_bar(mapping = aes(y = ..prop.., group = Conclusion),
           position = position_dodge(preserve = "single"), color="black") +
  facet_grid(Picture ~ ., labeller = label_both) +
  ggtitle("Based on this testimony, how would you rate your understanding of the 
          method described for the bullet matching algorithm?") +
  scale_fill_manual(values = c("#FF8E00", "grey20", "#037AC7")) +
  theme_bw()+
  theme(axis.title.x = element_blank())+ 
  ylab("Proportion")+
  scale_x_discrete(labels=c('1\nI\nunderstood\nnothing', '2', '3', '4',
                            '5\nI\nunderstood\neverything'))

```

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

algunder_subset <- filter(algorithm_results, algunder %in% c('2.0', '3\nI\nunderstood\nabout\nhalf\nof\nthe\nmethod', '4.0', '5\nI\nunderstood\neverything'))

algunder_subset$algunder <- factor(algunder_subset$algunder, levels=c('2.0', '3\nI\nunderstood\nabout\nhalf\nof\nthe\nmethod', '4.0', '5\nI\nunderstood\neverything'))

vglm_algunder<- vglm(algunder~Conclusion+Picture, algunder_subset, 
                    family = cumulative(link = "logit", parallel = TRUE, reverse = TRUE))

coef(summary(vglm_algunder))

vglm_algunder2<- vglm(algunder~Conclusion+Picture, algunder_subset,
                     family = cumulative(link = "logit", parallel = FALSE, reverse = TRUE))

# lrtest(vglm_algunder2,vglm_algunder)

```


## Strength

### Strength of Evidence against Cole

Figure \ref{fig:strength} shows the participants' percieved strength of evidence against the defendant.
All categories are considered in this case.

```{r strength, fig.cap= "Histogram of perceived strength of evidence against the defendant", fig.width= 7, fig.height= 4, fig.align= "center", echo= FALSE, message= FALSE, warning= FALSE}


merged_results$strength=factor(merged_results$strength, 
                               levels =c("1 <br/> Not at all strong","2.0", "3.0",
                                         "4.0","5 <br/> Moderately strong","6.0", 
                                         "7.0", "8.0", "9 <br/> Extremely strong"))

levels(merged_results$strength) <-
  gsub("<br/>", "", levels(merged_results$strength))

levels(merged_results$strength) <-
  gsub("  ", " ", levels(merged_results$strength))

levels(merged_results$strength) <-
  gsub(" ", "\n", levels(merged_results$strength))

ggplot(subset(merged_results, !is.na(strength)), aes(x=strength, fill = Algorithm))+
  geom_bar(mapping = aes(y = ..prop.., group = Algorithm), 
           position=position_dodge(preserve = "single"), color="black") +
  facet_grid(Conclusion~., labeller =  purrr::partial(label_both, sep = ":\n")) +
  ggtitle("How strong would you say the case against the defendant is?") +
  scale_fill_manual(values = c("grey20", "tomato")) +
  theme_bw()+
  ylab("Proportion")+
  theme(axis.title.x = element_blank())+ 
  scale_x_discrete(labels=c('1\nNot\nat all\nstrong', '2', '3', '4', 
                            '5\nModerately\nstrong', '6', '7', '8', '9\nExtremely\nstrong'))


```

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

vglm_strength<- vglm(strength~Conclusion+Picture+Algorithm, merged_results, 
                    family = cumulative(link = "logit", parallel = TRUE, reverse = TRUE))

coef(summary(vglm_strength))

vglm_strength2<- vglm(strength~Conclusion+Picture+Algorithm, merged_results,
                     family = cumulative(link = "logit", parallel = FALSE, reverse = TRUE))

# lrtest(vglm_strength2,vglm_strength)

```

### Strength of Evidence Against the Gun

In this case, there are not enough observations for the algorithm condition in each level for non-parallel odds (Figure \ref{fig:gunstrength}).
Therefore, only parallel odds are computed.
All levels are included in the model.

```{r gunstrength, fig.cap= "Histogram of perceived strength of evidence against the gun", fig.width= 7, fig.height= 4, fig.align= "center", echo= FALSE, message= FALSE, warning= FALSE}


merged_results$gunstrength=factor(merged_results$gunstrength, 
                                  levels =c("1 <br/> Not at all strong","2.0", "3.0",
                                            "4.0","5 <br/> Moderately strong","6.0", 
                                            "7.0", "8.0", "9 <br/> Extremely strong"))

levels(merged_results$gunstrength) <-
  gsub("<br/>", "", levels(merged_results$gunstrength))

levels(merged_results$gunstrength) <-
  gsub("  ", " ", levels(merged_results$gunstrength))

levels(merged_results$gunstrength) <-
  gsub(" ", "\n", levels(merged_results$gunstrength))

ggplot(subset(merged_results, !is.na(gunstrength)), aes(x=gunstrength, fill = Algorithm))+
  geom_bar(mapping = aes(y = ..prop.., group = Algorithm), 
           position=position_dodge(preserve = "single"), color="black") +
  facet_grid(Conclusion~., labeller = purrr::partial(label_both, sep = ":\n")) +
  ggtitle("How strong is the evidence that the defendant's gun was used to fire 
          the shot in the convenience store, in your opinion?") +
  scale_fill_manual("Algorithm", values = c("grey20", "tomato")) +
  theme_bw()+
  ylab("Proportion")+
  theme(axis.title.x = element_blank(),
        legend.position = c(1, 1), legend.justification = c(1, 1), 
        legend.direction = "horizontal", 
        legend.background = element_rect(fill = "transparent"))+
  scale_x_discrete(labels=c('1\nNot at\nall strong', '2', '3', '4', 
                            '5\nModerately\nstrong', '6', '7', '8', '9\nExtremely\nstrong'))


```

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

vglm_gunstrength<- vglm(gunstrength~Conclusion+Picture+Algorithm, merged_results, 
                    family = cumulative(link = "logit", parallel = TRUE, reverse = TRUE))

coef(summary(vglm_gunstrength))

vglm_gunstrength2<- vglm(gunstrength~Conclusion+Picture+Algorithm, merged_results,
                     family = cumulative(link = "logit", parallel = FALSE, reverse = TRUE))

# Can't compute - there are not enough algorithm observations in the non-parallel condition
# lrtest(vglm_gunstrength2,vglm_gunstrength)

```


## Mistakes

Figure \ref{fig:mistakes} shows the perceived frequency that firearms examiners make mistakes.
As can be seen, most individuals selected "Rarely", with little variations between conditions.
Very few people selected the extreme values of the scale - "Never" and "Usually".
Only the values of "Rarely", "Occasionally", and "Sometimes" are considered in the analysis.

```{r mistakes, fig.cap= "Histogram of perceived frequency of mistakes made by firearms examiners", fig.width= 7, fig.height= 3, fig.align= "center", echo= FALSE, message= FALSE, warning= FALSE}


merged_results$mistakes=factor(merged_results$mistakes, levels =c("Never","Rarely", "Occasionally","Sometimes", "Frequently", "Usually", "Always"))


ggplot(subset(merged_results, !is.na(mistakes)), aes(x=mistakes, fill = Conclusion))+
  geom_bar(mapping = aes(y = ..prop.., group = Conclusion), 
           position=position_dodge(preserve = "single"), color="black") +
  facet_grid(Picture~., labeller = label_both) +
  ggtitle("How often do firearms examiners make mistakes when determining 
          whether bullets were fired through the same gun?") +
  scale_fill_manual(values = c("#FF8E00", "grey20", "#037AC7")) +
  theme_bw()+
  ylab("Proportion")+
  theme(axis.title.x = element_blank())


```


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

mistake_subset <- filter(merged_results, mistakes %in% c('Rarely', 'Occasionally', 'Sometimes'))

mistake_subset$mistakes <- factor(mistake_subset$mistakes, levels=c('Rarely', 'Occasionally', 'Sometimes'))

vglm_mistakes<- vglm(mistakes~Conclusion+Picture+Algorithm, mistake_subset, 
                    family = cumulative(link = "logit", parallel = TRUE, reverse = TRUE))

coef(summary(vglm_mistakes))


# vglm_mistakes2<- vglm(mistakes~Conclusion+Picture+Algorithm, mistake_subset,
#                      family = cumulative(link = "logit", parallel = FALSE, reverse = TRUE))
# 
# lrtest(vglm_mistakes2,vglm_mistakes)

```


# Beta Distributions

## Probability

The graphs of probability are shown in Figure \ref{fig:prob}.
In both cases, the only significant effect is that of conclusion.
This analysis uses the 'gam' function from the 'mgcv' package.

```{r prob, fig.cap= "Probability the gun was used in the crime, or that Cole committed the crime. Black lines indicate bullet match scores for the algorithm.", out.width= "400px", fig.align= "center", echo= FALSE, eval= TRUE, warning= FALSE, message= FALSE, dev='png'}


# Putting gun probability and defendant probability in a single column
merged_probcombine<- merged_results %>% 
  dplyr::select(Algorithm, gunprob, Conclusion, probability) %>%
  pivot_longer(cols=c(gunprob, probability), names_to = "Type", values_to="Probability")

merged_probcombine$Type <- factor(merged_probcombine$Type, levels = c("probability", "gunprob"), 
                  labels = c("Cole", "Gun"))

ggplot(merged_probcombine, aes(x = Probability, fill = Algorithm)) +
  geom_density(alpha = 0.55, color=NA)+ #|> partition(vars(Algorithm)) |> blend("multiply") +
  geom_vline(data = filter(merged_probcombine, Conclusion == "NoMatch"),
             aes(xintercept = 3.4)) +
  geom_vline(data = filter(merged_probcombine, Conclusion == "Inconclusive"),
             aes(xintercept = 34)) +
  geom_vline(data = filter(merged_probcombine, Conclusion == "Match"),
             aes(xintercept = 98.9)) +
  ggtitle("Probability __ was Involved in/Committed the Crime") +
  facet_grid(Type ~ Conclusion, labeller = label_both) +
  scale_fill_manual(values = c("grey20", "tomato"))+
  theme_bw()

```

### Probability Cole Committed the Crime

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

merged_results$prob_trunc <- merged_results$probability
merged_results[merged_results$probability==0,]$prob_trunc <- 0.001
merged_results[merged_results$probability==100,]$prob_trunc <- 99.99

merged_results$sim_prob <- merged_results$prob_trunc/100

betaprob<-gam(sim_prob ~ Conclusion*Picture*Algorithm, family=betar(link="logit"),data=merged_results)

joint_tests(betaprob)
```

### Probability the Gun was Involved in the Crime

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

merged_results$gunprob_trunc <- merged_results$gunprob
merged_results[merged_results$gunprob==0,]$gunprob_trunc <- 0.001
merged_results[merged_results$gunprob==100,]$gunprob_trunc <- 99.99

merged_results$gun_prob <- merged_results$gunprob_trunc/100

betagunprob<-gam(gun_prob ~ Conclusion*Picture*Algorithm, family=betar(link="logit"),data=merged_results)


joint_tests(betagunprob)
```

# Binomial Responses

## Do you think guns leave unique markings on discharged bullets/casings?

Responses were recorded as in a yes/no format.
`r sum(merged_results$unique=="No")` individuals indicated that they did not think think firearms left unique markings, out of `r dim(merged_results)[1]` total responses.

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

merged_results$unique_num <- NA
merged_results[merged_results$unique=="No",]$unique_num <- 0
merged_results[merged_results$unique=="Yes",]$unique_num <- 1

olr_unique<- glm(unique_num~Conclusion+Picture+Algorithm, family="binomial", data=merged_results)

anova(olr_unique, test="Chisq")
```


## Conviction

Individuals were given the following question: 
"The State has the burden of proving beyond a reasonable doubt that the defendant is the person who committed the alleged crime. 
If you are not convinced beyond a reasonable doubt that the defendant is the person who committed the alleged crime, you must find the defendant not guilty.
Would you convict this defendant, based on the evidence that you have heard?".
Results are shown in Figure \ref{fig:probguilt}.

```{r probguilt, fig.cap= "Probabilities based on whether the participants thought the defendant was guilty", fig.width= 7, fig.height= 5, fig.align= "center", echo= FALSE, eval= TRUE, warning= FALSE, message= FALSE}

probplot <-
  ggplot(merged_results, aes(x = Conclusion, y = probability, fill = guilty)) +
  geom_point(
    position = position_jitterdodge(
      jitter.width = 0.2,
      jitter.height = 0.4,
      dodge.width = 1
    ),
    size = 0.5
  ) +
  geom_boxplot(position = position_dodge(1),
               alpha = 0.5,
               outlier.shape = NA) +
  ggtitle("Probability Cole\nCommited the Crime") +
  scale_fill_manual(values = c("grey20", "plum1"), name="Convict?")+
  theme_bw()+
  theme(legend.position = "none")

gunplot <-
  ggplot(merged_results, aes(x = Conclusion, y = gunprob, fill = guilty)) +
  geom_point(
    position = position_jitterdodge(
      jitter.width = 0.2,
      jitter.height = 0.4,
      dodge.width = 1
    ),
    size = 0.5
  ) +
  geom_boxplot(position = position_dodge(1),
               alpha = 0.5,
               outlier.shape = NA) +
  ggtitle("Probability Gun\nwas used in the Crime") +
  ylab("Probability")+
  scale_fill_manual(values = c("grey20", "plum1"), name="Convict?")+
  theme_bw()

probplot + gunplot

```

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

merged_results$guilt_num <- NA
merged_results[merged_results$guilty=="No",]$guilt_num <- 0
merged_results[merged_results$guilty=="Yes",]$guilt_num <- 1

olr_guilt<- glm(guilt_num~Conclusion*Picture*Algorithm, family="binomial", data=merged_results)

anova(olr_guilt, test="Chisq")

```
