---
title: "Lab 5"
output:
  html_document:
    theme: readable
    df_print: paged
    highlight: tango
    toc: yes
    toc_float: no
---

```{r setup, include=FALSE}
knitr::opts_chunk$set( echo=T, message=F, warning=F, image.width=10)
```










# Load Data

```{r}
library( dplyr )
library( scales )
library( stargazer )
library( pander )

# STARGAZER OUTPUT
#
# Use:
#
#   s.type="text"
#
# while running chunks interactively
# to see table output.
# This sets it to "html"
# when knitting the file.
s.type="html"
```


```{r, eval=F, echo=F}
# run this chunk to preview stargazer tables prior to knitting
s.type="text"
```


```{r}
URL <- "https://github.com/DS4PS/cpp-524-sum-2021/blob/main/labs/data/counterfactuals.csv?raw=true"
d <- read.csv( URL )

# set factor levels
d$group <- factor( d$group, 
                   levels=c("high.ses","treatment","control"))
```





```{r, fig.height=8, fig.width=8, echo=F}

t <- tapply( d$ability, list(d$group,d$time), mean )

x <- c(42,50,64,78)
plot.new()
plot.window( xlim=c(40,90), ylim=c(-2,5) )
abline( h=-2:5, lty=2, col="gray70" )
abline( h=-1.5:4.5, lty=3, col="gray85" )
points( x, t[1,], type="b", pch=19, col="steelblue", cex=3 )
points( x, t[2,], type="b", pch=19, col="darkred", cex=3 )
points( x, t[3,], type="b", pch=19, col="steelblue", cex=3 )
axis( side=1 )
axis( side=2, las=1 )
title( xlab="Age (months)", ylab="Ability (logits)" )
title( main="Average Outcomes of Study Groups" )

text( 80, t[1,4], "High SES", col="steelblue", pos=4 )
text( 80, t[2,4], "Treatment Group", col="darkred", pos=4 )
text( 80, t[3,4], "Control Group", col="steelblue", pos=4 )

segments( x0=50, y0=-2, y1=3.2, lty=3, col="darkred" )
text( 50, 3.7, "Start of \nTreatment \nPeriods", col="darkred" )

```





# Reflexive Pre-Post Estimator (T2-T1)



```{r, fig.height=8, fig.width=8, echo=F}

t <- tapply( d$ability, list(d$group,d$time), mean )

x <- c(42,50,64,78)
plot.new()
plot.window( xlim=c(40,90), ylim=c(-2,5) )

segments( x0=50, y0=-2, y1=3.1, lty=3, col="gray50" )
text( 50, 3.5, "Start of \nTreatment \nPeriods", col="gray50" )

abline( h=-2:5, lty=2, col="gray70" )
abline( h=-1.5:4.5, lty=3, col="gray85" )
points( x, t[1,], type="b", pch=19, col="steelblue", cex=3 )
points( x, t[2,], type="b", pch=19, col="steelblue", cex=3 )
points( x, t[3,], type="b", pch=19, col="steelblue", cex=3 )
axis( side=1 )
axis( side=2, las=1 )
title( xlab="Age (months)", ylab="Ability (logits)" )
title( main="Average Outcomes of Study Groups" )

text( 80, t[1,4], "High SES", col="steelblue", pos=4 )
text( 80, t[2,4], "Treatment Group", col="steelblue", pos=4 )
text( 80, t[3,4], "Control Group", col="steelblue", pos=4 )



points( x[2:3], t[2,2:3], 
        type="b", pch=19, col="darkred", cex=4, lwd=2, lty=3 )
```


```{r, results="asis"}
dm <- filter( d, 
              group %in% c("treatment") &
              time %in% c("time1","time2") )

dm$post.dummy <- ifelse( dm$time=="time2", 1, 0 )

m <- lm( ability ~ post.dummy, data=dm )

stargazer( m, type=s.type, 
           omit.stat=c("f","ser","adj.rsq"),
           intercept.top=TRUE, intercept.bottom=FALSE,
           digits=2 )
```

## Question 1a

> What is the average score for kids in the treatment group in period 1 of the study? What is the average score for the same kids in period 2?

*Note, you can check your answers against the group means in Table 2 above.*

**ANSWER:**



## Question 1b

> Explain what it means when b0 is statistically significant in the reflexive model. Explain what it means when b1 is statistically significant.

**ANSWER:**


## Question 1c

> What is the effect size according to this model?

**ANSWER:**


<br>
<hr>
<br>


Test for zero trend assumption: C1=C2

