---
title: "Homework 08 Solutions"
date: "11/21/2019"
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(knitr)
library(pander)
library(nlme)
library(ggplot2)
library(dplyr)
library(scales)
```

# 2 (a)
We first fit a model regression Sales on all of the covariates except Advert and Lag1Advert.
```{r 2a}
bookstore <- read.table("data/bookstore.txt", header=TRUE)
out1 <- lm(Sales ~ . -Advert - Lag1Advert, data=bookstore)
pander(summary(out1), caption = "Problem 2a: Summary of ordinary linear model.")
```

The diagnostic plots for our first model are shown in Figure \ref{f:2a_diagnostic}.  Our standard set of diagnostic plots, the first four plots, don't look too bad although there are some deviations from normality in the QQ plot and maybe an upward trend in the Scale-Location plot.  We see in the ACF plot that there is, statistically significant, auto-correlation in the lag1 term for the residuals.  We also plot the residuals against the Time variable, which could conceivably be used to address the auto correlation but we do not observe any trends in this plot.  We could estimate a model to account for the observed auto correlation but we will defer that until the next part since it is hypothesize that advertising in both the current and previous month may be important variables, and could conceivably account for the auto correlation.  We do note that the current model explains a large portion of the observed variation and is straight forward to interpret relying only on the month of the year and the Time (cumulative count of months).

```{r 2a_diagnostic, fig.cap="Problem 2a: Diagnostic plots for ordinary linear model.\\label{f:2a_diagnostic}", fig.height=4}
par(mfrow = c(2, 3),  mar = c(4, 5, 2, 1))
plot(out1)
acf(resid(out1), main = "Auto Correlation in Residuals")
plot(bookstore$Time, rstandard(out1), xlab = "Time", ylab = "Standardized Residuals")
abline(h =c(-2, 0, 2), lty = c(2,1,2))
```

# 2 (b)
We try adding the Advert and Lag1Advert 
```{r 2b}
out2 <- lm(Sales ~ ., data=bookstore)
pander(summary(out2), caption = "Problem 2b: Summary of ordinary linear model.")
```

We can see from the model outputs that neither Advert or Lag1Advert is statistically significant.  Conducting an ANOVA test (results below) we see that we fail to reject the null hypothesis that the model without either Advert or Lag1Advert performs as well as the larger model.  Examining the diagnostic plots for the larger model in Figure \ref{f:2b_diagnostic} we also the same issues as with the smaller model, namely: autocorrelation in the residuals and a slight trend in the scale-location plot.  Given the diagnostics we would prefer the model from part a) but might choose several methods for improving it.  We could try fitting a more general model to allow for the autocorrelation in the residuals.  Alternatively, we could remove or transform the Time variable yielding a smaller.  Since this variable is a just a linear trend it is possible that this could address the autocorrelation in the residuals observed in both models.

```{r 2b_diagnostic, fig.cap="Problem 2b: Diagnostic plots for ordinary linear model.\\label{f:2b_diagnostic}", fig.height=4}
par(mfrow = c(2, 3),  mar = c(4, 5, 2, 1))
plot(out2)
acf(resid(out2), main = "Auto Correlation in Residuals")
plot(bookstore$Time, rstandard(out2), xlab = "Time", ylab = "Standardized Residuals")
abline(h =c(-2, 0, 2), lty = c(2,1,2))
anova(out1, out2)
```

# 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("data/CarlsenQ.txt", header=TRUE)
m1 <- lm(Sales ~ . -Case - Time, data=carlsen)
pander(summary(m1), caption = "Problem 3a: Summary of ordinary linear model.")
```

```{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 = "")
```

# 3 (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)
```

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

# 3 (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 3c, 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 = "Problem 3d: 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 = "")
```

\newpage
# 4 (a)

We see in Figure \ref{f:4a_plot} that plotting the mean for each state is helpful in summarizing the overall income level within each stage but does not give an indication of the relationship between income and high school graduation rates.  Relative to the other plots Figure \ref{f:4a_plot} is easy to interpret but is probably not the best way to convey the information it is displaying.

