---
title: "Homework 06 Solutions"
date: "2022-10-25"
output: 
  pdf_document:
    keep_tex: true
    fig_caption: true
  header-includes:
  - \usepackage{amsmath}
  - \usepackage{amsthm}
  - \usepackage{hyperref}
  - \usepackage{xcolor}
  - \DeclareMathOperator{\Cov}{Cov}
  - \DeclareMathOperator{\Corr}{Corr}
  - \DeclareMathOperator{\Var}{Var}
  - \DeclareMathOperator{\upmodels}{\perp \!\!\! \perp}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, include = TRUE, message = FALSE)
knitr::opts_chunk$set(fig.width=6, fig.height=5, fig.align="center") 
# library(knitr)
```


\def\red#1{\textcolor{red}{#1}}
\def\blue#1{\textcolor{blue}{#1}}
\let\bar\overline
\let\ul\underline

\frenchspacing


\begin{center}
\Large
36-617: Applied Linear Models\\
Fall 2022\\
Solutions
\end{center}


```{r}

## library(arm)  ## includes lme4, MASS, Matrix

library(ggplot2); theme_set(theme_bw())
library(gridExtra) ## to arrange ggplots...

## library(GGally) ## for ggpairs...


library(leaps) ## regsubsets(), summary(), coef()

## library(car)   ## subsets(), mmps(), vif(), etc.

## library(marginalmodelplots) ## because mmps doesn't work for glm's...
## library(DHARMa)  ## for better glm residual plots...

## library(tidyverse) ## at last...

## library(foreign)  ## for "read.dta" for problem 3

library(ISLR2)

library(splines)

library(boot)

library(mgcv)

set.seed(1000000)  ## repeatable "random" selection of CV subsets...

```


# Problem 1: ISLR \#9, p. 324.

This question uses the variables dis (the weighted mean of distances to five Boston employment centers) and nox (nitrogen oxides concentration in parts per 10 million) from the Boston data. We will treat dis as the predictor and nox as the response.

## Problem 1(a) 

Use the poly() function to fit a cubic polynomial regression to
predict nox using dis. Report the regression output, and plot
the resulting data and polynomial fits.

\em

```{r}
data(Boston)

str(Boston)

## help(Boston)   -- to get the definitions of the variables

lm.1 <- glm(nox ~ poly(dis,3), data=Boston)
summary(lm.1)

plot(nox ~ dis, data=Boston)
fits <- predict(lm.1)
ord <- order(Boston$dis)
lines(fits[ord] ~ dis[ord], data=Boston)
```

\rm

## Problem 1(b) 

Plot the polynomial fits for a range of different polynomial
degrees (say, from 1 to 10), and report the associated residual
sum of squares.

\em

```{r fig.height=9}
par(mfrow=c(4,3))
for (i in 1:10) {
  lm.fit <- glm(nox ~ poly(dis,i), data=Boston)
  rss <- sum(residuals(lm.fit)^2)
  plot(nox ~ dis, data=Boston, col="grey",
       main=paste0("Polynomial of degree ",i,"\nRSS = ",round(rss,2)))
  fits <- predict(lm.fit)
  ord <- order(Boston$dis)
  lines(fits[ord] ~ dis[ord], data=Boston)
}
```

\rm

## Problem 1(c) 

Perform cross-validation or another approach to select the optimal degree for the polynomial, and explain your results.

\em

```{r}
## From help(cv.glm, package=boot) we see that 
## (1) cv.glm uses mean squared prediction error as the cross-validation criterion
## (2) We want $delta[1] as the raw cross-validation mean squared prediction error

## First we try LOOCV...

cv.err <- NULL
for (i in 1:10) {
  lm.fit <- glm(nox ~ poly(dis,i), data=Boston)
  cv.err <- c(cv.err,cv.glm(Boston,lm.fit)$delta[1])  
}
cv.err
cat(paste("Optimal LOOCV polynomial degree =",which(cv.err==min(cv.err)),"\n"))

## Now let's try 10-fold cross-validation...