```{r, fig.height=6, fig.width=7, echo=F}

x <- c(42,50,64,78)
plot.new()
plot.window( xlim=c(40,90), ylim=c(-2,5) )

segments( x0=50, y0=-2, y1=3.1, lty=3, col="gray50" )
text( 50, 3.5, "Start of \nTreatment \nPeriods", col="gray50" )

abline( h=-2:5, lty=2, col="gray70" )
abline( h=-1.5:4.5, lty=3, col="gray85" )
points( x, t[1,], type="b", pch=19, col="steelblue", cex=3 )
points( x, t[2,], type="b", pch=19, col="steelblue", cex=3 )
points( x, t[3,], type="b", pch=19, col="steelblue", cex=3 )
axis( side=1 )
axis( side=2, las=1 )
title( xlab="Age (months)", ylab="Ability (logits)" )
title( main="Average Outcomes of Study Groups" )

text( 80, t[1,4], "High SES", col="steelblue", pos=4 )
text( 80, t[2,4], "Treatment Group", col="steelblue", pos=4 )
text( 80, t[3,4], "Control Group", col="steelblue", pos=4 )

# segments( x0=50, y0=t[1,2], y1=3.2, lty=3, col="steelblue" )
# text( 50, 3.7, "Start of \nTreatment \nPeriods", col="steelblue" )

points( x[2:3], t[3,2:3], 
        type="b", pch=19, col="darkred", cex=4, lwd=2, lty=3 )
```

```{r, results="asis"}
dm <- filter( d, 
              group %in% c("control") &
              time %in% c("time1","time2") )

dm$post.dummy <- ifelse( dm$time=="time2", 1, 0 )

m <- lm( ability ~ post.dummy, data=dm )

stargazer( m, type=s.type, 
           omit.stat=c("f","ser","adj.rsq"),
           intercept.top=TRUE, intercept.bottom=FALSE,
           digits=2 )
```


## Question 2a

> What is the average score for kids in the control group in period 1 of the study? What is the average score for the same kids in period 2?

**ANSWER:**




## Question 2b

> Which coefficient represents the test for whether we observe a secular trend in student achievement gains independent of the treatment? What is the decision rule?

**ANSWER:**



## Question 2c

> What does this model tell us about the appropriateness of the reflexive model?

**ANSWER:**


<br>
<hr>
<br>



# Post-Test Only Estimator (T2-C2)


```{r, fig.height=8, fig.width=8, echo=F}

t <- tapply( d$ability, list(d$group,d$time), mean )

x <- c(42,50,64,78)
plot.new()
plot.window( xlim=c(40,90), ylim=c(-2,5) )

segments( x0=50, y0=-2, y1=3.1, lty=3, col="gray50" )
text( 50, 3.5, "Start of \nTreatment \nPeriods", col="gray50" )

abline( h=-2:5, lty=2, col="gray70" )
abline( h=-1.5:4.5, lty=3, col="gray85" )
points( x, t[1,], type="b", pch=19, col="steelblue", cex=3 )
points( x, t[2,], type="b", pch=19, col="steelblue", cex=3 )
points( x, t[3,], type="b", pch=19, col="steelblue", cex=3 )
axis( side=1 )
axis( side=2, las=1 )
title( xlab="Age (months)", ylab="Ability (logits)" )

text( 80, t[1,4], "High SES", col="steelblue", pos=4 )
text( 80, t[2,4], "Treatment Group", col="steelblue", pos=4 )
text( 80, t[3,4], "Control Group", col="steelblue", pos=4 )

# segments( x0=50, y0=t[1,2], y1=3.2, lty=3, col="steelblue" )
# text( 50, 3.7, "Start of \nTreatment \nPeriods", col="steelblue" )

points( c(64,64), t[2:3,3], 
        type="b", pch=19, col="darkred", cex=4, lwd=2, lty=3 )
```



```{r, results="asis"}
dm <- filter( d, 
              group %in% c("treatment","control") &
              time=="time2" )

dm$treat.dummy <- ifelse( dm$group=="treatment", 1, 0 )

m <- lm( ability ~ treat.dummy, data=dm )

stargazer( m, type=s.type, 
           omit.stat=c("f","ser","adj.rsq"),
           intercept.top=TRUE, intercept.bottom=FALSE,
           digits=2 )
```


## Question 3a

> What is the average score for kids in the control group in the study? What is the average score for the kids in the treatment group? 

*Note, you can check your answers against the group means in Table 2 above.*

**ANSWER:**


## Question 3b

