fp_essentials.qmd

---
title: "Final Report Essentials"
author: "Jeremy Springman"

institute: "University of Pennsylvania"
date: last-modified
toc: true

bibliography: references.bib

format: 
  html:
    self-contained: true

editor: source
---



```{r }
#| include: false
#| warning: false
#| message: false

library(ggplot2)
library(readr)
library(ggdag)
library(tidyverse)
library(gt)
library(modelsummary)

# read-in data
# cdat = read_csv("https://raw.githubusercontent.com/jrspringman/psci3200-globaldev/main/workshops/cambodia_ngos/data.csv")
# dat = read_csv("https://raw.githubusercontent.com/jrspringman/psci3200-globaldev/main/workshops/aau_survey/clean_endline_did.csv" ) %>%

#dat = read_csv("https://raw.githubusercontent.com/jrspringman/psci3200-globaldev/main/workshops/aau_survey/clean_endline_did.csv" ) %>%
dat = read_csv(here::here("workshops/aau_survey/clean_endline_did.csv" )) %>%
    # clean home region variable
  mutate(q8_baseline = ifelse(q8_baseline == "Southern Nations, Nationalities, and Peoples Region", "SNNPR", q8_baseline), 
         q8_baseline = str_remove(q8_baseline, " Region"))
# create color palette for plotting
palette = MetBrewer::met.brewer(name = "Cross")

```

# Index Variables

## Additive Scale

```{r}
#| echo: false
#| warning: false

# Contacted gov't official
# Signed a petition

## Find participation measures that are based on a count
add_ecols = grep("q13_4_1$|q13_5_1$", names(dat), value = T)

dat = dat %>%
  mutate(add_participation_end =  rowSums(across(add_ecols) ) )
  #mutate(add_participation_end =  rowSums(across(add_ecols), na.rm = T) )

#dat$add_participation_end = dat$q13_4_1 + dat$q13_5_1

```


```{r}
#| echo: false
#| warning: false

# participation
participate = dat %>% 
  select(`Additive Index` = add_participation_end, 
         `Signed petition` = q13_5_1, 
         `Contact official` = q13_4_1) %>% 
  drop_na() %>% 
  pivot_longer(everything()) %>% 
  group_by(name, value) %>% 
  tally() %>% 
  mutate(pct = n/sum(n))



ggplot(participate, aes(x = value, y = pct)) + 
  geom_col(fill = palette[9]) + 
  facet_wrap(vars(name), scales = "free") + 
  scale_y_continuous(labels = scales::percent) +
  labs(y = "Percent of respondents", x = NULL, 
       title = "How do youth participate?", 
       subtitle = glue::glue("Number of respondents = {scales::comma(sum(participate$n)/7)}")) 

```



## Averaged Z-Scores


```{r}
#| echo: false
#| warning: false


## Find participation measures that are based on likert
# baseline
bcols = grep("^q13_.*_baseline$", names(dat), value = T)
dat[, paste0(bcols, "_st")] = dat[, bcols]
bcols = paste0(bcols,"_st")

# endline
ecols = grep("^q13_[1-7]_\\d$", names(dat), value = T)
dat[, paste0(ecols, "_st")] = dat[, ecols]
ecols = paste0(ecols,"_st")


# Create treatment variable
dat = dat %>% mutate(moved = case_when(q8_baseline == "Addis Ababa" ~ 0, TRUE ~ 1) )

# clean q13_
levels = c("Never", "Once or Twice", "More than twice", "More than 5 times", 
           "More than 10 times")
dat = dat %>% 
  mutate(across(c(bcols), 
                .fns = ~ factor(.x, levels = levels)))

# Create z-score function from Kling, Liberman, and Katz (2007)
z_score = function(x, y){
  # calculate mean and sd of control group
  c_mean = mean( as.numeric( unlist(x[, y])) , na.rm = T)
  c_sd = sd( as.numeric( unlist(x[, y])) , na.rm = T)
  # subtract control group mean; divide by control group SD
  ( as.numeric(x[, y, drop = TRUE]) - c_mean) / c_sd
}

# calculate z-scores
for (i in c(bcols, ecols)) {
  dat[,i] = z_score(dat, i)
}

dat = dat %>% 
  rowwise() %>% 
  mutate( z_participation_end = mean(c_across(all_of(bcols)), na.rm = TRUE)) %>% 
  mutate( z_participation_base = mean(c_across(all_of(ecols)), na.rm = TRUE)) %>%
  ungroup()

```


