---
title: 'Code for QSS tidyverse Chapter 7: Uncertainty'
author: "Kosuke Imai and Nora Webb Williams"
date: "First Printing"
output:
  pdf_document: default
always_allow_html: true
---

# Uncertainty

```{r setup7, include = FALSE, purl=FALSE, eval = T}
library(knitr)
library(kableExtra)
library(tidyverse)
opts_chunk$set(tidy = FALSE,
        fig.align = "center",
        background = "#EEEEEE",
        fig.width = 7,
        fig.height = 7,
        out.width = ".7\\linewidth",
        out.height = ".7\\linewidth",
        cache = TRUE)
options(width = 60)
set.seed(12345)
```

## Estimation

### Unbiasedness and Consistency 

```{r echo = FALSE, warning=FALSE, message=FALSE}
library(tidyverse)
set.seed(1234)
```

```{r}
## simulation parameters
n <- 100 # sample size
mu0 <- 0 # mean of Y_i(0) [not treated]
sd0 <- 1 # standard deviation of Y_i(0)
mu1 <- 1 # mean of Y_i(1) [treated]
sd1 <- 1 # standard deviation of Y_i(1)

## generate a sample as a tibble
smpl <- tibble(id = seq_len(n),
               # Y if not treated
               Y0 = rnorm(n, mean = mu0, sd = sd0),
               # Y if treated
               Y1 = rnorm(n, mean = mu1, sd = sd1),
               # individual treatment effect
               tau = Y1 - Y0)

## true value of the sample average treatment effect
SATE <- smpl %>% select(tau) %>% summarize(SATE = mean(tau)) %>% pull()
SATE
```

```{r message=FALSE}
sim_treat <- function(smpl) {
  n <- nrow(smpl)
  # indexes of obs receiving treatment (randomly assign to half)
  idx <- sample(seq_len(n), floor(nrow(smpl) / 2), replace = FALSE)
  # "treat" variable is 1 for those receiving treatment, else 0
  smpl[["treat"]] <- as.integer(seq_len(nrow(smpl)) %in% idx)
  smpl %>%
    # what outcome we observe for unit based on treatment assignment
    mutate(Y_obs = if_else(treat == 1, Y1, Y0)) %>%
    group_by(treat) %>%
    summarize(Y_obs = mean(Y_obs)) %>%
    pivot_wider(names_from = treat, 
                values_from = Y_obs) %>%
    rename(Y1_mean = `1`, Y0_mean = `0`) %>%
    mutate(diff_mean = Y1_mean - Y0_mean,
           est_error = diff_mean - SATE)
}
## show the results of the function on the data
## values will differ each time it is run
sim_treat(smpl)
```

```{r, message=FALSE, cache=TRUE}
## number of simulations
sims <- 500
## run the created function sims times
sate_sims <- map_df(seq_len(sims), ~ sim_treat(smpl))
## what is the distribution of the error? 
summary(sate_sims$est_error)
```

```{r, message=FALSE, warning=FALSE, cache=TRUE}
PATE <- mu1 - mu0
PATE
## Update the function for PATE instead of SATE
sim_pate <- function(n, mu0, mu1, sd0, sd1) {
  smpl <- tibble(Y0 = rnorm(n, mean = mu0, sd = sd0),
                 Y1 = rnorm(n, mean = mu1, sd = sd1),
                 tau = Y1 - Y0)
  # indexes of obs receiving treatment
  idx <- sample(seq_len(n), floor(nrow(smpl) / 2), replace = FALSE)
  # treat variable are those receiving treatment, else 0
  smpl[["treat"]] <- as.integer(seq_len(nrow(smpl)) %in% idx)
  smpl %>%
    mutate(Y_obs = if_else(treat == 1L, Y1, Y0)) %>%
    group_by(treat) %>%
    summarize(Y_obs = mean(Y_obs)) %>%
    pivot_wider(names_from = treat, 
                values_from = Y_obs) %>%
    rename(Y1_mean = `1`, Y0_mean = `0`) %>%
    mutate(diff_mean = Y1_mean - Y0_mean,
           est_error = diff_mean - PATE)
}
## number of simulations
sims <- 500
## run the created function sims times
## input values are defined above
pate_sims <-  map_df(seq_len(sims), ~ sim_pate(n, mu0, mu1, sd0, sd1))
## what is the distribution of the error? 
summary(pate_sims$est_error)
```

