diff --git a/README.md b/README.md index 6959271..fc7fac6 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ -# Theology, Religious Studies, and Data Science +# Hacking Religion -An open textbook introducing data science to religious studies. +An open textbook introducing data science to religious studies (or vice versa) You can view a live demo of the book here: [https://kidwellj.github.io/hacking_religion_textbook/intro.html] @@ -25,6 +25,10 @@ Directory structure includes: 3. change to the `hacking_religion` subdirectory and run `quarto preview` 4. alternatively you can render a copy of the book using `quarto render` +# Cookbook + +There is a companion repository which contains recipes which will replicate the example data used in the book. This can be found here: [https://github.com/kidwellj/hacking_religion_cookbook] + # Copyright Content here, unless otherwise indicated are copyright by Jeremy H. Kidwell. Please re-use them as they are covered by Creative Commons Attribution 4.0 International Licence (CC BY 4.0). diff --git a/docs/Hacking-Religion--TRS---Data-Science-in-Action.pdf b/docs/Hacking-Religion--TRS---Data-Science-in-Action.pdf index b9bee7a..199fd32 100644 Binary files a/docs/Hacking-Religion--TRS---Data-Science-in-Action.pdf and b/docs/Hacking-Religion--TRS---Data-Science-in-Action.pdf differ diff --git a/docs/chapter_1.html b/docs/chapter_1.html index 7f22509..d6bd28b 100644 --- a/docs/chapter_1.html +++ b/docs/chapter_1.html @@ -252,37 +252,31 @@ div.csl-indent {

In the example below, we’re going to read in data from a comma separated value file (“csv”) which has rows of information on separate lines in a text file with each column separated by a comma. This is one of the standard plain text file formats. R has a function you can use to import this efficiently called “read.csv”. Each line of code in R usually starts with the object, and then follows with instructions on what we’re going to put inside it, where that comes from, and how to format it:

setwd("/Users/kidwellj/gits/hacking_religion_textbook/hacking_religion")
-library(here) # much better way to manage working paths in R across multiple instances
-
-
here() starts at /Users/kidwellj/gits/hacking_religion_textbook
-
-
library(tidyverse)
-
-
-- Attaching core tidyverse packages ------------------------ tidyverse 2.0.0 --
-v dplyr     1.1.3     v readr     2.1.4
-v forcats   1.0.0     v stringr   1.5.0
-v ggplot2   3.4.3     v tibble    3.2.1
-v lubridate 1.9.3     v tidyr     1.3.0
-v purrr     1.0.2
-
-
-
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
-x dplyr::filter() masks stats::filter()
-x dplyr::lag()    masks stats::lag()
-i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
-
-
here::i_am("chapter_1.qmd")
+library(here) |> suppressPackageStartupMessages() +library(tidyverse) |> suppressPackageStartupMessages() +here::i_am("chapter_1.qmd")
here() starts at /Users/kidwellj/gits/hacking_religion_textbook/hacking_religion
-
uk_census_2021_religion <- read.csv(here("example_data", "census2021-ts030-rgn.csv"))
+
# Set up local workspace:
+if (dir.exists("data") == FALSE) {
+  dir.create("data") 
+}
+if (dir.exists("figures") == FALSE) {
+  dir.create("figures") 
+}
+if (dir.exists("derivedData") == FALSE) {
+  dir.create("derivedData")
+}
+
+uk_census_2021_religion <- read.csv(here("example_data", "census2021-ts030-rgn.csv"))

2.2 Examining data:

What’s in the table? You can take a quick look at either the top of the data frame, or the bottom using one of the following commands:

-
head(uk_census_2021_religion)
+
head(uk_census_2021_religion)
                 geography   total no_religion christian buddhist  hindu jewish
 1               North East 2647012     1058122   1343948     7026  10924   4389
@@ -302,7 +296,7 @@ i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all

This is actually a fairly ugly table, so I’ll use an R tool called kable to give you prettier tables in the future, like this:

-
knitr::kable(head(uk_census_2021_religion))
+
knitr::kable(head(uk_census_2021_religion))
@@ -418,7 +412,7 @@ i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all

You can see how I’ve nested the previous command inside the kable command. For reference, in some cases when you’re working with really complex scripts with many different libraries and functions, they may end up with functions that have the same name. You can specify the library where the function is meant to come from by preceding it with :: as we’ve done knitr:: above. The same kind of output can be gotten using tail:

-
knitr::kable(tail(uk_census_2021_religion))
+
knitr::kable(tail(uk_census_2021_religion))
@@ -546,7 +540,7 @@ i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all

The first thing you’re going to want to do is to take a smaller subset of a large data set, either by filtering out certain columns or rows. Now let’s say we want to just work with the data from the West Midlands, and we’d like to omit some of the columns. We can choose a specific range of columns using select, like this:

You can use the filter command to do this. To give an example, filter can pick a single row in the following way:

-
uk_census_2021_religion_wmids <- uk_census_2021_religion %>% filter(geography=="West Midlands")
+
uk_census_2021_religion_wmids <- uk_census_2021_religion %>% filter(geography=="West Midlands")

Now we’ll use select in a different way to narrow our data to specific columns that are needed (no totals!).

Some readers will want to pause here and check out Hadley Wickham’s “R For Data Science” book, in the section, “Data visualisation” to get a fuller explanation of how to explore your data.
@@ -556,15 +550,15 @@ i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all

2.4 Making your first data visulation: the humble bar chart

We’ve got a nice lean set of data, so now it’s time to visualise this. We’ll start by making a pie chart:

-
uk_census_2021_religion_wmids <- uk_census_2021_religion_wmids %>% select(no_religion:no_response)
-uk_census_2021_religion_wmids <- gather(uk_census_2021_religion_wmids)
+
uk_census_2021_religion_wmids <- uk_census_2021_religion_wmids %>% select(no_religion:no_response)
+uk_census_2021_religion_wmids <- gather(uk_census_2021_religion_wmids)

There are two basic ways to do visualisations in R. You can work with basic functions in R, often called “base R” or you can work with an alternative library called ggplot:

2.4.1 Base R

-
df <- uk_census_2021_religion_wmids[order(uk_census_2021_religion_wmids$value,decreasing = TRUE),]
-barplot(height=df$value, names=df$key)
+
df <- uk_census_2021_religion_wmids[order(uk_census_2021_religion_wmids$value,decreasing = TRUE),]
+barplot(height=df$value, names=df$key)

@@ -573,48 +567,48 @@ i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all

2.4.2 GGPlot

-
ggplot(uk_census_2021_religion_wmids, aes(x = key, y = value)) +
-  geom_bar(stat = "identity")
+
ggplot(uk_census_2021_religion_wmids, aes(x = key, y = value)) +
+  geom_bar(stat = "identity")
-
2
+
2
-We’ll re-order the column by size. +We’ll re-order the column by size.

-
2ggplot(uk_census_2021_religion_wmids, aes(x= reorder(key,-value),value)) + geom_bar(stat ="identity")
+
2ggplot(uk_census_2021_religion_wmids, aes(x= reorder(key,-value),value)) + geom_bar(stat ="identity")

Let’s assume we’re working with a data set that doesn’t include a “totals” column and that we might want to get sums for each column. This is pretty easy to do in R:

