% pdflatex mult_regSL.tex ; pdflatex mult_regHO.tex
% HANDOUTS - use this line
%\documentclass[handout,mathserif]{beamer}
\usepackage{graphics,amsmath,amsfonts,amssymb}
\usepackage[T1]{fontenc}
\usepackage[numbers, sort&compress]{natbib}
\usepackage[utf8]{inputenc}
\usepackage{float, rotating, subfigure}
%\usepackage{epsfig}
\usepackage{color}
\usepackage{xcolor}
%\usepackage[usenames,dvipsnames,svgnames,table]{xcolor}
\usepackage{hyperref}
\usepackage{multicol}
\usepackage{hyperref}
\usepackage{bm}

\usepackage{manfnt}% for \textdbend

\makeatletter
\makeatother

\definecolor{links}{HTML}{2A1B81}
\hypersetup{colorlinks,linkcolor=,urlcolor=links}
%\hypersetup{pdfpagemode=FullScreen}

\definecolor{OliveGreen}{cmyk}{0.64,0,0.95,0.40}
\definecolor{MyGrey}{cmyk}{0.2,0.2,0.2,0.5}
%\usecolortheme[named=MyGrey]{structure} 
%\usecolortheme[named=OliveGreen]{structure}
\usetheme{CambridgeUS}
%\usecolortheme{CambridgeUS}
%\usetheme{Boadilla}
\newcommand{\ggr}{\color{red}}

%\pgfpagesuselayout{4 on 1}[a4paper,border shrink = 5mm, landscape]
\beamertemplatenavigationsymbolsempty % removes pdf-content, for printing


%\setbeamercolor{block}{fg=Red}
%\setbeamercolor{block}{bg=gray}
\setbeamercolor{math text}{fg=OliveGreen}


\title[Stat. modelling]{Statistical modelling: theory and practice}
\subtitle{Multilple linear regression }
\author{Gilles Guillot}
\institute{gigu@imm.dtu.dk}

\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
 \maketitle
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
 \tableofcontents
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Introductory example}


%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}{Example I: The Pine Processionary Caterpillar}

\hspace{3.5cm}\includegraphics[height=3.5cm]{./fig/pin.jpeg}
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\hspace{3cm}\includegraphics[height=5cm]{./fig/nest.jpg}
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\hspace{4cm}\includegraphics[height=5cm]{./fig/cat.jpg}
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\hspace{3cm}\includegraphics[height=4cm]{./fig/proc.jpg}
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{}
{Harmfulness still debated but
  \begin{itemize}
  \item Weaken the trees
  \item Cause strong irritating reactions to humans and pets
  \end{itemize}
}
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{What  environmental factors favour their development?}
\pause Experimental design: 33 forest plots (of 500 sq meters each) in the South of France\\
\pause Response variable: average number of nests per tree in the plot (last column of the file)\\
\pause Explanatory variables (col. 1-10):
\begin{itemize}
\pause \item Altitude (meters) 
\pause \item Slope (degrees)
\pause \item Height of the tree at the centre of the plot
\pause \item Diameter of the tree at the centre of the plot
\pause \item Index of tree density in the plot
\pause \item Orientation (from 1 if southbound to 2)
\pause \item Height of the dominant tree (meters)
\pause \item Number of vegetation strata
\pause \item Mixing index of vegetation cover: from 1 (not mixed) to 2 (mixed)
\end{itemize}
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}{Questions brought up by the caterpillar data}

 \begin{itemize}
  \pause \item How do local environmental conditions affect caterpillar spread?
\pause \item Where should caterpillars be expected?
  \end{itemize}

Data available as \href{http://www2.imm.dtu.dk/courses/02418/slides/w3_multiple_reg/data/caterpillar.dat}{text file caterpillar.dat} on course web page.
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Example II: relation between salary and experience}

A available as \href{http://www.stat.tamu.edu/~sheather/book/docs/datasets/profsalary.txt}{text file profsalary.txt} on Sheather's book web site.

\begin{itemize}
\pause \item EDA and simple model fit with R...

\pause  \item A simple regression does not capture the whole pattern

\pause  \item A model along the line of $$y = a_0 + a_1 x + a_2 x + \varepsilon$$ would be more accurate
 \end{itemize}
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Definition}

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[allowframebreaks]{The multiple regression model}

