---
title: "Homework 07 Solutions"
date: "11/09/2018"
output: 
  pdf_document:
    keep_tex: true
fig_caption: true
header-includes:
  - \usepackage{amsmath}
  - \usepackage{amsthm}
  - \DeclareMathOperator{\Cov}{Cov}
  - \DeclareMathOperator{\Corr}{Corr}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(cache = FALSE, echo = FALSE, include = TRUE, message = FALSE)
options(scipen = 1, digits = 2)

library(arm)
library(knitr)
library(nlme)
library(pander)
```

# 1.
## (a). 
Summaries and diagnostic plots for both models are below. Since the variance of a Poisson random variable is equal to the mean, the raw residual plot and scale-location plot are not very informative, so we will focus on the other two diagnostic plots.

For model1, the model with no interaction term, the residuals appear to be approximately normally distributed (Figure \ref{f:1a_diagnostic1}). All the points have relatively high leverage, and points 5, 9, and 24 have been labeled as influential points that be worth investigating further. Many of these points have standardized Pearson residuals with absolute values between 2 and 4, so there are appear to be a lot of outliers. The binned residual plot doesn't give any particular cause for concern (Figure \ref{f:1a_binned1}).

For model2, the model with the interaction term, the residuals also appear to be 
approximately normally distributed (Figure \ref{f:1a_diagnostic2}). Leverage is not shown directly on the leverage plot in this case, since R has chosen to plot the factors of the `Wool` variable on the x-axis instead, but there are no points that are marked as influential. There are still a lot of points that look like outliers, with very large standardized Pearson residual values. The binned residual plot now has values very close to 0, though all the values are slightly negative.

**Summary of Model 1 with no interaction:**

```{r 1a}
data(warpbreaks)
wool1 <- glm(breaks ~ wool + tension, data = warpbreaks, family = poisson)
wool2 <- glm(breaks ~ wool * tension, data = warpbreaks, family = poisson)
pander(summary(wool1), caption = "Problem 1(a): Summary of wool model without interaction term.")
```

```{r 1a_plot1, fig.height=4, fig.cap="Problem 1(a): Diagnostic plots for model without interaction term.\\label{f:1a_diagnostic1}"}
par(mfrow = c(2, 2),  mar = c(4, 4.5, 2, 1))
plot(wool1)
```

```{r 1a_binnedplot1, fig.cap="Problem 1(a): Binned plot of residuals vs. fitted values for model without interaction term.\\label{f:1a_binned1}"}
binnedplot(predict(wool1), resid(wool1))
```

**Summary of Model 2 with interaction term:**

```{r 1a_summary2}
pander(summary(wool2), caption = "Problem 1(a): Summary of wool model with interaction term.")
```

```{r 1a_plot2, fig.cap="Problem 1(a): Diagnostic plots for model with interaction term.\\label{f:1a_diagnostic2}"}
par(mfrow = c(2, 2),  mar = c(4, 5, 2, 1))
plot(wool2)
```

```{r 1a_binnedplot2, fig.cap="Problem 1(a): Binned plot of residuals vs. fitted values for model with interaction term.\\label{1a_binned2}"}
binnedplot(predict(wool2), resid(wool2))
```

## (b) i.
The AICs and BICs are as follows:

```{r, include=TRUE}
ICs <- data.frame(AIC = c(AIC(wool1), AIC(wool2)),
                  BIC = c(BIC(wool1), BIC(wool2)))
row.names(ICs) <- c("Model 1 (without interaction)",
                    "Model 2 (with interaction)")
