\documentclass{chapman}

%%% copy Sweave.sty definitions

%%% keeps `sweave' from adding `\usepackage{Sweave}': DO NOT REMOVE
%\usepackage{Sweave} 


\RequirePackage[T1]{fontenc}
\RequirePackage{graphicx,ae,fancyvrb}
\IfFileExists{upquote.sty}{\RequirePackage{upquote}}{}
\usepackage{relsize}

\DefineVerbatimEnvironment{Sinput}{Verbatim}{}
\DefineVerbatimEnvironment{Soutput}{Verbatim}{fontfamily=courier,
                                              fontshape=it,
                                              fontsize=\relsize{-1}}
\DefineVerbatimEnvironment{Scode}{Verbatim}{}
\newenvironment{Schunk}{}{}

%%% environment for raw output
\newcommand{\SchunkRaw}{\renewenvironment{Schunk}{}{}
    \DefineVerbatimEnvironment{Soutput}{Verbatim}{fontfamily=courier,
                                                  fontshape=it,
                                                  fontsize=\small}
    \rawSinput
}

%%% environment for labeled output
\newcommand{\nextcaption}{}
\newcommand{\SchunkLabel}{
  \renewenvironment{Schunk}{\begin{figure}[ht] }{\caption{\nextcaption}
  \end{figure} }
  \DefineVerbatimEnvironment{Sinput}{Verbatim}{frame = topline}
  \DefineVerbatimEnvironment{Soutput}{Verbatim}{frame = bottomline, 
                                                samepage = true,
                                                fontfamily=courier,
                                                fontshape=it,
                                                fontsize=\relsize{-1}}
}


%%% S code with line numbers
\DefineVerbatimEnvironment{Sinput}
{Verbatim}
{
%%  numbers=left
}

\newcommand{\numberSinput}{
    \DefineVerbatimEnvironment{Sinput}{Verbatim}{numbers=left}
}
\newcommand{\rawSinput}{
    \DefineVerbatimEnvironment{Sinput}{Verbatim}{}
}


