The following code chunks load the relevant libraries, and create subset indices for Belfast and Derry. Click the Show buttons to the right to view the relevant code for each section.

# Loading Libraries and Converting Data
library(tmaptools)
library(GISTools)
library(tidyverse)
library(sp)
library(spdep)
library(tmap)
library(dplyr)
library(OpenStreetMap)
library(grid)
library(readxl)
library(corrplot)
library(ggcorrplot)
library(bslib)
library(rJava)
setwd("C:/Users/eflaherty/iCloudDrive/Papers/NI census/Census Boundaries/Super Data Zone")
ni_sdv <- st_read("C:/Users/eflaherty/iCloudDrive/Papers/NI census/Census Boundaries/Super Data Zone/sdz_merge.shp")
ni_sf = st_as_sf(ni_sdv)
ni_sp <- as(ni_sf, "Spatial")

We also include some code to generate some OpenStreetMap background layers for context.

# Creating Subset Indices for Belfast and Derry
index.belfast <- c(161:335)
ni_belfast.sub <- ni_sdv[index.belfast,]
index.derry <- c(397:435)
ni_derry.sub <- ni_sdv[index.derry,]

# We will also create an OpenStreetMap layer for context using the OpenStreetMap package. The new OSM layers will be called 'osm_belfast, osm_derry, and osm_ni'.
library(OpenStreetMap)

osm_belfast <- read_osm(ni_belfast.sub)
tm_shape(osm_belfast) + tm_rgb()

osm_derry <- read_osm(ni_derry.sub)
tm_shape(osm_derry) + tm_rgb()

osm_ni <- read_osm(ni_sdv)
tm_shape(osm_ni) + tm_rgb()

Introduction

The figures and tables below explore patterns of non-response to the religion question in the census of population in Northern Ireland. The decennial census waves of interest are those of 2001, 2011, and 2021. The analysis below focuses on 2021 for now, pending further preparation of data from preceding years. Our preceding paper (Coulter, Shirlow, and Flaherty 2023) examined the ‘bio-politics’ of the Northern Irish Census, and in doing so, suggested an anomaly in the presumed count of a Catholic majority. Leaving aside the accuracy of the final recorded count (which we do not contest), we focus instead on the processes that may have led to an increase in those now opting to report their religion, who may previously have refused. There is precedent for this in Northern Ireland, with its history of organised census boycotts, particularly in 1971 on the back of the Civil Rights Movement, and again in 1981 following hunger strikes among Republican prisoners. The effects of these boycotts are still being uncovered through historical research (Cooley 2024). We suggest that something of this may be detectable in the data, more specifically the geography of religion and non-response, and the characteristics of those areas where changes in response patterns may be strongest.

Those reporting as ‘no religion’ are predominantly (ethnic) Protestants, as is visible in the pattern of responses. In Belfast, the distribution of these responses is characteristically West-East, and those areas with the lowest rates of ‘no religion’ are also those where more current Catholics reside. Those recorded as ‘not stated’ are more diffuse, and these are the rates of interest. One working hypothesis suggests that the effect may be limited to specific areas where compliance with past census boycotts may have been highest (acknowledging that this is difficult in itself to infer). If this is true, and this is reflected in the data on births, then the increase in Catholics in 2021 is more to do with greater affirmative response to the religion question, net of Catholic births and Protestant out-migration. The specific characteristics of areas with the following characteristics needs to be investigated further: highly religious areas (high levels of affirmative reporting of a specific religion), those with high levels of specific ethnic identification (i.e. British/Irish), and those with high levels of non-response. This analysis should proceed as follows:

  1. Analysis of the general geography of non-response across NI, and within Belfast and Derry.
  2. Local analysis of specific clusters or outliers.
  3. Cluster analysis to include as input conditions: religion (possibly ethno-national identity to counter colinearlity), non-response, and socio-economic conditions (education levels).
  4. Time-series comparison of census waves using NISRA Grid Square data.
  5. Geographically weighted modelling.

Basic Plots for ‘Religion Not Stated’

