\documentclass[dvipsnames, pdflatex,slidecentered]{beamer} %add red inside []
\usepackage{xcolor,colortbl}
\usepackage{amsmath,epsfig,latexsym}
\usepackage{hyperref}

%  options  [handout]

%  There is a VERY rich set of possible
%  styles of presentations and color themes.  See the
%  beamer documentation for a full list of possibilities.
\usetheme{default}  % JuanLesPins  Rochester
                   %  Berkeley  Palo Alto    with sidebar and top
                   % Goettingen  Marburg     with sidebar
                   %  Copenhagen  Luebeck  Warsaw
                   % Antibes with overlays
\useinnertheme{rounded}
\usecolortheme{rose} %crane %lily
%\usefonttheme{structuresmallcapsserif}

\usepackage{pgfarrows}
%  Don't need to load the pgf package, but it has
%  has itself some packages you might want, such as
%  pgfarrows,,pgfnodes,pgfautomata,pgfheaps
%  See the pgf documentation.


\usepackage{amsmath,amssymb}

% \beamertemplatetransparentcovereddynamic
                % overlays that are upcoming are transparent, in
                % a manner that depends upon how far ahead  they are.



%\usepackage{SweaveSlides}
\title[Stat 849'F06]{Stat 849: Linear regression diagnostics}
\author[S. Kele\c{s}]{S\"und\"uz Kele\c{s}}
\institute[UW]{Department of Statistics\\
 Department of Biostatistics and Medical Informatics\\
University of Wisconsin, Madison}
\date[]{}


%\mode<presentation>
%{ \usetheme{boxes} }



\begin{document}

\frame{\titlepage}

\begin{frame}
\frametitle{Regression Diagnostics}

\begin{itemize}
\item Checking the Gauss-Markov assumptions
on the error vector using graphical and numerical methods.

\item Checking the normality assumption.

\item Weighted least squares: Heteroscedascity.

\item  Outlier and influential point detection.

\end{itemize}

\end{frame}


\begin{frame}
\frametitle{Assumptions on the error vector}
\begin{eqnarray*}
E[\epsilon] & = & 0\\
var(\epsilon) & =  & \sigma^2 I_n \\
\epsilon &\sim& \mathcal{N}(0, \sigma^2 I_n)
\end{eqnarray*}

\pause  First two are the G-M assumptions under which $\hat{\beta}$
is BLUE.

\end{frame}

\begin{frame}
\frametitle{Residuals}
\[\hat{\epsilon}_i=y_i-\hat{y}_i, \quad i=1, \cdots, n.\]
Residuals are useful for answering questions such as

\begin{itemize}
\item <+-> The regression function is not linear.

\item <+-> The error terms do not have a constant variance.

\item <+-> The model fits pretty well but there are a few outliers.

\item <+-> The error terms are not normally distributed.

\item <+-> One or several important predictors have been omitted from
the model.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Informal diagnostic plots of residuals}

\begin{itemize}
\item <+-> residuals vs predictor variable,

\item <+-> absolute or squared residuals vs predictor variable,

\item <+-> residuals vs fitted values,

\item <+-> residuals vs time or other sequence,

\item <+-> residuals vs omitted predictor variables,

\item <+-> box plot of residuals,

\item <+-> normal probability plot (quantile plot) of residuals.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Nonlinearity of the regression function}
\begin{figure}
\centering
\includegraphics[scale=0.3]{resid_nonlin}
\caption{\textit{Relation between maps distributed and bus
ridership in eight test cities.} Whether a linear regression
function is appropriate for the data being analyzed can be studied
from a \textit{residual plot against the predictor variable} or
equivalently from a \textit{residual plot against the fitted
values}.}
\end{figure}

\end{frame}

\begin{frame}
\frametitle{Checking $var(\epsilon)  =   \sigma^2 I_n$}
\begin{itemize}
\item Plot residuals $\hat{\epsilon}$ versus fitted values
$\hat{y}$ \textit{or} residuals $\hat{\epsilon}$ versus predictor
variable $X$.

\end{itemize}

\begin{figure}
\centering
\includegraphics[scale=0.3]{resid_hete}
\caption{Heteroscasdicity?.}
\end{figure}
\end{frame}