### Standard Error

```{r include=FALSE, purl=FALSE}
## setting the ggplot theme
theme_set(theme_classic(base_size = 12))
```

```{r}
ggplot(pate_sims, aes(x = diff_mean, y = ..density..)) +
  geom_histogram(binwidth = 0.1) +
  geom_vline(xintercept = PATE, color = "black", size = 0.5) +
  ggtitle("Sampling distribution") +
  labs(x = "Difference-in-means estimator", y = "Density")
```

```{r}
## the standard deviation of the difference in means
pate_sims %>% 
  select(diff_mean) %>% 
  summarize(sd = sd(diff_mean))
```


```{r}
RMSE <- pate_sims %>% 
  summarize(rmse = sqrt(mean(est_error)^2))
RMSE
```

```{r message=FALSE, warning=FALSE}
## PATE simulation with standard error
sim_pate_se <- function(n, mu0, mu1, sd0, sd1) {
  # PATE - difference in means
  PATE <- mu1 - mu0
  # sample
  smpl <- tibble(Y0 = rnorm(n, mean = mu0, sd = sd0),
                 Y1 = rnorm(n, mean = mu1, sd = sd1),
                 tau = Y1 - Y0)
  # indexes of obs receiving treatment
  idx <- sample(seq_len(n), floor(nrow(smpl) / 2), replace = FALSE)
  # treat variable are those receiving treatment, else 0
  smpl[["treat"]] <- as.integer(seq_len(nrow(smpl)) %in% idx)
  # sample
  smpl %>%
    mutate(Y_obs = if_else(treat == 1, Y1, Y0)) %>%
    group_by(treat) %>%
    summarize(mean = mean(Y_obs),
              var = var(Y_obs),
              nobs = n()) %>%
    summarize(diff_mean = diff(mean),
              se = sqrt(sum(var / nobs)),
              est_error = diff_mean - PATE)
}
## test a single simulation
sim_pate_se(n, mu0, mu1, sd0, sd1)
## run 500 times
sims <- 500
pate_sims_se <- map_df(seq_len(sims), ~ sim_pate_se(n, mu0, mu1, sd0, sd1))

## standard deviation of difference-in-means
## and mean of standard errors
sd_se <- pate_sims_se %>% 
  summarize(sd = sd(diff_mean),
            mean_se = mean(se))
sd_se
```

### Confidence Interval

```{r}
## set the sample size
n <- 1000
## set the point estimate
x_bar <- 0.6
## calculate the standard error
se <- sqrt(x_bar * (1 - x_bar) / n)
## set the desired Confidence levels
levels <- c(0.99, 0.95, 0.90)
## build a tibble to calculate the ci at each level
tibble(level = levels) %>%
  mutate(
    ci_lower = x_bar - qnorm(1 - (1 - level) / 2) * se,
    ci_upper = x_bar + qnorm(1 - (1 - level) / 2) * se
  )

```


```{r}
## initial confidence level
level <- 0.95
## CI at that level for the PATE simulations with standard errors
pate_sims_ci <- pate_sims_se %>%
  mutate(ci_lower = diff_mean - qnorm(1 - (1 - level) / 2) * se,
         ci_upper = diff_mean + qnorm(1 - (1 - level) / 2) * se,
         includes_pate = PATE > ci_lower & PATE < ci_upper)
## view a subset of the CIs
glimpse(pate_sims_ci)

## compute the rate of PATE coverage 
pate_sims_ci %>% 
  summarize(coverage = mean(includes_pate))
```

