\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 and Ingeborg Hem, Department of Mathematical Sciences, NTNU},
            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{longtable,booktabs}
\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 8: Generalized linear mixed effects models (GLMM)}
  \author{Mette Langaas and Ingeborg Hem, Department of Mathematical Sciences,
NTNU}
    \preauthor{\centering\large\emph}
  \postauthor{\par}
      \predate{\centering\large\emph}
  \postdate{\par}
    \date{15.11.2018 {[}PL{]}, 16.11.2018 {[}IL{]}}


\begin{document}
\maketitle

{
\setcounter{tocdepth}{2}
\tableofcontents
}
(Latest changes: 16.11: IL added links, 13.11.2018 - first version).

\section{Overview}\label{overview}

Aim: Present methods for analysing correlated responses in a generalized
linear models setting - LMM meets GLM to become \emph{GLMM}.

Also here (as in module LMM) we will only consider \emph{two-level
models} and in particular focus on \emph{random intercept models} for
binomial and Poisson responses. Emphasis will be on understanding.

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

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

\begin{itemize}
\tightlist
\item
  Textbook: Fahrmeir et al (2013): Chapter 7.5: 389-394, 7.7.
\item
  \href{https://www.math.ntnu.no/emner/TMA4315/2018h/M8PL.pdf}{Classnotes
  (15.11.2018)}
\end{itemize}

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

\subsection{Topics}\label{topics}

\begin{itemize}
\tightlist
\item
  beaches example - revisited
\item
  notation
\item
  the generalized linear mixed effect model (three ingredients)
\item
  the GLMM with random intercept
\item
  the marginal model
\item
  parameter estimation and Laplace approximation
\item
  summing up: what do we need to know about the GLMM?
\item
  \protect\hyperlink{additional}{additional info on different software
  (not on reading list)}
\end{itemize}

Jump to \protect\hyperlink{interactive}{interactive}.

\section{Beaches example - revisited}\label{beaches-example---revisited}

This example is taken from Zuur et al. (2009, chapter 5, pages 101-142),
and data are referred to as RIKZ. Data were collected on nine different
beaches in the Netherlands with the aim to investigate if there is a
relationship between the

\begin{itemize}
\tightlist
\item
  richness of species (number of species observed) and
\item
  NAP: the height of the sampling station compared to mean tidal level.
\end{itemize}

Data: 45 observations, taken at 9 beaches:

\begin{itemize}
\tightlist
\item
  beach: the beach that the samples were taken, for each beach 5
  different samples were taken.
\end{itemize}

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

\includegraphics{8GLMM_files/figure-latex/unnamed-chunk-1-1.pdf}

\small

We want to use the richness of species as the response in a regression
model. In Module LMM we assumed that the response could be viewed as a
normal variable, but now we rather assume that the response comes from a
\emph{Poisson distribution}.

\textbf{Q}: Why do we want to assume a Poisson distribution? \normalsize

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

\textbf{A}:

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}

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

\begin{itemize}
\tightlist
\item
  The number of events in a Poisson process follows a Poisson
  distribution.
\end{itemize}

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

We go ahead and fit a Poission GLM with log-link and NAP as covariate.
In Module LMM we talked about possible solutions to how to handle
multiple observations for each beach and came up with:

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

\textbf{Possible solutions:}

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  Use all 45 observations in a Poisson regression - with a common
  intercept and linear effect in NAP. \textbf{Problem:} violation of
  assumption of independent observations lead to wrong estimates for
  variances of parameter estimates.
\item
  Add beach as covariate to Poisson regression - then we estimate one
  regression coefficient for each beach (intercepts) in addition to the
  linear effect in NAP. \textbf{Problem:} why do we want to add 8 extra
  parameters to estimate (why 8?) values for the 9 beaches? Loss of
  power and what would we use the beach estimates for?
\item
  Add beach as a \emph{random} covariate to the Poisson regression: this
  is called \emph{random intercept models}. \textbf{Problem:} new stuff
  - slightly complicated. We do this because beaches not of interest in
  themselves, only random sample from population of beaches, and
  therefore we only need to account for beaches, not estimate separate
  parameters.
\end{enumerate}

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

But - first: what do we remember about GLMs (and Poission regression
with log-link in particular)?

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

\textbf{Assumptions for GLM model (Poisson): }

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  Random component from exponential family: We have for our beaches:
  \(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 {\boldsymbol \beta}\).
\item
  Log link \[\eta_i=\ln(\lambda_i)=g(\lambda_i)\] and (inverse thereof)
  response function \[\lambda_i=\exp(\eta_i)\] Canonical link for
  Poisson is log.
\end{enumerate}

Parameter estimation is performed using the method of maximum
likelihood, numerically using Fisher scoring.

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

\textbf{Solution 1:} all observations in one GLM together - standard
errors not correct since assumption of independent obsevations is
violated.

Remember: For interpretation. 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 \emph{factor} \(\exp(\beta_1)\). Not so easy to compare to the normal
case. \footnotesize

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{fitall =}\StringTok{ }\KeywordTok{glm}\NormalTok{(Richness }\OperatorTok{~}\StringTok{ }\NormalTok{NAP, }\DataTypeTok{data =}\NormalTok{ RIKZ, }\DataTypeTok{family =} \KeywordTok{poisson}\NormalTok{(}\DataTypeTok{link =}\NormalTok{ log))}
\KeywordTok{summary}\NormalTok{(fitall)}
\NormalTok{fitall}\OperatorTok{$}\NormalTok{coefficients[}\DecValTok{2}\NormalTok{]}
\KeywordTok{exp}\NormalTok{(fitall}\OperatorTok{$}\NormalTok{coefficients[}\DecValTok{2}\NormalTok{])}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## glm(formula = Richness ~ NAP, family = poisson(link = log), data = RIKZ)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2029  -1.2432  -0.9199   0.3943   4.3256  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.79100    0.06329  28.297  < 2e-16 ***
## NAP         -0.55597    0.07163  -7.762 8.39e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 179.75  on 44  degrees of freedom
## Residual deviance: 113.18  on 43  degrees of freedom
## AIC: 259.18
## 
## Number of Fisher Scoring iterations: 5
## 
##        NAP 
## -0.5559735 
##       NAP 
## 0.5735136
\end{verbatim}

\normalsize

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

\small

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(ggplot2)}
\KeywordTok{library}\NormalTok{(ggpubr)}
\NormalTok{df=}\KeywordTok{data.frame}\NormalTok{(}\StringTok{"Y"}\NormalTok{=RIKZ}\OperatorTok{$}\NormalTok{Richness,}\StringTok{"x"}\NormalTok{=RIKZ}\OperatorTok{$}\NormalTok{NAP,}\StringTok{"fitted"}\NormalTok{=fitall}\OperatorTok{$}\NormalTok{fitted.values,}\StringTok{"dres"}\NormalTok{=}\KeywordTok{residuals}\NormalTok{(fitall,}\DataTypeTok{type=}\StringTok{"deviance"}\NormalTok{),}\StringTok{"pres"}\NormalTok{=}\KeywordTok{residuals}\NormalTok{(fitall,}\DataTypeTok{type=}\StringTok{"pearson"}\NormalTok{),}\DataTypeTok{catBeach=}\KeywordTok{as.factor}\NormalTok{(RIKZ}\OperatorTok{$}\NormalTok{Beach))}
\NormalTok{gg1=}\KeywordTok{ggplot}\NormalTok{(df)}\OperatorTok{+}\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{x,}\DataTypeTok{y=}\NormalTok{Y,}\DataTypeTok{colour=}\NormalTok{catBeach))}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{x,}\DataTypeTok{y=}\NormalTok{fitted))}\OperatorTok{+}\KeywordTok{xlab}\NormalTok{(}\StringTok{"NAP"}\NormalTok{)}\OperatorTok{+}\KeywordTok{ylab}\NormalTok{(}\StringTok{"Richness"}\NormalTok{)}
\NormalTok{gg2=}\KeywordTok{ggplot}\NormalTok{(df)}\OperatorTok{+}\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{fitted,}\DataTypeTok{y=}\NormalTok{dres,}\DataTypeTok{color=}\NormalTok{catBeach))}\OperatorTok{+}\KeywordTok{xlab}\NormalTok{(}\StringTok{"fitted"}\NormalTok{)}\OperatorTok{+}\KeywordTok{ylab}\NormalTok{(}\StringTok{"deviance residuals"}\NormalTok{)}
\KeywordTok{ggarrange}\NormalTok{(gg1,gg2)  }
\end{Highlighting}
\end{Shaded}

\includegraphics{8GLMM_files/figure-latex/unnamed-chunk-3-1.pdf}

\normalsize
**Q\textbf{: What is the formula for the black line on the left panel?
}A**: see classnotes.

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

\textbf{Solution 2:} fixed effects for each beach - many estimates not
so much of interest when population is in focus. \footnotesize

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{RIKZ}\OperatorTok{$}\NormalTok{beachfactor =}\StringTok{ }\KeywordTok{as.factor}\NormalTok{(RIKZ}\OperatorTok{$}\NormalTok{Beach)}
\NormalTok{fitbeach =}\StringTok{ }\KeywordTok{glm}\NormalTok{(Richness }\OperatorTok{~}\StringTok{ }\NormalTok{NAP }\OperatorTok{+}\StringTok{ }\NormalTok{beachfactor, }\DataTypeTok{data =}\NormalTok{ RIKZ, }\DataTypeTok{contrasts =} \KeywordTok{list}\NormalTok{(}\DataTypeTok{beachfactor =} \StringTok{"contr.sum"}\NormalTok{), }
    \DataTypeTok{family =} \KeywordTok{poisson}\NormalTok{(}\DataTypeTok{link =} \StringTok{"log"}\NormalTok{))}