-
1uk_census_2021_religion_totals <- uk_census_2021_religion %>% select(no_religion:no_response)
-uk_census_2021_religion_totals <- uk_census_2021_religion_totals %>%
-2   summarise(across(everything(), ~ sum(., na.rm = TRUE)))
-3uk_census_2021_religion_totals <- gather(uk_census_2021_religion_totals)
-4ggplot(uk_census_2021_religion_totals, aes(x= reorder(key,-value),value)) + geom_bar(stat ="identity")
+
1uk_census_2021_religion_totals <- uk_census_2021_religion %>% select(no_religion:no_response)
+uk_census_2021_religion_totals <- uk_census_2021_religion_totals %>%
+2   summarise(across(everything(), ~ sum(., na.rm = TRUE)))
+3uk_census_2021_religion_totals <- gather(uk_census_2021_religion_totals)
+4ggplot(uk_census_2021_religion_totals, aes(x= reorder(key,-value),value)) + geom_bar(stat ="identity")
-
1
+
1
-First, remove the column with region names and the totals for the regions as we want just integer data. +First, remove the column with region names and the totals for the regions as we want just integer data.
-
2
+
2
-Second calculate the totals. In this example we use the tidyverse library dplyr(), but you can also do this using base R with colsums() like this: uk_census_2021_religion_totals <- colSums(uk_census_2021_religion_totals, na.rm = TRUE). The downside with base R is that you’ll also need to convert the result into a dataframe for ggplot like this: uk_census_2021_religion_totals <- as.data.frame(uk_census_2021_religion_totals) +Second calculate the totals. In this example we use the tidyverse library dplyr(), but you can also do this using base R with colsums() like this: uk_census_2021_religion_totals <- colSums(uk_census_2021_religion_totals, na.rm = TRUE). The downside with base R is that you’ll also need to convert the result into a dataframe for ggplot like this: uk_census_2021_religion_totals <- as.data.frame(uk_census_2021_religion_totals)
-
3
+
3
-In order to visualise this data using ggplot, we need to shift this data from wide to long format. This is a quick job using gather() +In order to visualise this data using ggplot, we need to shift this data from wide to long format. This is a quick job using gather()
-
4
+
4
-Now plot it out and have a look! +Now plot it out and have a look!
@@ -624,17 +618,17 @@ i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all

You might have noticed that these two dataframes give us somewhat different results. But with data science, it’s much more interesting to compare these two side-by-side in a visualisation. We can join these two dataframes and plot the bars side by side using bind() - which can be done by columns with cbind() and rows using rbind():

-
uk_census_2021_religion_merged <- rbind(uk_census_2021_religion_totals, uk_census_2021_religion_wmids)
+
uk_census_2021_religion_merged <- rbind(uk_census_2021_religion_totals, uk_census_2021_religion_wmids)

Do you notice there’s going to be a problem here? How can we tell one set from the other? We need to add in something idenfiable first! This isn’t too hard to do as we can simply create a new column for each with identifiable information before we bind them:

-
uk_census_2021_religion_totals$dataset <- c("totals")
-uk_census_2021_religion_wmids$dataset <- c("wmids")
-uk_census_2021_religion_merged <- rbind(uk_census_2021_religion_totals, uk_census_2021_religion_wmids)
+
uk_census_2021_religion_totals$dataset <- c("totals")
+uk_census_2021_religion_wmids$dataset <- c("wmids")
+uk_census_2021_religion_merged <- rbind(uk_census_2021_religion_totals, uk_census_2021_religion_wmids)

Now we’re ready to plot out our data as a grouped barplot:

-
ggplot(uk_census_2021_religion_merged, aes(fill=dataset, x= reorder(key,-value), value)) + geom_bar(position="dodge", stat ="identity")
+
ggplot(uk_census_2021_religion_merged, aes(fill=dataset, x= reorder(key,-value), value)) + geom_bar(position="dodge", stat ="identity")

@@ -643,12 +637,12 @@ i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all

If you inspect our chart, you can see that we’re getting closer, but it’s not really that helpful to compare the totals. What we need to do is get percentages that can be compared side by side. This is easy to do using another dplyr feature mutate:

It’s worth noting that an alternative approach is to leave the numbers intact and simply label them differently so they render as percentages on your charts. You can do this with the `scales() library and the label_percent() function. The downside of this approach is that it won’t transfer to tables if you make them.
-
uk_census_2021_religion_totals <- uk_census_2021_religion_totals %>% 
-  dplyr::mutate(perc = scales::percent(value / sum(value), accuracy = 0.1, trim = FALSE))
-uk_census_2021_religion_wmids <- uk_census_2021_religion_wmids %>% 
-  dplyr::mutate(perc = scales::percent(value / sum(value), accuracy = 0.1, trim = FALSE))
-uk_census_2021_religion_merged <- rbind(uk_census_2021_religion_totals, uk_census_2021_religion_wmids)
-ggplot(uk_census_2021_religion_merged, aes(fill=dataset, x=key, y=perc)) + geom_bar(position="dodge", stat ="identity")
+
uk_census_2021_religion_totals <- uk_census_2021_religion_totals %>% 
+  dplyr::mutate(perc = scales::percent(value / sum(value), accuracy = 0.1, trim = FALSE))
+uk_census_2021_religion_wmids <- uk_census_2021_religion_wmids %>% 
+  dplyr::mutate(perc = scales::percent(value / sum(value), accuracy = 0.1, trim = FALSE))
+uk_census_2021_religion_merged <- rbind(uk_census_2021_religion_totals, uk_census_2021_religion_wmids)
+ggplot(uk_census_2021_religion_merged, aes(fill=dataset, x=key, y=perc)) + geom_bar(position="dodge", stat ="identity")

@@ -657,7 +651,7 @@ i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all

We’ll draw this data into comparison with later sets in the next chapter. But the one glaring issue which remains for our chart is that it’s lacking in really any aesthetic refinements. This is where ggplot really shines as a tool as you can add all sorts of things.

These are basically just added to our ggplot code. So, for example, let’s say we want to improve the colours used for our bars. You can specify the formatting for the fill on the scale using scale_fill_brewer. This uses a particular tool (and a personal favourite of mine) called colorbrewer. Part of my appreciation of this tool is that you can pick colours which are not just visually pleasing, and produce useful contrast / complementary schemes, but you can also work proactively to accommodate colourblindness. Working with colour schemes which can be divergent in a visually obvious way will be even more important when we work on geospatial data and maps in a later chapter.

-
ggplot(uk_census_2021_religion_merged, aes(fill=dataset, x=key, y=perc)) + geom_bar(position="dodge", stat ="identity") + scale_fill_brewer(palette = "Set1")
+
ggplot(uk_census_2021_religion_merged, aes(fill=dataset, x=key, y=perc)) + geom_bar(position="dodge", stat ="identity") + scale_fill_brewer(palette = "Set1")

@@ -665,15 +659,15 @@ i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all

We might also want to add a border to our bars to make them more visually striking (notice the addition of color to the geom_bar below. I’ve also added reorder() to the x value to sort descending from the largest to smallest.

You can find more information about reordering ggplots on the R Graph gallery.
-
uk_census_2021_religion_merged$dataset <- factor(uk_census_2021_religion_merged$dataset, levels = c('wmids', 'totals'))
-ggplot(uk_census_2021_religion_merged, aes(fill=fct_reorder(dataset, value), x=reorder(key,-value),value, y=perc)) + geom_bar(position="dodge", stat ="identity", colour = "black") + scale_fill_brewer(palette = "Set1")
+
uk_census_2021_religion_merged$dataset <- factor(uk_census_2021_religion_merged$dataset, levels = c('wmids', 'totals'))
+ggplot(uk_census_2021_religion_merged, aes(fill=fct_reorder(dataset, value), x=reorder(key,-value),value, y=perc)) + geom_bar(position="dodge", stat ="identity", colour = "black") + scale_fill_brewer(palette = "Set1")

We can fine tune a few other visual features here as well, like adding a title with ggtitle and a them with some prettier fonts with theme_ipsum() (which requires the hrbrthemes() library). We can also remove the x and y axis labels (not the data labels, which are rather important).

-
ggplot(uk_census_2021_religion_merged, aes(fill=fct_reorder(dataset, value), x=reorder(key,-value),value, y=perc)) + geom_bar(position="dodge", stat ="identity", colour = "black") + scale_fill_brewer(palette = "Set1") + ggtitle("Religious Affiliation in the UK: 2021") + xlab("") + ylab("")
+
ggplot(uk_census_2021_religion_merged, aes(fill=fct_reorder(dataset, value), x=reorder(key,-value),value, y=perc)) + geom_bar(position="dodge", stat ="identity", colour = "black") + scale_fill_brewer(palette = "Set1") + ggtitle("Religious Affiliation in the UK: 2021") + xlab("") + ylab("")

@@ -687,7 +681,7 @@ i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all

Finally, researchers have also found that when people are asked to mark their religious affiliation, sometimes they can prefer to mark more than one answer. A person might consider themselves to be “Muslim” but also “Spiritual but not religious” preferring the combination of those identities. It is also the case that respondents can identify with more unexpected hybrid religious identities, such as “Christian” and “Hindu”. The census only allows respondents to tick a single box for the religion category. It is worth noting that, in contrast, the responses for ethnicity allow for combinations. Given that this is the case, it’s impossible to know which way a person went at the fork in the road as they were forced to choose just one half of this kind of hybrid identity. Finally, it is interesting to wonder exactly what it means for a person when they tick a box like this. Is it because they attend synagogue on a weekly basis? Some persons would consider weekly attendance at workship a prerequisite for membership in a group, but others would not. Indeed we can infer from surveys and research which aims to track rates of participation in weekly worship that many people who tick boxes for particular religious identities on the census have never attended a worship service at all.

What does this mean for our results? Are they completely unreliable and invalid? I don’t think this is the case or that taking a clear-eyed look at the force and stability of our underlying data should be cause for despair. Instead, the most appropriate response is humility. Someone has made a statement which is recorded in the census, of this we can be sure. They felt it to be an accurate response on some level based on the information they had at the time. And with regard to the census, it is a massive, almost completely population level, sample so there is additional validity there. The easiest way to represent all this reality in the form of speaking truthfully about our data is to acknowledge that however valid it may seem, it is nonetheless a snapshot. For this reason, I would always advise that the best title for a chart is one which specifies the data set. We should also probably do something different with those non-responses:

-
ggplot(uk_census_2021_religion_merged, aes(fill=fct_reorder(dataset, value), x=reorder(key,-value),value, y=perc)) + geom_bar(position="dodge", stat ="identity", colour = "black") + scale_fill_brewer(palette = "Set1") + ggtitle("Religious Affiliation in the 2021 Census of England and Wales") + xlab("") + ylab("")
+
ggplot(uk_census_2021_religion_merged, aes(fill=fct_reorder(dataset, value), x=reorder(key,-value),value, y=perc)) + geom_bar(position="dodge", stat ="identity", colour = "black") + scale_fill_brewer(palette = "Set1") + ggtitle("Religious Affiliation in the 2021 Census of England and Wales") + xlab("") + ylab("")

@@ -713,204 +707,169 @@ What is Nomis?

For the UK, census data is made available for programmatic research like this via an organisation called NOMIS. Luckily for us, there is an R library you can use to access nomis directly which greatly simplifies the process of pulling data down from the platform. It’s worth noting that if you’re not in the UK, there are similar options for other countries. Nearly every R textbook I’ve ever seen works with USA census data, so you’ll find plenty of documentation available on the tools you can use for US Census data. Similarly for the EU, Canada, Austrailia etc.

-

Here’s the process to identify a dataset within the nomis platform:

-
-
# Process to explore nomis() data for specific datasets
-library(nomisr)
-# temporarily commenting out until renv can be implemented and runtime errors in other environments avoided:
-#religion_search <- nomis_search(name = "*Religion*")
-#religion_measures <- nomis_get_metadata("ST104", "measures")
-#tibble::glimpse(religion_measures)
-#religion_geography <- nomis_get_metadata("NM_529_1", "geography", "TYPE")
-
+

If you want to draw some data from the nomis platform yourself in R, have a look at the companion cookbook repository.

-
library(nomisr)
-# Get table of Census 2011 religion data from nomis
-# temporarily commenting out until renv can be implemented and runtime errors in other environments avoided:
-#z <- nomis_get_data(id = "NM_529_1", time = "latest", geography = "TYPE499", measures=c(20301))
-#saveRDS(z, file = "z.rds")
-z <- readRDS(file = (here("example_data", "z.rds")))
-
-# Filter down to simplified dataset with England / Wales and percentages without totals
-uk_census_2011_religion <- filter(z, GEOGRAPHY_NAME=="England and Wales" & RURAL_URBAN_NAME=="Total" & C_RELPUK11_NAME != "All categories: Religion")
-# Drop unnecessary columns
-uk_census_2011_religion <- select(uk_census_2011_religion, C_RELPUK11_NAME, OBS_VALUE)
-# Plot results
-plot1 <- ggplot(uk_census_2011_religion, aes(x = C_RELPUK11_NAME, y = OBS_VALUE)) + geom_bar(stat = "identity") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
-# ggsave(filename = "plot.png", plot = plot1)
-
-# grab daata from nomis for 2001 census religion / ethnicity
-
-# temporarily commenting out until renv can be implemented and runtime errors in other environments avoided:
-#z0 <- nomis_get_data(id = "NM_1872_1", time = "latest", geography = "TYPE499", measures=c(20100))
-#saveRDS(z0, file = "z0.rds")
-z0 <- readRDS(file = (here("example_data", "z0.rds")))
-
-# select relevant columns
-uk_census_2001_religion_ethnicity <- select(z0, GEOGRAPHY_NAME, C_RELPUK11_NAME, C_ETHHUK11_NAME, OBS_VALUE)
-# Filter down to simplified dataset with England / Wales and percentages without totals
-uk_census_2001_religion_ethnicity <- filter(uk_census_2001_religion_ethnicity, GEOGRAPHY_NAME=="England and Wales" & C_RELPUK11_NAME != "All categories: Religion")
-# Simplify data to only include general totals and omit subcategories
-uk_census_2001_religion_ethnicity <- uk_census_2001_religion_ethnicity %>% filter(grepl('Total', C_ETHHUK11_NAME))
-
-# grab data from nomis for 2011 census religion / ethnicity table
-# commenting out nomis_get temporarily until I can get renv working properly here
-#z1 <- nomis_get_data(id = "NM_659_1", time = "latest", geography = "TYPE499", measures=c(20100))
-#saveRDS(z1, file = "z1.rds")
-z1 <- readRDS(file = (here("example_data", "z1.rds")))
-
-# select relevant columns
-uk_census_2011_religion_ethnicity <- select(z1, GEOGRAPHY_NAME, C_RELPUK11_NAME, C_ETHPUK11_NAME, OBS_VALUE)
-# Filter down to simplified dataset with England / Wales and percentages without totals
-uk_census_2011_religion_ethnicity <- filter(uk_census_2011_religion_ethnicity, GEOGRAPHY_NAME=="England and Wales" & C_RELPUK11_NAME != "All categories: Religion" & C_ETHPUK11_NAME != "All categories: Ethnic group")
-# Simplify data to only include general totals and omit subcategories
-uk_census_2011_religion_ethnicity <- uk_census_2011_religion_ethnicity %>% filter(grepl('Total', C_ETHPUK11_NAME))
-
-# grab data from nomis for 2021 census religion / ethnicity table
-#z2 <- nomis_get_data(id = "NM_2131_1", time = "latest", geography = "TYPE499", measures=c(20100))
-#saveRDS(z2, file = "z2.rds")
-z2 <- readRDS(file = (here("example_data", "z2.rds")))
-
-# select relevant columns
-uk_census_2021_religion_ethnicity <- select(z2, GEOGRAPHY_NAME, C2021_RELIGION_10_NAME, C2021_ETH_8_NAME, OBS_VALUE)
-# Filter down to simplified dataset with England / Wales and percentages without totals
-uk_census_2021_religion_ethnicity <- filter(uk_census_2021_religion_ethnicity, GEOGRAPHY_NAME=="England and Wales" & C2021_RELIGION_10_NAME != "Total" & C2021_ETH_8_NAME != "Total")
-# 2021 census includes white sub-groups so we need to omit those so we just have totals:
-uk_census_2021_religion_ethnicity <- filter(uk_census_2021_religion_ethnicity, C2021_ETH_8_NAME != "White: English, Welsh, Scottish, Northern Irish or British" & C2021_ETH_8_NAME != "White: Irish" & C2021_ETH_8_NAME != "White: Gypsy or Irish Traveller, Roma or Other White")
-
-ggplot(uk_census_2011_religion_ethnicity, aes(fill=C_ETHPUK11_NAME, x=C_RELPUK11_NAME, y=OBS_VALUE)) + geom_bar(position="dodge", stat ="identity", colour = "black") + scale_fill_brewer(palette = "Set1") + ggtitle("Religious Affiliation in the 2021 Census of England and Wales") + xlab("") + ylab("") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
+
# Get table of Census 2011 religion data from nomis
+# Note: for reproducible code used to generate the dataset used in the book, see the cookbook here: 
+z <- readRDS(file = (here("example_data", "z.rds")))
+
+# Filter down to simplified dataset with England / Wales and percentages without totals
+uk_census_2011_religion <- filter(z, GEOGRAPHY_NAME=="England and Wales" & RURAL_URBAN_NAME=="Total" & C_RELPUK11_NAME != "All categories: Religion")
+# Drop unnecessary columns
+uk_census_2011_religion <- select(uk_census_2011_religion, C_RELPUK11_NAME, OBS_VALUE)
+# Plot results
+plot1 <- ggplot(uk_census_2011_religion, aes(x = C_RELPUK11_NAME, y = OBS_VALUE)) + geom_bar(stat = "identity") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
+# ggsave(filename = "plot.png", plot = plot1)
+
+# grab daata from nomis for 2001 census religion / ethnicity
+z0 <- readRDS(file = (here("example_data", "z0.rds")))
+
+# select relevant columns
+uk_census_2001_religion_ethnicity <- select(z0, GEOGRAPHY_NAME, C_RELPUK11_NAME, C_ETHHUK11_NAME, OBS_VALUE)
+# Filter down to simplified dataset with England / Wales and percentages without totals
+uk_census_2001_religion_ethnicity <- filter(uk_census_2001_religion_ethnicity, GEOGRAPHY_NAME=="England and Wales" & C_RELPUK11_NAME != "All categories: Religion")
+# Simplify data to only include general totals and omit subcategories
+uk_census_2001_religion_ethnicity <- uk_census_2001_religion_ethnicity %>% filter(grepl('Total', C_ETHHUK11_NAME))
+
+# grab data from nomis for 2011 census religion / ethnicity table
+z1 <- readRDS(file = (here("example_data", "z1.rds")))
+
+# select relevant columns
+uk_census_2011_religion_ethnicity <- select(z1, GEOGRAPHY_NAME, C_RELPUK11_NAME, C_ETHPUK11_NAME, OBS_VALUE)
+# Filter down to simplified dataset with England / Wales and percentages without totals
+uk_census_2011_religion_ethnicity <- filter(uk_census_2011_religion_ethnicity, GEOGRAPHY_NAME=="England and Wales" & C_RELPUK11_NAME != "All categories: Religion" & C_ETHPUK11_NAME != "All categories: Ethnic group")
+# Simplify data to only include general totals and omit subcategories
+uk_census_2011_religion_ethnicity <- uk_census_2011_religion_ethnicity %>% filter(grepl('Total', C_ETHPUK11_NAME))
+
+# grab data from nomis for 2021 census religion / ethnicity table
+z2 <- readRDS(file = (here("example_data", "z2.rds")))
+
+# select relevant columns
+uk_census_2021_religion_ethnicity <- select(z2, GEOGRAPHY_NAME, C2021_RELIGION_10_NAME, C2021_ETH_8_NAME, OBS_VALUE)
+# Filter down to simplified dataset with England / Wales and percentages without totals
+uk_census_2021_religion_ethnicity <- filter(uk_census_2021_religion_ethnicity, GEOGRAPHY_NAME=="England and Wales" & C2021_RELIGION_10_NAME != "Total" & C2021_ETH_8_NAME != "Total")
+# 2021 census includes white sub-groups so we need to omit those so we just have totals:
+uk_census_2021_religion_ethnicity <- filter(uk_census_2021_religion_ethnicity, C2021_ETH_8_NAME != "White: English, Welsh, Scottish, Northern Irish or British" & C2021_ETH_8_NAME != "White: Irish" & C2021_ETH_8_NAME != "White: Gypsy or Irish Traveller, Roma or Other White")
+
+ggplot(uk_census_2011_religion_ethnicity, aes(fill=C_ETHPUK11_NAME, x=C_RELPUK11_NAME, y=OBS_VALUE)) + geom_bar(position="dodge", stat ="identity", colour = "black") + scale_fill_brewer(palette = "Set1") + ggtitle("Religious Affiliation in the 2021 Census of England and Wales") + xlab("") + ylab("") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
-

+

The trouble with using grouped bars here, as you can see, is that there are quite sharp disparities which make it hard to compare in meaningful ways. We could use logarithmic rather than linear scaling as an option, but this is hard for many general public audiences to apprecaite without guidance. One alternative quick fix is to extract data from “white” respondents which can then be placed in a separate chart with a different scale.

-
# Filter down to simplified dataset with England / Wales and percentages without totals
-uk_census_2011_religion_ethnicity_white <- filter(uk_census_2011_religion_ethnicity, C_ETHPUK11_NAME == "White: Total")
-uk_census_2011_religion_ethnicity_nonwhite <- filter(uk_census_2011_religion_ethnicity, C_ETHPUK11_NAME != "White: Total")
-
-ggplot(uk_census_2011_religion_ethnicity_nonwhite, aes(fill=C_ETHPUK11_NAME, x=C_RELPUK11_NAME, y=OBS_VALUE)) + geom_bar(position="dodge", stat ="identity", colour = "black") + scale_fill_brewer(palette = "Set1") + ggtitle("Religious Affiliation in the 2021 Census of England and Wales") + xlab("") + ylab("") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
+
# Filter down to simplified dataset with England / Wales and percentages without totals
+uk_census_2011_religion_ethnicity_white <- filter(uk_census_2011_religion_ethnicity, C_ETHPUK11_NAME == "White: Total")
+uk_census_2011_religion_ethnicity_nonwhite <- filter(uk_census_2011_religion_ethnicity, C_ETHPUK11_NAME != "White: Total")
+
+ggplot(uk_census_2011_religion_ethnicity_nonwhite, aes(fill=C_ETHPUK11_NAME, x=C_RELPUK11_NAME, y=OBS_VALUE)) + geom_bar(position="dodge", stat ="identity", colour = "black") + scale_fill_brewer(palette = "Set1") + ggtitle("Religious Affiliation in the 2021 Census of England and Wales") + xlab("") + ylab("") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
-

+

This still doesn’t quite render with as much visual clarity and communication as I’d like. For a better look, we can use a technique in R called “faceting” to create a series of small charts which can be viewed alongside one another.

-
ggplot(uk_census_2011_religion_ethnicity_nonwhite, aes(x=C_RELPUK11_NAME, y=OBS_VALUE)) + geom_bar(position="dodge", stat ="identity", colour = "black") + facet_wrap(~C_ETHPUK11_NAME, ncol = 2) + scale_fill_brewer(palette = "Set1") + ggtitle("Religious Affiliation in the 2011 Census of England and Wales") + xlab("") + ylab("") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
+
ggplot(uk_census_2011_religion_ethnicity_nonwhite, aes(x=C_RELPUK11_NAME, y=OBS_VALUE)) + geom_bar(position="dodge", stat ="identity", colour = "black") + facet_wrap(~C_ETHPUK11_NAME, ncol = 2) + scale_fill_brewer(palette = "Set1") + ggtitle("Religious Affiliation in the 2011 Census of England and Wales") + xlab("") + ylab("") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
-

+

For our finale chart, I’d like to take the faceted chart we’ve just done, and add in totals for the previous two census years (2001 and 2011) so we can see how trends are changing in terms of religious affiliation within ethnic self-identification categories. We’ll draw on some techniques we’re already developed above using rbind() to connect up each of these charts (after we’ve added a column identifying each chart by the census year). We will also need to use one new technique to change the wording of ethnic categories as this isn’t consistent from one census to the next and ggplot will struggle to chart things if the terms being used are exactly the same. We’ll use mutate() again to accomplish this with some slightly different code.

-
# First add column to each dataframe so we don't lose track of the census it comes from:
-uk_census_2001_religion_ethnicity$dataset <- c("2001")
-uk_census_2011_religion_ethnicity$dataset <- c("2011")
-uk_census_2021_religion_ethnicity$dataset <- c("2021")
-
-# Let's tidy the names of each column:
-
-names(uk_census_2001_religion_ethnicity) <- c("Geography", "Religion", "Ethnicity", "Value", "Year")
-names(uk_census_2011_religion_ethnicity) <- c("Geography", "Religion", "Ethnicity", "Value", "Year")
-names(uk_census_2021_religion_ethnicity) <- c("Geography", "Religion", "Ethnicity", "Value", "Year")
-
-# Next we need to change the terms using mutate()
-uk_census_2001_religion_ethnicity <- uk_census_2001_religion_ethnicity %>% 
-  mutate(Ethnicity = str_replace_all(Ethnicity, 
-            pattern = "^White: Total$", replacement = "White")) %>%
-  mutate(Ethnicity = str_replace_all(Ethnicity, 
-            pattern = "^Mixed: Total$", replacement = "Mixed")) %>%
-  mutate(Ethnicity = str_replace_all(Ethnicity, 
-            pattern = "^Asian: Total$", replacement = "Asian")) %>%
-  mutate(Ethnicity = str_replace_all(Ethnicity, 
-            pattern = "^Black or Black British: Total$", replacement = "Black")) %>%
-  mutate(Ethnicity = str_replace_all(Ethnicity, 
-            pattern = "^Chinese or Other ethnic group: Total$", replacement = "Other"))
-  
-uk_census_2011_religion_ethnicity <- uk_census_2011_religion_ethnicity %>% 
-  mutate(Ethnicity = str_replace_all(Ethnicity, 
-            pattern = "^White: Total$", replacement = "White")) %>%
-  mutate(Ethnicity = str_replace_all(Ethnicity, 
-            pattern = "^Mixed/multiple ethnic group: Total$", replacement = "Mixed")) %>%
-  mutate(Ethnicity = str_replace_all(Ethnicity, 
-            pattern = "^Asian/Asian British: Total$", replacement = "Asian")) %>%
-  mutate(Ethnicity = str_replace_all(Ethnicity, 
-            pattern = "^Black/African/Caribbean/Black British: Total$", replacement = "Black")) %>%
-  mutate(Ethnicity = str_replace_all(Ethnicity, 
-            pattern = "^Other ethnic group: Total$", replacement = "Other"))
-
-uk_census_2021_religion_ethnicity <- uk_census_2021_religion_ethnicity %>% 
-  mutate(Ethnicity = str_replace_all(Ethnicity, 
-            pattern = "^White: Total$", replacement = "White")) %>%
-  mutate(Ethnicity = str_replace_all(Ethnicity, 
-            pattern = "^Mixed or Multiple ethnic groups$", replacement = "Mixed")) %>%
-  mutate(Ethnicity = str_replace_all(Ethnicity, 
-            pattern = "^Asian, Asian British or Asian Welsh$", replacement = "Asian")) %>%
-  mutate(Ethnicity = str_replace_all(Ethnicity, 
-            pattern = "^Black, Black British, Black Welsh, Caribbean or African$", replacement = "Black")) %>%
-  mutate(Ethnicity = str_replace_all(Ethnicity, 
-            pattern = "^Other ethnic group$", replacement = "Other"))
-
-# Now let's merge the tables:
-
-uk_census_merged_religion_ethnicity <- rbind(uk_census_2021_religion_ethnicity, uk_census_2011_religion_ethnicity)
-
-uk_census_merged_religion_ethnicity <- rbind(uk_census_merged_religion_ethnicity, uk_census_2001_religion_ethnicity)
-
-# As above, we'll split out non-white and white:
-
-uk_census_merged_religion_ethnicity_nonwhite <- filter(uk_census_merged_religion_ethnicity, Ethnicity != "White")
-
-# Time to plot!
-
-ggplot(uk_census_merged_religion_ethnicity_nonwhite, aes(fill=Year, x=Religion, y=Value)) + geom_bar(position="dodge", stat ="identity", colour = "black") + facet_wrap(~Ethnicity, ncol = 2) + scale_fill_brewer(palette = "Set1") + ggtitle("Religious Affiliation in the 2001-2021 Census of England and Wales") + xlab("") + ylab("") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
+
# First add column to each dataframe so we don't lose track of the census it comes from:
+uk_census_2001_religion_ethnicity$dataset <- c("2001")
+uk_census_2011_religion_ethnicity$dataset <- c("2011")
+uk_census_2021_religion_ethnicity$dataset <- c("2021")
+
+# Let's tidy the names of each column:
+
+names(uk_census_2001_religion_ethnicity) <- c("Geography", "Religion", "Ethnicity", "Value", "Year")
+names(uk_census_2011_religion_ethnicity) <- c("Geography", "Religion", "Ethnicity", "Value", "Year")
+names(uk_census_2021_religion_ethnicity) <- c("Geography", "Religion", "Ethnicity", "Value", "Year")
+
+# Next we need to change the terms using mutate()
+uk_census_2001_religion_ethnicity <- uk_census_2001_religion_ethnicity %>% 
+  mutate(Ethnicity = str_replace_all(Ethnicity, 
+            pattern = "^White: Total$", replacement = "White")) %>%
+  mutate(Ethnicity = str_replace_all(Ethnicity, 
+            pattern = "^Mixed: Total$", replacement = "Mixed")) %>%
+  mutate(Ethnicity = str_replace_all(Ethnicity, 
+            pattern = "^Asian: Total$", replacement = "Asian")) %>%
+  mutate(Ethnicity = str_replace_all(Ethnicity, 
+            pattern = "^Black or Black British: Total$", replacement = "Black")) %>%
+  mutate(Ethnicity = str_replace_all(Ethnicity, 
+            pattern = "^Chinese or Other ethnic group: Total$", replacement = "Other"))
+  
+uk_census_2011_religion_ethnicity <- uk_census_2011_religion_ethnicity %>% 
+  mutate(Ethnicity = str_replace_all(Ethnicity, 
+            pattern = "^White: Total$", replacement = "White")) %>%
+  mutate(Ethnicity = str_replace_all(Ethnicity, 
+            pattern = "^Mixed/multiple ethnic group: Total$", replacement = "Mixed")) %>%
+  mutate(Ethnicity = str_replace_all(Ethnicity, 
+            pattern = "^Asian/Asian British: Total$", replacement = "Asian")) %>%
+  mutate(Ethnicity = str_replace_all(Ethnicity, 
+            pattern = "^Black/African/Caribbean/Black British: Total$", replacement = "Black")) %>%
+  mutate(Ethnicity = str_replace_all(Ethnicity, 
+            pattern = "^Other ethnic group: Total$", replacement = "Other"))
+
+uk_census_2021_religion_ethnicity <- uk_census_2021_religion_ethnicity %>% 
+  mutate(Ethnicity = str_replace_all(Ethnicity, 
+            pattern = "^White: Total$", replacement = "White")) %>%
+  mutate(Ethnicity = str_replace_all(Ethnicity, 
+            pattern = "^Mixed or Multiple ethnic groups$", replacement = "Mixed")) %>%
+  mutate(Ethnicity = str_replace_all(Ethnicity, 
+            pattern = "^Asian, Asian British or Asian Welsh$", replacement = "Asian")) %>%
+  mutate(Ethnicity = str_replace_all(Ethnicity, 
+            pattern = "^Black, Black British, Black Welsh, Caribbean or African$", replacement = "Black")) %>%
+  mutate(Ethnicity = str_replace_all(Ethnicity, 
+            pattern = "^Other ethnic group$", replacement = "Other"))
+
+# Now let's merge the tables:
+
+uk_census_merged_religion_ethnicity <- rbind(uk_census_2021_religion_ethnicity, uk_census_2011_religion_ethnicity)
+
+uk_census_merged_religion_ethnicity <- rbind(uk_census_merged_religion_ethnicity, uk_census_2001_religion_ethnicity)
+
+# As above, we'll split out non-white and white:
+
+uk_census_merged_religion_ethnicity_nonwhite <- filter(uk_census_merged_religion_ethnicity, Ethnicity != "White")
+
+# Time to plot!
+
+ggplot(uk_census_merged_religion_ethnicity_nonwhite, aes(fill=Year, x=Religion, y=Value)) + geom_bar(position="dodge", stat ="identity", colour = "black") + facet_wrap(~Ethnicity, ncol = 2) + scale_fill_brewer(palette = "Set1") + ggtitle("Religious Affiliation in the 2001-2021 Census of England and Wales") + xlab("") + ylab("") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
-

+

There are a few formatting issues which remain. Our y-axis number labels are in scientific format which isn’t really very easy to read. You can use the very powerful and flexible scales() library to bring in some more readable formatting of numbers in a variety of places in R including in ggplot visualizations.

-
library(scales)
-
-

-Attaching package: 'scales'
-
-
-
The following object is masked from 'package:purrr':
-
-    discard
-
-
-
The following object is masked from 'package:readr':
-
-    col_factor
-
-
ggplot(uk_census_merged_religion_ethnicity_nonwhite, aes(fill=Year, x=Religion, y=Value)) + geom_bar(position="dodge", stat ="identity", colour = "black") + facet_wrap(~Ethnicity, ncol = 2) + scale_fill_brewer(palette = "Set1") + scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6), breaks = breaks_extended(8)) + ggtitle("Religious Affiliation in the 2001-2021 Census of England and Wales") + xlab("") + ylab("") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
+
library(scales) |> suppressPackageStartupMessages()
+ggplot(uk_census_merged_religion_ethnicity_nonwhite, aes(fill=Year, x=Religion, y=Value)) + geom_bar(position="dodge", stat ="identity", colour = "black") + facet_wrap(~Ethnicity, ncol = 2) + scale_fill_brewer(palette = "Set1") + scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6), breaks = breaks_extended(8)) + ggtitle("Religious Affiliation in the 2001-2021 Census of England and Wales") + xlab("") + ylab("") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
-

+

-
# https://ggplot2-book.org/scales-position#sec-position-continuous-breaks
+
# https://ggplot2-book.org/scales-position#sec-position-continuous-breaks

This chart shows an increase in almost every category, though it’s a bit hard to read in some cases. However, this information is based on the increase in raw numbers. It’s possbile that numbers may be going up, but in some cases the percentage share for a particular category has actually gone down. Let’s transform and visualise our data as percentages to see what kind of trends we can actually isolate:

-
uk_census_merged_religion_ethnicity <- uk_census_merged_religion_ethnicity %>%
-  group_by(Ethnicity, Year) %>%
-  dplyr::mutate(Percent = Value/sum(Value))
-
-ggplot(uk_census_merged_religion_ethnicity, aes(fill=Year, x=Religion, y=Percent)) + geom_bar(position="dodge", stat ="identity", colour = "black") + facet_wrap(~Ethnicity, scales="free_x") + scale_fill_brewer(palette = "Set1") + scale_y_continuous(labels = scales::percent) + ggtitle("Religious Affiliation in the 2001-2021 Census of England and Wales") + xlab("") + ylab("") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
+
uk_census_merged_religion_ethnicity <- uk_census_merged_religion_ethnicity %>%
+  group_by(Ethnicity, Year) %>%
+  dplyr::mutate(Percent = Value/sum(Value))
+
+ggplot(uk_census_merged_religion_ethnicity, aes(fill=Year, x=Religion, y=Percent)) + geom_bar(position="dodge", stat ="identity", colour = "black") + facet_wrap(~Ethnicity, scales="free_x") + scale_fill_brewer(palette = "Set1") + scale_y_continuous(labels = scales::percent) + ggtitle("Religious Affiliation in the 2001-2021 Census of England and Wales") + xlab("") + ylab("") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
-

+

Now you can see why this shift is important - the visualisation tells a completely different story in some cases across the two different charts. In the first, working off raw numbers we see a net increase in Christianity across all categories. But if we take into account the fact that the overall share of population is growing for each of these groups, their actual composition is changing in a different direction. The proportion of each group is declining across the three census periods (albeit with an exception for the “Other” category from 2011 to 2021).

To highlight a few features of this final plot, I’ve used a specific feature within facet_wrap scales = "free_x" to let each of the individual facets adjust the total range on the x-axis. Since we’re looking at trends here and not absolute values, having correspondence across scales isn’t important and this makes for something a bit more visually tidy. I’ve also shifted the code for scale_y_continuous to render values as percentages (rather than millions).

In case you want to print this plot out and hang it on your wall, you can use the ggsave tool to render the chart as an image file:

-
plot1 <- ggplot(uk_census_merged_religion_ethnicity, aes(fill=Year, x=Religion, y=Percent)) + geom_bar(position="dodge", stat ="identity", colour = "black") + facet_wrap(~Ethnicity, scales="free_x") + scale_fill_brewer(palette = "Set1") + scale_y_continuous(labels = scales::percent) + ggtitle("Religious Affiliation in the 2001-2021 Census of England and Wales") + xlab("") + ylab("") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
-
-ggsave("chart.png", plot=plot1, width = 8, height = 10, units=c("in"))
+
plot1 <- ggplot(uk_census_merged_religion_ethnicity, aes(fill=Year, x=Religion, y=Percent)) + geom_bar(position="dodge", stat ="identity", colour = "black") + facet_wrap(~Ethnicity, scales="free_x") + scale_fill_brewer(palette = "Set1") + scale_y_continuous(labels = scales::percent) + ggtitle("Religious Affiliation in the 2001-2021 Census of England and Wales") + xlab("") + ylab("") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
+
+ggsave("chart.png", plot=plot1, width = 8, height = 10, units=c("in"))
diff --git a/docs/chapter_1_files/figure-html/unnamed-chunk-18-1.png b/docs/chapter_1_files/figure-html/unnamed-chunk-18-1.png index ad999cc..b08d52c 100644 Binary files a/docs/chapter_1_files/figure-html/unnamed-chunk-18-1.png and b/docs/chapter_1_files/figure-html/unnamed-chunk-18-1.png differ diff --git a/docs/chapter_1_files/figure-html/unnamed-chunk-19-1.png b/docs/chapter_1_files/figure-html/unnamed-chunk-19-1.png index b08d52c..94a034b 100644 Binary files a/docs/chapter_1_files/figure-html/unnamed-chunk-19-1.png and b/docs/chapter_1_files/figure-html/unnamed-chunk-19-1.png differ diff --git a/docs/chapter_1_files/figure-html/unnamed-chunk-20-1.png b/docs/chapter_1_files/figure-html/unnamed-chunk-20-1.png index 94a034b..bce4f73 100644 Binary files a/docs/chapter_1_files/figure-html/unnamed-chunk-20-1.png and b/docs/chapter_1_files/figure-html/unnamed-chunk-20-1.png differ diff --git a/docs/chapter_1_files/figure-html/unnamed-chunk-21-1.png b/docs/chapter_1_files/figure-html/unnamed-chunk-21-1.png index bce4f73..62392dc 100644 Binary files a/docs/chapter_1_files/figure-html/unnamed-chunk-21-1.png and b/docs/chapter_1_files/figure-html/unnamed-chunk-21-1.png differ diff --git a/docs/chapter_1_files/figure-html/unnamed-chunk-22-1.png b/docs/chapter_1_files/figure-html/unnamed-chunk-22-1.png index 62392dc..3e227c8 100644 Binary files a/docs/chapter_1_files/figure-html/unnamed-chunk-22-1.png and b/docs/chapter_1_files/figure-html/unnamed-chunk-22-1.png differ diff --git a/docs/chapter_2.html b/docs/chapter_2.html index daa0f14..df94999 100644 --- a/docs/chapter_2.html +++ b/docs/chapter_2.html @@ -355,15 +355,15 @@ So who’s religious?
1
-First we generate new a dataframe with sums per category and +First we generate new a dataframe with sums per category and
2
-…sort in descending order +…sort in descending order
3
-Then we add new column with percentages based on the sums you’ve just generated +Then we add new column with percentages based on the sums you’ve just generated
diff --git a/docs/chapter_3.html b/docs/chapter_3.html index 8b773e1..b7e4e32 100644 --- a/docs/chapter_3.html +++ b/docs/chapter_3.html @@ -20,6 +20,40 @@ ul.task-list li input[type="checkbox"] { margin: 0 0.8em 0.2em -1em; /* quarto-specific, see https://github.com/quarto-dev/quarto-cli/issues/4556 */ vertical-align: middle; } +/* CSS for syntax highlighting */ +pre > code.sourceCode { white-space: pre; position: relative; } +pre > code.sourceCode > span { display: inline-block; line-height: 1.25; } +pre > code.sourceCode > span:empty { height: 1.2em; } +.sourceCode { overflow: visible; } +code.sourceCode > span { color: inherit; text-decoration: inherit; } +div.sourceCode { margin: 1em 0; } +pre.sourceCode { margin: 0; } +@media screen { +div.sourceCode { overflow: auto; } +} +@media print { +pre > code.sourceCode { white-space: pre-wrap; } +pre > code.sourceCode > span { text-indent: -5em; padding-left: 5em; } +} +pre.numberSource code + { counter-reset: source-line 0; } +pre.numberSource code > span + { position: relative; left: -4em; counter-increment: source-line; } +pre.numberSource code > span > a:first-child::before + { content: counter(source-line); + position: relative; left: -1em; text-align: right; vertical-align: baseline; + border: none; display: inline-block; + -webkit-touch-callout: none; -webkit-user-select: none; + -khtml-user-select: none; -moz-user-select: none; + -ms-user-select: none; user-select: none; + padding: 0 4px; width: 4em; + } +pre.numberSource { margin-left: 3em; padding-left: 4px; } +div.sourceCode + { } +@media screen { +pre > code.sourceCode > span > a:first-child::before { text-decoration: underline; } +} /* CSS for citations */ div.csl-bib-body { } div.csl-entry { @@ -176,7 +210,9 @@ div.csl-indent {

Table of contents

@@ -200,8 +236,196 @@ div.csl-indent { +

Until recently, most data science books didn’t have a section on geospatial data. It was considered a specialist form of research best left to GIS technicians who tended to use proprietary tools like ArcGIS. This has changed significantly in the past five years, but you’ll still be hard pressed to find an introduction to the subject which strays very far from a few simple data sets (mostly of the USA) and relatively uncomplicated geospatial operations. I actually first began learning R, back in 2013, right when open source geospatial research tools were beginning to be developed with quite a lot more energy and geospatial data is my personal favourite data science playground, so in this book we’re going to go much deeper than is usual. There are also good reasons to take things a few steps further in the particular forms of data and inquiry that religion takes us into.

+

Recommend https://r-spatial.org/book/

+

Geospatial data is, in the most basic form, working with maps. This means that most of your data can be a quite simple dataframe, e.g. just a list of names or categories associated with a set of X and Y coordinates. Once you have a set of items, however, things get interesting very quickly, as you can layer data sets on top of one another. We’re going to begin this chapter by developing a geolocated data set of churches in the UK. This information is readily and freely available online thanks to the UK Ordnance Survey, a quasi-governmental agency which maintains the various (now digital) maps of Britain. Lucky for us, the Ordnance Survey has an open data product that anyone can use!

+

Before we begin, there are some key things we should note about geospatial data. Geospatial data tends to fall into one of two kinds: points and polygons. Points can be any kind of feature: a house, a church, a pub, someone’s favourite ancient oak tree, or some kind of sacred relic. Polygons tend to be associated with wider areas, and as such can be used to describe large features, e.g. an Ocean, a local authority, or a mountain, or also demographic features, like a census Output Area with associated census summaries. Points are very simple data representations, an X and Y coordinate. Polygons are much more complex, often containing dozens or even thousands of points.

+

The most complex aspect of point data relates to the ways that coordinates are encoded, as they will aways need to be associated with a coordinate reference system (CRS) which describes how they are situated with respect to the planet earth. The most common CRS is the WGS, though for our data sets we’ll also come into contact with the BGS, a specifically British coordinate reference system. There are dozens of CRS, usually mapping onto a specific geographical region. Bearing in mind the way that you need to use a CRS to understand how coordinates can be associated with specific parts of the earth, you can see how this is a bit like survey data, where you also need a “codebook” to understand what the specific response values map onto, e.g. a “1” means “strongly agree” and so on. We also saw, in a previous chapter, how some forms of data have the codebook already baked into the factor data, simplifying the process of interpreting the data. In a similar way, some types of geospatial data sets can also come with CRS “baked in” while we’ll need to define CRS for other types. Here are some of the most common types of geospatial data files:

+

CSV: “comma separated values” a plain text file containing various coordinates Shapefile: a legacy file format, often still in use, but being replaced by others for a variety of good reasons. For more on this see [http://switchfromshapefile.org/] Geopackage: one of the more recent ways of packaging up geospatial data. Geopackages can contain a wide variety of different data and are easily portable. GeoJSON: a file format commonly used in other forms of coding, the “JSON” (an acronym for JavaScript Object Notation) is meant to be easily interchangeable across various platforms. GeoJSON is an augmented version of JSON data with coordinates added in.

+

Now that you have a sense of some of the basic aspects of geospatial data, let’s dive in and do a bit of learning in action.

+
+

5 Administrative shapes - the UK

+

A good starting point is to aquire some adminstrative data. This is a way of referring to political boundaries, whether country borders or those of a local authority or some other administrative unit. For our purposes, we’re going to import several different types of administrative boundary which will be used at different points in our visualisations below. It’s worth noting that the data we use here was prepared to support the 2011 census, and make use of the shapefile format.

+
+
library(sf) |> suppressPackageStartupMessages()
+library(here)  |> suppressPackageStartupMessages()
+library(tidyverse)  |> suppressPackageStartupMessages()
+setwd("/Users/kidwellj/gits/hacking_religion_textbook/hacking_religion")
+here::i_am("chapter_3.qmd")
+
+
here() starts at /Users/kidwellj/gits/hacking_religion_textbook/hacking_religion
+
+
# Download administrative boundaries for whole UK at country level
+if (file.exists(here("data", "infuse_uk_2011_clipped.shp")) == FALSE) {
+download.file("https://borders.ukdataservice.ac.uk/ukborders/easy_download/prebuilt/shape/infuse_uk_2011_clipped.zip", destfile = "data/infuse_uk_2011_clipped.zip")
+unzip("data/infuse_uk_2011_clipped.zip", exdir = "data")
+}
+uk_countries <- st_read(here("data", "infuse_uk_2011_clipped.shp"))
+
+
Reading layer `infuse_uk_2011_clipped' from data source 
+  `/Users/kidwellj/gits/hacking_religion_textbook/hacking_religion/data/infuse_uk_2011_clipped.shp' 
+  using driver `ESRI Shapefile'
+Simple feature collection with 1 feature and 3 fields
+Geometry type: MULTIPOLYGON
+Dimension:     XY
+Bounding box:  xmin: -69.1254 ymin: 5337.9 xmax: 655604.7 ymax: 1220302
+Projected CRS: OSGB36 / British National Grid
+
+
# Download administrative boundaries for whole UK at regions level
+if (file.exists(here("data", "infuse_rgn_2011_clipped.shp")) == FALSE) {
+download.file("https://borders.ukdataservice.ac.uk/ukborders/easy_download/prebuilt/shape/infuse_rgn_2011_clipped.zip", destfile = "data/infuse_rgn_2011_clipped.zip")
+unzip("data/infuse_rgn_2011_clipped.zip", exdir = "data")
+}
+uk_rgn <- st_read(here("data", "infuse_rgn_2011_clipped.shp"))
+
+
Reading layer `infuse_rgn_2011_clipped' from data source 
+  `/Users/kidwellj/gits/hacking_religion_textbook/hacking_religion/data/infuse_rgn_2011_clipped.shp' 
+  using driver `ESRI Shapefile'
+Simple feature collection with 9 features and 2 fields
+Geometry type: MULTIPOLYGON
+Dimension:     XY
+Bounding box:  xmin: 82672 ymin: 5337.9 xmax: 655604.7 ymax: 657534.1
+Projected CRS: OSGB36 / British National Grid
+
+
# Download administrative boundaries for whole UK at local authority level
+if (file.exists(here("data", "infuse_dist_lyr_2011_clipped.shp")) == FALSE) {
+download.file("https://borders.ukdataservice.ac.uk/ukborders/easy_download/prebuilt/shape/infuse_dist_lyr_2011_clipped.zip", destfile = "data/infuse_dist_lyr_2011_clipped.zip")
+unzip("data/infuse_dist_lyr_2011_clipped.zip", exdir = "data")
+}
+local_authorities <- st_read(here("data", "infuse_dist_lyr_2011_clipped.shp"))
+
+
Reading layer `infuse_dist_lyr_2011_clipped' from data source 
+  `/Users/kidwellj/gits/hacking_religion_textbook/hacking_religion/data/infuse_dist_lyr_2011_clipped.shp' 
+  using driver `ESRI Shapefile'
+Simple feature collection with 404 features and 3 fields
+Geometry type: MULTIPOLYGON
+Dimension:     XY
+Bounding box:  xmin: -69.1254 ymin: 5337.9 xmax: 655604.7 ymax: 1220302
+Projected CRS: OSGB36 / British National Grid
+
+
# Download building outlines for whole UK
+if (file.exists(here("data", "infuse_dist_lyr_2011_simplified_100m_buildings_simplified.gpkg")) == FALSE) {
+  download.file("https://zenodo.org/record/6395804/files/infuse_dist_lyr_2011_simplified_100m_buildings_overlay_simplified.gpkg", destfile = here("data", "infuse_dist_lyr_2011_simplified_100m_buildings_simplified.gpkg"))}
+local_authorities_buildings_clip <- st_read(here("data", "infuse_dist_lyr_2011_simplified_100m_buildings_simplified.gpkg"))
+
+
Reading layer `infuse_dist_lyr_2011_simplified_100_buildings_overlay_simplified' from data source `/Users/kidwellj/gits/hacking_religion_textbook/hacking_religion/data/infuse_dist_lyr_2011_simplified_100m_buildings_simplified.gpkg' 
+  using driver `GPKG'
+Simple feature collection with 403 features and 0 fields
+Geometry type: MULTIPOLYGON
+Dimension:     XY
+Bounding box:  xmin: -69.1254 ymin: 5524.797 xmax: 655986.4 ymax: 1219597
+Projected CRS: OSGB36 / British National Grid
+
+
+

Before we move on, let’s plot a simple map and have a look at one of our administrative layers. We can use ggplot with a new type of shape geom_sf() to plot the contents of a geospatial data file with polygons which is loaded as a simplefeature in R.

+
+
ggplot(uk_countries) +
+  geom_sf()
+
+

+
+
+
+
+

6 Load in Ordnance Survey OpenMap Points Data

+
+
# Note: for more advanced reproducible scripts which demonstrate how these data surces have been 
+# obtained, see the companion cookbook here: https://github.com/kidwellj/hacking_religion_cookbook/blob/main/ordnance_survey.R
+
+os_openmap_pow <- st_read(here("data", "os_openmap_pow.gpkg"))
+
+
Reading layer `os_openmap_pow' from data source 
+  `/Users/kidwellj/gits/hacking_religion_textbook/hacking_religion/data/os_openmap_pow.gpkg' 
+  using driver `GPKG'
+Simple feature collection with 48759 features and 5 fields
+Geometry type: POLYGON
+Dimension:     XY
+Bounding box:  xmin: 64594.12 ymin: 8287.54 xmax: 655238.1 ymax: 1214662
+Projected CRS: OSGB36 / British National Grid
+
+
ggplot(os_openmap_pow) +
+  geom_sf()
+
+