```{r}
pate_sims_coverage <- function(.data, level = 0.95) {
  mutate(.data,
         ci_lower = diff_mean - qnorm(1 - (1 - level) / 2) * se,
         ci_upper = diff_mean + qnorm(1 - (1 - level) / 2) * se,
         includes_pate = PATE > ci_lower & PATE < ci_upper) %>%
    summarize(coverage = mean(includes_pate))
}

pate_sims_coverage(pate_sims_se, level = 0.95)
pate_sims_coverage(pate_sims_se, level = 0.99)
pate_sims_coverage(pate_sims_se, level = 0.90)

```

```{r}
## Function to test if CI contains true parameter value
binom_ci_contains <- function(n, p, alpha = 0.05) {
  x <- rbinom(n, size = 1, prob = p)
  x_bar <- mean(x)
  se <- sqrt(x_bar * (1 - x_bar) / n)
  ci_lower <- x_bar - qnorm(1 - alpha / 2) * se
  ci_upper <- x_bar + qnorm(1 - alpha / 2) * se
  (ci_lower <= p) & (p <= ci_upper)
}
## Demonstrate the function
p <- 0.6 # true parameter value
n <- 10
binom_ci_contains(n = n, p = p, alpha = 0.05)
```

```{r}
## Show coverage by taking the average of the logical result for each sim
mean(map_lgl(seq_len(sims), ~ binom_ci_contains(n, p)))
```

```{r}
## Function to calculate CI coverage while varying number of simulations
binom_ci_coverage <- function(n, p, sims) {
  mean(map_lgl(seq_len(sims), ~ binom_ci_contains(n, p)))
}
## Apply the function to a range of simulations values
tibble(n = c(10, 100, 1000)) %>%
  mutate(coverage = map_dbl(n, binom_ci_coverage, 
                            p = p, 
                            sims = sims))
```

### Margin of Error and Sample Size Calculation in Polls 

```{r}
## First, write a function to find population proportion give MoE
moe_pop_prop <- function(MoE) {
  tibble(p = seq(from = 0.01, to = 0.99, by = 0.01),
         n = 1.96 ^ 2 * p * (1 - p) / MoE ^ 2,
         MoE = MoE)
}
glimpse(moe_pop_prop(0.01))

## Then use map_df to call the function for a range of MoEs
MoE <- c(0.01, 0.03, 0.05)
props <- map_df(MoE, moe_pop_prop)

## plot the results
ggplot(props, aes(x = p, y = n, color = factor(MoE))) +
  geom_line(aes(linetype = factor(MoE))) +
  labs(color = "Margin of error",
       linetype = "Margin of error",
       x = "Population proportion",
       y = "Sample size") +
  theme(legend.position = "bottom")
```


```{r, message=FALSE, warning=FALSE}
## load required library
library(lubridate)
## load final vote shares
data("pres08", package = "qss")
## load polling data
data("polls08", package = "qss")
## set the election date
ELECTION_DATE <- ymd(20081104)
## Add days to election variable
polls08 <- polls08 %>%
  mutate(DaysToElection = as.integer(ELECTION_DATE - middate))
## Calculate mean of latest polls by state
poll_pred <- polls08 %>%
  group_by(state) %>%
  # latest polls in the state
  filter(DaysToElection == min(DaysToElection)) %>%
  # take mean of latest polls and convert from 0-100 to 0-1
  summarize(Obama = mean(Obama) / 100)
## Add confidence intervals
## sample size (assumed)
sample_size <- 1000
# confidence level
alpha <- 0.05
## Add the CIs and se
poll_pred <- poll_pred %>%
  mutate(se = sqrt(Obama * (1 - Obama) / sample_size),
         ci_lwr = Obama + qnorm(alpha / 2) * se,
         ci_upr = Obama + qnorm(1 - alpha / 2) * se)
## Add actual outcome
## And check if coverage includes the actual result
poll_pred <-left_join(poll_pred,
            select(pres08, state, actual = Obama),
            by = "state") %>%
  mutate(actual = actual / 100,
         covers = (ci_lwr <= actual) & (actual <= ci_upr))
## Check the results
glimpse(poll_pred)
```

