---
title: "TMA4315 Generalized linear models H2018"
subtitle: "Module 6: Categorical regression"
author: "Mette Langaas, Department of Mathematical Sciences, NTNU, with contibutions from Ingeborg G. Hem"
date: "25.10.2017 [PL], 26.10.2017 [IL]"
output: #3rd letter intentation hierarchy
  # html_document:
  #   toc: true
  #   toc_float: true
  #   toc_depth: 2
 # pdf_document:
 #   toc: true
 #   toc_depth: 2
 #   keep_tex: yes
  beamer_presentation:
    keep_tex: yes
    fig_caption: false
    latex_engine: xelatex
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,tidy=TRUE,message=FALSE,warning=FALSE,results="hold")
showsol<-FALSE
```

(Latest changes: 24.10: added IL + corrected typos, 23.10.2018: first version)

# Overview

## Learning material
This topic is _new_ on the reading list this year.

* Textbook: Fahrmeir et al (2013): Chapter 6 (not p 344-345 nominal models and latent utility models, not 6.3.2 Sequential model, and  not category specific varables on page 344-3458).
* [Classnotes 25.10.2018](https://www.math.ntnu.no/emner/TMA4315/2018h/M62018.pdf) 

---

## Topics

* multinomial random component 
* nominal vs. ordinal response
* ungrouped and grouped data
* multivariate exponential family

* nominal response and logit models

* ordinal reponse and logit models - based on a latent model

* likelihood inference

<!-- * conditional poisson data -->
<!-- * chisquare test for contingency tables -->

Jump to [interactive](#interactive).

---

# Categorical random component 

We consider a situation where our random variable (response) is given as one of $c+1$ possible categories (where we will look at category $c+1$ as the reference category).

The categories will either be 

* Unordered: _nominal response variable_. Example: food types in alligator example.
* Ordered: _ordered response variable_. Example: degrees of mental impairment.

---

Assumptions:

* _Independent_ observation pairs $({\bf Y}_i,{\bf x}_i)$.
* $\pi_{ir}$: probability that the response is category $r$ for subject $i$. 
* $\sum_{s=1}^{c+1}\pi_{is}=1$ for all $i$, so that $\pi_{i,c+1}=1-\sum_{s=1}^{c}\pi_{is}$. So, we have $c$ probabilities to estimate.
* Further, the covariate vector ${\bf x}_i$ consists of the same measurements for each response category (that is, not different covariate types that are measured for each response category - which in our textbook is written as _independent of the response category_).


---

When coding the response variable we use a dummy variable coding with $c$ elements (the $c+1$ category is the reference level). This means that if we have that $\pi_{ir}=1$ then ${\bf y}_i=(0,0,\ldots,0,1,0,\ldots,0)$ with a value of $1$ in the $r$th element of ${\bf y}_i$. If observation $i$ comes from category $c+1$ we have ${\bf y}_i=(0,0,\ldots,0)$.

---

## Categorical regression

is modelling and estimating the probabilites $\pi_{ir}=P(Y_i=r)=P(Y_{ir}=1)$ as a function of the covariates ${\bf x}_i$.
The modelling is done differently for nominal (unordered) and ordered categories, but both rely upon the multinomial distribution.

For unordered categories, a Poisson distribution can also be used

---

## The multinomial distribution

Probability mass function for one observation:

$$f({\bf y})=\pi_1^{y_{1}} \pi_2^{y_{2}} \cdots \pi_c^{y_c} (1-\pi_1-\pi_2-\cdots-\pi_c)^{1-y_1-y_2-\cdots-y_c}$$
where then ${\bf y}=(y_1,y_2,\ldots,y_c)$ and $y_r=1$ if the observation comes from the $r$th category.

---

If we then have $m$ independent trials then ${\bf y}=(y_1,y_2,\ldots,y_c)$ is summed over our $m$ responses, so that $y_r$ is the number of observations where the response is from the $r$th category.

$$f({\bf y})=\frac{m!}{y_1!\cdots y_c! (m-y_1-\cdots -y_c)!}\pi_1^{y_{1}} \pi_2^{y_{2}} \cdots \pi_c^{y_c}(1-\pi_1-\pi_2-\cdots-\pi_c)^{m-y_1-y_2-\cdots-y_c}$$

---

The mean and the covariance matrix of the random vector ${\bf Y}$ are given by:

$$\text{E}({\bf Y})=m\boldsymbol{\pi}=\begin{pmatrix}m\pi_1\\m\pi_2\\ \vdots \\ m\pi_c\end{pmatrix}$$
$$\text{Cov}({\bf Y})=m \begin{pmatrix}\pi_1(1-\pi_1)& -\pi_1\pi_2&\cdots &-\pi_1\pi_c\\
-\pi_2\pi_1& \pi_2(1-\pi_2)&\cdots &-\pi_2\pi_c\\
\vdots & \vdots & \ddots & \vdots\\
-\pi_c\pi_1& -\pi_c\pi_2&\cdots & \pi_c(1-\pi_c)\end{pmatrix}$$

**Q**: what about $\text{E}(Y_{c+1})$ and $\text{Cov}(Y_1,Y_{c+1})$?

---

Finally, if we look at $\bar{Y}_r=\frac{1}{m}Y_r$ then $\bar{\bf Y}=\frac{1}{m}{\bf Y}$ follows a scaled multinomial distribution $\bar{\bf Y}\sim \frac{1}{m}M(m,{\boldsymbol \pi})$ with $\text{E}(\bar{\bf Y})={\boldsymbol \pi}$ and $\text{Cov}(\bar{\bf Y})=\frac{1}{m^2}\text{Cov}({\bf Y})$.

# Data

### Ungrouped data

$${\bf Y}=\begin{pmatrix} Y_{11}& Y_{12} & \cdots & Y_{1c}\\
Y_{21}& Y_{22} & \cdots & Y_{2c}\\
\vdots & \vdots & \ddots & \vdots\\
Y_{n1}& Y_{n2} & \cdots & Y_{nc}\\
\end{pmatrix}$$

and ${\bf X}$ is an $n \times p$ matrix as usual.


---

### Grouped data

As for the binomial case we look at the number of occurences with a group - that is, one covariate pattern. 

$${\bf Y}=\begin{pmatrix} Y_{11}& Y_{12} & \cdots & Y_{1c}\\
Y_{21}& Y_{22} & \cdots & Y_{2c}\\
\vdots & \vdots & \ddots & \vdots\\
Y_{G1}& Y_{G2} & \cdots & Y_{Gc}\\
\end{pmatrix}$$

The notation here is that we have $n_i$ observation for each covariate pattern (group) $i$ for $i=1,\ldots,G$. This will replace the $m$ used for the multinomial distribution above.

---

# Regression with nominal responses

nominal=unordered

Agresti (2015, p203): "The model treats the response variable as nominal scale in the following sense: if the model holds and the outcome categories are permuted in any way, the model still holds with the corresponding permuatation of the effects."

---

This is a generalization of the binary logit model with $P(Y=1)$ vs $P(Y=0)$, to $c$ models of $\pi_{ir}$ vs $\pi_{i,c+1}$ for $r=1,\ldots,c$.

The models can be written using log ratios:
$$ \ln(\frac{\pi_{ir}}{\pi_{i,c+1}})={\bf x}_i^T {\boldsymbol \beta}_r$$

Remark: ${\boldsymbol \beta}_r$ is the $p\times 1$ coefficient vector for the $r$th response 

Using this we may also look at the log ratio for any two probabilites $\pi_{ia}$ and $\pi_{ib}$:

$$\ln(\frac{\pi_{ia}}{\pi_{ib}})=\ln(\frac{\pi_{ia}}{\pi_{i,c+1}})-\ln(\frac{\pi_{ib}}{\pi_{i,c+1}})={\bf x}_i^T ({\boldsymbol \beta}_a-{\boldsymbol \beta}_b)$$

---

Alternatively, we may write out the model for the probabilites:

$$P(Y_i=r)=\pi_{ir}=\frac{\exp({\bf x}_i^T {\boldsymbol\beta}_r)}{1+\sum_{s=1}^{c}\exp({\bf x}_i^T {\boldsymbol\beta}_s)}$$

$$P(Y_i=c+1)=\pi_{i,c+1}=1-\pi_{i1}-\cdots \pi_{ic}=\frac{1}{1+\sum_{s=1}^{c}\exp({\bf x}_i^T {\boldsymbol\beta}_s)}$$

---

## Multivariate GLM

This is a multivariate GLM and the multinomial distribution is _a multivariate exponential family_.

$$f({\bf y}_i,{\boldsymbol \theta}_i,\phi)=\exp(\frac{{\bf y}_i^T{\boldsymbol \theta}_i-b({\boldsymbol \theta}_i)w_i}{\phi}+c({\bf y}_i,\phi,w_i))$$

where ${\boldsymbol \theta}$ has dimension $c$.

---

### Multivariate GLM-set-up

1. ${\bf Y}_i$ is multinomial with 

$${\boldsymbol \mu}_i=\text{E}({\bf Y}_i)=
{\boldsymbol \pi}_i=\begin{pmatrix} \pi_{i1}\\ \pi_{i2}\\ \vdots \\ \pi_{i,c+1}\end{pmatrix}$$
Remark: if grouped data we instead look at $\bar{\bf Y}_i\sim \frac{1}{n_i}M(n_i,\pi_i)$ so that the mean is ${\boldsymbol \pi}_i$

---

2. Linear predictor is now a $c \times 1$ vector: 

$${\boldsymbol \eta}_i=\begin{pmatrix} \eta_{i1}\\ \eta_{i2}\\ \vdots \\ \eta_{i,c}\end{pmatrix}=
\begin{pmatrix} {\bf x}_i^T{\boldsymbol \beta}_{1}\\ {\bf x}_i^T{\boldsymbol \beta}_{2}\\ \vdots \\ {\bf x}_i^T{\boldsymbol \beta}_{c} \end{pmatrix}$$

---

3. Link functions ($c$ of those): ${\bf g}(\boldsymbol{\mu}_i)={\boldsymbol \eta}_i$ where for the _nominal logit data model_ element $r$ (for $r=1,\ldots, c$) of ${\bf g}$ is

$$ g_{r}(\boldsymbol{\mu}_i)=\ln(\frac{\mu_{ir}}{1-\mu_{i1}-\cdots-\mu_{ic}})=
\ln(\frac{\pi_{ir}}{1-\pi_{i1}-\cdots-\pi_{ic}})$$

We also define response functions $(\bf h)$ with elements $h_r$ given by
$\pi_{ir}=h_r(\eta_{i1},\eta_{i2},\ldots,\eta_{ic})$, and we have for the _nominal data model_

$$\pi_{ir}=h_r((\eta_{i1},\eta_{i2},\ldots,\eta_{ic})=\frac{\exp(\eta_{ir})}{1+\sum_{s=1}^c \exp(\eta_{is})}$$



---

It turns out that the reference category logits are the canonical links for the multinomial distribution GLM. 

In this case, as for the univariate exponential family GLM the loglikelihood is concave with an unique maximum (if it exists) and the expected and observed Fisher information matrices are equal. 

As before, we find maximum likelihood parameter estimates from the Fisher scoring or Newton Raphson method.

Remember: now we have $p\times c$ parameters to estimate — $p$ for each category $c$. All of these coefficients may either be put into a long vector (length $p\cdot c$) — which might be easiest to understand for the estimation, or into a matrix of dimension $p \times c$ — might be easier for viewing.

---

### Likelihood
(grouped data)

With the notation that ${\boldsymbol \beta}$ is a long vector with the coefficients for the $c$ categories stacked upon eachother. 


$$L({\boldsymbol \beta})=\Pi_{i=1}^G f({\bf y}_i\mid {\boldsymbol \pi})$$
where $f$ is the multinomial distribution function.

---

### Loglikelihood

$$ l({\boldsymbol \beta})\propto\sum_{i=1}^G \sum_{s=1}^{c+1} y_{is}\ln(\pi_{is})$$

<!-- \sum_{i=1}^n (\sum_{j=1}^c y_{ij}\ln(\frac{\pi_{ij}}{\pi_{i,c+1}})+\ln\pi_{i,c+1})$$ -->
where we remember that $y_{i,c+1}=n_i-y_{i1}-\cdots-y_{ic}$, and $1-\pi_{i1}-\cdots \pi_{ic}$.

(This formula is also correct for the ordinal model of the next section.)
General formulas for the score function and expected Fisher information matrix follow later.

---

### Deviance

The derivation used for the binary GLM model generalizes directly ot the multinomial GLM. The fitted probabilities are $\hat\pi_{ij}$ (group $i$ and category $j$) and the saturated model (grouped data) is $n_i \tilde\pi_{ij}=y_{ij}$. 
$$D=2 \sum_{i=1}^G \sum_{s=1}^{c+1} y_{is}\ln(\frac{{y}_{is}}{n_i \hat{\pi}_{is}})$$

The asymptotic distribution is as before $\chi^2$ with "the number of groups times number of categories minus 1 (Gc)" minus "the number of covariates (cp)", giving $Gc-cp=c(G-p)$ degrees of freedom.

---

The deviance can be used for model check with grouped data ($G$ groups with $n_i$ observations), but can be used to compare nested unsaturated models also for individual (ungrouped) data, with again an asymptotic $\chi^2$ distribution with the difference of number of parameters between the two models.

This formula is also correct for the ordinal model of the next section, except that the number of parameters estimated differ.

---

## Alligators example 
Example and data are taken from Agresti (2015, pages 217-219).

Research question: what is the factors influencing the primary food choice of alligators?

Data are from 219 captured alligators from four lakes in Florida, where the stomack contents of the alligators were investigated. The weight of different types of food was measured, and then the primary food choice (highest weight) was noted. The primary choice is given as y1:y5 below. In addition the size of the alligator (non-adult or adult) was registered.

---

* lake: each of the 4 lakes in Florida (1:4)
* size: non-adult=the size of the alligator (0: 2.3 meters or smaller) and adult=(1: larger than 2.3 meters)
* y1: fish
* y2: inverterbrate
* y3: reptile
* y4: bird
* y5: other

These data are grouped, and we let y1:fish be the reference category.

---

\small

```{r}
# data from Agresti (2015), section 6, with use of the VGAM packages
data="http://www.stat.ufl.edu/~aa/glm/data/Alligators.dat"
ali = read.table(data, header = T)
ali
attach(ali)
```
---
  
  
```{r}
y.data=cbind(y2,y3,y4,y5,y1)
y.data
dim(y.data)
x.data=model.matrix(~size+factor(lake),data=ali)
x.data
dim(x.data)
```


---

```{r}
# We use library VGAM:
library(VGAM)

