--- 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) ) ```