```{r}
#| echo: false
#| warning: false


# raw
participate = dat %>% 
  select(`Attend meeting` = q13_1_baseline, 
         `Protest` = q13_2_baseline, 
         `NGO event` = q13_3_baseline, 
         `Signed petition` = q13_4_baseline, 
         `Contact official` = q13_5_baseline, 
         `Contact student representative` = q13_6_baseline, 
         `Contact an NGO` = q13_7_baseline) %>% 
  drop_na() %>% 
  pivot_longer(everything()) %>% 
  group_by(name, value) %>% 
  tally() %>% 
  mutate(pct = n/sum(n)) %>% 
  mutate(value = fct_relevel(value, "Never", after = Inf), 
         value = fct_rev(value))

ggplot(participate , aes(x = value, y = pct)) + 
  geom_col(fill = palette[9]) + 
  facet_wrap(vars(name), scales = "free", ncol = 2) + 
  scale_y_continuous(labels = scales::percent) +
  labs(y = "Percent of respondents", x = NULL, 
       title = "How do youth participate?", 
       subtitle = glue::glue("Number of respondents = {scales::comma(sum(participate$n)/7)}")) + 
  scale_x_discrete(labels = scales::label_wrap(10))

```


```{r}
#| echo: false
#| warning: false

# standardized
participate = dat %>% 
  select(`Z-Score Index` = z_participation_base,
         `Attend meeting` = q13_1_baseline_st, 
         `Protest` = q13_2_baseline_st, 
         `NGO event` = q13_3_baseline_st, 
         `Signed petition` = q13_4_baseline_st, 
         `Contact official` = q13_5_baseline_st, 
         `Contact student representative` = q13_6_baseline_st, 
         `Contact an NGO` = q13_7_baseline_st) %>% 
  drop_na() %>% 
  pivot_longer(everything()) 


ggplot(participate , aes(x = value )) + 
  geom_histogram(aes(y = after_stat(count / sum(count))), 
                 binwidth=.5, fill = palette[9]) +
  #geom_col(fill = palette[9]) + 
  facet_wrap(vars(name),  ncol = 2) + 
  scale_y_continuous(labels = scales::percent) +
  labs(y = "Percent of respondents", x = NULL, 
       title = "How do youth participate?")


```




# Interaction Terms

## Differences across groups

### Binary

```{r}
#| echo: false
#| warning: false

table(dat$class_year)

dat = dat %>% mutate(year = as.numeric(case_when(class_year == "Year I" ~ 0,
                                                 class_year == "Year II" ~ NA,
                                                 class_year == "Year III" ~ 1 )) )
table(dat$year)

models <- list()
models[['Bivariate']] = lm(z_participation_base ~ moved, dat)
models[['Multivariate']] = lm(z_participation_base ~ moved + year, dat)
models[['Interaction']] = lm(z_participation_base ~ moved + year + moved*year, dat)

modelsummary(
  models,
  estimate  = "{estimate}{stars} ({std.error})",
             statistic = NULL,
  gof_omit = 'IC|RMSE|Log|F|R2$|Std.')

```


```{r}
#| echo: false
#| warning: false
#| 
# Predict values based on the model
predictions <- predict(lm(z_participation_base ~ moved * year, dat), newdata = dat)

# Create a dataframe with predictions
prediction_df <- data.frame(year = dat$year, moved = dat$moved, prediction = predictions)

# Plot the interaction
ggplot(prediction_df, aes(x = year, y = prediction, color = factor(moved))) +
  geom_line() +
  geom_point() +
  scale_x_continuous(breaks = c(0, 1), labels=c("Y1", "Y3")) +
  labs(x = "Year", y = "Predicted z_participation_base", color = "Moved") +
  ggtitle("Interaction Plot of z_participation_base by Year and Moved")

```

### Continuous-ish

```{r}
#| echo: false
#| warning: false

dat = dat %>% mutate(year = as.numeric(case_when(class_year == "Year I" ~ 0,
                                                 class_year == "Year II" ~ 1,
                                                 class_year == "Year III" ~ 2 )) )
table(dat$year)

models <- list()
models[['Bivariate']] = lm(z_participation_base ~ moved, dat)
models[['Multivariate']] = lm(z_participation_base ~ moved + year, dat)
models[['Interaction']] = lm(z_participation_base ~ moved + year + moved*year, dat)

modelsummary(
  models,
  estimate  = "{estimate}{stars} ({std.error})",
             statistic = NULL,
  gof_omit = 'IC|RMSE|Log|F|R2$|Std.')

```


```{r}
#| echo: false
#| warning: false
# Predict values based on the model
predictions <- predict(lm(z_participation_base ~ moved * year, dat), newdata = dat)

# Create a dataframe with predictions
prediction_df <- data.frame(year = dat$year, moved = dat$moved, prediction = predictions)

# Plot the interaction
ggplot(prediction_df, aes(x = year, y = prediction, color = factor(moved))) +
  geom_line() +
  geom_point() +
  scale_x_continuous(breaks = c(0, 1, 2), labels=c("Y1", "Y2", "Y3")) +
  labs(x = "Year", y = "Predicted z_participation_base", color = "Moved") +
  ggtitle("Interaction Plot of z_participation_base by Year and Moved")

```

### More complex