# We fit a multinomial logit model with fish (y1) as the 
#  reference category:
fit.main=vglm(cbind(y2,y3,y4,y5,y1)~
                size+factor(lake),
         family=multinomial, data=ali)
# summary(fit.main)
pchisq(deviance(fit.main),df.residual(fit.main),lower.tail=FALSE)
```

---

\normalsize

**Q**: 

* Why is the number of degrees of freedom for the residual deviance 12? Hint: there are 8 covariate patterns, and we have 5 reponse categories.
* How can you interpret the coefficient for inverterbrate (y2) and size? Hint: we have y2,y3,y4,y5 as 1:4.

\small

```{r}
exp(coefficients(fit.main))
```

\normalsize

---

Testing out other models, and comparing with LRT-test - by using deviances for different models.

```{r}
# Fit model with only lake:
fit.lake=vglm(cbind(y2,y3,y4,y5,y1)~factor(lake),
              family=multinomial, data=ali)
# Test effect of size (no anova command is available)
(G2=deviance(fit.lake)-deviance(fit.main))
(df.diff=df.residual(fit.lake)-df.residual(fit.main))
1-pchisq(G2,df.diff)
# Size has a significant effect
```

---

```{r}
# Fit model with only size:
fit.size=vglm(cbind(y2,y3,y4,y5,y1)~size,family=multinomial, data=ali)

