Analysis.Rmd

---
title: "Analysis"
output: html_document
date: '2022-07-06'
---

# Setup 

```{r, results='hide'}
# Installs any necessary packages if not already installed 
for(
  pkg in c(
    "tidyverse", "ipumsr", "sf", "terra", "ggspatial", "gtools",
    "srvyr", "survey", "lme4", "broom.mixed", "broom", "remotes"
  )
){
  if(!require(pkg, quietly = TRUE, character.only = TRUE)){
    install.packages(pkg)
  }
}

# Installs a development version that plots raster data from `terra`
# See https://github.com/paleolimbot/ggspatial/
if(packageVersion("ggspatial") != '1.1.5.9000'){
  remotes::install_local("cellar/ggspatial_1.1.5.9000.tar.gz")
}
```


```{r, results='hide'}
# IPUMS & Spatial Data Tools
library(tidyverse)
library(ipumsr) 
library(sf)
library(terra)
library(ggspatial)
```

```{r, results='hide'}
# Analysis Tools 
library(srvyr)
library(survey)
library(lme4)
library(broom.mixed)
library(broom)
```


# 1 - IPUMS PMA Data 

## Load IPUMS PMA data into R

```{r, results='hide'}
pma <- read_ipums_micro(
  ddi = "data/pma_00136.xml", # Use your own file name here 
  data = "data/pma_00136.dat.gz" # Use your own file name here 
)
```

```{r}
pma
```

## Data format 

```{r}
pma %>% count(RESULTFQ_1, RESULTFQ_2)

pma <- pma %>% filter(RESULTFQ_1 == 1 & RESULTFQ_2 == 1)
```


# 2- PMA GPS data 

## Load PMA GPS data into R 

```{r}
gps <- read_csv("data/pma_gps.csv")

gps
```

## GPS data as a Simple Features Object

```{r}
gps <- gps %>% 
  st_as_sf(
    coords = c("GPSLONG", "GPSLAT"), 
    crs = 4326
  ) 

gps
```

## Adding a Shapefile

```{r}
shape <- st_read("data/geobf") %>% select(ADMIN_NAME)

shape
```

## Mapping GPS coordinates with a Shapefile 

```{r}
ggplot() + 
  layer_spatial(gps) + 
  theme_minimal() 

ggplot() + 
  layer_spatial(gps) + 
  layer_spatial(shape, alpha = 0) + 
  theme_minimal() 
```

## Projection to Meters 

```{r}
gps <- gps %>% st_transform(crs = 32630)
shape <- shape %>% st_transform(crs = 32630)
```

## Creating Buffers 

```{r, results='hide'}
gps_buffer <- gps %>% st_buffer(5000)
```

## Intersecting Regional Boundaries

```{r}
test_ea <- 854131009 

ggplot() + 
  layer_spatial(gps_buffer %>% filter(EAID == test_ea), alpha = 0) + 
  layer_spatial(gps %>% filter(EAID == test_ea), color = "red") + 
  annotation_spatial(shape, alpha = 0) + 
  theme_minimal()

gps_buffer <- gps_buffer %>% 
  st_intersection(shape) %>% 
  st_filter(gps)

ggplot() + 
  layer_spatial(gps_buffer %>% filter(EAID == test_ea), alpha = 0) + 
  layer_spatial(gps %>% filter(EAID == test_ea), color = "red") + 
  annotation_spatial(shape, alpha = 0) + 
  theme_minimal()
```

## Back to arc-degrees

```{r}
gps_buffer <- gps_buffer %>% st_transform(crs = 4326)
shape <- shape %>% st_transform(crs = 4326)
ggplot() + 
  layer_spatial(shape, alpha = 0) + 
  layer_spatial(gps_buffer, alpha = 0) + 
  theme_minimal()
```


# 3 - Mapping Birth Outcomes

## Birth outcomes for individuals

```{r}
pma %>% count(PANELBIRTH_2)

pma %>% count(PANELBIRTH_2 == 1)
```

## Birth outcomes by EA 

```{r}
ea_summary <- pma %>% 
  as_survey_design(weight = PANELWEIGHT) %>% 
  group_by(ea = EAID_1, urban = URBAN == 1) %>% 
  summarise(birth_prop = survey_mean(PANELBIRTH_2 == 1, vartype = NULL)) %>% 
  ungroup() %>% 
  mutate(
    ntile = ntile(birth_prop, 2),
    many_births = ntile == 2
  )

ea_summary
```

## Merging GPS data 

```{r}
ea_summary_gps <- full_join(
  gps_buffer %>% select(ea = EAID), 
  ea_summary, 
  by = "ea"
)

ea_summary_gps
```