tm_shape(ni_sdv) +
  tm_fill("per_nsta", title = "Rlgn Not Stated (2021)", palette = "Blues", 
  breaks = c(0, round(quantileCuts(ni_sdv$per_nsta, 6), 2))) +
  tm_borders(lty = "solid", col = "black") +
  tm_style("natural", bg.color = "grey90") + 
  tm_borders(lwd = 1) +
  tm_scale_bar(width = 0.22, position = c(0.05, 0.02)) +
  tm_compass(position = c(0.9, 0.07)) +
  tm_layout(legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","top"),
            legend.format = list(text.separator = "-"),
            legend.bg.alpha = 1,
            title = "Northern Ireland",
            title.size = 1,
            title.position = c(0.35, "top"))

Figure 1. Distribution of ‘Religion Not Stated’ (NI 2021)

tm_shape(ni_belfast.sub) +
  tm_fill("per_nsta", alpha = 0.9, title = "Rlgn Not Stated (2021)", palette = "Blues", 
          breaks = c(0, round(quantileCuts(ni_belfast.sub$per_nsta, 6), 2))) +
  tm_borders(lty = "solid", col = "black") +
  tm_style("natural", bg.color = "grey90") + 
  tm_borders(lwd = 1) +
  tm_scale_bar(width = 0.22) +
  tm_compass(position = c(0.8, 0.09)) +
  tm_layout(legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","top"),
            legend.format = list(text.separator = "-"),
            legend.bg.alpha = 1,
            title = "Belfast",
            title.size = 1,
            title.position = c("0.7", "top")) +
  tm_shape(osm_belfast) + tm_rgb(alpha = 0.4)

Figure 2. Distribution of ‘Religion Not Stated’ (Belfast 2021)

tm_shape(ni_derry.sub) +
  tm_fill("per_nsta", alpha = 0.9, title = "Rlgn Not Stated (2021)", palette = "Blues", 
          breaks = c(0, round(quantileCuts(ni_derry.sub$per_nsta, 6), 2))) +
  tm_borders(lty = "solid", col = "black") +
  tm_style("natural", bg.color = "grey90") + 
  tm_borders(lwd = 1) +
  tm_scale_bar(width = 0.22) +
  tm_compass(position = c(0.8, 0.09)) +
  tm_layout(legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","top"),
            legend.format = list(text.separator = "-"),
            legend.bg.alpha = 1,
            title = "Derry",
            title.size = 1,
            title.position = c(0.6, "top")) +
  tm_shape(osm_derry) + tm_rgb(alpha = 0.4)

Figure 3. Distribution of ‘Religion Not Stated’ (Derry 2021)

Visualisations

# assign some variables
p.cat <- ni_sp$per_cat
p.nrel <- ni_sp$per_nrel
p.nsta <- ni_sp$per_nsta
# bind these together
df <- data.frame(p.cat, p.nrel, p.nsta)
# fit regressions
mod.1.2 <- lm(p.nrel ~ p.cat, data = df)
mod.2.2 <- lm(p.nsta ~ p.cat, data = df)
# Plot regression results
p1 <- ggplot(df,aes(p.cat, p.nrel))+
  #stat_summary(fun.data=mean_cl_normal) +
  geom_smooth(method='lm') +
  geom_point(alpha = 0.2) +
  xlab("Percentage Catholic") +
  ylab("Percentage No Religion") +
  labs(title="No Rlgn by Catholic (NI 2021)")

# Plot regression results
p2 <- ggplot(df,aes(p.cat, p.nsta))+
  #stat_summary(fun.data=mean_cl_normal) +
  geom_smooth(method='lm') +
  geom_point(alpha = 0.2) +
  xlab("Percentage Catholic") +
  ylab("Percentage Rlgn Not Stated") +
  labs(title="Not Stated by Catholic (NI 2021)")

# Belfast only
index.belfast <- c(161:335)
ni_belfast.sub <- ni_sdv[index.belfast,]
# sp conversion
bf_sf = st_as_sf(ni_belfast.sub)
bf_sp <- as(bf_sf, "Spatial")

# assign some variables
pbf.cat <- bf_sp$per_cat
pbf.nrel <- bf_sp$per_nrel
pbf.nsta <- bf_sp$per_nsta
# bind these together
df.belfast <- data.frame(pbf.cat, pbf.nrel, pbf.nsta)
# fit regressions
mod.1.3 <- lm(pbf.nrel ~ pbf.cat, data = df)
mod.2.3 <- lm(pbf.nsta ~ pbf.cat, data = df)

