Wolfe Lab Team Meeting (SISMID Learnings)

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document.

The blocks in the light gray-ish color are code chunks. Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

The Basics and Basic Maps in R

Maps of locations linked to a table of attributes.

Geographic information system (GIS) “a technology designed to capture, store, manipulate, analyze, and visualize georeferenced data” (Goodchild, Parks, and Steyaert 1993). GIS is a database system containing locations for every value and allowing operations (search, sorting, etc) based on locations as well as attributes.

GIS Operation 1: Layering

Datasets ~ “Layers” (like a cake). Examples: site/case locations (points), road locations (lines), population levels (areas), satellite image(raster)

GIS Operation 2: Buffering

Find areas within a user specified distance.This can be done with points, lines, and areas. (ex: you can show like what falls within a 25 mi of a WWTP site & see if the buffers overlap).

GIS Operation 3: Joining What you have in spatial joins: attribute table linked to a map, 2nd table of data over the same features, common identifier for both What you want doing a spatial join: add (join) attributes in 2nd table to first table by linking tables based on common attributes SO you need a 1-to-1 correspondence

There are different softwares you can use for mapping and spatial analysis. From what I understand, there are statistical limitations with the number and specificity of toolboxes in GIS.

ArcGIS (need money): GIS/mapping software by Esri, geodatabases (ArcGIS’s native container/folder/database format), shapefiles (older too but still widely used vector data format) that come in a collection of files (.shp, .dbf, .shx, etc.)

SAS (need money): Statistical software heavily used in emory coursework/cdc/coporate/etc. SAS datasets (.sas7bdat)

R: Free statistical programming free language, objects (which is everything in R - vectors, lists, models, plots, functions), dataframes with rows as observations and columns as variables that can be numeric, characters etc. (tables/data structure in R )

Source: Emory SISMID 2026

Again, R is open-source. There are many old and new spatial anlysis packages available. sp package is the former standard for spatial processing. sf package is the new standard.

Installing Packages

We first need to install R packages that will contain tools, functions, datasets that expand the capabilities of base R programming language. Documentation for packages can be found online.

# Install missing packages if you need them. Don't if you got them already (check packages tab at the bottom right). #
#install.packages <- (c("ggplot2", "devtools", "dplyr", "stringr"))

#older/standard map packages
#install.packages(c("maps","mapdata"))

#install.packages('ggmap')

Loading Libraries

Now the next step is to load packages in r using the library () function.

library(ggplot2) #provides the map_data() function that turns a series of points along an outline into a data fram of those points
library(ggmap)

## ℹ Google's Terms of Service: <https://mapsplatform.google.com>
##   Stadia Maps' Terms of Service: <https://stadiamaps.com/terms-of-service>
##   OpenStreetMap's Tile Usage Policy: <https://operations.osmfoundation.org/policies/tiles>
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.

library(maps) #has lots of outlines for continents, countries, states, and counties
library(mapdata) #higher-resolution outlines 

Maps package maps with ggplot

usa <- map_data ("usa")


#doing some inspection
dim(usa) #get dimensions of an R object

## [1] 7243    6

head(usa) #view first 6 rows

##        long      lat group order region subregion
## 1 -101.4078 29.74224     1     1   main      <NA>
## 2 -101.3906 29.74224     1     2   main      <NA>
## 3 -101.3620 29.65056     1     3   main      <NA>
## 4 -101.3505 29.63911     1     4   main      <NA>
## 5 -101.3219 29.63338     1     5   main      <NA>
## 6 -101.3047 29.64484     1     6   main      <NA>

tail(usa) #view last 6 rows

##           long      lat group order         region subregion
## 7238 -122.6187 48.37482    10  7247 whidbey island      <NA>
## 7239 -122.6359 48.35764    10  7248 whidbey island      <NA>
## 7240 -122.6703 48.31180    10  7249 whidbey island      <NA>
## 7241 -122.7218 48.23732    10  7250 whidbey island      <NA>
## 7242 -122.7104 48.21440    10  7251 whidbey island      <NA>
## 7243 -122.6703 48.17429    10  7252 whidbey island      <NA>

Simple black map

usa <- map_data ("usa")

ggplot() + geom_polygon(data = usa, aes(x=long, y = lat, group = group)) + 
  coord_fixed (1.3)

Simple map with no fill using fill = NA

ggplot() + 
  geom_polygon(data = usa, aes(x=long, y = lat, group = group), fill = NA, color = "purple") + coord_fixed (1.3)

Simple map with fill and outline

gg1 <- ggplot() + 
  geom_polygon(data = usa, aes(x=long, y = lat, group = group), fill = "purple", color = "orange") + coord_fixed (1.3)

#group here in gg1 is super important for the next step (lifts the pen!)

gg1

#ga_wwscan <- data.frame(
    #long = c(), #combining longitude into a single vector
    #lat = c(), #combining latitude into a single vector
    #names = c("BC", "CC", "SC", "JC", "LR", "RMC", "SR", "UC"), 
    #stringsAsFactors = FALSE) #ensures text column doesn't get converted into a category when readind building this dataframe

#gg1 +
  #geom_point(data = ga_wwscan, aes(x=long, y = lat), color = "black", size = 5) + geom_point (data = ga_wwscan,   aes(x=long, y = lat), color = "white", size = 4)

US maps with states

states <- map_data("state")

dim(states)

## [1] 15537     6

head(states)

##        long      lat group order  region subregion
## 1 -87.46201 30.38968     1     1 alabama      <NA>
## 2 -87.48493 30.37249     1     2 alabama      <NA>
## 3 -87.52503 30.37249     1     3 alabama      <NA>
## 4 -87.53076 30.33239     1     4 alabama      <NA>
## 5 -87.57087 30.32665     1     5 alabama      <NA>
## 6 -87.58806 30.32665     1     6 alabama      <NA>

tail(states)