+
+
+

It’s worth noting that the way that you load geospatial data in R has changed quite dramatically since 2020 with the introduction of the simplefeature class in R. Much of the documentation you will come across “out there” will make reference to a set of functions which are now deprecated.

+

Let’s use that data we’ve just loaded to make our first map:

+
+
# Generate choropleth map of respondent locations
+# using temporary palette here so that 0s are white
+library(tmap) |> suppressPackageStartupMessages()
+# palette <- c(white, "#a8ddb5", "#43a2ca")
+map1 <- tm_shape(local_authorities) + 
+#  tm_fill(col = "surveys_count", , palette = palette, title = "Concentration of survey respondents") +
+  tm_borders(alpha=.5, lwd=0.1) +
+  # for intermediate polygon geometries
+  # tm_shape(local_authorities) +
+  # tm_borders(lwd=0.6) +
+  # for dots from original dataset
+  # tm_dots("red", size = .05, alpha = .4) +
+  tm_scale_bar(position = c("right", "bottom")) +
+  tm_style("gray") +
+  tm_credits("Data: UK Data Service (OGL)\n& Jeremy H. Kidwell,\nGraphic is CC-by-SA 4.0", 
+             size = 0.4, 
+             position = c("left", "bottom"),
+             just = c("left", "bottom"),
+             align = "left") +
+  tm_layout(asp = NA,
+            frame = FALSE, 
+            title = "Figure 1a", 
+            title.size = .7,
+            legend.title.size = .7,
+            inner.margins = c(0.1, 0.1, 0.05, 0.05)
+  )
+
+map1
+
+

