\documentclass[]{article}
\usepackage{lmodern}
\usepackage{amssymb,amsmath}
\usepackage{ifxetex,ifluatex}
\usepackage{fixltx2e} % provides \textsubscript
\ifnum 0\ifxetex 1\fi\ifluatex 1\fi=0 % if pdftex
  \usepackage[T1]{fontenc}
  \usepackage[utf8]{inputenc}
\else % if luatex or xelatex
  \ifxetex
    \usepackage{mathspec}
  \else
    \usepackage{fontspec}
  \fi
  \defaultfontfeatures{Ligatures=TeX,Scale=MatchLowercase}
\fi
% use upquote if available, for straight quotes in verbatim environments
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
% use microtype if available
\IfFileExists{microtype.sty}{%
\usepackage{microtype}
\UseMicrotypeSet[protrusion]{basicmath} % disable protrusion for tt fonts
}{}
\usepackage[margin=1in]{geometry}
\usepackage{hyperref}
\hypersetup{unicode=true,
            pdftitle={TMA4315 Generalized linear models H2018},
            pdfauthor={Mette Langaas, Department of Mathematical Sciences, NTNU -- with contributions from Ingeborg Hem},
            pdfborder={0 0 0},
            breaklinks=true}
\urlstyle{same}  % don't use monospace font for urls
\usepackage{color}
\usepackage{fancyvrb}
\newcommand{\VerbBar}{|}
\newcommand{\VERB}{\Verb[commandchars=\\\{\}]}
\DefineVerbatimEnvironment{Highlighting}{Verbatim}{commandchars=\\\{\}}
% Add ',fontsize=\small' for more characters per line
\usepackage{framed}
\definecolor{shadecolor}{RGB}{248,248,248}
\newenvironment{Shaded}{\begin{snugshade}}{\end{snugshade}}
\newcommand{\KeywordTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{\textbf{#1}}}
\newcommand{\DataTypeTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{#1}}
\newcommand{\DecValTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\BaseNTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\FloatTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\ConstantTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\CharTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\SpecialCharTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\StringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\VerbatimStringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\SpecialStringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\ImportTok}[1]{#1}
\newcommand{\CommentTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textit{#1}}}
\newcommand{\DocumentationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\AnnotationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\CommentVarTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\OtherTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{#1}}
\newcommand{\FunctionTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\VariableTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\ControlFlowTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{\textbf{#1}}}
\newcommand{\OperatorTok}[1]{\textcolor[rgb]{0.81,0.36,0.00}{\textbf{#1}}}
\newcommand{\BuiltInTok}[1]{#1}
\newcommand{\ExtensionTok}[1]{#1}
\newcommand{\PreprocessorTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textit{#1}}}
\newcommand{\AttributeTok}[1]{\textcolor[rgb]{0.77,0.63,0.00}{#1}}
\newcommand{\RegionMarkerTok}[1]{#1}
\newcommand{\InformationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\WarningTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\AlertTok}[1]{\textcolor[rgb]{0.94,0.16,0.16}{#1}}
\newcommand{\ErrorTok}[1]{\textcolor[rgb]{0.64,0.00,0.00}{\textbf{#1}}}
\newcommand{\NormalTok}[1]{#1}
\usepackage{graphicx,grffile}
\makeatletter
\def\maxwidth{\ifdim\Gin@nat@width>\linewidth\linewidth\else\Gin@nat@width\fi}
\def\maxheight{\ifdim\Gin@nat@height>\textheight\textheight\else\Gin@nat@height\fi}
\makeatother
% Scale images if necessary, so that they will not overflow the page
% margins by default, and it is still possible to overwrite the defaults
% using explicit options in \includegraphics[width, height, ...]{}
\setkeys{Gin}{width=\maxwidth,height=\maxheight,keepaspectratio}
\IfFileExists{parskip.sty}{%
\usepackage{parskip}
}{% else
\setlength{\parindent}{0pt}
\setlength{\parskip}{6pt plus 2pt minus 1pt}
}
\setlength{\emergencystretch}{3em}  % prevent overfull lines
\providecommand{\tightlist}{%
  \setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}
\setcounter{secnumdepth}{0}
% Redefines (sub)paragraphs to behave more like sections
\ifx\paragraph\undefined\else
\let\oldparagraph\paragraph
\renewcommand{\paragraph}[1]{\oldparagraph{#1}\mbox{}}
\fi
\ifx\subparagraph\undefined\else
\let\oldsubparagraph\subparagraph
\renewcommand{\subparagraph}[1]{\oldsubparagraph{#1}\mbox{}}
\fi

%%% Use protect on footnotes to avoid problems with footnotes in titles
\let\rmarkdownfootnote\footnote%
\def\footnote{\protect\rmarkdownfootnote}

%%% Change title format to be more compact
\usepackage{titling}

% Create subtitle command for use in maketitle
\newcommand{\subtitle}[1]{
  \posttitle{
    \begin{center}\large#1\end{center}
    }
}

\setlength{\droptitle}{-2em}

  \title{TMA4315 Generalized linear models H2018}
    \pretitle{\vspace{\droptitle}\centering\huge}
  \posttitle{\par}
  \subtitle{Module 4: Count and continuous positive response data (Poisson and gamma
regression)}
  \author{Mette Langaas, Department of Mathematical Sciences, NTNU -- with
contributions from Ingeborg Hem}
    \preauthor{\centering\large\emph}
  \postauthor{\par}
      \predate{\centering\large\emph}
  \postdate{\par}
    \date{27.09.2018 and 04.10.2018 {[}PL{]}, 28.09.2018 and 05.10.2018 {[}IL{]}}


\begin{document}
\maketitle

{
\setcounter{tocdepth}{2}
\tableofcontents
}
(Latest changes: 06.10: solutions added. 01.10: small changes for second
week. 27.09: added one Problem for ILw1, moved stuff to w2, added a few
dimensions to score test.)

\section{Overview}\label{overview}

\subsection{Learning material}\label{learning-material}

\begin{itemize}
\tightlist
\item
  Textbook: Fahrmeir et al (2013): Chapter 5.2, 5.3.
\item
  \href{https://www.math.ntnu.no/emner/TMA4315/2018h/TMA4315M4H20180927.pdf}{Classnotes
  27.09.2018}
\item
  \href{https://www.math.ntnu.no/emner/TMA4315/2018h/TMA4315M4H20181004.pdf}{Classnotes
  04.10.2018}
\end{itemize}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Topics}\label{topics}

\subsubsection{\texorpdfstring{\protect\hyperlink{firstweek}{First
week}}{First week}}\label{first-week}

\begin{itemize}
\tightlist
\item
  examples of count data
\item
  the Poisson distribution
\item
  regression with count data
\item
  Poisson regression with log-link
\item
  parameter estimation (ML): log-likelhood, score vector, information
  matrix to give iterative calculations
\item
  asymptotic MLE properties
\item
  confidence intervals and hypothesis tests (Wald, score and LRT)
\end{itemize}

Jump to \protect\hyperlink{interactivew1}{interactive (week 1)}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsubsection{\texorpdfstring{\protect\hyperlink{secondweek}{Second
week}}{Second week}}\label{second-week}

\begin{itemize}
\item
  Count data with Poisson regression (continued)
\item
  deviance, model fit and model choice
\item
  overdispersion
\item
  rate models and offset
\item
  Modelling continuous response data: lognormal and gamma
\item
  the gamma distribution
\item
  the gamma GLM model
\item
  gamma likelihood and derivations thereof
\item
  dispersion parameter: scaled and unscaled deviance
\end{itemize}

Jump to \protect\hyperlink{secondweek}{week 2} and
\protect\hyperlink{interactivew2}{interactive (week 2)}.

\textbf{FIRST WEEK}

\section{Examples of count data}\label{examples-of-count-data}

\begin{itemize}
\tightlist
\item
  the number of automobile thefts pr city worldwide
\item
  the number of UFO sightings around the world
\item
  the number of visits at web pages
\item
  the number of male crabs (satellites) residing nearby a female crab
\item
  the number of goals by the home team and the number of goals for the
  away team in soccer
\item
  the number of newspapers sold at newsagents
\end{itemize}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Sales of newspapers}\label{sales-of-newspapers}

This is a short description of a project run at the
\href{www.nr.no}{Norwegian Computing Centre} a few years ago.

The aim of the project was to provide a statistical model to predict the
number of newspapers to be sold at each of 11 thousand outlets all over
Norway for a given day (``tomorrow'' - maybe scaled to match front page
issues). And then based on the prediction to decide on how many
newspapers to be delivered to each outlet in order to optimize the
overall profit (lost sales if outlet are sold out, return costs if
unsold papers).

Response data: number of newspapers (delivered) sold at each outlet.
Covariate data: type of outlet, but mainly calendar information=
weekday, month, season, public holidays, winter/autumn/easter/xmas,
\ldots{}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Female crabs with
satellites}\label{female-crabs-with-satellites}

The example is taken from Agresti (1996): ``An Introduction to
Categorical Data Analysis'', and the data is from a study of nesting
horseshoe crabs (J. Brockmann, Ethology 1996)

First, the study objects were female crabs (horseshoe crabs). Each
female crab had a male crab attached to her in her nest. The objective
of the study was to investigate factors that affect whether the female
crab had any other males, called satellites, residing near her. The
following covariates were collected for 173 female crabs:

\begin{itemize}
\tightlist
\item
  \texttt{C}: the color of the female crab (1=light medium, 2=medium,
  3=dark medium, 4=dark)
\item
  \texttt{S}: spine condition (1=both good, 2=one worn or broken, 3=both
  worn or broken)
\item
  \texttt{W}: width of carapace (cm)
\item
  \texttt{Wt}: weight (kg)
\end{itemize}

The response was the number of satellites, \texttt{Sa}= male crabs
residing nearby.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(ggplot2)}
\KeywordTok{library}\NormalTok{(GGally)}
\NormalTok{crab =}\StringTok{ }\KeywordTok{read.table}\NormalTok{(}\StringTok{"https://www.math.ntnu.no/emner/TMA4315/2018h/crab.txt"}\NormalTok{)}
\KeywordTok{colnames}\NormalTok{(crab) =}\StringTok{ }\KeywordTok{c}\NormalTok{(}\StringTok{"Obs"}\NormalTok{, }\StringTok{"C"}\NormalTok{, }\StringTok{"S"}\NormalTok{, }\StringTok{"W"}\NormalTok{, }\StringTok{"Wt"}\NormalTok{, }\StringTok{"Sa"}\NormalTok{)}
\NormalTok{crab =}\StringTok{ }\NormalTok{crab[, }\OperatorTok{-}\DecValTok{1}\NormalTok{]  }\CommentTok{#remove column with Obs}
\NormalTok{crab}\OperatorTok{$}\NormalTok{C =}\StringTok{ }\KeywordTok{as.factor}\NormalTok{(crab}\OperatorTok{$}\NormalTok{C)}
\NormalTok{crab}\OperatorTok{$}\NormalTok{S =}\StringTok{ }\KeywordTok{as.factor}\NormalTok{(crab}\OperatorTok{$}\NormalTok{S)}
\KeywordTok{ggpairs}\NormalTok{(crab)}
\end{Highlighting}
\end{Shaded}

\includegraphics{4PoissonGamma_files/figure-latex/unnamed-chunk-1-1.pdf}

\textbf{Q}: Discuss what you see. Any potential covariates to influence
\texttt{Sa}? Which distribution can \texttt{Sa} have?

\begin{verbatim}
## **A:** Seems that `C`, `W` and `Wt` might be correlated with `Sa`. Also `W` and `Wt` are very correlated. Difficult to see the distribution of `Sa` since plotting a Poisson variable with different means might not be easy to recognize.
\end{verbatim}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\section{Modelling counts with the Poisson
distribution}\label{modelling-counts-with-the-poisson-distribution}

\subsection{The Poisson process}\label{the-poisson-process}

We observe events that may occur within a time interval or a region.

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
  The number of events occuring within a time interval or a region, is
  independent of the number of events that occurs in any other disjoint
  (non-overlapping) time interval or region.
\item
  The probability that a single event occurs within a small time
  interval or region, is proportional to the length of the interval or
  the size of the region.
\item
  The probability that more than one event may occur within a small time
  interval or region is negligable.
\end{enumerate}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

When all of these three properties are funfilled we have a \emph{Poisson
process}. This leads to three distributions

\begin{itemize}
\tightlist
\item
  The number of events in a Poisson process follows a Poisson
  distribution.
\item
  Time between two events in a Poisson process follows an exponential
  distribution.
\item
  Time between many events in a Poisson process follows a gamma
  distribution.
\end{itemize}

We will first study the Poisson distribution - and link it to a
regression setting.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{The Poisson distribution}\label{the-poisson-distribution}

We study a Poisson process within a time interval or a region of
specified size. Then, the number of events, \(Y\), will follow a
\emph{Poisson distribution} with parameter \(\lambda\) \[
f(y)=\frac{\lambda^y}{y!}e^{-\lambda} \text{ for } y=0,1,2,…
\] Here the parameter \(\lambda\) is the proportionality factor in the
requirement 2 (above) for the Poisson process. Another popular
parameterization is \(\mu\), or given some interval \(\lambda t\), but
we will stick with \(\lambda\). In R we calculate the Poisson point
probabilities using \texttt{dpois}.

If you want to see how this distribution function is derived form the
binomial distribution you may watch this video:
\href{https://mediasite.ntnu.no/Mediasite/Play/5517061a4ada489eb5eb8d90431923fe1d?catalog=0fce6173-7a98-4db7-84b7-50cba3a3a341}{Poisson
process and distribution}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\textbf{Cumulative distribution (cdf):} The cumulative distribution is
calculated by summing \(F(y)=\sum_{t\le y} f(t)\), and we might
calculate the Poisson cdf in R with \texttt{ppois}.

\textbf{Expected value and variance:} Let \(Y\) follow a Poisson
distribution with parameter \(\lambda\). Then
\[\text{E}(Y)=\lambda \text{ and } \text{Var}(Y)=\lambda\]

\textbf{Proof:} \[
E(y)=\sum_{y=0}^{\infty} y \frac{\lambda^y}{y!}e^{-\lambda}=
\sum_{y=1}^{\infty} y \frac{\lambda^y}{y!}e^{-\lambda}=\sum_{y=1}^{\infty} \frac{\lambda \lambda^{y-1}}{(y-1)!}e^{-\lambda}=\lambda \sum_{z=0}^{\infty} \frac{\lambda^{z}}{z!}e^{-\lambda}=\lambda
\] In the first transition we use that the term \(y=0\) gives no
contribution to the sum. Then we cancel out \(y\) in the numerator with
the first term of \(y!\) in the denominator. Then we let \(z=y-1\), and
finally we use that the sum of the Poisson distribution for all possible
outcomes equal 1.

To calculate the variance we first use that
\[ \text{Var}(Y)=\text{E}(Y^2)-(\text{E}(Y))^2\] and then that
\[Y^2=Y(Y-1)+Y\] so that
\[\text{Var}(Y)=\text{E}(Y(Y-1))+\text{E}(Y)-(\text{E}(Y))^2\] The
reason for this is that we may use the same type of trick as for
\(\text{E}(Y)\) since the sum is over all values of the distribution
function.
\[ \text{E}(Y(Y-1))=\sum_{y=0}^{\infty} y(y-1) \frac{\lambda^y}{y!}e^{-\lambda}=\sum_{y=2}^{\infty} y(y-1) \frac{\lambda^y}{y!}e^{-\lambda}\]
\[=\sum_{y=2}^{\infty} \frac{\lambda^{y-2}\lambda^2}{(y-2)!}e^{-\lambda}=\lambda^2\sum_{y=0}^{\infty}\frac{\lambda^{y}}{y!}e^{-\lambda}=\lambda^2\]

Putting it all togehter
\[\text{Var}(Y)=\lambda^2 +\lambda - \lambda^2=\lambda\] which was to be
shown.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Properties of the Poisson
distribution}\label{properties-of-the-poisson-distribution}

\begin{itemize}
\tightlist
\item
  A sum of \(n\) independent Poisson distributed random variables,
  \(Y_i\) with means \(\lambda_i\) are Poisson distributed with mean
  \(\sum_{i=1}^n \lambda_i\).
\item
  When the mean increases the Poisson distribution becomes more and more
  symmetric and for large \(\lambda\) the Poisson distribution can be
  approximated by a normal distribution.
\end{itemize}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Exponential family}\label{exponential-family}

In Module 1 we introduced distributions of the \(Y_i\), that could be
written in the form of a \emph{univariate exponential family}
\[ f(y_i\mid \theta_i)=\exp \left( \frac{y_i \theta_i-b(\theta_i)}{\phi}\cdot w_i + c(y_i, \phi, w_i) \right) \]
where we said that

\begin{itemize}
\item
  \(\theta_i\) is called the canonical parameter and is a parameter of
  interest
\item
  \(\phi\) is called a nuisance parameter (and is not of interest to
  us=therefore a nuisance (plage))
\item
  \(w_i\) is a weight function, in most cases \(w_i=1\)
\item
  \(b\) and \(c\) are known functions.
\end{itemize}

It can be shown that \(\text{E}(Y_i)=b'(\theta_i)\) and
\(\text{Var}(Y_i)=b''(\theta_i)\cdot \frac{\phi_i}{w_i}\), see
\href{https://www.math.ntnu.no/emner/TMA4315/2017h/M5ExpFamProofEVar.pdf}{derivation
from Module 1}.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

In Module 1 we found that the Poisson distribution
\(Y_i\sim \text{Poisson}(\lambda_i)\) is an exponential family
\href{https://www.math.ntnu.no/emner/TMA4315/2017h/Module1ExponentialFamily.pdf}{derivation
from Module 1},

and that

\begin{itemize}
\tightlist
\item
  \(\theta_i=\ln(\lambda_i)\) is the canonical parameter
\item
  \(\phi=1\), no nuisance
\item
  \(w_i=1\)
\item
  \(b(\theta_i)=\exp(\theta)\)
\item
  \(\mu_i=\text{E}(Y_i)=\lambda_i\)
\end{itemize}

For a GLM with linear predictor \(\eta_i\) - to have a canonical link we
need \[\theta_i=\eta_i\] Since \(\eta_i=g(\mu_i)=g(\lambda_i)\) this
means to us that we need \[ g(\mu_i)=g(\lambda_i)=\theta_i\] saying that
with the Poisson the canonical link is \(\ln(\lambda_i)\).

\textbf{Q}: Why may we want to choose a canonical link?

\section{Regression with count data}\label{regression-with-count-data}

\subsection{Aim}\label{aim}

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  Construct a model to help understand the relationship between a count
  variable and one or many possible explanatory variables. The response
  measurements are counts.
\item
  Use the model for understanding what can explain count, and for
  prediction of counts.
\end{enumerate}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{What is we instead want to use the multiple linear
regression
model?}\label{what-is-we-instead-want-to-use-the-multiple-linear-regression-model}

When modelling the counts some times the normal approximation might be
used, especially when the counts are high. It is also possible to use a
transformation of the response to get a constant variance. The
transformation \(\sqrt{Y_i}\) give an approximate constant variance, but
then it is not clear if there then is a linear relationship between
\(\text{E}(\sqrt{Y_i})\) and the covariates.

In general, with count data, we instead use a Possion GLM regression.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{The linear Poisson model}\label{the-linear-poisson-model}

It is possible to construct a \emph{linear Poisson model} where we have
the direct relationship \[\lambda_i=\eta_i\] between the mean of the
Poisson distribution and the linear predictor.

This means that the covariates have an additive effect on the rate
\(\lambda_i\).

However, since the rate \(\lambda_i\) can not be negative, then to use
this model restrictions on the parameter space of the \(\beta\) is
needed.

We will not use the linear Poisson model, but instead the
\emph{log-linear Poisson model}.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{The log-linear Poisson
model}\label{the-log-linear-poisson-model}

\textbf{Assumptions: }

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  \(Y_i \sim \text{Poisson}(\lambda_i)\), with
  \(\text{E}(Y_i)=\lambda_i\), and \(\text{Var}(Y_i)=\lambda_i\).
\item
  Linear predictor: \(\eta_i={\bf x}_i^T \beta\).
\item
  Log link \[\eta_i=\ln(\lambda_i)=g(\lambda_i)\] and (inverse thereof)
  response function \[\lambda_i=\exp(\eta_i)\]
\end{enumerate}

Assumptions 1 and 3 above can be written as
\[Y_i \sim \text{Poisson}(\exp(\eta_i)), \text{ }i=1,\ldots,n\]

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Interpreting parameters in the log-linear Poisson
model}\label{interpreting-parameters-in-the-log-linear-poisson-model}

In the log-linear model the mean, \(\text{E}(Y_i)=\lambda_i\) satisfy an
exponential relationship to covariates
\[ \lambda_i=\exp(\eta_i)=\exp({\bf x}_i^T \beta)=\exp(\beta_0)\cdot \exp(\beta_1)^{x_{i1}}
\cdots \exp(\beta_k)^{x_{ik}}.\]

Let us look in detail at \(\beta_1\) with covariate \(x_{i1}\) for
observation \(i\).

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
  If \(x_{i1}\) increases by one unit to \(x_{i1}+1\) then the mean
  \(\text{E}(Y_i)\) will in our model change by a factor
  \(\exp(\beta_1)\).
\item
  If \(\beta_1\)=0 then \(\exp(\beta_1)=1\), so that a change in
  \(x_{i1}\) does not change \(\text{E}(Y_i)\).
\item
  If \(\beta_1<0\) then \(\exp(\beta_1)<1\) so if \(x_{i1}\) increase
  then \(\text{E}(Y_i)\) decrease.
\item
  If \(\beta_1>0\) then \(\exp(\beta_1)>1\) so if \(x_{i1}\) increase
  then \(\text{E}(Y_i)\) increase.
\end{enumerate}

Thus, the covariates have a multiplicative effect on the rate
\(\lambda_i\).

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Example: interpreting parameters for the female crabs with
satellites}\label{example-interpreting-parameters-for-the-female-crabs-with-satellites}

We fit a log-linear model to \texttt{Sa}, assuming the number of
satellites follows a Poisson distribution with log-link.

\textbf{Q:}

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
  First the model is fitted with intercept only. What do we assume then?
  Interpret the fit.
\item
  Then width \texttt{W} is added as a covariate in the log-linear model.
  What happens if the width increase by one unit (cm)?
\item
  What is the predicted number of satellites for the average width?
\end{enumerate}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{model1 =}\StringTok{ }\KeywordTok{glm}\NormalTok{(Sa }\OperatorTok{~}\StringTok{ }\DecValTok{1}\NormalTok{, }\DataTypeTok{family =} \KeywordTok{poisson}\NormalTok{(}\DataTypeTok{link =}\NormalTok{ log), }\DataTypeTok{data =}\NormalTok{ crab)}
\KeywordTok{cat}\NormalTok{(}\StringTok{"Intercept only}\CharTok{\textbackslash{}n}\StringTok{"}\NormalTok{)}
\KeywordTok{print}\NormalTok{(model1}\OperatorTok{$}\NormalTok{coefficients)}
\KeywordTok{cat}\NormalTok{(}\StringTok{"Intercept only, exp}\CharTok{\textbackslash{}n}\StringTok{"}\NormalTok{)}
\KeywordTok{exp}\NormalTok{(model1}\OperatorTok{$}\NormalTok{coefficients)}
\KeywordTok{print}\NormalTok{(}\KeywordTok{mean}\NormalTok{(crab}\OperatorTok{$}\NormalTok{Sa))}
\NormalTok{model2 =}\StringTok{ }\KeywordTok{glm}\NormalTok{(Sa }\OperatorTok{~}\StringTok{ }\NormalTok{W, }\DataTypeTok{family =} \KeywordTok{poisson}\NormalTok{(}\DataTypeTok{link =}\NormalTok{ log), }\DataTypeTok{data =}\NormalTok{ crab)}
\KeywordTok{cat}\NormalTok{(}\StringTok{"Intercept + W}\CharTok{\textbackslash{}n}\StringTok{"}\NormalTok{)}
\KeywordTok{print}\NormalTok{(model2}\OperatorTok{$}\NormalTok{coefficients)}
\KeywordTok{cat}\NormalTok{(}\StringTok{"Intercept+W, exp}\CharTok{\textbackslash{}n}\StringTok{"}\NormalTok{)}
\KeywordTok{exp}\NormalTok{(model2}\OperatorTok{$}\NormalTok{coefficients)}
\KeywordTok{cat}\NormalTok{(}\StringTok{"summary of width}\CharTok{\textbackslash{}n}\StringTok{"}\NormalTok{)}
\KeywordTok{summary}\NormalTok{(crab}\OperatorTok{$}\NormalTok{W)}
\KeywordTok{cat}\NormalTok{(}\StringTok{"what is this?}\CharTok{\textbackslash{}n}\StringTok{"}\NormalTok{)}
\KeywordTok{print}\NormalTok{(}\KeywordTok{exp}\NormalTok{(model2}\OperatorTok{$}\NormalTok{coefficients[}\DecValTok{1}\NormalTok{] }\OperatorTok{+}\StringTok{ }\NormalTok{model2}\OperatorTok{$}\NormalTok{coefficients[}\DecValTok{2}\NormalTok{] }\OperatorTok{*}\StringTok{ }\KeywordTok{mean}\NormalTok{(crab}\OperatorTok{$}\NormalTok{W)))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Intercept only
## (Intercept) 
##    1.071267 
## Intercept only, exp
## (Intercept) 
##    2.919075 
## [1] 2.919075
## Intercept + W
## (Intercept)           W 
##  -3.3047572   0.1640451 
## Intercept+W, exp
## (Intercept)           W 
##  0.03670812  1.17826744 
## summary of width
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    21.0    24.9    26.1    26.3    27.7    33.5 
## what is this?
## (Intercept) 
##    2.744061
\end{verbatim}

\section{Parameter estimation with maximum
likelihood}\label{parameter-estimation-with-maximum-likelihood}

Our parameter of interest is the vector \(\beta\) of regression
coefficients, and we have no nuisance parameters. We would like to
estimate \(\beta\) from maximizing the likelihood - the presentation
here is essentially the same as for Module 3: Binary regression - with
``Poisson and log'' instead of ``Bernoulli and logit''. And, also here
we will not have a closed form solution for \(\hat{\beta}\) (except for
a few special cases).

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{\texorpdfstring{Likelihood
\(L(\beta)\)}{Likelihood L(\textbackslash{}beta)}}\label{likelihood-lbeta}

We assume that pairs of covariates and response are measured
independently of each other: \(({\bf x}_i,Y_i)\), and \(Y_i\) follows
the distribution specified above, and \({\bf x}_i\) is fixed.

\[L(\beta)=\prod_{i=1}^n L_i(\beta)=\prod_{i=1}^n f(y_i; \beta)=\prod_{i=1}^n\frac{\lambda_i^{y_i}}{y_i!}\exp(-\lambda_i)\]

\textbf{Note:} still a slight misuse of notation - where is \(\beta\)?

\textbf{A:} \(\eta_i=\ln(\lambda_i)={\bf x}_i^T\beta\), so replace
\(\lambda_i\) with \(\exp({\bf x}_i^T \beta)\).

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{\texorpdfstring{Loglikelihood
\(l(\beta)\)}{Loglikelihood l(\textbackslash{}beta)}}\label{loglikelihood-lbeta}

The log-likelihood is just the natural log of the likelihood, and we
work with the log-likelihood because this makes the mathematics simpler
- since we work with exponential families. The main aim with the
likelihood is to maximize it to find the maximum likelihood estimate,
and since the log is a monotone function the maximum of the
log-likelihood will be in the same place as the maximum of the
likelihood.

\[l(\beta)=\ln L(\beta)=\sum_{i=1}^n \ln L_i(\beta)=\sum_{i=1}^n l_i(\beta)=\sum_{i=1}^n [y_i \ln(\lambda_i)-\lambda_i-\ln(y!)]\]
Observe that the log-likelihood is a sum of invidual contributions for
each observation pair \(i\). We often omit the last term since it is not
a function of model parameters, only data.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

If we want a function of \(\eta_i=\ln(\lambda_i)\) or \(\beta\):
\[l(\beta)=\sum_{i=1}^n[y_i \eta_i-\exp(\eta_i)+C_i]=\sum_{i=1}^n y_i{\bf x}_i^T\beta-\sum_{i=1}^n \exp({\bf x}_i^T\beta) +C\]

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{\texorpdfstring{Score function
\(s(\beta)\)}{Score function s(\textbackslash{}beta)}}\label{score-function-sbeta}

The score function is a \(p\times 1\) vector, \(s(\beta)\), with the
partial derivatives of the log-likelihood with respect to the \(p\)
elements of the \(\beta\) vector. Remember, the score function is linear
in the individual contributions:
\[s(\beta)=\frac{\partial l(\beta)}{\partial \beta}=
\sum_{i=1}^n \frac{\partial l_i(\beta)}{\partial \beta}=
\sum_{i=1}^n s_i(\beta)\]

We work with \(l_i(\beta)=l_i(\pi_i(\eta_i(\beta)))\) and use the chain
rule to find \(s_i(\beta)\).

\[s_i(\beta)=\frac{\partial l_i(\beta)}{\partial \beta}=\frac{\partial l_i(\beta)}{\partial \eta_i}\cdot \frac{\partial \eta_i}{\partial \beta}=\frac{\partial [y_i\eta_i-\exp(\eta_i)+C_i]}{\partial \eta_i}\cdot \frac{\partial [{\bf x}_i^T\beta ]}{\partial \beta}\]
\[=[y_i-\exp(\eta_i)]\cdot {\bf x}_i=(y_i-\lambda_i){\bf x}_i \] See
Module 3 for rules for partial derivatives of scalar wrt vector.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

The score function is given as:

\[s(\beta)=\sum_{i=1}^n s_i(\beta)=\sum_{i=1}^n (y_i-\lambda_i){\bf x}_i\]
Again, observe that
\(\text{E}(s_i(\beta))=\text{E}((Y_i-\lambda_i){\bf x}_i)=0\) because
\(\text{E}(Y_i)=\lambda_i\), and thus also \(\text{E}(s(\beta))=0\).

To find the maximum likelihood estimate \(\hat{\beta}\) we solve the set
of \(p\) non-linear equations: \[s(\hat{\beta})=0\] And, as before we do
that using the Newton-Raphson or Fisher Scoring iterative methods, so we
need the derivative of the score vector (our Fisher information).

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{\texorpdfstring{The expected Fisher information matrix
\(F(\beta)\)}{The expected Fisher information matrix F(\textbackslash{}beta)}}\label{the-expected-fisher-information-matrix-fbeta}

We saw in Module 3 that the expected Fisher information matrix,
\(F(\beta)\) is equal the covariance matrix of the score function.

\begin{align} F(\beta) &= \text{Cov}(s(\beta)) =   \sum_{i=1}^n \text{Cov}(s_i(\beta)) \\
&= \sum_{i=1}^n E\left[\Big(s_i(\beta) - E(s_i(\beta))\Big)\Big(s_i(\beta)-E(s_i(\beta))\Big)^T\right] \\
&= \sum_{i=1}^n E(s_i(\beta)s_i(\beta)^T) = \sum_{i=1}^n F_i(\beta) 
\end{align}

where it is used that the responses \(Y_i\) and \(Y_j\) are independent,
and that \(E(s_i(\beta)) = 0 \ \forall i\).

Remember that \(s_i(\beta)=(Y_i-\lambda_i){\bf x}_i\), then:

\[ F_i(\beta) = E(s_i(\beta)s_i(\beta)^T) = E((Y_i-\lambda_i){\bf x}_i(Y_i-\lambda_i){\bf x}_i^T) = {\bf x}_i{\bf x}_i^T E((Y_i-\lambda_i)^2) = {\bf x}_i{\bf x}_i^T \lambda_i \]
where \(E((Y_i-\lambda_i)^2)=\text{Var}(Y_i)=\lambda\) is the variance
of \(Y_i\). Thus

\[F(\beta) = \sum_{i=1}^n {\bf x}_i{\bf x}_i^T \lambda_i.\]

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{\texorpdfstring{Observed Fisher information matrix
\(H(\beta)\)}{Observed Fisher information matrix H(\textbackslash{}beta)}}\label{observed-fisher-information-matrix-hbeta}

We do not really need the observed version of the Fisher information
matrix, and since we use canonical link \(H(\beta)=F(\beta)\) - so we
already have it.

But, for completeness, we add the direct derivation of \(H(\beta)\).

\[H(\beta) = -\frac{\partial^2l(\beta)}{\partial\beta\partial\beta^T} = -\frac{\partial s(\beta)}{\partial\beta^T} = \frac{\partial}{\partial\beta^T}\left[\sum_{i=1}^n (\lambda_i-y_i){\bf x}_i \right] \]
\[=\frac{\partial}{\partial\beta^T}\left[\sum_{i=1}^n (\exp(\eta_i)-y_i){\bf x}_i \right]\]

because \(s(\beta) = \sum_{i=1}^n (y_i-\lambda_i){\bf x}_i\) and hence
\(-s(\beta) = \sum_{i=1}^n (\lambda_i-y_i){\bf x}_i\). Note that
\(\lambda_i = \exp(\eta_i)\).

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\[H(\beta) = \sum_{i=1}^n \frac{\partial}{\partial\beta^T}[{\bf x}_i\lambda_i-{\bf x}_iy_i] = \sum_{i=1}^n \frac{\partial}{\partial\beta^T}{\bf x}_i\lambda_i = \sum_{i=1}^n {\bf x}_i \frac{\partial \lambda_i}{\partial \eta_i} \frac{\partial \eta_i}{\partial \beta^T} \]

Use that

\[ \frac{\partial \eta_i}{\partial \beta^T}=\frac{\partial {\bf x}_i^T\beta}{\partial \beta^T} = \left(\frac{\partial {\bf x}_i^T\beta}{\partial \beta}\right)^T = {\bf x}_i^T \]

and

\[ \frac{\partial \lambda_i}{\partial \eta_i} =  \frac{\partial\exp(\eta_i)}{\partial \eta_i} = \exp(\eta_i)=\lambda_i\]

And thus

\[H(\beta) =\sum_{i=1}^n {\bf x}_i {\bf x}_i^T\lambda_i.\] Note that the
observed and the expected Fisher information matrix are equal (see below
- canonical link - that this is a general finding).

\section{Parameter estimation - in
practice}\label{parameter-estimation---in-practice}

To find the ML estimate \(\hat{\beta}\) we need to solve
\[s(\hat{\beta})=0\] We have that the score function for the log-linear
model is:
\[s(\beta)=\sum_{i=1}^n {\bf x}_i (y_i-\lambda_i)=\sum_{i=1}^n {\bf x}_i (y_i-\exp({\bf x}_i^T\beta)).\]
Observe that this is a non-linear function in \(\beta\), and has no
closed form solution (except for a few special cases).

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Fisher scoring}\label{fisher-scoring}

To solve this we use the Fisher scoring algorithm, were we at interation
\(t+1\) have

\[\beta^{(t+1)}=\beta^{(t)} + F(\beta^{(t)})^{-1} s(\beta^{(t)})\]

Remark: what do we need to do to use the Newton-Raphson method instead?
Well, replace \(F\) with \(H\), but for canonical link (which is the
log-link for the Poisson) \(F=H\).

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Requirements for
convergence}\label{requirements-for-convergence}

For the Fisher scoring algorithm the expected Fisher information matrix
\(F\) needs to be invertible, and analogous to what we saw in Module 3
this is possible if \(\lambda_i>0\) for all \(i\) and that the design
matrix has full rank (\(p\)). Since we have the log-link we have that
\(\lambda_i=\exp({\bf x}_i^T\beta)\) which is always positive, so all
good. Note, with the linear link \(\lambda_i=\eta_i\) this might be a
challenge, and restrictions on \(\beta\) must be set.

Again, it is possible that the algorithm does not converge. This may
happen for ``unfavorable'' data configurations (especially for small
samples). Accoring to our text book, Fahrmeir et al (2013), page 284,
the conditions for uniqueness and existence of ML estimators are very
complex, and the authors suggest that the GLM user instead checks for
convergence in practice by performing the iterations - also for the
Poisson log-linear model.

\section{Statistical inference}\label{statistical-inference}

\subsection{Asymptotic properties of ML
estimates}\label{asymptotic-properties-of-ml-estimates}

We repeat what we found for Module 3: Under some (weak) regularity
conditions:

Let \(\hat{\beta}\) be the maximum likelihood (ML) estimate in the GLM
model. As the total sample size increases, \(n\rightarrow \infty\):

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
  \(\hat{\beta}\) exists
\item
  \(\hat{\beta}\) is consistent (convergence in probability, yielding
  asymptotically unbiased estimator, variances goes towards 0)
\item
  \(\hat{\beta} \approx N_p(\beta,F^{-1}(\hat{\beta}))\)
\end{enumerate}

Observe that this means that asymptotically
\(\text{Cov}(\hat{\beta})=F^{-1}(\hat{\beta})\): the inverse of the
expected Fisher information matrix evaluated at the ML estimate.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

In our case we have
\[F(\beta)=\sum_{i=1}^n {\bf x}_i{\bf x}_i^T \lambda_i={\bf X}^T {\bf W} {\bf X},\]
where \({\bf W}=\text{diag}(\lambda_i)\). This means
\(\text{Cov}(\hat{\beta})=({\bf X}^T {\bf W} {\bf X})^{-1}\) (remember
that \(\hat{\beta}\) comes in with \(\hat{\lambda}_i\) in \({\bf W}\)).

Let \({\bf A}(\beta)=F^{-1}(\beta)\), and \(a_{jj}(\beta)\) is diagonal
element number \(j\).

For one element of the parameter vector:
\[ Z_j=\frac{\hat{\beta}_j-\beta_j}{\sqrt{a_{jj}(\hat{\beta})}}\] is
standard normal, which can be used to make confidence intervals - and
test hypotheses.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Confidence intervals}\label{confidence-intervals}

In addition to providing a parameter estimate for each element of our
parameter vector \(\beta\) we should also report a \((1-\alpha)100\)\%
confidence interval (CI) for each element.

Let \(z_{\alpha/2}\) be such that \(P(Z_j>z_{\alpha/2})=\alpha/2\). We
then use \[ P(-z_{\alpha/2}\le Z_j \le z_{\alpha/2})=1-\alpha\] insert
\(Z_j\) and solve for \(\beta_j\) to get
\[ P(\hat{\beta}_j-z_{\alpha/2}\sqrt{a_{jj}(\hat{\beta})}
\le \beta_j \le \hat{\beta}_j-z_{\alpha/2}\sqrt{a_{jj}(\hat{\beta})})=1-\alpha\]

A \((1-\alpha)\)\% CI for \(\beta_j\) is when we insert numerical values
for the upper and lower limits.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Example: Female crabs with
satellites}\label{example-female-crabs-with-satellites}

\textbf{Q:}

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
  Explain what is done in the \texttt{R}-print-out.
\item
  What if we instead want a CI for \(\beta_0+\beta_1 x_{i1}\)?
\item
  What if we instead want a CI for
  \(\lambda_i=\exp(\beta_0+\beta_1 x_{i1})\)?
\end{enumerate}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{model2 =}\StringTok{ }\KeywordTok{glm}\NormalTok{(Sa }\OperatorTok{~}\StringTok{ }\NormalTok{W, }\DataTypeTok{family =} \KeywordTok{poisson}\NormalTok{(}\DataTypeTok{link =}\NormalTok{ log), }\DataTypeTok{data =}\NormalTok{ crab)}
\KeywordTok{summary}\NormalTok{(model2)}
\KeywordTok{cat}\NormalTok{(}\StringTok{"lower}\CharTok{\textbackslash{}n}\StringTok{"}\NormalTok{)}
\NormalTok{lower =}\StringTok{ }\NormalTok{model2}\OperatorTok{$}\NormalTok{coefficients }\OperatorTok{-}\StringTok{ }\KeywordTok{qnorm}\NormalTok{(}\FloatTok{0.975}\NormalTok{) }\OperatorTok{*}\StringTok{ }\KeywordTok{sqrt}\NormalTok{(}\KeywordTok{diag}\NormalTok{(}\KeywordTok{vcov}\NormalTok{(model2)))}
\NormalTok{lower}
\KeywordTok{cat}\NormalTok{(}\StringTok{"upper}\CharTok{\textbackslash{}n}\StringTok{"}\NormalTok{)}
\NormalTok{upper =}\StringTok{ }\NormalTok{model2}\OperatorTok{$}\NormalTok{coefficients }\OperatorTok{+}\StringTok{ }\KeywordTok{qnorm}\NormalTok{(}\FloatTok{0.975}\NormalTok{) }\OperatorTok{*}\StringTok{ }\KeywordTok{sqrt}\NormalTok{(}\KeywordTok{diag}\NormalTok{(}\KeywordTok{vcov}\NormalTok{(model2)))}
\NormalTok{upper}
\KeywordTok{confint}\NormalTok{(model2)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## glm(formula = Sa ~ W, family = poisson(link = log), data = crab)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.8526  -1.9884  -0.4933   1.0970   4.9221  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.30476    0.54224  -6.095  1.1e-09 ***
## W            0.16405    0.01997   8.216  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 632.79  on 172  degrees of freedom
## Residual deviance: 567.88  on 171  degrees of freedom
## AIC: 927.18
## 
## Number of Fisher Scoring iterations: 6
## 
## lower
## (Intercept)           W 
##  -4.3675312   0.1249137 
## upper
## (Intercept)           W 
##  -2.2419833   0.2031764 
##                  2.5 %     97.5 %
## (Intercept) -4.3662326 -2.2406858
## W            0.1247244  0.2029871
\end{verbatim}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Hypothesis testing}\label{hypothesis-testing}

There are three methods that are mainly used for testing hypotheses in
GLMs - these are called Wald test, likelihood ratio test and score test.
In Module 3 we looked at the Wald and likelihood ratio test - and what
we found there still applies for the Poisson GLM. We only repeat the
previous findings, and then add an example with the crab data.

We will look at null hypotheses and alternative hypotheses that can be
written
\[H_0: {\bf C}{\bf \beta}={\bf d} \text{ vs. } {\bf C}{\bf \beta}\neq {\bf d} \]
by specifying \({\bf C}\) to be a \(r \times p\) matrix and \({\bf d}\)
to be a column vector of length \(r\),

and/or where we define

\begin{itemize}
\tightlist
\item
  A: the larger model and
\item
  B: the smaller model (under \(H_0\)), and the smaller model is nested
  within the larger model (that is, B is a submodel of A).
\end{itemize}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsubsection{The Wald test}\label{the-wald-test}

The Wald test statistic is given as:
\[w=({\bf C}\hat{{\bf \beta}}-{\bf d})^{\text T}[{\bf C}F^{-1}(\hat{\beta}){\bf C}^{\text T}]^{-1}({\bf C}\hat{{\bf \beta}}-{\bf d}) \]
and measures the distance between the estimate \({\bf C}\hat{\beta}\)
and the value under then null hypothesis \({\bf d}\), weighted by the
asymptotic covariance matrix of \({\bf C}\hat{\beta}\), and under the
null follows a \(\chi^2\) distribution with \(r\) degrees of freedom
(where \(r\) is the number of hypotheses tested).

\(P\)-values are calculated in the upper tail of the
\(\chi^2\)-distribution.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsubsection{The likelihood ratio
test}\label{the-likelihood-ratio-test}

Notation:

\begin{itemize}
\tightlist
\item
  A: the larger model and
\item
  B: the smaller model (under \(H_0\)), and the smaller model is nested
  within the larger model (that is, B is a submodel of A).
\end{itemize}

The likelihood ratio statistic is defined as
\[- 2\ln \lambda=-2(\ln L(\hat{\beta}_B)-\ln L(\hat{\beta}_A)) \] which
under the null is asymptotically \(\chi^2\)-distributed with degrees of
freedom equal the difference in the number of parameters in the large
and the small model.

Again, \(p\)-values are calculated in the upper tail of the
\(\chi^2\)-distribution.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsubsection{The score test}\label{the-score-test}

We continue to use the notation that A: the larger model (under \(H_1\))
and B: the smaller model (under \(H_0\)), and still assume that the
smaller model is nested within the larger model (that is, B is a
submodel of A).

The \emph{score statistics} is based on the \emph{score function}, and
measures the distance to the score function at the maximum likelihood
for model A (which is 0) and scales with the covariance to form the test
statistic.

\begin{itemize}
\tightlist
\item
  Under the null hypothesis investigated let
  \(\tilde{\boldsymbol{\beta}}\) be the ML estimate (that is, model B,
  the smaller model) - that means that this is a restricted ML estimate,
  and
\item
  under \(H_1\) we have the larger model (A) with maximum likelihood
  \(\hat{\boldsymbol{\beta}}\).
\end{itemize}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

The score statistics is:

\[ U=(s(\tilde{\boldsymbol{\beta}})-{\bf 0})^T{\bf F}^{-1}(\tilde{\boldsymbol{\beta}})(s(\tilde{\boldsymbol{\beta}})-{\bf 0})\]

Here \(s(\tilde{\boldsymbol{\beta}})\) represents a subvector of the
score function where only the elements that is in \(H_1\) and not in
\(H_0\) are present (so dimension is the difference in number of
parameters between A and B models), but the score function is evaluated
based on parameter estimates under \(H_0\) (i.e.~the value of
\(\hat{\lambda}\) in the score and expected Fisher information is based
on \(\tilde{\boldsymbol{\beta}}\).

To calculate \({\bf F}^{-1}(\tilde{\boldsymbol{\beta}})\) this is a
submatrix of the full inverted matrix (not invert just the submatrix).
The dimension of this matrix is ``difference in number of parameters
between A and B models'' squared.

When the null hypothesis is true \(U\) as an asymptotic
\(\chi^2\)-distribution with \(r\) degrees of freedom (difference in
number of estimated parameters between the large and small model).

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

Remark: In R the function \texttt{glm.scoretest} in the package
\texttt{statmod} calculates the score test for a GLM when the difference
between \(H_0\) and \(H_1\) is one parameter. The output from the
\texttt{glm.scoretest} is the square root of \(U\) and must be squared
and related to a \(\chi^2_1\) distribution (see example below).

The score test is very useful for special situations when the smaller
model is to be tested towards many larger models, because only the
smaller model has to be fitted.

The score test is perhaps the the most complex and least studied of the
three tests, and in this course the main focus will be on the Wald and
LRT tests. But, it is important for you to have heard of the score test,
because in special situation it may be the preferred test.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsubsection{Example: Female crabs with satellites - the different
tests.}\label{example-female-crabs-with-satellites---the-different-tests.}

We fit a model with two covariates, one is categorical and we use effect
coding. We want to test if we need to add these covariates.

\textbf{Q}: Comment on what you see in the print-out. Which hypotheses
are we testing? Compare the test statistics.

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(statmod)}
\NormalTok{model3 =}\StringTok{ }\KeywordTok{glm}\NormalTok{(Sa }\OperatorTok{~}\StringTok{ }\NormalTok{W }\OperatorTok{+}\StringTok{ }\NormalTok{C, }\DataTypeTok{family =} \KeywordTok{poisson}\NormalTok{(}\DataTypeTok{link =}\NormalTok{ log), }\DataTypeTok{data =}\NormalTok{ crab, }\DataTypeTok{contrasts =} \KeywordTok{list}\NormalTok{(}\DataTypeTok{C =} \StringTok{"contr.sum"}\NormalTok{))}
\KeywordTok{summary}\NormalTok{(model3)}

\CommentTok{# possible to use type III in Anova}
\KeywordTok{library}\NormalTok{(car)}
\KeywordTok{Anova}\NormalTok{(model3, }\DataTypeTok{type =} \StringTok{"III"}\NormalTok{, }\DataTypeTok{test.statistic =} \StringTok{"Wald"}\NormalTok{)  }\CommentTok{#Q:same as summary?}
\KeywordTok{Anova}\NormalTok{(model3, }\DataTypeTok{type =} \StringTok{"III"}\NormalTok{, }\DataTypeTok{test.statistic =} \StringTok{"LR"}\NormalTok{)  }\CommentTok{#same as comparing deviances}

\CommentTok{# compare to type I with anova - just to remember the difference between}
\CommentTok{# type I and III anova(model3,test='Chisq') anova(model3,test='LRT')}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## glm(formula = Sa ~ W + C, family = poisson(link = log), data = crab, 
##     contrasts = list(C = "contr.sum"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.0415  -1.9581  -0.5575   0.9830   4.7523  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.92089    0.56010  -5.215 1.84e-07 ***
## W            0.14934    0.02084   7.166 7.73e-13 ***
## C1           0.27085    0.11784   2.298   0.0215 *  
## C2           0.07117    0.07296   0.975   0.3294    
## C3          -0.16551    0.09316  -1.777   0.0756 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 632.79  on 172  degrees of freedom
## Residual deviance: 559.34  on 168  degrees of freedom
## AIC: 924.64
## 
## Number of Fisher Scoring iterations: 6
## 
## Analysis of Deviance Table (Type III tests)
## 
## Response: Sa
##             Df  Chisq Pr(>Chisq)    
## (Intercept)  1 27.196  1.839e-07 ***
## W            1 51.350  7.726e-13 ***
## C            3  8.518    0.03644 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Deviance Table (Type III tests)
## 
## Response: Sa
##   LR Chisq Df Pr(>Chisq)    
## W   49.794  1  1.707e-12 ***
## C    8.534  3    0.03618 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\end{verbatim}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\textbf{Optional:} We also include the results for the score test, but
to test an hypothesis for a factor (more than one column in our design
matrix) the \texttt{glm.scoretest} can not be used and then some manual
programming is done below. There exists other packages to perform score
tests for GLMs, but is not considered here.

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{# for score test we need ML under H0, and we need to fit two different H0}
\CommentTok{# models here when H0 is model with intercept and W:}
\NormalTok{model2a =}\StringTok{ }\KeywordTok{glm}\NormalTok{(Sa }\OperatorTok{~}\StringTok{ }\NormalTok{W, }\DataTypeTok{family =} \KeywordTok{poisson}\NormalTok{(}\DataTypeTok{link =}\NormalTok{ log), }\DataTypeTok{data =}\NormalTok{ crab)}
\CommentTok{# when H0 is model with intercept and C:}
\NormalTok{model2b =}\StringTok{ }\KeywordTok{glm}\NormalTok{(Sa }\OperatorTok{~}\StringTok{ }\NormalTok{C, }\DataTypeTok{family =} \KeywordTok{poisson}\NormalTok{(}\DataTypeTok{link =}\NormalTok{ log), }\DataTypeTok{data =}\NormalTok{ crab, }\DataTypeTok{contrasts =} \KeywordTok{list}\NormalTok{(}\DataTypeTok{C =} \StringTok{"contr.sum"}\NormalTok{))}

\CommentTok{# First H0 is intercept and C and H1 is to add W (only one parameter)}
\NormalTok{testobs3b =}\StringTok{ }\KeywordTok{glm.scoretest}\NormalTok{(model2b, }\DataTypeTok{x2 =}\NormalTok{ crab}\OperatorTok{$}\NormalTok{W)}\OperatorTok{^}\DecValTok{2}
\NormalTok{testobs3b}
\DecValTok{1} \OperatorTok{-}\StringTok{ }\KeywordTok{pchisq}\NormalTok{(testobs3b, }\DecValTok{1}\NormalTok{)}

\CommentTok{# alternatively, implement from scratch}
\NormalTok{X =}\StringTok{ }\KeywordTok{model.matrix}\NormalTok{(model3)}
\KeywordTok{colnames}\NormalTok{(X)}
\NormalTok{estlambda =}\StringTok{ }\NormalTok{model2b}\OperatorTok{$}\NormalTok{fitted.values  }\CommentTok{# the means under H0 (equals the variances)}
\NormalTok{W =}\StringTok{ }\KeywordTok{diag}\NormalTok{(estlambda)  }\CommentTok{# variance on diagonal, estimated under H0}
\NormalTok{XH0 =}\StringTok{ }\NormalTok{X[, }\OperatorTok{-}\DecValTok{2}\NormalTok{]  }\CommentTok{# model under H0}
\NormalTok{XH1 =}\StringTok{ }\NormalTok{X[, }\DecValTok{2}\NormalTok{]  }\CommentTok{# part of model under H1 but not under H0}
\NormalTok{scorefunH1 =}\StringTok{ }\KeywordTok{t}\NormalTok{(XH1) }\OperatorTok{%*%}\StringTok{ }\NormalTok{(crab}\OperatorTok{$}\NormalTok{Sa }\OperatorTok{-}\StringTok{ }\NormalTok{estlambda)}
\NormalTok{expFH1full =}\StringTok{ }\KeywordTok{t}\NormalTok{(X) }\OperatorTok{%*%}\StringTok{ }\NormalTok{W }\OperatorTok{%*%}\StringTok{ }\NormalTok{X}
\NormalTok{fullinv =}\StringTok{ }\KeywordTok{solve}\NormalTok{(expFH1full)}
\NormalTok{subfullinv =}\StringTok{ }\NormalTok{fullinv[}\DecValTok{2}\NormalTok{, }\DecValTok{2}\NormalTok{]}
\NormalTok{scoreteststat =}\StringTok{ }\KeywordTok{t}\NormalTok{(scorefunH1) }\OperatorTok{%*%}\StringTok{ }\NormalTok{subfullinv }\OperatorTok{%*%}\StringTok{ }\NormalTok{scorefunH1}
\NormalTok{scoreteststat}
\DecValTok{1} \OperatorTok{-}\StringTok{ }\KeywordTok{pchisq}\NormalTok{(scoreteststat, }\DecValTok{1}\NormalTok{)}
\CommentTok{# difference within numerical accuracy}

\CommentTok{# test if C should be included in a model with W}
\NormalTok{mm =}\StringTok{ }\KeywordTok{model.matrix}\NormalTok{(Sa }\OperatorTok{~}\StringTok{ }\NormalTok{C, }\DataTypeTok{data =}\NormalTok{ crab, , }\DataTypeTok{contrasts =} \KeywordTok{list}\NormalTok{(}\DataTypeTok{C =} \StringTok{"contr.sum"}\NormalTok{))}
\NormalTok{testobs3a =}\StringTok{ }\KeywordTok{glm.scoretest}\NormalTok{(model2a, }\DataTypeTok{x2 =}\NormalTok{ mm[, }\OperatorTok{-}\DecValTok{1}\NormalTok{])}\OperatorTok{^}\DecValTok{2}  \CommentTok{#only test one level at a time, not all together - can not be used in general}
\NormalTok{testobs3a}
\DecValTok{1} \OperatorTok{-}\StringTok{ }\KeywordTok{pchisq}\NormalTok{(testobs3a, }\DecValTok{1}\NormalTok{)}
\CommentTok{# need to use other implementation (like above)}
\NormalTok{X =}\StringTok{ }\KeywordTok{model.matrix}\NormalTok{(model3)}
\KeywordTok{colnames}\NormalTok{(X)}
\NormalTok{estlambda =}\StringTok{ }\NormalTok{model2a}\OperatorTok{$}\NormalTok{fitted.values  }\CommentTok{# the means under H0 (and variances since Poisson)}
\NormalTok{W =}\StringTok{ }\KeywordTok{diag}\NormalTok{(estlambda)  }\CommentTok{# variance on diagonal, estimated under H0}
\NormalTok{XH0 =}\StringTok{ }\NormalTok{X[, }\DecValTok{1}\OperatorTok{:}\DecValTok{2}\NormalTok{]  }\CommentTok{# model under H0}
\NormalTok{XH1 =}\StringTok{ }\NormalTok{X[, }\DecValTok{3}\OperatorTok{:}\DecValTok{5}\NormalTok{]  }\CommentTok{# part of model under H1 but not under H0}
\NormalTok{scorefunH1 =}\StringTok{ }\KeywordTok{t}\NormalTok{(XH1) }\OperatorTok{%*%}\StringTok{ }\NormalTok{(crab}\OperatorTok{$}\NormalTok{Sa }\OperatorTok{-}\StringTok{ }\NormalTok{estlambda)}
\NormalTok{expFH1full =}\StringTok{ }\KeywordTok{t}\NormalTok{(X) }\OperatorTok{%*%}\StringTok{ }\NormalTok{W }\OperatorTok{%*%}\StringTok{ }\NormalTok{X}
\NormalTok{fullinv =}\StringTok{ }\KeywordTok{solve}\NormalTok{(expFH1full)}
\NormalTok{subfullinv =}\StringTok{ }\NormalTok{fullinv[}\DecValTok{3}\OperatorTok{:}\DecValTok{5}\NormalTok{, }\DecValTok{3}\OperatorTok{:}\DecValTok{5}\NormalTok{]}
\NormalTok{scoreteststat =}\StringTok{ }\KeywordTok{t}\NormalTok{(scorefunH1) }\OperatorTok{%*%}\StringTok{ }\NormalTok{subfullinv }\OperatorTok{%*%}\StringTok{ }\NormalTok{scorefunH1}
\NormalTok{scoreteststat}
\DecValTok{1} \OperatorTok{-}\StringTok{ }\KeywordTok{pchisq}\NormalTok{(scoreteststat, }\DecValTok{3}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] 51.51423
## [1] 7.107648e-13
## [1] "(Intercept)" "W"           "C1"          "C2"          "C3"         
##          [,1]
## [1,] 51.51418
##              [,1]
## [1,] 7.107648e-13
##        C1        C2        C3 
## 4.2333184 2.6565569 0.5677206 
##         C1         C2         C3 
## 0.03963787 0.10312375 0.45116612 
## [1] "(Intercept)" "W"           "C1"          "C2"          "C3"         
##          [,1]
## [1,] 8.578756
##            [,1]
## [1,] 0.03544893
\end{verbatim}

\section{Interactive session - first
week}\label{interactive-session---first-week}

\subsection{Problem 1: Exam 2005 (Problem 1d-f - slightly modified) -
Female crabs and
satellites}\label{problem-1-exam-2005-problem-1d-f---slightly-modified---female-crabs-and-satellites}

(Only a subset of 20 crabs in the original data was used in the exam
problems, but we will use the full data set of 173 crabs - so results
will not be the same. Permitted aids for the exam was ``all printed and
handwritten material and all calculators'' - NB that is not the case for
2018!)

We assume that the number of satellites for each female crab follows a
Poisson distribution and want to perform a Poisson regression using a
log-link to find if there is a connection between the expected number of
satellites \texttt{Sa} and the width \texttt{W} and colour \texttt{C} of
the carapace of the female crab.

\begin{itemize}
\tightlist
\item
  \texttt{C}: the color of the female crab (1=light medium, 2=medium,
  3=dark medium, 4=dark)
\item
  \texttt{W}: width of carapace (cm)
\end{itemize}

The following model was fitted.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{crab <-}\StringTok{ }\KeywordTok{read.table}\NormalTok{(}\StringTok{"https://www.math.ntnu.no/emner/TMA4315/2018h/crab.txt"}\NormalTok{)}
\KeywordTok{colnames}\NormalTok{(crab) <-}\StringTok{ }\KeywordTok{c}\NormalTok{(}\StringTok{"Obs"}\NormalTok{, }\StringTok{"C"}\NormalTok{, }\StringTok{"S"}\NormalTok{, }\StringTok{"W"}\NormalTok{, }\StringTok{"Wt"}\NormalTok{, }\StringTok{"Sa"}\NormalTok{)}
\NormalTok{crab <-}\StringTok{ }\NormalTok{crab[,}\OperatorTok{-}\DecValTok{1}\NormalTok{] }\CommentTok{#remove column with Obs}
\NormalTok{crab}\OperatorTok{$}\NormalTok{C <-}\StringTok{ }\KeywordTok{as.factor}\NormalTok{(crab}\OperatorTok{$}\NormalTok{C)}
\NormalTok{modelEXAM <-}\StringTok{ }\KeywordTok{glm}\NormalTok{(Sa }\OperatorTok{~}\StringTok{ }\NormalTok{W }\OperatorTok{+}\StringTok{ }\NormalTok{C, }\DataTypeTok{family =} \KeywordTok{poisson}\NormalTok{(}\DataTypeTok{link =}\NormalTok{ log), }\DataTypeTok{data =}\NormalTok{ crab, }\DataTypeTok{contrasts =} \KeywordTok{list}\NormalTok{(}\DataTypeTok{C =} \StringTok{"contr.sum"}\NormalTok{))}
\KeywordTok{summary}\NormalTok{(modelEXAM)}
\KeywordTok{library}\NormalTok{(car)}
\KeywordTok{Anova}\NormalTok{(modelEXAM, }\DataTypeTok{type =} \StringTok{"III"}\NormalTok{, }\DataTypeTok{test.statistic =} \StringTok{"Wald"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## glm(formula = Sa ~ W + C, family = poisson(link = log), data = crab, 
##     contrasts = list(C = "contr.sum"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.0415  -1.9581  -0.5575   0.9830   4.7523  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.92089    0.56010  -5.215 1.84e-07 ***
## W            0.14934    0.02084   7.166 7.73e-13 ***
## C1           0.27085    0.11784   2.298   0.0215 *  
## C2           0.07117    0.07296   0.975   0.3294    
## C3          -0.16551    0.09316  -1.777   0.0756 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 632.79  on 172  degrees of freedom
## Residual deviance: 559.34  on 168  degrees of freedom
## AIC: 924.64
## 
## Number of Fisher Scoring iterations: 6
## 
## Analysis of Deviance Table (Type III tests)
## 
## Response: Sa
##             Df  Chisq Pr(>Chisq)    
## (Intercept)  1 27.196  1.839e-07 ***
## W            1 51.350  7.726e-13 ***
## C            3  8.518    0.03644 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\end{verbatim}

\textbf{a)} Write down the estimated model and evaluate the significance
of the estimated coefficients.

How would you proceed to evaluate the significance of the coefficient
for the colour C4. Explain.

Perform a test to evaluate the fit of the model. (In Problem 3d we
derive the mathematical formula for the test statistic this test is
based on.)

Use a 5\% significance level for your evaluations. What is your
conclusion?

\textbf{b)} Make a sketch \emph{by hand} illustrating the connection
between the expected number of satellites and the width of the carapace
for each of the four colours of the crab.

The estimated multiplicative change in the expected number of satellites
when the width of the carapace is changed by 1 cm is, according to the
model, independent of the colour. Explain why.

Also find a 95\% confidence interval for this change.

\textbf{c)} Let \(\hat{\eta}(x_{W})\) be the estimated linear predictor
when the width of the carapace is \(x_{W}\) and the crab is ``light
medium''. What is the distribution of \(\hat{\eta}(x_{W})\)? Which value
of \(x_W\) would give estimated mean number of satellites equal to 5?

Optional: Using R -- for this value of \(x_W\), construct a 95\%
confidence interval for the mean number of satellites.

Optional: Use \texttt{ggplot} to create (and improve) the sketch from
b).

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Problem 2: Exam 2017 (Problem 1) - Poisson
regression}\label{problem-2-exam-2017-problem-1---poisson-regression}

Consider a random variable \(Y\). In our course we have considered the
univariate exponential family having distribution (probability density
function for continuous variables and probability mass function for
discrete variables)
\[ f(y)=\exp \left( \frac{y \theta-b(\theta)}{\phi}w + c(y,\phi,w) \right)\]
where \(\theta\) is called the \emph{natural parameter} (or parameter of
interest) and \(\phi\) the \emph{dispersion parameter}.

The Poisson distribution is a discrete distribution with probability
mass function
\[ f(y)=\frac{\lambda^{y}}{y!}\exp(- \lambda), \text{ for } y=0,1,\ldots\]
where \(\lambda>0\).

\textbf{a)} (10 points)

\begin{itemize}
\tightlist
\item
  Show that the Poisson distribution is a univariate exponential family,
  and specify what the elements of the exponential family (\(\theta\),
  \(\phi\), \(b(\theta)\), \(w\), \(c(y,\phi,w)\)) are.
\item
  What is the connections between \(\text{E}(Y)\) and elements of the
  exponential family?
\item
  What is the connections between \(\text{Var}(Y)\) and elements of the
  exponential family?
\item
  Use these connections to derive the mean and variance for the Poisson
  distribution.
\item
  If the Poisson distribution is used as the distribution for the
  response in a generalized linear model, what is the \emph{canonical
  link} function?
\end{itemize}

\textbf{b)} (15 points)

We consider a Poisson regression with log link
\(\eta_i=g(\mu_i)=\ln(\mu_i)\), and linear predictor equal to
\(\eta_i={\bf x}_i^T \boldsymbol{\beta}\). Further, let \(p\) be the
number of regression parameters in \(\boldsymbol{\beta}\) (intercept
included). The response--covariate pairs \((Y_i, {\bf x}_i)\) are
independent for \(i=1,\ldots, n\).

\begin{itemize}
\tightlist
\item
  Does this set-up satisfy the \emph{requirements} of a GLM model?
  Explain.
\item
  Write down the log-likelihood.
\item
  From the log-likelihood, \emph{derive} the formula for the score
  function \({\bf s}(\boldsymbol{\beta})\) and the expected Fisher
  information matrix, \({\bf F}(\boldsymbol{\beta})\).
\item
  What are the dimensions of \({\bf s}(\boldsymbol{\beta})\) and
  \({\bf F}(\boldsymbol{\beta})\)?
\item
  How can \({\bf s}(\boldsymbol{\beta})\) and
  \({\bf F}(\boldsymbol{\beta})\) be used to arrive at a maximum
  likelihood estimate for \(\boldsymbol{\beta}\)?
\end{itemize}

\textbf{c)} (10 points)

We now look at a data set giving the number of species of tortoise on
the various Galapagos Islands (Data taken from the book ``Practical
Regression and Anova using R'' by Julian J. Faraway.).

The data set contains measurements on 30 islands, and we study the
following variables:

\begin{itemize}
\tightlist
\item
  \texttt{Species}: The number of species of tortoise found on the
  island.
\item
  \texttt{Area}: The area of the island (km\(^2\)).
\item
  \texttt{Elevation}: The highest elevation of the island (m).
\item
  \texttt{Nearest}: The distance from the nearest island (km).
\item
  \texttt{Scruz}: The distance from Santa Cruz island (km).
\item
  \texttt{Adjacent}: The area of the adjacent island (km\(^2\)).
\end{itemize}

We have fitted a Poisson regression with log link to \texttt{Species} as
response, and the other five variables are used as continuous
covariates. Print-out from the fitted model is given in below.

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(faraway)}
\NormalTok{fit <-}\StringTok{ }\KeywordTok{glm}\NormalTok{(Species }\OperatorTok{~}\StringTok{ }\NormalTok{Area }\OperatorTok{+}\StringTok{ }\NormalTok{Elevation }\OperatorTok{+}\StringTok{ }\NormalTok{Nearest }\OperatorTok{+}\StringTok{ }\NormalTok{Scruz }\OperatorTok{+}\StringTok{ }\NormalTok{Adjacent, }\DataTypeTok{data =}\NormalTok{ gala, }
    \DataTypeTok{family =} \KeywordTok{poisson}\NormalTok{(}\DataTypeTok{link =}\NormalTok{ log))}
\KeywordTok{summary}\NormalTok{(fit)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## glm(formula = Species ~ Area + Elevation + Nearest + Scruz + 
##     Adjacent, family = poisson(link = log), data = gala)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -8.2752  -4.4966  -0.9443   1.9168  10.1849  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.155e+00  5.175e-02  60.963  < 2e-16 ***
## Area        -5.799e-04  2.627e-05 -22.074  < 2e-16 ***
## Elevation    3.541e-03  8.741e-05  40.507  < 2e-16 ***
## Nearest      8.826e-03  1.821e-03   4.846 1.26e-06 ***
## Scruz       -5.709e-03  6.256e-04  -9.126  < 2e-16 ***
## Adjacent    -6.630e-04  2.933e-05 -22.608  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 3510.73  on 29  degrees of freedom
## Residual deviance:  716.85  on 24  degrees of freedom
## AIC: 889.68
## 
## Number of Fisher Scoring iterations: 5
\end{verbatim}

Let \(\boldsymbol{\beta}\) be a \(6 \times 1\) column vector with the
regression coefficients (intercept included), and let
\(\hat{\boldsymbol{\beta}}\) be the maximum likelihood estimator for
\(\boldsymbol{\beta}\).

\begin{itemize}
\tightlist
\item
  Write down the asymptotic distribution for
  \(\hat{\boldsymbol{\beta}}\), and specify how the covariance matrix
  for \(\hat{\boldsymbol{\beta}}\) is estimated.
\end{itemize}

We will focus on the effect of \texttt{Elevation}, and denote the
corresponding regression coefficient \(\beta_2\).

\begin{itemize}
\tightlist
\item
  Write down the maximum likelihood estimate for \(\beta_2\) in the
  print-out above.
\item
  How can you explain this value to a biologist interested in
  understanding the effect of \texttt{Elevation} on the number of
  species of tortoise found on the islands?
\item
  What is numerical value for the estimated standard deviation of
  \(\hat{\beta}_2\) given in above?
\item
  Construct an approximate 95\% confidence interval for \(\beta_2\).
\end{itemize}

\subsection{Problem 3: Exam December 2017 from UiO, Problem
1.}\label{problem-3-exam-december-2017-from-uio-problem-1.}

(written out - with small changes to fit the notation we use - from
\url{https://www.uio.no/studier/emner/matnat/math/STK3100/h17/stk3100-4100_2017_2eng.pdf})

Assume that the random variable \(Y\) is Poisson distributed with
probability mass function (pmf)

\[\text{P}(Y = y | \lambda)= \frac{\lambda^y}{y!} \exp(-\lambda), \ y = 0, 1, 2, \dots.\]

\textbf{a)} Show that the distribution of \(Y\) is an exponential
family, that is, show that the pmf can be written on the form
\[\exp\left\{\frac{y \theta - b(\theta)}{\phi}w + c(y, \phi, w)\right\},\]
and determine \(\theta\), \(\phi\), \(w\), \(b(\theta)\) and
\(c(y, \phi, w)\).

We then assume that \(Y_1, Y_2, \dots, Y_n\) are independent with the
pmf from a), and let \(\mu_i = \text{E}(Y_i)\), \(i = 1, \dots, n\).

\textbf{b)} Explain what we mean by a generalized linear model (GLM) for
\(Y_1, Y_2, \dots, Y_n\) with link function \(g(\cdot)\), and determine
the canonical link function.

\textbf{c)} Derive an expression for the log-likelihood function
\(l(\mathbf{\mu}; \mathbf{y})\), where
\(\mathbf{y} = (y_1, \dots, y_n)^T\) is the observed value of
\(\mathbf{Y} = (Y_1, Y_2, \dots, Y_n)^T\) and
\(\mathbf{\mu} = (\mu_1, \dots, \mu_n)^T\).

\textbf{d)} Explain what we mean by a saturated model and determine the
maximum of \(l(\mathbf{\mu}; \mathbf{y})\) for the saturated model.

\textbf{e)} Explain what we mean by the deviance
\(D(\mathbf{y}; \mathbf{\hat{\mu}})\) of a Poisson GLM, find an
expression for the deviance, and discuss how it may be used.

\textbf{ SECOND WEEK}

\section{Poisson regression for count
data}\label{poisson-regression-for-count-data}

What did we do last week?

Examples --- GLM model --- loglikelihood, score function and Fisher
information matrix --- asympototic results for
\(\hat{\boldsymbol{\beta}}\) and Wald, score and LRT.

\section{Residuals}\label{residuals}

Two types of residuals are popular: \emph{deviance} and \emph{Pearson}.

\subsection{Deviance residuals}\label{deviance-residuals}

Recall that the deviance was defined by relating a \emph{candidate}
model to a \emph{saturated model} and calculating the likelihood ratio
statistic with these two models.

\textbf{Saturated model:} If we were to provide a perfect fit to our
data then we would estimate the mean \(\lambda_i\) by the observed count
for observation \(i\). That is, \(\tilde{\lambda}_i=y_i\) Then
\(\tilde{\lambda}\) is an \(n\) dimensional column vector with the
elements \(\tilde{\lambda}_i\).

\textbf{Candidate model:} The model that we are investigated can be
thought of as a \emph{candidate} model. Then we maximize the likelihood
and get \(\hat{\beta}\) which through our linear predictor and link
function we turn into \(\hat{\lambda}_i\), also called \(\hat{y}_i\).
Then \(\hat{\lambda}\) is an \(n\)-dimensional column vector with the
elements \(\hat{\lambda}_i\), also called \(\hat{y}\).

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\[D=-2(\ln L(\text{candidate model})-\ln L(\text{saturated model}))=-2(l(\hat{\lambda})-l(\tilde{\lambda}))=
-2\sum_{i=1}^n(l_i(\hat{\lambda}_i)-l_i(\tilde{\lambda}_i))\] Inserting
the Poisson likelihood and \(\lambda_i\) estimates this gives:

\[ D=2 \sum_{i=1}^n [y_i\ln(\frac{y_i}{\hat{y}_i})-(y_i-\hat{y}_i)]\]
Where \(\hat{y}_i=\exp({\bf x}_i^T \hat{\beta})\). Verify this by
yourself.

The deviance residuals are given by a signed version of each element in
the sum for the deviance, that is
\[d_i=\text{sign}(y_i-\hat{y}_i)\cdot \left\{ 
2[y_i\ln(\frac{y_i}{\hat{y}_i})-(y_i-\hat{y}_i)]\right\}^{1/2}\] where
the term \(\text{sign}(y_i-\hat{y}_i)\) makes negative residuals
possible - and we get the same sign as the \emph{Pearson residuals}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Pearson residuals}\label{pearson-residuals}

The Pearson residuals are given as \[ r_i=\frac{o_i-e_i}{\sqrt{e_i}}\]
where \(o_i\) is the observed count for observation \(i\) and \(e_i\) is
the estimated expected count for observation \(i\). We have that
\(o_i=y_i\) and
\(e_i=\hat{y}_i=\hat{\lambda}_i=\exp({\bf x}_i^T\hat{\beta})\).

Remark: A standardized version scales the Pearson residuals with
\(\sqrt{1-h_{ii}}\) similar to the standardized residuals for the normal
model. Here \(h_{ii}\) is the diagonal element number \(i\) in the hat
matrix \({\bf H}={\bf X}({\bf X}^T{\bf X})^{-1}{\bf X}^T\).

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Plotting residuals}\label{plotting-residuals}

Deviance and Pearson residuals can be used for checking the fit of the
model, by plotting the residuals against fitted values and covariates.
Normality of residuals are also assumed, and can be checked using
qq-plots as for the MLR in Module 2.

Below - notice the trend in the residuals, this is due to the discrete
nature of the response. The plot with different shades of blue shows
that the structures are for equal values of \(y\).

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{model3 =}\StringTok{ }\KeywordTok{glm}\NormalTok{(Sa }\OperatorTok{~}\StringTok{ }\NormalTok{W }\OperatorTok{+}\StringTok{ }\NormalTok{C, }\DataTypeTok{family =} \KeywordTok{poisson}\NormalTok{(}\DataTypeTok{link =}\NormalTok{ log), }\DataTypeTok{data =}\NormalTok{ crab, }\DataTypeTok{contrasts =} \KeywordTok{list}\NormalTok{(}\DataTypeTok{C =} \StringTok{"contr.sum"}\NormalTok{, }
    \DataTypeTok{S =} \StringTok{"contr.sum"}\NormalTok{))}
\NormalTok{df =}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{Sa =}\NormalTok{ crab}\OperatorTok{$}\NormalTok{Sa, }\DataTypeTok{fitted =}\NormalTok{ model3}\OperatorTok{$}\NormalTok{fitted.values, }\DataTypeTok{dres =} \KeywordTok{residuals}\NormalTok{(model3, }
    \DataTypeTok{type =} \StringTok{"deviance"}\NormalTok{), }\DataTypeTok{pres =} \KeywordTok{residuals}\NormalTok{(model3, }\DataTypeTok{type =} \StringTok{"pearson"}\NormalTok{))}
\end{Highlighting}
\end{Shaded}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(ggplot2)}

\NormalTok{gg1 =}\StringTok{ }\KeywordTok{ggplot}\NormalTok{(df) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ fitted, }\DataTypeTok{y =}\NormalTok{ dres, }\DataTypeTok{color =} \StringTok{"deviance"}\NormalTok{)) }\OperatorTok{+}\StringTok{ }
\StringTok{    }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ fitted, }\DataTypeTok{y =}\NormalTok{ pres, }\DataTypeTok{color =} \StringTok{"pearson"}\NormalTok{))}
\NormalTok{gg1}
\end{Highlighting}
\end{Shaded}

\includegraphics{4PoissonGamma_files/figure-latex/unnamed-chunk-11-1.pdf}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{gg2 =}\StringTok{ }\KeywordTok{ggplot}\NormalTok{(df) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ fitted, }\DataTypeTok{y =}\NormalTok{ pres, }\DataTypeTok{color =}\NormalTok{ Sa))}
\NormalTok{gg2}
\end{Highlighting}
\end{Shaded}

\includegraphics{4PoissonGamma_files/figure-latex/unnamed-chunk-12-1.pdf}

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{# qqnorm(residuals(model3,type='deviance'))}
\CommentTok{# qqline(residuals(model3,type='deviance'))}
\CommentTok{# qqnorm(residuals(model3,type='pearson'))}
\CommentTok{# qqline(residuals(model3,type='pearson'))}
\end{Highlighting}
\end{Shaded}

Pearson residuals

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{dff =}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{devres =} \KeywordTok{residuals}\NormalTok{(model3, }\DataTypeTok{type =} \StringTok{"deviance"}\NormalTok{), }\DataTypeTok{pearsonres =} \KeywordTok{residuals}\NormalTok{(model3, }
    \DataTypeTok{type =} \StringTok{"pearson"}\NormalTok{))}
\KeywordTok{ggplot}\NormalTok{(dff, }\KeywordTok{aes}\NormalTok{(}\DataTypeTok{sample =}\NormalTok{ devres)) }\OperatorTok{+}\StringTok{ }\KeywordTok{stat_qq}\NormalTok{(}\DataTypeTok{pch =} \DecValTok{19}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_abline}\NormalTok{(}\DataTypeTok{intercept =} \DecValTok{0}\NormalTok{, }
    \DataTypeTok{slope =} \DecValTok{1}\NormalTok{, }\DataTypeTok{linetype =} \StringTok{"dotted"}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{labs}\NormalTok{(}\DataTypeTok{x =} \StringTok{"Theoretical quantiles"}\NormalTok{, }\DataTypeTok{y =} \StringTok{"Deviance residuals"}\NormalTok{, }
    \DataTypeTok{title =} \StringTok{"Normal Q-Q"}\NormalTok{, }\DataTypeTok{subtitle =} \KeywordTok{deparse}\NormalTok{(model3}\OperatorTok{$}\NormalTok{call))}
\end{Highlighting}
\end{Shaded}

\includegraphics{4PoissonGamma_files/figure-latex/unnamed-chunk-13-1.pdf}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{ggplot}\NormalTok{(dff, }\KeywordTok{aes}\NormalTok{(}\DataTypeTok{sample =}\NormalTok{ pearsonres)) }\OperatorTok{+}\StringTok{ }\KeywordTok{stat_qq}\NormalTok{(}\DataTypeTok{pch =} \DecValTok{19}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_abline}\NormalTok{(}\DataTypeTok{intercept =} \DecValTok{0}\NormalTok{, }
    \DataTypeTok{slope =} \DecValTok{1}\NormalTok{, }\DataTypeTok{linetype =} \StringTok{"dotted"}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{labs}\NormalTok{(}\DataTypeTok{x =} \StringTok{"Theoretical quantiles"}\NormalTok{, }\DataTypeTok{y =} \StringTok{"Pearson residuals"}\NormalTok{, }
    \DataTypeTok{title =} \StringTok{"Normal Q-Q"}\NormalTok{, }\DataTypeTok{subtitle =} \KeywordTok{deparse}\NormalTok{(model3}\OperatorTok{$}\NormalTok{call))}
\end{Highlighting}
\end{Shaded}

\includegraphics{4PoissonGamma_files/figure-latex/unnamed-chunk-14-1.pdf}

\section{Model assessment and model
choice}\label{model-assessment-and-model-choice}

The fit of the model can be assessed based on goodness of fit statistics
(and related tests) and by residual plots. Model choice can be made from
analysis of deviance, or by comparing the AIC for different models.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Deviance test}\label{deviance-test}

We may use the deviance test presented in Module 3 to test if the model
under study is preferred compared to the saturated model.

We may write the deviance test as a sum of the squared deviance
residuals.

\[ D=2 \sum_{i=1}^n [y_i\ln(\frac{y_i}{\hat{y}_i})-(y_i-\hat{y}_i)]\]

Remark: if \(\sum_{i=1}^n y_i=\sum_{i=1}^n \hat{y}_i\) then deviance
residuals will be equal to
\[ D=2 \sum_{i=1}^n y_i\ln(\frac{y_i}{\hat{y}_i})\]

\textbf{Q}: is this the case for the log-linear model?

\textbf{A}: yes, but only if an intercept is included in the model?

The deviance statistic might be approximately \(\chi^2_{n-p}\), at least
when the counts are high.

Note: see below for remark.

\textbf{Remember} that the likelihood ratio test can be performed using
the difference between two deviances.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Pearson test}\label{pearson-test}

The Pearson \(\chi^2\)-goodness of fit statistic is given as the sum of
the squared Pearson residuals

\[ X_P^2=\sum_{i=1}^n r_i^2=\sum_{i=1}^n  \frac{(y_j-\hat{y}_i)^2}{\hat{y}_i}\]
where \(\hat{y}_i=\hat{\lambda}_i=\exp({\bf x}_i^T\hat{\beta})\). The
Pearson \(\chi^2\) statistic is asymptotically equivalent to the
deviance statistic and thus is asymptotically \(\chi^2_{n-p}\).

\textbf{Remark}: See connection with Pearson residuals.

\textbf{Remark}: the asymptotic distribution of both statistics
(deviance and Pearson) are questionable when there are many low counts.
Agresti (1996, page 990) suggest analysing grouped data, for example by
grouping by width in the crab example.

\textbf{Remark}: The Pearson statistic is also used for testing
independence in contingency tables - we will do that in Compulsory
Exercise 2.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Example: goodness of fit with female
crabs}\label{example-goodness-of-fit-with-female-crabs}

\textbf{Q}: Comment on the print-out. Is this a good fit? What might a
bad fit be due to?

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{model3 =}\StringTok{ }\KeywordTok{glm}\NormalTok{(Sa }\OperatorTok{~}\StringTok{ }\NormalTok{W }\OperatorTok{+}\StringTok{ }\NormalTok{C, }\DataTypeTok{family =} \KeywordTok{poisson}\NormalTok{(}\DataTypeTok{link =}\NormalTok{ log), }\DataTypeTok{data =}\NormalTok{ crab, }\DataTypeTok{contrasts =} \KeywordTok{list}\NormalTok{(}\DataTypeTok{C =} \StringTok{"contr.sum"}\NormalTok{))}
\KeywordTok{summary}\NormalTok{(model3)}
\DecValTok{1} \OperatorTok{-}\StringTok{ }\KeywordTok{pchisq}\NormalTok{(model3}\OperatorTok{$}\NormalTok{deviance, model3}\OperatorTok{$}\NormalTok{df.residual)}
\NormalTok{Xp =}\StringTok{ }\KeywordTok{sum}\NormalTok{(}\KeywordTok{residuals}\NormalTok{(model3, }\DataTypeTok{type =} \StringTok{"pearson"}\NormalTok{)}\OperatorTok{^}\DecValTok{2}\NormalTok{)}
\NormalTok{Xp}
\DecValTok{1} \OperatorTok{-}\StringTok{ }\KeywordTok{pchisq}\NormalTok{(Xp, model3}\OperatorTok{$}\NormalTok{df.residual)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## glm(formula = Sa ~ W + C, family = poisson(link = log), data = crab, 
##     contrasts = list(C = "contr.sum"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.0415  -1.9581  -0.5575   0.9830   4.7523  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.92089    0.56010  -5.215 1.84e-07 ***
## W            0.14934    0.02084   7.166 7.73e-13 ***
## C1           0.27085    0.11784   2.298   0.0215 *  
## C2           0.07117    0.07296   0.975   0.3294    
## C3          -0.16551    0.09316  -1.777   0.0756 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 632.79  on 172  degrees of freedom
## Residual deviance: 559.34  on 168  degrees of freedom
## AIC: 924.64
## 
## Number of Fisher Scoring iterations: 6
## 
## [1] 0
## [1] 543.249
## [1] 0
\end{verbatim}

\begin{verbatim}
## **A:** Reject the null that the model is good. This may be due to wrong model or missing covariate, or overdispersion.
\end{verbatim}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{AIC}\label{aic}

Identical to Module 3 - we may use the Akaike informations criterion.
Let \(p\) be the number of regression parameters in our model.
\[\text{AIC} =-2 \cdot l(\hat{\beta})+2p\] A scaled version of AIC,
standardizing for sample size, is sometimes preferred. And, we may also
use the BIC, where \(2p\) is replaced by \(\log(n)\cdot p\).

\subsection{Analysis of deviance}\label{analysis-of-deviance}

Identical to Module 3 we may also sequentially compare models, and use
analysis of deviance for this.

\section{Overdispersion}\label{overdispersion}

Count data might show greater variability in the response counts than we
would expect if the response followed a Poisson distribution. This is
called \emph{overdispersion}.

Example: newspaper sales with tourist bus.

Our model states that the variance \(\text{Var}(Y_i)=\lambda_i\). If we
change the model to \(\text{Var}(Y_i)=\phi\lambda_i\) we may allow for
an increased variance due to heterogeniety among subjects.

Agresti (1996, page 92) explains: ``For our crab data set, what if
width, weight, colour and spine affect the number of satellites for a
female crab, and we only fitted a model with width as covariate. Then
the crabs with a certain width are a mixture of crabs of various
weights, colours and spine condition - that is, a mixture of several
Poisson populations, each with its own mean for the response. This
heterogeniety may give an overall response distribution where the
variance is greater than the standard Poisson variance.''

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

The overdispersion parameter can be estimated as the average Pearson
statistic or average deviance

\[\hat{\phi}_D = \frac{1}{n-p} D\]

where \(D\) is the deviance. Note that similarity to
\(\hat{\sigma^2} = 1/(n-p)\cdot\text{SSE}\) in the MLR. The
Cov\((\hat{\beta})\) can then be changed to
\(\hat{\phi}F^{-1}(\hat{\beta})\), so we multiply the standard error by
the square root of \(\hat{\phi}_D\).

Remark: We are now moving from likelihood to quasi-likelihood theory,
where only E\((Y_j)\) and Var\((Y_j)\) - and not the distribution of
\(Y_j\) - are used in the estimation.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{model.disp =}\StringTok{ }\KeywordTok{glm}\NormalTok{(Sa }\OperatorTok{~}\StringTok{ }\NormalTok{W, }\DataTypeTok{family =} \KeywordTok{quasipoisson}\NormalTok{(}\DataTypeTok{link =}\NormalTok{ log), }\DataTypeTok{data =}\NormalTok{ crab)}
\KeywordTok{summary.glm}\NormalTok{(model.disp)}
\KeywordTok{summary.glm}\NormalTok{(model.disp)}\OperatorTok{$}\NormalTok{dispersion}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## glm(formula = Sa ~ W, family = quasipoisson(link = log), data = crab)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.8526  -1.9884  -0.4933   1.0970   4.9221  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.30476    0.96729  -3.417 0.000793 ***
## W            0.16405    0.03562   4.606 7.99e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 3.182205)
## 
##     Null deviance: 632.79  on 172  degrees of freedom
## Residual deviance: 567.88  on 171  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 6
## 
## [1] 3.182205
\end{verbatim}

\section{Rate models and offset}\label{rate-models-and-offset}

In the Poisson process we might analyse an event that occurs within a
time interval or region in space, and therefore it is often of interest
to model the \emph{rate} at which events occur.

Examples:

\begin{itemize}
\tightlist
\item
  crime rates in cities
\item
  death rate for smokers vs.~non-smokers
\item
  rate of auto thefts in cities
\end{itemize}

Agresti (1996, page 86): ``what if we want to model the number of auto
thefts for a year in a sample of cities. We would make a rate for each
city by dividing the number of auto thefts by the population size of the
city. The model could then describe how this rate depends on
unemployment rate, median income, percentage of residents having
completed high school.''"

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

Now we don't want a model for \(Y_i\) but for \(Y_i/t_i\), where

\begin{itemize}
\tightlist
\item
  Let \(t_i\) denote the index (population size in the example)
  associated with observation \(i\).
\item
  We still assume that \(Y_i\) follows a Poisson distribution, but we
  now include the index in the modelling and focus on \(Y_i/t_i\).
\item
  The expected value of \(Y_i/t_i\) would then be
  \(\text{E}(Y_i)/t_i=\lambda_i/t_i\).
\end{itemize}

A log-linear model would be \[ \log(\lambda_i/t_i)={\bf x}_i^T \beta\]
We may equivalently write the model as
\[ \log(\lambda_i)-\log(t_i)={\bf x}_i^T \beta\] This adjustment term is
called an \emph{offset} and is a known quantity.

The expected number of outcomes will then satisfy
\[ \text{E}(Y_i)=\lambda_i=t_i \exp({\bf x}_i^T \beta).\]

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Example: British doctors and rate
models}\label{example-british-doctors-and-rate-models}

Count data - the number of times and event occurs - is common. In one
famous study British doctors were in 1951 sent a questionnaire about
whether they smoked tobacco - and later information about their death
were collected.

Research questions that were asked were: 1) Is the death rate higher for
smokers than for non-smokers? 2) If so, by how much? 3) And, how is this
related to age?

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(boot)}
\StringTok{`}\DataTypeTok{?}\StringTok{`}\NormalTok{(breslow)}
\CommentTok{# n=person-year, ns=smoker-years, age=midpoint 10 year age group, y=number}
\CommentTok{# of deaths due to cad, smoke=smoking status}
\KeywordTok{head}\NormalTok{(breslow, }\DataTypeTok{n =} \DecValTok{10}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##    age smoke     n   y    ns
## 1   40     0 18790   2     0
## 2   50     0 10673  12     0
## 3   60     0  5710  28     0
## 4   70     0  2585  28     0
## 5   80     0  1462  31     0
## 6   40     1 52407  32 52407
## 7   50     1 43248 104 43248
## 8   60     1 28612 206 28612
## 9   70     1 12663 186 12663
## 10  80     1  5317 102  5317
\end{verbatim}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

To investigate this we will look at different ways of relating the
expected number of deaths and the number of doctors at risk in the
observation period for each smoke and age group. The aim is to model the
rate of cardiovascular mortality.

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{# first age and smoke (but not interaction thereof)}
\NormalTok{fit1 <-}\StringTok{ }\KeywordTok{glm}\NormalTok{(y }\OperatorTok{~}\StringTok{ }\KeywordTok{factor}\NormalTok{(age) }\OperatorTok{+}\StringTok{ }\KeywordTok{factor}\NormalTok{(smoke), }\DataTypeTok{offset =} \KeywordTok{log}\NormalTok{(n), }\DataTypeTok{family =}\NormalTok{ poisson, }
    \DataTypeTok{data =}\NormalTok{ breslow)}
\KeywordTok{summary}\NormalTok{(fit1)}

\CommentTok{# do we need interaction?}
\NormalTok{fit2 <-}\StringTok{ }\KeywordTok{update}\NormalTok{(fit1, . }\OperatorTok{~}\StringTok{ }\NormalTok{. }\OperatorTok{+}\StringTok{ }\KeywordTok{factor}\NormalTok{(smoke) }\OperatorTok{*}\StringTok{ }\KeywordTok{factor}\NormalTok{(age))}
\KeywordTok{summary}\NormalTok{(fit2)}
\KeywordTok{anova}\NormalTok{(fit1, fit2, }\DataTypeTok{test =} \StringTok{"Chisq"}\NormalTok{)}
\CommentTok{# yes, significant interaction between age and smoking - how does this}
\CommentTok{# compare to a deviance test for fit1?}

\CommentTok{# reporting on final model - give rate in each possible group}
\KeywordTok{cbind}\NormalTok{(fit2}\OperatorTok{$}\NormalTok{fitted.values, breslow}\OperatorTok{$}\NormalTok{y)  }\CommentTok{#perfect fit since number of coeffs equal number of groups}
\KeywordTok{exp}\NormalTok{(}\KeywordTok{predict}\NormalTok{(fit2, }\DataTypeTok{type =} \StringTok{"link"}\NormalTok{))}
\KeywordTok{predict}\NormalTok{(fit2, }\DataTypeTok{type =} \StringTok{"response"}\NormalTok{)}

\CommentTok{# I want to see the estimated value of lambda for a population size of 1 and}
\CommentTok{# of 1000}
\KeywordTok{length}\NormalTok{(fit2}\OperatorTok{$}\NormalTok{coefficients)}
\CommentTok{# year 40 nonsmokers should only be the intercept}
\KeywordTok{exp}\NormalTok{(fit2}\OperatorTok{$}\NormalTok{coefficients[}\DecValTok{1}\NormalTok{])  }\CommentTok{# expected number of deaths pr individual in a population of 40-year olds who do not smoke}
\CommentTok{# pr 1000}
\KeywordTok{exp}\NormalTok{(fit2}\OperatorTok{$}\NormalTok{coefficients[}\DecValTok{1}\NormalTok{]) }\OperatorTok{*}\StringTok{ }\DecValTok{1000}
\CommentTok{# 80 year olds who smoke}
\KeywordTok{exp}\NormalTok{(}\KeywordTok{sum}\NormalTok{(fit2}\OperatorTok{$}\NormalTok{coefficients[}\KeywordTok{c}\NormalTok{(}\DecValTok{1}\NormalTok{, }\DecValTok{5}\NormalTok{, }\DecValTok{6}\NormalTok{, }\DecValTok{10}\NormalTok{)]))}
\CommentTok{# pr 1000}
\DecValTok{1000} \OperatorTok{*}\StringTok{ }\KeywordTok{exp}\NormalTok{(}\KeywordTok{sum}\NormalTok{(fit2}\OperatorTok{$}\NormalTok{coefficients[}\KeywordTok{c}\NormalTok{(}\DecValTok{1}\NormalTok{, }\DecValTok{5}\NormalTok{, }\DecValTok{6}\NormalTok{, }\DecValTok{10}\NormalTok{)]))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## glm(formula = y ~ factor(age) + factor(smoke), family = poisson, 
##     data = breslow, offset = log(n))
## 
## Deviance Residuals: 
##        1         2         3         4         5         6         7  
## -2.17978  -1.30800  -0.13791   0.22882   1.91902   0.90160   0.51038  
##        8         9        10  
##  0.05135  -0.08732  -0.91237  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -7.9193     0.1918 -41.298  < 2e-16 ***
## factor(age)50    1.4840     0.1951   7.606 2.82e-14 ***
## factor(age)60    2.6275     0.1837  14.301  < 2e-16 ***
## factor(age)70    3.3505     0.1848  18.131  < 2e-16 ***
## factor(age)80    3.7001     0.1922  19.249  < 2e-16 ***
## factor(smoke)1   0.3545     0.1074   3.302  0.00096 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 935.067  on 9  degrees of freedom
## Residual deviance:  12.132  on 4  degrees of freedom
## AIC: 79.2
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## Call:
## glm(formula = y ~ factor(age) + factor(smoke) + factor(age):factor(smoke), 
##     family = poisson, data = breslow, offset = log(n))
## 
## Deviance Residuals: 
##  [1]  0  0  0  0  0  0  0  0  0  0
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   -9.1479     0.7071 -12.937  < 2e-16 ***
## factor(age)50                  2.3574     0.7638   3.087  0.00203 ** 
## factor(age)60                  3.8302     0.7319   5.233 1.67e-07 ***
## factor(age)70                  4.6227     0.7319   6.316 2.69e-10 ***
## factor(age)80                  5.2944     0.7296   7.257 3.96e-13 ***
## factor(smoke)1                 1.7469     0.7289   2.397  0.01654 *  
## factor(age)50:factor(smoke)1  -0.9866     0.7901  -1.249  0.21174    
## factor(age)60:factor(smoke)1  -1.3628     0.7562  -1.802  0.07151 .  
## factor(age)70:factor(smoke)1  -1.4423     0.7565  -1.906  0.05659 .  
## factor(age)80:factor(smoke)1  -1.8470     0.7572  -2.439  0.01471 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 9.3507e+02  on 9  degrees of freedom
## Residual deviance: 2.9754e-14  on 0  degrees of freedom
## AIC: 75.068
## 
## Number of Fisher Scoring iterations: 3
## 
## Analysis of Deviance Table
## 
## Model 1: y ~ factor(age) + factor(smoke)
## Model 2: y ~ factor(age) + factor(smoke) + factor(age):factor(smoke)
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1         4     12.132                       
## 2         0      0.000  4   12.132  0.01639 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##    [,1] [,2]
## 1     2    2
## 2    12   12
## 3    28   28
## 4    28   28
## 5    31   31
## 6    32   32
## 7   104  104
## 8   206  206
## 9   186  186
## 10  102  102
##   1   2   3   4   5   6   7   8   9  10 
##   2  12  28  28  31  32 104 206 186 102 
##   1   2   3   4   5   6   7   8   9  10 
##   2  12  28  28  31  32 104 206 186 102 
## [1] 10
##  (Intercept) 
## 0.0001064396 
## (Intercept) 
##   0.1064396 
## [1] 0.01918375
## [1] 19.18375
\end{verbatim}

\section{Modelling continuous positive response
data}\label{modelling-continuous-positive-response-data}

\subsection{Examples of continuous positive
responses}\label{examples-of-continuous-positive-responses}

\begin{itemize}
\tightlist
\item
  Insurance: Claim sizes
\item
  Medicine: Time to blood coagulation (main example)
\item
  Biology: Time in various development stages for fruit fly
\item
  Meteorology: Amount of precipitation (interactive session - exam
  question 2012)
\end{itemize}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Models for continuous positive
responses}\label{models-for-continuous-positive-responses}

\begin{itemize}
\tightlist
\item
  Lognormal distribution on response
\item
  Gamma distribution on response
\item
  Inverse Gaussian distribution on response (we will not consider this
  here)
\end{itemize}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Time to blood coagulation}\label{time-to-blood-coagulation}

This data is described in McCullagh and Nelder (1989, page 300). The
data represents clotting time of blood (in seconds) \texttt{y} for
normal plasma diluted to nine different percentage concentrations
\texttt{u} with socalled prothrombin-free plasma. To induce the clotting
a chemical called thromboplasting was used, and in the experiment two
different lots of the chemical were used - denoted \texttt{lot}. Our aim
is to investigate the relationship between the clotting time and the
dilution percentage, and look at differences between the lots.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{clot =}\StringTok{ }\KeywordTok{read.table}\NormalTok{(}\StringTok{"https://www.math.ntnu.no/emner/TMA4315/2018h/clot.txt"}\NormalTok{, }\DataTypeTok{header =}\NormalTok{ T)}
\NormalTok{clot}\OperatorTok{$}\NormalTok{lot =}\StringTok{ }\KeywordTok{as.factor}\NormalTok{(clot}\OperatorTok{$}\NormalTok{lot)}
\KeywordTok{summary}\NormalTok{(clot)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##        u            time       lot  
##  Min.   :  5   Min.   : 12.0   1:9  
##  1st Qu.: 15   1st Qu.: 18.0   2:9  
##  Median : 30   Median : 23.0        
##  Mean   : 40   Mean   : 32.5        
##  3rd Qu.: 60   3rd Qu.: 35.0        
##  Max.   :100   Max.   :118.0
\end{verbatim}

\section{Lognormal distribution}\label{lognormal-distribution}

Let \(Y_i\) be the response on the original scale, where \(Y_i>0\).

Transform the response to a logaritmic scale: \(Y^*_i=\ln(Y_i)\). Then,
assume that transformed responses follow a normal distribution (or
follows approximately) and use ordinary MLR. This means we have a GLM
with normal response and identity link (on logarithmic scale of
reponse).

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
  \(Y^*_i \sim N(\mu^{*}_i,\sigma^{*2})\)
\item
  \(\eta_i={\bf x}_i^T\beta\)
\item
  \(\mu^{*}_i=\eta_i\) (identity link)
\end{enumerate}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

There are two ways of looking at this,

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
  either this is just a transformation to achieve approximate normality,
  or
\item
  we assume that the original data follows a lognormal distribution.
\end{enumerate}

In genomics one usually assume the former, and reports back results on
the exponential scale - just say that the mean of original data is
\(\exp(\mu^{*}_i)\).

However, if on instead assume that the original data really comes from a
lognormal distribution, then it can be shown that

\[\text{E}(Y_i)=\exp(\mu^{*}_i) \cdot \exp(\sigma^{*2}/2)\]
\[\text{Var}(Y_i) =\exp(\sigma^{*2} -1)\cdot \mu_i^2\]

i.e.~standard deviation proportional to expectation. That is in general
a put-off for using the lognormal model.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{orgmu =}\StringTok{ }\DecValTok{1}
\NormalTok{orgsd =}\StringTok{ }\FloatTok{0.3}
\KeywordTok{library}\NormalTok{(ggplot2)}
\NormalTok{xrange =}\StringTok{ }\KeywordTok{range}\NormalTok{(}\KeywordTok{rlnorm}\NormalTok{(}\DecValTok{1000}\NormalTok{, orgmu, orgsd))}
\KeywordTok{ggplot}\NormalTok{(}\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ xrange), }\KeywordTok{aes}\NormalTok{(xrange)) }\OperatorTok{+}\StringTok{ }\KeywordTok{xlab}\NormalTok{(}\KeywordTok{expression}\NormalTok{(x)) }\OperatorTok{+}\StringTok{ }\KeywordTok{stat_function}\NormalTok{(}\DataTypeTok{fun =}\NormalTok{ dlnorm, }
    \DataTypeTok{args =} \KeywordTok{list}\NormalTok{(}\DataTypeTok{meanlog =}\NormalTok{ orgmu, }\DataTypeTok{sdlog =}\NormalTok{ orgsd), }\DataTypeTok{geom =} \StringTok{"line"}\NormalTok{, }\DataTypeTok{colour =} \StringTok{"red"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\includegraphics{4PoissonGamma_files/figure-latex/unnamed-chunk-21-1.pdf}

\section{Gamma regression}\label{gamma-regression}

\subsection{The gamma distribution}\label{the-gamma-distribution}

We have seen that a gamma distributed variable may be the result of the
time between events in a Poisson process. The well known
\(\chi^2_{\delta}\)-distribution is a special case of the gamma
distribution (\(\frac{\nu}{\mu_i}=2\), \(\nu=\frac{\delta}{2}\)).

There are many parameterization for the gamma distribution, but we will
stick with the one used in our textbook (page 643):

\(Y_i \sim Ga(\mu_i,\nu)\) with density
\[ f(y_i)=\frac{1}{\Gamma(\nu)} (\frac{\nu}{\mu_i})^{\nu} y_i^{\nu-1}\exp(-\frac{\nu}{\mu_i}y_i) \text{ for }y_i>0\]

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{mu =}\StringTok{ }\DecValTok{1}
\NormalTok{nu =}\StringTok{ }\FloatTok{0.3}
\KeywordTok{library}\NormalTok{(ggplot2)}
\NormalTok{xrange =}\StringTok{ }\KeywordTok{range}\NormalTok{(}\KeywordTok{rgamma}\NormalTok{(}\DecValTok{1000}\NormalTok{, }\DataTypeTok{shape =}\NormalTok{ mu}\OperatorTok{/}\NormalTok{nu, }\DataTypeTok{scale =}\NormalTok{ nu))}
\KeywordTok{ggplot}\NormalTok{(}\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ xrange), }\KeywordTok{aes}\NormalTok{(xrange)) }\OperatorTok{+}\StringTok{ }\KeywordTok{xlab}\NormalTok{(}\KeywordTok{expression}\NormalTok{(x)) }\OperatorTok{+}\StringTok{ }\KeywordTok{stat_function}\NormalTok{(}\DataTypeTok{fun =}\NormalTok{ dgamma, }
    \DataTypeTok{args =} \KeywordTok{list}\NormalTok{(}\DataTypeTok{shape =}\NormalTok{ mu}\OperatorTok{/}\NormalTok{nu, }\DataTypeTok{scale =}\NormalTok{ nu), }\DataTypeTok{geom =} \StringTok{"line"}\NormalTok{, }\DataTypeTok{colour =} \StringTok{"red"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\includegraphics{4PoissonGamma_files/figure-latex/unnamed-chunk-22-1.pdf}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

We found in Module 1, see
\href{https://www.math.ntnu.no/emner/TMA4315/2017h/Module1ExponentialFamily.pdf}{exponential
family} that the gamma distribution is an exponential family, with

\begin{itemize}
\tightlist
\item
  \(\theta_i=-\frac{1}{\mu_i}\) is the canonical parameter
\item
  \(\phi=\frac{1}{\nu}\),
\item
  \(w_i=1\)
\item
  \(b(\theta_i)=-\ln(-\theta_i)\)
\item
  \(\text{E}(Y_i)=b'(\theta_i)=-\frac{1}{\theta_i}=\mu_i\)
\item
  \(\text{Var}(Y_i)=b''(\theta_i)\frac{\psi}{w_i}=\frac{\mu_i^2}{\nu}\)
\end{itemize}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

For a GLM model we have canonical link if \[\theta_i=\eta_i\] Since
\(\eta_i=g(\mu_i)\) this means to us that we need
\[\theta_i=g(\mu_i)=-\frac{1}{\mu_i}\] saying that with the canonical
link is \(-\frac{1}{\mu_i}\).

However, the most commonly used link is \(g(\mu_i)=\ln(\mu_i)\), and the
identity link is also used.

\textbf{Q:} Discuss the implications on \(\eta_i\) when using the
canonical link. What about using log-link?

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

Remark: often the inverse and not the negative inverse is used, and
since \[g(\mu_i)=-\frac{1}{\mu_i}={\bf x}_i^T\beta\] then
\[\frac{1}{\mu_i}=-{\bf x}_i^T\beta={\bf x}_i^T\beta^*\] where
\(\beta^*=-\beta\).

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Gamma GLM model}\label{gamma-glm-model}

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  \(Y_i \sim Ga(\mu_i,\nu)\)
\item
  \(\eta_i={\bf x}_i^T\beta\)
\item
  Popular link functions:
\end{enumerate}

\begin{itemize}
\tightlist
\item
  \(\eta_i=\mu_i\) (identity link)
\item
  \(\eta_i=\frac{1}{\mu_i}\) (inverse link)
\item
  \(\eta_i=\ln(\mu_i)\) (log-link)
\end{itemize}

\textbf{Remark:} In our model the parameter \(\mu_i\) varies with \(i\)
but \(\nu\) is the same for all observations.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsubsection{Example: Time to blood
coagulation}\label{example-time-to-blood-coagulation}

A simple model to start with is as follows (dosages often analysed on
log scale):

\tiny

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{fit1 =}\StringTok{ }\KeywordTok{glm}\NormalTok{(time }\OperatorTok{~}\StringTok{ }\NormalTok{lot }\OperatorTok{+}\StringTok{ }\KeywordTok{log}\NormalTok{(u), }\DataTypeTok{data =}\NormalTok{ clot, }\DataTypeTok{family =} \KeywordTok{Gamma}\NormalTok{(}\DataTypeTok{link =}\NormalTok{ log))}
\KeywordTok{summary}\NormalTok{(fit1)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## glm(formula = time ~ lot + log(u), family = Gamma(link = log), 
##     data = clot)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.17470  -0.11596  -0.04281   0.06919   0.27749  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.44660    0.13453   40.48  < 2e-16 ***
## lot2        -0.47034    0.07095   -6.63 8.02e-06 ***
## log(u)      -0.58476    0.03772  -15.50 1.22e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.02265072)
## 
##     Null deviance: 7.7087  on 17  degrees of freedom
## Residual deviance: 0.3211  on 15  degrees of freedom
## AIC: 104.28
## 
## Number of Fisher Scoring iterations: 5
\end{verbatim}

\normalsize

\textbf{Q}: describe what you see in the print-out.

\section{Gamma regression: likelihood and derivations
thereof}\label{gamma-regression-likelihood-and-derivations-thereof}

\textbf{Likelihood:}
\[L(\beta)=\prod_{i=1}^n \exp(-\frac{\nu y_i}{\mu_i}-\nu \ln \mu_i +\nu \ln \nu +(\nu-1)\ln y_i -\ln(\Gamma(\nu)))\]

\textbf{Log-likelihood:}
\[l(\beta)=\sum_{i=1}^n [-\frac{\nu y_i}{\mu_i}-\nu \ln \mu_i +\nu \ln \nu +(\nu-1)\ln y_i -\ln(\Gamma(\nu))]\]
Observe that we now- for the first time - have a nuisance parameter
\(\nu\) here.

To produce numerical estimates for the parameter of interest \(\beta\)
we may proceed to the score function, and solve using Newton Raphson or
Fisher scoring. If we do not have the canonical link the observed and
expected Fisher information matrix may not be equal.

What about \(\phi=1/\nu\)? Also estimated using maximum likelihood.

Further analyses: as before we use asymptotic distribution of parameter
estimates, and of Wald, LRT and score test.

\section{Scaled and unscaled
deviance}\label{scaled-and-unscaled-deviance}

We have defined the deviance as
\[D=-2(\ln L(\text{candidate model})-\ln L(\text{saturated model}))\]

This is often called the \emph{scaled deviance}.

The \emph{unscaled deviance} is then defined as \(\phi D\), but is sadly
sometimes also called the deviance - for example by \texttt{R}.

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
  For the normal model the
\end{enumerate}

\begin{itemize}
\tightlist
\item
  scaled deviance is
  \(D=\frac{1}{\sigma^2}\sum_{i=1}^n (y_i-\hat{\mu}_i)^2\), while
\item
  unscaled deviance is \(\phi D=\sum_{i=1}^n (y_i-\hat{\mu}_i)^2\)
\end{itemize}

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\setcounter{enumi}{1}
\item
  For the binomial and Poisson model \(\phi=1\) so the scaled and
  unscaled deviance are equal.
\item
  What about the Gamma model?
\end{enumerate}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

Some calculations - see IL week 2, problem 2: 1b.

\[D=\frac{-2 \sum_{i=1}^n [\ln (\frac{y_i}{\hat{\mu}_i})-\frac{y_i-\hat{\mu}_i}{\hat{\mu}_i}]}{\phi}\]
and unscaled as
\(\phi D=-2 \sum_{i=1}^n [\ln (\frac{y_i}{\hat{\mu}_i})-\frac{y_i-\hat{\mu}_i}{\hat{\mu}_i}]\).

Compare to print-out from R: the deviance in R is the \emph{unscaled
deviance}.

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{deviance}\NormalTok{(fit1)}
\NormalTok{nu1 =}\StringTok{ }\DecValTok{1}\OperatorTok{/}\KeywordTok{summary}\NormalTok{(fit1)}\OperatorTok{$}\NormalTok{dispersion}
\NormalTok{nu1}
\NormalTok{D =}\StringTok{ }\OperatorTok{-}\DecValTok{2} \OperatorTok{*}\StringTok{ }\NormalTok{nu1 }\OperatorTok{*}\StringTok{ }\KeywordTok{sum}\NormalTok{(}\KeywordTok{log}\NormalTok{(fit1}\OperatorTok{$}\NormalTok{y}\OperatorTok{/}\NormalTok{fit1}\OperatorTok{$}\NormalTok{fitted.values) }\OperatorTok{-}\StringTok{ }\NormalTok{((fit1}\OperatorTok{$}\NormalTok{y }\OperatorTok{-}\StringTok{ }\NormalTok{fit1}\OperatorTok{$}\NormalTok{fitted.values)}\OperatorTok{/}\NormalTok{fit1}\OperatorTok{$}\NormalTok{fitted.values))}
\NormalTok{D}
\KeywordTok{deviance}\NormalTok{(fit1) }\OperatorTok{*}\StringTok{ }\NormalTok{nu1}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] 0.3210963
## [1] 44.14871
## [1] 14.17599
## [1] 14.17599
\end{verbatim}

\section{Comparing models}\label{comparing-models}

\subsection{Comparing models based on
deviance}\label{comparing-models-based-on-deviance}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{fit2 =}\StringTok{ }\KeywordTok{glm}\NormalTok{(time }\OperatorTok{~}\StringTok{ }\NormalTok{lot }\OperatorTok{+}\StringTok{ }\KeywordTok{log}\NormalTok{(u) }\OperatorTok{+}\StringTok{ }\NormalTok{lot}\OperatorTok{:}\KeywordTok{log}\NormalTok{(u), }\DataTypeTok{data =}\NormalTok{ clot, }\DataTypeTok{family =} \KeywordTok{Gamma}\NormalTok{(}\DataTypeTok{link =}\NormalTok{ log))}
\KeywordTok{anova}\NormalTok{(fit1, fit2)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Analysis of Deviance Table
## 
## Model 1: time ~ lot + log(u)
## Model 2: time ~ lot + log(u) + lot:log(u)
##   Resid. Df Resid. Dev Df  Deviance
## 1        15    0.32110             
## 2        14    0.31576  1 0.0053352
\end{verbatim}

The deviance table does not include \(\phi\), so the unscaled deviance
is reported. If significance testing is done, the estimated \(\phi\)
from the largest model is used, and \(p\)-values are based on the scaled
deviance.

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{anova}\NormalTok{(fit1, fit2, }\DataTypeTok{test =} \StringTok{"Chisq"}\NormalTok{)}
\DecValTok{1} \OperatorTok{-}\StringTok{ }\KeywordTok{pchisq}\NormalTok{((}\KeywordTok{deviance}\NormalTok{(fit1) }\OperatorTok{-}\StringTok{ }\KeywordTok{deviance}\NormalTok{(fit2))}\OperatorTok{/}\KeywordTok{summary}\NormalTok{(fit2)}\OperatorTok{$}\NormalTok{dispersion, fit1}\OperatorTok{$}\NormalTok{df.residual }\OperatorTok{-}\StringTok{ }
\StringTok{    }\NormalTok{fit2}\OperatorTok{$}\NormalTok{df.residual)}
\KeywordTok{anova}\NormalTok{(fit1, fit2, }\DataTypeTok{test =} \StringTok{"F"}\NormalTok{)}
\DecValTok{1} \OperatorTok{-}\StringTok{ }\KeywordTok{pf}\NormalTok{((}\KeywordTok{deviance}\NormalTok{(fit1) }\OperatorTok{-}\StringTok{ }\KeywordTok{deviance}\NormalTok{(fit2))}\OperatorTok{/}\KeywordTok{summary}\NormalTok{(fit2)}\OperatorTok{$}\NormalTok{dispersion, fit1}\OperatorTok{$}\NormalTok{df.residual }\OperatorTok{-}\StringTok{ }
\StringTok{    }\NormalTok{fit2}\OperatorTok{$}\NormalTok{df.residual, fit2}\OperatorTok{$}\NormalTok{df.residual)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Analysis of Deviance Table
## 
## Model 1: time ~ lot + log(u)
## Model 2: time ~ lot + log(u) + lot:log(u)
##   Resid. Df Resid. Dev Df  Deviance Pr(>Chi)
## 1        15    0.32110                      
## 2        14    0.31576  1 0.0053352   0.6355
## [1] 0.6355477
## Analysis of Deviance Table
## 
## Model 1: time ~ lot + log(u)
## Model 2: time ~ lot + log(u) + lot:log(u)
##   Resid. Df Resid. Dev Df  Deviance      F Pr(>F)
## 1        15    0.32110                           
## 2        14    0.31576  1 0.0053352 0.2246 0.6429
## [1] 0.642854
\end{verbatim}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Comparing models based on
AIC}\label{comparing-models-based-on-aic}

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{AIC}\NormalTok{(fit1, fit2)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##      df      AIC
## fit1  4 104.2763
## fit2  5 105.9738
\end{verbatim}

\textbf{Q}: would you prefer \texttt{fit1} or \texttt{fit2}?

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

AIC can also be used when we compare models with different link
functions (models that are not nested).

The literature suggests to plot \(y_i\) vs.~each covariate to get a hint
about which link function or transformation to use.

\begin{itemize}
\tightlist
\item
  Identity: Plot of \(y_i\) vs \(x_i\) should be close to linear
\item
  \(\ln:\) Plot of \(\ln(y_i)\) vs \(x_i\) should be close to linear
\item
  Inverse (reciprocal): Plot of \(1/y_i\) vs \(x_i\) should be close to
  linear
\end{itemize}

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(ggplot2)}
\KeywordTok{library}\NormalTok{(ggpubr)}
\NormalTok{y =}\StringTok{ }\NormalTok{clot}\OperatorTok{$}\NormalTok{time}
\NormalTok{x =}\StringTok{ }\NormalTok{clot}\OperatorTok{$}\NormalTok{u}

\NormalTok{df =}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{y =}\NormalTok{ y, }\DataTypeTok{x =}\NormalTok{ x)}
\NormalTok{gg1 =}\StringTok{ }\KeywordTok{ggplot}\NormalTok{(df) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =} \KeywordTok{log}\NormalTok{(x), }\DataTypeTok{y =}\NormalTok{ y))}
\NormalTok{gg2 =}\StringTok{ }\KeywordTok{ggplot}\NormalTok{(df) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =} \KeywordTok{log}\NormalTok{(x), }\DataTypeTok{y =} \KeywordTok{log}\NormalTok{(y)))}
\NormalTok{gg3 =}\StringTok{ }\KeywordTok{ggplot}\NormalTok{(df) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =} \KeywordTok{log}\NormalTok{(x), }\DataTypeTok{y =} \DecValTok{1}\OperatorTok{/}\NormalTok{y))}
\NormalTok{gg4 =}\StringTok{ }\KeywordTok{ggplot}\NormalTok{(df) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =} \KeywordTok{sqrt}\NormalTok{(}\DecValTok{1}\OperatorTok{/}\NormalTok{x), }\DataTypeTok{y =} \KeywordTok{log}\NormalTok{(y)))}
\KeywordTok{ggarrange}\NormalTok{(gg1, gg2, gg3, gg4)}
\end{Highlighting}
\end{Shaded}

\includegraphics{4PoissonGamma_files/figure-latex/unnamed-chunk-28-1.pdf}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{fit4 =}\StringTok{ }\KeywordTok{glm}\NormalTok{(time }\OperatorTok{~}\StringTok{ }\NormalTok{lot }\OperatorTok{+}\StringTok{ }\KeywordTok{sqrt}\NormalTok{(}\DecValTok{1}\OperatorTok{/}\NormalTok{u), }\DataTypeTok{data =}\NormalTok{ clot, }\DataTypeTok{family =} \KeywordTok{Gamma}\NormalTok{(}\DataTypeTok{link =}\NormalTok{ log))}
\KeywordTok{AIC}\NormalTok{(fit1, fit4)}

\NormalTok{df4 =}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{fitted =}\NormalTok{ fit4}\OperatorTok{$}\NormalTok{fitted.values, }\DataTypeTok{dres =} \KeywordTok{residuals}\NormalTok{(fit4, }\DataTypeTok{type =} \StringTok{"deviance"}\NormalTok{))}
\NormalTok{gg4 =}\StringTok{ }\KeywordTok{ggplot}\NormalTok{(df4) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ fitted, }\DataTypeTok{y =}\NormalTok{ dres)) }\OperatorTok{+}\StringTok{ }\KeywordTok{scale_color_discrete}\NormalTok{(}\StringTok{""}\NormalTok{)}
\NormalTok{df1 =}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{fitted =}\NormalTok{ fit1}\OperatorTok{$}\NormalTok{fitted.values, }\DataTypeTok{dres =} \KeywordTok{residuals}\NormalTok{(fit1, }\DataTypeTok{type =} \StringTok{"deviance"}\NormalTok{))}
\NormalTok{gg1 =}\StringTok{ }\KeywordTok{ggplot}\NormalTok{(df1) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ fitted, }\DataTypeTok{y =}\NormalTok{ dres)) }\OperatorTok{+}\StringTok{ }\KeywordTok{scale_color_discrete}\NormalTok{(}\StringTok{""}\NormalTok{)}
\KeywordTok{ggarrange}\NormalTok{(gg1, gg4)}
\end{Highlighting}
\end{Shaded}

\includegraphics{4PoissonGamma_files/figure-latex/unnamed-chunk-28-2.pdf}

\begin{verbatim}
##      df       AIC
## fit1  4 104.27633
## fit4  4  45.01688
\end{verbatim}

\section{Interactive session - second
week}\label{interactive-session---second-week}

\subsection{Problem 1: Exam 2007 (Problem 1, a bit modified) - Smoking
and lung
cancer}\label{problem-1-exam-2007-problem-1-a-bit-modified---smoking-and-lung-cancer}

(Permitted aids for the exam was ``Tabeller og formler i statistikk'',
Matematisk formelsamling (Rottmann), one A5 sheet with your own
handwritten notes, and a simple calculator.)

The dataset given in \texttt{smoking.txt} consists of four variables:

\begin{itemize}
\tightlist
\item
  \texttt{deaths}: number of lung cancer deaths over a period of six
  years {[}remark: incorrectly 1 year in exam question{]}
\item
  \texttt{population}: the number of people {[}remark: incorrectly in
  100 000 people in exam question{]}
\item
  \texttt{age}: in five-year age groups
  (\texttt{40-44,\ 45-49,\ 50-54,\ 55-59,\ 60-64,\ 65-69,\ 70-74,\ 75-79,\ 80+})
\item
  \texttt{ageLevel}: age group as numbers from 1 to 9 (1 corresponds to
  40-44, 2 to 45-59, and so on)
\item
  \texttt{smoking\ status}: doesn't smoke (\texttt{no}), smokes cigars
  or pipe only (\texttt{cigarPipeOnly}), smokes cigarettes and cigar or
  pipe (\texttt{cigarettePlus}), and smokes cigarettes only
  (\texttt{cigaretteOnly})
\end{itemize}

You can look at the dataset here:
\url{https://www.math.ntnu.no/emner/TMA4315/2018h/smoking.txt}, and can
be found here as well:
\url{http://data.princeton.edu/wws509/datasets/\#smoking}

Remark: the data set is probably taken from
\url{https://www.jstor.org/stable/41983444?seq=1\#page_scan_tab_contents}
and there it is said the the persons under study are males who have
contributed in wars before 1956 and who answered a questionary. The
authors point out that the dataset is not representative for the whole
population.

We are interested in studying if the mortality rate due to lung cancer
(the number of deaths due to lung cancer per individual during six year)
controlled for age group and smoking status. Assume that the number of
deaths for each set of covariate values, \(Y_i\), can be considered
Poisson distributed, \(Y_i \sim \text{Poisson}(\lambda_i)\). We fit the
following model:

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{# load data:}
\NormalTok{smoking <-}\StringTok{ }\KeywordTok{read.table}\NormalTok{(}\DataTypeTok{file =} \StringTok{"https://www.math.ntnu.no/emner/TMA4315/2018h/smoking.txt"}\NormalTok{)}
\KeywordTok{head}\NormalTok{(smoking)}
\KeywordTok{nrow}\NormalTok{(smoking)}
\NormalTok{model1 <-}\StringTok{ }\KeywordTok{glm}\NormalTok{(dead }\OperatorTok{~}\StringTok{ }\NormalTok{age }\OperatorTok{+}\StringTok{ }\NormalTok{smoke, }\DataTypeTok{family =} \StringTok{"poisson"}\NormalTok{, }\DataTypeTok{data =}\NormalTok{ smoking, }\DataTypeTok{offset =} \KeywordTok{log}\NormalTok{(pop))}
\CommentTok{# note that the size of the population for each combination of the}
\CommentTok{# covaraites is the offset here}
\KeywordTok{summary}\NormalTok{(model1)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##   dead  pop   age ageLevel smoke
## 1   18  656 40-44        1    no
## 2   22  359 45-59        2    no
## 3   19  249 50-54        3    no
## 4   55  632 55-59        4    no
## 5  117 1067 60-64        5    no
## 6  170  897 65-69        6    no
## [1] 36
## 
## Call:
## glm(formula = dead ~ age + smoke, family = "poisson", data = smoking, 
##     offset = log(pop))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.06055  -0.54773   0.06431   0.29963   1.48348  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -3.63222    0.06783 -53.552  < 2e-16 ***
## age45-59             0.55388    0.07999   6.924 4.38e-12 ***
## age50-54             0.98039    0.07682  12.762  < 2e-16 ***
## age55-59             1.37946    0.06526  21.138  < 2e-16 ***
## age60-64             1.65423    0.06257  26.439  < 2e-16 ***
## age65-69             1.99817    0.06279  31.824  < 2e-16 ***
## age70-74             2.27141    0.06435  35.296  < 2e-16 ***
## age75-79             2.55858    0.06778  37.746  < 2e-16 ***
## age80+               2.84692    0.07242  39.310  < 2e-16 ***
## smokecigarretteOnly  0.36915    0.03791   9.737  < 2e-16 ***
## smokecigarrettePlus  0.17015    0.03643   4.671 3.00e-06 ***
## smokeno             -0.04781    0.04699  -1.017    0.309    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4055.984  on 35  degrees of freedom
## Residual deviance:   21.487  on 24  degrees of freedom
## AIC: 285.51
## 
## Number of Fisher Scoring iterations: 4
\end{verbatim}

\textbf{a)} What is an offset?\\
For 53-year old non-smokers, what is the estimated number of deaths due
to lung cancer (per person over 6 years)?\\
Why is the number of degrees of freedom for the deviance of this model
24? Does the model give a good fit?

\textbf{b)} Let \(\lambda(a,s)\) denote the expected number of lung
cancer deaths per person in age group \(a\) with smoking status \(s\).
For two different smoking statuses \(s_1\) and \(s_2\), define

\[r(a, s_1, s_2) = \frac{\lambda(a, s_1)}{\lambda(a, s_2)}.\]

Explain why \(r(a, s_1, s_2)\) does \textbf{not} vary as a function of
\(a\) in \texttt{model1}.\\
For \(s_1 =\) \texttt{cigarPipeOnly} and \(s_2 =\)
\texttt{cigaretteOnly}, find an estimate value for \(r(a, s_1, s_2)\)
and an approximate 90 \% confidence interval. Is there a significant
difference in the expected number of lung cancer deaths for individuals
that smoke cigarettes \texttt{cigaretteOnly} versus those that smoke
cigar/pipe \texttt{cigarPipeOnly}?

\textbf{c)} We will now consider two alternative models, \texttt{model2}
and \texttt{model3}:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{model2 <-}\StringTok{ }\KeywordTok{glm}\NormalTok{(dead }\OperatorTok{~}\StringTok{ }\NormalTok{smoke, }\DataTypeTok{family =} \StringTok{"poisson"}\NormalTok{, }\DataTypeTok{data =}\NormalTok{ smoking, }\DataTypeTok{offset =} \KeywordTok{log}\NormalTok{(pop))}
\NormalTok{model3 <-}\StringTok{ }\KeywordTok{glm}\NormalTok{(dead }\OperatorTok{~}\StringTok{ }\NormalTok{ageLevel }\OperatorTok{+}\StringTok{ }\NormalTok{smoke, }\DataTypeTok{family =} \StringTok{"poisson"}\NormalTok{, }\DataTypeTok{data =}\NormalTok{ smoking, }\DataTypeTok{offset =} \KeywordTok{log}\NormalTok{(pop))}

\KeywordTok{summary}\NormalTok{(model2)}
\KeywordTok{summary}\NormalTok{(model3)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## glm(formula = dead ~ smoke, family = "poisson", data = smoking, 
##     offset = log(pop))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -24.546   -5.892   -2.310    8.343   15.612  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -1.47319    0.03099 -47.532  < 2e-16 ***
## smokecigarretteOnly -0.31219    0.03576  -8.729  < 2e-16 ***
## smokecigarrettePlus -0.43013    0.03468 -12.402  < 2e-16 ***
## smokeno             -0.36678    0.04669  -7.855 3.98e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4056.0  on 35  degrees of freedom
## Residual deviance: 3910.7  on 32  degrees of freedom
## AIC: 4158.7
## 
## Number of Fisher Scoring iterations: 5
## 
## 
## Call:
## glm(formula = dead ~ ageLevel + smoke, family = "poisson", data = smoking, 
##     offset = log(pop))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.0918  -1.1673  -0.2755   0.7803   2.6364  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -3.705950   0.050717 -73.071  < 2e-16 ***
## ageLevel             0.333006   0.005591  59.559  < 2e-16 ***
## smokecigarretteOnly  0.405019   0.037463  10.811  < 2e-16 ***
## smokecigarrettePlus  0.203426   0.035996   5.651 1.59e-08 ***
## smokeno             -0.032927   0.046894  -0.702    0.483    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4055.984  on 35  degrees of freedom
## Residual deviance:   75.734  on 31  degrees of freedom
## AIC: 325.76
## 
## Number of Fisher Scoring iterations: 4
\end{verbatim}

Why does \texttt{model2} and \texttt{model3} have 32 and 31 degrees of
freedom, respectively? If we want to compare the three models
\texttt{model1}, \texttt{model2} and \texttt{model3}, which model would
you choose as the best? Justify your answer by formulating relevant
hypotheses and perform hypothesis tests, based on the print-outs above.
Hint: remember that to compare two models then one model needs to be
nested within the other model.

New: For \texttt{model3}: Plot the regression line for the expected
number of lung cancer deaths per person as a function of \texttt{age}
for each of the four different smoking levels in the same plot. First
use all coefficients, and then only the ones that are significant. Use
\texttt{ggplot}! Discuss what you see.

\subsection{Problem 2: TMA4315 Exam 2012, Problem 3: Precipitation in
Trondheim,
amount}\label{problem-2-tma4315-exam-2012-problem-3-precipitation-in-trondheim-amount}

Remark: the text is slightly modified from the original exam since we
parameterized the gamma as in our textbook.

Remark: Problems 1 (binomial) and 2 (multinomial) at the 2012 exam also
asked about percipitation.

We want to model the amount of daily precipitation given that it
\emph{is} precipitation, and denote this quantity \(Y\). It is common to
model \(Y\) as a gamma distributed random variable,
\(Y \sim Gamma(\nu,\mu)\), with density

\[ f_Y(y) = \frac{1}{\Gamma(\nu)} \left(\frac{\nu}{\mu}\right)^{\nu} y^{\nu-1}\exp\left(-\frac{\nu}{\mu} y \right) \]

In this problem we consider \(N\) observations, each gamma distributed
with \(Y_i \sim Gamma(\nu, \mu_i)\) (remark: common \(\nu\)). Here
\(\nu\) is considered to be a known nuisance parameter, and the
\(\mu_i\)s are unknown.

\textbf{a)} Show that the gamma distribution function is member of the
exponential family when \(\mu_i\) is the parameter of interest.\\
Use this to find expressions for the expected value and the variance of
\(Y_i\), in terms of \((\nu,\mu_i)\), and interpret \(\nu\).\\
New: What is the canonical link for the Gamma distribution?

Hint: if you want to focus on discussing - you may look at the solutions
from
\href{https://www.math.ntnu.no/emner/TMA4315/2017h/Module1ExponentialFamily.pdf}{Module
1 together}.

\textbf{b)} Explain what a saturated model is.\\
Set up the log-likelihood function expressed by \(\mu_i\), and use it to
find the maximum likelihood estimators for \(\mu_i\)-s of the saturated
model.\\
Find the deviance (based on all \(N\) observations).

Hint: Do you see directly that \(\hat{\mu}_i=y_i\)? If not, you may look
at the likelihood for one observation and solve that the derivative
equals 0.

\textbf{c)} We now want to construct a model for amount of precipitation
(given that there are occurrence) with precipitation forecast as
explanatory variable.\\
Let \(Y_i\) be amount of precipitation for day \(i\), and let \(x_i\) be
the precipitation forecast valid for day \(i\). Set up a GLM for this,
and argue for your choice of link function and linear predictor.

Hint: explain why this fits into our GLM-gamma framework and set up our
3 equation model (random, systematic, link).

\begin{center}\rule{0.5\linewidth}{\linethickness}\end{center}

\subsection{Problem 3: Taken from UiO, STK3100, 2015, problem
2}\label{problem-3-taken-from-uio-stk3100-2015-problem-2}

(For reference:
\href{http://www.uio.no/studier/emner/matnat/math/STK3100/oppgaver/stk3100-f15eng_final.pdf}{here
is the original exam}).

This is a problem on the logistic regression.

Do not look at the dataset before the end of the exercise! You should
solve the exercise without using \texttt{R}, just as if you were at an
exam.

In this problem you shall consider data of survivals from a study of
treatment for breast cancer. The response is the number that survived
for three years. The covariates were the four factors

\begin{itemize}
\tightlist
\item
  \texttt{app}: appearance of tumor, two levels, 1 = malignant, 2 =
  benign
\item
  \texttt{infl}: inflammatory reaction, two levels, 1 = minimal, 2 =
  moderate or severe
\item
  \texttt{age}: age of patients, three levels, 1 = under 50, 2 = 50 to
  69, 3 = 70 or older
\item
  \texttt{country}: hospital of treatment, three levels, 1 = Japan, 2 =
  US, 3 = UK
\end{itemize}

The dataset we have used differs slightly from the one they used at UiO
(the number of survivals and non-survivals differ).

The number of survivors is modelled as a binomially distributed variable
using a canonical logit link. Dummy variable coding is used for all
factors, with the level coded as ``1'' as the reference category.

\textbf{a)} The output from fitting the model where only appearance and
country are used as covariates, i.e., a model with predictor on the form

\[ \eta = \beta_0 + \beta_1 \texttt{ fapp} + \beta_2 \texttt{ fcountry2} + \beta_3 \texttt{ fcountry3} \]

is displayed below. What is the interpretation of the estimate of the
coefficient of appearance, \texttt{fapp} (\texttt{f} means factor)?
Explain also how this coefficient can be expressed in terms of an odds
ratio.

Hint: first explain what is the predicted odds of survivial for country
\(j\) for a benign tumor, then the same for a malignant tumor- and then
make an odds ratio.

New: The main difference between our dataset and UiO's dataset is the
number of degrees of freedom for the null model; we have 34, they have
35. How many observations do we have in the dataset? And how many does
UiO have? Hint: All the covariates in the dataset are categorical, two
with 3 levels, and two with 2 levels. How many possible combinations of
observations does that make?

\begin{verbatim}
Call:
glm(formula = cbind(surv, nsurv) ~ fapp + fcountry, family = binomial, 
    data = brc)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8142  -0.7279   0.2147   0.7576   1.8715  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   1.0803     0.1656   6.522 6.93e-11 ***
fapp2         0.5157     0.1662   3.103 0.001913 ** 
fcountry2    -0.6589     0.1998  -3.298 0.000972 ***
fcountry3    -0.4945     0.2071  -2.388 0.016946 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 57.588  on 34  degrees of freedom
Residual deviance: 36.625  on 31  degrees of freedom
\end{verbatim}

\textbf{b)} The output below is an analysis of deviance for comparing
various model spesifications. Fill out the positions indicated
\texttt{i-iv}.\\
Note: \texttt{y\ \textasciitilde{}\ x1*x2} gives both the linear
components and the interaction, i.e., it is the same as
\texttt{y\ \textasciitilde{}\ x1\ +\ x2\ +\ x1:x2} (\texttt{x1:x2} gives
interaction only).

Remark: remember that \texttt{anova} gives the sequential comparison of
deviances - that is - comparing each model to the previous model.

\begin{verbatim}
Analysis of Deviance Table

Model 1: cbind(surv, nsurv) ~ fapp + fage + fcountry
Model 2: cbind(surv, nsurv) ~ fapp + fage + finfl + fcountry
Model 3: cbind(surv, nsurv) ~ fapp + finfl + fage * fcountry
Model 4: cbind(surv, nsurv) ~ fapp * finfl + fage * fcountry
Model 5: cbind(surv, nsurv) ~ fapp * finfl + fapp * fage + fage * fcountry
Model 6: cbind(surv, nsurv) ~ fapp * finfl * fage * fcountry
  Resid. Df Resid. Dev Df Deviance
1        29     33.102            
2         i     33.095  1   0.0065
3        24     25.674 ii   7.4210
4        23     25.504  1      iii
5        21     22.021  2   3.4823
6         0      0.000 iv  22.0214
\end{verbatim}

In the remaining parts of this problem we return to the model in part a)
and consider the hypothesis

\[ H_0 : \beta_2 + \beta_3 = -1 \text{ versus } H_1 : \beta_2 + \beta_3 \neq -1 \]

Note: what are you testing now?

\textbf{c)} The estimated covariance matrix between the estimators of
the coefficients \(\beta_2\) and \(\beta_3\) above is
\(\left(\begin{matrix} 0.040 & 0.021 \\ 0.021 & 0.043\end{matrix} \right)\).
Use a Wald test to test the null hypothesis above.

Note: remember - this was a written exam so you do this by hand!

\textbf{d)} Explain how the null hypothesis can be tested with a
likelihood ratio test by fitting two suitable models. No numerical
calculations are necessary, but it must be specified how the predictors
should be defined.\\
Remark: rather technical - and hint: offset term?

New: Do this test in \texttt{R}.
\href{https://www.math.ntnu.no/emner/TMA4315/2018h/breastcancer.txt}{Here
is the data set!}

\section{Work on your own: Exam
questions}\label{work-on-your-own-exam-questions}

\subsection{December 2013 (Essay exam)}\label{december-2013-essay-exam}

We will consider the following Poisson regression
\[Y_i \sim \text{Poisson}(\exp(\eta_i)), \text{ }i=1,\ldots,n\] where
the linear predictor is \(\eta_i={\bf x}_i^T\beta\). Here \({\bf x}_i\)
is a vector of the \(p\) covariates for the \(i\)th observation \(Y_i\)
and \(\beta\) is unknown \(p\) dimensional column vector of unknown
regression coefficients.

Write an introduction to Poisson regression and its practical usage, for
a student with a good background in statistics, but no knowledge about
Generalized Linear Models (GLM). Topics you may want to consider, are

\begin{itemize}
\tightlist
\item
  When to use it? Underlying assumptions.
\item
  Parameter estimation, limiting results for the MLE, Fisher information
  and observed Fisher information, confidence intervals and hypothesis
  testing.
\item
  Output analysis, residual plots and interpretation of results.
\item
  Deviance and its usage.
\item
  What do we do when a covariate is a factor, and should the results be
  interpreted?
\item
  The use of Poisson regression in the analysis of contingency tables.
\end{itemize}

\section{R packages}\label{r-packages}

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{install.packages}\NormalTok{(}\KeywordTok{c}\NormalTok{(}\StringTok{"tidyverse"}\NormalTok{, }\StringTok{"ggplot2"}\NormalTok{, }\StringTok{"statmod"}\NormalTok{, }\StringTok{"corrplot"}\NormalTok{, }\StringTok{"ggplot2"}\NormalTok{, }
    \StringTok{"GGally"}\NormalTok{, }\StringTok{"boot"}\NormalTok{))}
\end{Highlighting}
\end{Shaded}

\section{Further reading}\label{further-reading}

\begin{itemize}
\tightlist
\item
  A. Agresti (1996): ``An Introduction to Categorical Data Analysis''.
\item
  A. Agresti (2015): ``Foundations of Linear and Generalized Linear
  Models.'' Wiley.
\item
  A. J. Dobson and A. G. Barnett (2008): ``An Introduction to
  Generalized Linear Models'', Third edition.
\item
  J. Faraway (2015): ``Extending the Linear Model with R'', Second
  Edition. \url{http://www.maths.bath.ac.uk/~jjf23/ELM/}
\item
  P. McCullagh and J. A. Nelder (1989): ``Generalized Linear Models''.
  Second edition.
\end{itemize}


\end{document}
