---
title: "Lab 5: Practice Key"
subtitle: "Binary Models and Intro to tile"
author: Minji Jeong
date: \today
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, error=FALSE, warning=FALSE, message = FALSE)

# clean environment

rm(list=ls())

# Load packages
library(MASS)
library(simcf)
library(tile)
library(RColorBrewer)
source("source/theme_cavis.R") #ggplot theme from former TA Brian Leung; Download ZIP file if you want to use it
library(WhatIf)
library(tidyverse)  # Should come after MASS
```



# Visualization with Tile


Consider a new specification by adding the variable *ever married*.

We’ll estimate this model with `glm` and simulate expected values for new scenarios comparing married and non-married respondents.


```{r}
# Load data
data <- read.csv("https://faculty.washington.edu/cadolph/mle/nes00a.csv")

# Estimate logit model using glm()

# Set up model formula and model specific data frame
model2 <- vote00 ~ age + I(age^2) + hsdeg + coldeg + marriedo
mdata2 <- extractdata(model2, data, na.rm=TRUE)

# Run logit & extract results
logitM2 <- glm(model2, family=binomial(link="logit"), data=mdata2)
peM2 <- logitM2$coefficients  # point estimates
vcM2 <- vcov(logitM2)         # var-cov matrix

```




## Quantities of Interest using `simcf`


Now, we will compute quantities of interest (QoI) again. Specifically, we will calculate the **first differences** and **relative risks** between the scenarios of being married and not being married. This allows us to visualize the marginal effect of marriage on the expected probability and relative risk of voting.

When using `simcf` to compute first differences or relative risks with `simf`, we need to set up two parallel scenarios, `x` and `xpre`, in the `cfMake` object. The first, `x`, represents the current value, while the second, `xpre`, represents the prior value. For example, the first difference is calculated as `x - xpre`.



```{r}
# Simulate parameter distributions
sims <- 10000
simbetasM2 <- mvrnorm(n=sims, 
                      mu=peM2, 
                      Sigma=vcM2)


# Set up counterfactuals:  all ages, each of three educations
xhyp <- seq(18,97,1)
nscen <- length(xhyp)

marriedScen <- notmarrScen <- cfMake(formula=model2,
                                     data=mdata2,
                                     nscen=nscen)


for (i in 1:nscen) {
  
  # Married (loop over each age)
  # Note below the careful use of before scenarios (xpre) and after scenarios (x)
  # we will use the marriedScen counterfactuals in FDs and RRs as well as EVs
  marriedScen <- cfChange(xscen=marriedScen, 
                          covname="age", 
                          x = xhyp[i], 
                          xpre= xhyp[i],  
                          scen = i)
  
  marriedScen <- cfChange(xscen=marriedScen, 
                          covname="marriedo", 
                          x = 1,     # compare before and after
                          xpre= 0, 
                          scen = i)
  # We want to see the difference from not being married to married 
  # you can hold other things (hs, coll) at their means
  # we change age bc we are interested in age
  
  # Not married (loop over each age)
  notmarrScen <- cfChange(xscen=notmarrScen, 
                          covname="age",
                          x = xhyp[i], 
                          scen = i)
  
  notmarrScen <- cfChange(xscen=notmarrScen, 
                          covname="marriedo", 
                          x = 0, 
                          scen = i)
}

# Simulate expected probabilities for all scenarios
marriedSims <- logitsimev(marriedScen, simbetasM2, ci=0.95)
notmarrSims <- logitsimev(notmarrScen, simbetasM2, ci=0.95)

# Simulate first difference of voting with respct to marriage
marriedFD <- logitsimfd(marriedScen, simbetasM2, ci=0.95)

# Simulate relative risk of voting wrt marriage
marriedRR <- logitsimrr(marriedScen, simbetasM2, ci=0.95)

```



### Visualization: `tile`


Below we execute the plot using `tile`.

```{r, results='hide'}
## Make plots using tile

# Get 3 nice colors for traces
col <- brewer.pal(3,"Dark2")

# Set up lineplot traces of expected probabilities
marriedTrace <- lineplot(x=xhyp,
                         y=marriedSims$pe,
                         lower=marriedSims$lower,
                         upper=marriedSims$upper,
                         col=col[1],
                         extrapolate=list(data=mdata2[,2:ncol(mdata2)],
                                          cfact=marriedScen$x[,2:ncol(marriedScen$x)],
                                          omit.extrapolated=TRUE),
                         plot=1)

notmarrTrace <- lineplot(x=xhyp,
                         y=notmarrSims$pe,
                         lower=notmarrSims$lower,
                         upper=notmarrSims$upper,
                         col=col[2],
                         ci = list(mark="dashed"),
                         extrapolate=list(data=mdata2[,2:ncol(mdata2)],
                                          cfact=notmarrScen$x[,2:ncol(notmarrScen$x)],
                                         omit.extrapolated=TRUE),
                         plot=1)


# Set up traces with labels and legend
labelTrace <- textTile(labels=c("Currently Married", "Not Married"),
                       x=c( 35,    53),
                       y=c( 0.8,  0.56),
                       col=col,
                       plot=1)

legendTrace <- textTile(labels=c("Logit estimates:", "95% confidence", "interval is shaded"),
                        x=c(80, 80, 80),
                        y=c(0.2, 0.15, 0.10),
                        cex=0.9,
                        plot=1)

# Plot traces using tile
tile(marriedTrace,
     notmarrTrace,
     labelTrace,
     legendTrace,
     limits=c(18,94,0,1),
     xaxis=list(at=c(20,30,40,50,60,70,80,90), fontsize=11),
     yaxis=list(label.loc=-0.5, major=FALSE, fontsize=11),
     xaxistitle=list(labels="Age of Respondent", y=-0.5),
     yaxistitle=list(labels="Probability of Voting"),
     width=list(null=5,yaxistitle=4,yaxis.labelspace=-0.5)
     )




# Plot Relative Risk

# Set up lineplot trace of relative risk
marriedRRTrace <- lineplot(x=xhyp,
                         y=marriedRR$pe,
                         lower=marriedRR$lower,
                         upper=marriedRR$upper,
                         col=col[1],
                         extrapolate=list(data=mdata2[,2:ncol(mdata2)],
                                          cfact=marriedScen$x[,2:ncol(marriedScen$x)],
                                          omit.extrapolated=TRUE),
                         plot=1)


# Set up baseline: for relative risk, this is 1
baseline <- linesTile(x=c(18,94),
                      y=c(1,1),
                      plot=1)

# Set up traces with labels and legend
labelRRTrace <- textTile(labels=c("Married compared \n to Not Married"),
                         x=c( 55),
                         y=c( 1.25),
                         col=col[1],
                         plot=1)

legendRRTrace <- textTile(labels=c("Logit estimates:", "95% confidence", "interval is shaded"),
                          x=c(80, 80, 80),
                          y=c(0.98, 0.95, 0.92),
                          cex=0.9,
                          plot=1)

# Plot traces using tile
tile(marriedRRTrace,
     labelRRTrace,
     legendRRTrace,
     baseline,
     limits=c(18,94,0.9,1.5),
     xaxis=list(at=c(20,30,40,50,60,70,80,90), fontsize=11),
     yaxis=list(label.loc=-0.5, major=FALSE, fontsize=11),
     xaxistitle=list(labels="Age of Respondent", y=-0.5),
     yaxistitle=list(labels="Relative Risk of Voting"),
     width=list(null=5,yaxistitle=4,yaxis.labelspace=-0.5)
     )



```