kable(ICs)
```

Both metrics are lower for model 2, so that is the model we would choose.

## (b) ii.

```{r lrt}
T <- -2*(as.numeric(logLik(wool1)) - as.numeric(logLik(wool2)))
pval <- pchisq(T, 2, lower.tail = FALSE)
```

The log likelihoods are -243 for model 1 without the interaction and -228 for model 2 with the interaction. The difference in degrees of freedom is 2. This comes from the fact that the variable `wool` has two levels, which corresponds to a single dummy variable, and the variable `tension` has three levels, which corresponds to two dummy variables. The interaction between them therefore corresponds to `2 x 1 = 2` dummy variables.

We therefore compare the statistic $-2(-243-(-228)) = 28$ to a $\chi_{(2)}^2$ distribution. The pvalue is $`r pval`$, so we reject the null hypothesis and decide to include the interaction. The likelihood ratio test therefore selects the same model as the AIC and BIC.

## (c)

```{r 1c}
y <- warpbreaks$breaks
y_hat <- predict(wool2, type = "response")
r_pearson <- (y - y_hat)/sqrt(y_hat)
T = sum(r_pearson^2)
deg_freedom <- length(y) - length(coef(wool2))
pval <- pchisq(T, deg_freedom, lower.tail = FALSE)
```

The Pearson residuals are
\[
    r_i = \frac{y_i - \hat{y}_i}{\sqrt{\hat{y}_i}}
\]
If the model is correctly specified, then $\sum_{i=1}^n r_i$ should have approximately a $\chi_{n-(p+1)^2}$ distribution.

We have $\sum_{i=1}^n r_i = `r T`$. Comparing this to a $\chi^2$ distribution with $`r deg_freedom`$ degrees of freedom yields a p-value of $`r pval`$. We conclude therefore that the residuals are overdispersed, meaning that the variability of the responses is larger than what we would expect if they really followed a conditional Poisson distribution.

## (d)
**Summary of Model 3, the quasipoisson model:**

```{r 1d}
wool3 <- glm(breaks ~ wool * tension, data = warpbreaks, family = quasipoisson)
pander(summary(wool3), caption = "Problem 1(d): Summary of quasipoisson wool model with interaction term.")
```

The coefficient estimates are identical to those in Model 2, but the estimated standard errors are much larger, as expected. Consequently, the p-values are larger, though still significant at customary thresholds like 0.05 or 0.01. The dispersion parameter is estimated to be 3.8, whereas it was assumed to be 1 in Model 2.

# 2(a)
A summary of the model is below. Diagnostic plots are in Figure \ref{f:2a_diagnostic}, and the estimated autocorrelation between the residuals is plotted in Figure \ref{f:2a_acf}.

The residuals are very roughly normally distributed, but they appear to show an increasing trend relative to the fitted values. The scale-location plot also appears to show that the variance is increasing with fitted values, and the autocorrelation plot shows large autocorrelations, all of which suggest that the linear model with assumed independence between the residuals is not appropriate.

```{r 2a}
carlsen <- read.table("CarlsenQ.txt", header=TRUE)
m1 <- lm(Sales ~ . -Case - Time, data=carlsen)
pander(summary(m1), caption = "Problem 2a: Summary of ordinary linear model.")
# summary(m1)
```

```{r 2a_diagnostic, fig.cap="Problem 2a: Diagnostic plots for ordinary linear model.\\label{f:2a_diagnostic}", fig.height=4}
par(mfrow = c(2, 2),  mar = c(4, 5, 2, 1))
plot(m1)
```

```{r 2a_acf, fig.height=4, fig.align="center", fig.cap="Problem 2a: acf plot for ordinary linear model residuals.\\label{f:2a_acf}"}
acf(resid(m1), main = "")
```

\pagebreak

# 2(b)
Note that the model can be estimated with either Maximum Likelihood (ML) or Restricted Maximum Likelihood (REML). I have used the default, which is REML, but either one is fine. They will produce different estimates. The code should look like one of the following: 

```{r, echo=TRUE, eval=FALSE}
model <- gls(Sales ~ . -Case - Time, data=carlsen, correlation=corAR1())
model <- gls(Sales ~ . -Case - Time, data=carlsen, correlation=corAR1(),
             method = "ML")
