\documentclass{article}

\frenchspacing

\usepackage{fullpage,txfonts,graphicx,subfig}

\long\def\comment#1{}

\usepackage{Sweave}
\begin{document}
\input{hw03-solutions-concordance}

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

\begin{enumerate} % \itemsep=0pt\parsep=0pt\parskip=0pt


\item{} {\em Please do Sheather, Ch 3, p. 105, \#3, {\bf Part A}.  Remember that the data is in the ``0-textbook'' folder in the files area on Canvas for this class.}\label{3a}

\begin{enumerate}

\item {\em  Develop a simple linear regression model based on least squares that predicts
advertising revenue per page from circulation (i.e., feel free to transform
either the predictor or the response variable or both variables). Ensure that
you provide justification for your choice of model.}


First, a quick look at the data... \\
{\sf (You do not have to produce any EDA for your answert to part (a)).}

\begin{Schunk}
\begin{Sinput}
> magdata <- read.csv("AdRevenue.csv",header=T)
> str(magdata,width=72,strict.width = "cut")
\end{Sinput}
\begin{Soutput}
'data.frame':	70 obs. of  4 variables:
 $ Magazine                  : chr  "People" "Better Homes and Garden"..
 $ PARENT.COMPANY..SUBSIDIARY: chr  "Time Warner, (Time Inc.)" "Mered"..
 $ AdRevenue                 : num  233 397 286 877 304 ...
 $ Circulation               : num  3.75 7.64 4.07 32.7 3.21 ...
\end{Soutput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> par(mfrow=c(1,3))
> hist(magdata$AdRevenue,main="")
> hist(magdata$Circulation,main="")
> plot(AdRevenue ~ Circulation, data=magdata)
\end{Sinput}
\end{Schunk}

The plot is shown in Figure~\ref{f:initEDA}. We can see from the figure that there is severe right skewing in both the Circulation and AdRevenue variables.  

% \setkeys{Gin}{width=0.6\textwidth}
\begin{figure}[hb]
\begin{center}
\includegraphics[height=2in,width=4.5in]{hw03-solutions-fig-1.pdf}
\end{center}
\caption{Initial EDA.}\label{f:initEDA}
\end{figure}
% \setkeys{Gin}{width=0.8\textwidth}

Next, we fit a linear regression model to the untransformed variables: \\
{\sf (This is a good baseline model to fit when you are considering transformations, but you do not need to include it in your answer to part (a) for this assignment.)}

\clearpage

\begin{Schunk}
\begin{Sinput}
> lm.1 <- lm(AdRevenue ~ Circulation, data=magdata)
> summary(lm.1)
\end{Sinput}
\begin{Soutput}
Call:
lm(formula = AdRevenue ~ Circulation, data = magdata)

Residuals:
     Min       1Q   Median       3Q      Max 
-147.694  -22.939   -7.845   13.810  131.130 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  99.8095     5.8547   17.05   <2e-16 ***
Circulation  22.8534     0.9518   24.01   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 42.22 on 68 degrees of freedom
Multiple R-squared:  0.8945,	Adjusted R-squared:  0.8929 
F-statistic: 576.5 on 1 and 68 DF,  p-value: < 2.2e-16
\end{Soutput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> par(mfrow=c(2,2))
> plot(lm.1)
\end{Sinput}
\end{Schunk}
The diagnostic plots are shown in Figure~\ref{f:untrans}. We can see that the residuals are also right-skewed and have non-constant variance.

\setkeys{Gin}{width=0.6\textwidth}
\begin{figure}[hb]
\begin{center}
\includegraphics{hw03-solutions-fig-2}
\end{center}
\caption{Diagnostics for regression on untransformed variables.}\label{f:untrans}
\end{figure}
\setkeys{Gin}{width=0.8\textwidth}

We will try two transformations:\\
{\sf (You should have something like these two analyses in your answer to part (a).)}
\begin{itemize}
\item Just a log transformation on both variables (since logarithms reduce right-skew, and if it works it gives us an interpretable model);
\item The power transformations suggested by Box-Cox.
\end{itemize}

Here's the ``log everything model'':
\begin{Schunk}
\begin{Sinput}
> lm.2 <- lm(log(AdRevenue) ~ log(Circulation),data=magdata)
> summary(lm.2)
\end{Sinput}
\begin{Soutput}
Call:
lm(formula = log(AdRevenue) ~ log(Circulation), data = magdata)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.47022 -0.11142 -0.00532  0.10835  0.42705 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       4.67473    0.02525  185.16   <2e-16 ***
log(Circulation)  0.52876    0.02356   22.44   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1768 on 68 degrees of freedom
Multiple R-squared:  0.881,	Adjusted R-squared:  0.8793 
F-statistic: 503.6 on 1 and 68 DF,  p-value: < 2.2e-16
\end{Soutput}
\begin{Sinput}
> par(mfrow=c(2,2))
> plot(lm.2)
\end{Sinput}
\end{Schunk}
The diagnostic plots are shown in Figure \ref{sf:log-diag}, page \pageref{sf:log-diag}.  \\

For the Box-Cox transformations, I suggest
\begin{itemize}
\item First, find the best (rounded) Box-Cox power for $x$.
\item Then, for the model $y \sim $(transformed $x$), find the best (rounded) Box-Cox power for $y$.  
\end{itemize}
That way, you are using the information you have learned about $x$ to produce the best possible transformation of $y$.

First, the suggested transformation for $x$ = Circulation:
\begin{Schunk}
\begin{Sinput}
> library(car)
> with(magdata,boxCox(Circulation~1))
> with(magdata,powerTransform(Circulation~1)$roundlam)
\end{Sinput}
\begin{Soutput}
  Y1 
-0.5 
\end{Soutput}
\end{Schunk}
The profile likelihood is shown in Figure \ref{sf:bc-x}, page \pageref{sf:bc-x}.

Next, the Box-Cox transformation for $y$ = AdRevenue, using the transformed $x$, $1/\sqrt{\mbox{Circulation}}$:
\begin{Schunk}
\begin{Sinput}
> lm.3 <- lm(AdRevenue ~ I(Circulation^(-0.5)),data=magdata)
> with(magdata,boxCox(lm.3))
> with(magdata,powerTransform(lm.3)$roundlam)
\end{Sinput}
\begin{Soutput}
Y1 
-1 
\end{Soutput}
\end{Schunk}
The profile likelihood is shown in Figure \ref{sf:bc-y}.

\begin{figure}
\begin{center}
\subfloat[boxCox profile likelihood for $x\sim1$.  The ``rounded'' value of $\lambda$ is $-0.5$, so the transformation is $1/\sqrt{x}$.]{
\includegraphics[width=0.4\textwidth]{hw03-solutions-boxCox-L.pdf} % &
\label{sf:bc-x}
}
\hspace{4em}
\subfloat[boxCox profile likelihood for $y\sim 1/\sqrt{x}$. The ``rounded'' value of $\lambda$ is $-1$, so the transformation is $1/y$.]{\includegraphics[width=0.4\textwidth]{hw03-solutions-boxCox-R.pdf} 
\label{sf:bc-y}}
\end{center}
\caption{Selecting the boxCox power for $x\sim1$ and for $y\sim 1/\sqrt{x}$; $x$ = Circulation, $y$ = AdRevenue.}\label{f:bc-profile}
\end{figure}

So, the final model suggested by Box-Cox is $1/\mbox{(AdRevenue)} \sim 1/\sqrt{\mbox{(Circulation)}}$: 
\begin{Schunk}
\begin{Sinput}
> magdata$AdRevInv <- 1/magdata$AdRevenue
> magdata$InvSqrtCirc <- 1/sqrt(magdata$Circulation)
> lm.4 <- lm(AdRevInv ~ InvSqrtCirc,data=magdata)
> # the following caused an error in R, which is why I defined the variables above...
> # lm.4 <- lm(I(AdRevenue^(-1)) ~ I(Circulation^(-0.5)),data=magdata)
> summary(lm.4)
\end{Sinput}
\begin{Soutput}
Call:
lm(formula = AdRevInv ~ InvSqrtCirc, data = magdata)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0028448 -0.0008745 -0.0000689  0.0006133  0.0040733 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.0001662  0.0004000   0.416    0.679    
InvSqrtCirc 0.0091424  0.0004571  20.000   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.001223 on 68 degrees of freedom
Multiple R-squared:  0.8547,	Adjusted R-squared:  0.8526 
F-statistic:   400 on 1 and 68 DF,  p-value: < 2.2e-16
\end{Soutput}
\begin{Sinput}
> par(mfrow=c(2,2))
> plot(lm.4)
\end{Sinput}
\end{Schunk}
The diagnostic plots are shown in Figure \ref{sf:bc-diag}.

\begin{figure}
\begin{center}
\subfloat[Regression diagnostics for $\log y \sim \log x$.]{
\includegraphics[width=0.4\textwidth]{hw03-solutions-diag-L.pdf} 
\label{sf:log-diag}
}
\hspace{4em}
\subfloat[Regression diagnostics for $1/y \sim 1/\sqrt{x}$.]{
\includegraphics[width=0.4\textwidth]{hw03-solutions-diag-R.pdf}
\label{sf:bc-diag}
}
\end{center}
\caption{Comparing regression diagnostics for the model $\log y \sim \log x$, vs. the ``best'' boxCox model $1/y \sim 1/\sqrt{x}$; $x$ = Circulation, $y$ = AdRevenue.}\label{f:diag-compare}
\end{figure}

{\sf (Your choice of model, and justification, should be similar to the following.)}\\
Both models have high $R^2$'s: for lm.2 ($log y \sim log x$), $R^2 = 0.88$ and for lm.4  ($1/y \sim 1/\sqrt{x}$), $R^2 = 0.85$\footnote{Note that even though the original model lm.1 ($y\sim x$) had an even higher $R^2=0.89$, we do not seriously consider it since the regression diagnostic plots (Figure \ref{f:untrans}) are so bad!}.  Comparing coefficient estimates, 
\begin{Schunk}
\begin{Sinput}
> summary(lm.2)$coef
\end{Sinput}
\begin{Soutput}
                 Estimate Std. Error   t value     Pr(>|t|)
(Intercept)      4.674734 0.02524738 185.15717 1.168946e-93
log(Circulation) 0.528758 0.02356174  22.44138 3.754210e-33
\end{Soutput}
\begin{Sinput}
> summary(lm.4)$coef
\end{Sinput}
\begin{Soutput}
                Estimate   Std. Error    t value     Pr(>|t|)
(Intercept) 0.0001662359 0.0004000311  0.4155575 6.790423e-01
InvSqrtCirc 0.0091424191 0.0004571138 20.0003113 3.417303e-30
\end{Soutput}
\end{Schunk}
we see that the slope estimates for both models are highly significantly different from zero, so the regression output doesn't help us distinguish between the models very well.  The regression diagnostics for both models are also very similar (Figure \ref{f:diag-compare}): in both models the residuals are much more nearly normal, they have nearly constant variance, and there are almost no highly influential data points. \\  

Since both models fit similarly well, I prefer to use the model that is easier to talk about: $\log y \sim \log x$.  Referring to the coefficient tables above, we see that for each 1\% change in circulation, we can expect about a 0.53\% change in ad revenue. \\

{\sf (Here is where you identify which magazines are associated with unusual observations.)}\\
Referring back to Figure \ref{sf:log-diag}, we see that the most extreme residual outliers are observations 5, 60, and 64:

\clearpage

\begin{Schunk}
\begin{Sinput}
> with(magdata,
+      data.frame(Magazine,Circulation,AdRevenue,
+                 StdRes=rstandard(lm.2),leverage=hatvalues(lm.2))[c(5,60,64),])
\end{Sinput}
\begin{Soutput}
             Magazine Circulation AdRevenue    StdRes   leverage
5  Sports Illustrated       3.205   304.185  2.440209 0.02022629
60         Prevention       3.347   127.315 -2.668848 0.02115025
64      Cooking Light       1.717    89.153 -2.678798 0.01432293
\end{Soutput}
\begin{Sinput}
> p <- 1 # number of predictors: x only
> c(leverage.cutoff = 2*(p+1)/dim(magdata)[1])
\end{Sinput}
\begin{Soutput}
leverage.cutoff 
     0.05714286 
\end{Soutput}
\end{Schunk}
We see that {\em Sports Illustrated} ad revenue overperforms relative to its circulation, while both {\em Cooking Light} and {\em Prevention} (a health magazine) underperform.  None of these has very high leverage, however (using the rule of thumb that leverage above $4/n$ is ``high'').

\item {\em Find a 95\% prediction interval for the advertising revenue per page for
magazines with the following circulations:
\begin{enumerate}
\item 0.5 million
\item 20 million
\end{enumerate}
}

{\sf (Your answer should pretty much go like this.)}\\
This is slightly tricky, because we have to account for the transformation of $y$ in model lm.2 in producing the appropriate prediction interval:

\begin{enumerate}
\item 0.5 million
\begin{Schunk}
\begin{Sinput}
> print(pred.i <- predict(lm.2,newdata=data.frame(AdRevenue=0,Circulation=0.5),
+                         interval="prediction"))
\end{Sinput}
\begin{Soutput}
       fit      lwr    upr
1 4.308227 3.947855 4.6686
\end{Soutput}
\begin{Sinput}
> print(interval.i <- exp(pred.i[c(2,3)]))  # Have to un-do the logarithm
\end{Sinput}
\begin{Soutput}
[1]  51.82406 106.54846
\end{Soutput}
\end{Schunk}
So, we see that a magazine with circulation of half a million could expect ad revenue between \$51,820 and \$106,550.  

\item 20 million
\begin{Schunk}
\begin{Sinput}
> print(pred.ii <- predict(lm.2,newdata=data.frame(AdRevenue=0,Circulation=20),
+                         interval="prediction"))
\end{Sinput}
\begin{Soutput}
       fit      lwr      upr
1 6.258752 5.885815 6.631689
\end{Soutput}
\begin{Sinput}
> print(interval.ii <- exp(pred.ii[c(2,3)]))  # Have to un-do the logarithm
\end{Sinput}
\begin{Soutput}
[1] 359.8958 758.7626
\end{Soutput}
\end{Schunk}


So, a magazine with circulation of 20 million could expect ad revenue between \$359,900 and \$758,760.  
\end{enumerate}


\item {\em Describe any weaknesses in your model.}

{\sf (Your answer should pretty much go like this, though you might also reproduce the summary()'s here.)}\\
We already discussed the summary in part (a): it shows a high $R^2$ of 0.88, and a highly significant slope estimate (same with the $F$ statistic for overall fit).  Referring again to Figure \ref{sf:log-diag}:
\begin{itemize}
\item The residuals vs fitted plot doesn't really show any problems.  Except for a couple of outliers, that show up more clearly in the QQ plot, there really isn't any vertical pattern or curve that the residuals follow.
\item The QQ plot shows that the residuals are following the normal distribution fairly well.  There is still a bit of right skewing in the residuals, and two low outliers {\em Prevention} and {\em Cooking Light}.
\item The scale-location plot doesn't really show any serious problems, though there may be a bit of increasing variance as the predicted ad revenue increases.
\item The residuals vs leverage plot shows a couple of high-leverage points, but no really high Cook's Distance values, so these points are not very influential on the fit of the model lm.2; some more detail on the high Cook's Distance points is given in the R output below:\\
{\sf (You don't have to do the following for your answer, but it is interesting.)}
\begin{Schunk}
\begin{Sinput}
> res.lev <- data.frame(magdata$Magazine,StdRes=rstandard(lm.2),
+                       leverage=hatvalues(lm.2),Cooks.Dist=cooks.distance(lm.2))
> tail(res.lev[order(cooks.distance(lm.2)),])
\end{Sinput}
\begin{Soutput}
     magdata.Magazine     StdRes   leverage Cooks.Dist
64      Cooking Light -2.6787982 0.01432293 0.05213714
5  Sports Illustrated  2.4402089 0.02022629 0.06146309
20    Reader's Digest -1.3103123 0.06716783 0.06181267
49  AARP The Magazine -0.9279932 0.13138896 0.06513180
60         Prevention -2.6688479 0.02115025 0.07695149
4          Parade (1)  1.5938160 0.16374980 0.24870866
\end{Soutput}
\begin{Sinput}
> p <- 1 # number of predictors: x only
> c(leverage.cutoff = 2*(p+1)/dim(magdata)[1])
\end{Sinput}
\begin{Soutput}
leverage.cutoff 
     0.05714286 
\end{Soutput}
\end{Schunk}
\end{itemize}
Overall, the log-log model seems to fit the data well.

\end{enumerate}

\item{} {\em Please do Sheather, Ch 3, p. 105, \#3, {\bf Part B}.}\label{3b}

\begin{enumerate}

\item {\em Develop a polynomial regression model based on least squares that directly
predicts the effect on advertising revenue per page of an increase in circula-
tion of 1 million people (i.e., do not transform either the predictor nor the
response variable). Ensure that you provide detailed justification for your
choice of model. [Hint: Consider polynomial models of order up to 3.]}

Sheather suggests we try polynomials up to order 3. We will go even higher, to order 5, just to see what happens.  To save some space, I will just print out tables of estimated coefficients and standard errors, and $R^2$'s:\\
{\sf (You might arrive at the cubic model differently.)}
\begin{Schunk}
\begin{Sinput}
> lm.5 <- lm(AdRevenue ~ Circulation + I(Circulation^2) + I(Circulation^3) + 
+              I(Circulation^4) + I(Circulation^5), data=magdata)
> ## Note: you could get the same model with
> ## lm.5 <- lm(AdRevenue ~ poly(Circulation, degree=5, raw=T),data=magdata)
> summary(lm.5)$r.squared
\end{Sinput}
\begin{Soutput}
[1] 0.9366882
\end{Soutput}
\begin{Sinput}
> summary(lm.5)$coefficients
\end{Sinput}
\begin{Soutput}
                      Estimate   Std. Error    t value    Pr(>|t|)
(Intercept)      57.0382867767 1.737595e+01  3.2825999 0.001668725
Circulation      47.8288287740 2.213601e+01  2.1606799 0.034469648
I(Circulation^2)  1.3621060276 7.962225e+00  0.1710710 0.864707640
I(Circulation^3) -0.6557047334 9.924238e-01 -0.6607104 0.511169347
I(Circulation^4)  0.0370875190 4.489949e-02  0.8260120 0.411865983
I(Circulation^5) -0.0005798354 6.533235e-04 -0.8875166 0.378124078
\end{Soutput}
\end{Schunk}
If we remove some high powers of $x$, we may find a change in the significance of lower powers (this can happen because powers of $x$ can be {\em collinear}). So we will also try the model with powers up to order 4 and order 3 only. 
\begin{Schunk}
\begin{Sinput}
> lm.6 <- lm(AdRevenue ~ Circulation + I(Circulation^2) + I(Circulation^3) + 
+              I(Circulation^4), data=magdata)
> summary(lm.6)$r.squared
\end{Sinput}
\begin{Soutput}
[1] 0.935909
\end{Soutput}
\begin{Sinput}
> summary(lm.6)$coefficients
\end{Sinput}
\begin{Soutput}
                     Estimate   Std. Error   t value     Pr(>|t|)
(Intercept)      45.686665431 11.742208057  3.890807 2.377459e-04
Circulation      65.380572047  9.928749452  6.584976 9.325058e-09
I(Circulation^2) -5.499586772  1.900414956 -2.893887 5.173919e-03
I(Circulation^3)  0.220172602  0.104527351  2.106364 3.903839e-02
I(Circulation^4) -0.002733043  0.001694505 -1.612886 1.116144e-01
\end{Soutput}
\begin{Sinput}
> lm.7 <- lm(AdRevenue ~ Circulation + I(Circulation^2) + I(Circulation^3),
+            data=magdata)
> summary(lm.7)$r.squared
\end{Sinput}
\begin{Soutput}
[1] 0.933344
\end{Soutput}
\begin{Sinput}
> summary(lm.7)$coefficients
\end{Sinput}
\begin{Soutput}
                    Estimate  Std. Error   t value     Pr(>|t|)
(Intercept)      59.17036829 8.345045881  7.090478 1.118099e-09
Circulation      51.23581639 4.711234296 10.875243 2.334496e-16
I(Circulation^2) -2.50537894 0.411411261 -6.089719 6.476556e-08
I(Circulation^3)  0.05222479 0.009229702  5.658339 3.574381e-07
\end{Soutput}
\end{Schunk}
All three models have $R^2\approx 0.93$, but only model lm.7, with powers just up to order 3, has all of its $\hat\beta$'s significiantly different from zero.  Since this seems like the best model so far, we consider diagnostic plots for it:\\
{\sf (Always good to look at diagnostic plots!)}
\begin{Schunk}
\begin{Sinput}
> par(mfrow=c(2,2))
> plot(lm.7)
\end{Sinput}
\end{Schunk}
The plots appear in Figure \ref{f:poly3diag}.  They suggest that this model isn't doing as well\footnote{We could improve things a bit by fitting $\log y \sim x + x^2 + x^3$ rather than $y \sim x + x^2 + x^2$ but it is still not better than the log-log model we chose above (try it!).  Box-Cox on $y$ also does not help, and anyway these ideas go beyond what Sheather is asking for.} as the log-log model above!

\begin{figure}
\begin{center}
\includegraphics[width=0.5\textwidth]{hw03-solutions-polydiag}
\end{center}
\caption{Diagnostics for the 3$^{rd}$ order polynomial model.}\label{f:poly3diag}
\end{figure}




\item {\em Find a 95\% prediction interval for the advertising page cost for magazines
with the following circulations:
\begin{enumerate}
\item 0.5 million
\item 20 million
\end{enumerate}
}

{\sf (Your answer should pretty much go like this.)}\\
This goes just like the predictions we did earlier, except that now since there is no transformation on $y$, we can use the prediction intervals directly.

\clearpage

\begin{Schunk}
\begin{Sinput}
> print(interval.i <- predict(lm.7,newdata=data.frame(AdRevenue=0,Circulation=0.5),
+                         interval="prediction")[c(2,3)])
\end{Sinput}
\begin{Soutput}
[1]  14.92314 153.41378
\end{Soutput}
\begin{Sinput}
> print(interval.ii <- predict(lm.7,newdata=data.frame(AdRevenue=0,Circulation=20),
+                         interval="prediction")[c(2,3)])
\end{Sinput}
\begin{Soutput}
[1] 418.1790 580.8878
\end{Soutput}
\end{Schunk}

So from this model 
\begin{enumerate}
\item A magazine with a circulation of half a million could anticipate ad revenue between \$14,920 and \$153,410.
\item A magazine with a circulation of 20 million could anticipate ad revenue between \$418,180 and \$580,890.
\end{enumerate}

\item {\em Describe any weaknesses in your model.}

{\sf (Again, you might copy the model summary here, but otherwise your answer should be pretty similar to the below.)}\\
The summary of model lm.7 above shows an $R^2= 0.93$, and all significant predictors, which is great.  Referring to Figure \ref{f:poly3diag}, we see, however, that 
\begin{itemize}
\item The residuals vs fitted plot shows a lot of right-skew in the fitted values, pretty good symmetry of the residuals around zero, but some pretty large outliers.
\item The QQ plot shows a that both the right and left tails or the residual distribution are longer than the normal distribution's tails, with the right tail even longer than the left, and some very large outliers.
\item The scale-location plot doesn't show much stability around 1, but part of the problem may be that the fitted values are so skewed-right that there just isn't much data on the right side of the plot from which to make the (red) loess line. 
\item The residuals vs. leverage plot shows some very large outliers with large leverage values; these also have large Cook's distances, which suggests that they are influential on the fit of lm.7.  A few details on large outliers/leverage points are given below.
\begin{Schunk}
\begin{Sinput}
> res.lev <- data.frame(magdata$Magazine,StdRes=rstandard(lm.7),
+                       leverage=hatvalues(lm.7),Cooks.Dist=cooks.distance(lm.7))
> tail(res.lev[order(cooks.distance(lm.7)),],n=7)
\end{Sinput}
\begin{Soutput}
           magdata.Magazine    StdRes   leverage  Cooks.Dist
5        Sports Illustrated  3.125952 0.03102333  0.07821313
46         American Profile -1.509367 0.14002065  0.09273286
2  Better Homes and Gardens  2.179137 0.13037457  0.17797941
20          Reader's Digest -2.702626 0.17242129  0.38044641
49        AARP The Magazine -2.765478 0.47272305  1.71414898
8               USA Weekend  3.705447 0.47304388  3.08140011
4                Parade (1) -1.478471 0.99099778 60.15734624
\end{Soutput}
\begin{Sinput}
> p <-3 # Number of predictors: x, x^2, X^3
> c(leverage.cutoff = 2*(p+1)/dim(magdata)[1])
\end{Sinput}
\begin{Soutput}
leverage.cutoff 
      0.1142857 
\end{Soutput}
\end{Schunk}
Several magazines have outlying residuals or high leverage; three magazines have high Cook's distances, with higher influence on the fit of lm.7.  The last of these, {\em Parade}, has leverage near one ($h_{ii} = 0.99$!), and a whopping 60 for Cook's Distance!
\end{itemize}
Overall, this is not a great model. Although the regression output is good, the diagnostic plots reveal many problems with the fit.

\end{enumerate}

\clearpage

\item{} {\em Please do Sheather, Ch 3, p. 105, \#3, {\bf Part C}.}\label{3c}

\begin{enumerate}

\item {\em Compare the model in Part A with that in Part B. Decide which provides a
better model. Give reasons to justify your choice.}

{\sf (Your answer should go pretty much like this.}\\
All of the models have high $R^2$'s and significant predictors (except for some predictors in the higher order polynomial models) .  The differences really come in the diagnostic plots.  Referring to Figures \ref{f:diag-compare} and \ref{f:poly3diag}, as well as our summaries for Parts A(c) and B(c) above, we see that 
\begin{itemize}
\item The residual diagnostic plots for the polynomial regression model (Figure \ref{f:poly3diag} show severe model deficiencies, whereas the plots for either the log-log or Box-Cox models Figure \ref{f:diag-compare} show closer agreement with the regression assumptions.  So we should pick one of the models from Part A.
\item The two models from Part A perform about equally well, but the log-log model is easier to explain and talk about than the Box-Cox model.
\end{itemize}
For these reasons, I prefer the log-log model (lm.2) to all the others tried.

\item {\em Compare the prediction intervals in Part A with those in Part B. In each case,
decide which interval you would recommend. Give reasons to justify each
choice.}

{\sf (Your answer should be pretty much like this, although if you don't go into an explanation for the different interval widths, that's OK.)}\\
Here is a table comparing the prediction intervals:

\begin{center}
\let\mc\multicolumn
\begin{tabular}{rllll}
       & \mc{2}{c}{Circulation = 0.5 million} & \mc{2}{c}{Circulation = 20 million} \\
Model  & Low Endpoint & High Endpoint & Low Endpoint & High Endpoint\\
\hline
lm.2 (log-log)    & 51.82 & 106.55 & 359.9 & 758.76 \\
lm.7 (polynomial) & 14.92 & 153.41 & 418.18 & 580.89 
\end{tabular}
\end{center}
For the lower circulation, the log-log interval is narrower.  For the higher circulation, the polynomial interval is narrower.  We can say a little more: Both models give wider intervals for higher circulations; this is just because higher circulations are farther from the average circulation, and SE$_{pred}$($\hat y$) = $S\sqrt{1 + (x - \overline x)^2/SXX}$ increases with this distance.  The width of the log-log intervals grows faster because we also had to exponentiate the endpoints to get from log(\$\$) intervals back to just \$\$ intervals.  Calculations below:
\begin{Schunk}
\begin{Sinput}
> ## polynomial interval widths:
> data.frame(circ0.5mil=c(width = 153.41 - 14.92), circ20mil = c(width = 580.89 - 418.18))
\end{Sinput}
\begin{Soutput}
      circ0.5mil circ20mil
width     138.49    162.71
\end{Soutput}
\begin{Sinput}
> ## log-log interval widths:
> data.frame(circ0.5mil = c(width = 106.55 - 51.82), circ20mil = c(width = 758.76 - 359.9))
\end{Sinput}
\begin{Soutput}
      circ0.5mil circ20mil
width      54.73    398.86
\end{Soutput}
\end{Schunk}
I prefer the log-log intervals (even though the second one is quite wide) because they are based on a better-fitting model.

\end{enumerate}

\item{}\label{idmrad} {\em Write a brief IDMRAD paper based on your answers to problems \ref{3a}--\ref{3c}.\label{report}  Remember to label the {\bf Introduction}, {\bf Data}, {\bf Methods}, {\bf Results} and {\bf Discussion} and {\bf Technical Appendix} sections.  }

This report appears at the end of these solutions!

\clearpage

\comment{

Here are some additional instructions:
\begin{description}
\item[Title] Choose a brief title that gives the reader some idea of what the data are and what the result is.  
\item[Author] Your name, your institution (Department of Statistics and Data Science, Carnegie Mellon University), and your email address.
\item[Abstract] {\em We will skip the abstract for this exercise}.
\item[Introduction] This can largely parrot information in the problem statement before {\bf Part A} on p. 105 of Sheather.  The main research problem is to develop a good regresson model to predict gross advertising revenue per advertising page in 2006
(in thousands of dollars) from circulation (in millions).
\item[Data] First, say where the data came from (again, you can parrot part of the problem statement before {\bf Part A}). Then, provide (a) a table with the variable names and their definitions (including units of measurement), (b) a set of summary statistics for the quantitative variables (e.g. mean, SD, etc.), and (c) any exploratory plots that seem useful (histogram, boxplots,  scatter plot, etc.), but don't be repetitive (for example, for some distributions, histograms and boxplots tell the same story, so you wouldn't want to include both).  {\em Include at least one sentence about each of these three elements in the Data section}.
\item[Methods\footnotemark]\footnotetext{Note that I am citing parts of the Appendix for any reader who wants to check the technical details.} The methods section will read something like this (not much more than a paragraph or two, since this is a simple analysis): ``To build a good regression model for predicting gross advertising revenue per advertising page from circulation, we considered two approaches.  First we considered simple linear regression, with various transformations for the circulation variable (see Appendix, Part A).  The methods we used to find suitable transformations were {\em blah-blah-blah}.  The final transformation was adjusted to improve interpretation and clarity of communication.  Then we considered polynomial regression models (see Appendix, Part B), i.e., regressing gross advertising revenue on circulation, circulation$^2$, and so forth, up to the power {\em blah-blah}.   We selected the best model from both the transformed simple linear regression and the polynomial regression models by considering casewise residual diagnostic plots, {\em and blah-blah-blah} (see Appendix, Part C).'' 
\item[Results\footnotemark]\footnotetext{Again, it would be appropriate to cite parts of the Appendix for any reader who wishes to examine the technical details.}  This will be three short sections, perhaps just three paragraphs.  First, remind the reader what methods you used in Part A, give the best model from Part A, and provide a table of coefficient estimates and standard errors.  Second, same thing for Part B.  Finally, summarize your findings for Part C.  Normally I would {\em not} include regression diagnostic plots without a very good reason, but in this problem if there are a couple of plots that really explain why you chose the model you did, you can include it/them here.
\item[Discussion] Re-state the research question, and state the main result.  Interpret the result: what does it say about the relationship between gross advertising revenue per advertising page, and circulation?  Is there anything else the reader should know, or keep in mind when considering or using this result?  Finally, discuss any weakneses of your analyses, and suggest possible future work (could be work that fixes the weakness, or work that extends these results to a larger data set or a different setting, etc.)
\item[References] {\em We will skip the references for this exercise.}
\item[Technical Appendix] For the technical appendix, include the questions and the answers to problems \ref{3a}, \ref{3b} and \ref{3c} above.
\end{description} 

}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\item{} {\em [Based on Gelman \& Hill. Ch 3, \#1, p. 49] The file {\tt pyth.dat}, in the same folder as this hw,  contains outcome {\tt y} and inputs {\tt x1}, {\tt x2} for 40 data points, with a further 20 points with the inputs but no observed outcome (for this problem we will ignore these last 20 points). Save the file to your working directory and read it into R using the {\tt read.table()} function.}

\begin{enumerate}

\item\label{gh-a} {\em Fit the two models
\begin{eqnarray*}  
	{\bf M1:~} {\tt y} & = & \beta_0 + \beta_1 {\tt x1} + \varepsilon \\
	{\bf M2:~} {\tt y} & = & \beta_0 + \beta_1 {\tt x2} + \varepsilon 
\end{eqnarray*}
Which model provides a better fit for {\tt y}?  Why?}
\begin{Schunk}
\begin{Sinput}
> gh.data <- read.table("pyth.dat",header=T)
> gh.data <- gh.data[apply(gh.data,1,function(x) {!any(is.na(x))}),]
> str(gh.data)
\end{Sinput}
\begin{Soutput}
'data.frame':	40 obs. of  3 variables:
 $ y : num  15.68 6.18 18.1 9.07 17.97 ...
 $ x1: num  6.87 4.4 0.43 2.73 3.25 5.3 7.08 9.73 4.51 6.4 ...
 $ x2: num  14.09 4.35 18.09 8.65 17.68 ...
\end{Soutput}
\begin{Sinput}
> M1 <- lm(y ~ x1, data=gh.data)
> M2 <- lm(y ~ x2, data=gh.data)
> summary(M1)
\end{Sinput}
\begin{Soutput}
Call:
lm(formula = y ~ x1, data = gh.data)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.7409 -4.5056  0.7114  4.3739  7.7547 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  10.0633     1.5526   6.481 1.25e-07 ***
x1            0.6559     0.2499   2.625   0.0124 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.921 on 38 degrees of freedom
Multiple R-squared:  0.1535,	Adjusted R-squared:  0.1312 
F-statistic:  6.89 on 1 and 38 DF,  p-value: 0.01242
\end{Soutput}
\begin{Sinput}
> summary(M2)
\end{Sinput}
\begin{Soutput}
Call:
lm(formula = y ~ x2, data = gh.data)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.1751 -1.2352 -0.1867  1.0899  5.3755 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.78532    0.66037   5.732 1.33e-06 ***
x2           0.83223    0.05017  16.589  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.863 on 38 degrees of freedom
Multiple R-squared:  0.8787,	Adjusted R-squared:  0.8755 
F-statistic: 275.2 on 1 and 38 DF,  p-value: < 2.2e-16
\end{Soutput}
\begin{Sinput}
> ##
> par(mfrow=c(2,2))
> plot(M1)
> plot(M2)
\end{Sinput}
\end{Schunk}
The plots appear in Figure \ref{f:gh-a}.  The $R^2$ is much lower for M1 (0.1535) than for M2 (0.8787).  Neither set of residual diagnostic plots looks great: the residual vs fitted and scale-location plots somewhat favor M1, and the QQ plots somewhat favor M2. The Cook's distances are a bit better for M1 also.   For $R^2$ and normality of residuals, I prefer M2.

\begin{figure}
\begin{center}
\subfloat[Diagnostic plots for M1.]{
\includegraphics[width=0.4\textwidth]{hw03-solutions-M1-diag}
}
\hspace{\fill}
\subfloat[Diagnostic plots for M2.]{
\includegraphics[width=0.4\textwidth]{hw03-solutions-M2-diag}
}
\end{center}
\caption{Diagnostic plots for problem \ref{gh-a}.}\label{f:gh-a}
\end{figure}

\item\label{gh-b} {\em Construct new variables {\tt y2 = y$^2$}, {\tt x12 = x1$^2$}, and {\tt x22 = x2$^2$} and fit the models
\begin{eqnarray*}  
	{\bf M3:~} {\tt y2} & = & \beta_0 + \beta_1 {\tt x12} + \varepsilon \\
	{\bf M4:~} {\tt y2} & = & \beta_0 + \beta_1 {\tt x22} + \varepsilon 
\end{eqnarray*}
Compare the fits of these two models to the models in part (a).  Which fits best?  Why?}

\begin{Schunk}
\begin{Sinput}
> attach(gh.data)
> y2 <- y^2
> x12 <- x1^2
> x22 <- x2^2
> detach()
> M3 <- lm(y2 ~ x12,data=gh.data)
> M4 <- lm(y2 ~ x22,data=gh.data)
> summary(M3)
\end{Sinput}
\begin{Soutput}
Call:
lm(formula = y2 ~ x12, data = gh.data)

Residuals:
     Min       1Q   Median       3Q      Max 
-189.324 -125.674    4.988  131.052  214.089 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 161.7831    31.7873   5.090    1e-05 ***
x12           1.2971     0.6242   2.078   0.0445 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 131.1 on 38 degrees of freedom
Multiple R-squared:  0.102,	Adjusted R-squared:  0.07841 
F-statistic: 4.318 on 1 and 38 DF,  p-value: 0.04452
\end{Soutput}
\begin{Sinput}
> summary(M4)
\end{Sinput}
\begin{Soutput}
Call:
lm(formula = y2 ~ x22, data = gh.data)

Residuals:
    Min      1Q  Median      3Q     Max 
-41.280 -31.224  -7.463  25.422  59.571 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  35.1583     9.0306   3.893 0.000387 ***
x22           1.0198     0.0419  24.338  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 33.96 on 38 degrees of freedom
Multiple R-squared:  0.9397,	Adjusted R-squared:  0.9381 
F-statistic: 592.3 on 1 and 38 DF,  p-value: < 2.2e-16
\end{Soutput}
\begin{Sinput}
> par(mfrow=c(2,2))
> plot(M3)
> plot(M4)
\end{Sinput}
\end{Schunk}
The plots are in Figure \ref{f:gh-b}.  Model M4 has the highest $R^2$ (0.9397), and has residuals vs fitted and scale-location plots that are at least as good as any of the others; on the other hand, we seem to be losing normality of the residuals.  Nevertheless I prefer M4 so far.

\begin{figure}
\begin{center}
\subfloat[Diagnostic plots for M3.]{
\includegraphics[width=0.4\textwidth]{hw03-solutions-M3-diag}
}
\hspace{\fill}
\subfloat[Diagnostic plots for M4.]{
\includegraphics[width=0.4\textwidth]{hw03-solutions-M4-diag}
}
\end{center}
\caption{Diagnostic plots for problem \ref{gh-b}.}\label{f:gh-b}
\end{figure}



\item\label{gh-c} {\em To fit the model 
\[
	{\tt y2} ~ = ~ \beta_0 + \beta_1 {\tt x1} + \beta_2 {\tt x2} + \varepsilon ~,
\]
we just expand the R modeling language a little bit: \verb|y ~ x1 + x2|.  Fit both of the models
\begin{eqnarray*}
	{\bf M5:~} {\tt ~y} & = & \beta_0 + \beta_1 {\tt x1} + \beta_2 {\tt x2} + \varepsilon \\
	{\bf M6:~} {\tt y2} & = & \beta_0 + \beta_1 {\tt x12}  + \beta_2 {\tt x22} + \varepsilon 
\end{eqnarray*}
Compare these to the earlier models.  Which fits best?  Why?}

\begin{Schunk}
\begin{Sinput}
> M5 <- lm(y ~ x1 + x2, data=gh.data)
> M6 <- lm(y2 ~ x12 + x22, data=gh.data)
> summary(M5)
\end{Sinput}
\begin{Soutput}
Call:
lm(formula = y ~ x1 + x2, data = gh.data)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.9585 -0.5865 -0.3356  0.3973  2.8548 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.31513    0.38769   3.392  0.00166 ** 
x1           0.51481    0.04590  11.216 1.84e-13 ***
x2           0.80692    0.02434  33.148  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9 on 37 degrees of freedom
Multiple R-squared:  0.9724,	Adjusted R-squared:  0.9709 
F-statistic: 652.4 on 2 and 37 DF,  p-value: < 2.2e-16
\end{Soutput}
\begin{Sinput}
> summary(M6)
\end{Sinput}
\begin{Soutput}
Call:
lm(formula = y2 ~ x12 + x22, data = gh.data)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.26020 -0.05391 -0.00396  0.06367  0.35990 

Coefficients:
             Estimate Std. Error  t value Pr(>|t|)    
(Intercept) 0.0026691  0.0422669    0.063     0.95    
x12         0.9999672  0.0006419 1557.713   <2e-16 ***
x22         0.9998685  0.0001663 6011.909   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1344 on 37 degrees of freedom
Multiple R-squared:      1,	Adjusted R-squared:      1 
F-statistic: 2.013e+07 on 2 and 37 DF,  p-value: < 2.2e-16
\end{Soutput}
\begin{Sinput}
> par(mfrow=c(2,2))
> plot(M5)
> plot(M6)
\end{Sinput}
\end{Schunk}
The plots appear in Figure \ref{f:gh-c}.  \\

Putting both x1 and x2 in the model for y really improved the model: M5 has an $R^2$ of 0.9724, and both predictors are significant (have coefficient estimates significantly different from zero).  However the residual diagnostic plots don't look great; in particular it seems like the residuals have a lot of right skew, and leverage seems tio increase with the size of the standardized resduals. \\

M6 is really winning, though: $R^2=1$, and the QQ plot shows good agreement between the residuals and the normal distribution.  There seems to be some evidence for non-constant variance though: the residuals vs fitted plot fans out as the fitted valuess increase, and the scale-location plot tells a similar story.  On the other hand, only one data point seems to have a concerning Cook's distance.\\

Based on all of this I like M6 best.  Looking at the estimated coefficients for M6, I notice something interesting: $\hat\beta_0$ is indistinguishable from 0, and both $\hat\beta_1$ and $\hat\beta_2$ equal 1, to at least two decimal places (even if we compute the 95\% CI's!).

\begin{figure}
\begin{center}
\subfloat[Diagnostic plots for M5.]{
\includegraphics[width=0.4\textwidth]{hw03-solutions-M5-diag}
}
\hspace{\fill}
\subfloat[Diagnostic plots for M6.]{
\includegraphics[width=0.4\textwidth]{hw03-solutions-M6-diag}
}
\end{center}
\caption{Diagnostic plots for problem \ref{gh-c}.}\label{f:gh-c}
\end{figure}


\item\label{gh-d} {\em Can you find a simple, recognizable function {\tt x3 = (something involving both x1 and x2)}, so that 
\[
		{\bf M7:~} {\tt y} ~ = ~ \beta_0 + \beta_1 {\tt x3} + \varepsilon
\]
provides a fit comparable to the best fitting models above?  What is going on?}

In problem \ref{gh-c} we saw that the model M6 was very nearly
\[
    y^2 = x1^2 + x2^2 + \varepsilon
\]
If we ignore $\varepsilon$, take square roots, and put $\varepsilon$ and some ``unknown'' regression coefficients back in, we get a model like
\[
    y = \beta_0 + \beta_1\sqrt{x1^2 + x2^2} + \varepsilon
\]
i.e., $y$ is the distance to the origin from some points $(x1,x2)$ in Cartesian space.

Let's try fitting this model:
\begin{Schunk}
\begin{Sinput}
> M7 <- lm(y ~ I(sqrt(x1^2 + x2^2)),data=gh.data)
> summary(M7)
\end{Sinput}
\begin{Soutput}
Call:
lm(formula = y ~ I(sqrt(x1^2 + x2^2)), data = gh.data)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0083283 -0.0027000 -0.0007907  0.0031643  0.0089809 

Coefficients:
                      Estimate Std. Error  t value Pr(>|t|)    
(Intercept)          0.0018422  0.0019159    0.962    0.342    
I(sqrt(x1^2 + x2^2)) 0.9998313  0.0001316 7596.431   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.00434 on 38 degrees of freedom
Multiple R-squared:      1,	Adjusted R-squared:      1 
F-statistic: 5.771e+07 on 1 and 38 DF,  p-value: < 2.2e-16
\end{Soutput}
\begin{Sinput}
> par(mfrow=c(1,4))
> plot(M7)
\end{Sinput}
\end{Schunk}
\begin{figure}
\begin{center}
\def\bj{1.2in}
\begin{tabular}{cccc}
\includegraphics[height=\bj,width=\bj]{hw03-solutions-M7-diag1}
&
\includegraphics[height=\bj,width=\bj]{hw03-solutions-M7-diag2}
&
\includegraphics[height=\bj,width=\bj]{hw03-solutions-M7-diag3}
&
\includegraphics[height=\bj,width=\bj]{hw03-solutions-M7-diag5}
\end{tabular}
\end{center}
\caption{Diagnostic plots for model M7 (problem \ref{gh-d}).}\label{f:gh-d}
\end{figure}

The residual diagnostic plots are in Figure \ref{f:gh-d}.  This seems to confirm our suspicions!  $R^2=1$, the estimated regression coefficients are essentially $\hat\beta_0 = 0$ and $\hat\beta_1=1$, and the diagnostic plots looks great: 
\begin{itemize}
\item The residual vs fitted plot shows little vertical structure.
\item The QQ plot shows good adherence to normality.
\item The scale-location plot is consistent with constant-variance residuals.
\item None of the data points has Cook's distance above 0.5.
\end{itemize}


\end{enumerate}


\end{enumerate}


\clearpage

\setcounter{figure}{0} 
\setcounter{page}{1}

\begin{center}
{\Large\bf
Understanding the Relationship Between Circulation Size and Ad Revenue for a Selection of U.S. Consumer Magazines
}\\[2em]

Brian Junker, Department of Statistics and Data Science\\[2em]
brian@stat.cmu.edu
\end{center}

\section{Introduction}\label{s:intro}

The price of advertising (and hence revenue from advertising) is different from
one consumer magazine to another. Publishers of consumer magazines argue that magazines that reach more readers create more value for the advertiser. Thus, circulation is an important factor that affects revenue from advertising. In this report, we investigate the relationship between circulation and gross advertising revenue.  

In particular we will 
\begin{itemize}
\item Develop regression models to predict gross advertising revenue per advertising page in 2006 (in thousands of dollars) from circulation (in millions); and 
\item Illustrate the effect of circulation on ad revenue with two prediction intervals.
\end{itemize}

\section{Data}\label{s:data}

The data are for the top 70 US magazines ranked in terms of total gross advertising revenue in 2006.  The data were obtained from http://adage.com and are given in the file AdRevenue.csv which is available on the book web site. 

The variables in the data set are shown in Table \ref{t:variables}.

\begin{table}[h]
\begin{center}
\begin{tabular}{rl}
Variable & Definition \& Comments\\
\hline
Magazine & The name of each magazine for which data was collected\\
PARENT.COMPANY & The parent company or subsidiary which publishes this magazine\\
AdRevenue & The magazine's revenue per advertising page in 2006 (in thousands of dollars)\\
Circulation & The number of subscribers (in millions) to this magazine
\end{tabular}
\end{center}
\caption{Variable Definitions for the {\tt AdRevenue.csv} data set.}\label{t:variables}
\end{table}
Summary statistics for the two quanitative variables are given in Table \ref{t:sum-stat}.  Further EDA in Appendix \ref{a:eda} (page \pageref{a:eda} below) shows that both variables are substantially skewed right, but that an increasing relationship between the variables is plausible.
\begin{table}[h]
\begin{center}
\begin{Schunk}
\begin{Soutput}
              Min.   1st Qu.   Median       Mean   3rd Qu.    Max.
Circulation  0.331   0.99225   1.6755   3.118471   2.74325  32.700
AdRevenue   61.101 104.85050 133.7940 171.077200 179.39750 876.907
\end{Soutput}
\end{Schunk}
\end{center}
\caption{Summary Statistics for AdRevenue and Circulation.}\label{t:sum-stat}
\end{table}

\section{Methods}\label{s:methods}

To develop a regression model to predict AdRevenue from Circulation, we considered regression models using logarithmic and Box-Cox power transforms of the variables AdRevenue and Circulation (Appendix \ref{a:log-log}) as well as regression of AdRevenue on polynomial functions of Circulation (Appendix \ref{a:poly}), up to order 5.  We chose our final model based on a summary of each regression analysis and an examiniation of residual diagnostic plots.

With our final model, we calculated AdRevenue intervals in which we would expect to find 95\% of companies with circulations of 0.5 million subscribers and 20 million subscribers, respectively, accommodating transformation of variables, if any.

\section{Results}\label{s:results}

We considered regressions using the original variables AdRevenue and Circulation (details in Appendix \ref{a:log-log}, p. \pageref{aa:original}), logarithmic and power transformations (Appendix \ref{a:log-log}, pp. \pageref{aa:log-log}ff.) and polynomial regression using polynomials in Circulation of orders 3, 4 and 5 (Appendix \ref{a:poly}).  All approaches produced models with high $R^2$ values and highly significant predictors of AdRevenue. 

\subsection*{Logarithmic and Power Transformations}

Among models with logarithmic and power transformations, the models with the best residual diagnostic plots (Appendix \ref{a:log-log}, pages \pageref{aa:log-log-diag} and \pageref{aa:box-cox-diag}) were the following two models, shown with estimated regression coefficients:
\begin{equation}
      \log({\rm AdRevenue}) = 4.67 + 0.53\cdot\log({\rm Circulation}) + \varepsilon
\label{e:log-log}
\end{equation}
and
\begin{equation}
      1/{\rm (AdRevenue)} =  0.0002 + 0.0091\cdot1/\sqrt{\rm (Circulation)} + \varepsilon
\label{e:box-cox}
\end{equation}
Models (\ref{e:log-log}) and (\ref{e:box-cox})
had similar $R^2$ values (0.881 and 0.8547, respectively) and similarly good residual diagnostic plots. Because the log-log model is more easily interpreted in terms of percent-change, we prefer model (\ref{e:log-log}).  Table \ref{t:coef-table} gives the full table of estimated coefficients and standard errors for model (\ref{e:log-log}).

\begin{table}[hb]
\begin{center}
\begin{Schunk}
\begin{Soutput}
                 Estimate Std. Error   t value     Pr(>|t|)
(Intercept)      4.674734 0.02524738 185.15717 1.168946e-93
log(Circulation) 0.528758 0.02356174  22.44138 3.754210e-33
\end{Soutput}
\end{Schunk}
\end{center}
\caption{Estimated coefficients and standard errors for model (\ref{e:log-log}).}\label{t:coef-table}
\end{table}

\subsection*{Polynomial Regression Models}

We fitted polynomial regression models of order 3, 4 and 5, and found that the model of order 3, shown here with estimated regression coefficients, 
\begin{equation}
  {\rm (AdRevenue)} = 59.17 + 51.24\cdot{\rm (Circulation)} - 
  2.51\cdot{\rm (Circulation)}^2 + 0.05\cdot{\rm (Circulation)}^3 + \varepsilon
  \label{e:poly}
\end{equation}
was the most successful---all the predictors were significant, and $R^2\approx 0.93$.  Table \ref{t:poly-coef} gives the estimated coefficients and standard errors for this model.

\begin{table}
\begin{center}
\begin{Schunk}
\begin{Soutput}
                    Estimate  Std. Error   t value     Pr(>|t|)
(Intercept)      59.17036829 8.345045881  7.090478 1.118099e-09
Circulation      51.23581639 4.711234296 10.875243 2.334496e-16
I(Circulation^2) -2.50537894 0.411411261 -6.089719 6.476556e-08
I(Circulation^3)  0.05222479 0.009229702  5.658339 3.574381e-07
\end{Soutput}
\end{Schunk}
\end{center}
\caption{Estimated coefficients and standard errors for model (\ref{e:poly}).}\label{t:poly-coef}
\end{table}

Nevertheless, the residual diagnostic plots for this model (Appendix \ref{a:poly}, p. \pageref{aa:poly-diag}) did not look as good as the diagnostics for the log-log model (\ref{e:log-log}).


\subsection*{Final Model}

All three models (\ref{e:log-log}), (\ref{e:box-cox}) and (\ref{e:poly}) have high $R^2$'s and highly significant predictors.  Based on residual diagnostic plots (pages \pageref{aa:log-log-diag}, \pageref{aa:box-cox-diag} and \pageref{aa:poly-diag} in the Appendix), we can eliminate model (\ref{e:poly}), which does not follow the assumptions of the linear model as well as the other two.  Models (\ref{e:log-log}) and (\ref{e:box-cox}) have very similar residual diagnostics, so we are free to choose based on interpretability.   

Since the model with logarithms has a simpler interpretation (a 1\% change in Circulation is associated with an expected change of 0.53\% in AdRevenue), we chose model (\ref{e:log-log}) as our final model.  Figure \ref{f:regression-line} shows the fitted regression line under model (\ref{e:log-log}), laid over the raw data.  

\begin{figure}[h]
\begin{center}
\subfloat[Original scale.]{
\includegraphics[width=0.4\textwidth]{hw03-solutions-regline}
}
\hspace{\fill}
\subfloat[Log-log scale.]{
\includegraphics[width=0.4\textwidth]{hw03-solutions-logline}
}
\end{center}
\caption{The fitted model (\ref{e:log-log}) overlaid on the raw data.  Points in red correspond to magazines in Table \ref{t:outliers}.}\label{f:regression-line}
\end{figure}

The three magazines that followed our final model (\ref{e:log-log}) least well are listed in Table \ref{t:outliers}; nevertheless they did not influence the fit of the model very much.  We can see from the table that {\em Sports Illustrated} outperforms its predicted ad revenue, and both {\em Cooking Light} and {\em Prevention} underperform their predicted ad revenues.

\begin{table}
\begin{center}
\begin{Schunk}
\begin{Soutput}
             Magazine Circulation Predicted.AdRevenue Actual.AdRevenue
5  Sports Illustrated       3.205              198.46          304.185
60         Prevention       3.347              203.06          127.315
64      Cooking Light       1.717              142.68           89.153
\end{Soutput}
\end{Schunk}
\end{center}
\caption{Magazines with unusually high or low ad revenues (in thousands of dollars), given their circulation sizes (in millions of subscriptions), relative to their predicted ad revenue under model (\ref{e:log-log}).  These magazines are marked as red points in Figure \ref{f:regression-line}.}\label{t:outliers}
\end{table}

\subsection*{Prediction Intervals}

Using our final model model (\ref{e:log-log}), we predict 95\% of consumer magazines could expect ad revenues in the following intervals, based on their circulations:
\begin{itemize}
\item For a magazine with a circulation of 0.5 million subscriptions, the predicted interval for ad revenue is \$51,820 to \$106,550.
\item For a magazine with a circulation of 20 million subscriptions, the predicted interval for ad revenue is \$359,900 to \$758,760.
\end{itemize}

\section{Discussion}\label{s:discussion}

Among models we considered (power transformations in simple linear regression, as well as polynomial regression), we found that the model that fits the relationship between Circulation and AdRevenue best is a log-log model, shown here with estimated regression coefficients:
$$
      \log({\rm AdRevenue}) = 4.67 + 0.53\cdot\log({\rm Circulation}) + \varepsilon
      \eqno{(\ref{e:log-log})}
$$
The variable $\log({\rm Circulation})$ is a highly significant predictor of $\log({\rm AdRevenue})$; the variation in predicted \\ $\log({\rm AdRevenue})$ accounts for $R^2\cdot 100\% = 88.1\%$ of the variation in raw $\log({\rm AdRevenue})$.  The model can be interpreed as saying that we expect a 0.53\% change in Ad Revenue for every 1\% change in Circulation.  The relationship in model (\ref{e:log-log}) is illustrated in Figure \ref{f:regression-line} above.  

We also calculated intervals predicting a range of Ad Revenues for magazines with circulations of 0.5 million and 20 million subscriptions. As expected, the larger the circulation, the wider the range of possible ad revenues.

These calculations are helpful in determining whether particular magazines are over- or under-performing what we would expect, and we illustrated this with three magazines with the most unusual ad revenues for their circulation sizes, according to model \ref{e:log-log}; see Table \ref{t:outliers} above.

A key limitation of this work is that the data is quite old, from 2006.  This limits the generalizability of the results to the present time: the publishing industry has continued to undergo enormous upheavals due to competition from ``free'' content available on the internet; although we might still expect a log-log relationship to hold up with more current data, we would expect the estimated regression coeffcients (at least) to change.  

It also might be useful to have more than 70 magazines, especially to assess whether the relationship holds up for lower-circulation or lower-revenue magazines, and whether the relationship changes from one magazine genre (e.g. sports magazines) to another (e.g. health magazines).

\clearpage

\appendix

\renewcommand\thefigure{\thesection.\arabic{figure}} 

\section{Initial Look at the Data}\label{a:eda}

We begin by reading in the data and taking a quick look at it:
\begin{Schunk}
\begin{Sinput}
> magdata <- read.csv("AdRevenue.csv",header=T)
> str(magdata,width=72,strict.width = "cut")
\end{Sinput}
\begin{Soutput}
'data.frame':	70 obs. of  4 variables:
 $ Magazine                  : chr  "People" "Better Homes and Garden"..
 $ PARENT.COMPANY..SUBSIDIARY: chr  "Time Warner, (Time Inc.)" "Mered"..
 $ AdRevenue                 : num  233 397 286 877 304 ...
 $ Circulation               : num  3.75 7.64 4.07 32.7 3.21 ...
\end{Soutput}
\begin{Sinput}
> par(mfrow=c(1,3))
> hist(magdata$AdRevenue,main="")
> hist(magdata$Circulation,main="")
> plot(AdRevenue ~ Circulation, data=magdata)
\end{Sinput}
\end{Schunk}
\includegraphics[height=2in,width=4.5in]{hw03-solutions-fig-1.pdf}

\begin{Schunk}
\begin{Sinput}
> rbind(
+ Circulation = summary(magdata$Circulation),
+ AdRevenue = summary(magdata$AdRevenue)
+ )
\end{Sinput}
\begin{Soutput}
              Min.   1st Qu.   Median       Mean   3rd Qu.    Max.
Circulation  0.331   0.99225   1.6755   3.118471   2.74325  32.700
AdRevenue   61.101 104.85050 133.7940 171.077200 179.39750 876.907
\end{Soutput}
\end{Schunk}
We see from the plots that both AdRevenue and Circulation are highly skewed-right.  However, there does seem to be a linear relationship between these two variables.

\section{Simple Regression, Transformed Variables}\label{a:log-log}
\setcounter{figure}{0}

We tried 
\begin{itemize}
\item Simple regression on the original variables: AdRevenue $\sim$ Circulation.
\item Simple regression on the logs of the variables: log(AdRevenue) $\sim$ log(Circulation).
\item Simple regression with Box-Cox power transformations of the variables; this model turned out to be 1/AdRevenue $\sim$ 1/$\sqrt{\mbox{Circulation}}$.
\end{itemize}
The best model turned out to be the log-log model.  Some details of our analyses follow:

\subsection*{Original variables}\label{aa:original}

The regression output and residual diagnostic plots for the model AdRevenue $\sim$ Circulation are as follows:

\begin{Schunk}
\begin{Sinput}
> summary(lm.1 <- lm(AdRevenue ~ Circulation,data=magdata))
\end{Sinput}
\begin{Soutput}
Call:
lm(formula = AdRevenue ~ Circulation, data = magdata)

Residuals:
     Min       1Q   Median       3Q      Max 
-147.694  -22.939   -7.845   13.810  131.130 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  99.8095     5.8547   17.05   <2e-16 ***
Circulation  22.8534     0.9518   24.01   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 42.22 on 68 degrees of freedom
Multiple R-squared:  0.8945,	Adjusted R-squared:  0.8929 
F-statistic: 576.5 on 1 and 68 DF,  p-value: < 2.2e-16
\end{Soutput}
\begin{Sinput}
> par(mfrow=c(2,2))
> plot(lm.1)
\end{Sinput}
\end{Schunk}
\begin{center}
\includegraphics[width=0.5\textwidth]{hw03-solutions-fig-2}
\end{center}

Although $R^2=0.8945$ and Circulation is a highly statistically significant predictor, the residual diagnostic plots show skewing in the residuals, to go along with the skewing in AdRevenue and Circulation that we saw in the exploratory plots in Section \ref{a:eda}.

\subsection*{Log-transformed variables}\label{aa:log-log}

The regression output and residual diagnostic plots for the model log(AdRevenue) $\sim$ log(Circulation) are as follows:

\begin{Schunk}
\begin{Sinput}
> summary(lm.2 <- lm(log(AdRevenue) ~ log(Circulation),data=magdata))
\end{Sinput}
\begin{Soutput}
Call:
lm(formula = log(AdRevenue) ~ log(Circulation), data = magdata)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.47022 -0.11142 -0.00532  0.10835  0.42705 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       4.67473    0.02525  185.16   <2e-16 ***
log(Circulation)  0.52876    0.02356   22.44   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1768 on 68 degrees of freedom
Multiple R-squared:  0.881,	Adjusted R-squared:  0.8793 
F-statistic: 503.6 on 1 and 68 DF,  p-value: < 2.2e-16
\end{Soutput}
\begin{Sinput}
> par(mfrow=c(2,2))
> plot(lm.2)
\end{Sinput}
\end{Schunk}
\begin{center}
\includegraphics[width=0.5\textwidth]{hw03-solutions-diag-L}\label{aa:log-log-diag}
\end{center}
For this model, $R^2=0.881$ is still high, log(Circulation) is still a strong predictor, and the residual diagnostic plots look much better: residuals show no severe vertical patterns, they follow the normal distribution except for a small number of outliers, the location-scale plot shows at most mild violations of non-constant variance, and no data points with high Cook's distances.\\

We can look at the data points with the highest (but still not concerning) Cook's distances to see what ``extreme'' data looks like for this model; points with high leverage tend to have small residuals, and vice-versa:
\begin{Schunk}
\begin{Sinput}
> res.lev <- data.frame(Magazine=magdata$Magazine,StdRes=rstandard(lm.2),
+                       leverage=hatvalues(lm.2),Cooks.Dist=cooks.distance(lm.2))
> tail(res.lev[order(cooks.distance(lm.2)),],n=8)
\end{Sinput}
\begin{Soutput}
                   Magazine     StdRes   leverage Cooks.Dist
42             Country Home  2.3636899 0.01639929 0.04657547
2  Better Homes and Gardens  1.3575202 0.05146659 0.04999601
64            Cooking Light -2.6787982 0.01432293 0.05213714
5        Sports Illustrated  2.4402089 0.02022629 0.06146309
20          Reader's Digest -1.3103123 0.06716783 0.06181267
49        AARP The Magazine -0.9279932 0.13138896 0.06513180
60               Prevention -2.6688479 0.02115025 0.07695149
4                Parade (1)  1.5938160 0.16374980 0.24870866
\end{Soutput}
\begin{Sinput}
> p <- 1 # number of predictors: x only
> c(leverage.cutoff = 2*(p+1)/dim(magdata)[1])
\end{Sinput}
\begin{Soutput}
leverage.cutoff 
     0.05714286 
\end{Soutput}
\end{Schunk}
We see that none of these cases have Cook's distances exceeding 0.50.  The points that are identified as the three most extreme outliers in the QQ plot correspond to the magazines {\em Sports Illustrated}, which overperforms expectation, and {\em Cooking Light} and {\em Prevention}, which both underperform:
\begin{Schunk}
\begin{Sinput}
> data.frame(Magazine=magdata$Magazine,Circulation=magdata$Circulation,
+            Predicted.AdRevenue=round(exp(predict(lm.2)),2),
+            Actual.AdRevenue=magdata$AdRevenue)[c(5,60,64),]
\end{Sinput}
\begin{Soutput}
             Magazine Circulation Predicted.AdRevenue Actual.AdRevenue
5  Sports Illustrated       3.205              198.46          304.185
60         Prevention       3.347              203.06          127.315
64      Cooking Light       1.717              142.68           89.153
\end{Soutput}
\end{Schunk}
Finally, here are plots of the fitted regression line, overlaid on the raw data, in the original scale and in a log-log scale.  The points colored in red correspond to the ``outlier'' magazines in the table above.
\begin{Schunk}
\begin{Sinput}
> plot(AdRevenue ~ Circulation,data=magdata)
> regline <- function(x) {exp(4.67 + 0.53*log(x))}
> curve(regline,add=T)
> points(c(3.205,3.347,1.717),c(304.185,127.315,89.153),col="red",pch="*",cex=3)
\end{Sinput}
\end{Schunk}
\begin{Schunk}
\begin{Sinput}
> plot(AdRevenue ~ Circulation,data=magdata,log="xy")
> regline <- function(x) {exp(4.67 + 0.53*log(x))}
> curve(regline,add=T)
> points(c(3.205,3.347,1.717),c(304.185,127.315,89.153),col="red",pch="*",cex=3)
\end{Sinput}
\end{Schunk}
\begin{center}
\includegraphics[width=0.4\textwidth]{hw03-solutions-regline} 
\hspace{\fill}
\includegraphics[width=0.4\textwidth]{hw03-solutions-logline}
\end{center}

\subsection*{Box-Cox transformed variables}\label{aa:box-cox}

In order to find the Box-Cox transformations, we first find the best transform for $x$, and then using the transformed $x$, we find the best transform for $y$.

First, the suggested transformation for $x$ = Circulation:
\begin{Schunk}
\begin{Sinput}
> library(car)
> with(magdata,powerTransform(Circulation~1)$roundlam)
\end{Sinput}
\begin{Soutput}
  Y1 
-0.5 
\end{Soutput}
\begin{Sinput}
> with(magdata,boxCox(Circulation~1))
\end{Sinput}
\end{Schunk}
\begin{center}
\includegraphics[width=0.5\textwidth]{hw03-solutions-boxCox-L.pdf}
\end{center}
Then, the suggested transform for $y$ = AdRevenue, when regressing on $x=1/\sqrt{\rm Circulation}$:
\begin{Schunk}
\begin{Sinput}
> lm.3 <- lm(AdRevenue ~ I(Circulation^(-0.5)),data=magdata)
> with(magdata,powerTransform(lm.3)$roundlam)
\end{Sinput}
\begin{Soutput}
Y1 
-1 
\end{Soutput}
\begin{Sinput}
> with(magdata,boxCox(lm.3))
\end{Sinput}
\end{Schunk}
\begin{center}
\includegraphics[width=0.5\textwidth]{hw03-solutions-boxCox-R.pdf}
\end{center}
So our final Box-Cox model should be $1/{\rm AdRevenue} \sim 1/\sqrt{\rm Circulation}$:
\begin{Schunk}
\begin{Sinput}
> magdata$AdRevInv <- 1/magdata$AdRevenue
> magdata$InvSqrtCirc <- 1/sqrt(magdata$Circulation)
> lm.4 <- lm(AdRevInv ~ InvSqrtCirc,data=magdata)
> # the following caused an error in R, which is why I defined the variables above...
> # lm.4 <- lm(I(AdRevenue^(-1)) ~ I(Circulation^(-0.5)),data=magdata)
> summary(lm.4)
\end{Sinput}
\begin{Soutput}
Call:
lm(formula = AdRevInv ~ InvSqrtCirc, data = magdata)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0028448 -0.0008745 -0.0000689  0.0006133  0.0040733 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.0001662  0.0004000   0.416    0.679    
InvSqrtCirc 0.0091424  0.0004571  20.000   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.001223 on 68 degrees of freedom
Multiple R-squared:  0.8547,	Adjusted R-squared:  0.8526 
F-statistic:   400 on 1 and 68 DF,  p-value: < 2.2e-16
\end{Soutput}
\begin{Sinput}
> par(mfrow=c(2,2))
> plot(lm.4)
\end{Sinput}
\end{Schunk}
\begin{center}
\includegraphics[width=0.5\textwidth]{hw03-solutions-diag-R}\label{aa:box-cox-diag}
\end{center}
This model has an $R^2=0.8547$, nearly as good as the log-log model, and again (transformed) Circulation is a strong predictor of (transformed) AdRevenue.  The residual diagnostic plots are also very comparable to the corresponding plots for the log-log model.


\subsection*{Conclusions, Simple Regression}

Here's a brief comparison of the models ($x$=Circulation, $y$=AdRevenue):

\begin{center}
\begin{tabular}{rcllp{3in}}
         &       &       & Significant & Comments on \\
Strategy & Model & $R^2$ & Predictor? & Residual Diagnostics \\
\hline
No Transform &
$y\sim x$              & 0.8945 & yes & $x$, $y$ and residuals all skewed right; some severe outliers and larger Cook's distances.\\
log-log &
$\log(y) \sim \log(x)$ & 0.881  & yes & Assumptions of normality and constant variance for residuals approximately satisfied; few outliers; no large Cook's distances.\\
Box-Cox &
$1/y \sim 1/\sqrt{x}$ & 0.8547 & yes & Similar to log-log diagnostics.
\end{tabular}
\end{center}
Since there isn't much difference in terms of fit and residual diagnostics between the log-log and Box-Cox models, we should choose based on interpretability.  The log-log model has a simpler interpreation: a 1\% change in Circulation can be expected to produce a $\hat\beta_1=0.53$\% change in AdRevenue.  Therefore we prefer the log-log model.


\section{Polynomial Regression, Untransformed Variables}\label{a:poly}
\setcounter{figure}{0}

We tried polynomial models of order 5, 4, and 3.  To save space, we just quote the $R^2$ and coefficient tables for each model:
\begin{Schunk}
\begin{Sinput}
> lm.5 <- lm(AdRevenue ~ Circulation + I(Circulation^2) + I(Circulation^3) + 
+              I(Circulation^4) + I(Circulation^5), data=magdata)
> ## Note: you could get the same model with
> ## lm.5 <- lm(AdRevenue ~ poly(Circulation, degree=5, raw=T),data=magdata)
> summary(lm.5)$r.squared
\end{Sinput}
\begin{Soutput}
[1] 0.9366882
\end{Soutput}
\begin{Sinput}
> summary(lm.5)$coefficients
\end{Sinput}
\begin{Soutput}
                      Estimate   Std. Error    t value    Pr(>|t|)
(Intercept)      57.0382867767 1.737595e+01  3.2825999 0.001668725
Circulation      47.8288287740 2.213601e+01  2.1606799 0.034469648
I(Circulation^2)  1.3621060276 7.962225e+00  0.1710710 0.864707640
I(Circulation^3) -0.6557047334 9.924238e-01 -0.6607104 0.511169347
I(Circulation^4)  0.0370875190 4.489949e-02  0.8260120 0.411865983
I(Circulation^5) -0.0005798354 6.533235e-04 -0.8875166 0.378124078
\end{Soutput}
\begin{Sinput}
> lm.6 <- lm(AdRevenue ~ Circulation + I(Circulation^2) + I(Circulation^3) + 
+              I(Circulation^4), data=magdata)
> summary(lm.6)$r.squared
\end{Sinput}
\begin{Soutput}
[1] 0.935909
\end{Soutput}
\begin{Sinput}
> summary(lm.6)$coefficients
\end{Sinput}
\begin{Soutput}
                     Estimate   Std. Error   t value     Pr(>|t|)
(Intercept)      45.686665431 11.742208057  3.890807 2.377459e-04
Circulation      65.380572047  9.928749452  6.584976 9.325058e-09
I(Circulation^2) -5.499586772  1.900414956 -2.893887 5.173919e-03
I(Circulation^3)  0.220172602  0.104527351  2.106364 3.903839e-02
I(Circulation^4) -0.002733043  0.001694505 -1.612886 1.116144e-01
\end{Soutput}
\begin{Sinput}
> lm.7 <- lm(AdRevenue ~ Circulation + I(Circulation^2) + I(Circulation^3),
+            data=magdata)
> summary(lm.7)$r.squared
\end{Sinput}
\begin{Soutput}
[1] 0.933344
\end{Soutput}
\begin{Sinput}
> summary(lm.7)$coefficients
\end{Sinput}
\begin{Soutput}
                    Estimate  Std. Error   t value     Pr(>|t|)
(Intercept)      59.17036829 8.345045881  7.090478 1.118099e-09
Circulation      51.23581639 4.711234296 10.875243 2.334496e-16
I(Circulation^2) -2.50537894 0.411411261 -6.089719 6.476556e-08
I(Circulation^3)  0.05222479 0.009229702  5.658339 3.574381e-07
\end{Soutput}
\end{Schunk}
All three models have $R^2\approx 0.93$; in the models of order 4 and 5, the predictors (Circulation)$^4$ and (Circulation)$^5$ were not significant predictors; all predictors were significant for the model of order 3.  Therefore we concentrate on the model of order 3.  Here is the complete summary and residual diagnostics for that model:
\begin{Schunk}
\begin{Sinput}
> summary(lm.7)
\end{Sinput}
\begin{Soutput}
Call:
lm(formula = AdRevenue ~ Circulation + I(Circulation^2) + I(Circulation^3), 
    data = magdata)

Residuals:
   Min     1Q Median     3Q    Max 
-83.75 -13.56  -2.16  11.46 104.82 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      59.17037    8.34505   7.090 1.12e-09 ***
Circulation      51.23582    4.71123  10.875 2.33e-16 ***
I(Circulation^2) -2.50538    0.41141  -6.090 6.48e-08 ***
I(Circulation^3)  0.05223    0.00923   5.658 3.57e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 34.06 on 66 degrees of freedom
Multiple R-squared:  0.9333,	Adjusted R-squared:  0.9303 
F-statistic: 308.1 on 3 and 66 DF,  p-value: < 2.2e-16
\end{Soutput}
\begin{Sinput}
> par(mfrow=c(2,2))
> plot(lm.7)
\end{Sinput}
\end{Schunk}
\begin{center}
\includegraphics[width=0.5\textwidth]{hw03-solutions-polydiag}\label{aa:poly-diag}
\end{center}
This set of residual diagnostic plots does not look as good as either the log-log or Box-Cox residual diagnostics in Section \ref{a:log-log}.  Although $R^2$ is higher, we seem to be farther from the assumptions underlying regression here.

\section{Final Model and Predictions}\label{a:final}
\setcounter{figure}{0}

Comparing the regression output and residual diagnostic plots for the log-log model and polynomial model of order 3: 
\begin{itemize}
\item The log-log model has dramatically better residual diagnostic plots;
\item The log-log model has a significant predictor and an $R^2$ nearly as high as the polynomial model;
\item The log-log model has a simple interpretation: for every 1\% increase in Circulation, we can expect a 0.53\% increase in Ad Revenue.
\end{itemize}
For these reasons, we prefer the log-log model.\\

Here are 95\% prediction intervals for Ad Revenue, for a publication with circulation of 0.5 million and 20 million, respectively, from the log-log model. Note that we have to exponentiate the endpoints of the intervals, to ``undo'' the log transformation on AdRevenue.
\begin{Schunk}
\begin{Sinput}
> int.A.i <- round(exp(predict(lm.2,newdata=data.frame(AdRevenue=0,Circulation=0.5),
+                         interval="prediction")[c(2,3)]),2)
> int.A.ii <- round(exp(predict(lm.2,newdata=data.frame(AdRevenue=0,Circulation=20),
+                         interval="prediction")[c(2,3)]),2)
> data.frame("Lower Endpoint"=c("Circulation 0.5 Million"=int.A.i[1],
+                               "Circulation 20 Million"=int.A.ii[1]), 
+            "Upper Endpoint"=c("Circulation 0.5 Million"=int.A.i[2],
+                               "Circulation 20 Million"=int.A.ii[2]))
\end{Sinput}
\begin{Soutput}
                        Lower.Endpoint Upper.Endpoint
Circulation 0.5 Million          51.82         106.55
Circulation 20 Million          359.90         758.76
\end{Soutput}
\end{Schunk}

\end{document}