\begin{block}{Setting:}
 \begin{itemize}
  \pause \item $Y = (y_1,....,y_n)^t$ response variable observed on $n$ individuals\\
  \pause \item $X_1=x_{i1}=(x_{11},...,x_{n1})^t$ first explanatory variable 
  \pause \item $X_2=x_{i1}=(x_{12},...,x_{n2})^t$ second explanatory variable \\
 ...
 \pause \item $X_p=x_{i1}=(x_{1p},...,x_{np})^t$ p-th explanatory variable


  \end{itemize}
All variables are quantitative, preferably continuous
\end{block}


\begin{block}{Model assumptions:}
 \begin{itemize}
  \pause \item The explanatory  variables $x_{ij}$-s are observed and non random
  \pause \item The response variables $y_i$-s are observed and random 
\pause  \item $y_i$ relates linearly to the $x_{ij}$-s:
$$ y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + ... +  \beta_p x_{ip} + \varepsilon_i$$
\pause  \item The coefficients $\beta_0,...\beta_p$ are deterministic (and unknown)
\pause  \item The error terms are independent and indentically distributed (iid)\\
  as a ${\cal N}(0,\sigma^2_e)$
\pause  \item  The error variance $\sigma^2_e$ is deterministic (and unknown)
  \end{itemize}
\end{block}


\begin{block}{Multiple regression equation in matrix form:}

  \pause  Denoting 
 \begin{itemize}  
  \pause \item    $\bm{ \beta} = (\beta_0,\beta_1,...,\beta_p)^t$ 
  \pause \item    $\bm{X} = (\bm{1},X_1,...,X_p)$
 \pause \item    $E = (\varepsilon_1,...,\varepsilon_n)^t$
  \end{itemize}  

 
We have: 
\begin{eqnarray}
\pause Y & = & \bm{X} \bm{\beta} + E \nonumber \\
\pause  ( n \times 1) & & (n \times p+1) (p+1 \times 1) (n \times 1)  \nonumber
\end{eqnarray}
\end{block}
{\footnotesize \pause In the above, $\bm{X}$ is known as the {\em design} matrix. \\
 Stems from the situations where the scientist choses the $x_{ij}$ values in view of maximising 
the amount of information obtained for a given cost (design of experiment or DoE).}




\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Parameter estimation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[allowframebreaks]{Least square estimation}

\pause What parameters?\\
\pause $\bm{ \beta} = (\beta_0,\beta_1,...,\beta_p)^t$  and $\sigma^2_e$

\begin{block}{Least square principle}
We seek  $\bm{ \beta}$ that minimizes globally the squared errors 
$$[y_i - (\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + ... +  \beta_p x_{ip})]^2$$

Our estimate is
$$\hat{\bm{ \beta}} = 
Argmax_{\bm{ \beta}} \sum_i [y_i - (\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + ... +  \beta_p x_{ip})]^2
$$ 
\end{block}


\begin{block}{Geometric interpretation:}
 \begin{itemize}
  \pause \item $\beta_1 \bm{1} + \beta_1 X_1 + .... + \beta_p X_p$ 
is an element of  $Span( \bm{1},X_1,...,X_p)$
  \pause \item $\sum_i [y_i - (\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + ... +  \beta_p x_{ip})]^2 = ||Y-\bm{X} \bm{\beta}||^2$
 \pause \item $\bm{X} \hat{\bm{\beta}}$ is the orthogonal projection of $Y$ on  $Span( \bm{1},X_1,...,X_p)$
  \end{itemize}
\end{block}

\begin{block}{Deriving $\hat{\bm{ \beta}}$ explictly: the normal equations}
It can be shown that  $\hat{\bm{ \beta}}$ satisfies:
$$(\bm{X}^t\bm{X})\hat{\bm{ \beta}} = \bm{X}^t y $$
which is known as the {\em normal equations}. \\
\end{block}
\pause\begin{block}{The LSE estimate of $\bm{ \beta}$}
 If $(\bm{X}^t\bm{X})$ if of full rank, normal equations give:
$$ \hat{\bm{ \beta}} =(\bm{X}^t\bm{X})^{-1} \bm{X}^t y$$ 

$ \hat{\bm{ \beta}}$ defined above is an {\em unbiased} estimate of $\beta$.

\end{block}
\bigskip
Note: the LSE estimate exists if the columns of $X$  are linearly independent.


\pagebreak
\begin{block}{LSE estimate of the residual variance $\sigma^2$:}

$$S^2 = \frac{RSS}{n-p-1} = \frac{1}{n-p-1} \sum_i (\hat{y}_i-y_i)^2$$
 is an unbiased estimate of $\sigma^2$.
\end{block}


\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}{Model fitting in R}
\begin{block}{R function {\tt lm}:}
 \begin{itemize}
\pause \item {\tt lm(y \~{} 1 + x1 + x2 + x3)}\\

\pause \item {\tt lm(y \~{} x1 + x2 + x3) \# includes a constant term by default}\\

\pause \item {\tt lm(y \~{} 0 +x1 + x2 + x3) \# no constant term }\\


\pause \item {\tt lm(y \~{} .,data=mydata) \\
\# include all variables in the rhs   from data.frame mydata}

\pause\item   {\tt lm(y \~{} . - x2 ,data=mydata) \# drop variable x2}
\end{itemize}
 \end{block}

By default, the {\tt lm} function does much more than just parameter estimation.  


\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Testing}

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[allowframebreaks]{Testing global significance}

We are interested in testing whether the response $y$ is pure noise and does not relate 
linearly to any of the putative explanatory variables. 
This is a global test of significance and answers to the question ``does the model explain anything at all?''

\pause \begin{block}{Null hypothesis in global significance testing}
The null hypthesis is $H_0: \beta_1=\beta_2=....=\beta_p=0$ i.e 
 $H_0: Y = E$.
\end{block}

\pagebreak
We define 
 \begin{itemize}
  \pause \item $SST = \sum_i(y_i-\bar{y})^2$
  \pause \item $SSM = \sum_i (\hat{y}_i-\bar{y})^2$
 \pause \item  $SSE = \sum_i (\hat{y}_i-y_i)^2$
  \end{itemize}

If one of the parameters $\beta_j$ is $\neq0$, the model should ``do good'' with 
$SSM \approx SST$ and $SSE \approx 0$.

\pause \begin{block}{Test statistic in global significance testing}
Under $H_0$, $ \displaystyle  F = \frac{SSM/p}{SSE/(n-p-1)} \sim F(p,n-p-1)$.
\pause  A large value of $F$ suggests that the model 
``explains something'' and leads to reject $H_0$.
\end{block}

\pagebreak

In R, the p-value for the global F test can be found at the first line 
of the analysis of variance table ({\tt anova(lm(formula))}) or on the last line 
of {\tt summary(lm(formula))}.
\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[allowframebreaks]{Testing individual coefficients}
If $H_0$ is rejected in the global test above, we want to identify 
wich coefficient $\beta_j$ is non zero. \\
Mor generally, we can be interested in testing $H_0:\beta_j=\beta_j^0$ for some specific 
$\beta_j^0$ values. 

\pause \begin{block}{Test for $\beta_j$ }
Under $H_0:\beta_j=\beta_j^0$, 
$$\displaystyle T_j = \frac{\beta_j-\beta_j^0}{\widehat{sd(\hat{\beta_j})}} \sim St_{n-p-1}$$
\end{block}

\pagebreak
The previous test of $\beta_j=0$ is performed under the assumption that other parameters are non zero. It measures the improvement brought by $X_j$ in a model \underline{containing the other variables}.

\bigskip 

\pause The result of the test of $\beta_j=0$ depends on which other variables 
are included in the model. Cf. cubic term for salaries with and without quadratic term.
 

\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Variable selection}


%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[allowframebreaks]{Model assessment and variable selection}

\pause \begin{block}{Coefficient of determination $R^2$}
Same definition as for simple linear regression:
$$R^2 = 1- \frac{SSE}{SST}$$
\end{block}

\bigskip 
\pause {\ggr \textdbend} An important word of caution:
$$ R^2_{X_1,....,X_p} \leq R^2_{X_1,....,X_p,X_{p+1}} $$

Including an extra variable will always increase {\em automatically} $R^2$.

\pause  \begin{block}{Adjusted coefficient of determination $R^2_{adj}$:}
\pause  A coefficient that includes a penalty term in $p$
$$R^2_{adj} = 1- \frac{SSE/(n-p-1)}{SST/(n-1)}$$
\end{block}
\bigskip 
$R^2_{adj}$ makes it possible (in principle) to compare models of different dimensions.

\bigskip
\pause Other criteria for model selection include Akaike Information Criterion (AIC) and 
Bayeisan Information Criterion  (BIC). 
To be defined after chapters on likelihood and Bayesian theory. 


\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%




\section{References}
%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Suggested reading}
Sheather, chapter 5. \href{http://www.springerlink.com.globalproxy.cvt.dk/content/pj526800177g7186/fulltext.pdf}{Multiple linear regression}.
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Exercise: the caterpillar case}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[fragile,allowframebreaks]{Caterpillar case: suggested steps}


{\footnotesize
  \begin{itemize}
  \item Load the data and give the variables some names, e.g.
\begin{verbatim}
cat.dat <- read.table("...")

colnames(cat.dat) <- c("Altitude",
              "Slope",
              "Nb.pines.in.area",
              "Height.tree.center",
              "Diameter.tree.center",
              "Density.index",
              "Orientation",
              "Height.dominant.tree",
              "Nb.vegetation.strata",
              "Mixing.index",
              "Nb.nests")
\end{verbatim}

\item  Explore graphically the dataset (the R functions {\tt plot, pairs} and {\tt hist} are your best friends).\\
  You may want to print the correlation coefficients on the pairs plot with:\\
\begin{verbatim}
panel.cor <- function(x, y, digits=2, prefix="", cex.cor, ...)
     {
         usr <- par("usr"); on.exit(par(usr))
         par(usr = c(0, 1, 0, 1))
         r <- abs(cor(x, y))
         txt <- format(c(r, 0.123456789), digits=digits)[1]
         txt <- paste(prefix, txt, sep="")
         if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
         text(0.5, 0.5, txt, cex = cex.cor * r)
     }

pairs(cat.dat, lower.panel=panel.smooth, upper.panel=panel.cor)
\end{verbatim}
  \end{itemize}
}

%% \end{frame}


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% \begin{frame}[fragile]
%% \frametitle{}

{\footnotesize
\begin{itemize}
\item The raw correlation coefficient can be misleading as some of the explanatory variables are correlated with each other.
One can ``get rid'' of the effect of one variable using the partial correlation coefficient. It is defined as the correlation coefficent 
bewteen the residuals of a regression on the variable one wants to get rid of. It is obtained as:\\
\begin{verbatim}
pcor <- function(v1, v2, v3)
       {
       c12 <- cor(v1, v2)
       c23 <- cor(v2, v3)
       c13 <- cor(v1, v3)
       partial <- (c12-(c13*c23))/(sqrt(1-(c13^2)) * sqrt(1-(c23^2)))
       return(partial)
       }
\end{verbatim}
Use this to evaluate the correlation between the response variable and the explanatory 
variables given {\tt Nb.vegetation.strata}. Compare to the plain correlation. 
Which variable would you expect to be useless in a multiple regression once {\tt Nb.vegetation.strata} 
is included in the model?
\end{itemize}
}

%% \end{frame}



%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% \begin{frame}[fragile]
%% \frametitle{}
{\footnotesize
\begin{itemize}
\item Fit a linear model explaining the number of nests with all the other variables as predictors.\\
%Hints: {\tt lm.all = lm(formula = Nb.nests \~ \; . , data=cat.dat)}\\
Give the exact math equation that describes the model fitted by the R command above. What are the parameters involved? 
What are the main assumptions involved?
\item Can you reject with confidence the null hypothesis in a global test?
\item Plot the residuals againtst the predicted values. Does it seem reasonable to assume that the variance of the noise 
is the same for all observations?\\
Hints: The residuals are computed by the  {\tt lm } function and stored as a component of the returned list named residuals. 
The predicted values can be obtained by {\tt predict(lm.all)}
\item Can the model quality be improved by a regression of the logarithm of the number of nests (instead of the number of nests)?
\end{itemize}
}

%% \end{frame}


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% \begin{frame}[fragile]
%% \frametitle{}
{\small From now on, we work with the logarithm of the number of nests as response variable.\\}


{\footnotesize
\begin{itemize}
\item The last line printed by the R command {\tt summary(lm.all)} contains the p-value for the null hypothesis of  absence 
of correlation with any explanatoty variable. Can you accept safely this null hypothesis?
\item Plot the predicted values against the observed response values and compute their coefficient of correlation. 
Does this seem to be a useful model?
\item Print the esimated coefficients and associated p-values. What variables stand out as good or poor predictors? 
Compare this to a linear model with Density.index as single predictor. How do you explain the difference?
\end{itemize}
}

%% \end{frame}



%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% \begin{frame}[fragile]
%% \frametitle{}

{\tiny
\begin{itemize}
\item From the above, there are obsviously some variables that are useless in the model.
Compute leave-one-out prediction error with all variables then with all variables except Mixing.index. \\

Hints:\\
\begin{verbatim}
## leave one out prediction with full model
mse.1 <- 0
loo.pred1 <- numeric(nrow(cat.dat))
for(i in 1:nrow(cat.dat))
  {
    ## fitting a model with all variables
    ## to the dataset without obs. i
    lm.res.loo <- lm(formula=log(Nb.nests) ~. , data=cat.dat[-i,])
    ## predicting Log.nb.nests[i] with this model
    ## (just a linear combination with estimated coef.):
    loo.pred1[i] <- sum(c(1,as.numeric(cat.dat[i,1:10]))*lm.res.loo$coefficients)
    ## squre error in prediction
    mse.1 <- mse.1 + (log(cat.dat$Nb.nests[i]) - loo.pred1[i])^2
  }

mse.1 <- mse.1/(nrow(cat.dat)-1)
rmse.1 <- sqrt(mse.1)



\end{verbatim}
\end{itemize}
}

%% \end{frame}

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% \begin{frame}[fragile]
%% \frametitle{}

{\tiny

\begin{verbatim}
## leave one out prediction with all variables exept Mixing.index
mse.2 <- 0
loo.pred2 <- numeric(nrow(cat.dat))
for(i in 1:nrow(cat.dat))
  {
    ## fitting a model with all variables except Mixing.index
    ## to the dataset without obs. i
    lm.res.loo <- lm(formula=(log(Nb.nests) ~.  - Mixing.index), data=cat.dat[-i,])
    ## predicting Log.nb.nests[i] with this model
    ## (just a linear combination with estimated coef.):
    loo.pred2[i] <- sum(c(1,as.numeric(cat.dat[i,1:9]))*lm.res.loo$coefficients)
    ## squre error in prediction
    mse.2 <- mse.2 + (log(cat.dat$Nb.nests[i]) - loo.pred2[i])^2
  }
mse.2 <- mse.2/(nrow(cat.dat)-1)
rmse.2 <- sqrt(mse.2)

rmse.1
rmse.2

##
plot(log(cat.dat$Nb.nests),loo.pred1,xlab="Log nb  of nest",
     ylab="LOO prediction"); abline(0,1, col=3)
points(log(cat.dat$Nb.nests),loo.pred2,col=2,pch=2)


\end{verbatim}
On the basis of the above, would you include {\tt Mixing.index} in the model?

}



%% {\footnotesize
%% \begin{itemize}
%% %% \item We want now to perform variable selection. We will minimize the  Akaike information criterion (AIC) defined e.g. in [Bingham and Fry p. 140].\\
%% %% This is performed in R with the function {\tt step()} that takes the output of the {\tt lm ()} function as an argument.\\
%% %% What is the final model retained? 

%% \item What leave-one-out  Root Mean Square Error (RMSE) do you obtain with this ``best'' model?

%% \end{itemize}
%% }

\end{frame}




\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%





%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}


\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}[allowframebreaks]{}


\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
  \frametitle{}
  \begin{itemize}
  \pause \item
  \pause \item
  \end{itemize}
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{block}{}
\end{block}