##            long      lat group order  region subregion
## 15532 -106.3295 41.00659    63 15594 wyoming      <NA>
## 15533 -106.8566 41.01232    63 15595 wyoming      <NA>
## 15534 -107.3093 41.01805    63 15596 wyoming      <NA>
## 15535 -107.9223 41.01805    63 15597 wyoming      <NA>
## 15536 -109.0568 40.98940    63 15598 wyoming      <NA>
## 15537 -109.0511 40.99513    63 15599 wyoming      <NA>

#plot all states colored differently

ggplot(data = states) + 
  geom_polygon(aes(x = long, y = lat, fill = region, group = group), color = "white") +
  coord_fixed(1.3) + 
  guides (fill = FALSE) #no color legend

## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Zoom in at Georgia and Counties

Get GA data only in a dataframe

ga_df <- subset(states, region == "georgia")

head(ga_df)

##           long      lat group order  region subregion
## 2326 -80.89018 32.03980    10  2335 georgia      <NA>
## 2327 -80.85007 32.02834    10  2336 georgia      <NA>
## 2328 -80.84435 32.01688    10  2337 georgia      <NA>
## 2329 -80.85007 31.99396    10  2338 georgia      <NA>
## 2330 -80.88445 31.97105    10  2339 georgia      <NA>
## 2331 -80.93029 31.95959    10  2340 georgia      <NA>

tail(ga_df)

##           long      lat group order  region subregion
## 2701 -81.09072 32.09710    10  2710 georgia      <NA>
## 2702 -81.05061 32.09137    10  2711 georgia      <NA>
## 2703 -80.99905 32.08564    10  2712 georgia      <NA>
## 2704 -80.96467 32.06845    10  2713 georgia      <NA>
## 2705 -80.92456 32.05126    10  2714 georgia      <NA>
## 2706 -80.89018 32.03980    10  2715 georgia      <NA>

Get GA county lines

counties <-map_data("county")

ga_county <- subset(counties, region == "georgia")

head(ga_county)

##            long      lat group order  region subregion
## 14370 -82.44862 31.94813   358 14727 georgia   appling
## 14371 -82.42570 31.94813   358 14728 georgia   appling
## 14372 -82.40852 31.94240   358 14729 georgia   appling
## 14373 -82.39706 31.94240   358 14730 georgia   appling
## 14374 -82.38560 31.93094   358 14731 georgia   appling
## 14375 -82.35122 31.91948   358 14732 georgia   appling

tail(ga_county)

##            long      lat group order  region subregion
## 19481 -83.66902 31.31215   516 19996 georgia     worth
## 19482 -83.74351 31.30641   516 19997 georgia     worth
## 19483 -83.74351 31.30068   516 19998 georgia     worth
## 19484 -83.78362 31.30068   516 19999 georgia     worth
## 19485 -83.78934 31.31215   516 20000 georgia     worth
## 19486 -84.01853 31.31215   516 20001 georgia     worth

Plot GA

ga_base <-ggplot(data = ga_df, mapping = aes(x = long, y = lat, group = group)) + 
  coord_fixed(1.3) + 
  geom_polygon(color = "black", fill = "gray")
ga_base #+ theme_nothing()

Plot the county boundaries in white

ga_base +
  #theme_nothing() +
  geom_polygon(data = ga_county, fill = NA, color = "white") +
  geom_polygon(color = "black", fill = NA)   # state border on top

Plot the WWSCAN sites in GA

#ga_base +
  #geom_polygon(data = ga_county, fill = NA, color = "white") +
  #geom_polygon(color = "black", fill = NA) +  # state border on top
  #geom_point(data = ga_wwscan, aes(x = long, y = lat), inherit.aes = FALSE, color = "black", size = 5) +
  #geom_point(data = ga_wwscan, aes(x = long, y = lat), inherit.aes = FALSE, color = "white", size = 4)

##Switching Gears to Key Steps in a Spatial Analysis (Mini Project)

• Identify and acquire geospatial datasets.
• Prepare spatial data objects in R.
• Visualize the data with maps.
• Perform exploratory data analysis.
• Specify and fit regression models.
• Summarize results.

Goal: Identify your research question by picking an outcome of interest and some covariates with data available at the neighborhood/county/state-level etc.

Question: In 2019, to what extent do state-level poverty rates and low educational attainment (adults aged 25+ without a high school diploma) predict state-level variations in Hepatitis A rates?

According to the WHO/CDC, Hepatitis A risk factors include those experiencing homelessness, poor sanitation, lack of safe water, living in a household with an infected person, being a sexual partner of someone with acute hepatitis A infection, use of /injecting recreational drugs, MSM, and traveling to areas of high endemicity w/o being immunized.

How could it or should wastewater data be added to a question like this? Or are there better questions framed differently Note: Questions heavily depend on data availability.

1. Finding Geospatial Datasets

Example Data Sources (read their ‘about’ pages): 1. AltasPlus from the National Center for HIV, Viral Hepatitis, STD, and Tuberculosis Prevention https://gis.cdc.gov/grasp/nchhstpatlas/maps.html 2. Data Explorer from the National Environmental Public Health Tracking Network https://ephtracking.cdc.gov/DataExplorer/

library(tigris)

## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ lubridate 1.9.5     ✔ tibble    3.3.0
## ✔ purrr     1.2.2     ✔ tidyr     1.3.2
## ✔ readr     2.2.0

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ purrr::map()    masks maps::map()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(here)

## here() starts at /Users/kymoore/Downloads

library(tmap)

#replace file paths with yours after downloading the data as a csv from the website
#poverty <- read.csv("C:/Users/kmoor38/Downloads/MapData (1).csv", skip=9) #skips the first 9 rows
#hepa <- read.csv("C:/Users/kmoor38/Downloads/MapData.csv", skip=9)
#education <- read.csv("C:/Users/kmoor38/Downloads/MapData (2).csv", skip=9)