\begin{frame}
\frametitle{Checking $var(\epsilon)  =   \sigma^2 I_n$}
\begin{itemize}

\item Plotting the absolute values of the residuals or of the
squared residuals against the predictor variable $X$ or against the
fitted values $\hat{y}$ are also useful for diagnosing nonconstancy
of the error variance. Especially useful when there are not many
observations in the data set.
\end{itemize}

\end{frame}


\begin{frame}

\frametitle{How can we deal with these?}
\begin{itemize}
\item <+-> Weighted least squares (when the linearity assumption is ok, but the only concern is nonconstant
variance).

\item <+-> Transformation of the variables, i.e., response or the
predictors (nonlinearity and nonconstant variance might go
together).

\item <+-> Changing the linear model.
\end{itemize}
\end{frame}


\begin{frame}
\frametitle{Weighted Least Squares}

\begin{itemize}
\item <+-> Although the mean of the dependent variable might be a linear
function of the regressors, the variance of the error terms might
also depend on those same regressors.

\item <+-> Heteroscedasticity is often seen with \textit{aggregate data}.

\item <+-> Assume we are collecting data from $n$ groups, e.g., different
days, regions etc...

\item <+-> For individual $i$ in group $j$, we have the following linear
model

\[Y_{ij}=\beta_0+\beta_1X_{ij}+\epsilon_{ij}, \quad \mbox{var}(\epsilon_{ij})=\sigma^2.\]

\pause Now, assume  that we only observe
\textcolor{red}{\textit{group averages}} (Can you think of
examples??):

\[\bar{Y}_j=\frac{1}{n_j} \sum_{i=1}^{n_j} Y_{ij}, \quad
\bar{X}_j=\frac{1}{n_j} \sum_{i=1}^{n_j} X_{ij}.\]

\pause  Then,
\[\bar{Y}_j=\beta_0+\beta_1\bar{X}_j +\bar{\epsilon}_j, \quad
\mbox{var}(\bar{\epsilon}_j)=\sigma^2 \left ( \frac{1}{n_j} \right
)=\sigma^2 h_j.\]

\end{itemize}

\end{frame}

\begin{frame}
\frametitle{Apple trees} [Weisberg, 1985] Apple trees produce ``long
shoots'' (which may grow as much as 15 to 20 cm over a growing
season) as well as ``short shoots''. Samples of both of these types
of shoots were taken from trees in an orchard every few days during
a growing season (106 days), and they were then analyzed in a
laboratory. We will look at only long shoots. Among the measurements
taken was a count of the number of ``stem units'' on each shoot.
\pause

\begin{itemize}
\item day of the growing season ($\bar{X}_j$),

\item \# shoots sampled that day ($n_j$),

\item mean number of stem units among the sampled shoots
($\bar{Y}_j=1/n_j\sum_{i=1}^{n_j}Y_{ij}$),

\item standard deviation of the stem unit counts ($\hat{\sigma}
\sqrt{n_j}$).
\end{itemize}

\end{frame}

\begin{frame}
\textcolor{red}{Interested in modelling}
\[ \mbox{mean \# stem units} \sim \mbox{day of the growing
season}.\]

\textcolor{red}{General linear model:}

\[Y=X\beta+\epsilon, \quad E[\epsilon]=0, \quad
cov(\epsilon)=\sigma^2V=\Sigma.\]

This no longer satisfies the G-M conditions.

\begin{itemize}
\item <+-> What is the least squares estimator for $\beta$?
 \pause \[(X^TX)^{-1}X^TY\]

