%% modified from 201A-f08/lecture/regression.tex  for 201B 
%\documentclass[11pt]{article}
\documentclass[12pt]{article}   %% use 12pt \large for lectures
\usepackage{graphicx}
\usepackage{amssymb}

\textwidth = 6.5 in
\textheight = 9 in
\oddsidemargin = 0.0 in
\evensidemargin = 0.0 in
\topmargin = 0.0 in
\headheight = 0.0 in
\headsep = 0.0 in
\parskip = 2pt
%\parindent = 2em

\newcommand{\be}{\begin{enumerate}}
\newcommand{\ee}{\end{enumerate}}
\newcommand{\bi}{\begin{itemize}}
\newcommand{\ei}{\end{itemize}}
\newcommand{\bd}{\begin{description}}
\newcommand{\ed}{\end{description}}
\newcommand{\para}{\medskip\noindent\textbf}

\newcommand{\var}{\mbox{\rm var}}
\newcommand{\cov}{\mbox{\rm cov}}
\newcommand{\logit}{\mbox{\rm logit}}
\newcommand{\ilogit}{\mbox{\rm ilogit}}

\newcommand{\graph}[1]{\includegraphics{../graph/#1}}

\begin{document}

 \begin{flushleft}
 {\bf Stats 201B  (W14)} Regression Analysis: Model Building, Fitting, and Criticism %\hfill Hongquan Xu \\
\end{flushleft}

{\small
\bi
%\item  Textbook: Chapter 10 of  Montgomery (2005).   ``Design and Analysis of Experiments,"  Wiley.
\item  Textbook: J. J. Faraway (2005). ``Linear Models with R," Chapman \& Hall.
%\item  J. J. Faraway (2006). ``Extending the Linear Model with R," Chapman \& Hall.
\ei
}

\LARGE \Large % \large
\section*{Chapter 1. Introduction}

\para{Before you start}
\bi
\item Look before you leap!
\item Understand the physical problem and objective
\item Put the problem into stat. terms
\item Understand how the data were collected.
\bi
\item experimental data
\item observational data (survey data)
\item nonresponses/ missing values?
\item how are the data coded (for qualitative variables)?
\item possible errors?
\ei
\ei

\para{Initial data analysis}
\bi
\item numerical summaries (mean, sd, min, max, cor, etc.)
\item graphical summaries (boxplots, histograms, scatter plots, etc.)
\item look for outliers, data-entry errors, skewed or unusual distributions
\item cleaning and preparing data for analysis is an important part in practice
\ei

\para{Example} data(pima)
\bi
\item  Any missing or unusual observations?
\item  What does 0 represent for each variable?
\item  missing values are coded in various ways in practice!
\ei 


\para{Regression Analysis} models the relationship between
\bi
\item response $y$ (output or dependent variable)
\item and predictors $x_1, \ldots, x_k$ (input, independent or explanatory variables)
\ei

\para{Types of regressions}
\bi
\item simple regression: $k=1$
\item multiple regression: $k>1$
\item response $y$ is often continuous
\item predictors can be continuous or discrete (categorical)
\item ANOVA (analysis of variance): all predictors are qualitative
\item ANCOVA (analysis of covariance): a mix of continuous and discrete predictors
\item When response $y$ is discrete, {\bf logistic regression} or {\bf Poisson regression} could be used. 
\item multivariate multiple regression: multiple responses $y_1, \ldots, y_m$
\ei

\para{Objectives of regression analysis}
\bi
\item prediction of future observations
\item assessment of effects of, or relationship betweens $x$'s and $y$
\item description of data structure
\ei

\para{History}
\bi
\item Gauss developed {\bf least squares} methods in early 1800.
\item Galton coined the term "regression to mediocrity" in 1825 when he studied the heights of sons and fathers as
\[ \frac{y-\bar y}{SD_y} = r \frac{x - \bar x}{SD_x} \]
where $r$ is the correlation between $x$ and $y$.
\item phenomenon ({\bf regression effect}): sons of tall fathers tend to be tall but not as tall as their fathers while  sons of short fathers tend to be short but not as short as their fathers.
\ei

\para{Example} data(stat500)

\newpage
\section*{Chapter 2. Linear Model/Estimation}

\section*{1.   Linear  Model}  
Consider a general model
\[ Y=f(x_1, \ldots, x_k) + \epsilon \]
where $f()$ is an unknown function and $\epsilon$ is random error.  This model is too general and not useful at all.  A more practical approach is to restrict $f()$ to some parametrical form, say $f$ is a linear function of the predictors as follows.

\para{Linear Model:}
\[ y = \beta_0 + \beta_1 x_1 + \dots + \beta_k x_k +\epsilon,  \]
where $y$ is the response (dependent variable), $x_i$ are predictors (covariates, independent variables), $\beta_i$ are unknown regression coefficients, and $\epsilon$ is a random error.  Assume that $E(\epsilon)=0$ and $var(\epsilon)=\sigma^2$.
%\bi
%\item Assume that $\epsilon$ has mean 0 and variance $\sigma^2$, i.e., 
%\ei

\para{Remarks:} 
\be
\item In a linear model, the parameters $\beta_i$ enter linearly, the predictors do not have to be linear.
% Linear regression models include certain types of nonlinear models. For example, consider a second-order model 
\[ y = \beta_0 + \beta_1 x_1 + \beta_2 \log (x_2) + \beta_3 x_1 x_2  + \epsilon \]
%\[ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_1^2 + \beta_4 x_1 x_2 + \beta_5 x_2^2 + \epsilon .\]
%Let $x_3=x_1^2$, $x_4=x_1 x_2$ and $x_5=x_2^2$. Then the model becomes
%\[ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \beta_4 x_4 + \beta_5 x_5 + \epsilon,\]
is a linear  model.
\[ y = \beta_0 + \beta_1 x_1^{\beta_2} + \epsilon \]
is not a linear model.
\item linear models seem restrictive, but are actually very flexible, because predictors can be transformed and combined in many ways.
\item Truly nonlinear models are rarely absolutely necessary (unless supported by theory).
\ee


For data with $n$ observations $(y_i, x_{i1}, \ldots, x_{ik})$,  the model is:
\[ y_i = \beta_0 + \beta_1 x_{i1} + \dots + \beta_k x_{ik} +\epsilon_i , i=1, \dots, n. \]
\bi
\item Assume that $\epsilon_i$ are independent with mean 0 and common variance $\sigma^2$.
%\item $E(\epsilon)=0$ and $var(\epsilon)=\sigma^2$.
\ei

%It is usually assumed $\epsilon \sim N(0, \sigma^2 )$, then $y$ is also normally distributed, i.e., $y \sim N(\mu, \sigma^2 )$, where 
%\begin{eqnarray*}
%\mu = E(y) & = & \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p .
%\end{eqnarray*}

\para{In matrix form}, the model is
\[ {y = X} \beta+ \epsilon,  E( \epsilon)=0, cov( \epsilon)=\sigma^2 I  \]
where 
\[ y = \pmatrix{ y_1 \cr \vdots \cr y_n },  \ X = \left(\begin{array} {llll}
1 & x_{11} & \dots & x_{1k} \\
\vdots & \vdots & \vdots & \vdots \\
1 & x_{n1} & \dots & x_{nk} \\
\end{array}
    \right) ,
 \  \beta = \pmatrix{ \beta_0 \cr  \beta_1 \cr \vdots \cr \beta_k}, \ \epsilon=  \pmatrix{  \epsilon_1 \cr \vdots \cr \epsilon_n } \]
\bi
\item  $X$, an $n \times p$ matrix,  is called the {\bf model matrix}, where $p=k+1$ here.
\item  $\beta$  is a vector of unknown regression coefficients.
\item $\epsilon$ is a vector of random errors.
\ei


\section*{2. Estimation}
The \textbf{least squares estimate} of $\beta$ minimizes
\[  L = \sum_{i=1}^n \epsilon_i^2 = \epsilon^T \epsilon = {({y - X} \beta)}^T {({y - X} \beta)} \]
Let
\[ \frac{\partial L}{\partial \beta} = - 2 X^T ({y - X} \beta) = 0 \]
The {\bf normal equations} are  
\[ X^T X \beta = X^T y \]
The \textbf{least squares estimate} (LSE) is 
 \[ \hat{\beta}  = {(X^{T}X)}^{-1}X^Ty, \]
The {\bf predicted or fitted values} are 
\[ \hat{y} = X \hat{\beta} = X {(X^{T}X)}^{-1}X^Ty = H y,\]
where $H=  X {(X^{T}X)}^{-1}X^T$ is the {\bf hat matrix}. 
\bi
\item The hat matrix plays an important role in residual analysis.
\item Property: $H$ is idempotent, i.e., $H^2=H$.
\ei

The {\bf residuals} are 
\[  {\hat{\epsilon}} = y - X \hat{\beta} = y - H y = (I - H) y \]
The residual sum of squares (RSS) is 
\[ RSS = {\hat{\epsilon}}^T{\hat{\epsilon}} =  y^T (I - H) (I - H) y = y^T (I - H)y  \]

\para{An unbiased estimate of $\sigma^2$} is 
\[ \hat \sigma^2 = \frac{RSS}{n-p} = MSE, \]
where $n-p$ is the {\bf degrees of freedom} of the model.

The least squares estimate has a nice {\bf geometrical interpretation} (Fig. 2.1).

\vspace{2in}

\para{Properties of the least squares estimate $\hat{\beta}$}  % has a multivariate normal distribution with
\be
\item LSE $\hat{\beta}$ is an unbiased estimator of $\beta$, i.e., $E(\hat{\beta}) = \beta$.
\item The covariance depends on $\sigma$ and the model matrix,  $cov(\hat{\beta}) =  \sigma^2 {(X^{T} X)}^{-1}$.
\item LSE $\hat \beta$ is max. likelihood estimator (MLE) if $\epsilon$ has a normal distribution.
\item {\bf Gauss-Markov Theorem}:  If the model is correct and $E(\epsilon)=0$ and $cov(\epsilon)=\sigma^2 I$, then $\hat \beta$ is the best linear unbiased estimator (BLUE).
\ee

%\newpage
\para{Examples of calculating $\hat \beta$}
\be
\item mean only model:  $y_i = \mu + \epsilon_i$
	\[ X = 1_n,  \ \beta=\mu, \] \[ \hat \beta = {(X^{T}X)}^{-1}X^Ty = (1^T 1)^{-1} 1^T y = \frac{1}{n} 1^T y = \bar y \]
\item simple linear regression (one predictor): $ y_i = \beta_0 + \beta_1 x_i + \epsilon_i $
	\[ y=X \beta + \epsilon \mbox{ means }  \pmatrix{y_1 \cr \vdots \cr y_n} = \pmatrix{ 1 & x_1 \cr \vdots & \vdots \cr 1 & x_n} \pmatrix{ \beta_0 \cr \beta_1} + \pmatrix{\epsilon_1 \cr \vdots \cr \epsilon_n} \]
A simpler approach:
\[   y_i = \beta_0 + \beta_1 \bar x + \beta_1(x_i - \bar x) + \epsilon_i = \beta'_0 + \beta_1 (x_i -\bar x) + \epsilon_i \]
Then
\[  X = \pmatrix{  1 & x_1-\bar x \cr \vdots & \vdots \cr 1 & x_n -\bar x} ,  X^T X = \pmatrix{ n & 0 \cr 0 & \sum (x_i -\bar x)^2}  \]
\[ \pmatrix{ \hat\beta'_0 \cr \hat\beta_1} = (X^T X)^{-1}X^T y\] \[ =  \pmatrix{ \frac{1}{n} & 0 \cr 0 & \frac{1}{\sum (x_i -\bar x)^2}} \pmatrix{  1 & \cdots & 1 \cr x_1-\bar x & \cdots & x_n -\bar x} \pmatrix{y_1 \cr \vdots \cr y_n} \]
\[ \hat \beta'_0 = \frac{1}{n} \sum y_i = \bar y, \ \ \hat \beta_1 = \frac{\sum (x_i - \bar x) y_i}{\sum (x_i - \bar x)^2}, \ \  \hat\beta_0 = \bar y - \hat\beta_1 \bar  x \]

\ee

\para{Remarks:}
\bi
\item In general, it is quite complicated to compute the inverse of $X^T X$. 
\item However, the computation is easy if $X^T X$ is diagonal, which happens for designed experiments.
\item When $X^T X$ is diagonal, estimates of $\hat \beta_1, \ldots, \hat \beta_p$ are uncorrelated. (why?)
\ei

\section*{3. Goodness of fit}
The coefficient of multiple determination 
%\[ R^2 = \frac{SS_{reg}}{SS_T} = 1 - \frac{SS_E}{SS_T} = \frac{\sum (\hat y_i - \bar y)^2}{\sum (y_i -\bar y)^2} \]
\[ R^2 = \frac{\sum (\hat y_i - \bar y)^2}{\sum (y_i -\bar y)^2} = 1-  \frac{\sum (\hat y_i - y_i)^2}{\sum (y_i -\bar y)^2} = 1 - \frac{RSS}{Total SS} \]
is a measure of goodness of fit.  
\bi
\item It tells the percentage of variation in $y$ explained by the regression model.
\item $0 \le R^2 \le 1$.
\item Large $R^2$ (close to 1) indicates a good fit.
%\item Adding a new variable to a model will increase $R^2$.
\item A large $R^2$ does NOT necessarily imply that the regression model is a good one.
%\item The adjusted $R^2$ is another measure and preferred for variable selection.
%$$R^2_{adj} = 1 - \frac{SS_E/(n-p)}{SS_T/(n-1)} = 1- \frac{n-1}{n-p} (1- R^2) $$
\item Caution: If the model does not include the intercept, $R^2$ defined here is not meaningful.
\ei


\para{Example} data(gala)


\section*{4. Identifiability}
LSE normal equations:  
\[ X^TX \hat\beta = X^T y \]
where $X$ is an $n \times p$ matrix, $X^TX$ is a $p \times p$ matrix.
\bi
\item If $X^TX$ is singular ($rank(X^TX)<p$), the inverse does not exist and $\hat\beta$ is unidentifiable.
\item It happens when some columns of $X$ are linearly dependent.
\item Common in experimental design, e.g., 2-sample comparison experiment or one-way layout.
\[ y_{ij} = \mu + \alpha_i + \epsilon_{ij},  \ \ i=1,2; j=1,\ldots,m \]
%\item Software handle non-identifiability differently.
\item Caution: R fits the largest identifiable model automatically.
\item In practice, be careful about near un-identifiability (near collinearity).
\ei

\para{Example} data(gala)

\newpage
\para{\bf Gauss-Markov Theorem}: LSE $\hat \beta$ is BLUE.

\para{\bf Proof.} Let $A=(X^TX)^{-1}X^T$.  Then $AX= (X^TX)^{-1}X^T X=I$, 
\[ AA^T=(X^TX)^{-1}X^T X(X^TX)^{-1}=(X^TX)^{-1}. \]
LSE $\hat\beta=Ay$ is a linear combination of the random response vector $y$.
% and  $E(\hat\beta)=E(Ay)=AE(y)=AX\beta=\beta$.

Let $\hat\beta^*=By$ be any other linear transformation of $y$, where $B$ is a $p\times n$ matrix of fixed coefficients. Denote $D=B-A$. Then
\[ E(\hat\beta^*)=E(By)= BE(y)=BX\beta=(D+A)X\beta=DX\beta + \beta .\]
So $\hat\beta^*$ is  unbiased iff $DX=0$.
\[ V(\hat\beta^*)=V(By)= BV(y)B^T=B(\sigma^2I)B^T=\sigma^2 BB^T \]
\[ =\sigma^2 (A+D)(A+D)^T = \sigma^2 (AA^T+DD^T) \]
Note $AD^T=(X^TX)^{-1}X^TD^T=0$ and $D^TA=0$ since $DT=0$ for any unbiased estimator. So
\[ V(\hat\beta^*)= \sigma^2 ((X^TX)^{-1}+DD^T) \]
\[ V(\hat\beta^*)= V(\hat\beta)+\sigma^2 DD^T  \ge V(\hat\beta) \]
because $DD^T\ge 0$ is semi-positive definite (i.e., for any vector $a$, $a^TDD^Ta={\tilde{a}}^T {\tilde{a}}  = \sum {\tilde{a}_i}^2 \ge 0$).

Gauss-Markov Theorem implies that for any linear combination of $\beta$, say $a^T\beta$, LSE $a^T\hat\beta$ is BLUE.
\[ V(a^T\hat\beta)=a^TV(\hat\beta)a=\sigma^2 a^T(X^TX)^{-1} a \]
\[ V(a^T\hat\beta^*)=a^TV(\hat\beta^*)a= a^TV(\hat\beta)a +  \sigma^2 a^TDD^T a  \ge \sigma^2 a^T(X^TX)^{-1} a. \]

%\end{document}
\newpage
\section*{Chapter 3. Linear Model/Inference} %\para{Linear Model}:
\[ y=X \beta + \epsilon,  E(\epsilon)=0, cov(\epsilon)=\sigma^2 I \]
For statistical inference, it is further assumed that the errors follow a {\bf normal distribution}, i.e., $\epsilon \sim N(0, \sigma^2 I)$.

\section*{3.1 Hypothesis Tests} to compare two {\bf nested} linear models,
\bi
\item  Null $H_0$: smaller model $\omega$, $dim(\omega)=q$
\item  Alternative $H_1$: larger model $\Omega$, $dim(\Omega)=p$
\ei

\para{Basic idea:} 
\bi
\item Geometric representation (Fig. 3.1)
\vspace{2in}
\item reject $H_0$ if $RSS_{\omega} - RSS_{\Omega}$ is large.
\item or reject $H_0$ if the ratio $(RSS_{\omega} - RSS_{\Omega})/RSS_{\Omega}$ is large.
\ei
\para{The F test statistic}
\[ F = \frac{ (RSS_{\omega} - RSS_{\Omega})/(df_{\omega}-df_{\Omega})}{RSS_{\Omega}/df_{\Omega}} = \frac{ (RSS_{\omega} - RSS_{\Omega})/(p-q)}{RSS_{\Omega}/(n-p)}  \]
\bi
\item  This is indeed the maximum likelihood ratio test statistic.
\item When $\epsilon \sim N(0, \sigma^2 I)$, the $F$ statistic has an $F$ distribution with df $p-q$ and $n-p$, i.e., $F \sim F_{p-q, n-p}$, under the null model $\omega$.
\item Reject $H_0$ at level $\alpha$ if $F > F_{\alpha, p-q, n-p}$.  (upper $\alpha$th percentile, $F^{(\alpha)}_{p-q, n-p}$ )
\item P-value = $P(F_{p-q, n-p} > F)$.
\item works when $\omega$ is a  subset or subspace of $\Omega$.
\ei

\section*{3.2 Testing Examples}

\para{a. Test of all the predictors}.   Are any of the predictors useful? \\
Model: $y=X\beta + \epsilon$ or $ y = \beta_0 + \beta_1 x_1 + \dots + \beta_k x_k +\epsilon $
\begin{eqnarray*}
H_0 &:& \beta_1 =\ldots=\beta_k=0
\\ H_1 &:& \beta_j \ne 0 \mbox{ for at least one } j \ge 1
\end{eqnarray*}
Here $p=k+1$,  $q=1$,  $\omega$ is $y=\beta_0 + \epsilon$ and  $\hat\beta_0=\bar y$ under $\omega$. So
\begin{eqnarray*}
 RSS_{\omega} &=& (y-\bar y)^T(y-\bar y) = \sum(y_i -\bar y)^2 = SS_T  \mbox{ (or } TSS)
\\ RSS_{\Omega} &=& (y-X \hat \beta)^T(y-X \hat \beta) = SS_E \mbox{ (or } RSS)
\end{eqnarray*}
\para{ANOVA table} for testing all predictors
\begin{center}
\renewcommand{\arraystretch}{1.4}
\begin{tabular}{l |cccc}
Source & DF & SS & MS & F \\ \hline
Regression & $p-1$ & $SS_{reg}$ & $MS_{reg}=\frac{SS_{reg}}{p-1}$ & $F=\frac{MS_{reg}}{MS_E}$ \\
Error (or residual) & $n-p$ & $SS_E$ & $MS_E=\frac{SS_E}{n-p}$ & \\ \hline
Total & $n-1$ & $SS_T$ \\
\end{tabular}
\end{center}
P-value = $P(F_{p-1, n-p} > F)$

\para{Remarks}
\bi
\item A failure of reject the null is not the end of the analysis!  
\item ``Fail to reject" the null is not the same as ``accept" the null.
\item Rejecting the null does not imply the alternative model is the best model.
\ei


\para{Example}  data(savings)
%\begin{verbatim}
%library(faraway)
%data(savings)
%?savings
%savings
%g <- lm(sr ~ pop15 + pop75 + dpi + ddpi, savings)
%summary(g)
%g <- lm(sr ~ pop15 + pop75 + dpi + ddpi, savings)
%(tss <- sum((savings$sr-mean(savings$sr))^2))
%(rss <- deviance(g))
%df.residual(g)
%(fstat <- ((tss-rss)/4)/(rss/df.residual(g)))
%1-pf(fstat,4,df.residual(g))
%\end{verbatim}


 
 \para{b. Testing just one predictor}.   Can one particular predictor be dropped? \\
Model: $ y = \beta_0 + \beta_1 x_1 + \dots + \beta_k x_k +\epsilon $
\begin{eqnarray*}
H_0 : \beta_j=0, \ \  H_1:  \beta_j \ne 0 
\end{eqnarray*}
Here $p=k+1$,  $q=p-1=k$ and $\omega$ has one parameter less than $\Omega$.  So $RSS_{\omega}$ is obtained by fitting a regression without predictor $x_j$.  

The $F$ test statistic is
\[ F = \frac{ (RSS_{\omega} - RSS_{\Omega})/(p-q)}{RSS_{\Omega}/(n-p)} = \frac{ (RSS_{\omega} - RSS_{\Omega})}{RSS_{\Omega}/(n-p)}  \]
\[ \sim F_{1, n-p} \mbox{ under } H _0 \]
Alternatively, use a $t$-test
\[ t = \frac{\hat{\beta}_j}{se(\hat{\beta}_j)} =  \frac{\hat{\beta}_j}{\hat \sigma \sqrt{(X^TX)^{-1}_{jj}}} \sim t_{n-p} \mbox{ under } H_0 \]
% where $(X^TX)^{-1}_{jj}$ is the $j$th diagonal element of $(X^TX)^{-1}$ and $\hat \sigma  = \sqrt{RSS_{\Omega} /(n-p)}$.
\bi
\item Reject $H_0$ at level $\alpha$ if $|t| > t_{\alpha/2, n-p}$
\item P-value = $P( | t_{n-p}| \ge |t|)$.
\ei
\para{Remarks}
\bi
\item Fact: $F=t^2$ when testing one predictor only.
\item The result (of testing $\beta_j=0$) depends on other predictors $x_i$ in the model. (Why?)
%\item Permutation tests may be used without normality assumption. % The result is similar to the F test.
\ei

\para{Example}  data(savings)

%\begin{verbatim}
%g2 <- lm(sr ~ pop75 + dpi + ddpi, savings)
%(rss2 <- deviance(g2))
%(fstat <- (deviance(g2)-deviance(g))/(deviance(g)/df.residual(g)))
%1-pf(fstat,1,df.residual(g))
%sqrt(fstat)
%(tstat <- summary(g)$coef[2,3])
%2*(1-pt(sqrt(fstat),45))
%anova(g2,g)
%\end{verbatim}

 
 \para{c. Testing a pair of predictors}.   \\
 Suppose both $x_i$ and $x_j$ have p-values $>0.05$. Can we eliminate both from the model?
 
\[  \]
 
% \begin{verbatim}
%g3 <- lm(sr ~ pop15+dpi , savings)
%anova(g3,g)
%\end{verbatim}

%\vspace*{1in}

 \para{d. Testing a subspace}.   For model: 
 $$ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_k x_k +\epsilon $$
Can we combine two variables ($x_1+x_2$) to simplify the model as
\[ y = \beta_0 + \beta_{12} (x_1+x_2) + \beta_3 x_3 + \dots + \beta_{k} x_k +\epsilon  \]
The null hypothesis is a linear subspace of the original model.

What is the null hypothesis here?
\[   \]

\para{Example}  data(savings)

%\begin{verbatim}
%g <- lm(sr ~ .,savings)
%gr <- lm(sr ~ I(pop15+pop75)+dpi+ddpi,savings)
%anova(gr,g)
%gr <- lm(sr ~ pop15+pop75+dpi+offset(0.5*ddpi),savings)
%anova(gr,g)
%(tstat <- (0.409695-0.5)/0.196197)
%2*pt(tstat,45)
%tstat^2
%\end{verbatim}

\section*{3.3 Permutation Tests}
\bi
\item The F or t test assumes the errors are normally distributed.
\item A permutation test does not require the normality assumption.
\item Procedure:  
 \bi \item consider all $n!$ permutations of the response variable
  \item compute the F statistic for each permutation
  \item compute proportion of F statistics exceed the observed F statistic for the original response
  \ei
 \item The proportion is estimated by the p-value calculated in the usual way (assuming normal errors).
 \item When $n$ is large, do a random sample of the $n!$ permutations.
 \item To test one predictor ($x_i$): permute that predictor ($x_i$) rather than the response $y$.
\ei

\para{Example}  data(savings)



\section*{3.4 Confidence Intervals (CIs) for $\beta$} 

\para{The 100(1-$\alpha)$\%  CI} for $\beta_j$ is 
\[ \hat{\beta}_j \pm t_{\alpha/2, n-p} \cdot se(\hat{\beta}_j) = \hat{\beta_j} \pm t_{\alpha/2, n-p} \hat{\sigma}  \sqrt{ (X^T X)^{-1}_{jj}}   \]
\bi
\item The $100(1-\alpha)\%$ CI for $\beta_j$ does not contain 0 if and only if the $t$ test rejects $H_0: \beta_j=0$ (vs $H_1: \beta_j \ne 0$) at level $\alpha$.
\ei

\para{The 100(1-$\alpha)$\%  confidence region} for $\beta$ is 
\[ (\hat{\beta} - \beta)^T (X^TX)(\hat{\beta} - \beta) \le p {\hat \sigma}^2 F_{\alpha, p, n-p} \]
\bi
\item The confidence region  is ellipsoidally shaped.
\item The $100(1-\alpha)\%$ confidence region for $\beta$ does not contain the origin if and only if the $F$ test rejects $H_0: \beta=0$  at level $\alpha$.
\item CIs and confidence region may give conflict conclusions. Why?
\item The conclusion from confidence region is preferred to individual CIs.
\ei

\para{Example}  data(savings)

\bi
\item Distinguish the correlation between predictors, $cor(x_i, x_j)$, and correlation between estimators, $cor(\hat\beta_i, \hat\beta_j)$ 
\ei


\section*{3.5 CIs for Predictions}

For a new set of predictors $x_0$, the predicted response is $\hat y_0 = x_0^T \hat \beta$.
\[ var(\hat y_0) = var(x_0^T \hat \beta) = x_0^T cov(\hat \beta) x_0 = \sigma^2 x_0^T (X^TX)^{-1}x_0 \]
\bi
\item The 100(1-$\alpha)$\%  {\bf CI}  on the {\bf mean response} is 
\[ \hat y_0 \pm t_{\alpha/2, n-p} \hat \sigma \sqrt{x_0^T (X^TX)^{-1}x_0} \]
\item The 100(1-$\alpha)$\%  {\bf prediction interval}  for a {\bf future observation} is 
\[ \hat y_0 \pm t_{\alpha/2, n-p} \hat \sigma \sqrt{1+x_0^T (X^TX)^{-1}x_0} \]
\item The prediction intervals become wider as $x_0$ moves away from the data center.
\item Caution on extrapolation (when the new $x_0$ is outside the range of the original data)
\ei
%\para{In R}, ``confidence" interval stands for mean response while ``predition" interval for a future observation.

\para{Example}  data(gala)


\iffalse
\para{Testing linear combinations of the $\beta_j$'s}
For testing $H_0: a^T \beta= 0$, 
where $a^T \beta=\sum_{j=0}^{p} a_j \beta_j$ and $a=(a_0, \ldots, a_p)$,
the test statistic is 
\[ t= \frac{a^T \hat{\beta}}
{  \sqrt{ \hat{\sigma}^2  a^T{(X^{T} X)}^{-1} a} } , \]
which has a $t$ distribution with $N-p-1$ degrees of freedom under $H_0$.
\fi

\section*{3.6 Designed Experiments}
In a designed experiment, we have some control over $X$ and can make it {\bf orthogonal} or near {orthogonal}.  
Two important design features:
\bi
\item orthogonality
\item randomization
\ei
which allows us to make stronger conclusions from the analysis.

Suppose $X=(X_1 X_2)$  such that $X_1^T X_2 =0$ and $\beta=\pmatrix{ \beta_1 \cr \beta_2}$.
\[ y=X\beta  + \epsilon = X_1 \beta_1 + X_2 \beta_2 + \epsilon \]
\[ X^T X = \pmatrix{X_1^TX_1 & X_1^TX_2 \cr X_2^T X_1 & X_2^T X_2} = \pmatrix{X_1^TX_1 & 0 \cr 0 & X_2^T X_2}  \]
\[ \hat\beta = (X^TX)^{-1}X^Ty =  \pmatrix{X_1^TX_1 & 0 \cr 0 & X_2^T X_2}^{-1}  \pmatrix{X_1^T y  \cr  X_2^T y} \]
\[    \hat\beta_1 =(X_1^TX_1)^{-1}X_1^T y, \  \  \hat\beta_2= (X_2^T X_2)^{-1}X_2^T y \]
\bi
\item Estimate of $\beta_1$ does not depend on whether $X_2$ is in the model or not.
\item Does the significance of $\beta_1$ depend on whether $X_2$ is in the model or not?  Why?
\item Randomization provides protection against the influence of unknown {\bf lurking} variables 
\ei
 

\para{Example} data(odor) vs data(savings)
 
\section*{3.7 Observational Data}
\bi
\item Interpreting models built on observational data is problematic.
\[ \hat y = \hat \beta_0 + \hat\beta_1 x_1 + \hat \beta_2 x_2 \]
\item What does $\hat\beta_1$ mean? 
\item ``A unit change in $x_1$ will produce a change of $\hat\beta_1$ in the response"?
\item ``$\hat\beta_1$ is the effect of $x_1$ when all the other specified predictors are held constant"?
\item prediction is more stable than parameter estimation.
\ei

 \para{Example} data(savings)

\section*{3.8 Practical Difficulties}
\[ y=X\beta + \epsilon \]
We have gone over the linear model theory of estimation and inference.  What are the difficulties?


Albert Einstein: ``So far as theories of mathematics are about reality; they are not certain; so far as they are certain, they are not about reality."  
\bi
\item Nonrandom samples
\item Choice and range of predictors
\item Model misspecification. George Box: \\ {\bf ``All models are wrong but some are useful."} 

\item Publication and experimenter bias
\item Practical and statistical significance
\ei

%\end{document}
\newpage
\section*{Chapter 4. Linear Model/Diagnostics} % Residual Analysis: Assessment of Model Assumptions}
%\para{Linear Model}:
\[ Y=X \beta + \epsilon, \   \epsilon \sim N(0, \sigma^2 I) \] % E(\epsilon)=0, cov(\epsilon)=\sigma^2 I \]
\para{Regression Diagnostics:} Use residuals to check assumptions %:  $\epsilon \sim N(0, \sigma^2 I)$.
%Use residuals to check % the structure of the model and the model assumptions. 
\bi
\item Errors:  independent,  equal variance, normally distributed?
\item Model structure:  $E(y)=X \beta $?
\item Unusual observations: outliers, influential cases?
\ei

\section*{4.1 Checking Error Assumptions}
\para{Residuals}.  %  Recall that  $ \hat y = H y,  H = X(X^TX)^{-1}X^T, cov(y) = cov(\epsilon) = \sigma^2 I$
$  \hat \epsilon = y -\hat y = (I-H) y ,   $
\[ \cov(\hat y, \hat \epsilon) = \cov(Hy, (I-H) y) = H \cov(y) (I-H)^T= \sigma^2 H (I-H) =0 \]
\bi
\item $\hat y$ and $\hat \epsilon$ are uncorrelated if errors are independent.
\[ \cov(\hat \epsilon) = \cov((I-H) y) = (I-H) \cov(y) (I-H)^T  = \sigma^2 (I-H) \]
\item  $\hat \epsilon_i$ may be correlated even if  $\epsilon_i$ are independent.
\item  $\hat \epsilon_i$ may not have equal variance even if $\epsilon_i$ do.
\item but we still use residuals as the impact is usually small.
%\[ var (\hat \epsilon_i) = \sigma^2 (1-h_{ii}) \]
% where $h_{ii}$, the $i$th diagonal element of $H$, is the {\bf leverage} of the $i$th case  
%	\bi
%	\item A large leverage $h_{ii}$ will make $var(\hat\epsilon_i)$  small.
%	\item Points far from the center of the $x$ space have large leverages.
%	\ei
\ei


\para{Residual plots}
\bi
\item residuals vs. fitted:  plot $\hat\epsilon_i$ (or $r_i$ defined in section 4.2) vs. $\hat{y}_i$
\item residuals vs. predictors:  plot $\hat\epsilon_i$ (or $r_i$) vs. $x_{ij}$ for each $j$
\item check model structure, constant variance, nonlinearity, outliers (see Fig. 4.1)
\ei

\para{Example} data(savings)

\para{How to deal with nonconstant variance?}
\bi
\item Use weighted least squares (Section 6.1) if we know the form of nonconstant variance.
\item Transform the response
\ei

\para{How to choose a function $h()$ so that $\var(h(y))$ is constant?} 
\[ h(y) = h(Ey) +(y-Ey) h'(Ey) + \cdots \]
\[ \var(h(y)) = h'(Ey)^2 \var(y) + \cdots \]
Ignoring higher order terms,  for $\var(h(y))$ to be constant, make
\[ h'(Ey) \propto (\var(y))^{-1/2} \]
\[ h(y) = \int (\var(y))^{-1/2} dy = \int \frac{1}{sd(y)} dy \]

Two commonly used transformations:
\bi
\item If $\var(y)=\var(\epsilon) \propto (Ey)^2$, then $h(y)=\log(y)$
\item If $\var(y)=\var(\epsilon) \propto (Ey)$, then $h(y)=\sqrt{y}$
%\item 
\ei

\para{Example} data(gala)

\para{Normality?  Normal Q-Q plot:}  
\bi 
\item plot  sorted $\hat\epsilon_i$ (or $r_i$) against standard normal quantiles $\Phi^{-1}(\frac{i}{n+1})$ for $i=1, \ldots, n$, where $\Phi$ is the cdf for $N(0,1)$.
\item points shall close to a straight line if the model is correct and errors are $N(0, \sigma^2 I)$.
\ei

\para{What to do when non-normality is found?}
\bi
\item For short-tailed distribution, not a serious problem and can be ignored.
\item For skewed errors, try to transform the response.
\item For long-tailed errors, need to do something ... 
\ei


\para{Correlated errors?} 
\bi
\item For temporal or spatial data, errors are often correlated.
\item graphical check: 
	\bi \item plot $\hat\epsilon$ against time (or run order)
	\item  plot $\hat\epsilon_i$ against $\hat\epsilon_{i-1}$
	\ei
\item Durbin-Watson test % (dwtest() in library(lmtest) )
\[ DW= \frac{\sum_{i=2}^n (\hat\epsilon_i  - \hat\epsilon_{i-1})^2}{\sum_{i=2}^n \hat\epsilon_i^2} \]

\ei

\para{Example} data(airquality)


\section*{4.2 Finding unusual observations}
\para{Leverage}: $h_i=h_{ii}$, the $i$th diagonal element of $H=X(X^TX)^{-1}X^T$
%\item $\hat y$ and $\hat \epsilon$ are uncorrelated if errors are independent.
\[ \cov(\hat \epsilon) =  \sigma^2 (I-H) \mbox{ so } \var (\hat \epsilon_i) = \sigma^2 (1-h_{ii}) \]
% where $h_{ii}$, the $i$th diagonal element of $H$, is the {\bf leverage} of the $i$th case  
	\bi
	\item A large leverage $h_{ii}$ will make $\var(\hat\epsilon_i)$  small.
	\item Points far from the center of the $x$ space have large leverages.
	\item $\sum_{i} h_{ii} = p$.
	\item Rule of thumb: look more closely if $h_{ii} > 2p/n$.
	\ei


\para{Standardized  Residuals}  (or internally studentized residuals):
\[ r_i =\frac{\hat \epsilon_i}{\hat\sigma \sqrt{ 1-h_{ii}}} \]

\bi
\item $r_i$ has approximately variance 1 when  the model is correct and assumptions hold.
\item $r_i$ may be preferred to $\hat\epsilon_i$ (e.g., in R)
%\item For many designed experiments, $h_{ii}$ are constant.
\ei

\para{Outliers?}  Which one is an outlier in Fig. 4.10, p. 66? \\
How to test whether case $i$ is a potential outlier?  
\bi
\item Remove it from the data and recompute the estimates. 
\item Let $\hat\beta_{(i)}$ and $\hat\sigma_{(i)}$ be the estimate of $\beta$ and $\sigma$ without case $i$. 
\item The fitted value for case $i$ is
\[ \hat y_{(i)} = x_i^T \hat\beta_{(i)} \]
\item If $y_i - \hat y_{(i)} $ is large, case $i$ is an outlier.
\[ \var(y_i - \hat y_{(i)}) = \var(y_i) + \var(\hat y_{(i)}) = \sigma^2 +  \sigma^2 x_i^T (X_{(i)}^T X_{(i)})^{-1} x_i \]
\ei

\para{Studentized  residuals}  (or jackknife residuals):
\[ t_i =\frac{y_i - \hat y_{(i)}}{\hat\sigma_{(i)} \sqrt{ 1 +   x_i^T (X_{(i)}^T X_{(i)})^{-1} x_i}} \]
\bi
\item $t_i \sim t_{n-p-1}$ if case $i$ is not an outlier. Why?
\item $\hat\sigma_{(i)}^2$ is an unbiased estimator of $\sigma^2$ and has df=$n-p-1$.
%which has a $t$ distribution with df=$n-p-1$ if case $i$ is not an outlier.
%\bi
%\item $t_i$ will differ little from $r_i$ if case $i$ is not an outlier.
%\item $t_i$ will differ significantly from $r_i$ when case $i$ is an outlier.
\item $t_i$ is a function of $r_i$ (so no need to fit $n$ regressions)
\[ t_i =\frac{\hat \epsilon_i}{\hat\sigma_{(i)} \sqrt{ 1-h_{ii}}} =  r_i \sqrt{\frac{n-p-1}{n-p-r_i^2}} \]
\item declare case $i$ as an outlier (with Bonferroni correction) if $|t_i| > t_{\alpha/(2n), n-p-1}$. 
\item may not work if there are two or more outliers next to each other.
\item An outlier in one model may not be an outlier in another model.
\ei

\para{What should be done about outliers?} Read page 68.

 \para{Example} data(savings) and data(star)


\para{Influential Cases and Cook's Distance}  \\
Let $\hat\beta_{(i)}$ be the estimate of $\beta$ without case $i$.  \\
Let $\hat y_{(i)} = X \hat\beta_{(i)}$ be the fitted values without case $i$.

\para{\bf Cook's distance} is
\[ D_i =  \frac{(\hat y_{(i)} - \hat y)^T (\hat y_{(i)} - \hat y)}{p \cdot \hat\sigma^2} = \frac{(\hat\beta_{(i)} - \hat\beta)^T (X^TX)(\hat\beta_{(i)} - \hat\beta)}{p \cdot \hat\sigma^2} \]
\bi
\item Large $D_i$ implies case $i$ is influential (the estimates and fitted values will be quite different with or without case $i$).
%\item $D_i$ is a function of $r_i$ and leverage $h_{ii}$: 
\[ D_i = \frac{ r_i^2}{p} \frac{h_{ii}}{1-h_{ii}} \]
\item A highly influential case must have a large leverage or a large standardized residual.
\item An influential case may, but not necessarily, be an outlier; see Fig. 4.10, p. 66.
\ei

\para{Example} data(savings)

\newpage
\section*{4.3 Checking the model structure}
\bi
\item residual plots ($\hat\epsilon$ vs. $\hat y$ and $x_i$) may reveal problems of model structure, but
\item other predictors can impact the relationship between $y$ and $x_i$
\ei

\para{Added variable plot (or partial regression plot)} can help isolate the effect of $x_i$ on $y$
\bi
\item regress $y$ on all $x$ except $x_i$, and get residuals $\hat \delta$.
\item regress $x_i$ on all $x$ except $x_i$, and get residuals $\hat \gamma$.
\item plot $\hat \delta$ against $\hat \gamma$.
\item fit a regression line to the plot. The slope is $\hat\beta_i$.
\ei

\para{Partial residual plot} 
\[ y - \sum_{j \ne i} x_j \hat\beta_j = \hat y + \hat \epsilon -\sum_{j \ne i} x_j \hat\beta_j = \hat \epsilon + x_i \hat\beta_i  \]
\bi
\item plot $\hat \epsilon + \hat\beta_i x_i$ against $x_i$.
\item The slope of the fitted line is also $\hat\beta_i$.
\item better for nonlinearity detection than added variable plot.
\ei

\para{Example} data(savings)


%\para{Note:} There are other names for the studentized and jackknife residuals. For example, they are called standardized and studentized residuals in R respectively; see commands \texttt{rstandard} and \texttt{rstudent}. 


\newpage
\section*{Chapter 5. Problems with the Predictors} 
\para{What will change if we scale the predictors or response?}
\bi
\item If  $x_i \longrightarrow (x_i+a)/b$, then $\hat\beta_i \longrightarrow b \hat\beta_i$
\item  If  $y \longrightarrow b y$, then $\hat\beta \longrightarrow b \hat\beta$ and $\hat\sigma \longrightarrow b \hat\sigma$
\item The t-tests, F-tests and $R^2$ are unchanged.
\item Scaling or standardizing the predictors and response makes comparisons simpler.
\ei

\para{Example} data(savings)

\para{Collinearity}
\bi
\item If $X^TX$ is close to singular, we have collinearity (or multicollinearity).
\item Collinearity causes serious problems with the estimation of $\beta$ and interpretation.
\ei

\para{Ways to detect collinearity}
\bi
\item Examine correlation matrix.
\item Regress $x_j$ on all other predictors, get $R_j^2$.  A large $R_j^2$ (close to one) indicates a problem.
\item Examine the eigenvalues of $X^TX$.  Let $\lambda_1 \ge \ldots \ge \lambda_p$ be the eigenvalues.
\item The {\bf condition number}
\[ \kappa = \sqrt{\frac{\lambda_1}{\lambda_p}} \]
where $\kappa \ge 30$ is considered large.
\ei

Let 
\[  S_{x_jx_j} = {\sum_i (x_{ij}-\bar x_j)^2}  \]
Then 
\[ \var({\hat\beta_j}) = \sigma^2 \left(\frac{1}{1-R_j^2}\right) \frac{1}{S_{x_jx_j}} \]
\bi
\item $\var({\hat\beta_j})$ will be large if $S_{x_jx_j}$ is small.
\item $\var({\hat\beta_j})$ will be large if $R_j^2$ is close to 1.
\item $\var({\hat\beta_j})$ will be minimized if $R_j^2=0$.
\item {\bf variance inflation factor}:  vif =  $\frac{1}{1-R_j^2}$
\item design strategy:  
 \bi \item spread $X$ to make $S_{x_jx_j}$ large.
 \item make predictors orthogonal to make $R_j^2=0$.
 \ei
\item one cure of collinearity is amputation, i.e., drop some (redundant) predictors.
\ei


\para{Example} data(seatpos)

%\end{document}
\newpage
\section*{Chapter 6. Problems with the Errors} 
\[ Y=X \beta + \epsilon, \   \epsilon \sim N(0, \sigma^2 I) \] % E(\epsilon)=0, cov(\epsilon)=\sigma^2 I \]
\bi
\item Assume errors are  independent,  equal variance, normally distributed
\item What to do if errors are not?
\ei

\section*{6.1 Generalized Least Squares} 
If we know the error variance structure, say,
\[\var(\epsilon)=\sigma^2 \Sigma = \sigma^2 S S^T, \]
where $\Sigma = S S^T$ is the Choleski decomposition and $S$ is a triangular matrix. 
\[ Y=X \beta + \epsilon \]
\[ S^{-1} Y= S^{-1}  X \beta +  S^{-1} \epsilon  \]
Consider the transformed variables and model, 
\[ Y'=S^{-1} Y, X'=S^{-1}  X, \epsilon' = S^{-1} \epsilon \]
\[ Y'=X' \beta + \epsilon' \]
\[ \var( \epsilon') = \var(S^{-1} \epsilon) = S^{-1} \var( \epsilon) S^{-T} = S^{-1} \sigma^2 S S^T S^{-T} =  \sigma^2 I \]
So the Generalized LS (GLS)  is an ordinary LS (OLS) for the transformed  model, 
\[ \hat\beta = ({X'}^TX')^{-1}{X'}^TY' = (X^T \Sigma^{-1} X)^{-1} X^T \Sigma^{-1}Y \]
\[ \var(\hat \beta) = \sigma^2  ({X'}^T { X'})^{-1} =  \sigma^2  ({X}^T \Sigma^{-1} X)^{-1} \]
\[ RSS = ({Y'} - {X'}\hat\beta)^T  ({Y'} - {X'}\hat\beta) = (Y-X\hat\beta)^T \Sigma^{-1} (Y-X\hat\beta) \]
Diagnostics should be applied to the residuals $\hat\epsilon' = S^{-1} \hat\epsilon$.
\bi
\item  Main problem in practice: $\Sigma$ is unknown 
\item need to estimate $\Sigma$ from data, not easy
\item need to guess the form of $\Sigma$ in many cases 
\item Two popular correlation structures:
\bi 
\item completely symmetrical (CS):  $\Sigma_{ij} = \rho$ for $i\ne j$
\item autoregressive AR(1) :  $\Sigma_{ij} = \rho^{|i-j|}$
\[ \epsilon_{i+1} = \rho  \epsilon_{i} + \delta_i \]
\ei
\ei 

\para{Example} data(longley)

\section*{6.2 Weighted Least Squares}

WLS is a special case of GLS,
\bi
\item errors are uncorrelated but have unequal variances
\[ \Sigma = diag(1/w_1, \ldots, 1/w_n) \]
\item so $S=diag(\sqrt{1/w_1}, \ldots, \sqrt{1/w_n})$ and 
\item regress $\sqrt{w_i} y_i$ on $\sqrt{w_i} x_i$ ( and $\sqrt{w_i} 1$)
\item E.g., if $\var(\epsilon_i) \propto x_i$ suggest $w_i=1/x_i$.
\item E.g., if $y_i$ is the average of $n_i$ observations, $\var(\epsilon_i)= \sigma^2/n_i$ suggest $w_i=n_i$.
\item use  $\sqrt{w_i} \hat\epsilon_i$ for diagnostics
\ei

\para{Example} data(fps)

\newpage
\section*{6.3 Testing for Lack of Fit}
\bi
\item If the model is correct, $\hat\sigma^2$ is an unbiased estimator of $\sigma^2$.
\item but we do not know  $\sigma^2$ and 
\item $\hat\sigma^2$ depends on the model.
\item If we have {\bf replicates} in the data, we can have an estimate of  $\sigma^2$ that does not depend on any model.
\ei

Let $y_{ij}$ be the $i$th observation in the group of replicates $j$.
\bi
\item ``pure error" estimate of $\sigma^2$ is $SS_{pe}/df_{pe}$, where
\[  SS_{pe} = \sum_{j} \sum_{i} (y_{ij}-\bar{y}_j)^2 \]
\[ df_{pe} = \sum_j (\#replicates -1) = n -\#groups \]
\item this is indeed the $\hat\sigma^2$ from the one-way ANOVA model with group as a factor.
\[ y_{ij} = \mu + \tau_j + \epsilon_{ij}  \]
\item Caution: whether the replicates are genuine?
\item Caution: avoid overfitting with complex model.
\ei

\para{Example} data(corrosion)

\newpage
\section*{6.4 Robust Regression}
\bi
\item LSE is clearly the best when the errors are normal, but 
\item other methods are preferred for long-tailed error distributions.
\ei

\para{M-Estimation} chooses $\beta$ to minimize
\[ \sum_{i=1}^n \rho(y_i - x_i^T \beta) \]
\bi
\item $\rho(x)=x^2$ is just least squares (LS)
\item $\rho(x)=|x|$ is called least absolute deviation (LAD) regression or $L_1$ regression
\item Huber's method is a compromise between LS and  LAD
\[ \rho(x) = \left\{ \begin{array}{cc}
x^2/2  &  \mbox{ if } |x|\le c   \\
c\, x - c^2/2  &   \mbox{ otherwise } 
\end{array}
\right.
\]
where $c$ should be a robust estimate of $\sigma$, e.g., a value proportional to the median of $\hat\epsilon$.
\ei

Robust regression is similar to WLS with weight function $w(u)=\rho'(u)/u$; see page 99.
\bi
\item LS: $w(u)$ is constant
\item LAD: $w(u)=1/|u|$
\item Huber: 
\[ w(u) = \left\{ \begin{array}{cc}
1  &  \mbox{ if } |u|\le c   \\
c/|u|  &   \mbox{ otherwise } 
\end{array}
\right.
\]
\item weights $w_i=w(u_i)$ depend on the residuals $u_i = y_i - \sum_{j=1}^p x_{ij} \hat\beta_j$
\item so use an iteratively reweighted least squares (IRWLS) approach to fit.
\ei
\[ \widehat{\var}(\hat\beta) = \hat\sigma^2 (X^TWX)^{-1} \]
where $W=diag(w_1, \ldots, w_n)$. 

\para{Example} data(gala)


\para{Least Trimmed Squares} (LTS)  minimizes
\[ \sum_{i=1}^q \hat\epsilon^2_{(1)}  \]
where $q$ is some number less than $n$ and
\bi
\item  $\hat\epsilon_{(1)} \le \ldots \le \hat\epsilon_{(n)}$ are ordered residuals.
\item  LTS is a {\bf resistant} regression method and can tolerate a large number of outliers.
\item default choice of $q$ in \verb|ltsreg| is $\lfloor n/2 \rfloor + \lfloor(p+1)/2 \rfloor$ where $\lfloor x \rfloor$ is the largest integer $\le x$.
\ei

\para{Bootstrap} is a general method to obtain standard errors or confidence intervals
\bi
\item Generate $\epsilon^*$ by {\bf sampling with replacement} from  $\hat\epsilon_{1}, \ldots, \ \hat\epsilon_{n}$
\item Form $y^*= X\hat\beta + \epsilon^*$
\item Compute $\hat\beta^*$ from $(X, y^*)$
\ei

Read summary on page 105-106.

\para{Example} data(gala) and data(star)

\newpage
\section*{Chapter 7  Transformation}
\section*{7.1  Transforming the Response}
\para{Box-Cox Method}  For $y>0$, consider  $y \longrightarrow g_{\lambda}(y)$
\[  g_{\lambda}(y) = \left\{ \begin{array}{cc}
\frac{y^{\lambda}-1}{\lambda}  &   \lambda \ne 0   \\
\log y &  \lambda = 0
\end{array}
\right. \]
\bi
\item choose $\lambda$ to  maximize the log-likelihood
\[ L(\lambda) = - \frac{n}{2}\log(RSS_{\lambda}/n) + (\lambda -1) \sum \log y_i \]
\item  A $100(1-\alpha)\%$ CI for $\lambda$ is 
\[ \lambda : L(\lambda) > L(\hat\lambda) - \frac{1}{2} {\chi_1^2}^{(1-\alpha)} \]
\item Box-Cox method gets upset by outliers 
\item If $\max_i y_i / \min_i y_i$ is small, the Box-Cox method will not have much real effect.
\item choose a convenient $\lambda$ for easy of interpretation
\item remember to transform back to the original scale.  For example,
\[ \log(y) =\beta_0 + \beta_1x+ \epsilon \]
\[ y = \exp(\beta_0 + \beta_1x) \exp(\epsilon) \]
the errors enter {\bf multiplicatively} not {\bf additively} as they usually do.
%\item alternative:  generalized least squares w/o transformation
\ei

\para{Example} data(savings)
%\end{document}

\newpage
\section*{7.2  Transforming the Predictors}
\bi
\item  replace $x$ with more than one term, say $f(x)+g(x)+\cdots$
\item add cross-product terms, say $x_i x_j$, etc.
\item fit different models in different regions, e.g., subset regression
\ei

\para{Broken Stick Regression}  is continuous but non-smooth % for one predictor $x$:
\[ y =\beta_0 + \beta_1 B_l(x) + \beta_2 B_r(x)  + \epsilon \]
where $B_l(x)$ and $B_r(x)$ are two hockey-stick functions:
\[
 B_l(x) =  \left\{ \begin{array}{cc}
c-x  &  \mbox{ if } x < c   \\
0 &   \mbox{ otherwise } 
\end{array}
\right.
\ 
B_r(x) =  \left\{ \begin{array}{cc}
x-c  &  \mbox{ if } x > c   \\
0 &   \mbox{ otherwise } 
\end{array}
\right.
 \]

\para{Polynomials} are flexible and smooth
\[ y= \beta_0 + \beta_1 x + \ldots + \beta_d x^d  + \epsilon \]
\bi
\item Two ways to choose $d$:  sequentially add or delete terms
\item Do no eliminate lower order terms from the model
\item Orthogonal polynomials are useful
\[ y= \beta_0 + \beta_1 \phi_1(x) + \ldots + \beta_d \phi_d(x)  + \epsilon \]
where $\phi_i(x)$ is a polynomials of order $i$, and orthogonal, i.e., $\sum_x \phi_i(x) \phi_j(x) = 0$ for $i \ne j$.
\item Response surface models for more than one predictors
\[ y= \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_{11} x_1^2  + \beta_{22} x_2^2 + \beta_{12} x_1 x_2+ \epsilon \]
%\item Caution: avoid overfitting
\ei

\para{Example} data(savings)

\para{Regression Splines} are smooth and have local influence property
\[ y= \beta_0 + \beta_1 \phi_1(x) + \ldots + \beta_d \phi_d(x)  + \epsilon \]
where $\phi_i(x)$ are B-spline basis functions on the interval $[a, b]$ 
\bi
\item requires knotpoints $t_1, \ldots, t_k$
\item A basis function is continuous and smooth (derivative $\phi'_i(x)$ is continuous).
\item A cubic B-spline basis function is nonzero on an interval defined by four successive knots and zero elsewhere.
\item A cubic B-spline basis function is a cubic polynomial for each subinterval between successive knots.
\item A basis function integrates to one over its support.
\item The broken stick regression is an example of linear splines.
\ei

\para{Example} simulated data
\[ y= \sin^3(2\pi x^3) + \epsilon, \ \ \epsilon \sim N(0, (0.1)^2) \]

\para{Remarks}
\bi
\item many possible ways to transform the predictors
\item many different models  with various complexity
\item Complex models may be good for prediction but difficult to interpret 
\item For small data sets or where the noise level is high, standard regression is most effective.
\ei

%\end{document}
\newpage
\section*{Chapter 8  Variable Selection}
Intended to select the ``best" subset of predictors.  Why?
\bi
\item To explain the data in the simplest way. ``The simplest is best".
\item To reduce the noise or variation of the estimation.
\item To avoid collinearity
\ei

\para{How to choose subset?}
\bi
\item {\bf Backward elimination}: Start with all predictors, get rid of the least significant one, repeat until some criterion is reached (i.e., all predictors are significant at a chosen level)

\item {\bf Forward selection}: Start with no predictor, add the most significant predictor, repeat until some criterion is reached (i.e., all predictors are not significant at a chosen level)

\item {\bf Stepwise selection}: combination of backward and forward selection.

\item {\bf Criterion-based procedures}:  search over various models (e.g., all subset regressions) and  choose a model based on a criterion (e.g., AIC, BIC, Adjusted $R^2$, $C_p$, etc.)
\ei

\para{Akaike Information Criterion (AIC)}
\[ AIC = -2\max \mbox{log-likelihood} + 2 p  \]
For regression, $-2\max \mbox{log-likelihood} = n \log(RSS/n) + constant$.

\para{Bayes Information Criterion (BIC)}
\[ BIC = -2\max \mbox{log-likelihood} +  p \log n  \]

\para{Adjusted $R^2$}
\[ R^2_a = 1 - \frac{RSS/(n-p)}{TSS/(n-1)} = 1 - \left( \frac{n-1}{n-p} \right) (1 -R^2) = 1 - \frac{\hat\sigma^2_{model}}{\hat\sigma^2_{null}} \]


\para{Mallow's $C_p$}
\[ C_p = \frac{RSS_p}{\hat\sigma^2} + 2p -n \]
\bi
\item $RSS_p$ is the RSS from the model with $p$ predictors
\item $\hat\sigma^2$ is from the model with all predictors
\item for the model with all predictors, $\hat\sigma^2 = RSS_p/(n-p)$ so $C_p=(n-p)+ 2p-n = p$
\item If a model fits poorly, then $RSS_p$ will be large and $C_p > p$.
\item choose models with small $p$ and $C_p$ around or less than $p$
\ei

\para{Example} data(state)

\section*{Chapter 9 Shrinkage Methods}
\bd
\item {9.1 Principal Component Regression}
\item {9.2 Partial Least Squares}
\item {9.3 Ridge Regression}  $ \hat\beta = (X^TX + \lambda I)^{-1} X^T y $
\ed

\section*{Chapter 10 Statistical Strategy and Model Uncertainty}

Diagnostics,  Transformation, Variable selection, Diagnostics, ...

\para{Chapter 11  A Complete Example}  % Insurance Redlining ---

\para{Chapter 12   Missing Data}: Imputation  % Insurance Redlining ---

\newpage
\section*{Chapter 13 Analysis of Covariance}
\bi
\item deal with a mixture of quantitative (numeric) and qualitative (categorical) predictors
\item use {\bf dummy variables} or {\bf contrast coding} for qualitative predictors
\item for a factor of $k$ groups (or levels), define $k$ dummy variables, $d_1, \ldots, d_k$
%\[ d_{ij} = \left\{ \begin{array}{cl} 1  &  \mbox{if  case $i$ has level $j$}    \\  0 &   \mbox{ otherwise }  \end{array} \right.  \mbox{ for $i=1,\ldots,n; j=1,\ldots, k$} \]
%\item Note: $d_{i1}+\cdots + d_{ik} = 1$
\[ d_{j} = \left\{ \begin{array}{cl} 1  &  \mbox{if  a case belongs to group $j$}    \\  0 &   \mbox{ otherwise }  \end{array} \right.  \]
\item Note: $d_{1}+\cdots + d_{k} = 1$
\item so only need $k-1$ dummy variables
\item or use any $k-1$ linear combinations of the $k$ dummy variables
\item The default choice in R is the {\bf treatment coding}, which use the last $k-1$ dummy variables, so the first group is the reference (or base) group.
\item The choice of coding does not affect the $R^2$, $\hat\sigma^2$ and overall $F$-statistic.
\item but it does affect the $\hat\beta$. 
\item The interpretation of the estimation depends on the coding.
\item interaction between a categorical and numeric variable is often of interest.
\ei

\newpage
\para{Example} data(sexab) \\
 $y=ptsd$ (Post-traumatic stress disorder on standard scale) \\
 $x=cpa$ (Childhood physical abuse on standard scale) \\
 $d=csa$ (Childhood sexual abuse - abused or not abused)
  \[  d =csa = \left\{ \begin{array}{cl} 1  &  \mbox{Abused}    \\  0 &   \mbox{ NotAbused }  \end{array} \right. \]
Consider 3 possible models
\begin{eqnarray*}
model  &1:&  y= \beta_0 + \beta_1 x + \epsilon \\
model  &2:& y= \beta_0 + \beta_1 x + \beta_2 d + \epsilon \\
model  &3:& y= \beta_0 + \beta_1 x + \beta_2 d + \beta_3 x \cdot d + \epsilon
\end{eqnarray*}
\bi
\item model 1: a single regression line  (csa has no effect on ptsd)
\item model 2:  two parallel regression lines with same slope \\ (the effect of csa does not depend on cpa)
\item model 3: two separate regression lines \\ (the effect of csa depends on cpa)
\ei
%In R,
%\begin{verbatim}
%ptsd ~ cpa
%ptsd ~ cpa + csa
%ptsd ~ cpa + csa + cpa:csa  (or ptsd ~ cpa * csa)
%\end{verbatim}

%\end{document}
\clearpage
 \begin{flushleft} \normalsize
 {\bf Stats 201B  (W13)} Regression Analysis: Model Building, Fitting, and Criticism %\hfill Hongquan Xu \\
\end{flushleft}

{\small
\bi
%\item  Textbook: Chapter 10 of  Montgomery (2005).   ``Design and Analysis of Experiments,"  Wiley.
%\item  Textbook: J. J. Faraway (2005). ``Linear Models with R," Chapman \& Hall.
\item  Textbook: J. J. Faraway (2006). ``Extending the Linear Model with R," Chapman \& Hall.
\ei
}

\section*{Chapter 2 Binomial Data} 
Linear model: 
\[ Y_i = \beta_0 + \beta_1 x_{i1} + \dots + \beta_k x_{ik} +\epsilon_i , \] %i=1, \dots, n. \]
where the error $\epsilon_i$ are independent and follow $N(0,\sigma^2)$, so
\bi
\item the $Y_i$ are independent and follow $N(\mu_i, \sigma^2)$ where
%\item so $y_i \sim N(\beta_0 + \beta_1 x_{i1} + \dots + \beta_k x_{ik},\sigma^2)$.
\[ \mu_i=E(Y_i) = \beta_0 + \beta_1 x_{i1} + \dots + \beta_k x_{ik} \]
%\item $E(\epsilon)=0$ and $var(\epsilon)=\sigma^2$.
\ei

\section*{2.2 Binomial Regression Model} 
Consider data 
\[ (y_i, n_i-y_i, x_{i1}, \ldots, x_{ik}), \mbox{ for $i=1, \ldots, n$} \]
where 
\bi
\item $y_i$ and $n_i-y_i$ are the number of ``successes'' and ``failures", respectively, out of $n_i$ independent trials
\item so the response has a {\bf binomial distribution}, $Y_i \sim B(n_i, p_i)$ 
\[  P(Y_i = y_i ) = {n_i \choose y_i} p_i^{y_i} (1-p_i)^{n_i-y_i} \]
\[ E(Y_i) = n_i p_i, \ \ \var(Y_i)=n_i p_i(1-p_i) \]
\item assume that the $Y_i$ are independent 
\item $x_{i1}, \ldots, x_{ik}$ are covariates
\item model the relationship between the $p_i$ and  predictors $x_{i1}, \ldots, x_{ik}$ as
\[ \eta_i = g(p_i) =  \beta_0 + \beta_1 x_{i1} + \dots + \beta_k x_{ik} \]
\item $\eta = g(p)$ is called the {\bf link function}
\item since $0\le p_i \le 1$, it is not appropriate to let $\eta_i=p_i$
\item {link function} $g$ is monotone and satisfies $0\le p=g^{-1}(\eta) \le 1$
\ei

%\newpage
\para{Common choices of link functions $\eta=g(p)$} 
\bi
\item Logit: $\eta = \log(p/(1-p))$  ({\bf logistic regression})
\item Probit: $\eta=\Phi^{-1}(p)$ where $\Phi(x)$ is the cdf of $N(0,1)$
\item Complementary log-log: $\eta=\log(-\log(1-p))$
\ei

\para{Inverse of link functions $p=g^{-1}(\eta)$} 
\bi
\item Inverse of Logit: $p = \frac{\exp(\eta)}{1+\exp(\eta)}$ 
\item Inverse of Probit: $p=\Phi(\eta)=P(Z \le \eta)$ where $Z \sim N(0,1)$
\item Inverse of Complementary log-log: $p = 1 - \exp(-\exp(\eta))$
\ei

\iffalse %
library(faraway)
par(mfrow=c(1,2))
p=seq(.01, .99, .01)
matplot(p, cbind(logit(p), qnorm(p),  log(-log(1-p))), type="l", ylab="g(p)", main="Link functions")
legend(0, 4, c("Logit", "Probit", "C-log-log"), lty=1:3, col=1:3, title="Link function")

eta = seq(-6,6, .1)
matplot(eta, cbind(ilogit(eta), pnorm(eta), 1-exp(-exp(eta))  ), type="l", ylab="p", main="Inverse of link functions")
legend(-5, 1, c("Logit", "Probit", "C-log-log"), lty=1:3, col=1:3, title="Link function")
\fi %

\para{Max Likelihood Estimator (MLE) $\hat\beta$}  max the log-likelihood 
%\[ l(\beta) = \sum_{i=1}^n \left[y_i \log \frac{p_i}{1-p_i} - n_i \log \frac{1}{1-p_i}+\log {n_i \choose y_i} \right] \]
\[ l(\beta) = \sum_{i=1}^n \left[y_i \log (p_i) + (n_i-y_i) \log ({1-p_i})+\log {n_i \choose y_i} \right] \]
For the logit link function, this is equivalent to
\[ l(\beta) = \sum_{i=1}^n \left[y_i \eta_i - n_i \log(1+\exp(\eta_i))+\log {n_i \choose y_i} \right] \]

\para{Example} data(orings)

\section*{2.3 Inference} 
To compare two models $\omega$ and $\Omega$ with $s$ and $l$ parameters, $s < l$, the likelihood ratio test (LRT) statistic is
\[  \Delta = 2\log \frac{L(\Omega)}{L(\omega)} = 2\log L(\Omega) - 2\log L(\omega) \]
\bi
\item  is approximately $\chi^2$ distributed with $df=l-s$ under the smaller model, if the $n_i$ are relatively large.
\item For a  $\chi^2$ r.v. with $df=d$, 
\[ E(\chi^2_d)= d, \  \var(\chi^2_d)=2d \]
\ei

When the larger model is saturated, $\hat p_i = y_i/n_i$, and the LRT statistic becomes
\[ D = 2 \sum_{i=1}^n \left[y_i \log \frac{y_i}{\hat y_i} + (n_i -y_i) \log \frac{n_i-y_i}{n_i -\hat y_i} \right] \]
where $\hat y_i$ are the fitted values from the smaller model.
\bi
\item The $D$ is called {\bf deviance} and measures how good the smaller model fits (compared to the saturated model).
\item The {deviance} is a measure of goodness of fit.
\item The deviance $D$ is approximately $\chi^2$ distributed with $df=n-s$, if the $n_i$ are relatively large (and the model is correct).
\item For a good approximation, it is often requested that $n_i \ge 5$.

\item {\bf Residual deviance} in R is the deviance for the current model.
\item {\bf Null deviance} is the deviance for a model with intercept only. 
\item approximate CI for $\beta_i$ is:   $\hat\beta_i \pm z_{\alpha/2} se(\hat\beta_i)$  or use \texttt{confint()}
\ei

\para{Example} data(orings)

\section*{2.5 Interpreting Odds} 
\bi
\item  Odds are used to represent chance in bets.
%\item A 3-1 {\bf against} bet would pay \$3 for every \$1 bet; 
\item Let $p$ be the probability of success and $o$ be the odds, 
\[  \frac{p}{1-p} = o \ \ \ \ p = \frac{o}{1+o} \]
%\item A 3-1 {\bf against} bet has $o=1/3$ and $p=1/4$;
\item A 3-1 {\bf on} bet has $o=3$ and $p=3/4$.
\item A 3-1 on bet would pay only \$1 for every \$3 bet.
\ei

\para{For logistic regression},
\[  \log(odds) = \log\left(\frac{p}{1-p} \right)=  \beta_0 + \beta_1 x_{1} + \dots + \beta_k x_{k} \]
\bi
\item $\beta_i$ has a simple interpretation in terms of log-odds
\item a unit increases in $x_1$ with all other predictors held fixed increases the log-odds of success by $\beta_1$ (or increases the odds of success by a factor of $\exp(\beta_1)$).
\ei

\para{Example} data(babyfood)


\section*{2.7 Choice of Link Function} 
\bi
\item when $p$ is moderate (not close to 0 or 1),  the link functions are similar.
\item Larger differences are apparent in the tails; see R plots.
\item but for very small $p$, one needs a very large amount of data to obtain just a few successes 
\ei


\para{Example} data(bliss)

%\section*{2.9 Goodness of Fit} 

\section*{2.10 Prediction and Effective Doses} 
\bi
\item To predict the outcome for given covariates, say $x_0$,
\[ \hat\eta = x_0^T \hat\beta \]
\[ \hat{\var}(\hat\eta) = x_0^T \hat{\cov}(\hat\beta) x_0 = x_0^T (X^T W X)^{-1} x_0 \]
where $\hat{\cov}(\hat\beta)=(X^T W X)^{-1}$ can be extracted using the \texttt{cov.unscaled} component of the model summary.

\item To get an answer in probability,  transform back using $p=g^{-1}(\eta)$
\ei

\para{Effective Doses}: to determine the $x$ value  for a  given $p$
\bi
\item For a simple logistic model
\[ \logit(p)= \log\left(\frac{p}{1-p} \right) = \beta_0 + \beta_1 x \]
the effective dose $x_p$  for prob. of success $p$ is
\[ x_p = \frac{\logit(p)-\beta_0}{\beta_1} \]
\item ED50 stands for {effective dose} for which there will be a 50\% chance of success;  $ED50= -\beta_0 / \beta_1$.
%\item LD50 stands for {\bf lethal dose} for which there will be a 50\% chance of killing a subject.
%\item when $p=1/2$, we get $ED50= -\beta_0 / \beta_1$.
\item To determine the standard error, use the delta method
\[ \var (h(\hat\theta)) \approx h'(\hat\theta)^T \var(\hat\theta)h'(\hat\theta) \]

\[ \theta =  \pmatrix{\beta_0 \cr \beta_1}, \  h(\theta) = \frac{\logit(p)-\beta_0}{\beta_1}, \ h'(\theta)  = \pmatrix{-1/\beta_1 \cr \beta_0/ \beta_1^2} \]
% = \pmatrix{\partial h(\theta)/\partial\beta_0 \cr \partial h(\theta)/\partial\beta_1}

%\item \textt{dose.p} 
\ei

\para{Example} data(bliss)

%\section*{2.11 Overdispersion} 

%\section*{Chapter 3 Count Regression} 
\section*{Chapter 3.1 Poisson Regression} 
If $Y$ has a Poisson distribution with mean $\mu > 0$, 
\[ P(Y=y) = \frac{\mu^y}{y!}e^{-\mu}, \ \ \ y=0, 1, 2, \ldots \]
\bi
\item $ E(Y)=\var(Y)=\mu$
\item Poisson distributions arise naturally when we count the number of (independent) events in a given time period.
%\ei
%\para{Poisson or Binomial?}
%\bi
\item Binomial $B(n,p)$ has a fixed total number of counts, which is $n$
\item If $n$ is large and $p$ is small, Poisson($\mu$) is a good approximation of $B(n,p)$ with $\mu=np$.
\ei

\para{Poisson regression model} (with link function $\eta=\log(\mu)$)
\[ \eta_i = \log(\mu_i) = x_i^T \beta \] 
The log-likelihood is
\[ l(\beta) = \sum_{i=1}^n (y_i \log(\mu_i) - \mu_i - \log(y_i!)) = \sum_{i=1}^n (y_i x_i^T \beta - \exp(x_i^T\beta) - \log(y_i!)) \]
Differentiating wrt $\beta_j$ gives the MLE $\hat\beta$ as the solution to
\[ \sum_{i=1}^n (y_i - \exp(x_i^T\hat\beta)) x_{ij} =0   \mbox{ for any } j \]
which can be written as
\[ X^T y = X^T \hat\mu,  \ \ \ \ \hat\mu=\exp(X\hat\beta) \]
\para{Deviance} (also known as the G-statistic)
\[ D = 2 \sum_{i=1}^n (y_i \log(y_i/\hat\mu_i) -  (y_i - \hat\mu_i)) \]
\para{Example} data(gala)
\clearpage

\para{Overdispersion}
\bi
\item  Poisson (or Binomial) distributions have only one parameter $\mu$ (or $p$) so it is not very flexible for empirical fitting purposes.
\item Overdispersion or  underdispersion can occur in Poisson (or Binomial) models; see text for possible reasons.
\item When overdispersion happens, the standard errors of the estimates are inflated.
\item One approach is to introduce a {\bf dispersion parameter} $\phi$, which can be estimated as 
\[  \hat\phi = \frac{X^2}{n-p} \]
where $X^2$ is the Pearson's $X^2$ statistic for goodness of fit.
\[ X^2 =\sum_{i=1}^n \frac{ (y_i - \hat \mu_i)^2}{\hat \mu_i} \]
\item For overdispersion: $\phi>1$;  for underdispersion: $\phi<1$.
%\item regular Poisson model, $\phi=1$.
\item When overdispersion is considered, an $F$-test rather than a $\chi^2$ test should be used.
\ei

\para{Example} data(gala)

\clearpage
\section*{Chapter 6 Generalized Linear Models} 

%Generalized Linear Model (GLM) includes LM, logisitic 

\section*{6.1 GLM Definition}
The distribution of response $Y$ is from the \textbf{exponential family}:
\[ f(y|\theta, \phi) = \exp\left[ \frac{y \theta - b(\theta)}{a(\phi)} + c(y, \phi) \right] \]
\bi
\item $\theta$ is called the {\bf canonical parameter} and represents the location
\item $\phi$ is called the {\bf dispersion parameter} and represents the scale.
\ei

\para{Examples of exponential family}
\bi
\item Normal or Gaussian:
\[ f(y|\theta, \phi) = \frac{1}{ \sqrt{2\pi} \sigma}  \exp\left[ -\frac{(y -\mu)^2}{2\sigma^2} \right]  \]
\[ =\exp\left[ \frac{y \mu - \mu^2/2}{\sigma^2} -\frac{1}{2} \left(\frac{y^2}{\sigma^2} + \log(2\pi\sigma^2)\right) \right] \]
so we can write $\theta=\mu, \phi=\sigma^2$, 
\[ a(\phi)=\phi, b(\theta)=\theta^2/2, c(y,\phi)=-(y^2/\phi+\log(2\pi\phi))/2 \]

\item Poisson:
\[ f(y|\theta, \phi) = \frac{\mu^y}{y!} e^{-\mu} = \exp[ y \log(\mu) -\mu - \log(y!) ] \]
%\[ =\exp\left[ \frac{y \mu - \mu^2/2}{\sigma^2} -\frac{1}{2} \left(\frac{y^2}{\sigma^2} + \log(2\pi\sigma^2)\right) \right] \]
so we can write $\theta=\log(\mu), \phi \equiv 1$, 
\[ a(\phi)=1, b(\theta)=\exp(\theta), c(y,\phi)=-\log(y!) \]

\item Binomial $B(m,p)$:  % replace \mu with p
\begin{eqnarray*}
 && f(y|\theta, \phi) = {m \choose y} p^y (1-p)^{m-y}  \\
&=& \exp \left[ y \log(p) +(m-y)\log(1-p)+\log {m \choose y} \right]  \\
 &=& \exp\left[ y \log \frac{p}{1- p} + m\log(1-p) +\log {m \choose y} \right] 
\end{eqnarray*}
so we can write $\theta=\log \frac{p}{1- p}, \phi \equiv1, a(\phi)=1$, 
\[ b(\theta)=-  m\log(1-p) = m\log(1+e^\theta)) , c(y,\phi)=\log{m \choose y} \]
\ei


\para{mean and variance of exponential family distributions}
\[ E(Y) = \mu = b'(\theta), \ \ \var(Y) = b''(\theta) a(\phi) \]
\bi
\item mean is a function of $\theta$ only;
\item variance is a product of functions of the location and the scale.
\item $V(\mu)=b''(\theta)$ is called the {\bf variance function} and describes how the variance relates to the mean.
\item For Gaussian case, $b''(\theta)=1$ and so the variance is independent of the mean.
\ei

\para{Linear predictor}
\[  \eta =  \beta_0 + \beta_1 x_{1} + \dots + \beta_k x_{k} = x^T \beta \]

\para{Link function}
\[ \eta = g(\mu) \] 
\bi
\item $g$ links mean response, $\mu=E(Y)$, to covariates through $\eta$
\item $g$ should be monotone, continuous and differentiable
\item {\bf canonical link} function satisfies 
\[ \eta=g(\mu)=\theta \longrightarrow g(b'(\theta)) = \theta \]
\begin{center}
\begin{tabular}{l l l}
Family & Canonical Link  & Variance function $b''(\theta)$ \\ \hline
Normal & $\eta=\mu$ & 1 \\
Poisson & $\eta=\log \mu$ & $\mu$ \\
%Binomial & $\eta=\log \frac{\mu}{1-\mu}$  & $\mu(1-\mu)$ % 
Binomial & $\eta=\log \frac{p}{1-p}$  & $mp(1-p)$
\end{tabular}
\end{center}
\item If a canonical link is used, $X^T Y$ is {\bf sufficient} for $\beta$.
\ei

\section*{6.2 Fitting a GLM}
One-step Taylor expansion 
\[ g(y) \approx  g(\mu) + (y-\mu) g'(\mu) \equiv z  \]
so 
\[ {\var}(z) = (g'(\mu))^2 \var(y) = (g'(\mu))^2 V(\mu) a(\phi) \] %\equiv 1/w \]
% \[ \widehat{\var}(z) = (g'(\mu))^2 V(\hat\mu) = 1/w \]

The following IRWLS procedure is used to fit a GLM:
\be
\item  Set initial estimates $\hat\mu_0$ and $\hat\eta_0=g(\hat\mu_0)$.
\item Form the ``adjusted dependent variable'' $z_0 = \hat\eta_0 + (y-\hat\mu_0) g'(\hat\mu_0)$.
\item Form the weights $w_0 = [ (g'(\hat\mu_0))^2 V(\hat\mu_0) ] ^{-1}$.
\item Re-estimate $\beta$ using WLS with response $z_0$ and weights $w_0$. %; $\hat\beta=(X^TWX)^{-1}X^TWz$
\item Update $\hat\eta_1 = x^T \hat\beta$ and $\hat\mu_1=g^{-1}(\hat\eta_1)$.
\item Iterate steps 2--5 until convergence.
\ee

Estimates of variance are obtained from:
\[ \widehat{\var}(\hat\beta) = (X^TWX)^{-1} \hat\phi \]


\para{Example} data(bliss)

For a binomial response, consider $p$ (rather than $\mu=mp$), 
\[ \eta=g(p) = \log \frac{p}{1-p},  g'(p) =  \frac{1}{p(1-p)}, V(p)= \frac{p(1-p)}{m}, w=mp(1-p) \]


\section*{6.3 Hypothesis Tests}
\bi
\item For a saturated (or full) model $\hat\mu=y$.
\item  Likelihood ratio test (LRT) statistic for comparing the current model with a saturated model is
\[ 2(l(y,\phi |y) - l(\hat\mu, \phi | y)) \]
\item when $a(\phi)=\phi/w_i$, this simplifies to
\[ D(y,\hat\mu)/\phi = \sum_{i} 2w_i (y_i(\tilde\theta_i - \hat\theta_i) - b(\tilde\theta_i)+b(\hat\theta_i))/\phi \]
where $\tilde\theta$ and $ \hat\theta$ are the estimates under the full and current model, respectively.
\item $D(y,\hat\mu)$ is called the {\bf deviance}.
\item $D(y,\hat\mu)/\phi$ is called the {\bf scaled deviance}.

\begin{center}
\begin{tabular}{l l }
GLM & Deviance $D(y,\hat\mu)$ \\ \hline
Normal & $\sum_i (y_i - \hat\mu_i)^2$ \\
Poisson & $2 \sum_i (y_i \log(y_i/\hat\mu_i) -  (y_i - \hat\mu_i))$ \\
%Binomial & $\eta=\log \frac{\mu}{1-\mu}$  & $\mu(1-\mu)$ % 
Binomial & $2 \sum_i \left[y_i \log \frac{y_i}{\hat \mu_i} + (n_i -y_i) \log \frac{n_i-y_i}{n_i -\hat \mu_i} \right]$
\end{tabular}
\end{center}
%\item For Gaussian distribution, $D(y,\hat\mu)=\sum_i (y_i - \hat\mu_i)^2$ is RSS.
\item An alternative to deviance is Pearson's $X^2$ statistic:
\[ X^2 = \sum_i \frac{(y_i - \hat\mu_i)^2}{V(\hat\mu_i)} \]
where $V(\hat\mu) = \var(\hat\mu)$.

\item The scaled deviance $D(y,\hat\mu)/\phi$ and Pearson's $X^2$ statistic are asymptotically $\chi^2$ distributed.
\item If we know $\phi$, we can test the goodness of fit with a $\chi^2$ test.
\item For nested models $\omega$ and $\Omega$,
$D_{\omega} - D_{\Omega}$ is asymptotically $\chi^2$.
\item If we know $\phi$, use $\chi^2$ test.
\item If we don't know $\phi$, use an $F$-statistic (which is approximately $F$ distributed)
\[ \frac{(D_{\omega} - D_{\Omega})/(df_{\omega} - df_{\Omega})}{\hat\phi} \]
where $\hat\phi = X^2/(n-p)$ is a good estimate of the dispersion.
\item For the Gaussian model, $\hat\phi = RSS_{\Omega}/df_{\Omega}$, the F test is exact.
\ei

\para{Example} data(bliss) 

\newpage
\section*{6.4 GLM Diagnostics}

\para{Residuals}
\bi
\item The {\bf response residual} $\hat\epsilon = y -\hat\mu$, which do not have constant variance.
\item {\bf Pearson residual}
\[  r_P = \frac{y -\hat\mu}{\sqrt{V(\hat\mu)}} \] 
where $V(\mu)\equiv b''(\theta)$. Note $\sum r_P^2 = X^2$.
\item Let deviance $D=\sum d_i = \sum r_D^2$. The {\bf deviance residual} is
\[ r_D =sign(y-\hat\mu) \sqrt{d_i} \] 
\item The {\bf working residual} is a by-product of the IRWLS fitting procedure.
\ei

\para{Leverage and influence}: 
\bi
\item  The {\bf hat matrix} is 
\[ H = W^{1/2} X(X^TWX)^{-1}X^TW^{1/2} \]
where $W=\mbox{diag}(w)$ and $w$ are the weights used in fitting. % the GLM.
\ei

\para{Standardized  residuals} :
\[ r_{SD} =\frac{r_D}{ \sqrt{ \hat\phi(1-h_{ii})}} \]
%\para{Studentized or Jackknife  residuals} :
%\[ sign(y-\hat\mu) \sqrt{ (1-h_{ii})r^2_{SD} + h_{ii} r^2_{SP}} \]
%where $r_{SP} = r_P/\sqrt{1-h_{ii}}$.
\para{\bf Cook's distance} is
\[ D_i =  \frac{(\hat\beta_{(i)} - \hat\beta)^T (X^TWX)(\hat\beta_{(i)} - \hat\beta)}{p \cdot \hat\phi} \]
\para{Example} data(bliss) 

\newpage
\para{Model diagnostics}
\bi
\item Constant variance is not assumed, so be careful on the choice of residual plots.
\item The response residuals do not have constant variance.
\item The deviance residuals are expected to have  constant variance.
\item Plot the deviance residuals against the fitted linear predictor $\hat\eta$.
\item however, residual plots are not helpful in some cases; see  p.127 for reasons.
\item Normality is not assumed except for the Gaussian.
\item Use half-normal plot to check {\bf unusual points} (such as outliers or influential cases)
\ei

\para{Difficulty:} How to judge the {\bf curvilinear} relationship between the response and the predictors? 
\bi
\item  Use the {\bf linearized response} 
\[ z = \eta +(y-\mu) g'(\mu) \]
\item plot $z$ versus linear  predictor $\hat\eta$ to check link function $\eta=g(\mu)$ 
 \item partial residual plot:  $z-\eta + \hat\beta_j x_j$ versus $x_j$
\ei

\para{Example}  data(gala)

\end{document}


\end{document}


 