```

i. A summary of this model is below. From this summary, the estimated value of $\rho$, the lag-1 autocorrelation parameter, is 0.87. (If you use ML for estimation, the estimate is 0.79).

ii. The coefficient values have changed in magnitude, though not in sign. The estimated standard errors in the AR1 model are much smaller than in the ordinary linear model, which is to be expected if the model is properly accounting for correlation among the residuals. As a result, both `Sun` and `Q4` are significant in the AR1 model but not in the ordinary linear model at conventional thresholds.

```{r 2b}
m2 <- gls(Sales ~ . -Case - Time, data=carlsen, correlation=corAR1())
cat("Problem 2b: Summary of AR1 model:")
summary(m2)
```

# 2(c)

```{r 2c}
T <- -2*(as.numeric(logLik(m1)) - as.numeric(logLik(m2)))
pval <- pchisq(T, df = 1, lower.tail = FALSE)
```

The log likelihoods are -131 for the ordinary linear model and -98 for the AR1 mode. The difference in degrees of freedom is 1, so we compare the statistic $-2(-131-(-98)) = 66$ to a $\chi_{(1)}^2$ distribution. The pvalue is $`r pval`$, so we reject the null hypothesis and conclude that the autocorrelation is not 0.

# 2(d)
A summary of the model is below. Diagnostic plots are in Figure \ref{f:2d_diagnostic}, and the autocorrelation plot of the residuals is in Figure \ref{f:2d_acf}.

The diagnostic plots suggest problems with the fit. The residuals are not very close to normally distributed, and the scale-location plot appears to show a relationship between the fitted values and the variance. The autocorrelation plot, however, shows much reduced autocorrelation among the residuals, as hoped for.

The likelihood of the model is `-130`, as compared to `-98` for the model in 2c. It appears to be more effective to estimate the autocorrelation simultaneously with estimating the model, as in 2c, rather than following a two-stage estimation procedure.

```{r 2d}
r <- resid(m1)
rho <- cor(r[-1], r[-length(r)])
Sigma <- diag(nrow(carlsen))
Sigma <- rho^abs(row(Sigma) - col(Sigma))
S <- chol(Sigma)
S_inv <- solve(t(S))

Xstar <- S_inv %*% model.matrix(m1)
ystar <- S_inv %*% carlsen$Sales

m3 <- lm(ystar ~ Xstar - 1)
pander(summary(m3), caption = "Problem 2d: Summary of ordinary linear model with transformed variables.")
```

```{r 2d_diagnostics, fig.cap="Problem 2d: diagnostic plots.\\label{f:2d_diagnostic}"}
par(mfrow = c(2, 2),  mar = c(4, 4.5, 2, 1))
plot(m3)
```

```{r 2d_acf, fig.cap="Problem 2d: acf plot for ordinary linear model residuals.\\label{f:2d_acf}"}
acf(resid(m3), main = "")
```

# 3.
(a) The treatment variable is no longer randomized in the way intended by the researcher, because individuals have selected their own level of treatment. This selection process is very likely confounded with the outcome variable, meaning that there is some unobserved variable that is correlated with both the treatment variable and the outcome.

For example, imagine that the therapy actually has no effect, but there is an "optimism" variable that determines both how many sessions a person attends and their emotional state at the end of the study. In other words, the people who engage the most with the treatment are those in the best state at the end of the study, but not because of the therapy. Then the number of sessions attended will strongly predict the outcome variable, but not because of an effect of the therapy.

An unobserved confounder could also obscure a real treatment effect. Suppose that instead of an optimism variable, there is a variable representing emotional wellbeing that is inversely correlated with the number of sessions attended. That is, those who start out the most depressed seek the most treatment, perhaps because they feel that they have the most to gain. Suppose the treatment *does* have a beneficial effect that increases with the number of sessions attended, and suppose that everyone seeks out just enough treatment to get them to roughly the same level of emotional wellbeing at the end of the study. Then a coefficient for number of treatment sessions could end up being 0 even though there is a treatment effect!

*In sum:* whenever there is noncompliance in an experiment, there is potential for confounding due to selection effects. Including a variable for the amount of treatment received *does not* address this and can lead to highly misleading results.

# 3(b)
This is perfectly set up for an Instrumental Variable analysis. The assignment variable (assignment to treatment or control) is the instrument, and the number of sessions attended is the treatment variable. The assignment variable almost certainly affects the treatment variable, but it should have no effect on the outcome through other pathways. (That is, the mere fact of having been assigned to the treatment or control group shouldn't affect participants' emotional wellbeing at the end of the study.) An experimental design with noncompliance is almost always a good candidate for an IV analysis.

An analysis based on propensity score matching won't work here unless we have variables with which to construct the propensity score. The propensity score is the probability of receiving treatment as a function of some set of covariates $X$. That is, the propensity score is a function $\pi(x) := P(A = 1|X = x)$, where $A = 1$ means that someone receives treatment. We can define it more broadly as the probability of receiving a particular amount of treatment $a$: $\pi_a(x) = P(A = a|X= x)$. We want as rich a set of covariates $X$ as possible in order to perform matching. If the investigator collected these prior to the start of the study, then this analysis may be possible.