# Plot regression results
p3 <- ggplot(df.belfast,aes(pbf.cat, pbf.nrel))+
  #stat_summary(fun.data=mean_cl_normal) +
  geom_smooth(method='lm') +
  geom_point(alpha = 0.2) +
  xlab("Percentage Catholic") +
  ylab("Percentage No Religion") +
  labs(title="No Rlgn by Catholic (Belfast 2021)")

# Plot regression results
p4 <- ggplot(df.belfast,aes(pbf.cat, pbf.nsta))+
  #stat_summary(fun.data=mean_cl_normal) +
  geom_smooth(method='lm') +
  geom_point(alpha = 0.2) +
  xlab("Percentage Catholic") +
  ylab("Percentage Rlgn Not Stated") +
  labs(title="Not Stated by Catholic (Belfast 2021)")

grid.newpage()
# set up the layout
pushViewport(viewport(layout=grid.layout(2,2)))
# plot using the print command
print(p1, vp=viewport(layout.pos.row = 1, layout.pos.col = 1, height = 5))
print(p2, vp=viewport(layout.pos.row = 1, layout.pos.col = 2, height = 5))
print(p3, vp=viewport(layout.pos.row = 2, layout.pos.col = 1, height = 5))
print(p4, vp=viewport(layout.pos.row = 2, layout.pos.col = 2, height = 5))

Figure 4. No Religion and Not Stated by Catholic Population (2021)

sum_table <- read_excel("C:/Users/eflaherty/iCloudDrive/Papers/NI census/Census Boundaries/Super Data Zone/summary_table_1.xlsx")
knitr::kable(sum_table,
 caption = "Summary Statistics (NI 2021")