\KeywordTok{summary}\NormalTok{(fitbeach)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## glm(formula = Richness ~ NAP + beachfactor, family = poisson(link = "log"), 
##     data = RIKZ, contrasts = list(beachfactor = "contr.sum"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.8849  -0.5240  -0.1194   0.3805   2.6037  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.64239    0.07356  22.326  < 2e-16 ***
## NAP          -0.49467    0.07642  -6.473 9.61e-11 ***
## beachfactor1  0.45928    0.14943   3.074 0.002115 ** 
## beachfactor2  0.88808    0.13527   6.565 5.20e-11 ***
## beachfactor3 -0.49422    0.22710  -2.176 0.029540 *  
## beachfactor4 -0.59982    0.26525  -2.261 0.023740 *  
## beachfactor5  0.58824    0.16469   3.572 0.000355 ***
## beachfactor6 -0.27433    0.21124  -1.299 0.194068    
## beachfactor7 -0.47708    0.28070  -1.700 0.089202 .  
## beachfactor8 -0.16553    0.21057  -0.786 0.431786    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 179.753  on 44  degrees of freedom
## Residual deviance:  47.073  on 35  degrees of freedom
## AIC: 209.07
## 
## Number of Fisher Scoring iterations: 5
\end{verbatim}

\normalsize

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

\small

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(ggplot2)}
\KeywordTok{library}\NormalTok{(ggpubr)}
\NormalTok{df=}\KeywordTok{data.frame}\NormalTok{(}\StringTok{"Y"}\NormalTok{=RIKZ}\OperatorTok{$}\NormalTok{Richness,}\StringTok{"x"}\NormalTok{=RIKZ}\OperatorTok{$}\NormalTok{NAP,}\StringTok{"fitted"}\NormalTok{=fitbeach}\OperatorTok{$}\NormalTok{fitted.values,}\StringTok{"dres"}\NormalTok{=}\KeywordTok{residuals}\NormalTok{(fitbeach,}\DataTypeTok{type=}\StringTok{"deviance"}\NormalTok{),}\StringTok{"pres"}\NormalTok{=}\KeywordTok{residuals}\NormalTok{(fitbeach,}\DataTypeTok{type=}\StringTok{"pearson"}\NormalTok{),}\DataTypeTok{catBeach=}\KeywordTok{as.factor}\NormalTok{(RIKZ}\OperatorTok{$}\NormalTok{Beach))}
\NormalTok{gg1=}\KeywordTok{ggplot}\NormalTok{(df)}\OperatorTok{+}\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{x,}\DataTypeTok{y=}\NormalTok{Y,}\DataTypeTok{colour=}\NormalTok{catBeach))}\OperatorTok{+}\KeywordTok{geom_line}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{x,}\DataTypeTok{y=}\NormalTok{fitted))}\OperatorTok{+}\KeywordTok{xlab}\NormalTok{(}\StringTok{"NAP"}\NormalTok{)}\OperatorTok{+}\KeywordTok{ylab}\NormalTok{(}\StringTok{"Richness"}\NormalTok{)}\OperatorTok{+}\KeywordTok{facet_wrap}\NormalTok{(}\OperatorTok{~}\NormalTok{catBeach)}
\NormalTok{gg2=}\KeywordTok{ggplot}\NormalTok{(df)}\OperatorTok{+}\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x=}\NormalTok{fitted,}\DataTypeTok{y=}\NormalTok{dres,}\DataTypeTok{color=}\NormalTok{catBeach))}\OperatorTok{+}\KeywordTok{xlab}\NormalTok{(}\StringTok{"fitted"}\NormalTok{)}\OperatorTok{+}\KeywordTok{ylab}\NormalTok{(}\StringTok{"deviance residuals"}\NormalTok{)}
\KeywordTok{ggarrange}\NormalTok{(gg1,gg2)  }
\end{Highlighting}
\end{Shaded}

\includegraphics{8GLMM_files/figure-latex/unnamed-chunk-5-1.pdf}
\normalsize

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

\section{Generalized linear mixed effects
models}\label{generalized-linear-mixed-effects-models}

\subsection{Notation}\label{notation}

combining GLM (Module 5) and LMM (Module 7)

We have \(n_i\) repeated observations (\(j=1,\ldots, n_i\)) from each of
\(i=1,\ldots,m\) clusters (or individuals).

\begin{itemize}
\tightlist
\item
  Responses: \(Y_{i1},Y_{i2},\ldots, Y_{in_i}\) (e.g.~species richness
  at beach \(i\) sample \(j\))
\item
  Covariates: \({\bf x}_{i1},{\bf x}_{i1}, \ldots, {\bf x}_{in_i}\)
  (e.g.~NAP for beach \(i\) sample \(j\))
\end{itemize}

The covariates \({\bf x}_{ij}\) are \(p\times 1\) vectors (as before,
\(k\) covariates and one intercept so \(p=k+1\)).

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

\subsection{Repetition GLM}\label{repetition-glm}

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
  \textbf{Random component} - exponential family (now \(y_{ij}\))
\end{enumerate}

\[ f(y_{ij}\mid \theta_{ij})=\exp \left( \frac{y_{ij} \theta_{ij}-b(\theta_{ij})}{\phi}\cdot w_{ij} + c(y_{ij}, \phi, w_{ij}) \right) \]

\begin{itemize}
\item
  \(\theta_{ij}\) 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_{ij}\) is a weight function, in most cases \(w_{ij}=1\)
\item
  \(b\) and \(c\) are known functions.
\end{itemize}

We looked in detail into the binomial, Poisson and gamma distribution
(in addition to the normal).

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

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\setcounter{enumi}{1}
\tightlist
\item
  \textbf{Linear predictor} (previously ``systematic component''):
\end{enumerate}

\[\eta_{ij}={\bf x}_{ij}^T {\boldsymbol \beta}\]

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\setcounter{enumi}{2}
\tightlist
\item
  \textbf{Link function \(g\)- and response function \(h\)}
\end{enumerate}

\[\eta_{ij}=g(\mu_{ij})\]

\[\mu_{ij}=h(\eta_{ij})\]

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

\subsection{New element}\label{new-element}

we add a random component (effect) to the \emph{linear predictor}!

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\setcounter{enumi}{1}
\tightlist
\item
  \[\eta_{ij}={\bf x}_{ij}^T {\boldsymbol \beta}+ {\bf u}^T_{ij} \gamma_i\]
  and we make the same assumptions as for the LMM
\end{enumerate}

\[ \gamma_i \sim N({\bf 0},{\bf Q})\] and independent from covariates.
Observe: we still assume normality for the random effects!

In addition, the mean \(\mu_{ij}=\text{E}(Y_{ij}\mid \gamma_i)\) is now
to be seen to be conditional on the random effect \(\gamma_i\). This
means that the specified response distribution is conditional on the
random effect: \(f(y_{ij}\mid \gamma_i)\).

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

\subsection{GLMM model assumptions}\label{glmm-model-assumptions}

\subsubsection{Distributional
assumptions}\label{distributional-assumptions}

(previously ``1. Random component from exponential family'')

\begin{itemize}
\tightlist
\item
  Given the random effects \(\gamma_i\) and the covariates
  \({\bf x}_{ij}\) and \({\bf u}_{ij}\)
\item
  the responses \(Y_{ij}\) are \emph{conditionally independent} and
\item
  the \emph{conditional distribution} \(f(y_{ij}\mid \gamma_i)\) belongs
  to and exponential family.
\end{itemize}

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

\subsubsection{Structural assumptions}\label{structural-assumptions}

(previously ``2. Systematic component'' and ``3. Link function'')

\begin{itemize}
\tightlist
\item
  The \emph{conditional mean} \(\mu_{ij}=\text{E}(Y_{ij}\mid \gamma_i)\)
  is linked to
\item
  the \emph{linear predictor}
  \(\eta_{ij}={\bf x}_{ij}^T {\boldsymbol \beta}+ {\bf u}_{ij} \gamma_i\)
\item
  through the link function \[\eta_{ij}=g(\mu_{ij})\]
\item
  or equivalently through the response function
  \[\mu_{ij}=h(\eta_{ij})\]
\item
  where \(h^{-1}=g\) (inverse functions).
\end{itemize}

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

\subsubsection{Distributional assumptions for random
effects}\label{distributional-assumptions-for-random-effects}

\begin{itemize}
\tightlist
\item
  The random effects \(\gamma_i\), \(i=1,\ldots,m\) are independent and
  identically distributed \[ \gamma_i \sim N({\bf 0}, {\bf Q})\]
\item
  where the covariance matrix \({\bf Q}\) is a \((q+1)\times (q+1)\)
  positive definite.
\item
  An important special case is
  \({\bf Q}=diag(\tau_0^2, \tau_1^2,\ldots,\tau_q^2)\).
\end{itemize}

\textbf{Q}: Any questions? Notice that we still have normal distribution
for the random effects!

\textbf{Remark}: for LMMs we used that
\(\varepsilon \sim N({\bf 0}, \sigma^2 {\bf I})\), but said that it was
possible to let the error terms be correlated
\(\varepsilon \sim N({\bf 0}, {\bf R})\) where \({\bf R}\) could har
non-zero off-diagonal elements. However, for GLMMs this is much more
complicated.

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

\subsection{Alternative two-step
formulation}\label{alternative-two-step-formulation}

Consider a random variable \(Y_{ij}\) and a random effect \(\gamma_i\)
for cluster \(i\).

Alternative two step formulation:

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
  \(\gamma_i \sim N({\bf 0},{\bf Q})\)
\item
  Given the realization of \(\gamma_i\), then the \(Y_{ij}\)s are
  independent and each have the distribution \(f(y_{ij}\mid \gamma_i)\)
  with mean \(\mu_{ij}=h(\eta_{ij})\) and linear predictor
  \(\eta_{ij}={\bf x}_{ij}^T {\boldsymbol \beta}+ {\bf u}^T_{ij} \gamma_i\).
\end{enumerate}

Remark: this conditional formulation is ``easy'' to write down, but to
do inference we need the marginal distribution of all \(Y_{ij}\)s
together to construct the likelihood. We will soon see that this is
hard.

\textbf{Q:} How would you proceed to generate a GLMM data set -
e.g.~based on the beaches example?

\textbf{A:} see class notes

\section{GLMM random intercept model}\label{glmm-random-intercept-model}

\textbf{Distributional assumptions:} we will focus on binomial and
Poisson responses.

\textbf{Structural assumptions:} The random effect is added to the
linear predictor

\[\eta_{ij}={\bf x}_{ij}^T {\boldsymbol \beta}+ \gamma_{0i}\]

and we only consider canonical links, so logit for binomial and log for
Poisson.

\textbf{Distributional assumptions for random effects}

The \(\gamma_{0i}\), \(i=1,\ldots,m\) are independent and identically
distributed \[\gamma_{0i}\sim N(0,\tau_0^2)\]

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

\subsection{Poisson random intercept
model}\label{poisson-random-intercept-model}

\textbf{Distributional assumptions:} The conditional distribution of the
response \(Y_{ij}\) is Poisson with mean \(\lambda_{ij}\) \[
f(y_{ij}\mid \gamma_i)=\frac{\lambda_{ij}^{y_{ij}}}{y_{ij}!}\exp(-\lambda_{ij}) \text{ for } y=0,1,2,\ldots
\] This is conditional on \(\lambda_{ij}\), which means that it is
really conditional on the linear predictor - since the mean is a
function of the linear predictor, and in the linear predictor we have
fixed and random effects. So, conditional on the values for the fixed
and random effects.

Then the observations \(Y_{ij}\) are \emph{conditionally} independent
for all \(i\) and \(j\) - but, not \emph{marginally} independent.

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

\textbf{Structural assumptions:}
\[\eta_{ij}={\bf x}_{ij}^T {\boldsymbol \beta}+ \gamma_{0i}\]

\[\ln(\lambda_{ij})=\eta_{ij} \text{ or } \lambda_{ij}=\exp(\eta_{ij})\]

\textbf{Distributional assumptions for random effects}

The \(\gamma_{0i}\), \(i=1,\ldots,m\) are independent and identically
distributed \[\gamma_{0i}\sim N(0,\tau_0^2)\]