```{r}
## Plot the results
ggplot(poll_pred, aes(x = actual, y = Obama)) +
  geom_abline(intercept = 0, slope = 1, color = "black", size = 0.5) +
  geom_pointrange(aes(ymin = ci_lwr, ymax = ci_upr),
                  shape = 1) +
  scale_y_continuous("Poll prediction", limits = c(0, 1)) +
  scale_x_continuous("Obama's vote share", limits = c(0, 1)) +
  scale_color_discrete("CI includes result?") +
  coord_fixed()
```

```{r}
poll_pred %>%
  summarize(mean(covers))
```

```{r}
poll_pred <- poll_pred %>%
  # calculate the bias
  mutate(bias = Obama - actual) %>%
  # bias corrected prediction, se, and CI
  mutate(Obama_bc = Obama - mean(bias),
         se_bc = sqrt(Obama_bc * (1 - Obama_bc) / sample_size),
         ci_lwr_bc = Obama_bc + qnorm(alpha / 2) * se_bc,
         ci_upr_bc = Obama_bc + qnorm(1 - alpha / 2) * se_bc,
         covers_bc = (ci_lwr_bc <= actual) & (actual <= ci_upr_bc))
## Updated coverage rate
poll_pred %>%
  summarize(mean(covers_bc))
```

### Analysis of Randomized Controlled Trials

```{r message=FALSE}
## load the data
data("STAR", package = "qss")
## Add meaningful labels to the classtype variable:
STAR <- STAR %>%
  mutate(classtype = factor(classtype,
                            labels = c("Small class", "Regular class",
                                       "Regular class with aide")))
## Summarize scores by classroom type:
classtype_means <- STAR %>%
  group_by(classtype) %>%
  summarize(g4reading = mean(g4reading, na.rm = TRUE))

## Plot the distribution of scores by classroom type for two classroom types
classtypes_used <- c("Small class", "Regular class")
ggplot(filter(STAR,
              classtype %in% classtypes_used,
              !is.na(g4reading)),
       aes(x = g4reading, y = ..density..)) +
  geom_histogram(binwidth = 20) +
  geom_vline(data = filter(classtype_means, classtype %in% classtypes_used),
             mapping = aes(xintercept = g4reading),
             color = "blue", size = 0.5) +
  facet_grid(classtype ~ .) +
  labs(x = "Fourth grade reading score", y = "Density")
```

```{r}
## alpha for 95% confidence
alpha <- 0.05
## calculate the mean, se, and CIs
star_estimates <- STAR %>%
  filter(!is.na(g4reading),
         classtype %in% classtypes_used) %>%
  group_by(classtype) %>%
  summarize(n = n(),
            est = mean(g4reading),
            se = sd(g4reading) / sqrt(n)) %>%
  mutate(lwr = est + qnorm(alpha / 2) * se,
         upr = est + qnorm(1 - alpha / 2) * se)

star_estimates
```

```{r}
## difference-in-means estimator
star_ate <- star_estimates %>%
  # ensure that it is ordered small then regular
  arrange(desc(classtype)) %>%
  summarize(
    se = sqrt(sum(se ^ 2)),
    est = diff(est)
  ) %>%
  mutate(ci_lwr = est + qnorm(alpha / 2) * se,
         ci_up = est + qnorm(1 - alpha / 2) * se)

star_ate
```

### Analysis based on Student's t Distribution

```{r}
## alpha for 95% confidence
alpha <- 0.05
## calculate the mean, se, and CIs
star_estimates_t <- STAR %>%
  filter(!is.na(g4reading),
         classtype %in% classtypes_used) %>%
  group_by(classtype) %>%
  summarize(n = n(),
            est = mean(g4reading),
            se = sd(g4reading) / sqrt(n)) %>%
  mutate(lwr = est + qt(alpha / 2, df = n - 1) * se,
         upr = est + qt(1 - alpha / 2, df = n - 1) * se)

star_estimates_t
## compare to original
star_estimates
```


```{r}
## compare reading scores between small and regular classes
reading_small <-  filter(STAR, classtype == "Small class")$g4reading
reading_reg <- filter(STAR, classtype == "Regular class")$g4reading

t_ci <- t.test(reading_small, reading_reg)
       
t_ci
```

## Hypothesis Testing