# Test effect of lake 
(G2=deviance(fit.size)-deviance(fit.main))
(df.diff=df.residual(fit.size)-df.residual(fit.main))
1-pchisq(G2,df.diff)
# Lake has a significant effect
```

---

**Q**: explain what is presented below, in particular "what is the probability that the main food source is fish given size=0 and lake=1"?

\small

```{r,eval=TRUE}
library(knitr)
# Fitted values for main effect model "fit.main":
fitted=data.frame(fitted(fit.main),lake=ali$lake,size=ali$size)
kable(fitted)
```

\normalsize

```{r,echo=FALSE,eval=FALSE}

ggplot(fitted, aes(x = lake, y = probability, colour = ses)) + geom_line() + facet_grid(variable ~., scales = "free")
```

\normalsize

---

# Regression with ordinal responses
(we will only consider cumulative models - and not sequential models)

An unobservable latent variable $U_i$ drives the observed category $Y_i$.

$$ Y_i=r \Leftrightarrow \theta_{r-1}\le U_i \le \theta_r$$
where these $\theta$s are our unobservable thresholds, and the thresholds are monotonely increasing, $-\infty=\theta_0< \theta_1 < \cdots < \theta_{c+1}=\infty$.

We further assume that the latent variables are dependent on our covariates through

$$ U_i=-{\bf x}_i^T{\boldsymbol \beta} + \varepsilon_i$$

where we have a new random variable that has cumulative distribution function (cdf) $F$. No intercept is included due to identifiability issuse (shift in intercept would produce the same effect as negative shift in threshold).

---

We get rid of the latent variable $U_i$ by considereing

$$P(Y_i\le r)=P(U_i\le \theta_r)=P(-{\bf x}_i^T {\boldsymbol \beta}+\varepsilon_i \le \theta_r)$$
$$=P(\varepsilon_i \le \theta_r+{\bf x}_i^T {\boldsymbol \beta})=F(\theta_r+{\bf x}_i^T {\boldsymbol \beta}) $$

Observe that the final expression does not include the latent variable $U_i$, but includes the unknown threshold and $k$ regression parameters.

---

Different choices of $F$ will give different models, and we will only consider $F$ to be the cdf for the logistic distribution. (Another popular choice is the cdf of the standard normal distribution.)

$$P(Y_i\le r)=\frac{\exp(\theta_r+{\bf x}_i^T {\boldsymbol \beta})}{1+\exp(\theta_r+{\bf x}_i^T {\boldsymbol \beta})}$$

which also can be written as

$$ \ln(\frac{P(Y_i\le r)}{P(Y_i> r)})=\theta_r+{\bf x}_i^T {\boldsymbol \beta}$$

---

Our model is a proportional odds model, in the sense that the cumulative odds are proportional across categories

$$\frac{\frac{P(Y\le r\mid {\bf x}_i)}{P(Y> r\mid {\bf x}_i)}}
{\frac{P(Y\le r\mid {\bf x}^*_i)}{P(Y> r\mid {\bf x}^*_i)}}= \exp(({\bf x}_i-{\bf x}^*_i)^T {\boldsymbol \beta})$$

Observe that this is independent of $r$.

---

### Response function

What is the response function here?

$$\pi_{i1}=F(\eta_{i1})$$
$$\pi_{ir}=F(\eta_{ir})-F(\eta_{i,r-1})$$
<!-- =h_r(\eta_{i1},\ldots, \eta{i}) -->

where $\eta_{ir}=\theta_r+{\bf x}_i^T{\boldsymbol \beta}$, and $F$ is the logistic cdf.

---

### Plotting the cumulative probabilites

We have five categories, where the fifth is the reference category. 

True parameters

* $\theta_1=0$, $\theta_2=1$, $\theta_3=4$ and $\theta_4=6$ and 
* one covariate with parameter $\beta=1$. 

The graph shows the cumulative probability $P(Y\le r)$ for $r=1$ (red), $r=2$ (organge), $r=3$ (green), $r=4$ (purple). 

Observe the parallell lines.

What would $P(Y\le 5)$ be? Why is this missing from the plot?

---

```{r,echo=showsol}
library(ggplot2)
#print(c(beta0, beta1))
xrange=c(-8,8)
beta0s=c(0,1,4,6)
beta1=1
etas1=beta1*seq(xrange[1],xrange[2],length=100)
#etas=data.frame(cbind(beta0s[1]+etas1,beta0s[2]+etas1,beta0s[3]+etas1,beta0s[4]+etas1))
#colnames(etas)=paste("etas",1:4,sep="")
#etarange=range(c(etas))
ggplot(data.frame(etas1=etas1), aes(x=etas1))+
  xlab(expression(eta))+ ylab("")+theme_minimal()+scale_x_continuous(limits = c(-12,10))+
  stat_function(fun=function(etas1) exp(beta0s[1]+etas1)/(1+exp(beta0s[1]+etas1)), geom="line", colour="red")+
  stat_function(fun=function(etas1) exp(beta0s[2]+etas1)/(1+exp(beta0s[2]+etas1)), geom="line", colour="orange")+
  stat_function(fun=function(etas1) exp(beta0s[3]+etas1)/(1+exp(beta0s[3]+etas1)), geom="line", colour="green")+
  stat_function(fun=function(etas1) exp(beta0s[4]+etas1)/(1+exp(beta0s[4]+etas1)), geom="line", colour="purple")