\textbf{Remark}: It is not necessary that \(n_i>1\) to use this model,
and this model can be used to take care of overdispersion (when all
\(n_i=1\)).

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

\subsection{Beach-example: Poisson GLMM with beach as random
intercept}\label{beach-example-poisson-glmm-with-beach-as-random-intercept}

(we will talk about parameter estimation soon)

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(}\StringTok{"AED"}\NormalTok{)}
\KeywordTok{data}\NormalTok{(RIKZ)}
\KeywordTok{library}\NormalTok{(lme4)}
\NormalTok{fitRI =}\StringTok{ }\KeywordTok{glmer}\NormalTok{(Richness }\OperatorTok{~}\StringTok{ }\NormalTok{NAP }\OperatorTok{+}\StringTok{ }\NormalTok{(}\DecValTok{1} \OperatorTok{|}\StringTok{ }\NormalTok{Beach), }\DataTypeTok{data =}\NormalTok{ RIKZ, }\DataTypeTok{family =} \KeywordTok{poisson}\NormalTok{(}\DataTypeTok{link =}\NormalTok{ log))}
\KeywordTok{summary}\NormalTok{(fitRI)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Richness ~ NAP + (1 | Beach)
##    Data: RIKZ
## 
##      AIC      BIC   logLik deviance df.resid 
##    220.8    226.2   -107.4    214.8       42 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9648 -0.6155 -0.2243  0.2236  3.1869 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Beach  (Intercept) 0.2249   0.4743  
## Number of obs: 45, groups:  Beach, 9
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.66233    0.17373   9.569  < 2e-16 ***
## NAP         -0.50389    0.07535  -6.687 2.28e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## NAP 0.013
\end{verbatim}

\normalsize

\textbf{Q}: Comment on what you see.

\textbf{A}: \texttt{glmer}, ``Laplace approximation'', we now have
\(p\)-values, and for the fixed effects we read z and not t. No estimate
for \(\sigma^2\) (or \texttt{residual} term in random effects part) -
since the model does not contain any dispersion parameter.

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

\tiny

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

\KeywordTok{plot_model}\NormalTok{(fitRI, }\DataTypeTok{type =} \StringTok{"re"}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{ylab}\NormalTok{(}\StringTok{"BLUP"}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{xlab}\NormalTok{(}\StringTok{"Beaches"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\includegraphics{8GLMM_files/figure-latex/unnamed-chunk-7-1.pdf}

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

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{df <-}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{Y =}\NormalTok{ RIKZ}\OperatorTok{$}\NormalTok{Richness, }\DataTypeTok{x =}\NormalTok{ RIKZ}\OperatorTok{$}\NormalTok{NAP, }\DataTypeTok{fitted =} \KeywordTok{fitted}\NormalTok{(fitRI), }\DataTypeTok{dres =} \KeywordTok{residuals}\NormalTok{(fitRI, }
    \DataTypeTok{scale =} \OtherTok{TRUE}\NormalTok{, }\DataTypeTok{type =} \StringTok{"deviance"}\NormalTok{), }\DataTypeTok{catBeach =} \KeywordTok{as.factor}\NormalTok{(RIKZ}\OperatorTok{$}\NormalTok{Beach))}
\KeywordTok{ggplot}\NormalTok{(df) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ x, }\DataTypeTok{y =}\NormalTok{ Y, }\DataTypeTok{col =}\NormalTok{ catBeach)) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_line}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ x, }
    \DataTypeTok{y =}\NormalTok{ fitted)) }\OperatorTok{+}\StringTok{ }\KeywordTok{xlab}\NormalTok{(}\StringTok{"NAP"}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{ylab}\NormalTok{(}\StringTok{"Richness"}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{facet_wrap}\NormalTok{(}\OperatorTok{~}\NormalTok{catBeach)}
\end{Highlighting}
\end{Shaded}

\includegraphics{8GLMM_files/figure-latex/unnamed-chunk-8-1.pdf}

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

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{xx <-}\StringTok{ }\KeywordTok{rep}\NormalTok{(}\KeywordTok{seq}\NormalTok{(}\KeywordTok{min}\NormalTok{(RIKZ}\OperatorTok{$}\NormalTok{NAP), }\KeywordTok{max}\NormalTok{(RIKZ}\OperatorTok{$}\NormalTok{NAP), }\DataTypeTok{length.out =} \DecValTok{100}\NormalTok{), }\DataTypeTok{each =} \DecValTok{9}\NormalTok{)}
\NormalTok{df <-}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ xx, }\DataTypeTok{eta =} \KeywordTok{coef}\NormalTok{(fitRI)}\OperatorTok{$}\NormalTok{Beach[, }\DecValTok{1}\NormalTok{] }\OperatorTok{+}\StringTok{ }\KeywordTok{coef}\NormalTok{(fitRI)}\OperatorTok{$}\NormalTok{Beach[, }
    \DecValTok{2}\NormalTok{] }\OperatorTok{*}\StringTok{ }\NormalTok{xx, }\DataTypeTok{Beach =} \KeywordTok{factor}\NormalTok{(}\KeywordTok{rep}\NormalTok{(}\DecValTok{1}\OperatorTok{:}\DecValTok{9}\NormalTok{, }\DataTypeTok{times =} \DecValTok{100}\NormalTok{)))}
\KeywordTok{ggplot}\NormalTok{(df, }\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ x, }\DataTypeTok{y =} \KeywordTok{exp}\NormalTok{(eta), }\DataTypeTok{col =}\NormalTok{ Beach)) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_line}\NormalTok{(}\DataTypeTok{size =} \DecValTok{1}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{labs}\NormalTok{(}\DataTypeTok{x =} \StringTok{"NAP"}\NormalTok{, }
    \DataTypeTok{y =} \StringTok{"Predicted incidents"}\NormalTok{, }\DataTypeTok{title =} \StringTok{"Predicted incidents of NAP on Richness"}\NormalTok{) }\OperatorTok{+}\StringTok{ }
\StringTok{    }\KeywordTok{facet_wrap}\NormalTok{(}\OperatorTok{~}\NormalTok{Beach)}
\end{Highlighting}
\end{Shaded}

\includegraphics{8GLMM_files/figure-latex/unnamed-chunk-9-1.pdf}

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

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{df <-}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{fitted =} \KeywordTok{fitted}\NormalTok{(fitRI), }\DataTypeTok{resid =} \KeywordTok{residuals}\NormalTok{(fitRI, }\DataTypeTok{scaled =} \OtherTok{TRUE}\NormalTok{, }
    \DataTypeTok{type =} \StringTok{"deviance"}\NormalTok{), }\DataTypeTok{Beach =} \KeywordTok{factor}\NormalTok{(RIKZ}\OperatorTok{$}\NormalTok{Beach))}