cv.err <- NULL
for (i in 1:10) {
  lm.fit <- glm(nox ~ poly(dis,i), data=Boston)
  cv.err <- c(cv.err,cv.glm(Boston,lm.fit, K=10)$delta[1])  
}
cv.err
cat(paste("Optimal 10-fold CV polynomial degree =",which(cv.err==min(cv.err)),"\n"))
```
Leave-one-out cross validation (LOOCV) does not depend on a random selection of folds; LOOCV uses all $n$ folds with training set of size $n-1$ and test set of size 1 (the one observation we leave out).  On the other hand $K$-fold cross-validation will depend somewhat on what random set of folds get selected.  

Here, LOOCV suggests that the optimal polynomial fit has degree 3, and 10-fold cross-validation suggests degree 4.  While LOOCV will always choose degree 3 for this data (it is not random), the degree chosen by 10-fold cross-validation may change, depending on the particular random folds that are chosen.  We will use degree 3 in the rest of this problem.

\rm

## Problem 1(d) 

Use the bs() function to fit a regression spline to predict nox
using dis. Report the output for the fit using four degrees of
freedom. How did you choose the knots? Plot the resulting fit.

\em 

```{r}

lm.2 <- glm(nox ~ bs(dis,df=4),data=Boston)
summary(lm.2)

plot(nox ~ dis, data=Boston)
fits <- predict(lm.2)
ord <- order(Boston$dis)
lines(fits[ord] ~ dis[ord], data=Boston)
```

According to help(bs,package=splines):

* The default degree of the spline is 3.

* When we specify "df=4", the bs() function chooses df $-$ degree $=$ $4-3=1$ internal knot, and also one knot each at the min and max of "dis", for a total of 3 knots.  The internal knot is at the median of "dis".

This means we are fitting one cubic between min(dis) and median(dis), and another cubic between median(dis) and max(dis).

\rm

## Problem 1(e) 

Now fit a regression spline for a range of degrees of freedom, and
plot the resulting fits and report the resulting RSS. Describe the
results obtained.

\em

```{r fig.height=9}
par(mfrow=c(4,3))
for (i in 1:12) {
  lm.fit <- glm(nox ~ bs(dis,df=i), data=Boston)
  rss <- sum(residuals(lm.fit)^2)
  plot(nox ~ dis, data=Boston, col="grey",
       main=paste0("Regression spline, df = ",i,"\nRSS = ",round(rss,2)))
  fits <- predict(lm.fit)
  ord <- order(Boston$dis)
  lines(fits[ord] ~ dis[ord], data=Boston)
}
```

From the warnings above and what we learned about the bs() function in part (d), when the formula "df $-$ degree" says we should use 0 knots or a negative number of knots, the bs() function just sticks a knot at the median again.  Effectively, when we specify a df that is too small, the bs() function just increases the df to the minimum acceptable, which is df=3 in our case.

\rm


## Problem 1(f) 

Perform cross-validation or another approach in order to select
the best degrees of freedom for a regression spline on this data.
Describe your results.

\em

```{r warning=FALSE}

## again, we can try LOOCV first...
cv.err <- NULL
for (i in 1:12) {
  lm.fit <- glm(nox ~ bs(dis,df=i), data=Boston)
  cv.err <- c(cv.err,cv.glm(Boston,lm.fit)$delta[1])  
}
cv.err
cat(paste("Optimal LOOCV df for bs() =",which(cv.err==min(cv.err)),"\n"))

## ...and then 10-fold cv...
cv.err <- NULL
for (i in 1:12) {
  lm.fit <- glm(nox ~ bs(dis,df=i), data=Boston)
  cv.err <- c(cv.err,cv.glm(Boston,lm.fit, K=10)$delta[1])  
}
cv.err
cat(paste("Optimal 10-fold CV df for bs() =",which(cv.err==min(cv.err)),"\n"))
```

The 10-fold cross-validation result is somewhat unstable here and depends on which random folds we get.  In different runs, I've gotten df=5, 6, 7, 8 and 10 as optimal.  On the other hand, the LOOCV result is stable and non-random (because it considers all $n=506$ folds of size $n-1$, leaving one observation out at a time as a "test dataset"), and the result is df=10.

As a final check I will compare AIC and BIC for the 3rd degree polynomial and several of the regression spline models:

```{r}
lm.poly3 <- lm.1
lm.bs3 <- lm.2
lm.bs5 <- glm(nox ~ bs(dis,df=5), data=Boston)
lm.bs6 <- glm(nox ~ bs(dis,df=6), data=Boston)
lm.bs7 <- glm(nox ~ bs(dis,df=7), data=Boston)
lm.bs8 <- glm(nox ~ bs(dis,df=8), data=Boston)
lm.bs10 <- glm(nox ~ bs(dis,df=10), data=Boston)

models <- list(lm.poly3,lm.bs3,lm.bs5,lm.bs6,lm.bs7,lm.bs8,lm.bs10)
names(models) <- strsplit("lm.poly3,lm.bs3,lm.bs5,lm.bs6,lm.bs7,lm.bs8,lm.bs10",
                          ",")[[1]]