\item <+->  What happens if we blindly use
$\hat{\beta}_{OLS}=(X'X)^{-1}X'Y$ as an estimator of $\beta$?

\item <+-> Statistical properties to check: \textit{unbiasedness},
\textit{efficiency}!

\item <+->  Can we find the unique BLUE for $\beta$ (when $X$ is full
rank)?

\end{itemize}
\pause \textbf{Consider transforming this linear model into a form
that we know about!}

\end{frame}


\begin{frame}
\begin{eqnarray*}
V &=& KK^T \quad \mbox{(Cholesky decomposition)}\\
Y^*&=&K^{-1}Y, \quad X^*=K^{-1}X\\
\end{eqnarray*}
Now multiply $Y=X\beta+\epsilon$, $E[\epsilon=0]$,
$\mbox{var}(\epsilon)=\sigma^2V$.
\begin{eqnarray*}
\underbrace{K^{-1}Y}_{Y^*}&=&\underbrace{K^{-1}X}_{X^*}\beta+K^{-1}+K^{-1}\epsilon, \quad \mbox{var}(K^{-1}\epsilon)=\sigma^2I_n\\
Y^* & = & X^* \beta+\epsilon^* \quad \mbox{G-M satisfied!}\\
0 & = & \frac{\partial}{\partial \beta}(Y^*-X^*\beta)^T(Y^*-X^*\beta)\\
&\Longrightarrow& \hat{\beta}=((X^*)^TX^*)^{-1}(X^*)^TY^*=(X^TV^{-1}X)^{-1}X^TV^{-1}Y\\
&\Longrightarrow& var(\hat{\beta})=\sigma^2 (X^tV^{-1}X)^{-1}.
\end{eqnarray*}
$\hat{\beta}$ is called the \textit{generalized least squares
estimator}, and it is called \textit{weighted least squares
estimator} when $V$ is diagonal.

\begin{center}
\textcolor{red}{How about the statistical properties of
$\hat{\beta}$?}
\end{center}

\end{frame}


\begin{frame}
\frametitle{Properties of the GLE $\hat{\beta}$?}

GLE $\hat{\beta}$ is simply the ordinary least squares estimator
for the transformed model, we would expect $\hat{\beta}$ to have
the same statistical properties of the ordinary least squares
estimator in this transformed model.

\begin{itemize}
\item  $E[\hat{\beta}] = \beta$ (OLS for the transformed model) and
$E[\hat{\beta}_{OLS}]=\beta$ (OLS for the original model).

\item $var(a'\hat{\beta}) \le var(a'\hat{\beta}_{OLS})$,  $\forall
a_{p \times 1}$ and $var(a'\hat{\beta}) \le var(a'\tilde{\beta})$
where $\tilde{\beta}$ is any unbiased estimator in the general
linear model.
\end{itemize}

So, $\hat{\beta}_{OLS}$ is \textcolor{red}{inefficient} but
\textcolor{red}{unbiased}. Can we use it for inference?

\end{frame}

\begin{frame}
\frametitle{Apple trees data}
\begin{figure}
\begin{center}
\includegraphics[scale = 0.4] {appleshoots_pairplot.pdf}
\end{center}
\end{figure}
\end{frame}


\begin{frame}[fragile]
\frametitle{R script for the apple trees example} \small
%sshoots = read.table("sshoots", header = F)
%lshoots = read.table("lshoots", header = F)
% all = rbind(sshoots, lshoots)
\begin{verbatim}
 all = read.table("apple_shoots.txt", header = F)
 names(all) = c("day", "total", "mc", "std")
 attach(all)
 lm.wls = lm(mc ~ day, data = all, weights = total)
 lm.ols = lm(mc ~ day, data = all)
 ? lm
 ...
  weights: an optional vector of weights to be used
          in the fitting process.  Should be 'NULL' or a
          numeric vector. If non-NULL, weighted least squares
          is used with weights 'weights' (that is, minimizing
          'sum(w*e^2)'); otherwise ordinary least
          squares is used.
\end{verbatim}
\end{frame}

\begin{frame}[fragile]
\footnotesize
\begin{verbatim}
 > summary(lm.ols)$coef
              Estimate  Std. Error  t value     Pr(>|t|)
 (Intercept) 9.7688811 0.337028691 28.98531 6.310058e-33
  day        0.1963376 0.005648492 34.75930 1.094398e-36

#Utilize total number of counts as weights.
 > summary(lm.wls)$coef
               Estimate  Std. Error  t value     Pr(>|t|)
 (Intercept) 10.0354917 0.414189300 24.22924 2.736109e-29
  day         0.1909468 0.006364845 30.00023 1.237558e-33

#Utilize 1/(standard errors^2) as weights
 lm.wls1 = lm(mc ~ day, data = all, weights = std^-2)
 > summary(lm.wls1)$coef
             Estimate  Std. Error  t value     Pr(>|t|)
 (Intercept) 9.8150263 0.271764428 36.11593 5.655701e-37
  day        0.1894794 0.004640411 40.83246 1.671095e-39
\end{verbatim}
\end{frame}


\begin{frame}[fragile]
\frametitle{Diagnostic plots for the apple trees example}
\vspace{-1.5cm}
\begin{columns}
\onslide  <+->
\begin{column}{4.5cm}
\begin{figure}
\begin{center}
\includegraphics[scale = 0.35]{lm_ols.pdf}
\caption{Ordinary least squares fit.}
\end{center}
\end{figure}
\end{column}
\onslide  <+->
\begin{column}{4.5cm}
\begin{figure}
\begin{center}
\includegraphics[scale = 0.35]{lm_wls_n.pdf}
\caption{Weighted least squares fit using $n$ as weights.}
\end{center}
\end{figure}
\end{column}
\end{columns}

\textcolor{red}{No substantial differences between the two fits!}
\end{frame}


\begin{frame}[fragile]
\frametitle{Diagnostic plots for the apple trees example}
\vspace{-1.5cm}
\begin{columns}
\onslide  <+->
\begin{column}{4.5cm}
\begin{figure}
\begin{center}
\includegraphics[scale = 0.35]{std_vs_mc}
\caption{Mean counts vs standard errors.}
\end{center}
\end{figure}
\end{column}
\onslide  <+->
\begin{column}{4.5cm}
\begin{figure}
\begin{center}
\includegraphics[scale = 0.35]{lm_wls_n.pdf}
\caption{Weighted least squares fit using
$\frac{1}{\hat{\sigma}_i^2}$ as weights.}
\end{center}
\end{figure}
\end{column}
\end{columns}

\pause \textcolor{red}{Which of three models is preferable?} \pause
\small
\begin{verbatim}
 > c(summary(lm.ols)$r.squared, summary(lm.wls)$r.squared,
     summary(lm.wls1)$r.squared)
 [1] 0.9602610 0.9473692 0.9714500
\end{verbatim}


\end{frame}

\begin{frame}
\frametitle{Checking normality of the error terms: $\epsilon \sim
\mathcal{N}(0, \sigma^2I_n)$}

\begin{itemize}

\item We observe estimates $\hat{\epsilon}$ of the error terms.

\item We can use these estimates to compare observed values to the theoretical
values.


\item This can be  done via a normal  probability plot or a
quantile-quantile (Q-Q) plot.

\item Each residual is plotted against its expected value under
normality (comparing empirical cumulative distribution with the
theoretical one). \textcolor{red}{A plot  that departs substantially
from linearity provides evidence against  normality of the error
distribution.}

\end{itemize}

\end{frame}

\begin{frame}
\frametitle{Q-Q plots}
 Quantile point $q_p$ for random variable $X$ is the
point such that
\[F_X(q_p)=P(X \le q_p)=p, \quad q_p=F^{-1}(p).\]

If $q_p^X$ and $q_p^Y$ are quantile functions of random variables
$X$ and $Y$, Q-Q plot of $X$ and $Y$ is the plots of $(q_p^X,
q_p^Y)$ for all $p$.


%If we standardize the residuals (dividing by $\hat{\sigma}$, then
%their theoretical distribution is standard normal.
\begin{enumerate}

\item Sort the residuals $\hat{\epsilon}_{(1)} \le
\hat{\epsilon}_{(2)}\le \cdots, \le \hat{\epsilon}_{(n)}$.

 \item
Compute $u_i = \Psi^{-1}\left ( \frac{i}{n} \right )$, where
$\Psi$ is the probability distribution function of the standard
normal distribution.

 \item Plot $\hat{\epsilon}_{(i)}$ against
$u_i$.
\end{enumerate}


%\begin{table}
%\begin{center}
%\begin{tabular}{llll}
% & (1) & (2) & (3) \\
% &&&Expected\\
% Run & Residual & Rank & value under \\
% $i$ & $\hat{\epsilon}_i$ & k& normality\\
% \hline
% 1 & 51.02 & 22 & 51.95\\
% 2 & -48.47 & 5 & -44.10 \\
% 3 & -19.88 & 10 & -14.76 \\
% $\cdots$ & $\cdots$ & $\cdots$ & $\cdots$ \\
% 23 & 38.83 & 19 &  31.05\\
% 24 & -5.98 & 13 & 0\\
% 25 & 10.72 & 17 & 19.93\\
%\end{tabular}
%\end{center}
%\end{table}


\end{frame}

\begin{frame}[fragile]
\frametitle{Example} \footnotesize
\begin{verbatim}
 set.seed(1)
 n = 300
 x = rnorm(n, 1, 3)
 y = 0.5 + 4*x + rnorm(n, 0, 1)
 lm1 = lm(y ~ x)
 > summary(lm1)
 Call: lm(formula = y ~ x)
 Residuals:
      Min       1Q   Median       3Q      Max
 -2.97733 -0.67295 -0.02005  0.70834  3.85267

 Coefficients:
             Estimate Std. Error t value Pr(>|t|)
 (Intercept)  0.47995    0.06459   7.431 1.15e-12 ***
               4.00851    0.02091 191.708  < 2e-16 ***
 ---
 Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

 Residual standard error: 1.045 on 298 degrees of freedom Multiple
 R-Squared: 0.992,      Adjusted R-squared: 0.9919 F-statistic:
 3.675e+04 on 1 and 298 DF,  p-value: < 2.2e-16
\end{verbatim}
\end{frame}

\begin{frame}[fragile]
\frametitle{Example} \footnotesize
\begin{verbatim}
 sigma = sqrt(sum(lm1$resid^2)/(n-2))
 > sigma
 [1] 1.045298

 sresid = sort(lm1$resid)
 tsresid1 = qnorm(c(1 : n) / n)
 postscript("qqplot.eps")
 qqnorm(lm1$resid)
 points(tsresid1, sresid, pch = 3, col = "red")
 qqline(lm1$resid)
 dev.off()
\end{verbatim}
\end{frame}

\begin{frame}
\begin{columns}
\begin{column}{4.5cm}
\begin{figure}
\centering \includegraphics[angle = 0, scale = 0.2]{qqplot}
\caption{$\hat{\epsilon}_{(i)}$ versus \textcolor{red}{$u_i =
\Psi^{-1}\left( \frac{i}{n}\right)$}.}
\end{figure}
\end{column}
\begin{column}{4.5cm}
\begin{figure}
\centering \includegraphics[angle = 0, scale = 0.2]{hist}
\caption{Histogram of residuals is also useful for checking
normality.}
\end{figure}
\end{column}
\end{columns}
\end{frame}

\begin{frame}
\begin{figure}
\centering \includegraphics[angle = 0, scale = 0.4]{nonnormal}
\caption{Examples of deviation from normality.}
\end{figure}

Possible solutions: Transform $Y$, i.e., use $\log{Y}$, $\exp{Y}$ as
the response variable.
\end{frame}

%\begin{frame}
%\begin{figure}
%\centering \includegraphics[angle = 0, scale = 0.3]{manyqqplot}
%\caption{}
%\end{figure}
%\end{frame}


\begin{frame}
\frametitle{Are the errors independent?}

In general no way to check this assumption from the plots alone,
need to know how the data were collected.

%If the data were gathered in certain order, could plot residuals
%against this order (e.g., different times of the day, year etc...)

%If parts of data were collected by different people, could make
%separate box-plots of the residuals for each person to see if
%there are some relationship in the errors based on the person who
%gathered the data.

Whenever the data is obtained in a time sequence or some other
type of sequence, such as for adjacent geographic areas, it is
useful to prepare a \textcolor{red}{sequence plot of the
residuals}.

\pause
\begin{figure}
\centering \includegraphics[angle = 0, scale = 0.3]{nonindepres}
\caption{Nonindependence of error terms.}
\end{figure}

\end{frame}


\begin{frame}[fragile]
\frametitle{Omission of important predictor variables}

Residuals should also be plotted against variables omitted from
the model that might have important effects on the response.

\pause
\textbf{E.g.}

 Consider three explanatory variables $X_1$, $X_2$, $X_3$ where $X_2$ and $X_3$ are binary
variables.

 True model $Y \sim X_1 + X_2$.

 Fitted model $Y \sim X_1$.
 \footnotesize

\pause
\begin{verbatim}
X1 = runif(n, 20, 45)

X2 = rbinom(n, c(0, 1), 0.7)

X3 = rbinom(n, c(0, 1), 0.5)

Y = 0.4 + 0.1 * X1 + 2 * X2 + rnorm(n, 0, 0.5)
\end{verbatim}

\pause \normalsize Plot residuals from the fitted model $Y \sim X_1$
versus $X_2$ and $X_3$ to see whether either of these should be
included in the
model.\\
\textcolor{red}{Any ideas how these plots might look like?}
\end{frame}

\begin{frame}
\begin{figure}
\centering \includegraphics[angle = 90, scale = 0.3]{omittedvar}
\caption{\textit{Omission of important predictor variables}.
Residuals from the fitted model $Y \sim X_1$ versus $X_1$, $X_2$,
and $X_3$.}
\end{figure}

\end{frame}

\begin{frame}
\begin{figure}
\centering \includegraphics[angle = 90, scale = 0.25]{omitted}
\caption{\textit{Omission of important predictor variables}. $Y$:
blood pressure in the morning. $X_1$:  Calories at dinner. $X_2$:
Age range.}
\end{figure}
\end{frame}

\begin{frame}
\frametitle{Omission of important predictor variables}

Does not mean the original model is \textit{wrong}.

Implies that it can be improved by adding another predictor
variable.
\end{frame}


\begin{frame}
\frametitle{Added-variable plots}
\begin{itemize}
\item Residual plots may not properly show the nature of the
marginal effect of a predictor variable, given the other predictor
variables in a model.

\item \textit{\textcolor{red}{Added-variable plots}}, also called
\textit{partial regression plots} or \textit{adjusted variable
plots}, are refined residual plots that provide graphic
information about the marginal importance of a predictor variable
$X_j$, given the other predictor variables already in the model.

\item In addition, these plot might be useful for identifying the
nature of the marginal relation for a predictor variable in the
regression model.

\end{itemize}
\end{frame}

\begin{frame}
\frametitle{How to obtain added-variable plots?}
\begin{itemize}
\item Both the response variable $Y$ and the predictor variable
$X_j$ under consideration are regressed against the other
predicted variables in the model and the residuals are obtained
for each.

\item These residuals reflect the part of each variable that is
\textit{not linearly associated} with the other predicted
variables already in the regression model.


\item The plot of these residuals against each other (1) shows the
marginal importance of this variable in reducing the residual
variability and (2) may provide information about the nature of
the marginal regression relation for the predictor variable $X_j$
under consideration for possible inclusion in the regression
model.

\end{itemize}

\end{frame}

\begin{frame}[fragile]
Consider $Y, X_1, X_2$.

\begin{eqnarray*}
\hat{Y}(X_1) &=& \hat{\beta}_0 +\hat{\beta}_1X_1\\
\hat{\epsilon}(Y \mid X_1) &=& Y -\hat{Y}(X_1)\\
\hat{X}_2(X_1)&=&\hat{\beta}_0^* +\hat{\beta}_1^*X_1\\
\hat{\epsilon}(X_2 \mid X_1) &=& X_2 -\hat{X}_2(X_1)\\
\end{eqnarray*}

The added-variable plot for predictor variable $X_2$ consists of a
plot of the $Y$ residuals $\hat{\epsilon}(Y \mid X_1)$ against the
$X_2$ residuals $\hat{\epsilon}(X_2 \mid X_1)$.

\end{frame}

\begin{frame}
\begin{figure}
\centering \includegraphics[angle = 0, scale = 0.4]{addvar1}
\caption{\textit{Added-variable plot}. True model: $Y \sim
X_1+X_2$.}
\end{figure}
\end{frame}

\begin{frame}
\begin{figure}
\centering \includegraphics[angle = 0, scale = 0.4]{addvar2}
\caption{\textit{Added-variable plot}. True model: $Y \sim
X_1+X_2^2$.}
\end{figure}
\end{frame}

\begin{frame}
\begin{figure}
\centering \includegraphics[angle = 0, scale = 0.35]{addvar3}
\caption{\textit{Added-variable plot}. Question: How does
$\hat{\beta}_1$ from $Y \sim  X_1 + X_2$ relate to the slope of the
regression line on the right panel?}
\end{figure}

\end{frame}

\begin{frame}[fragile]
\footnotesize
\begin{verbatim}
lm0 = lm(Y ~ X1 + X2)

lm1 = lm(Y ~ X2)

lm2 = lm(X1 ~ X2)

 > lm0
 Call: lm(formula = Y ~ X1 + X2)
 Coefficients:
 (Intercept)           X1       X2
    -205.719        6.288    4.738

 > lm(lm1$resid ~ lm2$resid)
 Call: lm(formula = lm1$resid ~ lm2$resid)
  Coefficients: (Intercept)    lm2$resid
                -2.458e-15    6.288e+00
\end{verbatim}
\normalsize
\textbf{Exercise:} Can you argue this analytically?
%(\textit{Hint:} Think about orthogonality.)
\end{frame}

\begin{frame}
\frametitle{Notes on added-variable plots}

\begin{itemize}
\item An added variable plot only suggests the nature of the
functional relation in which a predictor variable should be added
to the regression model but does not provide an analytical
expression.

\item The relation shown is for predictor $X_j$ adjusted for the
other predictor variables in the regression model, not for $X_j$
directly.

\item These plots may not show the proper form of the marginal
effect of a predictor variable if the functional relations for some
or all of the variables already in the regression model are
misspecified.

\end{itemize}

\end{frame}



\begin{frame}
\frametitle{Outliers and Influential observations}

\begin{center}
\begin{itemize}
\item Regression outliers.

\item High leverage points.

\item Influential points.
\end{itemize}
\textcolor{red}{What are these and how can we measure them?}
\end{center}
\pause One should be cautious about unusual data points in linear
models since they can influence the results of the analysis, and
their presence might may be a signal that the model fails to capture
important characteristics of the data.


\end{frame}


\begin{frame}
\frametitle{Outliers}

\begin{itemize}
\item Outliers are extreme observations.


\item It is common practice to distinguish between two types of
outliers.


\item Outliers in the response variable  are called
\textcolor{red}{residual outliers.}

\item Outliers with respect to the predictors are called
\textcolor{red}{leverage points.} They can affect the regression
model, too. Their response variables need not be outliers. However,
they may almost uniquely determine regression coefficients. They may
also cause the standard errors of regression coefficients to be much
smaller than they would be if the observation were excluded.

\end{itemize}

\end{frame}

\begin{frame}
\frametitle{Outlier detection}

\begin{itemize}
\item  Residual outliers can be identified from (1) residual plots
against $X$ or $\hat{Y}$; (2) box plots, stem-and-leaf plots, and
dot plots of the residuals.


\item Plotting standardized residuals might be  helpful.

\item Standardized residuals:
$\hat{\epsilon_i}/\sqrt{\mbox{var}(\hat{\epsilon}_i)}$.

\[\mbox{var}(\hat{\epsilon})=?\]

\end{itemize}

\end{frame}


\begin{frame}
\frametitle{Standardized (internally studentized) residuals}
\begin{eqnarray*}
\mbox{var}(\hat{\epsilon}) &=& \mbox{var}(Y-\hat{Y})\\
&=& \mbox{var}(Y-PY)= \sigma^2(I-P),
\end{eqnarray*}
where $P=X(X^TX)^{-1}X^T$ is the projection matrix. It is also
called the \textcolor{red}{\textit{hat matrix H}}.

Then, $\mbox{var}(\hat{\epsilon}_i)=\sigma^2(1-h_{ii})$ where
$h_{ii}$ is the $i$-th element on the main diagonal of the hat
matrix, and the covariance between residuals $\hat{\epsilon}_i$
and $\hat{\epsilon}_j$ $(i \neq j)$ is
$\sigma^2(0-h_{ij})=-h_{ij}\sigma^2$, where $h_{ij}$ is the
element in the $i$-th row and $j$-th column of the hat matrix.


\textcolor{red}{Note:} This can only correct for the natural
non-constant variance in residuals when errors $\epsilon_i$ have
constant variance.
\end{frame}

\begin{frame}
\frametitle{Leverage}
\begin{itemize}
\item  $h_{ii}$ is called the \textcolor{red}{leverage} (in terms of
the $X$ values) of the $i$-th case.

\item Properties of $h_{ii}$: $0 \le h_{ii} \le 1$, $\sum_{i=1}^n
h_{ii} = p$.

\item The $h_{ii}$ values theoretically range from $1/n$ to $1$. Those
that exceed $2p/n$ are said to be large.

\item \textcolor{red}{Interpretation:} It is a measure of distance
between the $X$ values for the $i$-th case and the means of the $X$
values for all of the observations. It measures the degree of
conformity of a single observation to the linear pattern established
by the other $n-1$ observations.


\item  For a linear regression with model with one predictor, the
leverage associated with the specific observation $(x_i, y_i)$ is
\[h_{ii} =  \frac{1}{n}+
\frac{(x_i-\bar{x})^2}{\sum_{i=1}^n(x_i-\bar{x})^2}.\]

\end{itemize}

\end{frame}


%\begin{frame}

%All observations are not equally important in determining the
%estimated value $\hat{y}_i$.

%Leverage
%\end{frame}


\begin{frame}[fragile]
\small
\begin{verbatim}
 #Computing standardized  residuals
 xx = cbind(1, x)
 H = xx %*% solve((t(xx)%*%xx)) %*% t(xx) #Hat Matrix
 dhat = (1-diag(H))
 #Could also use lm.influence(lm1)$hat
 lm1.stdres0 = lm1$resid/sqrt((sum(lm1$resid^2)/(n-2 ))*dhat)
 #Two R functions that will give us standardized residuals
 lm1.stdres2 = rstandard(lm1)
 library(MASS)
 lm1.stdres1 = stdres(lm1)
\end{verbatim}

\end{frame}

\begin{frame}

\begin{figure}
\centering \includegraphics[angle = 0, scale = 0.3]{normqqplot}
\caption{\textit{Outlier detection.}}
\end{figure}

\end{frame}

\begin{frame}
\begin{figure}
\centering \includegraphics[angle = 90, scale = 0.3]{outlier}
\caption{\textit{Outlier detection}. \footnotesize What to do with
outliers? Should only discard these if there is a direct evidence
that it represents an error in recoding, a miscalculation, a
malfunctioning of equipment, or a similar type of circumstance.
Sometimes detection of an outlier itself might be of interest.}
\end{figure}



\end{frame}

\begin{frame}[fragile]
\footnotesize
\begin{verbatim}
#Model fit with all observations:
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.47995    0.06459   7.431 1.15e-12 *** x
             4.00851    0.02091 191.708  < 2e-16 ***
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

Residual standard error: 1.045 on 298 degrees of freedom Multiple
R-Squared: 0.992,      Adjusted R-squared: 0.9919 F-statistic:
3.675e+04 on 1 and 298 DF,  p-value: < 2.2e-16

#Model fit omitting the circled outlier point:
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.46061    0.06340   7.265  3.3e-12 *** x[-195]
             4.01431    0.02051 195.681  < 2e-16 ***
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

Residual standard error: 1.023 on 297 degrees of freedom Multiple
R-Squared: 0.9923,     Adjusted R-squared: 0.9923 F-statistic:
3.829e+04 on 1 and 297 DF,  p-value: < 2.2e-16

#Model fit omitting the highest leverage point:
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.47060    0.06499   7.241 3.83e-12 *** x[-232]
             4.01307    0.02122 189.098  < 2e-16 ***
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

Residual standard error: 1.044 on 297 degrees of freedom Multiple
R-Squared: 0.9918,     Adjusted R-squared: 0.9917 F-statistic:
3.576e+04 on 1 and 297 DF,  p-value: < 2.2e-16

\end{verbatim}

\end{frame}

\end{document}