### Lady Tasting Tea Experiment

```{r}
## Number of cups of tea
cups <- 4
## Number guessed correctly
k <- c(0, seq_len(cups))
## Calculate probability correct
true <-  tibble(correct = k * 2,
                n = choose(cups, k) * choose(cups, cups - k)) %>%
  mutate(prob = n / sum(n))
true
```

```{r, cache=TRUE}
## Number of simulations
sims <- 1000
## The lady's guess (fixed); M for milk first; T for tea first
guess <- tibble(guess = c("M", "T", "T", "M", "M", "T", "T", "M"))

## A function to randomize the tea and calculate correct guesses
randomize_tea <- function(df) {
  # randomize the order of teas
  assignment <- sample_frac(df, 1) %>%
    rename(actual = guess)
  bind_cols(df, assignment) %>%
    summarize(correct = sum(guess == actual))
}

## Run the function 1000 times
approx <-
  map_df(seq_len(sims), ~ randomize_tea(guess)) %>%
  count(correct) %>%
  mutate(prob = n / sum(n))

## Then merge with the analytical solution
results <- approx %>% 
  select(correct, prob_sim = prob) %>% 
  left_join(select(true, correct, prob_exact = prob),
          by = "correct") %>%
  mutate(diff = prob_sim - prob_exact)

results
```

### The General Framework

```{r}
## all correct
x <- matrix(c(4, 0, 0, 4), byrow = TRUE, ncol = 2, nrow = 2)
## six correct
y <- matrix(c(3, 1, 1, 3), byrow = TRUE, ncol = 2, nrow = 2)
## `M' milk first, `T' tea first
rownames(x) <- colnames(x) <- rownames(y) <- colnames(y) <- c("M", "T")
x
y
```

```{r}
## one-sided test for 8 correct guesses
fisher.test(x, alternative = "greater") 
## two-sided test for 6 correct guesses
fisher.test(y)
```

### One-sample Tests

```{r}
n <- 1018
x.bar <- 550 / n
se <- sqrt(0.5 * 0.5 / n) # standard deviation of sampling distribution
## upper red area in the figure
upper <- pnorm(x.bar, mean = 0.5, sd = se, lower.tail = FALSE)  
## lower red area in the figure; identical to the upper area
lower <- pnorm(0.5 - (x.bar - 0.5), mean = 0.5, sd = se)  
## two-side p-value
upper + lower
```

```{r}
2 * upper
```

```{r}
## one-sided p-value
upper
```


```{r}
z.score <- (x.bar - 0.5) / se
z.score
pnorm(z.score, lower.tail = FALSE) # one-sided p-value
2 * pnorm(z.score, lower.tail = FALSE) # two-sided p-value
```

```{r}
## 99% confidence interval contains 0.5
c(x.bar - qnorm(0.995) * se, x.bar + qnorm(0.995) * se)
## 95% confidence interval does not contain 0.5
c(x.bar - qnorm(0.975) * se, x.bar + qnorm(0.975) * se)
```

```{r}
## no continuity correction to get the same p-value as above
prop.test(550, n = n, p = 0.5, correct = FALSE)
## with continuity correction
prop.test(550, n = n, p = 0.5)
```

```{r}
prop.test(550, n = n, p = 0.5, conf.level = 0.99)
```

```{r}
# two-sided one-sample t-test
t.test(STAR$g4reading, mu = 710)
```

### Two-sample Tests

```{r}
star_ate %>%
  mutate(p_value_1sided = pnorm(-abs(est),
                                mean = 0, sd = se),
         p_value_2sided = 2 * pnorm(-abs(est), mean = 0,
                                    sd = se))

```

```{r}
## testing the null of zero average treatment effect
reading_small <- filter(STAR, classtype == "Small class")$g4reading
reading_reg <- filter(STAR, classtype == "Regular class")$g4reading
## t-test
t.test(reading_small,
       reading_reg)
```

```{r}
## load the data
data("resume", package = "qss")
## reshape the data
x <- resume %>%
  count(race, call) %>%
  pivot_wider(names_from = call, values_from = n) %>%
  ungroup()