poverty <- read.csv("/Users/kymoore/Downloads/MapData-3.csv", skip=9)
hepa <- read.csv("/Users/kymoore/Downloads/MapData-2.csv", skip=9)
education <- read.csv("/Users/kymoore/Downloads/MapData-4.csv", skip=9)

2. Data Wrangling

Construct data for analysis from downloads. Tasks:

-Identify the variables needed from the downloaded datasets.(Use FIPS for linkage/joining) -Check issues related to missing data. Helps to review documentation like data dictionaries. -Examine summary statistics of the variables and create some exploratory plots (boxplot, histogram, scatter plot).

library(dplyr)
library(stringr)

# Prepare Poverty Dataset (Select and Rename)
# Syntax: new_name = old_name
poverty_sub <- poverty %>%
  select(
    FIPS, 
    poverty_p = Percent,
    poverty_num = Numerator, 
    pov_den = Denominator
  )

# Prepare Education Dataset (Select and Rename)
education_sub <- education %>%
  select(
    FIPS, 
    education_p = Percent,
    education_num = Numerator, 
    education_den = Denominator
  )

# Merge everything into the Main Hep A Dataset
hepa_final <- hepa %>%
  left_join(poverty_sub, by = "FIPS") %>%
  left_join(education_sub, by = "FIPS")

hepa_final

