--- title: "Reproducing the Website Analyses" author: "" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Reproducing the Website Analyses} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", message = FALSE ) ``` ```{r libraries} library(dplyr) library(ggplot2) library(knitr) library(tidyr) library(salem) ``` The [Salem Witchcraft Website] (https://www2.tulane.edu/~salem/index.html) presents a number of data analyses and explains how to reproduce them using Excel. This vignette will document how to as closely as possible reproduce those analyses using R. In some cases the analyses will be extended in order to illustrate points in the website text that were not easily expressible in Excel. While the website has seven data sets, the package has only 5. This is because R makes it relatively simple to manipulate the data and also because this organization is more in line with how social scientists would likely create a data set today. As is always true in R, there are many different ways to do the same thing. If you are using this for teaching purposes you will likely want to use the approach that most closely matches what you are doing in your course and what your goals are for students to work on their own analyses. The focus of this document is on use of _ggplot2_ for graphing and in subsequent iterations will also illustrate use of the base R plotting system. There are some examples that use the tidyverse for analysis and some that use base R. However, even within these major groupings there are many ways to acheive the same goals. We have aimed to present ways to do this that are relatively accessible to people who have not used R extensively. Many of the analyses presented are graphical and so that will be the focus here. However some preprocessing will also be illustrated. This version does not yet include spatial analysis, but that is on the roadmap. For the ggplots this theme is being used in may cases. ```{r} salem_theme <- theme(axis.text.x = element_text(angle = 45, vjust = .5)) ``` ## Accusations and Executions The website includes displays of the distribution by month. #### Histogram and frequency table of accusations by month ```{r fig.width=5} ggplot(data = accused_witches) + aes(x = Month.of.Accusation.Name) + geom_bar() + labs(title = "Accusations by month", x = "Month") + scale_x_discrete(drop = FALSE) + salem_theme table(accused_witches$Month.of.Accusation.Name, useNA = "ifany", dnn = "Month") |> kable(caption = "Accusations by Month") ``` #### Cumulative accusations by month ```{r fig.width=5} ggplot(data = accused_witches) + aes(x = as.numeric(Month.of.Accusation.Name)) + stat_ecdf(na.rm = TRUE, geom = "line") + labs(title = "Cumulative accusations by month", x = "Month", y = "Cumulative Accusations") + scale_x_continuous(breaks = c(1:12), labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul","Aug", "Sept", "Oct", "Nov", "Dec") ) cumsum(table(accused_witches$Month.of.Accusation.Name)) |> kable() ``` An alternative approach using dplyr that provides some other options. This approach would also let you add other variables to the new monthly data, whether the number of executions or other data you collect. One advantage of this approach is that it is easier to create a Pareto Histogram as displayed on the website. These examples show three possible ways to display the Pareto histogram. As usual there are other approaches. ```{r message=FALSE, fig.width=5} monthly_accusations <- accused_witches |> filter(!is.na(Month.of.Accusation)) |> group_by(Month.of.Accusation.Name) |> summarize(number = n()) |> mutate(cumulative = cumsum(number) ) monthly_accusations |> kable(caption = "Accusations by Month") ggplot(data = monthly_accusations) + aes(x = Month.of.Accusation.Name, y = cumulative, group =1) + geom_col(aes(y = number)) + geom_line() + labs(title = "Frequency and Cumulative Frequency of Accusations by Month") + salem_theme monthly_accusations |> arrange(desc(number)) |> mutate(Month.of.Accusation.Name = factor(Month.of.Accusation.Name, levels=Month.of.Accusation.Name)) |> mutate(pareto_cumulative = cumsum(number) ) |> ggplot() + aes(x = Month.of.Accusation.Name, y = number) + geom_col() + labs(title = "Pareto Histogram of Accusations by Month") + salem_theme monthly_accusations |> arrange(desc(number)) |> mutate(Month.of.Accusation.Name = factor(Month.of.Accusation.Name, levels=Month.of.Accusation.Name)) |> mutate(pareto_cumulative = cumsum(number) ) |> ggplot() + aes(x = Month.of.Accusation.Name, y = pareto_cumulative) + geom_point( ) + geom_path(aes(y=pareto_cumulative, group=1), colour="blue", lty=3, linewidth= 0.9) + labs(title = "Pareto Cumulative Distribution of Accusations by Month") + salem_theme monthly_accusations |> arrange(desc(number)) |> mutate(Month.of.Accusation.Name = factor(Month.of.Accusation.Name, levels = Month.of.Accusation.Name)) |> mutate(pareto_cumulative = cumsum(number) ) |> ggplot() + aes(x = Month.of.Accusation.Name, y = number) + geom_col() + geom_path(aes(y=pareto_cumulative, group = 1), colour = "blue", lty = 3, linewidth = 0.9) + labs(title = "Pareto Histogram of Accusations by Month") + salem_theme ``` ### Executions Creating the bar chart of executions over time is similar to creating that for accusations. Additionally, these are shown as what is termed scatterplots but now would be called line graphs. Since most accused were not executed (and this is indicated by setting them to NA or missing) we'll remove the missing cases from the graph. ```{r fig.width=5} ggplot(data = na.omit(accused_witches)) + aes(x = Month.of.Execution.Name) + geom_bar() + labs(title = "Executions by month", x = "Month") + scale_x_discrete(drop = FALSE) + salem_theme table(accused_witches$Month.of.Execution.Name, useNA = "ifany", dnn = "Month") |> kable(caption = "Executions by Month") ``` There are many ways to explore the relationship between executions and accusations. On the website this is displayed using overlapping line plots of the number of executions in a month and the number of accusations. If teaching, it is important to point out that the execuations in a month were not generally of the same people who were accused in that month. To reproduce this, one approach is to again use summarize() to create a monthly data set. ```{r fig.width=6} monthly_accused <- accused_witches |> filter(!is.na(Month.of.Accusation)) |> group_by(Month.of.Accusation.Name) |> summarize(accusations = n()) |> rename(Month = Month.of.Accusation.Name) monthly_executed <- accused_witches |> filter(!is.na(Month.of.Execution)) |> group_by(Month.of.Execution.Name) |> summarize(executions = n()) |> rename(Month = Month.of.Execution.Name) monthly_data <- left_join(monthly_accused, monthly_executed) |> mutate(executions = ifelse(is.na(executions), 0, executions)) ggplot(monthly_data) + aes(x = Month, y = accusations, group = 1) + geom_path( color = "blue" ) + geom_line(aes(y = executions, x = Month, group = 1), color = "red") + labs(title = "Executions and Accusations by Month", subtitle = "Accusations in blue and executions in red") + salem_theme ``` ### Geographic Spread of Accusations The pivot table of month of accusation by town. However that table is huge and only the first seven columns are shown on the website. A more compact display is given below rather than using kable(). Also note that the row and column totals are added. ```{r} addmargins(table( accused_witches$Month.of.Accusation.Name, accused_witches$Residence, dnn = c("Month", "Residence") )) ``` The site discusses the "grand total" data (the row sum) as reflecting the varying intensity in the different towns. This can also be shown using a univariate table (or dplyr::summarize()) ```{r} table(accused_witches$Residence, dnn = "Residence") |> kable(caption = "Residences of Accused") accused_witches |> group_by(Residence) |> summarize(number = n()) |> kable(caption = "Residences of Accused") ``` The website contains a towns data set that essentially reproduces the large contingency table. This makes sense to do in Excel, but in R we do not need to makea separate data set. To make a histogram for each town we can select the town we want. This is illustrated below for Andover. Note that in the data set the white space from the beginning and ending of the place names has been removed. ```{r fig.width=5} accused_witches |> filter(Residence == "Andover") |> ggplot() + aes(x = Month.of.Accusation.Name) + geom_bar() + scale_x_discrete(drop=FALSE) + salem_theme ``` The website encourages examining the period before and after the convening of the Court of Oyer and Terminer in June. To help make this more visible we can add a vline. ```{r fig.width=5} accused_witches |> filter(Residence == "Andover") |> ggplot() + aes(x = Month.of.Accusation.Name) + geom_bar() + geom_vline(xintercept = 6.5, color = "red") + scale_x_discrete(drop = FALSE) + labs(title = "Accusations by Month", subtitle = "Red line indicates reconvening of the Court of Oyer and Terminer", x = "Month of accusation") + salem_theme ``` Alternatively, we can use faceting to show all of the histograms. ```{r fig.width=6, fig.height=8} accused_witches |> ggplot() + aes(x = Month.of.Accusation.Name) + geom_bar() + geom_vline(xintercept = 6.5, color = "red") + scale_x_discrete(drop = FALSE) + facet_wrap(facets = vars(Residence), ncol = 4) + labs(title = "Number of accusations by month", subtitle = "Red line indicates reconvening of the Court of Oyer and Terminer", x = "Month of Accusation") + theme(axis.text.x = element_text(angle = 90, vjust = .5)) ``` ### Maps We can use maps to represent the geographic distribution over time spatially. The salem_region data set is an sf (simple features) data set that provides geographic data for the towns in the three Massachuessets counties represented in our data: Essex, Middlesex and Sussex. To use this data you will need to install the sf package. The code at the start of the second code chunk tests for whether you have done this. Not all of the towns in these three counties had accusations: the majority did not. The data also include the total accused for each town and the monthly accused for the towns that ever had an accusation. However, these variables should _never_ be used outside of maps because the results will be incorrect. This is because many towns have multiple records. This is especially true if they include islands, which is common in this region of the state. These are based on modern municipal boundaries which are close to but not the same as the historical ones. There are many geographic analyses that can be done. Below a comparison of the distribution of whether towns had at least on accusation in each of four months is given. First, create variables indicating whether towns did or did not have an accusation in a month and remove the TOWN_LABELs for all the towns that never had an accusation. ```{r} newdata <- salem_region |> mutate(February.Any = February > 0, March.Any = March > 0, April.Any = April > 0, May.Any = May > 0, June.Any = June > 0, July.Any = July > 0, August.Any = August > 0, September.Any = September >0, October.Any = October > 0, November.Any = November > 0) newdata$TOWN_LABEL <- ifelse(newdata$n_accused == 0, NA, newdata$TOWN_LABEL) ``` ```{r} ``` Here `requireNamespace` is used to check whether the `sf` package is installed. That is because `sf` is only suggested and not required. It can be complex to install on some computers. ```{r eval=FALSE, fig.width=8, fig.height=8} if (requireNamespace("sf", quietly = TRUE)) { library(sf) p1 <- ggplot(newdata) p2 <- geom_sf_text(aes(label = TOWN_LABEL), color = "blue", size = 2, nudge_x = 5, nudge_y = 5, na.rm = TRUE) p3 <- scale_fill_manual(values = c( "grey", "red"), na.value = "white") } ``` ```{r eval=FALSE,fig.width=8, fig.height=8} if (requireNamespace("sf", quietly = TRUE)) { p1 + geom_sf(data = newdata, aes(geometry = newdata$geometry, fill = February.Any ), color = "black", size = .1) + p3 + p2 + labs(title = "Location of Accusations in February") } if (requireNamespace("sf", quietly = TRUE)) { library(sf) p1 + geom_sf(data = newdata, aes(geometry = geometry, fill = July.Any), color = "black", size = .1) + p3 + p2 + labs(title = "Location of Accusations in July") } ``` ```{r eval=FALSE,fig.width=8, fig.height=8} if (requireNamespace("sf", quietly = TRUE)) { p1 + geom_sf(data = newdata, aes(fill = August.Any, geometry = geometry), color = "black", size = .1) + p3 + p2 + labs(title = "Location of Accusations in August") } if (requireNamespace("sf", quietly = TRUE)) { p1 + geom_sf(data = newdata, aes(fill = November.Any, geometry = geometry), color = "black", size = .1) + p3 + p2 + labs(title = "Location of Accusations in November") } ``` ## Social Conflict The website explores a variety of dimensions of social conflict in the Salem region. One of the important social cleavages was between those who did or did not support Reverend Parris, who had been appointed as the minister for the new church in Salem Village. (This was the Puritan church, which was the official church in Massacheusetts.) The site has two data sets, one of pro-Parris individuals and one of anti-Parris individuals. For use in R these have been combined with an added variable identifying the position of the indivual. This data set is called parris_social. Also the original column label for the anti-Parris data was "[Church] Member" while for the pro-Parris data it was "Church-Member." In the parris_social data Church-Member is used. Using the variable `view` it can be subset to pro- and anti- groups. However it is not necessary to do this. A three-way table can be created, as one example. Two ways of doing this are shown. Marginals and proportions could be added. ```{r} table(parris_social$Sex, parris_social$Identification, parris_social$view, dnn = c("Sex", "Identification", "View")) ftable(parris_social$Sex, parris_social$Identification, parris_social$view, dnn = c("Sex", "Identification", "View")) ``` Treating position on Parris as the dependent varable, we might rearrange the tables. ```{r} table(parris_social$view, parris_social$Identification, parris_social$Sex, dnn = c("View", "Identification", "Sex")) ftable( parris_social$Identification, parris_social$view, parris_social$Sex, dnn = c("View", "Identification", "Sex")) ``` The parris_social data set focuses on those who signed pro- and anti-Parris petitions. In contrast, the salem_village data set includes individuals who did not sign petitions in addition to those who did. However, it is limited to those who paid taxes so many of the young men and women of all ages who are in the parris_social data are not included. This also includes data from multiple years. The basic table for this is the analysis of church membership and position on Parris. ```{r} table(salem_village$Petition, salem_village$Church.to.1696, dnn = c( "Membership", "View")) addmargins( table(salem_village$Petition, salem_village$Church.to.1696, dnn = c( "Membership", "View"))) ``` ### Political Power As explained on the website, the Village of Salem was governed by a Committee of Five in addition to the Town Meeting [see Salem Village Politics](https://www2.tulane.edu/~salem/Politics%20Intro.html). The data set that is part of this package called committee_list is structured differently than the original data. It is in "long" format with one row per person-year for each committee member for the years they served. This can be coverted to various other formats if desired. You can make a table see how the make up of the committee changed over the period from 1687 to 1698. ```{r} table(committee_list$Petition, committee_list$Year, dnn = c("Petition", "Year")) ``` The website shows a complex pivot table involving three variables. ```{r} table( committee_list$Petition, committee_list$Social, committee_list$Year ) ``` One way to explore the association is to cross tabulate the Petition and Social variables. However, this counts the same people multiple times if they served for multiple years. ```{r} table( committee_list$Petition, committee_list$Social ) ``` To adjust for this and get the data matching the website we need to remove those duplicates. Then the table matches that found at [Committee Social Report](https://www2.tulane.edu/~salem/Committeemen%20Social%20Report.html). ```{r} committee_list2 <- committee_list |> group_by(Committee.Members) |> summarize(Social = first(Social), Petition = first(Petition), Terms = n()) addmargins(table( committee_list2$Petition, committee_list2$Social )) ``` ### Wealth and Salem Village Conflict The tax_comparison data set is also set up in long format. This can then be used to create a comparison table similar to the table on [Tax Comparison](https://www2.tulane.edu/~salem/1690%20Tax%20Comparisons.html). One limitation is that R does not have a built in function for the mode. This is the most complex analysis on the website. Complications include people moving in and out of the data in different years (individuals appear between one and six times), the changing amount of taxes over the time period and name changes. The analysis looks at both group and individual level comparisons and individual level mobility. Here is the distribution of tax payers by year. ```{r} table(tax_comparison$Year, dnn = c("Year")) |> kable(caption = "Tax payers by year") ``` The next table summarizes descriptive statistics for all of the years. The website general shows one year at a time or in some cases compares two years. You can add `filter(Year = 1681)` to display only the 1681 data, changing the filter to select whatever years you want. For example to select all of the data for 1681 through 1694 one way would be to add `select(Year >= 1691 & Year <= 1694)`. Let's start by looking at the group level data. The first table is for everyone in a year. The second separates the data based on which petition the tax payer signed (or if they) did not sign. ```{r} tax_comparison |> group_by(Year) |> summarize(mean = round(mean(Tax), 1), median = round(median(Tax), 1), minimum = min(Tax), maximum = max(Tax), range = maximum - minimum, sum = sum(Tax), Count = n()) |> kable(caption = "Taxes paid by year") tax_comparison |> group_by(Year, Petition) |> summarize(mean = round(mean(Tax), 1), median = round(median(Tax), 1), minimum = min(Tax), maximum = max(Tax), range = maximum - minimum, sum = sum(Tax), Tax_payers = n()) |> kable(caption = "Taxes paid by year and petition signed") ``` To do the analysis comparing the median tax paid for the three groups we can graph the values over time. In this case, we will first make a data set of statistics to work with. ```{r fig.width=5} yearly_tax <- tax_comparison |> group_by(Year, Petition) |> summarize(mean = round(mean(Tax), 1), median = round(median(Tax), 1), minimum = min(Tax), maximum = max(Tax), range = maximum - minimum, sum = sum(Tax), Count = n()) ggplot(yearly_tax) + aes(x = Year, y = median, group = Petition, color = Petition) + geom_line() + labs(title = "Tax paid by petition", y = "Shillings", x = "Year") ``` The analysis of annual data on the website specifically focuses on thse who were Pro- and Anti- Parris based on petition signatures. As can be seen in the previous graph, those who did not sign petititions always had a lower median tax, but the other two groups were inconsistent. To address this the site presents the ratio of the medians of the two group. A ratio of 1.0 means that the median tax for the two groups was identical. A ratio of great than 1 means that the Pro-Parris group had a higher median. A ratio of less than one means that the Anti-Parris group had a higher median. ```{r fig.width=5} ratios <- yearly_tax |> filter(Petition %in% c("Pro-P", "Anti-P")) |> select(Year, Petition, median) |> pivot_wider(id_cols = c(Year), values_from = median, names_from = Petition) |> mutate(median_ratio = `Pro-P`/`Anti-P`) ratios |> kable() ggplot(ratios, aes(x = Year, y = median_ratio)) + geom_point() + geom_line(group = 1) + geom_hline(yintercept = 1, color = "blue") + labs(title = "Ratio of Pro- to Anti- Parris median tax", y = "Ratio", x = "Year") ``` The same analysis can be done but substituting mean for median. ```{r} ratios <- yearly_tax |> filter(Petition %in% c("Pro-P", "Anti-P")) |> select(Year, Petition, mean) |> pivot_wider(id_cols = c(Year), values_from = mean, names_from = Petition) |> mutate(mean_ratio = `Pro-P`/`Anti-P`) ratios |> kable() ggplot(ratios, aes(x = Year, y = mean_ratio)) + geom_point() + geom_line(group = 1) + geom_hline(yintercept = 1, color = "blue") + labs(title = "Ratio of Pro- to Anti- Parris mean tax", y = "Ratio", x = "Year") ``` The site also looks at individual trajectories using ranks and percentiles within each year, allowing analysis of whether an individual's financial position changed substantially over time. Individuals whose percentile changed by more than 10 are considered to have moved substantially. There are many ways to do this, and the method chosen will depend on your purposes. The example below is not necessarily the best or most efficient. To do the analysis we need to obtain the rank and percentile for each of the tax amounts within the years. Finally, we need to compare those. In this example we will focus on the [1681-1690 comparison shown](https://www2.tulane.edu/~salem/1681-90%20Rank.html) while providing a framework for the others. First create a wide data set and calculate ranks and cumulative proportions. This allows viewing and sorting of individuals for specific years. ```{r} persisters <- tax_comparison |> pivot_wider(id_cols = c(Name, Petition), names_from = Year, names_prefix = "X", values_from = Tax ) |> relocate(X1681, .before = X1690) ``` Next we want to compare our two selected years and determine the number who went up, down and stayed the same and how that relates to the petitions. This analysis is based on the empirical cumulative relative frequency. To reproduce the analysis on the website we use the dplyr min_rank function (which is equivalent to rank(ties.method = "min") in base R). This ranks all of the tax payments with gaps between the ties. So if two people are tied for the second rank, the person after them would be ranked fourth. Turn the rank into a percentile by subtracting 1 and dividing by the number of tax payers in that year. ```{r} comparison <- persisters |> mutate(across(starts_with("X"), list( rank = min_rank ) )) |> filter(!is.na(X1690) & !is.na(X1681)) |> mutate( X1690_ptile = 100 * (X1690_rank - 1)/100, X1681_ptile = 100* (X1681_rank - 1)/94, compare_1681_1690 = case_when( X1690_ptile -X1681_ptile >= 10 ~ "Up", X1690_ptile -X1681_ptile <= -10 ~ "Down", TRUE ~ "No Change" ) ) table(comparison$compare_1681_1690, comparison$Petition, dnn = c("Change", "Petition")) |> kable(caption = "Change in position 1681-1690") ``` ## Conclusion This has illustrated how to reproduce the analyses from the Salem Witchcraft Website. The raw data is available at the websie and at [the respository for this package](https://github.com/elinw/salem/tree/master/rawdata). The code used to transform the original data into the files presented here is also available at that location.