---
title: " Homework 2 Solution"
author: "STA2101F 2020"
date: ''
output:
  pdf_document: default
  html_document: default
fontsize: 12pt
---
```{r setup, include=FALSE}
# knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(
  echo = FALSE,
  message = FALSE,
  warning = FALSE,
  error = FALSE, 
  out.width = "100%",
  fig.width = 10,
  fig.height = 6.7, 
  fig.retina = 3,
  cache = FALSE)
```
```{r library}
library(faraway)
library(plotrix)
library(kableExtra)
library(xtable)
library(multcomp)
library(caret)
data(butterfat)
data(prostate)
```
## Due November 5 2020 11.59 pm

**Homework to be submitted through Quercus**

Please submit both the .Rmd and the .pdf file for this homework. (It is not possible to annotate .html files in Quercus.)  In the .pdf document you should just include properly formatted output, not unedited chunks of ##R output. See the solutions to HW 1 Q 4, for example.

You are welcome to discuss the questions with others, but the solutions and code must be written independently. 

\bigskip

1. *Two-way layout with replication*.  If we have a balanced two-factor experiment with replication, a suitable linear model might be 
   \begin{equation} y_{ijk} = \mu + \alpha_i + \beta_j + \gamma_{ij} + \epsilon_{ijk},
   \end{equation}
   where $i = 1, \dots, I; j = 1, \dots, J; k = 1, \dots, K$. In this model we separate the main effects of factors $A$ and $B$ from the interaction effect. If the observations were set out in a two-way layout, there would be $K$ observations in each cell.   
    (a) If we assume that $\epsilon_{ijk}$ are all independent and identically distributed $N(0, \sigma^2)$, argue that the variation in a single cell, $\Sigma_k (y_{ijk} - \bar y_{ij.})^2$ estimates $\sigma^2$ with $K-1$ degrees of freedom. From this conclude that $\Sigma_{ij} \Sigma_k (y_{ijk} - \bar y_{ij.})^2$ estimates $\sigma^2$ with $IJ(K-1)$  degrees of freedom. Suggest and construct a plot that might be used to check on the assumption of constant variance. \newline
    \textcolor{blue} {\textbf{Solution:} Since $\epsilon_{ijk} \overset{\text{iid}}{\sim} N(0, \sigma^2)$, for fixed $i$ and $j$:
\[ \frac{\sum_k (y_{ijk} - \bar{y}_{ij \cdot})^2}{\sigma^2} \sim \chi^2_{(K-1)} \quad \text{(using Cochran's theorem)}\]
For $i' \neq i$ and $j' \neq j$, the random variables $\frac{\sum_k (y_{ijk} - \bar{y}_{ij \cdot})^2}{\sigma^2}$ and $\frac{\sum_k (y_{i'j'k} - \bar{y}_{i'j' \cdot})^2}{\sigma^2}$ are independent. The sum of these random variables are the sum of two independent $\chi^2$ random variables with $K-1$ df each, which yields a $\chi^2$ with $2(K - 1)$ df. Repeating this by summing over $i$ and $j$:
\[ \sum_{ij} \frac{\sum_k (y_{ijk} - \bar{y}_{ij \cdot})^2}{\sigma^2} \sim \chi^2_{IJ(K-1)} \]
In Figure 1, the log of sample variance for each cell is plotted against the cell mean. A positive pattern can clearly be observed, which indicates a violation of constant variance.}
```{r, echo=FALSE, fig.cap = "Confidence intervals for the variance estimates"}
Sij <- aggregate(Butterfat ~ Breed + Age, butterfat, var)[,3]
ybarij <- aggregate(Butterfat ~ Breed + Age, butterfat, mean)[,3]

plot(ybarij, log(Sij), xlab="Mean", ylab="log sample variance")
```
    (b) The data set `butterfat` in `library(faraway)` is a two-way layout with $10$ observations per cell. Present a two-way table of cell means, showing as well the row and column means. Is there evidence of interaction between the two factors? Explain. \newline
    \textcolor{blue} {\textbf{Solution:} The following table shows the cell means of the butterfat content for the Breed-Age combinations as well as the row and column means:}
    ```{r, results='asis', echo=FALSE}
cellmean <- aggregate(Butterfat ~ Breed + Age, butterfat, mean)
rowmean <- aggregate(Butterfat ~ Breed, butterfat, mean)
colmean <- aggregate(Butterfat ~ Age, butterfat, mean)
table <- data.frame('2Years' = cellmean$Butterfat[cellmean$Age=='2year'], 
                    Mature = cellmean$Butterfat[cellmean$Age=='Mature'],
                    'Row Mean' = rowmean$Butterfat)
table <- rbind(table, c(colmean$Butterfat, mean(butterfat$Butterfat)))
rownames(table) <- c(as.character(rowmean$Breed), 'Column Means')
colnames(table) <- c('Two Years', 'Mature', 'Row Means')
kable(table, 'latex', digits = 3, booktabs = T, align = "c")%>%
  kable_styling(latex_options = "HOLD_position")
```
    \textcolor{blue} {The cow's age appears to have a positive effect on the average butterfat content for all breeds except the Canadian and Holstein-Fresian cows whose average butterfat content was almost equal or lower. These differential age effects suggests that the breed and age factors potentially interact in their effect on butterfat content of the milk.}
    (c) Fit the linear model of equation (1) to this data and give the analysis of variance table. Give an algebraic expression for the sum of squares due to interaction in the analysis of variance table. What is the estimate of $\sigma^2$? Does it have $IJ(K-1)$ degrees of freedom?^[if not you've made a mistake somewhere.] \newline
    \textcolor{blue} {\textbf{Solution:} Model (1) was fitted to the butterfat data. The resulting ANOVA table is show by the following:}
```{r, echo=FALSE, results='asis'}
model1 <- anova(lm(Butterfat ~ factor(Breed) + factor(Age) + factor(Breed)*factor(Age), data = butterfat))
model1[ ,3:5] <- round(model1[ ,3:5],3)
model1[4,4:5] <- ''
kable(model1, 'latex', digits = 3, booktabs = T, align = "c")%>%
  kable_styling(latex_options = "HOLD_position")
```
    \textcolor{blue} {The total sum of squares can be partitioned as follows:
    \begin{align*}
    \sum_{ijk}(y_{ijk} - \bar{y})^2 &= \sum_{ijk}\left((y_{ijk} - \bar{y}_{ij.}) + (\bar{y}_{ij.} - \bar{y}_{i..} - \bar{y}_{.j.} + \bar{y}) + (\bar{y}_{i..} - \bar{y}) + (\bar{y}_{.j.}-\bar{y})\right)^2 \\
    &= \underbrace{\sum_{ijk}(y_{ijk} - \bar{y}_{ij.})^2}_{\text{Residuals}} + \underbrace{K\sum_{ij}(\bar{y}_{ij.} - \bar{y}_{i..} - \bar{y}_{.j.} + \bar{y})^2}_{\text{Interaction}} + \\
    & \underbrace{KJ\sum_{i}(\bar{y}_{i..} - \bar{y})^2}_{\text{Main effect of }\alpha} + \underbrace{KI\sum_{j}(\bar{y}_{.j.} - \bar{y})^2}_{\text{Main effect of }\beta}.
    \end{align*}
    The summary output in R shows that the estimated variance is 0.1731 with 90 df. We can confirm that the df is the same as in (a): $IJ(K-1) = 5 \times 2 \times (10 - 1) = 90$.}
    (d) Provide at most three plots of the data and/or fitted model. What information is available from these plots? \newline
    \textcolor{blue} {\textbf{Solution:} Figure 2 - a) shows the standardized residuals plotted against the fitted values. The spread of the standardized residuals appears to increase over the fitted values. Figure 2 - b) shows the QQ-plot for the standardized residuals. The plots lie approximately on the diagonal line shown in red, which indicates that the normality assumption appears to be appropriate. Figure 2 - c) shows lines joining the average butterfat contents of two-year-old cows and those of mature cows for each breed. Since the lines are not parallel, there is a possibility of a Breed-Age interaction. }
```{r, echo=FALSE, fig.cap="Plots of the data and the fitted model"}
cellmean <- cellmean[order(cellmean$Age, cellmean$Breed), ]
cellmean$mature <- (cellmean$Age == 'Mature')*1
bread_kind <- unique(cellmean$Breed)

model1 <- lm(Butterfat ~ factor(Breed) + factor(Age) + factor(Breed)*factor(Age), data = butterfat)

par(mfrow=c(2,2))
plot(model1$fitted.values, scale(model1$residuals), main="a) Standardised residuals by fitted values", xlab="Fitted values", ylab="Standardized residuals")

qqnorm(scale(model1$residuals), main="b) Normal QQ-plot")
abline(coef=c(0, 1), col='red')

plot(cellmean$mature[cellmean$Breed == bread_kind[1]], cellmean$Butterfat[cellmean$Breed == bread_kind[1]], xaxt="n", xlab="Age", ylab="Butterfat content", xlim=c(-0.5, 1.5), ylim=c(min(cellmean$Butterfat)-0.5, max(cellmean$Butterfat+0.5)), type='b',
     main="c) Age and Breed effects on the butterfat content")
axis(1, c(0, 1), c("2 Years", "Mature"))
for(i in 2:length(bread_kind))
{
  lines(cellmean$mature[cellmean$Breed == bread_kind[i]], cellmean$Butterfat[cellmean$Breed == bread_kind[i]], type='b')
}
```
    (e) Is the best breed in terms of butterfat content clearly superior to the second best breed?  Why or why not? Explain your calculations. \newline
    \textcolor{blue} {\textbf{Solution:} Since the data showed showed a weak evidence for the interaction effects, we ignore the interaction effects for this question. The following table shows the sample means $\bar{y}_{i..}=\frac{\sum_j \sum_k y_{ijk}}{JK}$:}
```{r, echo=FALSE, results='asis'}
ybari <- aggregate(Butterfat ~ Breed, butterfat, mean)
ybari <- cbind(i=1:5, ybari)
kable(ybari, 'latex', digits = 3, booktabs = T, align = "c")%>%
  kable_styling(latex_options = "HOLD_position")   
```
    \textcolor{blue} {The two largest mean butterfat contents came from Jersey $(i = 5)$ and Guernsey $(i = 3)$ breeds. The null hypothesis is $H_0: \alpha_5 - \alpha_3 = 0$. Assuming that age and breed don't interact, an unbiased estimator of this difference is $\bar{Y}_{5..} - \bar{Y}_{3..}$. The variance for the difference is $Var(\bar{Y}_{5..} - \bar{Y}_{3..}) =  \frac{\sigma^2}{10}$. Therefore, an unbiased variance estimator is $\frac{MSE}{10}$. MSE can be computed as follows:
\[ MSE = \frac{SSE}{IJ(K-1)} = \frac{\sum_i \sum_k (y_{ik} - \bar{y}_{i \cdot})}{90}\]
From the ANOVA table in (b), $MSE = 0.173$. Based on this, we can use a T-test to test our hypothesis. The test statistic is: 
    \[T = \frac{\bar{y}_{5\cdot} - \bar{y}_{3\cdot}}{\sqrt{MSE/10}} = \frac{5.2925-4.9500}{\sqrt{0.173/10}}=2.6040\]
    Under $H_0$, the test statistic is distributed as $T(90)$. The corresponding p-value is 0.0103. We conclude that we have evidence that the two largest mean butterfat contents differ.}

   <!--  \textcolor{blue} {The breeds with the two largest mean butterfat contents are respectively Jersey and Guernsey. We formulate the following null hypothesis: \[H_0: \alpha_{\text{Jersey}} - \alpha_{\text{Guernsey}} = 0. \] -->
   <!-- To test the hypothesis, we can use a contrast which takes a value of 1 for Jersey, -1 for Guernsey and 0 for other breeds. Since there are 2 levels for the age variable, the difference at each level with the corresponding interaction must be equal under $H_0$. Therefore, we test for the difference of the differences in the mean butterfat content at the two levels of age. The output generated by R shows a p-value of 0.212. We conclude that we have evidence that the two largest mean butterfat contents are equal.} -->

\bigskip
  
2. (*FELM Exercise 2.2*) The dataset `wbca` in the library `faraway` comes from a study of breast cancer in Wisconsin (see `?wbca`). There are 681 cases of potentially cancerous tumours, of which 238 are actually malignant. Malignancy of a tumor is traditionally determined by an invasive surgical procedure.  The purpose of this study was to determine whether a new procedure called fine needle aspiration, which draws only a small sample of tissue, could be effective in determining tumor status.   
    (a) Fit a binomial regression with `Class` as the response and the other nine variables as predictors.  Report the residual deviance and associated degrees of freedom.  Can this information be used to determine if this model fits the data? Explain. \newline
    \textcolor{blue} {\textbf{Solution:} The residual deviance and associated df after fitting the model is shown in Table 4. The residual deviance is 89.46 with 671 df. Here, the residual deviance does not provide any useful information about the model fit. As we will show in Q3 - b) below, in the case of a binary response, the residual deviance is not a function of the response $y$ but the estimated regression parameters. }
```{r, echo=FALSE, results='asis'}
data(wbca)
model2 <- glm(Class~., data = wbca, family = binomial())
tb <- anova(model2)
tb <- round(tb, 4)
tb[is.na(tb)] <- ""
colnames(tb) <- c("Df (Deviance)", "Deviance", "Df (Residual Deviance)", "Residual Deviance")
kable(tb, 'latex', digits = 3, booktabs = T, align = "c")%>%
  kable_styling(latex_options = "HOLD_position")
```
    (b) Use the `step` function to choose final model. \newline
    \textcolor{blue} {\textbf{Solution:} The following table shows the summary of the final model selected from the backward stepwise regression:}
```{r, include=FALSE}
step_res <- step(model2)
```
```{r, echo=FALSE, results='asis'}
kable(summary(step_res)$coef, 'latex', digits = 3, booktabs = T, align = "c")%>%
  kable_styling(latex_options = "HOLD_position")
```
    (c) For a patient with the following measured values of the predictors:  `Adhes` = 1, `BNucl` = 1,  `Chrom` = 3,  `Epith` = 2,  `Mitos` = 1,  `NNucl` = 1,  `Thick` = 4,  `UShap` = 1,  `USize` = 1, predict the  `Class` of the tumor, using the full model and the final model. Give confidence intervals for your predictions, and explain how you obtained them. \newline
    \textcolor{blue} {\textbf{Solution:} On the logit scale, the prediction for the patient using the model selected from the backward stepwise regression is:
    \[ \hat{\beta}_{Adhes} \cdot 1 + \hat{\beta}_{BNucl} \cdot 1 + \hat{\beta}_{Chrom} \cdot 3 + \cdots + \hat{\beta}_{UShap} \cdot 1 = 4.834\]
    The variance of this prediction is:
    \begin{align*}
    \sigma^2_{\text{pred}} &= Var(\hat{\beta}_{Adhes} \cdot 1 + \hat{\beta}_{BNucl} \cdot 1 + \cdots + \hat{\beta}_{UShap} \cdot 1) \\
    &= Var(\hat{\beta}_{Adhes}) + Var(\hat{\beta}_{BNucl}) + \cdots + Var(\hat{\beta}_{USize}) + \\
    & 2 Cov(\hat{\beta}_{Adhes}, \hat{\beta}_{BNucl}) + 2 \times 3 Cov(\hat{\beta}_{Adhes}, \hat{\beta}_{Chrom}) \cdots + 2 \times 4 Cov(\hat{\beta}_{Thick}, \hat{\beta}_{UShap})
    \end{align*}
    An estimated variance-covariance matrix of $\hat{\boldsymbol{\beta}}$ can be obtained using an iterative method (See Q5 in Assignment 1). Based on the output from R, $\hat{\sigma^2}_{\text{pred}}= 0.3382$. Hence, a 95\% prediction intervals for this prediction on the logit scale is $4.834 \pm z_{0.975} \times \sqrt{0.3382}$ or $(3.695, 5.974)$. Since the logit function is strictly monotone, we apply its inverse (i.e. logistic function) to obtain the prediction interval on the original scale: $\left( \frac{1}{1+\exp(-3.695)}=0.9757, \frac{1}{1+\exp(-5.974)}=0.9975 \right )$. \newline
    Similarly, the same procedure can be repeated for the full model. The 95\% prediction interval is (0.9759, 0.9976).}
```{r, echo=FALSE, results='asis'}
newdata <- data.frame(Adhes=1, BNucl=1, Chrom=3, Epith=2, Mitos=1, NNucl=1, Thick=4,UShap=1, USize=1)
predict_full <- predict(model2, newdata = newdata, se.fit = T)
predict_step <- predict(step_res, newdata = newdata, se.fit = T)
predict_f_resp <- data.frame(estimate = exp(predict_full$fit)/(1+exp(predict_full$fit)), 
                            lower=exp(predict_full$fit-1.96*predict_full$se.fit)/(1+exp(predict_full$fit-1.96*predict_full$se.fit)),
                            upper=exp(predict_full$fit+1.96*predict_full$se.fit)/(1+exp(predict_full$fit+1.96*predict_full$se.fit)))
predict_s_resp <- data.frame(estimate = exp(predict_step$fit)/(1+exp(predict_step$fit)), 
                            lower=exp(predict_step$fit-1.96*predict_step$se.fit)/(1+exp(predict_step$fit-1.96*predict_step$se.fit)),
                            upper=exp(predict_step$fit+1.96*predict_step$se.fit)/(1+exp(predict_step$fit+1.96*predict_step$se.fit)))
```
    (d) Suppose that a cancer is classified as benign (`Class` = 1) if $p(\hat\beta) > 0.5$, and otherwise is classified as malignant.  Give a table showing the misclassification errors if this method is applied to the current data, using the final model from (b).  Give a similar table but using the cutoff $0.9$ instead. \newline
    \textcolor{blue} {\textbf{Solution}: The following table shows the misclassification error with a cut-off of 0.5 where the column represents the actual data and the row, the predicted outcome:}
```{r, echo=FALSE, results='asis'}
reference_0.5 <- as.factor(ifelse(wbca$Class > 0.5, "Benign", "Malignant"))
predicted <- predict(step_res, newdata = wbca, type = 'response')
predicted_0.5 <- as.factor(ifelse(predicted > 0.5, "Benign", "Malignant"))
cm0.5 <- confusionMatrix(predicted_0.5, reference_0.5)

reference_0.9 <- as.factor(ifelse(wbca$Class > 0.9, "Benign", "Malignant"))
predicted_0.9 <- as.factor(ifelse(predicted > 0.9, "Benign", "Malignant"))
cm0.9 <- confusionMatrix(predicted_0.9, reference_0.9)

cm0.5 <- cm0.5$table/681
cm0.9 <- cm0.9$table/681

kable(cm0.5, 'latex', digits = 6, booktabs = T, align = "c")%>%
  kable_styling(latex_options = "HOLD_position")
```
    \textcolor{blue} {The following table shows the misclassification error with a cut-off of 0.9:}
```{r, echo=FALSE, results='asis'}
kable(cm0.9, 'latex', digits = 6, booktabs = T, align = "c")%>%
  kable_styling(latex_options = "HOLD_position")
```
    (e)  The error estimates in (d) are overly optimistic, because they are computed on the same data that was used to fit the model.  A better strategy is to fit the model on training data and assess the errors on test data.  Choose a random 1/3 of the data to hold out as a test set, and find the best fitting model on the remainder.  Construct the tables of misclassification errors for the cutoffs 0.5 and 0.9, and compare to the results in (d).
 *Note: If you have issues with convergence on the randomly selected training set, you could try following Faraway's suggestion of choosing each 3rd case to construct the test set.*   \newline
```{r, include=FALSE}
idx <- seq(from=1, to=681, by=3)
train <- wbca[-idx,]
test <- wbca[idx,]

model3 <- glm(Class~., data = train, family = binomial())
final_train <- step(model3)

reference_0.5_test <- as.factor(ifelse(test$Class > 0.5, "Benign", "Malignant"))
predicted_test <- predict(model3, newdata = test, type = 'response')
predicted_0.5_test <- as.factor(ifelse(predicted_test > 0.5, "Benign", "Malignant"))
cm0.5_test <- confusionMatrix(predicted_0.5_test, reference_0.5_test)

reference_0.9_test <- as.factor(ifelse(test$Class > 0.9, "Benign", "Malignant"))
predicted_0.9_test <- as.factor(ifelse(predicted_test > 0.9, "Benign", "Malignant"))
cm0.9_test <- confusionMatrix(predicted_0.9_test, reference_0.9_test)

cm0.5_test <- cm0.5_test$table/227
cm0.9_test <- cm0.9_test$table/227
```
    \textcolor{blue} {\textbf{Solution}: With a random test dataset, the following table shows the misclassification error matrix with a cut-off of 0.5 where the column represents the actual data and the row, the predicted outcome:}
```{r, echo=FALSE, results='asis'}
kable(cm0.5_test, 'latex', digits = 6, booktabs = T, align = "c")%>%
  kable_styling(latex_options = "HOLD_position")
```
    \textcolor{blue} {The misclassification error matrix with a cut-off of 0.9 is shows as follows:}
```{r, echo=FALSE, results='asis'}
kable(cm0.9_test, 'latex', digits = 6, booktabs = T, align = "c")%>%
  kable_styling(latex_options = "HOLD_position")
```
    \textcolor{blue} {As expected, the error estimates, the off-digonal elements in the matrices above, increase when we use a random test dataset compared with the results in (d).}
    
\bigskip

3. *SM, Exercise 10.4.1* Suppose $y_1, \dots, y_n$ are independent Binomial$(n_i, p_i)$ random variables, with $p_i = \exp(x_i^{T}\beta)/\{1+\exp(x_i^{T}\beta)\}, i = 1, \dots, n$. 
    (a) Derive expressions for the log-likelihood function, the equations defining the maximum likelihood estimator of $\beta$, and the residual deviance. \newline
    \textcolor{blue} {\textbf{Solution}: The kernel of the likelihood function for our model is:
    \[ L(\boldsymbol{\beta}) \propto \prod_{i=1}^np_i^{y_i}(1-p_i)^{n_i-y_i} \]
    By taking the log, we obtain the log-likelihood:
    \[ l(\boldsymbol{\beta}) = \sum_{i=1}^ny_i\log p_i + (n_i - y_i)\log(1 - p_i) \]
    Taking the first derivative and setting it equal to 0, we obtain the score equations:
    \[ \begin{cases} 
    \frac{\partial l(\boldsymbol{\beta})}{\partial \beta_0} = \sum_{i=1}^n y_i(1-p_i) - (n_i - y_i) p_i = 0 \\
        \frac{\partial l(\boldsymbol{\beta})}{\partial \beta_j} = \sum_{i=1}^n x_{ij} y_i(1-p_i) - x_{ij} (n_i - y_i) p_i = 0 \quad (j=1, \cdots, q)
    \end{cases}\]
    It can be shown that the log-logistic function is concave. Since the log-likelihood is sum of log-logistic functions, it is itself concave. Hence, the score equations have a unique solution which is the MLE. \newline
    To obtain the residual deviance $D$, we derive the log-likelihood for the saturated model:
    \[ l_S(\boldsymbol{\beta}) = \sum_{i=1}^n y_i\log \left(\frac{y_i}{n_i}\right) + (n_i - y_i)\log \left(\frac{n_i - y_i}{n_i}\right) \]
    The log-likelihood for the fitted model is:
    \[ l_M(\boldsymbol{\beta}) = \sum_{i=1}^n y_i\log \hat{p}_i + (n_i - y_i)\log(1 - \hat{p}_i) \]
    Then, the residual deviance $D$ is:
    \begin{align*} D &= -2 (l_M(\boldsymbol{\beta}) - l_S(\boldsymbol{\beta})) \\
    &= -2 \sum_{i=1}^n y_i\log \hat{p}_i + (n_i - y_i)\log(1 - \hat{p}_i) - y_i\log \left(\frac{y_i}{n_i}\right) - (n_i - y_i)\log \left(\frac{n_i - y_i}{n_i}\right) \\
    &= 2  \sum_{i=1}^n y_i \log\left(\frac{y_i}{n_i \hat{p}_i}\right) + (n_i - y_i)\log\left(\frac{n_i - y_i}{n_i - n_i\hat{p}_i}\right)
    \end{align*}}
    (b) Use the approximation $\log (1+x) \simeq x - { x^2/2}$ to show that the residual deviance is approximately 
$$
X^2 = \sum_{i=1}^n \frac{(y_i - n_i\hat p_i)^2}{n_i\hat p_i(1-\hat p_i)}. $$  \newline
    \textcolor{blue} {\textbf{Solution}: From (a), the residual deviance $D$ is:
    \begin{align*} D &= 2  \sum_{i=1}^n y_i \log\left(\frac{y_i}{n_i \hat{p}_i}\right) + (n_i - y_i)\log\left(\frac{n_i - y_i}{n_i - n_i\hat{p}_i}\right) \\
    &= 2  \sum_{i=1}^n y_i \log\left(1 + \frac{y_i - n_i \hat{p}_i}{n_i \hat{p}_i} \right) + (n_i - y_i)\log\left(1 + \frac{n_i\hat{p}_i- y_i}{n_i - n_i\hat{p}_i}\right)
    \end{align*}
    Since $\hat{p}_i \rightarrow p_i$, $\frac{y_i - n_i \hat{p}_i}{n_i \hat{p}_i}$ and $\frac{n_i\hat{p}_i- y_i}{n_i - n_i\hat{p}_i}$ are small. Based on this, we can use the approximation $\log (1+x) \simeq x - { x^2/2}$ when $x \simeq 0$:
    \begin{align*} D &\simeq 2\sum_{i=1}^n \left[y_i\left(\frac{y_i - n_i\hat{p}_i}{n_i\hat{p}_i}\right) + (n_i - y_i)\left( \frac{n_i\hat{p}_i-y_i}{n_i-n_i\hat{p}_i}\right)\right] - \\
    & \sum_{i=1}^n\left[y_i\left(\frac{y_i - n_i\hat{p}_i}{n_i\hat{p}_i}\right)^2 +(n_i - y_i)\left( \frac{n_i\hat{p}_i-y_i}{n_i-n_i\hat{p}_i}\right)^2\right] \\
    &= 2\sum_{i=1}^n\frac{(y_i-n_i\hat{p}_i)^2}{n_i\hat{p}_i(1-\hat{p}_i)} - \\
    & \sum_{i=1}^n\left[(y_i - n_i\hat{p}_i +n_i\hat{p}_i) \left(\frac{y_i - n_i\hat{p}_i}{n_i\hat{p}_i}\right)^2 + (n_i - y_i + n_i\hat{p}_i - n_i\hat{p}_i)\left( \frac{n_i\hat{p}_i-y_i}{n_i-n_i\hat{p}_i}\right)^2\right]\\
    &= \sum_{i=1}^n\frac{(y_i-n_i\hat{p}_i)^2}{n_i\hat{p}_i(1-\hat{p}_i)} - \sum_{i=1}^n\left[\frac{(y_i - n_i\hat{p}_i)^3}{n_i^2\hat{p}_i^2} + \frac{(n_i\hat{p}_i - y_i)^3}{(n_i-n_i\hat{p}_i)^2}\right].
    \end{align*}
    Assuming the third order terms $\frac{(y_i - n_i\hat{p}_i)^3}{n_i^2\hat{p}_i^2} + \frac{(n_i\hat{p}_i - y_i)^3}{(n_i-n_i\hat{p}_i)^2}$ are small and $\simeq 0$, the residual deviance is approximately 
    \[ D \simeq X^2 = \sum_{i=1}^n\frac{(y_i-n_i\hat{p}_i)^2}{n_i\hat{p}_i(1-\hat{p}_i)} \]}
    (c) Show that if $n_i \equiv 1$, that the deviance is a function only of $\hat\beta$, and not otherwise a function of $y$. Show that if $n_i \equiv 1$ and $p_i \equiv p$, that $X^2 = n$. \textcolor{blue}{\it It is worth looking at Davison's version of this question in SM p.497.}  \newline
    \textcolor{blue} {\textbf{Solution}: With $n_i \equiv 1$, the score equations in (b) becomes the following:
    \[ \begin{cases} 
    \frac{\partial l(\boldsymbol{\beta})}{\partial \beta_0} = \sum_{i=1}^n y_i(1-p_i) - (1 - y_i) p_i = 0 \\
        \frac{\partial l(\boldsymbol{\beta})}{\partial \beta_j} = \sum_{i=1}^n x_{ij} y_i(1-p_i) - x_{ij} (1 - y_i) p_i = 0 \quad (j=1, \cdots, p)
    \end{cases}\]
    With some simplification:
    \[ \begin{cases} 
    \sum_{i=1}^n p_i = \sum_{i=1}^n y_i \\
    \sum_{i=1}^n x_{ij} p_i = \sum_{i=1}^n x_{ij} y_i \quad (j=1, \cdots, p)
    \end{cases}\]
    We will use these conditions for the proof. \newline
    The deviance $D$ is:
    \begin{align*}
    D &= -2\sum_{i=1}^n y_i\log \hat{p}_i + (1 - y_i)\log(1 - \hat{p}_i) \\
    &= -2\sum_{i=1}^n\log\left(\frac{\hat{p}_i}{1-\hat{p}_i}\right)y_i + \log(1 - \hat{p}_i)\\
    &= -2\sum_{i=1}^n \hat{\beta}_0 y_i + \hat{\beta}_1 x_{i1} y_i + \cdots + \hat{\beta}_p x_{ip} y_i  + \log(1 - \hat{p}_i)\\
    \end{align*}
    Using the conditions from the score equations, this becomes:
    \[ D = -2\sum_{i=1}^n \hat{\beta}_0 \hat{p}_i + \hat{\beta}_1 x_{i1} \hat{p}_i + \cdots + \hat{\beta}_p x_{ip} \hat{p}_i  + \log(1 - \hat{p}_i) \]
    This shows that the deviance is no longer a function of $y$ but $\hat{\boldsymbol{\beta}}$. \newline
    When $n_i\equiv 1$ and $p_i \equiv p$, the conditions from the score equations lead to $\hat{p} = \bar{y}$. Using this, the residual deviance is approximately:
    \begin{align*}
    X^2 &= \sum_{i=1}^n\frac{(y_i-n_i\hat{p}_i)^2}{n_i\hat{p}_i(1-\hat{p}_i)} \\ 
    &= \sum_{i=1}^n\frac{(y_i-\bar{y})^2}{\bar{y}(1-\bar{y})} \\
    &= \frac{\sum_{i=1}^n y_i^2 - n\bar{y}^2}{\bar{y}(1-\bar{y})}\\
    &= \frac{n\bar{y} - n\bar{y}^2}{\bar{y}(1-\bar{y})} \\
    & = n.
    \end{align*}}
    (d) Explain why the results in (c) mean that neither the residual deviance nor Pearson's $X^2$ can be used to assess goodness of fit with binary data.   \newline
    \textcolor{blue} {A goodness of fit requires a comparison of observed outcome and expected outcome. However, for a binary outcome, the observed data do not enter in the equation of residual deviance or Pearson's $\chi^2$ residuals. Therefore, both are not useful when the outcome is binary.}

\bigskip 

4. This question refers to the article "Clinical characteristics and predictors of outcomes of hospitalized patients with Coronavirus Disease 2019 in a multiethnic London National Health Service Trust: a retrospective cohort study" by Perez-Guzman et al. (2020) in *Clinical Infectious Diseases*. The article and supplementary material are posted on my web page for this class. 
    (a) What is a "retrospective cohort study", and why is this study described that way? \newline
    \textcolor{blue}{I was looking for an indication of "back in time". A cohort study records data on a group of individuals. Several definitions online and elsewhere refer to "treatment" being compared to "control", but that is not really the case here. The main point is that they assess vital status (on May 6), and related this to various covariates measured on admission.}
    (b) What is the "Elixhauser comorbidity score"? \newline
    \textcolor{blue}{This is an overall measure of health conditions; instead of trying to assess the effect of each one, an index measure like this can be used for data sets that aren't very large.}
    (c) What was the primary outcome of the study? \newline
    \textcolor{blue}{Vital status on May 6. (Some people interpreted "outcome" to mean the results of the study. But the authors do use the phrase "primary outcome" in their paper.)}
    (d) The methods section refers to "unadjusted" and "adjusted" logistic regression. What is the difference between these? \newline
    \textcolor{blue}{Unadjusted means a simple logistic regression with just one covariate. There are so many covariates, that the authors first did this to try to see which ones to keep in a multiple regression model. }
    (e) There are 9 models reported in Table S4 of the Supplementary Material. Which of these models ended up being summarized in  Table 2 of the main paper?  \newline
    \textcolor{blue}{I could not figure this out, so I gave everyone full marks for this. The descriptions are very confusing, and the numbers don't seem to match very well.}
    (f) The first 22 or 23 lines of Table 2 show one column of estimates called "unadjusted", and a second column "adjusted for age." As the authors note, "when accounting for age, only diabetes and chronic kidney disease remained statistically significant" (p.6). Why do you think this is?  \newline
    \textcolor{blue}{I think because most of the other measures are confounded with age. Probably diabetes and CKD are also confounded with age, but their effect on COVID survival is so large that it's important even after adjusting for age.}
    (g) On p.7 the authors state that they performed "unadjusted and adjusted logistic regression analyses to assess the odds of death of BAME groups compared to white patients". In your own words, what did these analyses show? \newline
    \textcolor{blue}{Once adjusted for age and co-morbidities, death rates are higher for Black patients. From another point of view the Black patients had similar mortality to whites, but they were younger on admission.}
    (h) In the discussion section, the authors say "our findings .. sugest that important biological differences may be potentially driving differences in COVID-19 hospitalization outcomes by ethnicity, which cannot solely be explained by socioeconomic differences". What is the basis for this claim?  \newline
    \textcolor{blue}{I'm not sure what the basis for this was. It might be their finding that Black patients did seem younger and with fewer co-morbidities, but had worse outcomes. In their introduction they mention that other studies that have attributed worse outcomes in BAME groups to socio-economic status did not have detailed information on their symptoms or underlying conditions when they arrived at the hospital. }
    (i) Summarize in your own words what you think are the two most important findings in this paper. \newline
    \textcolor{blue}{That Black, and to some extent Asian, patients have poorer outcomes after adjusting for co-morbidities and age seems very important. There is a huge effect of gender, which they don't discuss, perhaps because they thought it was well known. This paper had the most confusing tables ever, imho. } 

\bigskip

5. *Bonus Question* 
 The article by Roozenbeek et al. (2020) *RSOS* was discussed in class on October 15. The authors carried out several statistical analyses, as described in Section 2.3: Pearson's correlation coefficent; one-way analysis of variance; ordingary least squares regression;  two logistic regressions, and an OLS linear regression.
    (a) Choose two of these analyses and try to reproduce them using the authors' data.^[You are welcome to use their code or your own; I found it necessary to read the data using `read.csv`, for example, not `read.xls` as they recommend.] 
    (b) Try to reproduce one of the plots in the paper or the supplementary document.
    (c) In the supplementary information, Figure S1 shows residual plots from OLS regression. Why do you think there is an apparent linear boundary on the residuals vs fitted plot? 
    (d) In the paper (Footnote 10 on p.8), the authors refer to a "robust standard error regression". What method did they use for this? (It might be obvious from their code what `R` package they used, but I am looking for some detail on what statistical method they used.) 
    (e) The headline from  this paper stated "Poor numerical literacy linked to great susceptibility to Covid-19 fake news". Do you think this headline is supported by the analyses in the paper? Why or why not? 
    

<!-- ## Including Plots -->

<!-- You can also embed plots, for example: -->

<!-- ```{r pressure, echo=FALSE} -->
<!-- plot(pressure) -->
<!-- ``` -->

<!-- Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot. -->