##       Indicator            Geography FIPS Year      Age.Group
## 1  Hepatitis  A              Alabama    1 2019 All age groups
## 2  Hepatitis  A               Alaska    2 2019 All age groups
## 3  Hepatitis  A              Arizona    4 2019 All age groups
## 4  Hepatitis  A             Arkansas    5 2019 All age groups
## 5  Hepatitis  A           California    6 2019 All age groups
## 6  Hepatitis  A             Colorado    8 2019 All age groups
## 7  Hepatitis  A          Connecticut    9 2019 All age groups
## 8  Hepatitis  A             Delaware   10 2019 All age groups
## 9  Hepatitis  A District of Columbia   11 2019 All age groups
## 10 Hepatitis  A              Florida   12 2019 All age groups
## 11 Hepatitis  A              Georgia   13 2019 All age groups
## 12 Hepatitis  A               Hawaii   15 2019 All age groups
## 13 Hepatitis  A                Idaho   16 2019 All age groups
## 14 Hepatitis  A             Illinois   17 2019 All age groups
## 15 Hepatitis  A              Indiana   18 2019 All age groups
## 16 Hepatitis  A                 Iowa   19 2019 All age groups
## 17 Hepatitis  A               Kansas   20 2019 All age groups
## 18 Hepatitis  A             Kentucky   21 2019 All age groups
## 19 Hepatitis  A            Louisiana   22 2019 All age groups
## 20 Hepatitis  A                Maine   23 2019 All age groups
## 21 Hepatitis  A             Maryland   24 2019 All age groups
## 22 Hepatitis  A        Massachusetts   25 2019 All age groups
## 23 Hepatitis  A             Michigan   26 2019 All age groups
## 24 Hepatitis  A            Minnesota   27 2019 All age groups
## 25 Hepatitis  A          Mississippi   28 2019 All age groups
## 26 Hepatitis  A             Missouri   29 2019 All age groups
## 27 Hepatitis  A              Montana   30 2019 All age groups
## 28 Hepatitis  A             Nebraska   31 2019 All age groups
## 29 Hepatitis  A               Nevada   32 2019 All age groups
## 30 Hepatitis  A        New Hampshire   33 2019 All age groups
## 31 Hepatitis  A           New Jersey   34 2019 All age groups
## 32 Hepatitis  A           New Mexico   35 2019 All age groups
## 33 Hepatitis  A             New York   36 2019 All age groups
## 34 Hepatitis  A       North Carolina   37 2019 All age groups
## 35 Hepatitis  A         North Dakota   38 2019 All age groups
## 36 Hepatitis  A                 Ohio   39 2019 All age groups
## 37 Hepatitis  A             Oklahoma   40 2019 All age groups
## 38 Hepatitis  A               Oregon   41 2019 All age groups
## 39 Hepatitis  A         Pennsylvania   42 2019 All age groups
## 40 Hepatitis  A         Rhode Island   44 2019 All age groups
## 41 Hepatitis  A       South Carolina   45 2019 All age groups
## 42 Hepatitis  A         South Dakota   46 2019 All age groups
## 43 Hepatitis  A            Tennessee   47 2019 All age groups
## 44 Hepatitis  A                Texas   48 2019 All age groups
## 45 Hepatitis  A                 Utah   49 2019 All age groups
## 46 Hepatitis  A              Vermont   50 2019 All age groups
## 47 Hepatitis  A             Virginia   51 2019 All age groups
## 48 Hepatitis  A           Washington   53 2019 All age groups
## 49 Hepatitis  A        West Virginia   54 2019 All age groups
## 50 Hepatitis  A            Wisconsin   55 2019 All age groups
## 51 Hepatitis  A              Wyoming   56 2019 All age groups
## 52 Hepatitis  A       American Samoa   60 2019 All age groups
## 53 Hepatitis  A                 Guam   66 2019 All age groups
## 54 Hepatitis  A     Marshall Islands   68 2019 All age groups
## 55 Hepatitis  A   N. Mariana Islands   69 2019 All age groups
## 56 Hepatitis  A                Palau   70 2019 All age groups
## 57 Hepatitis  A          Puerto Rico   72 2019 All age groups
## 58 Hepatitis  A  U.S. Virgin Islands   78 2019 All age groups
##           Race.Ethnicity        Sex            Country.of.Birth
## 1  All races/ethnicities Both sexes All transmission categories
## 2  All races/ethnicities Both sexes All transmission categories
## 3  All races/ethnicities Both sexes All transmission categories
## 4  All races/ethnicities Both sexes All transmission categories
## 5  All races/ethnicities Both sexes All transmission categories
## 6  All races/ethnicities Both sexes All transmission categories
## 7  All races/ethnicities Both sexes All transmission categories
## 8  All races/ethnicities Both sexes All transmission categories
## 9  All races/ethnicities Both sexes All transmission categories
## 10 All races/ethnicities Both sexes All transmission categories
## 11 All races/ethnicities Both sexes All transmission categories
## 12 All races/ethnicities Both sexes All transmission categories
## 13 All races/ethnicities Both sexes All transmission categories
## 14 All races/ethnicities Both sexes All transmission categories
## 15 All races/ethnicities Both sexes All transmission categories
## 16 All races/ethnicities Both sexes All transmission categories
## 17 All races/ethnicities Both sexes All transmission categories
## 18 All races/ethnicities Both sexes All transmission categories
## 19 All races/ethnicities Both sexes All transmission categories
## 20 All races/ethnicities Both sexes All transmission categories
## 21 All races/ethnicities Both sexes All transmission categories
## 22 All races/ethnicities Both sexes All transmission categories
## 23 All races/ethnicities Both sexes All transmission categories
## 24 All races/ethnicities Both sexes All transmission categories
## 25 All races/ethnicities Both sexes All transmission categories
## 26 All races/ethnicities Both sexes All transmission categories
## 27 All races/ethnicities Both sexes All transmission categories
## 28 All races/ethnicities Both sexes All transmission categories
## 29 All races/ethnicities Both sexes All transmission categories
## 30 All races/ethnicities Both sexes All transmission categories
## 31 All races/ethnicities Both sexes All transmission categories
## 32 All races/ethnicities Both sexes All transmission categories
## 33 All races/ethnicities Both sexes All transmission categories
## 34 All races/ethnicities Both sexes All transmission categories
## 35 All races/ethnicities Both sexes All transmission categories
## 36 All races/ethnicities Both sexes All transmission categories
## 37 All races/ethnicities Both sexes All transmission categories
## 38 All races/ethnicities Both sexes All transmission categories
## 39 All races/ethnicities Both sexes All transmission categories
## 40 All races/ethnicities Both sexes All transmission categories
## 41 All races/ethnicities Both sexes All transmission categories
## 42 All races/ethnicities Both sexes All transmission categories
## 43 All races/ethnicities Both sexes All transmission categories
## 44 All races/ethnicities Both sexes All transmission categories
## 45 All races/ethnicities Both sexes All transmission categories
## 46 All races/ethnicities Both sexes All transmission categories
## 47 All races/ethnicities Both sexes All transmission categories
## 48 All races/ethnicities Both sexes All transmission categories
## 49 All races/ethnicities Both sexes All transmission categories
## 50 All races/ethnicities Both sexes All transmission categories
## 51 All races/ethnicities Both sexes All transmission categories
## 52 All races/ethnicities Both sexes All transmission categories
## 53 All races/ethnicities Both sexes All transmission categories
## 54 All races/ethnicities Both sexes All transmission categories
## 55 All races/ethnicities Both sexes All transmission categories
## 56 All races/ethnicities Both sexes All transmission categories
## 57 All races/ethnicities Both sexes All transmission categories
## 58 All races/ethnicities Both sexes All transmission categories
##                 Cases    Rate.per.100000 Population          poverty_p
## 1                 242                4.9  4,903,185               11.2
## 2                   2                0.3    731,545                6.6
## 3                 584                8.0  7,278,717                9.5
## 4                 203                6.7  3,017,804               11.7
## 5                 256                0.6 39,512,223                8.2
## 6                 333                5.8  5,758,736                5.8
## 7                  17                0.5  3,565,287                6.8
## 8                  36                3.7    973,764                7.6
## 9                  15                2.1    705,749                8.8
## 10              3,392               15.8 21,477,737                8.7
## 11                844                7.9 10,617,423                9.6
## 12                  1                0.1  1,415,872                6.4
## 13                 75                4.2  1,787,065                7.4
## 14                185                1.5 12,671,821                7.9
## 15              1,398               20.8  6,732,219                7.8
## 16                  9                0.3  3,155,070                7.3
## 17                 10                0.3  2,913,314                7.5
## 18              1,318               29.5  4,467,673               11.8
## 19                687               14.8  4,648,794               14.3
## 20                 45                3.3  1,344,212                6.5
## 21                 88                1.5  6,045,680                5.8
## 22                204                3.0  6,892,503                  6
## 23                 70                0.7  9,986,857                8.8
## 24                 76                1.3  5,639,632                5.2
## 25                128                4.3  2,976,149               14.8
## 26                359                5.8  6,137,428                8.8
## 27                 15                1.4  1,068,778                7.9
## 28                 15                0.8  1,934,408                6.2
## 29                102                3.3  3,080,156                8.7
## 30                309               22.7  1,359,711                4.2
## 31                610                6.9  8,882,190                6.2
## 32                104                5.0  2,096,829               13.7
## 33                391                2.0 19,453,561                9.3
## 34                154                1.5 10,488,084                9.6
## 35                  4                0.5    762,062                6.5
## 36              1,802               15.4 11,689,100                9.2
## 37                 13                0.3  3,956,971               10.8
## 38                 27                0.6  4,217,737                6.6
## 39                696                5.4 12,801,989                8.2
## 40                  6                0.6  1,059,361                6.6
## 41                662               12.9  5,148,714                9.8
## 42                  8                0.9    884,659                7.3
## 43              2,160               31.6  6,829,174                 10
## 44                159                0.5 28,995,881               10.5
## 45                 20                0.6  3,205,958                6.1
## 46                 12                1.9    623,989                5.8
## 47                309                3.6  8,535,519                6.5
## 48                181                2.4  7,614,893                6.2
## 49                467               26.1  1,792,147               11.1
## 50                 31                0.5  5,822,434                6.2
## 51                 12                2.1    578,759                6.6
## 52 Data not available Data not available       <NA> Data not available
## 53 Data not available Data not available       <NA> Data not available
## 54 Data not available Data not available       <NA> Data not available
## 55 Data not available Data not available       <NA> Data not available
## 56 Data not available Data not available       <NA> Data not available
## 57 Data not available Data not available       <NA>               39.5
## 58 Data not available Data not available       <NA> Data not available
##           poverty_num   pov_den        education_p      education_num
## 1             138,643 1,237,883               12.9            433,073
## 2              10,767   163,134                6.4             31,090
## 3             165,367 1,740,704               12.4            612,998
## 4              88,371   755,305               12.5            254,993
## 5             735,757 8,972,643                 16          4,301,513
## 6              82,467 1,421,844                7.6            302,220
## 7              60,242   885,911                9.3            231,839
## 8              18,436   242,584                9.7             66,522
## 9              10,998   124,978                8.1             41,168
## 10            442,245 5,083,272               11.6          1,794,635
## 11            245,322 2,555,440               12.1            856,025
## 12             20,237   316,206                7.6             75,978
## 13             33,224   448,979                8.5             99,906
## 14            241,666 3,059,067               10.2            882,637
## 15            127,860 1,639,230               10.4            468,083
## 16             59,043   808,802                7.4            156,481
## 17             55,012   733,489                8.2            157,757
## 18            133,255 1,129,276               12.8            390,653
## 19            159,278 1,113,831                 14            440,844
## 20             22,288   342,894                6.8             67,314
## 21             84,800 1,462,070                9.6            402,152
## 22             99,900 1,665,005                8.7            420,941
## 23            219,696 2,496,550                8.6            596,287
## 24             72,884 1,401,623                6.4            245,690
## 25            106,321   718,382               14.7            290,917
## 26            136,052 1,546,045                9.3            391,815
## 27             21,252   269,009                5.8             42,682
## 28             30,567   493,013                  8            101,272
## 29             62,967   723,755               13.1            280,915
## 30             14,549   346,413                6.7             65,672
## 31            138,976 2,241,555                9.7            599,619
## 32             68,084   496,961               14.1            201,328
## 33            431,181 4,636,356               12.4          1,694,230
## 34            252,515 2,630,365               11.4            818,924
## 35             12,318   189,509                6.5             33,019
## 36            270,717 2,942,581                9.2            740,847
## 37            105,381   975,754               11.6            303,266
## 38             67,425 1,021,588                8.6            258,399
## 39            264,714 3,228,224                  9            811,832
## 40             16,499   249,987               10.7             80,207
## 41            126,060 1,286,326               11.7            415,678
## 42             16,349   223,964                7.9             46,294
## 43            172,790 1,727,895                 12            562,817
## 44            720,052 6,857,641               15.4          2,882,388
## 45             46,419   760,968                  7            133,060
## 46              9,057   156,153                6.9             30,858
## 47            136,160 2,094,763                 10            584,712
## 48            116,740 1,882,896                8.3            440,074
## 49             52,222   470,472               12.9            164,805
## 50             91,147 1,470,116                7.2            290,862
## 51             10,089   152,859                5.5             21,250
## 52 Data not available      <NA> Data not available Data not available
## 53 Data not available      <NA> Data not available Data not available
## 54 Data not available      <NA> Data not available Data not available
## 55 Data not available      <NA> Data not available Data not available
## 56 Data not available      <NA> Data not available Data not available
## 57            305,650   773,798               21.2            488,408
## 58 Data not available      <NA> Data not available Data not available
##    education_den
## 1      3,360,058
## 2        484,058
## 3      4,944,540
## 4      2,036,456
## 5     26,937,872
## 6      3,974,943
## 7      2,496,420
## 8        687,311
## 9        505,145
## 10    15,484,502
## 11     7,080,222
## 12       996,668
## 13     1,170,997
## 14     8,694,694
## 15     4,502,015
## 16     2,123,004
## 17     1,918,081
## 18     3,048,442
## 19     3,140,201
## 20       991,152
## 21     4,183,858
## 22     4,850,576
## 23     6,894,627
## 24     3,847,212
## 25     1,979,664
## 26     4,206,162
## 27       741,950
## 28     1,271,770
## 29     2,136,466
## 30       979,750
## 31     6,191,229
## 32     1,425,988
## 33    13,664,734
## 34     7,187,077
## 35       504,375
## 36     8,049,805
## 37     2,619,244
## 38     2,988,118
## 39     9,028,036
## 40       746,952
## 41     3,563,204
## 42       588,029
## 43     4,693,962
## 44    18,772,550
## 45     1,911,592
## 46       445,558
## 47     5,872,757
## 48     5,290,324
## 49     1,281,086
## 50     4,015,285
## 51       389,847
## 52          <NA>
## 53          <NA>
## 54          <NA>
## 55          <NA>
## 56          <NA>
## 57     2,301,735
## 58          <NA>