comparison <- data.frame(
  logLik = sapply(models,logLik),
  df = sapply(models, function(x) attributes(logLik(x))$df),
  AIC = sapply(models,AIC),
  BIC = sapply(models,BIC)
)
rownames(comparison) <- names(models)

round(comparison,2)
```

We see that AIC likes the regression spline with df $=$ 10 best, and BIC likes the df $=$ 5 regression spline model best.  Referring back to the regression spline plots, the curves for df $=$ 5,6 or 7 seem preferable to the somewhat wiggly df $=$ 8 or 10 plots.  (This is my usual "interpretability" bias: smooth curves are generally easier to tell a story about than wiggly ones.  If your goal is prediction, however, you may indeed prefer the df $=$ 10 plot, which did better with AIC and LOOCV.)


\rm


# Problem 2: ISLR \#10, pp. 324--325. 

This question relates to the College data set.

## Problem 2(a) 

Split the data into a training set and a test set. Using out-of-state
tuition as the response and the other variables as the predictors,
perform forward stepwise selection on the training set in order
to identify a satisfactory model that uses just a subset of the
predictors.

\em

```{r}
data(College)
str(College)
## help(College)  -- to get the definitions of the variables
```
We have $n=777<1000$ which is a smallish sample size for cross-validation, so I will use a 50-50 split (see lecture 8 from class).  

```{r}
n <- dim(College)[1] ## should be n=777
train <- sample(1:n,size=floor(n/2))
train <- ifelse(1:n %in% train,T,F)
test <- !train
train.data <- College[train,]
test.data <- College[test,]

n.train <- dim(train.data)[1]
n.test <- dim(test.data)[1]


forward <- regsubsets(Outstate ~ ., data=train.data,method="forward",nvmax=18)

tmp <- summary(forward)
p.max <- dim(train.data)[2]
p <- 1:(p.max-1)

results <- data.frame(tmp$which,rss=tmp$rss,adjr2=tmp$adjr2,
                    bic=n.train*log(tmp$rss)+log(n.train)*(p+2),
                    aic=n.train*log(tmp$rss)+2*(p+2),
                    caic=n.train*log(tmp$rss) + 2*(p+2) +
                      2*(p+2)*(p+1)/(n.train-p-1))

## results

names(results)
names(results)[2] <- "Private"  

minimize <- function(res,col) {
  obj <- res[,col]
  k <- (1:length(obj))[obj==min(obj)]
  return(res[k,])
}


sel.BIC <- minimize(results,"bic")[2:p.max]
sel.AIC <- minimize(results,"aic")[2:p.max]

fla.BIC <- paste0("Outstate ~ ",paste(names(sel.BIC)[unlist(sel.BIC)],collapse=" + "))
fla.AIC <- paste0("Outstate ~ ",paste(names(sel.AIC)[unlist(sel.AIC)],collapse=" + "))

fwd.BIC <- lm(fla.BIC, data=train.data)
fwd.AIC <- lm(fla.AIC, data=train.data)
                  
```

Just for fun, I compare the best AIC and BIC models with a partial F test, using the anova() function.  

From the anova() output below we see that the AIC model adds the variables Top25perc, Books and GradRate to the BIC model.

Unsurprisingly, AIC and the F test prefer the more complex model for this split of the data.  However, I will work with the simpler BiC model in what follows.


```{r}
anova(fwd.BIC,fwd.AIC)
```

Note that the details of the results above depend strongly on the particular split of the data into training and test subsets.  A different split, with a different set.seed() seed, could produce different (and possibly non-nested) "best" AIC and BIC models; the result of the anova() test (when the models turn out to be nested) can also change.

\rm

## Problem 2(b) 

Fit a GAM on the training data, using out-of-state tuition as
the response and the features selected in the previous step as
the predictors. Plot the results, and explain your findings.

\em

Since "Private" is in the model, but "Private" is a 0/1 variable, it doesn't get a transformation.  But all the other variables are quantitative and so get a transformation
```{r}

## Note that I am using cubic smoothing splines (bs='cr' in the code below) 
## to be consistent with what we did in class, but you could equally well 
## use mgcv's default "thin plate" splines...

fla.gam <- paste0("Outstate ~ Private + s(",
                  paste(names(sel.BIC)[-1][unlist(sel.BIC)[-1]],
                  collapse=",bs='cr') + s("),",bs='cr')")

gam.1 <- gam(formula(fla.gam), data=train.data,method="GCV.Cp")