+
+
# save image
+tmap_save(map1, here("figures", "map.png"), width=1920, height=1080, asp=0)
+
+
Map saved to /Users/kidwellj/gits/hacking_religion_textbook/hacking_religion/figures/map.png
+
+
+
Resolution: 1920 by 1080 pixels
+
+
+
Size: 6.4 by 3.6 inches (300 dpi)
+
+
+
+
# subsetting ordnance survey openmap data for measuring clusters and proximity
+
+os_openmap_important_buildings <- st_read(here("data", "os_openmap_important_buildings.gpkg"))
+
+
Reading layer `important_buildings' from data source 
+  `/Users/kidwellj/gits/hacking_religion_textbook/hacking_religion/data/os_openmap_important_buildings.gpkg' 
+  using driver `GPKG'
+Simple feature collection with 229800 features and 5 fields
+Geometry type: POLYGON
+Dimension:     XY
+Bounding box:  xmin: 64594.12 ymin: 8125.44 xmax: 655500.5 ymax: 1214662
+Projected CRS: OSGB36 / British National Grid
+
+
# add pubs, check_cashing, pawnbrokers, SSSI
+## subsets
+
+
+
# OSM data
+
+# Note: for more advanced reproducible scripts which demonstrate how these data surces have been 
+# obtained, see the companion cookbook here: https://github.com/kidwellj/hacking_religion_cookbook/blob/main/ordnance_survey.R
+
+
+# osm_uk_points <- st_read(system.file(here("data", "pow_osm.gpkg", package = "spData")))
+# vector_filepath = system.file("data/osm-gb-2018Aug29_pow_osm.pbf", package = "sf")
+# osm_uk_points = st_read(vector_filepath)
+

Guides to geographies: https://rconsortium.github.io/censusguide/ https://ocsi.uk/2019/03/18/lsoas-leps-and-lookups-a-beginners-guide-to-statistical-geographies/

Extact places of worship from Ordnance survey open data set Calculate proximity to pubs

+

References