x
## run the test on the relevant columns
prop.test(as.matrix(select(x, -race)), alternative = "greater")
```

```{r}
## sample size
n0 <- sum(resume$race == "black")
n1 <- sum(resume$race == "white")

## sample proportions
p <- mean(resume$call) # overall
p0 <- mean(filter(resume, race == "black")$call)
p1 <- mean(filter(resume, race == "white")$call)

## point estimate
est <- p1 - p0
est

## standard error
se <- sqrt(p * (1 - p) * (1 / n0 + 1 / n1))
se

## z-statistic
zstat <- est / se
zstat

## one-sided p-value
pnorm(-abs(zstat))
```

```{r}
prop.test(as.matrix(select(x, -race)), alternative = "greater", correct = FALSE)
```

### Pitfalls of Hypothesis Testing

### Power Analysis

```{r}
## set the parameters
n <- 250
p.star <- 0.48 # data generating process
p <- 0.5 # null value
alpha <- 0.05
## critical value
cr.value <- qnorm(1 - alpha / 2) 
## standard errors under the hypothetical data generating process
se.star <- sqrt(p.star * (1 - p.star) / n)
## standard error under the null
se <- sqrt(p * (1 - p) / n)  
## power
pnorm(p - cr.value * se, mean = p.star, sd = se.star) + 
  pnorm(p + cr.value * se, mean = p.star, sd = se.star, lower.tail = FALSE)
```

```{r}
## parameters
n1 <- 500
n0 <- 500
p1.star <- 0.05
p0.star <- 0.1
```

```{r}
## overall call back rate as a weighted average
p <- (n1 * p1.star + n0 * p0.star) / (n1 + n0) 
## standard error under the null
se <- sqrt(p * (1 - p) * (1 / n1 + 1 / n0)) 
## standard error under the hypothetical data generating process
se.star <- sqrt(p1.star * (1 - p1.star) / n1 
                + p0.star * (1 - p0.star) / n0)
```

```{r}
pnorm(-cr.value * se, mean = p1.star - p0.star, sd = se.star) +
  pnorm(cr.value * se, mean = p1.star - p0.star, sd = se.star, 
          lower.tail = FALSE)
```

```{r}
power.prop.test(n = 500, p1 = 0.05, p2 = 0.1, sig.level = 0.05)
```

```{r}
power.prop.test(p1 = 0.05, p2 = 0.1, sig.level = 0.05, power = 0.9)
```

```{r}
power.t.test(n = 100, delta = 0.25, sd = 1, type = "one.sample")
```


```{r}
power.t.test(power = 0.9, delta = 0.25, sd = 1, type = "one.sample")
```


```{r}
power.t.test(delta = 0.25, sd = 1, type = "two.sample", 
             alternative = "one.sided", power = 0.9)
```

## Linear Regression Model with Uncertainty

### Linear Regression as a Generative Model 

```{r}
## load the data
data("minwage", package = "qss")
## compute proportion of full employment before minimum wage increase
## same thing after minimum wage increase
## an indicator for NJ: 1 if it's located in NJ and 0 if in PA
minwage <- 
  mutate(minwage,
         fullPropBefore = fullBefore / (fullBefore + partBefore),
         fullPropAfter = fullAfter / (fullAfter + partAfter),
         NJ = if_else(location == "PA", 0 , 1))
```


```{r}
fit_minwage <- lm(fullPropAfter ~ -1 + NJ + fullPropBefore +
                    wageBefore + chain, data = minwage)
## regression result
fit_minwage
```

```{r}
## with intercept
fit_minwage1 <- lm(fullPropAfter ~ NJ + fullPropBefore +
                     wageBefore + chain, data = minwage)
fit_minwage1
```


```{r, message=FALSE, warning=FALSE}
## load the required package
library(modelr)
## Generate prediction from the first model
pred_1 <-minwage %>% 
  slice(1) %>% 
  add_predictions(fit_minwage) %>%
  select(pred) %>% 
  mutate(model = "fit_minwage")