# Remove punctuation
hepa_final_1 <- hepa_final %>% 
  # 1. Clean and convert numeric columns
  mutate(across(
    c(Cases, Rate.per.100000, poverty_p, poverty_num, education_p, education_num), 
    ~ {
      x <- na_if(.x, "Data not available")
      x <- str_remove_all(x, "[%, ]")
      as.numeric(x)
    }
  )) %>% 
  # 2. Convert fips to numeric and keep only states + DC (< 60)
  filter(as.numeric(FIPS) < 60)

# Replaces with NA
#hepa_final_1 <- hepa_final %>% 
#  mutate(across(c(Cases, Rate.per.100000, poverty_p, poverty_num, education_p, education_num), ~as.numeric(na_if(.x, "Data not available"))))


hist(hepa_final_1$Cases)

hist(hepa_final_1$Rate.per.100000)

#right skewed


summary(hepa_final_1$Cases)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0    16.0   128.0   369.5   375.0  3392.0

boxplot(hepa_final_1$Cases)

#outliers


plot(x = hepa_final_1$education_p, y = hepa_final_1$Rate.per.100000,
     main = "Low Educational Attainment vs. HepA Rate",
     xlab = "% of Population 25 years and older w/o HS diploma",
     ylab = "HepA Rate per 100,000",
     col = "black",          # Change point color
     pch = 19)              # Change point shape to solid circle