```


---

## Mental health data example

Example and data are taken from Agresti (2015, pages 219-223).

Research question: understand mental health issues.

The data comes from a random sample of size 40 of adult residents of Alachua County, Florida, USA. 

* Mental impairment $Y$: 1=well, 2=mild symptom formation, 3=moderate symptom formation, 4=impaired. 
* Life event index ($x_1$): compsite measure of the number and severity of important life events within the last three years (birth, new job, divorce, death in the family, ...)
* SES ($x_2$): socioeconomic index, 1=high, 0=low. 

These data are ungrouped (but could be grouped). In the original study several other explanatory variables were studied.

---


```{r}
# Read mental health data from the web:
library(knitr)
data="http://www.stat.ufl.edu/~aa/glm/data/Mental.dat"
mental = read.table(data, header = T)
colnames(mental)
apply(mental,2,table)
#kable(mental)
```

---

```{r}
library(VGAM)
# We fit a cumulative logit model with main effects of "ses" and "life":
fit.imp=vglm(impair~life+ses,family=cumulative(parallel=T),data=mental)
# parallell=T gives proportional odds structure - only intercepts differ
summary(fit.imp)
```

---

The ML fit for this model can be written as

$$ \text{logit}(\hat{P}(y_i\le r))=\hat{\theta}_r+0.319 x_{i1}+1.111 x_{i2}$$

**Q**: give an interpretation of this model!

Remember:

* Life event index ($x_1$): compsite measure of the number and severity of important life events within the last three years (birth, new job, divorce, death in the family, ...)
* SES ($x_2$): socioeconomic index, 1=high, 0=low. 

---

**Q**: How can you interpret the last line below? Why is it exp(CI(beta)) and not CI(exp(beta))?

```{r}
exp(confint(fit.imp))
```

---

**Q**: How are these predictions calculated? What is the interpretation?

```{r,eval=TRUE}
fitted=data.frame(fitted(fit.imp),ses=mental$ses, life=mental$life)
fitted[c(6,18,10),] #0,7 not fitted