```{r}
#| echo: false
#| warning: false

dat = dat %>% mutate(year1 = as.numeric(case_when(class_year == "Year I" ~ 1,
                                                 class_year == "Year II" ~ 0,
                                                 class_year == "Year III" ~ 0 )),
                     year2 = as.numeric(case_when(class_year == "Year I" ~ 0,
                                                 class_year == "Year II" ~ 1,
                                                 class_year == "Year III" ~ 0 )),
                     year3 = as.numeric(case_when(class_year == "Year I" ~ 0,
                                                 class_year == "Year II" ~ 0,
                                                 class_year == "Year III" ~ 1 )),)
table(dat$year)

models <- list()
models[['Bivariate']] = lm(z_participation_base ~ moved, dat)
models[['Multivariate']] = lm(z_participation_base ~ moved + year2 + year3, dat)
models[['Interaction']] = lm(z_participation_base ~ moved + year2 + year2*moved + year3 + year3*moved , dat)

modelsummary(
  models,
  estimate  = "{estimate}{stars} ({std.error})",
             statistic = NULL,
  gof_omit = 'IC|RMSE|Log|F|R2$|Std.')

```


```{r}
#| echo: false
#| warning: false

# Assuming dat already has year2 and year3 as binary indicators
# If not, create them based on your criteria for what constitutes year2 and year3

# Model including interactions
model <- lm(z_participation_base ~ moved + year2 + year2*moved + year3 + year3*moved, dat)

# Generate new data for predictions
# This data should cover all combinations of moved, year2, and year3
newdata <- expand.grid(moved = c(0, 1), year2 = c(0, 1), year3 = c(0, 1))

# Exclude impossible combinations (e.g., year2 and year3 both being 1)
newdata <- newdata[!with(newdata, year2 == 1 & year3 == 1),]

# Predict values
predictions <- predict(model, newdata = newdata)

# Add predictions to the newdata dataframe
newdata$prediction <- predictions

# Convert year2 and year3 back to a single year variable for plotting
newdata$year <- with(newdata, ifelse(year2 == 1, "1", ifelse(year3 == 1, "2", "0")))

# Plot
library(ggplot2)
ggplot(newdata, aes(x = as.numeric(year), y = prediction, color = factor(moved))) +
  geom_point() +
  geom_line() +
  scale_x_continuous(breaks = c(0, 1, 2), labels=c("Y1", "Y2", "Y3")) +
  labs(x = "Year", y = "Predicted z_participation_base", color = "Moved") +
  ggtitle("Interaction Plot of z_participation_base by Year and Moved")


```


## Differences over time

### Cross-Sectional

```{r}
#| echo: false
#| warning: false

regd = dat %>% select(z_participation_end, z_participation_base, moved, response_id ) %>%
  pivot_longer(cols = c(z_participation_end, z_participation_base),
               names_to = "time",
               values_to = "z_participation") %>%
  mutate(time = case_when(time == "z_participation_end" ~ 1, 
                          TRUE ~ 0))

models <- list()
models[['Bivariate']] = lm(z_participation ~ moved, regd)
models[['Multivariate']] = lm(z_participation ~ moved + time, regd)
models[['Interaction']] = lm(z_participation ~ moved + time + moved*time, regd)

modelsummary(
  models,
  estimate  = "{estimate}{stars} ({std.error})",
             statistic = NULL,
  gof_omit = 'IC|RMSE|Log|F|R2$|Std.')

```



```{r}
#| echo: false
#| warning: false
# Predict values based on the model
predictions <- predict(lm(z_participation ~ moved*time, regd), newdata = regd)

# Create a dataframe with predictions
prediction_df <- data.frame(time = regd$time, moved = regd$moved, prediction = predictions)

# Plot the interaction
ggplot(prediction_df, aes(x = time, y = prediction, color = factor(moved))) +
  geom_line() +
  geom_point() +
  scale_x_continuous(breaks = c(0, 1)) +
  labs(x = "Year", y = "Predicted z_participation", color = "Moved") +
  ggtitle("Interaction Plot of z_participation_base by Year and Moved")

```


### Fixed Effects


```{r}
#| echo: false
#| warning: false

models <- list()
models[['Simple']] = lm(z_participation ~ time , regd)
models[['Fixed Effects']] = lm(z_participation ~ time + factor(response_id), regd)
models[['Interaction']] = lm(z_participation ~ moved + time + moved*time + factor(response_id), regd)

modelsummary(
  models,
  coef_omit = "factor\\(.*",
  estimate  = "{estimate}{stars} ({std.error})",
             statistic = NULL,
  gof_omit = 'IC|RMSE|Log|F|R2$|Std.')


```



```{r}
#| echo: false
#| warning: false
#| 

# Predict values based on the model
predictions <- predict(lm(z_participation ~ moved + time + moved*time + factor(response_id), regd), newdata = regd)

# Create a dataframe with predictions
prediction_df <- data.frame(time = regd$time, moved = regd$moved, prediction = predictions) %>%
  group_by(time, moved) %>%
  summarise(prediction = mean(prediction))

# Plot the interaction
ggplot(prediction_df, aes(x = time, y = prediction, color = factor(moved))) +
  geom_line() +
  geom_point() +
  scale_x_continuous(breaks = c(0, 1)) +
  labs(x = "Year", y = "Predicted z_participation", color = "Moved") +
  ggtitle("Interaction Plot of z_participation_base by Year and Moved")

```