---
title: "Pygmalion"
output: html_document
---

```{r setup, include=TRUE, message=FALSE}
knitr::opts_chunk$set(echo = TRUE)
require(Sleuth3)
require(ggplot2)
require(tigerstats)
require(xtable)
require(lattice)

pyg = case1302
```

First, let's visualize the data. We need to use `aggregate` since we have two control observations for each company

```{r pygvis}
pyg_mean <- aggregate(pyg$Score, by=list(pyg$Company, pyg$Treat), mean)
#pyg_mean <- pyg_mean[order(pyg_mean$Company),]
colnames(pyg_mean) <- colnames(pyg)
barchart(Score~Company, groups=Treat, data=pyg_mean, auto.key=list(columns=2))
```

This presents an exaggerated picture of what's actually going on... Let's make the origin be 0 in that barchart

```{r}
barchart(Score~Company, groups=Treat, data=pyg_mean, auto.key=list(columns=2), origin=0)
```


Another way to visualize the data is using interaction plots.

```{r}
with(pyg, interaction.plot(x.factor=Company, trace.factor=Treat, Score, col=c("red", "blue")))
with(pyg, interaction.plot(x.factor=Treat, trace.factor=Company, Score, col=c("red", "blue")))
```

It looks like for most companies, there is a roughly constant Pygmalion effect.

Let's now run the regression for model with the interactions

```{r}
fit_saturated <- lm(Score~Treat+Company+Treat*Company, data=pyg)
summary(fit_saturated)
```


Let's look at the residual plot:
```{r}
plot(fit_saturated,1)
```

It looks like the variance is not really constant, with larger residuals for smaller scores. (Why would that be?)

It's hard to come up with a data transformation that would fix that, however. What would a log transformation working imply?


Let's now perform a partial F-test:
```{r}
anova(fit_saturated)
```


```{r}
fit_reduced <- lm(Score~Treat+Company, data=pyg)
summary(fit_reduced)
anova(fit_reduced)
plot(fit_reduced, 1)
```

Note that the estimate for the effect is quite different from the full fit. Deciding which estimate to use has to do with whether we have a prior reason to believe in there being interactions. In this case, we probably don't, and there was no evidence in the data that there are interactions, so it's reasonable to leave them out. 

Visualize the Pygmalion effect using boxplots:

```{r}
fit1 <- lm(Score ~ Company, data= pyg)
pyg$resids <- resid(fit1)
ggplot(pyg, aes(x= Treat, y= resids)) + geom_boxplot()
```