```{r 4a_plot, fig.height=7.5, fig.cap="Problem 4a: HS graduation rate vs. county level per capita income with state level mean per capita income line.\\label{f:4a_plot}", fig.pos='b'}
dat <- read.table("data/cdi.dat")

mean.df <- group_by(dat, state) %>%
    summarise(per.cap.income = mean(per.cap.income))

ggplot(dat) +
    geom_point(aes(pct.hs.grad, per.cap.income)) +
    geom_hline(data = mean.df, aes(yintercept = per.cap.income), linetype = 2) +
    facet_wrap(~state, ncol = 5) +
    scale_y_continuous(label = comma) +
    xlab("Percentage High School Graduates") + ylab("Per Capitca Income ($)") +
    theme_minimal()
```

# 4 (b)

Figure \ref{f:4b_plot} allows us to examine if the relationship between high school graduation rate and income are similar to what is observed at a national level but it is hard to see what the state level trends are.  Relative to the others plots Figure \ref{f:4b_plot} is best suited the answering the question of ``Does state X look like nation as a whole?'' instead of conveying state level results or tends.

```{r 4b_plot, fig.height=7.5, fig.pos='h!', fig.cap="Problem 4b: HS graduation rate vs. county level per capita income with national level regression line.\\label{f:4b_plot}"}
fit <- lm(per.cap.income ~ pct.hs.grad, data = dat)
ggplot(dat, aes(pct.hs.grad, per.cap.income)) +
    geom_point() +
    geom_abline(intercept = fit[["coefficients"]][1],
                slope = fit[["coefficients"]][2],
                linetype = 2) +
    facet_wrap(~state, ncol = 5) +
    scale_y_continuous(label = comma) +
    xlab("Percentage High School Graduates") + ylab("Per Capitca Income ($)") +
    theme_minimal()
```

# 4 (c)

Figure \ref{f:4c_plot} allows us to examine if the relationship between high school graduation rate and income are similar to what is observed at a national level while adjusting for mean income in the state but makes it harder to see how the overall state income level compares to the nation as a whole.  Relative to the others plots Figure \ref{f:4b_plot} is best suited the answering the question of ``Does the relationship between high gradation and income in state X look like nation as a whole?'' but does not convey much about the mean income level as done in Figure \ref{f:4a_plot} and Figure \ref{f:4b_plot}.

```{r 4c_plot, fig.height=7.5, fig.cap="Problem 4c: HS graduation rate vs. county level per capita income with state level intercept and national slope line.\\label{f:4c_plot}"}
fit2 <- lm(per.cap.income ~ -1 + state + pct.hs.grad, data = dat)
slopes.df <- data.frame(state = substr(names(fit2[["coefficients"]][1:48]), 6, 7),
                       state.intercept = fit2[["coefficients"]][1:48],
                       slope = fit2[["coefficients"]][49])

ggplot(dat, aes(pct.hs.grad, per.cap.income)) +
    geom_point() +
    geom_abline(data = slopes.df, aes(intercept = state.intercept, slope = slope), linetype = 2) +
    facet_wrap(~state, ncol = 5) +
    scale_y_continuous(label = comma) +
    xlab("Percentage High School Graduates") + ylab("Per Capitca Income ($)") +
    theme_minimal()
```

# 4 (d)

We see in Figure \ref{f:4d_plot} that the state level regression lines match the observed data, as closely as possible, but we are unable to fit them for all states since some states contain only a single county.  In contrast to the other plots there is much greater variance in the fitted regression lines perhaps better describing each state but Figure \ref{f:4d_plot} fails to convey any country level trends.

```{r 4d_plot, fig.height=7.5, fig.cap="Problem 4d: HS graduation rate vs. county level per capita income with state level regression line.\\label{f:4d_plot}", warning=FALSE}
state_lm <- function(df){
    mod <- lm(per.cap.income ~ pct.hs.grad, data = df)
    return(data.frame(intercept = mod[["coefficients"]][1], slope = mod[["coefficients"]][2]))
}
lm.df <- group_by(dat, state) %>% do(., state_lm(.))

ggplot(dat, aes(pct.hs.grad, per.cap.income)) +
    geom_point() +
    geom_abline(data = lm.df, aes(intercept = intercept, slope = slope), linetype = 2) +
    facet_wrap(~state, ncol = 5) +
    scale_y_continuous(label = comma) +
    xlab("Percentage High School Graduates") + ylab("Per Capitca Income ($)") +
    theme_minimal()

```