> What is the effect size identified by the model?


## Question 3c

> What is the identifying assumption of this model? Or stated differently, what must be true in order for the post-test only estimator to be appropriate?

**ANSWER:**

## Question 3d

> According to the model below is the assumption met? How can you tell?

**ANSWER:**


```{r, results="asis"}
dm <- filter( d, 
              group %in% c("treatment","control") &
              time=="time1" )

dm$treat.dummy <- ifelse( dm$group=="treatment", 1, 0 )

m <- lm( ability ~ treat.dummy, data=dm )

stargazer( m, type=s.type, 
           omit.stat=c("f","ser","adj.rsq"),
           intercept.top=TRUE, intercept.bottom=FALSE,
           digits=2 )
```

<br>
<hr>
<br>



# Diff-in-Diff Estimator [ gains - trend ]

Total Gains: T2-T1  
Trend: C2-C1  
DID Estimator: [ gains - trend ] = [ (T2-T1) - (C2-C1) ]   


```{r, fig.height=8, fig.width=8, echo=F}

t <- tapply( d$ability, list(d$group,d$time), mean )

x <- c(42,50,64,78)
plot.new()
plot.window( xlim=c(40,90), ylim=c(-2,5) )

segments( x0=50, y0=-2, y1=3.1, lty=3, col="gray50" )
text( 50, 3.5, "Start of \nTreatment \nPeriods", col="gray50" )

abline( h=-2:5, lty=2, col="gray70" )
abline( h=-1.5:4.5, lty=3, col="gray85" )
points( x, t[1,], type="b", pch=19, col="steelblue", cex=3 )
points( x, t[2,], type="b", pch=19, col="steelblue", cex=3 )
points( x, t[3,], type="b", pch=19, col="steelblue", cex=3 )
axis( side=1 )
axis( side=2, las=1 )
title( xlab="Age (months)", ylab="Ability (logits)" )

text( 80, t[1,4], "High SES", col="steelblue", pos=4 )
text( 80, t[2,4], "Treatment Group", col="steelblue", pos=4 )
text( 80, t[3,4], "Control Group", col="steelblue", pos=4 )

# segments( x0=50, y0=t[1,2], y1=3.2, lty=3, col="steelblue" )
# text( 50, 3.7, "Start of \nTreatment \nPeriods", col="steelblue" )

points( x[2:3], t[2,2:3], 
        type="b", pch=19, col="darkred", cex=4, lwd=2, lty=3 )
points( x[2:3], t[3,2:3], 
        type="b", pch=19, col="darkred", cex=4, lwd=2, lty=3 )
```



```{r, results="asis"}
dm <- filter( d, group %in% c("treatment","control") &
                time %in% c("time1","time2") )

dm$treat.dummy <- ifelse( dm$group=="treatment", 1, 0 )
dm$post.dummy  <- ifelse( dm$time=="time2", 1, 0 )

dm$treat.post.dummy <- dm$treat.dummy * dm$post.dummy

m <- lm( ability ~ treat.dummy + post.dummy + treat.post.dummy,
         data=dm)

stargazer( m, type=s.type, 
           omit.stat=c("f","ser","adj.rsq"),
           intercept.top=TRUE, intercept.bottom=FALSE,
           digits=2 )
```


```{r}
b0 <- m$coefficients[1] %>% as.numeric() %>% round(2)
b1 <- m$coefficients[2] %>% as.numeric() %>% round(2)
b2 <- m$coefficients[3] %>% as.numeric() %>% round(2)
b3 <- m$coefficients[4] %>% as.numeric() %>% round(2)

b0                 # C1
b0 + b1            # T1
b0 + b2            # C2
b0 + b1 + b2 + b3  # T2
b0 + b1 + b2       # CF

CF <- b0 + b1 + b2
T2 <- b0 + b1 + b2 + b3   
T2-CF
```


## Question 4a

> Are the treatment and control groups equivalent prior to the intervention? How do you know?

**ANSWER:**


## Question 4b

> Do we observe secular trends (gains independent of the treatment)? How do you know?

**ANSWER:**


## Question 4c

> What is the effect size in this model (gains from the treatment)?

**ANSWER:**


## Question 4d

> What does statistical significance of b3 represent? In other words, which contrast is being tested?

**ANSWER:**



<br>
<hr>
<br>


## Question 5a

> Do the reflexive and diff-in-diff models generate the same results (approximately)?  Why?

**ANSWER:**


## Question 5b

> Do the post-test only and diff-in-diff models generate the same results (approximately)? Why?

**ANSWER:**