Summary Statistics (NI 2021
Variable Min Max Mean Median
% Catholic 2.02 95.02 40.48 34.04
% Protestant 1.21 83.08 38.30 42.88
% Not Stated 0.21 14.62 1.62 1.41
% No Religion 2.44 43.24 18.19 17.31
% Irish 0.94 83.07 27.74 19.42
% British 2.11 68.72 32.83 35.90
% Northern Irish 6.88 39.46 19.76 20.00
% Level 4 + 8.15 60.49 25.38 24.41

Spatial Lag Plots

# List details of neighbour connections (Queen)
ni.nb <- poly2nb(ni_sp)
# Convert neighbour list to listw object (queen's case)
ni_nb_lw <- nb2listw(ni.nb)
ni_sp$per_nsta.lagged.means <- lag.listw(ni_nb_lw,ni_sp$per_nsta)
tm_shape(ni_sp) + tm_polygons(col= 'per_nsta.lagged.means',title= '% Not Stated.') +
  tm_scale_bar(width = 0.22, position = c("left", "bottom")) +
  tm_compass(position = c(0.07, 0.08)) +
  tm_layout(legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","top"),
            legend.format = list(text.separator = "-"),
            legend.bg.alpha = 1,
            title = "N.I. Spatial Lag Plot (2021)",
            title.size = 1,
            title.position = c(0.2, "top"))

Spatial Lag Map for ‘Religion Not Stated’ (NI, 2021)

## 
##  Moran I test under randomisation
## 
## data:  ni_sp$per_nsta  
## weights: ni_nb_lw    
## 
## Moran I statistic standard deviate = 11.49, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.2407434697     -0.0011778563      0.0004433108

Moran’s I Table

# List details of neighbour connections (Queen)
ni.belfast <- poly2nb(ni_belfast.sub)
# Convert neighbour list to listw object (queen's case)
ni_belfast_lw <- nb2listw(ni.belfast)
ni_belfast.sub$per_nsta.lagged.means <- lag.listw(ni_belfast_lw,ni_belfast.sub$per_nsta)
tm_shape(ni_belfast.sub) + tm_polygons(alpha = 0.9, col= 'per_nsta.lagged.means',title= '% Not Stated.') +
  tm_scale_bar(width = 0.22) +
  tm_compass(position = c(0.8, 0.09)) +
  tm_layout(legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","top"),
            legend.format = list(text.separator = "-"),
            legend.bg.alpha = 1,
            title = "Belfast Spatial Lag Plot (2021)",
            title.size = 1,
            title.position = c(0.5, "top")) +
   tm_shape(osm_belfast) + tm_rgb(alpha = 0.4)

Spatial Lag Map for ‘Religion Not Stated’ (Belfast, 2021)

## 
##  Moran I test under randomisation
## 
## data:  bf_sp$per_nsta  
## weights: ni_belfast_lw    
## 
## Moran I statistic standard deviate = 4.6229, p-value = 1.892e-06
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.187585571      -0.005747126       0.001748992

Moran’s I Table

with(data.frame(ni_sp), {
  plot(per_nsta,per_nsta.lagged.means, asp= 1, xlim=range(0,5), ylim=range(0,6), xlab="% religion not stated", ylab="Lagged mean")
  abline(a= 0, b= 1)
  abline(v=mean(per_nsta), lty= 2)
  abline(h=mean(per_nsta.lagged.means), lty= 2)
})

Lagged Mean Plot (NI 2021)

with(data.frame(ni_belfast.sub), {
  plot(per_nsta,per_nsta.lagged.means, asp= 1, xlim=range(0, 5), ylim=range(0, 6), 
       xlab="% religion not stated", ylab="Lagged mean")
  abline(a= 0, b= 1)
  abline(v=mean(per_nsta), lty= 2)
  abline(h=mean(per_nsta.lagged.means), lty= 2)
})

Lagged Mean Plot (Belfast 2021)

Local analysis plots

# Belfast - generate neighbour list
bf.nb <- poly2nb(ni_belfast.sub)
# Convert neighbour list to listw object (queen's case)
bf_nb_lw <- nb2listw(bf.nb)

# Compute the local Moran's I for the Belfast layer
ni_belfast.sub$lI <- localmoran(ni_belfast.sub$per_nsta,bf_nb_lw)[, 1]

# Compute and filter for all locations with values >=1
lisabelfast <- ni_belfast.sub[ni_belfast.sub$lI >= 1, ]
lisabelfast_sf = st_as_sf(lisabelfast)
lisabelfast_sp <- as(lisabelfast_sf, "Spatial")

# Draw the map with the 'lisabelfast' layer called for name information
tm_shape(ni_belfast.sub) +
  tm_polygons(col= 'lI', title= "Local Moran's I",legend.format=list(flag= "+"), border.alpha = 0.3) +
  tm_scale_bar(width=0.15) +
  tm_compass(position = c(0.8, 0.09)) +
  tm_style('col_blind') + tm_scale_bar(width= 0.15) +
  tm_layout(legend.position = c("right", "top"),
            legend.text.size= .6) +
  tm_shape(lisabelfast_sp) + tm_text("SDZ2021_nm", col="black", size = 0.8) +
  tm_shape(osm_belfast) + tm_rgb(alpha = 0.3)

The above map identifies and labels those areas indicative of high-value clustering using Local Moran’s I

# Compute the local p-values for the Belfast layer
ni_belfast.sub$pval <- localmoran(ni_belfast.sub$per_nsta,bf_nb_lw)[, 1]

# Compute and filter for all locations with values >0.05
lisapbelfast <- ni_belfast.sub[ni_belfast.sub$pval >= 0.05, ]
lisapbelfast_sf = st_as_sf(lisapbelfast)
lisapbelfast_sp <- as(lisapbelfast_sf, "Spatial")

# Map the p-value
tm_shape(ni_belfast.sub) +
tm_polygons(col= 'pval', title = "p-value", breaks= c(0, 0.01, 0.05, 0.10, 1), border.alpha = 0.3) +
  tm_scale_bar(width=0.15) +
  tm_compass(position = c(0.8, 0.09)) +
  tm_layout(legend.position = c("right", "top")) +
   tm_shape(lisapbelfast_sp) + tm_text("SDZ2021_nm", col="black", size = 0.4) +
   tm_shape(osm_belfast) + tm_rgb(alpha = 0.3)

The above plots statistically significant locations at .05 level.

The plots below explore the relationship between secularism and education more closely, first for the whole sample, and then by local authority district.

# Whole region
plot2.1 <- ggplot(data = ni_sf,
            mapping = aes(x = per_l4,
                          y = per_b_non)) +
  labs(x = "% level 4 education and above", y = "% brought up with no religion",
         title = "Secularism and Education (Northern Ireland)",
         subtitle = "Data points are SDZs",
         caption = "Source: Census 2021")
plot2.1 + geom_point(color = "black", alpha = 0.5, size = 1) + 
  geom_smooth(color = "blue", method = "lm") + geom_smooth(color = "orange", method = "gam")

# Split by LGD
plot2.2 <- ggplot(data = ni_sf,
            mapping = aes(x = per_l4,
                          y = per_b_non)) +
  labs(x = "% level 4 education and above", y = "% brought up with no religion",
         title = "Secularism and Education (Northern Ireland)",
         subtitle = "Data points are SDZs",
         caption = "Source: Census 2021")
plot2.2 + geom_point(color = "black", alpha = 0.5, size = 1) + 
  geom_smooth(color = "blue", method = "lm", size = 0.5) + geom_smooth(color = "orange", method = "gam", size = 0.5) +
  facet_wrap(~ LGD2014_nm, ncol=4) +
  theme(
    axis.text = element_text( size = 7),
    axis.text.x = element_text( size = 7),
    axis.title = element_text( size = 7, face = "bold"),
    strip.text = element_text(size = 7)
  )

Distribution of those raised with no religion

tm_shape(ni_belfast.sub) +
  tm_fill("per_b_non", alpha = 0.95,  title = "% Brought up No Religion (2021)", palette = "Blues", 
          breaks = c(0, round(quantileCuts(ni_belfast.sub$per_b_non, 6), 2))) +
  tm_borders(lty = "solid", col = "black") +
  tm_style("natural", bg.color = "grey90") + 
  tm_borders(lwd = 1) +
  tm_scale_bar(width = 0.22) +
  tm_compass(position = c(0.8, 0.09)) +
  tm_layout(legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","top"),
            legend.format = list(text.separator = "-"),
            legend.bg.alpha = 1,
            title = "Belfast",
            title.size = 1,
            title.position = c(0.6, "top")) +
   tm_shape(osm_belfast) + tm_rgb(alpha=.5)

Next we measure and plot correlations between variables in our dataset.

# First read in data as excel sheet only
sdz_master_data_entry_2021 <- read_excel("C:/Users/eflaherty/iCloudDrive/Papers/NI census/sdz_master_data_entry_2021_transferonly.xlsx")
names(sdz_master_data_entry_2021)
##  [1] "SDZ2021_cd"  "SDZ2021_nm"  "population"  "catholic"    "presby"     
##  [6] "cofireland"  "methodist"   "oth_cst"     "protestant"  "oth_rel"    
## [11] "no_relig"    "not_stated"  "per_cat"     "per_pro"     "per_nrel"   
## [16] "per_nsta"    "brou_clic"   "brou_pro"    "brou_other"  "brou_none"  
## [21] "per_b_cat"   "per_b_pro"   "per_b_oth"   "per_b_non"   "d_rate_cat" 
## [26] "d_rate_prot" "d_rate_oth"  "d_rate_none" "diff_cath"   "diff_prot"  
## [31] "diff_other"  "diff_none"   "per_bcath"   "per_bpro"    "per_both"   
## [36] "per_bnon"    "per_tism"    "per_iri"     "per_brit"    "per_ni"     
## [41] "per_l4"
# Then subset with only the variables of interest
ni_corr_subset <- sdz_master_data_entry_2021[, c("per_cat", "per_pro", "per_nrel", "per_nsta", "per_b_oth", "per_b_non", "d_rate_none", "per_iri", "per_ni")]
names(ni_corr_subset)
## [1] "per_cat"     "per_pro"     "per_nrel"    "per_nsta"    "per_b_oth"  
## [6] "per_b_non"   "d_rate_none" "per_iri"     "per_ni"
# Try the ggcorrplot method also - first generate a matrix for r...
corr_ni_r <- round(cor(ni_corr_subset), 1)
head(corr_ni_r[, 1:6])
##           per_cat per_pro per_nrel per_nsta per_b_oth per_b_non
## per_cat       1.0    -0.9     -0.8     -0.1      -0.2      -0.8
## per_pro      -0.9     1.0      0.6     -0.1      -0.1       0.5
## per_nrel     -0.8     0.6      1.0      0.1       0.4       1.0
## per_nsta     -0.1    -0.1      0.1      1.0       0.4       0.2
## per_b_oth    -0.2    -0.1      0.4      0.4       1.0       0.4
## per_b_non    -0.8     0.5      1.0      0.2       0.4       1.0
# ...and p-values
p.corr_ni_r <- cor_pmat(ni_corr_subset)
head(p.corr_ni_r[, 1:4])
##                 per_cat      per_pro      per_nrel     per_nsta
## per_cat    0.000000e+00 0.000000e+00 6.130451e-178 1.300118e-01
## per_pro    0.000000e+00 0.000000e+00  1.738867e-68 4.438416e-02
## per_nrel  6.130451e-178 1.738867e-68  0.000000e+00 1.184672e-04
## per_nsta   1.300118e-01 4.438416e-02  1.184672e-04 0.000000e+00
## per_b_oth  5.005472e-06 5.330226e-02  7.683950e-34 8.310815e-39
## per_b_non 7.416313e-167 6.781139e-66  0.000000e+00 2.045151e-11
# Then plot:
ggcorrplot(corr_ni_r, hc.order = TRUE, type = "lower",
     outline.col = "white",  lab = TRUE, method = "square", p.mat = p.corr_ni_r)

Let’s also try some preliminary modelling. For now this is principally bivariate linear, with some controls presented visually.

Preliminary Analysis of Autism Data

# Whole region
plot2.3 <- ggplot(data = ni_sf,
            mapping = aes(x = per_l4,
                          y = per_tism)) +
  labs(x = "% level 4 education and above", y = "% reporting Autism / Aspergers",
         title = "Autism and Education (Northern Ireland, 2021)",
         subtitle = "r = -.49",
         caption = "Source: Census 2021")
plot2.3 + geom_point(color = "black", alpha = 0.5, size = 1) + 
  geom_smooth(color = "blue", method = "lm") + geom_smooth(color = "orange", method = "gam")

# Split by LGD
plot2.4 <- ggplot(data = ni_sf,
            mapping = aes(x = per_l4,
                          y = per_tism)) +
  labs(x = "% level 4 education and above", y = "% reporting Autism / Aspergers",
         title = "Autism and Education by Local Government District (Northern Ireland, 2021)",
         subtitle = "Data points are SDZs",
         caption = "Source: Census 2021")
plot2.4 + geom_point(color = "black", alpha = 0.5, size = 1) + 
  geom_smooth(color = "blue", method = "lm", size = 0.5) + geom_smooth(color = "orange", method = "gam", size = 0.5) +
  facet_wrap(~ LGD2014_nm, ncol=4) +
  theme(
    axis.text = element_text( size = 7),
    axis.text.x = element_text( size = 7),
    axis.title = element_text( size = 7, face = "bold"),
    strip.text = element_text(size = 7),
  )

What does this look like in urban areas? Let’s start with Belfast, where the relationship is especailyl strong. The difference in correlation between education levels and autism identification for the whole-sample (r = -.49) and Belfast (r = -.79) is striking. Isolating the scatterplot for Belfast we see the following:

plot2.5 <- ggplot(data = bf_sf,
            mapping = aes(x = per_l4,
                          y = per_tism)) +
  labs(x = "% level 4 education and above", y = "% reporting Autism / Aspergers",
         title = "Autism and Education (Belfast, 2021)",
         subtitle = "r = -.79",
         caption = "Source: Census 2021")
plot2.5 + geom_point(color = "black", alpha = 0.5, size = 1) + 
  geom_smooth(color = "blue", method = "lm") + geom_smooth(color = "orange", method = "gam")

tm_shape(ni_belfast.sub) +
  tm_fill("per_tism", title = "Reported Autism / Aspergers (2021)", palette = "Blues", 
          breaks = c(0, round(quantileCuts(ni_belfast.sub$per_tism, 6), 2))) +
  tm_borders(lty = "solid", col = "black") +
  tm_style("natural", bg.color = "grey90") + 
  tm_borders(lwd = 1) +
  tm_scale_bar(width = 0.22) +
  tm_compass(position = c(0.8, 0.09)) +
  tm_layout(legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","top"),
            legend.format = list(text.separator = "-"),
            legend.bg.alpha = 1,
            title = "Belfast",
            title.size = 1,
            title.position = c(0.6, "top")) +
 tm_shape(osm_belfast) + tm_rgb(alpha=.5)

Finally, examine the correlation coefficients between education and autism for Belfast and NI

tism_edu_bf <- cor(ni_belfast.sub$per_tism, ni_belfast.sub$per_l4, method = 'pearson')
tism_edu_bf
## [1] -0.7855008
tism_edu_ni <- cor(ni_sf$per_tism, ni_sf$per_l4, method = 'pearson')
tism_edu_ni
## [1] -0.4911214