xs=cbind(c(2,7,2,7),c(0,0,1,1))
coeff=coefficients(fit.imp)
linpreds=cbind(coeff[1]+xs%*%coeff[4:5],
      coeff[2]+xs%*%coeff[4:5],
      coeff[3]+xs%*%coeff[4:5])
(cprobs=exp(linpreds)/(1+exp(linpreds)))
(pprobs=cbind(cprobs[,1],cprobs[,2]-cprobs[,1],cprobs[,3]-cprobs[,2],1-cprobs[,3]))

```

---

**Q**: What do you see here, and what is the formula for this matrix?

```{r}
vcov(fit.imp)
```

---

```{r}
# We consider a model with interaction between "ses" and "life":
fit.int=vglm(impair~life+ses+life:ses,family=cumulative(parallel=T),data=mental)
summary(fit.int)
```

---

```{r}
# And test if there is a significant effect of interaction:
G2=deviance(fit.imp)-deviance(fit.int)
df.diff=df.residual(fit.imp)-df.residual(fit.int)
1-pchisq(G2,df.diff)
# The effect of interaction is not significant
```

```{r}
# We consider a model where the effect of the covariates may 
# differ between the cumulative logits - so not parallell lines for the cdfs
fit.nopar=vglm(impair~life+ses,family=cumulative,data=mental)
summary(fit.nopar)