<br>
<hr>
<br>



# Alternative Diff-in-Diff Estimator 

```{r, fig.height=8, fig.width=8, echo=F}

x <- c(42,50,64,78)
plot.new()
plot.window( xlim=c(40,90), ylim=c(-2,5) )

segments( x0=50, y0=-2, y1=3.1, lty=3, col="gray50" )
text( 50, 3.5, "Start of \nTreatment \nPeriods", col="gray50" )

abline( h=-2:5, lty=2, col="gray70" )
abline( h=-1.5:4.5, lty=3, col="gray85" )
points( x, t[1,], type="b", pch=19, col="steelblue", cex=3 )
points( x, t[2,], type="b", pch=19, col="steelblue", cex=3 )
points( x, t[3,], type="b", pch=19, col="steelblue", cex=3 )
axis( side=1 )
axis( side=2, las=1 )
title( xlab="Age (months)", ylab="Ability (logits)" )

text( 80, t[1,4], "High SES", col="steelblue", pos=4 )
text( 80, t[2,4], "Treatment Group", col="steelblue", pos=4 )
text( 80, t[3,4], "Control Group", col="steelblue", pos=4 )

# segments( x0=50, y0=t[1,2], y1=3.2, lty=3, col="steelblue" )
# text( 50, 3.7, "Start of \nTreatment \nPeriods", col="steelblue" )

points( x[2:3], t[2,2:3], 
        type="b", pch=19, col="darkred", cex=4, lwd=2, lty=3 )
points( x[2:3], t[1,2:3], 
        type="b", pch=19, col="darkred", cex=4, lwd=2, lty=3 )
```

```{r, results="asis"}
dm <- filter( d, group %in% c("treatment","high.ses") &
                time %in% c("time1","time2") )

dm$treat.dummy <- ifelse( dm$group=="treatment", 1, 0 )
dm$pre.dummy   <- ifelse( dm$time=="time1", 1, 0 )
dm$post.dummy  <- ifelse( dm$time=="time2", 1, 0 )

dm$treat.post.dummy <- dm$treat.dummy * dm$post.dummy 

m <- lm( ability ~ treat.dummy + post.dummy + treat.post.dummy,
         data=dm)

stargazer( m, type=s.type,
           omit.stat=c("f","ser","adj.rsq"),
           intercept.top=TRUE, intercept.bottom=FALSE,
           digits=2 )
```


## Question 6a

> Are the pre-treatment differences (C1=T1?) different in this model versus the previous diff-in-diff? Why or why not?

**ANSWER:**



## Question 6b

> Does the diff-in-diff model require that study groups are equivalent prior to treatment to generate valid results?

**ANSWER:**


## Question 6c

> Is the secular trend identified by this model different from the previous diff-in-diff (approximately)? Why or why not?

**ANSWER:**



## Question 6d

> The treatment effects from this model are approximately the same as the previous diff-in-diff model, even though they use very different comparison groups. Why does this model still work using the high SES group?

**ANSWER:**



## Question 6e

> What is the identification assumption of the diff-in-diff model?

**ANSWER:**



<br>
<hr>
<br>

# Parallel Lines Test

Does the high SES group model secular trend appropriately?

Test: are the study group trend lines parallel prior to the intervention?

```{r, fig.height=8, fig.width=8, echo=F}
x <- c(42,50,64,78)
plot.new()
plot.window( xlim=c(40,90), ylim=c(-2,5) )

segments( x0=50, y0=-2, y1=3.1, lty=3, col="gray50" )
text( 50, 3.5, "Start of \nTreatment \nPeriods", col="gray50" )

abline( h=-2:5, lty=2, col="gray70" )
abline( h=-1.5:4.5, lty=3, col="gray85" )
points( x, t[1,], type="b", pch=19, col="steelblue", cex=3 )
points( x, t[2,], type="b", pch=19, col="steelblue", cex=3 )
points( x, t[3,], type="b", pch=19, col="steelblue", cex=3 )
axis( side=1 )
axis( side=2, las=1 )
title( xlab="Age (months)", ylab="Ability (logits)" )

text( 80, t[1,4], "High SES", col="steelblue", pos=4 )
text( 80, t[2,4], "Treatment Group", col="steelblue", pos=4 )
text( 80, t[3,4], "Control Group", col="steelblue", pos=4 )

# segments( x0=50, y0=t[1,2], y1=3.2, lty=3, col="steelblue" )
# text( 50, 3.7, "Start of \nTreatment \nPeriods", col="steelblue" )

points( x[1:2], t[2,1:2], 
        type="b", pch=19, col="darkred", cex=4, lwd=2, lty=3 )
points( x[1:2], t[1,1:2], 
        type="b", pch=19, col="darkred", cex=4, lwd=2, lty=3 )
```