summary(gam.1)
```


```{r fig.height=9}
plot(gam.1,pages=1)
```
It appears that PhD and perc.alumni not need transformations, while Accept and Expend clearly do.  In particular it looks like 

* Out of state tuition increases with institutional expenditures per student, up to a point at which it levels off (perhaps this is the highest tuition that the market will bear)

* Out of state tuition also increases with the number of students accepted (which seems a bit suprising to me) but the influence of more students accepted, on out of state tuition, lessens after roughly Accept=4000.

\rm


## Problem 2(c) 

Evaluate the model obtained on the test set, and explain the
results obtained.

\em

First we do mean squared prediction error for the two models
```{r}
mse.BIC <- mean((test.data$Outstate - predict(fwd.BIC,newdata=test.data))^2)
mse.gam <- mean((test.data$Outstate - predict(gam.1,newdata=test.data))^2)

mse.results <- data.frame("Prediction MSE"=c(fwd.BIC=mse.BIC,gam.1=mse.gam))
round(mse.results)
```
Now we do a goodness-of-fit comparison:
```{r}

test.BIC <- lm(formula(fla.BIC),data=test.data)
test.gam <- gam(formula(fla.gam),data=test.data)
models <- list(test.BIC,test.gam)
names(models) <- strsplit("test.BIC,test.gam",",")[[1]]

comparison <- data.frame(
  logLik = sapply(models,logLik),
  df = sapply(models, function(x) attributes(logLik(x))$df),
  AIC = sapply(models,AIC),
  BIC = sapply(models,BIC)
)
rownames(comparison) <- names(models)

round(comparison,2)