# Formula predicts Rate from education:
abline(lm(Rate.per.100000 ~ education_p, data = hepa_final_1), col = "red", lwd = 2)

plot(x = hepa_final_1$poverty_p, y = hepa_final_1$Rate.per.100000,
     main = "Poverty vs. HepA Rate",
     xlab = "% of Households living below the federal poverty level",
     ylab = "HepA Rate per 100,000",
     col = "black",          # Change point color
     pch = 19)              # Change point shape to solid circle

# Optional: Add a linear regression trend line
#abline(lm(poverty_p ~ Rate.per.100000, data = hepa_final_1), col = "red", lwd = 2)
abline(lm(Rate.per.100000 ~ poverty_p, data = hepa_final_1), col = "red", lwd = 2)

### Calculating the standardized differences
hepa_final_1 <- hepa_final_1 %>%
  mutate(
    # Standardized Hepatitis A Rate (Z-score)
    rate_z = (Rate.per.100000 - mean(Rate.per.100000, na.rm = TRUE)) / sd(Rate.per.100000, na.rm = TRUE),
    
     # Standardized Education
    edu_z = (education_p - mean(education_p, na.rm = TRUE)) / sd(education_p, na.rm = TRUE),
    
    # Standardized Poverty
    poverty_z = (poverty_p - mean(poverty_p, na.rm = TRUE)) / sd(poverty_p, na.rm = TRUE),
    
    # If you are taking the difference between two standardized proportions (e.g., Poverty vs Education):
    poverty_edu_diff_z = poverty_z - edu_z
  )

hist(hepa_final_1$rate_z)

hist(hepa_final_1$edu_z)

hist(hepa_final_1$poverty_z)

hist(hepa_final_1$poverty_edu_diff_z)

# 1. Set up a grid of 1 row and 3 columns
par(mfrow = c(2, 2))

# 2. Draw each histogram
hist(hepa_final_1$rate_z, 
     main = "Standardized Hepatitis A Rate", 
     xlab = "Z-score", 
     col = "lightblue")

hist(hepa_final_1$poverty_z, 
     main = "Standardized Poverty", 
     xlab = "Z-score", 
     col = "salmon")

hist(hepa_final_1$edu_z, 
     main = "Standardized Education", 
     xlab = "Z-score", 
     col = "lightgreen")

hist(hepa_final_1$poverty_edu_diff_z, 
     main = "Standardized Covariate Differences", 
     xlab = "Z-score", 
     col = "purple")

# 3. Reset plotting grid back to default (1x1)
par(mfrow = c(1, 1))

3. Creating Geospatial Datasets and Mapping

-Download shapefiles for your study region -Create maps showing your outcome and exposures of interests

library (tigris) #downloads shapefiles from the US Census
options(tigris_use_cache = TRUE)

# Load spatial geometry for US states and shift AK/HI
us_shifted <- states(cb = TRUE, resolution = "20m", year = 2020) %>% 
  filter(as.numeric(STATEFP) < 60) %>%  # Keep only states + DC
  shift_geometry()

# Plot using ggplot2
ggplot(us_shifted) +
  geom_sf(fill = "grey80", color = "white", size = 0.2) +
  theme_void() + # Removes empty margins completely 
  theme(legend.position = "none") # Drops the legend box

# Standardize FIPS to a 2-digit character vector with leading zeros
hepa_final_1 <- hepa_final_1 %>%
  mutate(FIPS = str_pad(as.character(FIPS), width = 2, side = "left", pad = "0"))

# Join your clean spatial geometry with your non-spatial dataset
# Note: Ensure state FIPS codes match 
hepa_sf <- us_shifted %>% 
  left_join(hepa_final_1, by = c("STATEFP" = "FIPS"))

Now to map the outcome, using the mapping package tmap’s routines tm_shape and tm_fill. tmap builds the maps as objects making it easier to customize parts of the map (legends, credits, and the layout) and to build multi-map figures.

hepa_map <- tm_shape(hepa_sf) + 
  tm_fill('Rate.per.100000', 
style='quantile', 
palette='BuPu', 
title='HepA Rates \n Per 100,000 people') + 
  tm_borders(alpha=0.7) + 
  tm_credits('Quantile (Equal-Frequency) Class Intervals', 
             position=c('RIGHT', 'BOTTOM')) + 
  tm_layout(main.title="HepA Rates in the US",
            inner.margins = c(0.1, 0.1, 0.05, 0.05), 
            main.title.size=1.2, legend.title.size=0.5,  
            legend.text.size=0.5)

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_fill()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_borders()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').
## [v3->v4] `tm_borders()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

hepa_map

## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "BuPu" is named
## "brewer.bu_pu"
## Multiple palettes called "bu_pu" found: "brewer.bu_pu", "matplotlib.bu_pu". The first one, "brewer.bu_pu", is returned.

---

Next, map standardized log poverty.

poverty_map <- tm_shape(hepa_sf) + 
  tm_fill('poverty_z', 
style='quantile', 
palette='BuPu', 
title='Standardized Log \n Poverty (%)') + # "\n" moves text to the next line 
  tm_borders(alpha=0.2) + 
  tm_credits('Quantile (Equal-Frequency) Class Intervals', 
             position=c('RIGHT', 'BOTTOM')) + 
  tm_layout(main.title="% of Households living below the federal poverty level in the US",
            inner.margins = c(0.1, 0.1, 0.05, 0.05), 
            main.title.size=1.2, legend.title.size=0.5,  
            legend.text.size=0.5)

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_fill()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_borders()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