\KeywordTok{ggplot}\NormalTok{(df, }\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =} \KeywordTok{log}\NormalTok{(fitted), }\DataTypeTok{y =}\NormalTok{ resid)) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_hline}\NormalTok{(}\DataTypeTok{yintercept =} \DecValTok{0}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{col =}\NormalTok{ Beach)) }\OperatorTok{+}\StringTok{ }
\StringTok{    }\KeywordTok{geom_smooth}\NormalTok{(}\DataTypeTok{method =} \StringTok{"loess"}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{labs}\NormalTok{(}\DataTypeTok{x =} \StringTok{"Log-predicted values"}\NormalTok{, }\DataTypeTok{y =} \StringTok{"Deviance residuals"}\NormalTok{, }
    \DataTypeTok{title =} \StringTok{"Residuals plot"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\includegraphics{8GLMM_files/figure-latex/unnamed-chunk-10-1.pdf}

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

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{df <-}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{NAP =}\NormalTok{ RIKZ}\OperatorTok{$}\NormalTok{NAP, }\DataTypeTok{resid =} \KeywordTok{residuals}\NormalTok{(fitRI, }\DataTypeTok{scaled =} \OtherTok{TRUE}\NormalTok{, }\DataTypeTok{type =} \StringTok{"deviance"}\NormalTok{), }
    \DataTypeTok{Beach =} \KeywordTok{factor}\NormalTok{(RIKZ}\OperatorTok{$}\NormalTok{Beach))}
\KeywordTok{ggplot}\NormalTok{(df, }\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ NAP, }\DataTypeTok{y =}\NormalTok{ resid)) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_hline}\NormalTok{(}\DataTypeTok{yintercept =} \DecValTok{0}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{col =}\NormalTok{ Beach)) }\OperatorTok{+}\StringTok{ }
\StringTok{    }\KeywordTok{geom_smooth}\NormalTok{(}\DataTypeTok{method =} \StringTok{"loess"}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{labs}\NormalTok{(}\DataTypeTok{x =} \StringTok{"NAP"}\NormalTok{, }\DataTypeTok{y =} \StringTok{"Residuals"}\NormalTok{, }\DataTypeTok{title =} \StringTok{"Relationship between predictor and residuals"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\includegraphics{8GLMM_files/figure-latex/unnamed-chunk-11-1.pdf}

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

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{df <-}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{Beach2 =}\NormalTok{ RIKZ}\OperatorTok{$}\NormalTok{Beach, }\DataTypeTok{resid =} \KeywordTok{residuals}\NormalTok{(fitRI, }\DataTypeTok{scaled =} \OtherTok{TRUE}\NormalTok{, }
    \DataTypeTok{type =} \StringTok{"deviance"}\NormalTok{), }\DataTypeTok{Beach =} \KeywordTok{factor}\NormalTok{(RIKZ}\OperatorTok{$}\NormalTok{Beach))}
\KeywordTok{ggplot}\NormalTok{(df, }\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ Beach2, }\DataTypeTok{y =}\NormalTok{ resid)) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_hline}\NormalTok{(}\DataTypeTok{yintercept =} \DecValTok{0}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{col =}\NormalTok{ Beach)) }\OperatorTok{+}\StringTok{ }
\StringTok{    }\KeywordTok{geom_smooth}\NormalTok{(}\DataTypeTok{method =} \StringTok{"loess"}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{labs}\NormalTok{(}\DataTypeTok{x =} \StringTok{"Beach"}\NormalTok{, }\DataTypeTok{y =} \StringTok{"Residuals"}\NormalTok{, }\DataTypeTok{title =} \StringTok{"Relationship between Beach and residuals"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\includegraphics{8GLMM_files/figure-latex/unnamed-chunk-12-1.pdf}

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

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{plot_model}\NormalTok{(fitRI, }\DataTypeTok{type =} \StringTok{"diag"}\NormalTok{)}\OperatorTok{$}\NormalTok{Beach}
\end{Highlighting}
\end{Shaded}

\includegraphics{8GLMM_files/figure-latex/unnamed-chunk-13-1.pdf}

\normalsize

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

\subsection{Simulation example}\label{simulation-example}

Aim: to show the difference between random intercept and random slope
log Poisson models, and to compare the cluster curves, mean of cluster
curves and the marginal mean

This example a slightly modified version (added line and changed to
ggplot) of an example in the Lecture notes by Magne Aldrin, Norwegian
Computing Centre (which is way the notation differs slightly from what
we have used).

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

The model we use is
\[\begin{pmatrix} \gamma_{0i}\\ \gamma_{1i} \end{pmatrix} \sim N \left( \begin{pmatrix}0\\ 0 \end{pmatrix},\begin{pmatrix} 0.2^2 & 0 \\ 0 &0.07^2 \end{pmatrix} \right)\]
and then \[\eta_{ij} = \beta_0+\beta_1 x_{ij} + \gamma_{0i} \] for the
random intercept model and
\[\eta_{ij} = \beta_0+\beta_1 x_{ij} + \gamma_{0i}+\gamma_{1i}x_{ij} \]
for the random slope model.

Values chosen are \(\beta_0=-2\) and \(\beta_1=0.3\).

Log-link is used so \(\eta_{ij}=\ln(\mu_{ij})\).

Only data for the random intercept and slope are simulated, while
parameter estimates are considered known. The number of clusters is
\(m=1000\) and \(n_i = 201\) value of \(x\) from \(-10\) to \(10\) are
studied.

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

Below 30 random lines are plotted together (in different colours). The
black dashed line is \(\exp(\beta_0+\beta_i x)\) (denoted ``not mean''
in the plot) and the black solid line is the population mean (really,
the average of the \(m\) curves).

(Q: what would happen if we simulate data and estimate parameters
instead of only plotting means?)

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

\small
\includegraphics{8GLMM_files/figure-latex/unnamed-chunk-14-1.pdf}
\normalsize

\textbf{Q}: comment on what you see. Would you expect that the solid and
dashed curve were identical?

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

\subsection{Binomial random intercept
model}\label{binomial-random-intercept-model}

We will study this model in the interactive session.

\subsubsection{Simulation example}\label{simulation-example-1}

Aim: to show the difference between random intercept and random slope
logit models, and to compare the mean cluster curve with the population
curve (as for the Poisson model).

This example a slightly modified version (added line and changed to
ggplot) of an example in the Lecture notes by Magne Aldrin, Norwegian
Computing Centre (which is way the notation differs slightly from what
we have used).

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

The model we use is
\[\begin{pmatrix} \gamma_{0i}\\ \gamma_{1i} \end{pmatrix} \sim N \left( \begin{pmatrix}0\\ 0 \end{pmatrix},\begin{pmatrix} 4 & 0 \\ 0 &1 \end{pmatrix} \right)\]
and then \[\eta_{ij} = \beta_0+\beta_1 x_{ij} + \gamma_{0i} \] for the
random intercept model and
\[\eta_{ij} = \beta_0+\beta_1 x_{ij} + \gamma_{0i}+\gamma_{1i}x_{ij} \]
for the random slope model where \(\beta_0=3\) and \(\beta_1=1\).
Logit-link is used so \(\eta_{ij}=\ln(\frac{\pi_{ij}}{1-\pi_{ij}})\).

Only data for the random intercept and slope are simulated, while
parameter estimates are considered known. The number of clusters is
\(m=1000\) and \(n_i = 201\) value of \(x\) from \(-10\) to \(10\) are
studied.

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

Below 30 random lines are plotted together (in different colours). The
black dashed line is
\(\exp(\beta_0+\beta_i x)/(1+\exp(\beta_0+\beta_i x))\) (denoted ``not
mean'' in the plot) and the black solid line is the population mean
(really, the average of the \(m\) curves).

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

\small
\includegraphics{8GLMM_files/figure-latex/unnamed-chunk-16-1.pdf}

\normalsize

\textbf{Q}: comment on what you see.

\section{Conditional and marginal
model}\label{conditional-and-marginal-model}

To do parameter estimation in the LMM we started with the marginal model
for \(Y_{ij}\).

\[f(y_{ij})=\int_{\gamma_i} f(y_{ij}\mid \gamma_i) f(\gamma_i) d\gamma_i\]

where \[f(\gamma_i) \text{ is } N({\bf 0},{\bf Q})\] and
\(f(y_{ij}\mid \gamma_i)\) might be the binomial, Poisson, normal,
gamma, \ldots{} distribution - chosen as an exponential family.

It turns out that this integral can only be written out analytically in
special cases.

In the special case that \(f(y_{ij}\mid \gamma_i)\) is normal we saw in
Module LMM that the marginal distribution was also multivariate normal.

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

Another special case is the log-linear Poisson random intercept where
all \(n_i=1\). See pages 393-394 in our textbook to see the solution
(that can be used to handle overdispersion in Poisson regression).

Remark: for the random intercept Poisson GLMM with log-link the expected
marginal mean is (Agresti, 2015, page 309)
\[ \text{E}(Y_{ij})=\text{E}[\text{E}(Y_{ij}\mid \gamma_{0i})]=\exp({\bf x}_{ij}^T{\boldsymbol \beta}+\tau_0^2/2)\]
The solid line in the plot to the left above gives this curve, while the
dashed is \(\exp({\bf x}_{ij}^T{\boldsymbol \beta})\). In Agresti (2015,
page 310) the variance of the marginal distribution is given as
\[\text{Var}(Y_{ij})=\cdots = \text{E}(Y_{ij})+[\text{E}(Y_{ij})^2(\exp(\tau_0^2)-1)]\]
so we see that the variance is larger than that of a Poisson
distribution.

Remark: the conditional mean is
\(\exp({\bf x}_{ij}^T{\boldsymbol \beta}+\gamma_{0i})\), while the
marginal mean is \(\exp({\bf x}_{ij}^T{\boldsymbol \beta}+\tau_0^2/2)\).

(If the link function chosen had been the identity link then the
marginal mean would have been \({\bf x}_{ij}^T{\boldsymbol \beta}\).)

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

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

The most popular method for parameter estimation for regression is
maximum likelihood. The likelihood of the data is given by the marginal
distribution of all \(Y_{ij}\)s jointly. However, this can not be found
analytically (as we saw above), so we must resort to numerial methods to
evaluate the likelihood of a general GLMM.

The parameters we want to estimate are the fixed effects
\({\boldsymbol \beta}\)s, and the parameters in \({\bf Q}\) for the
random effects - denoted by \(\vartheta\).

The contribution from observations in the \(i\)th cluster is
\[f({\bf y}_i\mid {\boldsymbol \beta},\vartheta)=\int_{{\boldsymbol \gamma}_i} f({\bf y}_i \mid {\boldsymbol \gamma}_i,{\boldsymbol \beta}) f({\boldsymbol \gamma}_i \mid {\bf Q}) d{\boldsymbol \gamma}_i
=\int_{{\boldsymbol \gamma}_i} \prod_{j=1}^{n_i}f(y_{ij} \mid {\boldsymbol \gamma}_i,{\boldsymbol \beta}) f({\boldsymbol \gamma}_i \mid {\bf Q}) d{\boldsymbol \gamma}_i\]

And the likelihood is then (since the clusters are independent)
\[L({\boldsymbol \beta},\vartheta)=\prod_{i=1}^m f(\bf{y}_i\mid {\boldsymbol \beta},\vartheta)\]

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

To arrive at parameter estimates we need to maximize the likelihood with
respect to \({\boldsymbol \beta}\) and \(\vartheta\), which is a
complicated numerical problem.

There are several solution to this, and the \texttt{glmer} function in
\texttt{lme4} uses the Laplace approximation. Other solutions are
adaptive Gaussian quadratures, penalized quasilikelihood and various
Bayesian methods (study the posterior distribution of the parameters).
This is an active area of research.

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

\subsection{Gauss-Hermite Quadrature}\label{gauss-hermite-quadrature}

If the dimension of \({\boldsymbol \gamma}_i\) is small, for example our
random intercept model, then the Gauss-Hermite quadrature may be used to
approximate the integral above. This will be of the form

\[\int_{-\infty}^{\infty} h(\gamma)\exp(-\gamma^2)d\gamma\] which can be
written as a sum of weights and quadrature points that are roots of
Hermite polynomials. This can then be maximized using Newton-Raphson
methods.

When the dimension of \({\boldsymbol \gamma}\) gets more than two, other
methods are used.

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

\subsection{Laplace approximation (for the interested
student)}\label{laplace-approximation-for-the-interested-student}

The notation in this section is taken from Kristensen et al (2016):
``TMB: Automatic Differentiation and Laplace Approximation'',
\emph{Journal of Statistical Software}, 70, 5.
\url{doi:10.18637/jss.v070.i05}, which is the background for the
\texttt{glmmTMB} package. The notation is different from what is used in
this module, but this reflects that often different notation is used for
the same models.

Let \(f(u,\theta)\) be the \emph{joint log-likelihood of the data and
the random effects}. This will in our case be related to
\(\ln[f({\bf y}_i \mid \gamma_i,{\boldsymbol \beta}) f(\gamma_i \mid {\bf Q})]\).

Observe that \(f(u,\theta)\) depends on the unknown random effects \(u\)
and the parameters in \(\theta\) (here we have the fixed effects and the
parameters in \({\bf Q}\) for the random effects).

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

The maximum likelihood estimate for \(\theta\) maximized
\[ L(\theta)=\int_{u} \exp(f(u,\theta))du\] with respect to \(\theta\).
Observe that the random effects \(u\) are now integrated out and the
marginal likelihood \(L(\theta)\) is the likelihood of the data as a
function of the parameters only.

Let \(\hat{u}(\theta)\) be the minimizer of \(f(u,\theta)\), and
\(H(\theta)\) is the Hessian (matrix of partial second derivatives) of
\(f(u,\theta)\) evaluated at \(\hat{u}(\theta)\).

The Laplace approximation for the marginal likelihood \(L(\theta)\) is
\[ L^{*}(\theta)=(2\pi)^{n/2} \text{det}(H(\theta))^{-1/2} \exp(f(\hat{u},\theta))\]

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

Some more details: Let \(g(u)=\ln(f(y\mid u)f(u))\), and look at the
second order Taylor expansion around \(\hat{u}(\theta)\):
\[g(u) \approx \tilde{g}(u)=g(\hat{u})+(u-\hat{u})g`(\hat{u})+\frac{1}{2}(u-\hat{u})^2 g''(\hat{u})=g(\hat{u})-\frac{1}{2}(u-\hat{u})^2 (-g''(\hat{u}))\]
This means that \(\exp(\tilde{g}(u))\) is proportional to the normal
density with mean \(\mu=\hat{u}\) and variance
\(\sigma^2=-1/g''(\hat{u})\). Putting this back into our integral
\[L(\theta)=\int_{u} \exp(g(u))du \approx\int_{u} \exp(\tilde{g}(u))du=\exp(g(\hat{u}))\int_u \exp(-\frac{1}{2 \sigma^2}(u-\mu)^2)du=\exp(g(\hat{u})\sqrt{2\pi \sigma^2}\]
This last part was taken from
\url{http://people.math.aau.dk/~rw/Undervisning/Topics/Handouts/6.hand.pdf}.

\section{What have we not covered?}\label{what-have-we-not-covered}

\begin{itemize}
\tightlist
\item
  We have skipped a lot of technicalities about the parameter
  estimation.
\item
  How to predict the random effects \(\gamma_i\): ``the conditional
  modes'' of the random effects.
\item
  What is fitted values, and how are residuals calculated.
\item
  How to test hypotheses?
\item
  AIC to compare models (very similar to LMM).
\item
  And, surely, much more.
\end{itemize}

\section{Summing up - what have we learned about the
GLMM?}\label{summing-up---what-have-we-learned-about-the-glmm}

\begin{itemize}
\tightlist
\item
  The GLMM can be formulated using three ingredients:

  \begin{itemize}
  \tightlist
  \item
    distributional assumption: \(f(y_{ij}\mid \gamma_i)\) from
    exponential family
  \item
    structural assumptions: linear predictor
    \(\eta_{ij}={\bf x}_{ij}^T {\boldsymbol \beta}+{\bf u}_{ij}\gamma_i\)
    and link function \(\eta_{ij}=g(\mu_{ij})\) where
    \(\mu_{ij}=\text{E}(Y_{ij}\mid \gamma_i)\).
  \item
    distributional assumptions for the random effects:
    \(\gamma_i\sim N({\bf 0},{\bf Q})\).
  \end{itemize}
\item
  The GLMM likelihood function is expressed as an integral with respect
  to the random effects and does (in general) not have a closed form
  solution.
\item
  Numerical approximation methods need to be used to find parameter
  estimates, and one possibility is to use the Laplace approximation,
  but many competing method exists.
\end{itemize}

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

\begin{itemize}
\tightlist
\item
  Three \texttt{R} packages that can be investigated is
  \texttt{lme4}(function \texttt{glmer}) and \texttt{glmmTMB} (template
  model builder), and the NTNU-flagship \texttt{inla}. How to use these
  three packages on a simulated data set (binary data, logit link,
  random intercept and slope) is shown in the end of the module page
  (NOT on our reading list but for the interested student).
\end{itemize}

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

This is the last interactive session of the course.

\begin{itemize}
\item
  First hour: Problem 1
\item
  15 min before end first hour (or if you finish earlier): answer the
  following questions in your group (togehter) and fill in
  \href{https://goo.gl/forms/lcW6AXDgcXpIMHsG3}{this Google form}
\item
  What was the top 5 (or so) concepts in this course that you feel that
  you need to work more to understand? (input to summing-up last
  lecture)
\item
  You are here today - you are the students that come to the IL, even
  the last one. Why do you attend the IL? What is good about the IL?
  What can be challenging?
\item
  If you want:
  \href{https://www.math.ntnu.no/emner/TMA4315/2018h/Summingup2018draft.pdf}{here
  is the draft for the summing-up} (module 6 on categorical missing and
  some QC needed) - any comments?
\item
  Second hour: continue the work above, or work on (and get help) with
  Compulsory 3!
\end{itemize}

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

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

\textbf{a)} Show that the binomial distribution is a member of the
exponential family of distributions. What do we mean with canonical
link, and which advantages do we get when using canonical link? What is
the canonical link for the binomial distibution?

Remark: this has been done several times before. Instead of doing the
maths on the board, focus on the concepts. Also discuss this
relationship:
\(\theta_i\leftrightarrow \mu_i\leftrightarrow \eta_i \leftrightarrow \mathbf{\beta}\)
and fill in what is over the arrows, and the missing arrow for the
canonical link.

\textbf{b)} \emph{This excercise is modified!}

Assume the following model:

\[y_{ij} \sim \text{Bin}(1, \pi_{ij})\]

\[ g(\pi_{ij}) = \beta_0 + \beta_1x_{1,ij} + \beta_2x_{2,i} \]

where \(g(\cdot)\) is a suitable link function and all \(y_{ij}\) are
independent. What kind of model is this?\\
The table below shows AIC-values for three different link functions
often used with binary response data:

\begin{longtable}[]{@{}lr@{}}
\toprule
Link function & AIC\tabularnewline
\midrule
\endhead
probit & 57.51571\tabularnewline
cloglog & 58.30338\tabularnewline
logit & 57.69237\tabularnewline
\bottomrule
\end{longtable}

Explain what AIC is, and argue why it is reasonable to use such criteria
in this situation (instead of other tests from this course). Based on
these three values, which link-function do you prefer? Why?

Remark: what is the take home message here? Hint: nested models?

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

We will, as we did a few weeks back, look at a fish dataset consisting
of observations of cod from the year 2000, but we are now interested in
the age of each fish, rather than the weight. The dataset is from
Havforskningsinstituttet in Bergen, and we use only a subset of a large
dataset with information about fish in Barentshavet. The following
variables are available on a random sample within each catch (a draft of
trawl (trål)):

\begin{itemize}
\tightlist
\item
  \texttt{length}: the length of fish (in cm)
\item
  \texttt{weight}: the weight of fish (in grams)
\item
  \texttt{age}: the age of fish (in years)
\item
  \texttt{haulsize}: the size of the total catch (in ton (1000 kg))
\end{itemize}

Let \(i\) be the index for catch, and \(j\) be the index for an
individual fish in a given catch.

Age is a categorical variable between 2 and 12, but we create a new
variable \(A_{ij}\):

\[ A_{ij} = \begin{cases} 1 \text{ if age}_{ij} > 6 \\ 0 \text{ else} \end{cases} \]

which is a binary variable, and we use this for the response.

\textbf{Remark: before we used weight as response (normal), but now we
use dichotomized age as a binary response}.

We look at the following model:

\[ A_{ij}|\gamma_{0i} \sim \text{Bin}(1, \pi_{ij}) \]
\[ g(\pi_{ij}) = \beta_0 + \beta_1\log(\texttt{length}_{ij}) + \beta_2\log(\texttt{haulsize}_i) + \gamma_{0i} \]

\[ \gamma_{0i} \sim N(0, \tau_0^2) \]

where \(g(\cdot)\) is the logit link function, and all \(A_{ij}\) are
conditional independent (given \(\gamma_i\)).

\textbf{c1)} \emph{New:} What kind of model is this? What is
\(\gamma_{0i}\)? Compare to the linear mixed effects model- similarities
and differences.

Below you find an excerpt of this model:

\begin{verbatim}
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: A ~ log(length) + log(haulsize) + (1 | haul)
##    Data: fish
## 
##      AIC      BIC   logLik deviance df.resid 
##   2524.7   2551.2  -1258.4   2516.7     5483 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -20.037  -0.176  -0.040   0.138  40.513 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  haul   (Intercept) 1.921    1.386   
## Number of obs: 5487, groups:  haul, 63
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -89.3202     3.0256 -29.522   <2e-16 ***
## log(length)    20.7783     0.7080  29.348   <2e-16 ***
## log(haulsize)  -0.2396     0.1400  -1.711    0.087 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg(ln)
## log(length) -0.997       
## log(haulsz)  0.021  0.017
\end{verbatim}

\textbf{c2)} \emph{New:} Explain what it means that the model is fitted
by maximum likelihood. What does Laplace approximation has to do with
the model fit?

\textbf{Skip d)- not in our reading list to test random effects - so too
advanced. } Instead discuss why this is difficult.

\textbf{d)} The log-likelihood value for the corresponding model without
any random effects (which corresponds to a model with \(\tau_0^2 = 0\))
is -1498.04. Use this to perform a likelihood ratio test with
\(H_0: \tau_0^2 = 0\). Calculate the corresponding p-value and draw a
conclusion.

\textbf{e)} The Havforskningsinstituttet believes that the haulsize is
important in the modelling of age and length of the fish (read: as a
fixed effect). Based on the excerpts above (on dichotomized age) and
below (on weight), what do you think about this?

\begin{verbatim}
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(weight) ~ log(length) + log(haulsize) + (1 | haul)
##    Data: fish
## 
## REML criterion at convergence: -10188.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1204 -0.6107 -0.0314  0.5813  5.1458 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  haul     (Intercept) 0.003030 0.05505 
##  Residual             0.008594 0.09270 
## Number of obs: 5433, groups:  haul, 63
## 
## Fixed effects:
##                Estimate Std. Error  t value
## (Intercept)   -4.286029   0.041554 -103.143
## log(length)    2.901037   0.009703  298.984
## log(haulsize)  0.006169   0.005328    1.158
## 
## Correlation of Fixed Effects:
##             (Intr) lg(ln)
## log(length) -0.980       
## log(haulsz)  0.047  0.056
\end{verbatim}

\textbf{Remark}: It is in general a very bad idea to do a
dichotomization of a continuous variable - like what is done here for
transforming the age into \(A_{ij}\), because then information is lost.
In a hypothesis testing set-up this will in general give substantial
loss of power to detect the effect of a covariate.

\section{\texorpdfstring{Software for GLMM: demonstration of
\texttt{glmer}, \texttt{glmmTMB} and \texttt{inla} for analysing
GLMMs}{Software for GLMM: demonstration of glmer, glmmTMB and inla for analysing GLMMs}}\label{software-for-glmm-demonstration-of-glmer-glmmtmb-and-inla-for-analysing-glmms}

This is \textbf{NOT} on the reading list, but an extra service to the
interested students.

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"arm"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"reshape2"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"sp"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"glmmTBM"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"INLA"}\NormalTok{, }\DataTypeTok{repos=}\StringTok{"https://inla.r-inla-download.org/R/stable"}\NormalTok{, }\DataTypeTok{dep=}\OtherTok{TRUE}\NormalTok{)}
\ErrorTok{\}}
\end{Highlighting}
\end{Shaded}

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(ggplot2)}
\KeywordTok{library}\NormalTok{(ggpubr)}
\KeywordTok{library}\NormalTok{(lme4)}
\KeywordTok{library}\NormalTok{(glmmTMB)}
\KeywordTok{library}\NormalTok{(INLA)}
\KeywordTok{library}\NormalTok{(reshape2)}
\KeywordTok{library}\NormalTok{(arm)}
\end{Highlighting}
\end{Shaded}

For this demonstration, we will simulate a dataset with one fixed effect
(x) and random intercept and random slope (always nice to know the true
values for parameters to compare with estimates).

\textbf{Q}: What is the model that we simulate from?

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{set.seed}\NormalTok{(}\DecValTok{90}\NormalTok{) }\CommentTok{# to get reproduciability}

\NormalTok{no.gr <-}\StringTok{ }\DecValTok{15}
\NormalTok{x1 <-}\StringTok{ }\KeywordTok{runif}\NormalTok{(}\DecValTok{200}\NormalTok{, }\DecValTok{0}\NormalTok{, }\DecValTok{1}\NormalTok{)}

\NormalTok{y <-}\StringTok{ }\KeywordTok{matrix}\NormalTok{(}\OtherTok{NA}\NormalTok{, }\DataTypeTok{nrow =} \DecValTok{200}\NormalTok{, }\DataTypeTok{ncol =} \DecValTok{15}\NormalTok{)}

\NormalTok{Q <-}\StringTok{ }\KeywordTok{matrix}\NormalTok{(}\KeywordTok{c}\NormalTok{(}\FloatTok{0.5}\NormalTok{, }\FloatTok{0.3}\NormalTok{, }\FloatTok{0.3}\NormalTok{, }\FloatTok{0.5}\NormalTok{), }\DataTypeTok{nrow =} \DecValTok{2}\NormalTok{)}
\NormalTok{random_both <-}\StringTok{ }\KeywordTok{mvrnorm}\NormalTok{(no.gr, }\KeywordTok{c}\NormalTok{(}\DecValTok{0}\NormalTok{,}\DecValTok{0}\NormalTok{), Q)}
\NormalTok{random_int <-}\StringTok{ }\NormalTok{random_both[,}\DecValTok{1}\NormalTok{]; random_slope <-}\StringTok{ }\NormalTok{random_both[,}\DecValTok{2}\NormalTok{]}

\CommentTok{# to test out}
\CommentTok{#random_int <- runif(15, -2, 2) # when we want to simulate from a model that does not fit exactly}
\CommentTok{#random_slope <- runif(15, -2, 2)}

\NormalTok{beta_}\DecValTok{0}\NormalTok{ <-}\StringTok{ }\OperatorTok{-}\FloatTok{0.45}
\NormalTok{beta_}\DecValTok{1}\NormalTok{ <-}\StringTok{ }\FloatTok{0.76}

\ControlFlowTok{for}\NormalTok{ (i }\ControlFlowTok{in} \DecValTok{1}\OperatorTok{:}\DecValTok{15}\NormalTok{)\{}
\NormalTok{  y[,i] <-}\StringTok{ }\NormalTok{beta_}\DecValTok{0} \OperatorTok{+}\StringTok{ }\NormalTok{beta_}\DecValTok{1}\OperatorTok{*}\NormalTok{x1 }\OperatorTok{+}\StringTok{ }\NormalTok{random_int[i] }\OperatorTok{+}\StringTok{ }\NormalTok{random_slope[i]}\OperatorTok{*}\NormalTok{x1 }\OperatorTok{+}\StringTok{ }\KeywordTok{rnorm}\NormalTok{(}\DecValTok{200}\NormalTok{, }\DecValTok{0}\NormalTok{, }\DecValTok{1}\NormalTok{)}
\NormalTok{\}}

\NormalTok{y_binom <-}\StringTok{ }\KeywordTok{as.numeric}\NormalTok{(}\KeywordTok{c}\NormalTok{(y) }\OperatorTok{>}\StringTok{ }\DecValTok{0}\NormalTok{)}

\KeywordTok{ggplot}\NormalTok{(}\KeywordTok{cbind}\NormalTok{(}\KeywordTok{melt}\NormalTok{(}\KeywordTok{data.frame}\NormalTok{(y)), }\DataTypeTok{x =} \KeywordTok{rep}\NormalTok{(x1, }\DecValTok{15}\NormalTok{), }\DataTypeTok{y =} \KeywordTok{as.factor}\NormalTok{(y_binom)), }
       \DataTypeTok{mapping =} \KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ x, }\DataTypeTok{y =}\NormalTok{ value)) }\OperatorTok{+}
\StringTok{  }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{col =}\NormalTok{ y)) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_smooth}\NormalTok{(}\DataTypeTok{method =} \StringTok{"lm"}\NormalTok{, }\DataTypeTok{col =} \StringTok{"black"}\NormalTok{, }\DataTypeTok{se =} \OtherTok{FALSE}\NormalTok{) }\OperatorTok{+}
\StringTok{  }\KeywordTok{facet_wrap}\NormalTok{(}\OperatorTok{~}\StringTok{ }\NormalTok{variable) }\OperatorTok{+}\StringTok{ }\KeywordTok{labs}\NormalTok{(}\DataTypeTok{x =} \StringTok{"x"}\NormalTok{, }\DataTypeTok{y =} \StringTok{"y"}\NormalTok{) }
\end{Highlighting}
\end{Shaded}

\includegraphics{8GLMM_files/figure-latex/unnamed-chunk-22-1.pdf}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{mydata <-}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{y =}\NormalTok{ y_binom, }\DataTypeTok{x =}\NormalTok{ x1, }\DataTypeTok{group =} \KeywordTok{rep}\NormalTok{(}\DecValTok{1}\OperatorTok{:}\DecValTok{15}\NormalTok{, }\DataTypeTok{each =} \DecValTok{200}\NormalTok{))}

\KeywordTok{head}\NormalTok{(mydata)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##   y         x group
## 1 0 0.5307603     1
## 2 0 0.8850698     1
## 3 0 0.4399639     1
## 4 1 0.9422908     1
## 5 1 0.1883275     1
## 6 0 0.6555310     1
\end{verbatim}

Use a binary regression with logit link to model the probability. That
is, we have one covariate and make the following assumptions.

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  Distributional assumptions:
  \(Y_{ij}\mid \gamma_i \sim \text{Bin}(n_{i}, \pi_{ij}) \text{ for } i = 1, \dots, 15, \ j = 1, \dots, n_i = 200, \ N = \sum_i n_i = 3000\).
\item
  Structural assumptions: The linear predictor (with random effects)
  with one fixed covariate:
  \[\eta_{ij} = \beta_0+\beta_1 x_{ij} + \gamma_{0i}+\gamma_{1i}x_{ij} \]
  where we will choose \(x_{ij}=u_{ij}\). The link function is:
  \[\eta_{ij} = \ln \left(\frac{\pi_i}{1-\pi_i}\right)\]
\item
  Distributional assumptions for random effects
  \[\begin{pmatrix} \gamma_{0i}\\ \gamma_{1i} \end{pmatrix} \sim N \left( \begin{pmatrix}0\\ 0 \end{pmatrix},\begin{pmatrix} \tau_0^2 & \tau_{01} \\ \tau_{01} &\tau_1^2 \end{pmatrix} \right)\]
\end{enumerate}

Now we will look at three common ways of analysing such models. For each
type we will fit the model and look at results.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{fit_glmer <-}\StringTok{ }\KeywordTok{glmer}\NormalTok{(y }\OperatorTok{~}\StringTok{ }\NormalTok{x }\OperatorTok{+}\StringTok{ }\NormalTok{(}\DecValTok{1} \OperatorTok{+}\StringTok{ }\NormalTok{x }\OperatorTok{|}\StringTok{ }\NormalTok{group), }\DataTypeTok{data =}\NormalTok{ mydata, }\DataTypeTok{family =} \KeywordTok{binomial}\NormalTok{(}\DataTypeTok{link =} \StringTok{"logit"}\NormalTok{))}

\NormalTok{fit_glmmTMB <-}\StringTok{ }\KeywordTok{glmmTMB}\NormalTok{(y }\OperatorTok{~}\StringTok{ }\NormalTok{x }\OperatorTok{+}\StringTok{ }\NormalTok{(}\DecValTok{1} \OperatorTok{+}\StringTok{ }\NormalTok{x }\OperatorTok{|}\StringTok{ }\NormalTok{group), }\DataTypeTok{data =}\NormalTok{ mydata, }\DataTypeTok{family =} \KeywordTok{binomial}\NormalTok{(}\DataTypeTok{link =} \StringTok{"logit"}\NormalTok{))}

\NormalTok{n.gr <-}\StringTok{ }\KeywordTok{max}\NormalTok{(mydata}\OperatorTok{$}\NormalTok{group)}
\NormalTok{mydata}\OperatorTok{$}\NormalTok{group2 <-}\StringTok{ }\NormalTok{mydata}\OperatorTok{$}\NormalTok{group }\OperatorTok{+}\StringTok{ }\NormalTok{n.gr}
\NormalTok{fit_inla <-}\StringTok{ }\KeywordTok{inla}\NormalTok{(y }\OperatorTok{~}\StringTok{ }\NormalTok{x }\OperatorTok{+}\StringTok{ }\KeywordTok{f}\NormalTok{(group, }\DataTypeTok{model =} \StringTok{"iid2d"}\NormalTok{, }\DataTypeTok{n =} \DecValTok{2}\OperatorTok{*}\NormalTok{n.gr) }\OperatorTok{+}\StringTok{ }\KeywordTok{f}\NormalTok{(group2, x, }\DataTypeTok{copy =} \StringTok{"group"}\NormalTok{),}
                 \DataTypeTok{data =}\NormalTok{ mydata, }\DataTypeTok{family =} \StringTok{"binomial"}\NormalTok{,}
                 \DataTypeTok{control.compute =} \KeywordTok{list}\NormalTok{(}\DataTypeTok{dic =} \OtherTok{TRUE}\NormalTok{),}
                 \DataTypeTok{control.family =} \KeywordTok{list}\NormalTok{(}\DataTypeTok{link =} \StringTok{"logit"}\NormalTok{))}
\end{Highlighting}
\end{Shaded}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{coefdf <-}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{mean =} 
                       \KeywordTok{c}\NormalTok{(}\KeywordTok{summary}\NormalTok{(fit_glmer)}\OperatorTok{$}\NormalTok{coefficients[,}\DecValTok{1}\NormalTok{], }
                         \KeywordTok{summary}\NormalTok{(fit_glmmTMB)}\OperatorTok{$}\NormalTok{coefficients}\OperatorTok{$}\NormalTok{cond[,}\DecValTok{1}\NormalTok{], }
\NormalTok{                         fit_inla}\OperatorTok{$}\NormalTok{summary.fixed[,}\DecValTok{1}\NormalTok{]),}
                     \DataTypeTok{sd =} 
                       \KeywordTok{c}\NormalTok{(}\KeywordTok{summary}\NormalTok{(fit_glmer)}\OperatorTok{$}\NormalTok{coefficients[,}\DecValTok{2}\NormalTok{], }
                         \KeywordTok{summary}\NormalTok{(fit_glmmTMB)}\OperatorTok{$}\NormalTok{coefficients}\OperatorTok{$}\NormalTok{cond[,}\DecValTok{2}\NormalTok{], }
\NormalTok{                         fit_inla}\OperatorTok{$}\NormalTok{summary.fixed[,}\DecValTok{2}\NormalTok{]),}
                     \DataTypeTok{mod =} \KeywordTok{rep}\NormalTok{(}\KeywordTok{c}\NormalTok{(}\StringTok{"glmer"}\NormalTok{, }\StringTok{"glmmTMB"}\NormalTok{, }\StringTok{"INLA"}\NormalTok{), }\DataTypeTok{each =} \KeywordTok{length}\NormalTok{(fit_inla}\OperatorTok{$}\NormalTok{names.fixed)),}
                     \DataTypeTok{par =} \KeywordTok{rep}\NormalTok{(fit_inla}\OperatorTok{$}\NormalTok{names.fixed, }\DecValTok{3}\NormalTok{))}

\NormalTok{true_frame <-}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{par =}\NormalTok{ fit_inla}\OperatorTok{$}\NormalTok{names.fixed, }\DataTypeTok{beta =} \KeywordTok{c}\NormalTok{(beta_}\DecValTok{0}\NormalTok{, beta_}\DecValTok{1}\NormalTok{))}

\NormalTok{critval <-}\StringTok{ }\KeywordTok{qnorm}\NormalTok{(}\FloatTok{0.025}\NormalTok{, }\DataTypeTok{lower.tail =} \OtherTok{FALSE}\NormalTok{)}
\KeywordTok{ggplot}\NormalTok{(coefdf) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ mod, }\DataTypeTok{y =}\NormalTok{ mean)) }\OperatorTok{+}\StringTok{ }
\StringTok{  }\KeywordTok{geom_hline}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{yintercept =} \DecValTok{0}\NormalTok{), }\DataTypeTok{col =} \StringTok{"#D55E00"}\NormalTok{) }\OperatorTok{+}
\StringTok{  }\KeywordTok{geom_hline}\NormalTok{(}\DataTypeTok{data =}\NormalTok{ true_frame, }\KeywordTok{aes}\NormalTok{(}\DataTypeTok{yintercept =}\NormalTok{ beta), }\DataTypeTok{col =} \StringTok{"forestgreen"}\NormalTok{, }\DataTypeTok{size =} \DecValTok{1}\NormalTok{) }\OperatorTok{+}
\StringTok{  }\KeywordTok{geom_errorbar}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ mod, }\DataTypeTok{ymin =}\NormalTok{ mean}\OperatorTok{-}\NormalTok{critval}\OperatorTok{*}\NormalTok{sd, }\DataTypeTok{ymax =}\NormalTok{ mean}\OperatorTok{+}\NormalTok{critval}\OperatorTok{*}\NormalTok{sd)) }\OperatorTok{+}
\StringTok{  }\KeywordTok{facet_wrap}\NormalTok{(}\OperatorTok{~}\StringTok{ }\NormalTok{par, }\DataTypeTok{scales =} \StringTok{"free_y"}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{labs}\NormalTok{ (}\DataTypeTok{x =} \StringTok{""}\NormalTok{, }\DataTypeTok{y =} \StringTok{""}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\includegraphics{8GLMM_files/figure-latex/unnamed-chunk-24-1.pdf}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{randdf_intercept <-}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{mean =} 
                      \KeywordTok{c}\NormalTok{(}\KeywordTok{ranef}\NormalTok{(fit_glmer)}\OperatorTok{$}\NormalTok{group[,}\DecValTok{1}\NormalTok{],}
                        \KeywordTok{ranef}\NormalTok{(fit_glmmTMB)}\OperatorTok{$}\NormalTok{cond}\OperatorTok{$}\NormalTok{group[,}\DecValTok{1}\NormalTok{],}
\NormalTok{                        fit_inla}\OperatorTok{$}\NormalTok{summary.random}\OperatorTok{$}\NormalTok{group}\OperatorTok{$}\NormalTok{mean[}\KeywordTok{unique}\NormalTok{(mydata}\OperatorTok{$}\NormalTok{group)],}
\NormalTok{                        random_int),}
                    \DataTypeTok{mod =} \KeywordTok{rep}\NormalTok{(}\KeywordTok{c}\NormalTok{(}\StringTok{"glmer"}\NormalTok{, }\StringTok{"glmmTMB"}\NormalTok{, }\StringTok{"INLA"}\NormalTok{, }\StringTok{"true"}\NormalTok{), }\DataTypeTok{each =} \KeywordTok{length}\NormalTok{(}\KeywordTok{unique}\NormalTok{(mydata}\OperatorTok{$}\NormalTok{group))),}
                    \DataTypeTok{x =} \KeywordTok{paste0}\NormalTok{(}\StringTok{"y"}\NormalTok{, }\KeywordTok{sprintf}\NormalTok{(}\StringTok{"%02d"}\NormalTok{, }\KeywordTok{rep}\NormalTok{(}\DecValTok{1}\OperatorTok{:}\KeywordTok{length}\NormalTok{(}\KeywordTok{unique}\NormalTok{(mydata}\OperatorTok{$}\NormalTok{group)), }\DecValTok{4}\NormalTok{))))}

\NormalTok{critval <-}\StringTok{ }\KeywordTok{qnorm}\NormalTok{(}\FloatTok{0.025}\NormalTok{, }\DataTypeTok{lower.tail =} \OtherTok{FALSE}\NormalTok{)}
\NormalTok{randdf_intercept}\OperatorTok{$}\NormalTok{low <-}\StringTok{ }\KeywordTok{c}\NormalTok{(}\KeywordTok{se.ranef}\NormalTok{(fit_glmer)}\OperatorTok{$}\NormalTok{group[,}\DecValTok{1}\NormalTok{]}\OperatorTok{*}\NormalTok{(}\OperatorTok{-}\NormalTok{critval) }\OperatorTok{+}\StringTok{ }\KeywordTok{ranef}\NormalTok{(fit_glmer)}\OperatorTok{$}\NormalTok{group[,}\DecValTok{1}\NormalTok{],}
                          \KeywordTok{ranef}\NormalTok{(fit_glmmTMB)}\OperatorTok{$}\NormalTok{cond}\OperatorTok{$}\NormalTok{group[,}\DecValTok{1}\NormalTok{], }\CommentTok{# do not know how to get these values for glmmTMB}
\NormalTok{                          fit_inla}\OperatorTok{$}\NormalTok{summary.random}\OperatorTok{$}\NormalTok{group}\OperatorTok{$}\StringTok{`}\DataTypeTok{0.025quant}\StringTok{`}\NormalTok{[}\KeywordTok{unique}\NormalTok{(mydata}\OperatorTok{$}\NormalTok{group)],}
                          \KeywordTok{rep}\NormalTok{(}\DecValTok{0}\NormalTok{, }\DecValTok{15}\NormalTok{))}
\NormalTok{randdf_intercept}\OperatorTok{$}\NormalTok{high <-}\StringTok{ }\KeywordTok{c}\NormalTok{(}\KeywordTok{se.ranef}\NormalTok{(fit_glmer)}\OperatorTok{$}\NormalTok{group[,}\DecValTok{1}\NormalTok{]}\OperatorTok{*}\NormalTok{(critval) }\OperatorTok{+}\StringTok{ }\KeywordTok{ranef}\NormalTok{(fit_glmer)}\OperatorTok{$}\NormalTok{group[,}\DecValTok{1}\NormalTok{],}
                          \KeywordTok{ranef}\NormalTok{(fit_glmmTMB)}\OperatorTok{$}\NormalTok{cond}\OperatorTok{$}\NormalTok{group[,}\DecValTok{1}\NormalTok{], }\CommentTok{# do not know how to get these values for glmmTMB}
\NormalTok{                          fit_inla}\OperatorTok{$}\NormalTok{summary.random}\OperatorTok{$}\NormalTok{group}\OperatorTok{$}\StringTok{`}\DataTypeTok{0.975quant}\StringTok{`}\NormalTok{[}\KeywordTok{unique}\NormalTok{(mydata}\OperatorTok{$}\NormalTok{group)],}
                          \KeywordTok{rep}\NormalTok{(}\DecValTok{0}\NormalTok{, }\DecValTok{15}\NormalTok{))}
\NormalTok{randdf_intercept}\OperatorTok{$}\NormalTok{x2 <-}\StringTok{ }\KeywordTok{rep}\NormalTok{(}\DecValTok{1}\OperatorTok{:}\DecValTok{15}\NormalTok{, }\DecValTok{4}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{rep}\NormalTok{(}\KeywordTok{c}\NormalTok{(}\OperatorTok{-}\FloatTok{0.15}\NormalTok{,}\OperatorTok{-}\FloatTok{0.05}\NormalTok{,}\FloatTok{0.05}\NormalTok{,}\FloatTok{0.15}\NormalTok{), }\DataTypeTok{each =} \DecValTok{15}\NormalTok{)}

\KeywordTok{ggplot}\NormalTok{(randdf_intercept) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ mean, }\DataTypeTok{y =}\NormalTok{ x2, }\DataTypeTok{col =}\NormalTok{ mod)) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_vline}\NormalTok{(}\DataTypeTok{xintercept =} \DecValTok{0}\NormalTok{, }\DataTypeTok{lty =} \DecValTok{2}\NormalTok{) }\OperatorTok{+}
\StringTok{  }\KeywordTok{geom_segment}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ low, }\DataTypeTok{xend =}\NormalTok{ high, }\DataTypeTok{y =}\NormalTok{ x2, }\DataTypeTok{yend =}\NormalTok{ x2, }\DataTypeTok{col =}\NormalTok{ mod)) }\OperatorTok{+}
\StringTok{  }\KeywordTok{labs}\NormalTok{(}\DataTypeTok{x =} \StringTok{"random intercept"}\NormalTok{, }\DataTypeTok{y =} \StringTok{""}\NormalTok{) }\OperatorTok{+}\StringTok{ }
\StringTok{  }\KeywordTok{scale_y_continuous}\NormalTok{(}\DataTypeTok{breaks =} \DecValTok{1}\OperatorTok{:}\DecValTok{15}\NormalTok{, }\DataTypeTok{labels =} \KeywordTok{paste0}\NormalTok{(}\StringTok{"y"}\NormalTok{, }\KeywordTok{sprintf}\NormalTok{(}\StringTok{"%02d"}\NormalTok{, }\KeywordTok{rep}\NormalTok{(}\DecValTok{1}\OperatorTok{:}\KeywordTok{length}\NormalTok{(}\KeywordTok{unique}\NormalTok{(mydata}\OperatorTok{$}\NormalTok{group))))))}
\end{Highlighting}
\end{Shaded}

\includegraphics{8GLMM_files/figure-latex/unnamed-chunk-24-2.pdf}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{randdf_slope <-}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{mean =}
                      \KeywordTok{c}\NormalTok{(}\KeywordTok{ranef}\NormalTok{(fit_glmer)}\OperatorTok{$}\NormalTok{group[,}\DecValTok{2}\NormalTok{],}
                        \KeywordTok{ranef}\NormalTok{(fit_glmmTMB)}\OperatorTok{$}\NormalTok{cond}\OperatorTok{$}\NormalTok{group[,}\DecValTok{2}\NormalTok{],}
\NormalTok{                        fit_inla}\OperatorTok{$}\NormalTok{summary.random}\OperatorTok{$}\NormalTok{group}\OperatorTok{$}\NormalTok{mean[}\KeywordTok{unique}\NormalTok{(mydata}\OperatorTok{$}\NormalTok{group2)],}
\NormalTok{                        random_slope),}
                    \DataTypeTok{mod =} \KeywordTok{rep}\NormalTok{(}\KeywordTok{c}\NormalTok{(}\StringTok{"glmer"}\NormalTok{, }\StringTok{"glmmTMB"}\NormalTok{, }\StringTok{"INLA"}\NormalTok{, }\StringTok{"true"}\NormalTok{), }\DataTypeTok{each =} \KeywordTok{length}\NormalTok{(}\KeywordTok{unique}\NormalTok{(mydata}\OperatorTok{$}\NormalTok{group))),}
                    \DataTypeTok{x =} \KeywordTok{paste0}\NormalTok{(}\StringTok{"y"}\NormalTok{, }\KeywordTok{sprintf}\NormalTok{(}\StringTok{"%02d"}\NormalTok{, }\KeywordTok{rep}\NormalTok{(}\DecValTok{1}\OperatorTok{:}\KeywordTok{length}\NormalTok{(}\KeywordTok{unique}\NormalTok{(mydata}\OperatorTok{$}\NormalTok{group)), }\DecValTok{4}\NormalTok{))))}


\NormalTok{critval <-}\StringTok{ }\KeywordTok{qnorm}\NormalTok{(}\FloatTok{0.025}\NormalTok{, }\DataTypeTok{lower.tail =} \OtherTok{FALSE}\NormalTok{)}
\NormalTok{randdf_slope}\OperatorTok{$}\NormalTok{low <-}\StringTok{ }\KeywordTok{c}\NormalTok{(}\KeywordTok{se.ranef}\NormalTok{(fit_glmer)}\OperatorTok{$}\NormalTok{group[,}\DecValTok{2}\NormalTok{]}\OperatorTok{*}\NormalTok{(}\OperatorTok{-}\NormalTok{critval) }\OperatorTok{+}\StringTok{ }\KeywordTok{ranef}\NormalTok{(fit_glmer)}\OperatorTok{$}\NormalTok{group[,}\DecValTok{2}\NormalTok{],}
                          \KeywordTok{ranef}\NormalTok{(fit_glmmTMB)}\OperatorTok{$}\NormalTok{cond}\OperatorTok{$}\NormalTok{group[,}\DecValTok{2}\NormalTok{], }\CommentTok{# do not know how to get these values for glmmTMB}
\NormalTok{                          fit_inla}\OperatorTok{$}\NormalTok{summary.random}\OperatorTok{$}\NormalTok{group}\OperatorTok{$}\StringTok{`}\DataTypeTok{0.025quant}\StringTok{`}\NormalTok{[}\KeywordTok{unique}\NormalTok{(mydata}\OperatorTok{$}\NormalTok{group2)],}
                          \KeywordTok{rep}\NormalTok{(}\DecValTok{0}\NormalTok{, }\DecValTok{15}\NormalTok{))}
\NormalTok{randdf_slope}\OperatorTok{$}\NormalTok{high <-}\StringTok{ }\KeywordTok{c}\NormalTok{(}\KeywordTok{se.ranef}\NormalTok{(fit_glmer)}\OperatorTok{$}\NormalTok{group[,}\DecValTok{2}\NormalTok{]}\OperatorTok{*}\NormalTok{(critval) }\OperatorTok{+}\StringTok{ }\KeywordTok{ranef}\NormalTok{(fit_glmer)}\OperatorTok{$}\NormalTok{group[,}\DecValTok{2}\NormalTok{],}
                          \KeywordTok{ranef}\NormalTok{(fit_glmmTMB)}\OperatorTok{$}\NormalTok{cond}\OperatorTok{$}\NormalTok{group[,}\DecValTok{2}\NormalTok{], }\CommentTok{# do not know how to get these values for glmmTMB}
\NormalTok{                          fit_inla}\OperatorTok{$}\NormalTok{summary.random}\OperatorTok{$}\NormalTok{group}\OperatorTok{$}\StringTok{`}\DataTypeTok{0.975quant}\StringTok{`}\NormalTok{[}\KeywordTok{unique}\NormalTok{(mydata}\OperatorTok{$}\NormalTok{group2)],}
                          \KeywordTok{rep}\NormalTok{(}\DecValTok{0}\NormalTok{, }\DecValTok{15}\NormalTok{))}
\NormalTok{randdf_slope}\OperatorTok{$}\NormalTok{x2 <-}\StringTok{ }\KeywordTok{rep}\NormalTok{(}\DecValTok{1}\OperatorTok{:}\DecValTok{15}\NormalTok{, }\DecValTok{4}\NormalTok{) }\OperatorTok{+}\StringTok{ }\KeywordTok{rep}\NormalTok{(}\KeywordTok{c}\NormalTok{(}\OperatorTok{-}\FloatTok{0.15}\NormalTok{,}\OperatorTok{-}\FloatTok{0.05}\NormalTok{,}\FloatTok{0.05}\NormalTok{,}\FloatTok{0.15}\NormalTok{), }\DataTypeTok{each =} \DecValTok{15}\NormalTok{)}

\KeywordTok{ggplot}\NormalTok{(randdf_slope) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ mean, }\DataTypeTok{y =}\NormalTok{ x2, }\DataTypeTok{col =}\NormalTok{ mod)) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_vline}\NormalTok{(}\DataTypeTok{xintercept =} \DecValTok{0}\NormalTok{, }\DataTypeTok{lty =} \DecValTok{2}\NormalTok{) }\OperatorTok{+}
\StringTok{  }\KeywordTok{geom_segment}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ low, }\DataTypeTok{xend =}\NormalTok{ high, }\DataTypeTok{y =}\NormalTok{ x2, }\DataTypeTok{yend =}\NormalTok{ x2, }\DataTypeTok{col =}\NormalTok{ mod)) }\OperatorTok{+}
\StringTok{  }\KeywordTok{labs}\NormalTok{(}\DataTypeTok{x =} \StringTok{"random slope"}\NormalTok{, }\DataTypeTok{y =} \StringTok{""}\NormalTok{) }\OperatorTok{+}\StringTok{ }
\StringTok{  }\KeywordTok{scale_y_continuous}\NormalTok{(}\DataTypeTok{breaks =} \DecValTok{1}\OperatorTok{:}\DecValTok{15}\NormalTok{, }\DataTypeTok{labels =} \KeywordTok{paste0}\NormalTok{(}\StringTok{"y"}\NormalTok{, }\KeywordTok{sprintf}\NormalTok{(}\StringTok{"%02d"}\NormalTok{, }\KeywordTok{rep}\NormalTok{(}\DecValTok{1}\OperatorTok{:}\KeywordTok{length}\NormalTok{(}\KeywordTok{unique}\NormalTok{(mydata}\OperatorTok{$}\NormalTok{group))))))}
\end{Highlighting}
\end{Shaded}

\includegraphics{8GLMM_files/figure-latex/unnamed-chunk-24-3.pdf}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{fittedvalues <-}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(}\DataTypeTok{mean =} 
                             \KeywordTok{c}\NormalTok{(}\KeywordTok{fitted.values}\NormalTok{(fit_glmer),}
                               \KeywordTok{fitted.values}\NormalTok{(fit_glmmTMB),}
\NormalTok{                               fit_inla}\OperatorTok{$}\NormalTok{summary.fitted.values}\OperatorTok{$}\NormalTok{mean),}
                           \DataTypeTok{mod =} \KeywordTok{rep}\NormalTok{(}\KeywordTok{c}\NormalTok{(}\StringTok{"glmer"}\NormalTok{, }\StringTok{"glmmTMB"}\NormalTok{, }\StringTok{"INLA"}\NormalTok{), }\DataTypeTok{each =} \KeywordTok{nrow}\NormalTok{(mydata)),}
                           \DataTypeTok{x =} \KeywordTok{rep}\NormalTok{(mydata}\OperatorTok{$}\NormalTok{x, }\DecValTok{3}\NormalTok{),}
                           \DataTypeTok{group =} \KeywordTok{rep}\NormalTok{(mydata}\OperatorTok{$}\NormalTok{group,}\DecValTok{3}\NormalTok{))}

\KeywordTok{ggplot}\NormalTok{(fittedvalues) }\OperatorTok{+}\StringTok{ }\KeywordTok{geom_count}\NormalTok{(}\DataTypeTok{data =}\NormalTok{ mydata, }\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ y}\OperatorTok{*}\FloatTok{0.8}\OperatorTok{+}\FloatTok{0.1}\NormalTok{, }\DataTypeTok{y =}\NormalTok{ y}\OperatorTok{*}\FloatTok{0.8}\OperatorTok{+}\FloatTok{0.1}\NormalTok{), }\DataTypeTok{col =} \KeywordTok{grey}\NormalTok{(}\FloatTok{0.4}\NormalTok{)) }\OperatorTok{+}\StringTok{ }\KeywordTok{scale_size_area}\NormalTok{() }\OperatorTok{+}
\StringTok{  }\KeywordTok{geom_point}\NormalTok{(}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{x =}\NormalTok{ x, }\DataTypeTok{y =}\NormalTok{ mean, }\DataTypeTok{col =}\NormalTok{ mod, }\DataTypeTok{pch =}\NormalTok{ mod, }\DataTypeTok{size =} \FloatTok{1.5}\NormalTok{)) }\OperatorTok{+}\StringTok{ }
\StringTok{  }\KeywordTok{facet_wrap}\NormalTok{(}\OperatorTok{~}\StringTok{ }\NormalTok{group) }\OperatorTok{+}\StringTok{ }\KeywordTok{labs}\NormalTok{(}\DataTypeTok{x =} \StringTok{"x"}\NormalTok{, }\DataTypeTok{y =} \KeywordTok{expression}\NormalTok{(}\KeywordTok{hat}\NormalTok{(y)))}
\end{Highlighting}
\end{Shaded}

\includegraphics{8GLMM_files/figure-latex/unnamed-chunk-24-4.pdf}

The results are quite similar (as approximation is necessary in the
computations, some differences are expected). We can make more
complicated models, and (especially from INLA) we also get much more
information about the model, but this is way out of the scope of this
course. Models like this require computer intenisve methods, and if you
are interested in learning more about this, you can take the course
TMA4300 - Computer Intensive Statistical Methods in the spring semester.

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

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"arm"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"reshape2"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"sp"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"glmmTBM"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"INLA"}\NormalTok{, }\DataTypeTok{repos =} \StringTok{"https://inla.r-inla-download.org/R/stable"}\NormalTok{, }
    \DataTypeTok{dep =} \OtherTok{TRUE}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"lme4"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"devtools"}\NormalTok{)}
\KeywordTok{library}\NormalTok{(devtools)}
\KeywordTok{install_github}\NormalTok{(}\StringTok{"romunov/AED"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"sjPlot"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"sjmisc"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

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

\begin{itemize}
\tightlist
\item
  \href{http://avesbiodiv.mncn.csic.es/estadistica/curso2011/regm26.pdf}{Bolker
  et al. (2008):\emph{Generalized linear mixed models: a practical guide
  for ecology and evolution}}
\item
  Zuur et al. (2009): ``Mixed Effects Models and Extensions in Ecology
  with R'', chapter 13 (pages 323-341). Available as free ebook from
  Springer for NTNU students. More explanations and less mathematics
  than Fahrmeir et al (2013), more focus on understanding.
  \href{https://link.springer.com/chapter/10.1007/978-0-387-87458-6_13}{Link
  to ebook Chapter 13}
\item
  Agresti (2015): Foundations of linear and generalized linear models,
  Chapter 9.
\end{itemize}


\end{document}
