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

```{r message = F, warning = F}
library(tidyverse)
library(pander)
library(nlme)
library(janitor)
```


# Problem 1

```{r}
data_mat <- matrix(c(
  1, 300, 0, 0, 4, 6,
  2, 300, 1, 0, 4, 6,
  3, 500, 0, 1, 4, 6,
  4, 500, 1, 1, 4, 6,
  5, 200, 0, 0, 10, 12,
  6, 200, 1, 0, 10, 12,
  7, 200, 0, 1, 10, 12,
  8, 200, 1, 1, 10, 12), byrow = TRUE, nrow = 8)

colnames(data_mat) <- c("Category", "num", "x",
                        "T", "y_0", "y_1")
df <- data.frame(data_mat)
```

## (a)


```{r}
N <- sum(df$num)
ace <- 1/N*sum(df$num * (df$y_1 - df$y_0))
ace
```

We can actually directly see that for each individual, the difference between $y^1$ and $y^0$ is 2 (which means the average is also this value), but the above equation also shows us the the ACE is `r ace`.

## (b)

If we only view $x$ as the only pre-treatment predictor (as suggested in the beginning of the problem), then we can observe that there $T$ and $x$ appear to be independent - aka $P(T=t,X=x) = P(T=t)P(X=x)$, this is explored specifically in part (d), but the below tables show this through showing that $P(T=t) = P(T=t|x = x)$ for $x \in \{0,1\}$. 


```{r warning = F}
tab_info <- df %>% 
  mutate(T = paste0("T",T)) %>%
  group_by(x,T) %>% summarise(count = sum(num)) %>%
  tidyr::pivot_wider(id_cols = x, names_from = T, values_from = count) %>%
  janitor::adorn_totals(where = c("row", "col")) 

tab_info

tab_info %>%
  janitor::adorn_percentages()
```
With our full fiew of the data we can see that $(y^0, y^1)$ pairing actually differ between the first four categories and the last four categories. Given that $x$ doesn't capture this different it wouldn't be too crazy to suggest that there are actually more pre-treatment predictors of the outcome. If one just defines another grouping variable `x2` to be $(1,1,1,1,0,0,0,0)$ to capture this we can see that `x2` and `T` are not independent of each other. None-the-less this goes beyond the problem description.

## (c)

```{r}
N_T <- sum(df$num[df$T == 1])
N_not_T <- sum(df$num[df$T == 0])

1/N_T * sum((df$num * df$y_1)[df$T == 1]) - 
  1/N_not_T * sum((df$num * df$y_0)[df$T == 0])
```

If there was no other confounders we'd expect this value to approximate the ACE. Interestingly it does not, but that seems to be related to the second commentary we had in part (b), which showed that the proportion those in the treatment group where `x2 = 1` was larger than the proportion of those in the treatment group where `x2 = 0`, and this relateds to the smaller value we see for the treatment effect as those in the `x2 = 1` group had smaller outcomes under both treatments than the `x2 = 0` group.


## (d)

```{r}
prob_T <-sum(df$num * (df$T == 1)) / sum(df$num)

prob_T
```

```{r}
prob_T_given_x1 <- sum(df$num * (df$x == 1) * (df$T == 1)) / 
  sum(df$num * (df$x == 1))

prob_T_given_x1
```

We find that, in the above population, $P(T=1) =$ `r round(prob_T,4)`, which is the same as $P(T=1|x=1)$.

## (e)

Below are 2 approaches to estimate the desired coefficient and se. Note that the coefficient values for the intercept and $T$ directly related to parts of the calculation in part (c). The intercept being the `mean(y[T=0])` and the coefficient of $T$ being the difference between `mean(y[T=1]) - mean(y[T=0])`.

### Weighted regression approach

A weighted regression solution is faster that expanding out the full dataset, but requires a correct ion for the model's degrees of freedom. See this [blogpost](https://www.r-bloggers.com/2015/09/linear-models-with-weighted-observations/) for more comments.

```{r}
df <- df %>% 
  mutate(y = ifelse(T, y_1, y_0))

wls <- lm(y ~ T + x, data = df, weights = num)
wls$df.residual <- sum(df$num) - length(wls$coefficients) 
#^ correction for number of observations / degrees of freedom

summ_w <- summary(wls)
summ_w
```

```{r}
coef(summ_w)[2,]
```

### Individualized approach

Another way would just be to create an "individual" data frame (where each row represents one individual). Below is a messy way to create the desired data.

```{r}
df2 <- data.frame(x = unlist(sapply(1:nrow(df), 
                                    function(idx){rep(df$x[idx], 
                                                      each = df$num[idx])})),
                  T = unlist(sapply(1:nrow(df), 
                                    function(idx){rep(df$T[idx], 
                                                      each = df$num[idx])})),
                  y = unlist(sapply(1:nrow(df), 
                                    function(idx){rep(df$y[idx], 
                                                      each = df$num[idx])})))

ls_all_obs <- lm(y ~ T + x, data = df2)
summ <- summary(ls_all_obs)
summ
```



```{r}
coef(summ)[2,]
```

# Problem 2.

## (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.

## (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.


# Problem 3
## (a)
A summary of the model is below. Diagnostic plots are in Figure \ref{f:3a_diagnostic}, and the estimated autocorrelation between the residuals is plotted in Figure \ref{f:3a_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 3a}
carlsen <- read.table("CarlsenQ.txt", header=TRUE)
m1 <- lm(Sales ~ . -Case - Time, data=carlsen)
pander(summary(m1), caption = "Problem 3a: Summary of ordinary linear model.")
# summary(m1)
```

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

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

\pagebreak

## (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 3b}
m2 <- gls(Sales ~ . -Case - Time, data=carlsen, correlation=corAR1())
cat("Problem 3b: Summary of AR1 model:"); summary(m2)
```

## (c)

```{r 3c}
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.

## (d)
A summary of the model is below. Diagnostic plots are in Figure \ref{f:3d_diagnostic}, and the autocorrelation plot of the residuals is in Figure \ref{f:3d_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 3c. 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 3d}
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 = paste("Problem 2d: Summary of ordinary linear model",
                       "with transformed variables."))
```

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

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