## Generate prediction from the second model
## then add predictions from first to compare
pred_compare <-minwage %>% 
  slice(1) %>% 
  add_predictions(fit_minwage1) %>%
  select(pred) %>% 
  mutate(model = "fit_minwage1") %>% 
  bind_rows(pred_1)

pred_compare  
```
                  
### Unbiasedness of Estimated Coefficients

### Standard Errors of Estimated Coefficients

### Inference about Coefficients

```{r message=FALSE, warning=FALSE}
library(broom)
## load the data
data("women", package = "qss")
## fit the model
fit_women <- lm(water ~ reserved, data = women)
## view the coefficients 
summary(fit_women)
tidy(fit_women)
```

```{r}
## display confidence intervals
tidy(fit_women, conf.int = TRUE)
```

```{r}
summary(fit_minwage)
tidy(fit_minwage, conf.int = TRUE)
```

### Inference about Predictions

```{r}
## load the data and subset them into two parties
data("MPs", package = "qss")
MPs_labour <- filter(MPs, party == "labour")
MPs_tory <- filter(MPs, party == "tory")

## two regressions for labour: negative and positive margin
labour_fit1 <- lm(ln.net ~ margin, data = filter(MPs_labour, margin < 0))
labour_fit2 <- lm(ln.net ~ margin, data = filter(MPs_labour, margin > 0))

## two regressions for tory: negative and positive margin
tory_fit1 <- lm(ln.net ~ margin, data = filter(MPs_tory, margin < 0))
tory_fit2 <- lm(ln.net ~ margin, data = filter(MPs_tory, margin > 0))
```

```{r}
## tory party: prediction at the threshold
tory_y0 <- augment(tory_fit1, newdata = tibble(margin = 0), 
          interval = "confidence",
          conf.level = 0.95)

tory_y0
tory_y1 <- augment(tory_fit2, newdata = tibble(margin = 0), 
                   interval = "confidence")
tory_y1
```

```{r}
## create data with the ranges of "margin" less than zero
y1_range <- data_grid(filter(MPs_tory, margin <= 0), margin)
## add predictions and CIs
tory_y0 <- augment(tory_fit1, newdata = y1_range, 
                   interval = "confidence")
## create the data with the ranges of "margin" greater than zero
y2_range <- data_grid(filter(MPs_tory, margin >= 0), margin)
## add predictions and CIs
tory_y1 <- augment(tory_fit2, newdata = y2_range, 
                   interval = "confidence")
```

```{r include=FALSE, purl=FALSE}
## setting the ggplot theme
theme_set(theme_classic(base_size = 12))
```

```{r}
ggplot() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  # plot losers
  geom_ribbon(aes(x = margin, ymin = .lower, ymax = .upper),
              data = tory_y0, alpha = 0.3) +
  geom_line(aes(x = margin, y = .fitted), data = tory_y0) +
  # plot winners
  geom_ribbon(aes(x = margin, ymin = .lower, ymax = .upper),
              data = tory_y1, alpha = 0.3) +
  geom_line(aes(x = margin, y = .fitted), data = tory_y1) +
  xlim(-.5,.25) +
  labs(x = "Margin of victory", y = "log net wealth")
```


```{r}
## predictions at threshold with SEs
tory_y0 <- augment(tory_fit1, newdata = tibble(margin = 0), 
          interval = "confidence",
          se_fit = TRUE)

tory_y0
tory_y1 <- augment(tory_fit2, newdata = tibble(margin = 0), 
                   interval = "confidence",
                   se_fit = TRUE)
tory_y1
```

```{r}
summary(tory_fit2)
```

```{r}
# standard error
se_diff <- sqrt(tory_y0$.se.fit ^ 2 + tory_y1$.se.fit ^ 2)
se_diff
# point estimate
diff_est <- tory_y1$.fitted - tory_y0$.fitted
diff_est
# confidence interval
CI <- c(diff_est - se_diff * qnorm(0.975),
        diff_est + se_diff * qnorm(0.975))
CI
# hypothesis test
z_score <- diff_est / se_diff
# two sided p value
p_value <- 2 * pnorm(abs(z_score), lower.tail = FALSE)
p_value
```