poverty_map

## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "BuPu" is named
## "brewer.bu_pu"
## Multiple palettes called "bu_pu" found: "brewer.bu_pu", "matplotlib.bu_pu". The first one, "brewer.bu_pu", is returned.
## 
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

---

Plot standardized log educational attainment.

edu_map <- tm_shape(hepa_sf) + 
  tm_fill('edu_z', 
style='quantile', 
palette='BuPu', 
title='Standardized Log \n Low Educational Attainment') + #reverse coding?? 
  tm_borders(alpha=0.2) + 
  tm_credits('Quantile (Equal-Frequency) Class Intervals', 
             position=c('RIGHT', 'BOTTOM')) + 
  tm_layout(main.title="% of Population 25 years and older w/o HS diploma",
            inner.margins = c(0.1, 0.1, 0.05, 0.05), 
            main.title.size=1.2, legend.title.size=0.5,  
            legend.text.size=0.5)

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_fill()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_borders()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

edu_map

## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "BuPu" is named
## "brewer.bu_pu"
## Multiple palettes called "bu_pu" found: "brewer.bu_pu", "matplotlib.bu_pu". The first one, "brewer.bu_pu", is returned.

4. Spatial Exploratory Analysis

What should we use as the neighborhood structure to induce spatial correlation? Tasks: -Construct an adjacency matrix for the analysis. -Is there evidence of spatial auto-correlation in the outcome and covariates?

Defining Spatial Proximity

#install.packages('spdep')
library (spdep)

## Loading required package: spData

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

## Loading required package: sf

## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

#Drop missing data
hepa_sf = subset (hepa_sf, !is.na (hepa_sf$Rate.per.100000))

#Because our shapefile here is un-projected, distance operations can be difficult. Alternatively, we can project the data to, for example, the Lambert Conformal Conic North America.
# 1. Project the data to Lambert Conformal Conic North America (ESRI:102004)
# (Recommended for accurate distance and spatial operations across the US)
hepa_sf_proj = st_transform(hepa_sf, crs = "ESRI:102004")
center_proj = st_centroid(hepa_sf_proj)
#Create 1-st order adjacency and distance buffer (100km)
#Plot the projected geometry
#plot (st_geometry(hepa_sf_proj))
nb.1st = poly2nb (hepa_sf_proj)
nb.buffer = dnearneigh(center_proj, d1=0, d2 = 25*1000)
#nb2 = poly2nb (hepa_sf_proj, queen=FALSE) 
#take the union of neighbors
nb = union.nb (nb.1st, nb.buffer) #take the union of neighbors


#A lot of good summaries about the neighbor structure
summary (nb)

## Neighbour list object:
## Number of regions: 51 
## Number of nonzero links: 218 
## Percentage nonzero weights: 8.381392 
## Average number of links: 4.27451 
## 2 regions with no links:
## 50, 51
## 3 disjoint connected subgraphs
## Link number distribution:
## 
##  0  1  2  3  4  5  6  7  8 
##  2  1  5  9 10 10 10  2  2 
## 1 least connected region:
## 37 with 1 link
## 2 most connected regions:
## 6 18 with 8 links

The function nb2mat () then creates a proximity/weight matrix

#Create weight matrix
W= nb2listw(nb, style="W", zero.policy = TRUE)
#W = nb2mat (nb, style = "B", zero.policy = TRUE) #Stopped code - Hawaii (or states separated by water/projected geometry boundaries) do not share a physical border with any other state SO add zero.policy

W

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 51 
## Number of nonzero links: 218 
## Percentage nonzero weights: 8.381392 
## Average number of links: 4.27451 
## 2 regions with no links:
## 50, 51
## 3 disjoint connected subgraphs
## 
## Weights style: W 
## Weights constants summary:
##    n   nn S0       S1       S2
## W 49 2401 49 24.57728 202.0046

#center_proj = st_centroid(hepa_sf_proj) ##Extract centroids of polygons
#center = st_transform(center_proj, crs = 4326) #Re-project back to WGS 84
#center = st_coordinates (center)                 

#plot (st_geometry(hepa_sf))
#plot (nb, center, add = TRUE, col = "blue")
#plot (nb2, center, add = TRUE, col = "red", lwd = 1)
#title (main="Blue lines = additional neighbours by queen's case")

## Using kth-nearest distance
#nb.k1 = knn2nb (knearneigh(center_proj,k=1))
#nb.k2 = knn2nb (knearneigh(center_proj,k=2), row.names=row.names(center))

#plot (st_geometry(hepa_sf))
#plot (nb.k2, center, add = TRUE, col = "red", lwd = 1)
#plot (nb.k1, center, add = TRUE, col = "blue", lwd = 1)
#title (main="Blue lines = additional second nearest-neighbour")

## Using buffer distance
#nb.buffer1 = dnearneigh(center_proj, d1=0, d2 = 25*1000)
#nb.buffer2 = dnearneigh(center_proj, d1=0, d2 = 40*1000)

#plot (st_geometry(hepa_sf))
#plot (nb.buffer2, center, add = TRUE, col = "red", lwd = 1)
#plot (nb.buffer1, center, add = TRUE, col = "blue", lwd =2)
#legend ("topright", legend =c("<25 km", "< 40 km"), col= c("blue", "red"), pch = 16)

Monte Carlo-based Moran’s I

moran.mc(hepa_sf$Rate.per.100000, W, nsim = 10000)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  hepa_sf$Rate.per.100000 
## weights: W  
## number of simulations + 1: 10001 
## 
## statistic = 0.19811, observed rank = 9814, p-value = 0.0187
## alternative hypothesis: greater