```
All three comparisons, prediction MSE, AIC, and BIC, prefer the gam model in this example.  

Explanation: The extra flexibility of the nonparametric transformations of the predictors provide better predictions, reducing the mean squared prediction error by (`r mse.results[2,1]`)/(`r mse.results[1,1]`) = `r round(mse.results[2,1]/mse.results[1,1],2)`, as well as dramatic improvements in fit as reflected by reductions of 60--100 in AIC and BIC. 

\rm


## Problem 2(d) 

For which variables, if any, is there evidence of a non-linear
relationship with the response?


\em

From the transformation plots above, it seems clear that Accept and Expand have a non-linear relationship with the response.  There really isn't a clear case for any of the others, though one could make an argument for Enroll and Personal.

One could more formally assess whether a nonlinear term is needed for Personal, for example, by comparing the model with s(Personal) as a predictor, with the model that just uses Personal as a predictor.

Once again, these results may vary, depending on what split into training and test data sets you are working with.

\rm


# Problem 3 (Based on Gelman \& Hill, Chapter 9, \#4).  

The table below describes a hypothetical experiment on 2400 persons.  Each row of the table specifies a category of person, as defined by his or her pre-treatment predictor $x$, treatment indicator $T$, and potential outcomes $y^0$ and $y^1$.  (For simplicity, we assume---unrealistically---that all people in this experiment fit into one of these eight categories).
\begin{center}
\begin{tabular}{cccccc}
Category & \# persons in category & $x$ & $T$ & $y^0$ & $y^1$ \\
\hline
\hline
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 \\
\hline
\end{tabular}
\end{center}
In making this table we are assuming omniscience, so that we know both $y^0$ and $y^1$ for all observations.  But the (non-omniscient) investigator woud only observe $x$, $T$, and $y=y^T$ for each unit.  (For example, a person in category 1 would have $x=0$, $T=0$ and $y=4$, and a person in category 3 would have $x=0$, $T=1$, and $y=6$.)

\em

Let's start by putting the data above into R...
```{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)
```

\rm

## Problem 3(a)

What is the true ACE (average causal effect), if we could observe $y^1$ and $y^0$ for every person in this population of 2400 persons?

\em

```{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`.


\rm

## Problem 3(b)

Another population quantity is the mean of $y$ for those who received the treatment, minus the mean of $y$ for those who did not.  What is the relationship between this quantity and the quantity you calculated in part (a)?

\em

```{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!  Below in part (d) we explore the data further and discover that there is another confounding variable not accounted for in the data.



\rm


## Problem 3(c) 

Suppose we draw a person randomly from this population or 2400 people.  What is the probability that $T=1$ for this person?  If I tell you that $x=1$ for this person, does that change the probability that $T=1$?  I.e., is $P[T=1|x=1]=P[T=1]$?

\em

```{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]$.


\rm

## Problem 3(d)

Is it plausible to believe this data came from a randomized experiment?  Defend your answer.

\em

From part (c) we know that $P[T=1]=P[T=1|x=1]$, and more generally, $P[T=t,x=x] = P[T=t]P[x=x]$, i.e. $T$ and $X$ are independent.  So, if we think $x$ is the only pre-treatment covariate (as suggested in the beginning of the problem), it is plausible that this is data from a randomized experiemnt (since randomization forces $T$ to be indpendent of any pre-treatment covariates). 

However, looking back at the data table, it seems that there is a jump of 6 in the level of response (for both $y^0$ and $y^1$) going from the first four categories, vs. the last four categories.  For example, if $y$ is measuring a health outcome after a hospital stay, and the first four categories represent people who would go to hospital \#1 and the second four categories represent people who would go to hospital \#2, it appears that the first group of people are generally less healthy than the second group.

If we define another grouping variable 

```{r}
df$x2 <- c(1,1,1,1,0,0,0,0)
```

then (i) we can see that `x2` and `T` are not independent: `P[T=1|x2=1] = 1000/1600 = ` `r 1000/1600` $\neq$ `P[T=1|x2=0] = 400/800 = ` `r 400/800`, and (ii) in the regression `y ~ T + x2`, the coefficient $\hat\beta_T$ would plausibly be an unbiased estimate of ACE.  The simpler regression `y ~ T` (which is essentially what we were doing in part (b) above) would produce a biased $\hat\beta_T$. 
```{r}

```

\rm

## Problem 3(e)

For this hypothetical data, figure out (by hand, or by using R) the estimate and standard error of the coefficient of $T$ in a regression of $y$ on $T$ and $x$.  Is this a reasonable way to estimate the true ACE?

\em

Below are two 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 (b). 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])`.

### Splitting the data into individual cases

```{r}

new.df <- data.frame(0,0,0,0,0,0,0)  ## a dummy line so we can rbind to it
names(new.df) <- names(df)
for (i in 1:8) {
  for (j in 1:df[i,"num"]) {
    new.df <- rbind(new.df,df[i,])
  }
}
new.df <- new.df[-1,]  ## remove the dummy line

new.df$y <- ifelse(new.df$T==0,new.df$y_0,new.df$y_1)

summary(lm.splitting <- lm(y ~ T + x,data=new.df))
```

### Weighted regression approach

A weighted regression solution is faster than expanding out the full dataset, 
```{r}
df$y <- ifelse(df$T==0,df$y_0,df$y_1)

lm.weighted <- lm(y ~ T + x, data = df, weights = num)
summary(lm.weighted)  ## uncorrected summary
```
but requires a correction for the model's residual degrees of freedom. See this [blogpost](https://www.r-bloggers.com/2015/09/linear-models-with-weighted-observations/) for more comments.


The correct residual degrees of freedom are
```{r}
df.res <- sum(df$num) - length(lm.weighted$coefficients)
df.res
```
and then new $s^2$ and $SE(\hat\beta_T)$ can be calculated with this new residual df.

For now we will just observe that the point estimates are the same:
Note that we get the same estimate either way:
```{r}
round(rbind("b_T from lm.splitting"=summary(lm.splitting)$coef[2,1],
      "b_T from lm.weighted"=summary(lm.weighted)$coef[2,1]),2)
```

This is *not* a reasonable way to estimate the ACE: Since we are not accounting for the variable `x2`, we get the same biased estimate as in part (b).

If we include `x2` also, we can get a better answer:
```{r}
corrected.lm.splitting <- lm(y ~ T + x + x2,data=new.df)
summary(corrected.lm.splitting)
```
We can see that the coefficient $\hat\beta_T$ is now the correct value, 2, to many decimal places (and with a very tight confidence interval!).

\rm

# Problem 4 (Based on Gelman \& Hill, Chapter 9, \#6).  

You are consulting for a researcher who has performed a randomized trial where the treatment was a series of 26 weekly therapy sessions, the control was no therapy, and the outcome was self-report of emotional state one year later. However, most people in the treatment group did not attend every therapy session. In fact there was a good deal of variation in the number of therapy sessions actually attended. The researcher is concerned that her results represent ``watered down'' estimates because of this variation and suggests adding in another predictor to the model: number of therapy sessions attended. 

## Problem 4(a)

What would you advise her about her suggestion?  Carefully justify your answer in terms of our discussion about controlled experiments, observational studies, confounders, and so forth.

\em

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.


\rm

## Problem 4(b)

Can you suggest another kind of analysis that might give her a better estimate of the treatment effect?  Carefully explain why this might work.  \em (Hint: Think about the more sophisticated analyses presented in lecture 16 in class.)

\rm

\em

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.



\rm