## Mapping births by EA

```{r, fig.align='center'}
ggplot() + 
  layer_spatial(ea_summary_gps, aes(fill = many_births)) + 
  layer_spatial(shape, alpha = 0) + 
  theme_minimal()
```


# 4 - CHIRPS: Annual rainfall summary 

## Load Raster Data into R

```{r}
june5_2020 <- rast("data/chirps/20200605.tif")

june5_2020
```

## June 5, 2020

```{r, `code-line-numbers`="5-7", fig.align='center'}
ggplot() + 
  layer_spatial(mask(june5_2020, vect(shape), touches = FALSE)) + 
  layer_spatial(shape, alpha = 0) + 
  theme_minimal() + 
  scale_fill_gradient2(
    low = "transparent", 
    high = "#4375B7", 
    na.value = "transparent"
  ) + 
  labs(fill = "Rainfall (mm)")
```

## Aggregation 

```{r, `code-line-numbers`="3-10"}
years <- map(2001:2020, ~{
  list.files("data/chirps/", pattern = paste0("^", .x), full.names = TRUE)
})

years <- set_names(years, 2001:2020) 
```

## Raster Layers 

```{r}
years <- years %>% map(~.x %>% rast)

years$`2020`
```

## Seasonal Rainfall Accumulation

```{r}
years$`2020` %>% sum()

chirps_seasonal_sum <- map(years, ~.x %>% sum)
```

## One Layer per Year 

```{r}
chirps_seasonal_sum <- rast(chirps_seasonal_sum)

chirps_seasonal_sum
```

## How rainy was 2020? 

```{r}
chrips_avg <- mean(chirps_seasonal_sum)
chirps_sd <- stdev(chirps_seasonal_sum)
chirps_z <- (chirps_seasonal_sum - chrips_avg) / chirps_sd

chirps_z
```

## 2020 Seasonal Rainfall Accumulation 

```{r}
ggplot() + 
  layer_spatial(mask(chirps_z$`2020`, vect(shape), touches = FALSE)) + 
  layer_spatial(shape, alpha = 0) +
  theme_minimal() + 
  scale_fill_gradient2(
    low = "#D72B22", 
    high = "#4375B7", 
    na.value = "transparent"
  ) + 
  labs(fill = "Z-score (2001-2020)")
```


# 5 - Bringing it All Together 

## Spatially Weighted Averages

```{r}
test_buffer <- ea_summary_gps %>% 
  filter(ea == test_ea)

ggplot() + 
  layer_spatial(
    chirps_z$`2020` %>% 
      crop(test_buffer, snap = "out")
  ) + 
  layer_spatial(test_buffer, alpha = 0) + 
  theme_minimal() + 
  scale_fill_gradient2(
    high = "#4375B7", 
    limits = c(0, 2)
  )

test_z <- chirps_z$`2020` %>%
  extract(vect(test_buffer), weights = TRUE)

test_z %>% summarise(wtd_mean = weighted.mean(`2020`, weight))
```

## Enumeration Area Z-scores

```{r}
ea_summary_gps <- ea_summary_gps %>% rowid_to_column("ID")

ea_summary_z <- chirps_z$`2020` %>% 
  extract(vect(ea_summary_gps), weights = TRUE) %>% 
  tibble() %>% 
  group_by(ID) %>% 
  summarise(z = weighted.mean(`2020`, weight)) 

ea_summary_z

ea_summary_z <- ea_summary_gps %>% 
  st_drop_geometry() %>% 
  select(ID, ea) %>% 
  full_join(ea_summary_z, by = "ID")

ea_summary_z
```


# 6 - Analysis

## Join Z-scores with Original PMA Data 

```{r}
pma_chirps <- pma %>% 
  transmute(
    ea = EAID_1, 
    wt = PANELWEIGHT,
    birth = PANELBIRTH_2 == 1,
    urban = URBAN == 1
  ) %>% 
  left_join(ea_summary_z, by = "ea") 

pma_chirps
```

## Binomial model with cluster-robust SEs

```{r}
# Weighted with cluster-robust SEs
pma_chirps %>%
  as_survey_design(weight = wt, id = ea) %>%
  svyglm(birth ~ z*urban, design = ., family = "quasibinomial") %>%
  tidy(exp = TRUE, conf.int = TRUE) %>%
  mutate(sig = gtools::stars.pval(p.value))
```

## Multilevel model with EA random effects 

```{r}
glmer(
  birth ~ z*urban + (1 | ea), 
  data = pma_chirps,
  family = "binomial") %>%
  tidy(exp = TRUE, conf.int = TRUE) %>%
  mutate(sig = gtools::stars.pval(p.value))
```