moran.mc(hepa_sf$poverty_z, W, nsim = 10000) #look at the p value

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  hepa_sf$poverty_z 
## weights: W  
## number of simulations + 1: 10001 
## 
## statistic = 0.50058, observed rank = 10001, p-value = 9.999e-05
## alternative hypothesis: greater

moran.mc(hepa_sf$edu_z, W, nsim = 10000) #look at the p value

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  hepa_sf$edu_z 
## weights: W  
## number of simulations + 1: 10001 
## 
## statistic = 0.52672, observed rank = 10001, p-value = 9.999e-05
## alternative hypothesis: greater

What Moran’s I values mean:

I = 1 → perfect clustering (similar values always neighbor each other) I = 0 → completely random spatial pattern I = -1 → perfect dispersion (opposites always neighbor each other)

The big takeaway: Poverty and education are strongly spatially clustered, which justifies using a spatial model (like CAR/INLA) rather than a standard Poisson regression — because nearby states aren’t independent of each other.

5. Spatial Regression Modeling

Goal: Use conditional autoregressive models to estimate the association between the outcome and covariates.

-Describe or write down the statistical model. -Fit the model using INLA and summarize results. -Compare model to a coventional model without accounting for spatial dependence.

#For state 𝑠, let 𝑌𝑠 be the number of Hepatitis A cases, 𝑋𝑠 be poverty or low educational attainment, and 𝑃𝑠 be the population size. We consider several Poisson log-linear models of the following form.

#𝑌𝑠log(𝜇𝑠)∼Poisson(𝜇𝑠)=log(𝑃𝑠)+𝛽0+𝛽1𝑋𝑠+𝜃𝑠
#The standard Poisson model assumes no spatial variation (𝜃𝑠=0 for all 𝑠).

#The independent random intercept model assumes that 𝜃𝑠∼𝑁(0,𝜏2).

#The spatial random intercept model (under the improper CAR model formulation) assumes that 𝜃𝑠=𝜂𝑠+𝛾𝑠, where 𝜂𝑠∼𝑁(0,𝜏2) and 𝛾𝑠 is a conditional autoregressive model.

library (INLA)

## Warning: package 'INLA' was built under R version 4.5.3

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

## 

# 1. Strip commas and convert population to numeric
hepa_sf <- hepa_sf %>%
  mutate(
    Population = as.numeric(str_remove_all(Population, ",")),
    log_pop = log(Population)
  )

#Some data wrangling first
hepa_sf$areal_ID = 1:nrow (hepa_sf) #Create a county ID from 1 to 226
hepa_sf$log_pop = log (hepa_sf$Population)
hepa_sf$poverty_std = hepa_sf$poverty_p / 0.1 #standardized by the interquartile range (0.1)

#INLA needs a specific summary of adjacency structure
nb2INLA ("adj.txt", nb) #write out a file
G <- inla.read.graph(filename = "adj.txt") #read in to get INLA's graph format

#Fit glm
fit.glm = glm (Cases~poverty_std+offset (log_pop), dat = hepa_sf, family = "poisson")

#Fit exchangeable (random effect model)
#Note: "E" here is the offset and it does not need to be log-transformed
fit.exch = inla (Cases~1+f(areal_ID)+poverty_std, E= hepa_sf$Population, family = "poisson", data = hepa_sf,
                 control.compute = list(dic = TRUE, waic = TRUE, return.marginals.predictor=TRUE))

#Fit improper CAR model (Besag)
fit.iCAR = inla (Cases~1+f(areal_ID,  model = "besag", graph = G)+poverty_std, E = hepa_sf$Population, family = "poisson", data = hepa_sf,control.compute = list(dic = TRUE, waic = TRUE))


library(knitr)

# 1. Extract GLM results
glm_sum <- summary(fit.glm)$coefficients["poverty_std", ]
glm_beta <- glm_sum["Estimate"]
glm_se   <- glm_sum["Std. Error"]
glm_rr   <- exp(glm_beta)
glm_lower <- exp(glm_beta - 1.96 * glm_se)
glm_upper <- exp(glm_beta + 1.96 * glm_se)

# 2. Extract Independent Random Effect (Exchangeable) results
exch_beta  <- fit.exch$summary.fixed["poverty_std", "mean"]
exch_se    <- fit.exch$summary.fixed["poverty_std", "sd"]
exch_rr    <- exp(exch_beta)
exch_lower <- exp(fit.exch$summary.fixed["poverty_std", "0.025quant"])
exch_upper <- exp(fit.exch$summary.fixed["poverty_std", "0.975quant"])

# 3. Extract Spatial (iCAR) results
icar_beta  <- fit.iCAR$summary.fixed["poverty_std", "mean"]
icar_se    <- fit.iCAR$summary.fixed["poverty_std", "sd"]
icar_rr    <- exp(icar_beta)
icar_lower <- exp(fit.iCAR$summary.fixed["poverty_std", "0.025quant"])
icar_upper <- exp(fit.iCAR$summary.fixed["poverty_std", "0.975quant"])

# 4. Combine into a clean Data Frame
results_table <- data.frame(
  Model = c("GLM", "Ind Random Effect", "Spatial"),
  beta = round(c(glm_beta, exch_beta, icar_beta), 2),
  `Std Err` = round(c(glm_se, exch_se, icar_se), 2),
  RR = round(c(glm_rr, exch_rr, icar_rr), 2),
  `95% Interval` = sprintf("(%.2f, %.2f)", 
                           c(glm_lower, exch_lower, icar_lower), 
                           c(glm_upper, exch_upper, icar_upper)),
  check.names = FALSE
)

# Display formatted table in R Console
kable(results_table, align = "lcccc")

Model	beta	Std Err	RR	95% Interval
GLM	0.02	0.00	1.02	(1.02, 1.02)
Ind Random Effect	0.02	0.01	1.02	(1.01, 1.04)
Spatial	0.01	0.01	1.01	(1.00, 1.03)