%%% R / System symbols
\newcommand{\R}{\textsf{R}}
\newcommand{\rR}{{R}}
\renewcommand{\S}{\textsf{S}}
\newcommand{\SPLUS}{\textsf{S-PLUS}}
\newcommand{\rSPLUS}{{S-PLUS}}
\newcommand{\SPSS}{\textsf{SPSS}}
\newcommand{\EXCEL}{\textsf{Excel}}
\newcommand{\ACCESS}{\textsf{Access}}
\newcommand{\SQL}{\textsf{SQL}}
%%\newcommand{\Rpackage}[1]{\hbox{\rm\textit{#1}}}
%%\newcommand{\Robject}[1]{\hbox{\rm\texttt{#1}}}
%%\newcommand{\Rclass}[1]{\hbox{\rm\textit{#1}}}
%%\newcommand{\Rcmd}[1]{\hbox{\rm\texttt{#1}}}
\newcommand{\Rpackage}[1]{\index{#1 package@{\fontseries{b}\selectfont #1} package} {\fontseries{b}\selectfont #1}}
\newcommand{\rpackage}[1]{{\fontseries{b}\selectfont #1}}
\newcommand{\Robject}[1]{\texttt{#1}}
\newcommand{\Rclass}[1]{\index{#1 class@\textit{#1} class}\textit{#1}}
\newcommand{\Rcmd}[1]{\index{#1 function@\texttt{#1} function}\texttt{#1}}
\newcommand{\Roperator}[1]{\texttt{#1}}
\newcommand{\Rarg}[1]{\texttt{#1}}
\newcommand{\Rlevel}[1]{\texttt{#1}}


%%% other symbols
\newcommand{\file}[1]{\hbox{\rm\texttt{#1}}}
%%\newcommand{\stress}[1]{\index{#1}\textit{#1}} 
\newcommand{\stress}[1]{\textit{#1}} 
\newcommand{\booktitle}[1]{\textit{#1}} %%'

%%% Math symbols
\usepackage{amstext}
\usepackage{amsmath}
\newcommand{\E}{\mathsf{E}}   
\newcommand{\Var}{\mathsf{Var}}   
\newcommand{\Cov}{\mathsf{Cov}}   
\newcommand{\Cor}{\mathsf{Cor}}   
\newcommand{\x}{\mathbf{x}}   
\newcommand{\y}{\mathbf{y}}   
\renewcommand{\a}{\mathbf{a}}
\newcommand{\W}{\mathbf{W}}   
\newcommand{\C}{\mathbf{C}}   
\renewcommand{\H}{\mathbf{H}}   
\newcommand{\X}{\mathbf{X}}   
\newcommand{\B}{\mathbf{B}}   
\newcommand{\V}{\mathbf{V}}   
\newcommand{\I}{\mathbf{I}}   
\newcommand{\D}{\mathbf{D}}   
\newcommand{\bS}{\mathbf{S}}   
\newcommand{\N}{\mathcal{N}}   
\renewcommand{\P}{\mathsf{P}}   
\newcommand{\K}{\mathbf{K}}
\newcommand{\m}{\mathbf{m}}
\newcommand{\argmin}{\operatorname{argmin}\displaylimits}
\newcommand{\argmax}{\operatorname{argmax}\displaylimits}
%%% links
\usepackage{hyperref}

\hypersetup{%
  pdftitle = {A Handbook of Statistical Analyses Using R},
  pdfsubject = {Book},
  pdfauthor = {Brian S. Everitt and Torsten Hothorn},
  colorlinks = {black},
  linkcolor = {black},
  citecolor = {black},
  urlcolor = {black},
  hyperindex = {true},
  linktocpage = {true},
}


%%% captions & tables
%% <FIXME>: conflics with figure definition in chapman.cls
%%\usepackage[format=hang,margin=10pt,labelfont=bf]{caption}
%% </FIMXE>
\usepackage{longtable}
\usepackage[figuresright]{rotating}

%%% R symbol in chapter 1
\usepackage{wrapfig}

%%% Bibliography
\usepackage[round,comma]{natbib}
\renewcommand{\refname}{References \addcontentsline{toc}{chapter}{References}}
\citeindexfalse

%%% texi2dvi complains that \newblock is undefined, hm...
\def\newblock{\hskip .11em plus .33em minus .07em}

%%% Example sections
\newcounter{exercise}[chapter]
\setcounter{exercise}{0}
\newcommand{\exercise}{\item{\stepcounter{exercise} Ex.
                       \arabic{chapter}.\arabic{exercise} }}


%% URLs
\newcommand{\curl}[1]{\begin{center} \url{#1} \end{center}}

%%% for manual corrections
%\renewcommand{\baselinestretch}{2}

%%% plot sizes
\setkeys{Gin}{width=0.95\textwidth}

%%% color
\usepackage{color}

%%% hyphenations
\hyphenation{drop-out}
\hyphenation{mar-gi-nal}

%%% new bidirectional quotes need 
\usepackage[utf8]{inputenc}
\begin{document}

%% Title page

\title{A Handbook of Statistical Analyses Using \R{} --- 2nd Edition}

\author{Brian S. Everitt and Torsten Hothorn}

\maketitle
%%\VignetteIndexEntry{Chapter Logistic Regression and Generalised Linear Models}
\setcounter{chapter}{6}


\SweaveOpts{prefix.string=figures/HSAUR,eps=FALSE,keep.source=TRUE} 

<<setup, echo = FALSE, results = hide>>=
rm(list = ls())
s <- search()[-1]
s <- s[-match(c("package:base", "package:stats", "package:graphics", "package:grDevices",
                "package:utils", "package:datasets", "package:methods", "Autoloads"), s)]
if (length(s) > 0) sapply(s, detach, character.only = TRUE)
if (!file.exists("tables")) dir.create("tables")
if (!file.exists("figures")) dir.create("figures")
set.seed(290875)
options(prompt = "R> ", continue = "+  ",
    width = 63, # digits = 4, 
    show.signif.stars = FALSE,
    SweaveHooks = list(leftpar = function() 
        par(mai = par("mai") * c(1, 1.05, 1, 1)),
        bigleftpar = function()
        par(mai = par("mai") * c(1, 1.7, 1, 1))))
HSAURpkg <- require("HSAUR2")
if (!HSAURpkg) stop("cannot load package ", sQuote("HSAUR2"))
rm(HSAURpkg)
 ### </FIXME> hm, R-2.4.0 --vanilla seems to need this
a <- Sys.setlocale("LC_ALL", "C")
 ### </FIXME>
book <- TRUE
refs <- cbind(c("AItR", "DAGD", "SI", "CI", "ANOVA", "MLR", "GLM", 
                "DE", "RP", "GAM", "SA", "ALDI", "ALDII", "SIMC", "MA", "PCA", 
                "MDS", "CA"), 1:18)
ch <- function(x) {
    ch <- refs[which(refs[,1] == x),]
    if (book) {
        return(paste("Chapter~\\\\ref{", ch[1], "}", sep = ""))
    } else {
        return(paste("Chapter~", ch[2], sep = ""))
    }
}
if (file.exists("deparse.R"))
    source("deparse.R")

setHook(packageEvent("lattice", "attach"), function(...) {
    lattice.options(default.theme = 
        function()
            standard.theme("pdf", color = FALSE))
    })
@

\pagestyle{headings}
<<singlebook, echo = FALSE>>=
book <- FALSE
@

\chapter[Logistic Regression and Generalised Linear Models]{Logistic
Regression and Generalised Linear Models: Blood Screening, Women's Role in %'
Society, Colonic Polyps, and Driving and Back Pain\label{GLM}}

\section{Introduction}



\section{Logistic Regression and Generalised Linear Models}


\section{Analysis Using \R{}}

\subsection{ESR and Plasma Proteins}


\begin{figure}
\begin{center}
<<GLM-plasma-plot, echo = TRUE, fig = TRUE, height = 4>>=
data("plasma", package = "HSAUR2")
layout(matrix(1:2, ncol = 2))
cdplot(ESR ~ fibrinogen, data = plasma)
cdplot(ESR ~ globulin, data = plasma)
@
\caption{Conditional density plots of 
 the erythrocyte sedimentation rate (ESR) given fibrinogen and globulin.
 \label{GLM:plasma1}}
\end{center}
\end{figure}

We can now fit a logistic regression model to the data using 
the \Rcmd{glm} function. We start with a model that includes 
only a single explanatory variable, \Robject{fibrinogen}. The 
code to fit the model is
<<GLM-plasma-fit1, echo = TRUE>>=
plasma_glm_1 <- glm(ESR ~ fibrinogen, data = plasma, 
                    family = binomial())
@
The formula implicitly defines a parameter for the global mean (the
intercept term) as discussed in \Sexpr{ch("ANOVA")} and \Sexpr{ch("MLR")}. The
distribution of the response is defined by the \Robject{family} argument, a
binomial distribution in our case.
\index{family argument@\Rcmd{family} argument}
\index{Binomial distribution}
(The default link function when the binomial family is requested 
is the logistic function.)



\renewcommand{\nextcaption}{\R{} output of the \Robject{summary} 
                            method for the logistic regression model fitted
                            to ESR and fibrigonen.
                            \label{GLM-plasma-summary-1}}
\SchunkLabel
<<GLM-plasma-summary-1, echo = TRUE>>=
summary(plasma_glm_1)
@
\SchunkRaw


From the results in Figure~\ref{GLM-plasma-summary-1} we 
see that the regression 
coefficient for fibrinogen is significant at the $5\%$ level. 
An increase 
of one unit in this variable increases the log-odds in favour 
of an ESR value greater than $20$ by an estimated
$\Sexpr{round(coef(plasma_glm_1)["fibrinogen"], 2)}$ with 95\% 
confidence interval 
<<GLM-plasma-confint, echo = FALSE, results = hide>>=
ci <- confint(plasma_glm_1)["fibrinogen",]
@
<<GLM-plasma-confint, echo = TRUE, results = hide>>=
confint(plasma_glm_1, parm = "fibrinogen")
@
<<GLM-plasma-confint, echo = FALSE>>=
print(ci) 
@
These values are more helpful 
if converted to the corresponding values for the odds themselves 
by exponentiating the estimate
<<GLM-plasma-exp, echo = TRUE>>=
exp(coef(plasma_glm_1)["fibrinogen"])
@
and the confidence interval 
<<GLM-plasma-exp-ci, echo = FALSE, results = hide>>=
ci <- exp(confint(plasma_glm_1, parm = "fibrinogen"))
@
<<GLM-plasma-exp-ci, echo = TRUE, results = hide>>=
exp(confint(plasma_glm_1, parm = "fibrinogen"))
@
<<GLM-plasma-exp-ci, echo = FALSE>>=
print(ci)
@
The confidence interval is very wide because there are few observations 
overall and very few where the ESR value is greater than $20$. 
Nevertheless it seems likely that increased values of fibrinogen 
lead to a greater probability of an ESR value greater than $20$.

We can now fit a logistic regression model that includes both 
explanatory variables using the code
<<GLM-plasma-fit2, echo = TRUE>>=
plasma_glm_2 <- glm(ESR ~ fibrinogen +  globulin, 
    data = plasma, family = binomial())
@
and the output of the \Rcmd{summary} method is shown in Figure
\ref{GLM-plasma-summary-2}.

\renewcommand{\nextcaption}{\R{} output of the \Robject{summary} 
                            method for the logistic regression model fitted
                            to ESR and both globulin and fibrinogen.
                            \label{GLM-plasma-summary-2}}
\SchunkLabel
<<GLM-plasma-summary-2, echo = TRUE>>=
summary(plasma_glm_2)
@
\SchunkRaw
<<GLM-plasma-anova-hide, echo = FALSE, results = hide>>=
plasma_anova <- anova(plasma_glm_1, plasma_glm_2, test = "Chisq")
@

The coefficient for gamma 
globulin is not significantly different from zero. Subtracting 
the residual deviance of the second model from the corresponding 
value for the first model we get a value of
$\Sexpr{round(plasma_anova$Deviance[2], 2)}$. Tested using a  
$\chi^2$-distribution with a single degree of freedom this is not significant 
at the 5\% level and so we conclude that gamma globulin is not 
associated with ESR level. In \R{}, the task of comparing the two nested
models can be performed using the \Rcmd{anova} function
<<GLM-plasma-anova, echo = TRUE>>=
anova(plasma_glm_1, plasma_glm_2, test = "Chisq")
@
Nevertheless we shall use the predicted 
values from the second model and plot them against the values 
of \stress{both} explanatory variables using a \stress{bubbleplot} to illustrate the use of the  \Rcmd{symbols} function. 
\index{Bubbleplot}
The estimated conditional probability of a ESR value larger $20$ for all
observations can be computed, following formula (\ref{GLM:logitexp}), by
<<GLM-plasma-predict, echo = TRUE>>=
prob <- predict(plasma_glm_2, type = "response")
@
and now we can assign a larger circle to observations with larger probability
as shown in Figure~\ref{GLM:bubble}. The plot clearly 
shows the increasing probability of an ESR value above $20$ (larger 
circles) as the values of fibrinogen, and to a lesser extent, 
gamma globulin, increase.

\begin{figure}
\begin{center}
<<GLM-plasma-bubble, echo = TRUE, fig = TRUE>>=
plot(globulin ~ fibrinogen, data = plasma, xlim = c(2, 6), 
     ylim = c(25, 55), pch = ".")
symbols(plasma$fibrinogen, plasma$globulin, circles = prob,
        add = TRUE)
@
\caption{Bubbleplot of fitted values for a logistic regression model
fitted to the \Robject{plasma} data. \label{GLM:bubble}}     
\end{center}
\end{figure}


\subsection{Women's Role in Society} %'

Originally the data in Table~\ref{GLM-womensrole-tab} 
would have been in a completely 
equivalent form to the data in Table~\ref{GLM-plasma-tab} 
data, but here 
the individual observations have been grouped into counts 
of numbers of agreements and disagreements for the two explanatory 
variables, \Robject{gender} and \Robject{education}.
To fit a logistic 
regression model to such grouped data using the \Rcmd{glm} function 
we need to specify the number of agreements and disagreements 
as a two-column matrix on the left hand side of the model formula. 
We first fit a model that includes the two explanatory variables 
using the code
<<GLM-womensrole-fit1, echo = TRUE>>=
data("womensrole", package = "HSAUR2")
fm1 <- cbind(agree, disagree) ~ gender + education
womensrole_glm_1 <- glm(fm1, data = womensrole, 
                        family = binomial())
@

\renewcommand{\nextcaption}{\R{} output of the \Robject{summary}
                            method for the logistic regression model fitted
                            to the \Robject{womensrole} data.
                            \label{GLM-womensrole-summary-1}}
\SchunkLabel
<<GLM-womensrole-summary-1, echo = TRUE>>=
summary(womensrole_glm_1)
@
\SchunkRaw

From the \Rcmd{summary} output in Figure~\ref{GLM-womensrole-summary-1} 
it appears that education 
has a highly significant part to play in predicting whether a respondent 
will agree with the statement read to them, but the respondent's %'
gender is apparently unimportant. As years of education increase 
the probability of agreeing with the statement declines. 
We now are going to construct a plot comparing the observed proportions of
agreeing with those fitted by our fitted model. Because we will reuse this
plot for another fitted object later on, we define a function which plots
years of education against some fitted probabilities, e.g., 
<<GLM-womensrole-probfit, echo = TRUE>>=
role.fitted1 <- predict(womensrole_glm_1, type = "response")
@
and labels each observation with the person's gender: %%'
\numberSinput
<<GLM-plot-setup, echo = TRUE>>=
myplot <- function(role.fitted)  {
    f <- womensrole$gender == "Female"
    plot(womensrole$education, role.fitted, type = "n", 
         ylab = "Probability of agreeing",
         xlab = "Education", ylim = c(0,1))
    lines(womensrole$education[!f], role.fitted[!f], lty = 1)
    lines(womensrole$education[f], role.fitted[f], lty = 2)  
    lgtxt <- c("Fitted (Males)", "Fitted (Females)")
    legend("topright", lgtxt, lty = 1:2, bty = "n")
    y <-  womensrole$agree / (womensrole$agree + 
                              womensrole$disagree)
    text(womensrole$education, y, ifelse(f, "\\VE", "\\MA"), 
         family = "HersheySerif", cex = 1.25)
} 
@
\rawSinput

\begin{figure}
\begin{center}
<<GLM-role-fitted1, echo = TRUE, fig = TRUE>>=
myplot(role.fitted1)
@
\caption{Fitted (from \Robject{womensrole\_glm\_1}) and observed probabilities of
         agreeing for the \Robject{womensrole} data. \label{GLM-role1plot}}
\end{center}
\end{figure}
In lines 3--5 of function \Rcmd{myplot}, an empty scatterplot of education and
fitted probabilities (\Rcmd{type = "n"}) is set up, basically to set the scene
for the following plotting actions. Then, two lines are drawn (using function
\Rcmd{lines} in lines 6 and 7), one for males (with line type 1) 
and one for females (with line type 2, i.e., a dashed line), where
the logical vector \Robject{f} describes both genders. In line 9 a legend
is added. Finally, in lines 12 and 13 we plot `observed' values, i.e., 
the frequencies of agreeing in each of the groups (\Robject{y} as computed
in lines 10 and 11) and use the Venus and Mars symbols to indicate gender.

The two curves 
for males and females in Figure~\ref{GLM-role1plot} 
are almost the same reflecting the non-significant 
value of the regression coefficient for gender in
\Robject{womensrole\_glm\_1}. But the observed values plotted on 
Figure~\ref{GLM-role1plot} suggest that 
there might be an interaction of education and gender, a possibility 
that can be investigated by applying a further logistic regression 
model using
\index{Interaction}
<<GLM-womensrole-fit2, echo = TRUE>>=   
fm2 <- cbind(agree,disagree) ~ gender * education
womensrole_glm_2 <- glm(fm2, data = womensrole, 
                        family = binomial())
@
The \Robject{gender} and \Robject{education}
interaction term is seen to be highly significant, as can be seen from the
\Rcmd{summary} output in Figure~\ref{GLM-womensrole-summary-2}.
\renewcommand{\nextcaption}{\R{} output of the \Robject{summary}
                            method for the logistic regression model fitted
                            to the \Robject{womensrole} data.
                            \label{GLM-womensrole-summary-2}}
\SchunkLabel
<<GLM-womensrole-summary-2, echo = TRUE>>=
summary(womensrole_glm_2)
@
\SchunkRaw


\begin{figure}
\begin{center}
<<GLM-role-fitted2, echo = TRUE, fig = TRUE>>=
role.fitted2 <- predict(womensrole_glm_2, type = "response")
myplot(role.fitted2)
@
\caption{Fitted (from \Robject{womensrole\_glm\_2}) and observed probabilities of
agreeing for the \Robject{womensrole} data. \label{GLM-role2plot}}
\end{center}
\end{figure}


We can obtain a plot of deviance residuals plotted against 
fitted values using the following code above Figure~\ref{GLM:devplot}.
\begin{figure}
\begin{center}
<<GLM-role-plot2, echo = TRUE, fig = TRUE>>=
res <- residuals(womensrole_glm_2, type = "deviance")
plot(predict(womensrole_glm_2), res,
     xlab="Fitted values", ylab = "Residuals", 
     ylim = max(abs(res)) * c(-1,1))
abline(h = 0, lty = 2)
@
\caption{Plot of deviance residuals from logistic
         regression model fitted to the \Robject{womensrole} data.
         \label{GLM:devplot}}
\end{center}
\end{figure}
The residuals fall into a horizontal band between $-2$ and $2$. 
This pattern does not suggest  a poor fit for any particular observation 
or subset of observations.


 
\subsection{Colonic Polyps}

The data on colonic polyps in Table~\ref{GLM-polyps-tab} 
involves \stress{count}
data. We could try to model this using multiple regression 
but there are two problems. The first is that a response that 
is a count can take only positive values, and secondly such a 
variable is unlikely to have a normal distribution. Instead we 
will apply a GLM with a log link function, ensuring that fitted 
values are positive, and a Poisson error distribution, i.e., 
\index{Poisson error distribution}
\index{Poisson regression}
\begin{eqnarray*}
\P(y) = \frac{e^{-\lambda}\lambda^y}{y!}.
\end{eqnarray*}
This type of GLM is often known as \stress{Poisson regression}. 
We can apply the model using
<<GLM-polyps-fit1, echo = TRUE>>=
data("polyps", package = "HSAUR2")
polyps_glm_1 <- glm(number ~ treat + age, data = polyps, 
                    family = poisson())
@
(The default link function when the Poisson family is requested 
is the log function.)

\renewcommand{\nextcaption}{\R{} output of the \Robject{summary}
                            method for the Poisson regression model fitted
                            to the \Robject{polyps} data.
                            \label{GLM-polyps-summary-1}}
\SchunkLabel
<<GLM-polyps-summary-1, echo = TRUE>>=
summary(polyps_glm_1)
@
\SchunkRaw


We can deal with overdispersion by using a procedure known 
as \stress{quasi-likelihood}, 
\index{Quasi-likelihood}
which allows the estimation of 
model parameters without fully knowing the error distribution 
of the response variable. \cite{HSAUR:McCullaghNelder1989} give full 
details of the quasi-likelihood approach. In many respects it 
simply allows for the estimation of $\phi$ 
from the data rather than defining it to be unity for the 
binomial and Poisson distributions. We can apply quasi-likelihood 
estimation to the colonic polyps data using the following \R{} code
<<GLM-polyp-quasi, echo = TRUE>>=
polyps_glm_2 <- glm(number ~ treat + age, data = polyps, 
                    family = quasipoisson())
summary(polyps_glm_2)
@
The regression coefficients 
for both explanatory variables remain significant but their estimated 
standard errors are now much greater than the values given in 
Figure~\ref{GLM-polyps-summary-1}. A possible reason for 
overdispersion in these data 
is that polyps do not occur independently of one another, but 
instead may `cluster' together. %'
\index{Overdispersion|)}


\subsection{Driving and Back Pain}

A frequently used design in medicine is the matched case-control study in which 
each patient suffering from a particular condition of interest included in the 
study is matched to one or more people without the condition.  The most commonly 
used matching variables are age, ethnic group, mental status etc.  A design with $m$ 
controls per case is known as a $1:m$ matched study.  In many cases $m$ will be one, 
and it is the $1:1$ matched study that we shall concentrate on here where we analyse 
the data on low back pain given in Table~\ref{GLM-backpain-tab}.  
To begin we shall describe the form of the logistic model appropriate for 
case-control studies in the simplest case where there is only one binary explanatory 
variable.

With matched pairs data the form of the logistic model involves the probability, $\varphi$, 
that in matched pair number $i$, for a given value of the explanatory variable 
the member of the pair is a case.  Specifically the model is
\begin{eqnarray*}
\text{logit}(\varphi_i) = \alpha_i + \beta x.
\end{eqnarray*}
The odds that a subject with $x=1$ is a case equals $\exp(\beta)$ times the odds 
that a subject with $x=0$ is a case.

The model generalises to the situation where there are $q$ explanatory variables as
\begin{eqnarray*}
\text{logit}(\varphi_i) = \alpha_i + \beta_1 x_1 + \beta_2 x_2 + \dots \beta_q x_q.
\end{eqnarray*}

Typically one $x$ is an explanatory variable of real interest, such as past 
exposure to a risk factor, with the others being used as a form of statistical 
control in addition to the variables already controlled by virtue of using them 
to form matched pairs. This is the case in our back pain example where it is the 
effect of car driving on lower back pain that is of most interest.

The problem with the model above is that the number of parameters increases at 
the same rate as the sample size with the consequence that maximum likelihood 
estimation is no longer viable.  We can overcome this problem if we regard 
the parameters $\alpha_i$ as of little interest and so are willing to forgo 
their estimation.  If we do, we can then create a \stress{conditional likelihood function} 
that will yield maximum likelihood estimators of the coefficients, $\beta_1, \dots, \beta_q$, 
that are consistent and asymptotically normally distributed. The mathematics behind 
this are described in \cite{HSAUR:Collett2003}.

The model can be fitted using the \Rcmd{clogit} function from
package \Rpackage{survival}; the results are shown in Figure~\ref{GLM-backpain-print}.
<<GLM-backpain-clogit, echo = TRUE>>=
library("survival")
backpain_glm <- clogit(I(status == "case") ~ 
    driver + suburban + strata(ID), data = backpain)
@
The response has to be a logical (\Rcmd{TRUE} for cases) and the
\Rcmd{strata} command specifies the matched pairs.

\renewcommand{\nextcaption}{\R{} output of the \Robject{print}
                            method for the conditional logistic regression model fitted
                            to the \Robject{backpain} data.
                            \label{GLM-backpain-print}}
\SchunkLabel
<<GLM-backpain-print, echo = TRUE>>=
print(backpain_glm)
@
\SchunkRaw

The estimate of the odds ratio of a herniated disc occurring
in a driver relative to a nondriver is $\Sexpr{round(exp(coef(backpain_glm)[1]),2)}$
with a $95\%$ confidence interval of 
$\Sexpr{paste("(", paste(round(exp(confint(backpain_glm)[1,]), 2), collapse = ","),")", sep = "")}$.
Conditional on residence we can say that the risk of a herniated disc
occurring in a driver is about twice that of a nondriver.
There is no evidence that where a person lives affects the risk
of lower back pain.

\section{Summary}

Generalised linear models provide a very powerful and flexible framework
for the application of regression models to a variety of non-normal
response variables, for example, logistic regression to binary responses   
and Poisson regression to count data.


\bibliographystyle{LaTeXBibTeX/refstyle}
\bibliography{LaTeXBibTeX/HSAUR}   
\end{document}