# The change in the deviance compared to the model "fit.imp" 
# is 99.0979-96.7486=2.3493 with df.diff=115-111=4, which is not significant

# So model "fit.imp"seems fine.
```

---

## Why not use MLR instead of ordinal regression?
Based on Agresti (2015, p 214-216)

To use MLR the ordinal categories need to be replaced with numerical values, and we then need to assume a normal error structure. The following are questions to be answered and possible limitation to be assumed for using MLR instead of ordinal regression:

* how to translate ordered categories into numerical scores?
* is it better with an ordinal variable with some range than a single numerical number?
* MLR will not give probabilities for each response category
* variability in the response may be dependent on the category, for MLR we assume homoscedasticity


# Likelihood inference

We use the notation that ${\boldsymbol \beta}$ is a long vector with all regression parameters.
The content of this vector is slightly different for our two models, with intercept and $k$ covariate effects for each response category for the nominal model - and with $c$ thresholds but the same $k$-dimensional ${\boldsymbol \beta}$ vector for all categories.

Full matrix versions (over all $i$) can be found in our textbook, page 345-346.

## Loglikelihood

We have seen that the loglikelihood is:
$$ l({\boldsymbol \beta})\propto\sum_{i=1}^n \sum_{s=1}^{c+1} y_{is}\ln(\pi_{is})$$

<!-- \sum_{i=1}^n (\sum_{j=1}^c y_{ij}\ln(\frac{\pi_{ij}}{\pi_{i,c+1}})+\ln\pi_{i,c+1})$$ -->
where we remember that $y_{i,c+1}=n_i-y_{i1}-\cdots-y_{ic}$, and $1-\pi_{i1}-\cdots \pi_{ic}$.

---

## Design matrix and coefficient vector
The design matrix ${\bf X}$ and coefficient vector are different for our nominal logit model and our ordinal cumulative model.

### Nominal logit model

$${\bf X}_i=\text{diag}({\bf x}_i^T)=\begin{pmatrix}{\bf x}_i^T & 0 & \cdots & 0\\
0 & {\bf x}_i^T  & \cdots & 0\\
\vdots & \vdots & \ddots & \vdots\\
0 & 0 & \cdots & {\bf x}_i^T\\
\end{pmatrix}$$

where the $0$s are $1\times p$ vectors. The dimension of the design matrix for covariate pattern $i$ is $c \times c\cdot p$.

---

The vector of coefficients has dimension $c\cdot p \times 1$.
$${\boldsymbol \beta}=\begin{pmatrix}{\boldsymbol \beta}_1\\{\boldsymbol \beta}_2\\\vdots\\{\boldsymbol \beta}_c\\ \end{pmatrix}$$

---

### Ordinal cumulative model

$${\bf X}_i=\begin{pmatrix}1 & 0 & \cdots & 0 & {\bf x}_i^T\\
0 & 1 & \cdots & 0 & {\bf x}_i^T\\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & 1 & {\bf x}_i^T\\
\end{pmatrix}$$
he dimension of the design matrix for covariate pattern $i$ is $c \times (c+k)$

The vector of coefficients has dimension $(c+ k) \times 1$ (where $p=k+1$), and now the thresholds replace the intercept and are put first in the vector, and the effects of the covariates are the same for all categories. 
$${\boldsymbol \beta}=\begin{pmatrix}\theta_1\\ \theta_2 \\ \vdots \\ \theta_c \\{\boldsymbol \beta}\\ \end{pmatrix}$$

---

## Score function

$${\bf s}({\boldsymbol \beta})=\sum_{i=1}^G {\bf X}_i^T {\bf D}_i {\boldsymbol \Sigma}_i^{-1}({\bf y}_i-n_i{\boldsymbol \pi}_i)$$

where 

* ${\bf D}_i=\frac{\partial {\boldsymbol h}({\boldsymbol \eta})}{\partial {\boldsymbol \eta}}\lvert_{{\boldsymbol \eta}={\boldsymbol \eta}_i}$ has dimension $c \times c$
* ${\boldsymbol \Sigma}_i=\text{Cov}({\bf Y}_i)$ 

---

## Fisher information
The dimension of the matrix is $cp\times cp$ for the nominal case and $(c+k) \times (c+k)$ for the ordinal case studied.

$$F({\boldsymbol \beta})=\sum_{i=1}^G {\bf X}_i^T {\bf W}_i {\bf X}_i$$

where ${\bf W}_i$ is given as ${\bf D}_i{\boldsymbol \Sigma}_i^{-1}{\bf D}_i^T$.

---

## Finding the ML estimate

As in modules 1-5 we find the ML estimate by the Fisher scoring or Newton Raphson method. 

## Asymptotic distribution

As in modules 1-5 the ML estimator $\hat{\boldsymbol \beta}$ asymptotically follows a multivariate normal distribution with unbiased mean and covariance matrix given by the inverse of the expected Fisher information matrix.

---

# Summing up


---

# R packages

```{r, eval=FALSE}
install.packages(c("VGAM", 
  "ggplot2", 
  "statmod",
  "knitr"))
```

# Further reading

* A. Agresti (2015): "Foundations of Linear and Generalized Linear Models." Chapter 6. Wiley.