```{r, results="asis"}
dm <- filter( d, group %in% c("treatment","high.ses") &
                time %in% c("time0","time1") )

dm$treat.dummy <- ifelse( dm$group=="treatment", 1, 0 )
dm$pre.dummy   <- ifelse( dm$time=="time0", 1, 0 )
dm$post.dummy  <- ifelse( dm$time=="time1", 1, 0 )

dm$treat.post <- dm$treat.dummy * dm$post.dummy 

m <- lm( ability ~ treat.dummy + post.dummy + treat.post,
         data=dm)

stargazer( m, type=s.type,
           omit.stat=c("f","ser","adj.rsq"),
           intercept.top=TRUE, intercept.bottom=FALSE,
           digits=2 )
```

Does the control group model secular trend appropriately?

Test: are the study group trend lines parallel prior to the intervention? 

```{r, fig.height=8, fig.width=8, echo=F}
x <- c(42,50,64,78)
plot.new()
plot.window( xlim=c(40,90), ylim=c(-2,5) )

segments( x0=50, y0=-2, y1=3.1, lty=3, col="gray50" )
text( 50, 3.5, "Start of \nTreatment \nPeriods", col="gray50" )

abline( h=-2:5, lty=2, col="gray70" )
abline( h=-1.5:4.5, lty=3, col="gray85" )
points( x, t[1,], type="b", pch=19, col="steelblue", cex=3 )
points( x, t[2,], type="b", pch=19, col="steelblue", cex=3 )
points( x, t[3,], type="b", pch=19, col="steelblue", cex=3 )
axis( side=1 )
axis( side=2, las=1 )
title( xlab="Age (months)", ylab="Ability (logits)" )

text( 80, t[1,4], "High SES", col="steelblue", pos=4 )
text( 80, t[2,4], "Treatment Group", col="steelblue", pos=4 )
text( 80, t[3,4], "Control Group", col="steelblue", pos=4 )

# segments( x0=50, y0=t[1,2], y1=3.2, lty=3, col="steelblue" )
# text( 50, 3.7, "Start of \nTreatment \nPeriods", col="steelblue" )

points( x[1:2], t[2,1:2], 
        type="b", pch=19, col="darkred", cex=4, lwd=2, lty=3 )
points( x[1:2], t[3,1:2], 
        type="b", pch=19, col="darkred", cex=4, lwd=2, lty=3 )
```

```{r, results="asis"}
dm <- filter( d, group %in% c("treatment","control") &
                time %in% c("time0","time1") )

dm$treat.dummy <- ifelse( dm$group=="treatment", 1, 0 )
dm$post.dummy  <- ifelse( dm$time=="time1", 1, 0 )

dm$treat.post <- dm$treat.dummy * dm$post.dummy 

m <- lm( ability ~ treat.dummy + post.dummy + treat.post,
         data=dm)

stargazer( m, type=s.type,
           omit.stat=c("f","ser","adj.rsq"),
           intercept.top=TRUE, intercept.bottom=FALSE,
           digits=2 )
```

## Question 7

> Which coefficient captures the parallel lines test? Do we want it to be significant or not?

**ANSWER:**


























<br>

-----

<br>



```{css, echo=F}
body{
     font-family:system-ui,-apple-system,"Segoe UI",Roboto,Helvetica,Arial,sans-serif;
     font-size:calc(1.5em + 0.25vw);
     font-weight:300;line-height:1.65;
     -webkit-font-smoothing:antialiased;
     -moz-osx-font-smoothing:grayscale;
     margin-left:20%;
     margin-right:20%} 
     

h1, h2, h3, h4 { color: #995c00; }

h2 { margin-top:80px; }



.footer {
  background-color:#726e6e;
  height:340px;
  color:white;
  padding: 20px 3px 20px 3px;
  margin:0px;
  line-height: normal;
}

.footer a{ color:orange; text-decoration:bold !important; } 
 
 
 
 table{
   border-spacing:1px;
   margin-top:80px;
   margin-bottom:100px !important;
   margin-left: auto;
   margin-right: auto;
   align:center} 


 
td{ padding: 6px 10px 6px 10px } 

th{ text-align: left; } 

blockquote {
    margin-top:80px;
    border-left: 5px solid #995c00;
    color: #995c00;
}
```