heatherturner.net/talks/uRos2021   @HeathrTurnr

Overview

• Introduction to Generalized Nonlinear Models and the gnm package

• Case study of a structured main effect:

  Modelling the Consequences of Social Mobility with the Diagonal Reference Model

• Case study of a structured interaction:

  Modelling Mortality Trends with the Lee-Carter Model

heatherturner.net/talks/uRos2021   @HeathrTurnr

Linear models

$$\widehat{y_i}={\beta }_0+{\beta }_1{x}_i+{\beta }_2{z}_i+{ϵ}_i$$

$$E\left(y_i\right)={\beta }_0+\exp \left({\theta }_1\right)x_i+{\beta }_2{x}_i^2$$

In general

$$E\left(y_i\right)={\eta }_i\left(\beta \right)=\text{linear function of unknown parameters}$$

Also assumes variance essentially constant: $$\text{Var}\left(y_i\right)=\varphi a_i$$ with $a_i$ known (often $a_i\equiv 1$).

heatherturner.net/talks/uRos2021   @HeathrTurnr

Generalized linear models

Problems with linear models in many applications:

• range of $y$ is restricted (e.g., $y$ is a count, or is binary, or is a duration)
• effects are not additive
• variance depends on mean (e.g., large mean $\Rightarrow$ large variance)

Generalized linear models specify a non-linear link function and variance function to allow for such things, while maintaining the simple interpretation of linear models.

heatherturner.net/talks/uRos2021   @HeathrTurnr

GLM structure

A GLM is made up of a linear predictor

$$\eta ={\beta }_0+{\beta }_1{x}_1+...+{\beta }_p{x}_p$$

a link function that describes how $E\left(Y\right)=\mu$ depends on $\eta$

$$g\left(\mu \right)=\eta$$

a variance function that describes how $\text{Var}\left(Y\right)$ depends on $\mu$

$$\text{Var}\left(Y\right)=\varphi a_{i}V\left(\mu \right)$$

where the dispersion parameter $\varphi$ is a constant and $a_i$ are known values.

heatherturner.net/talks/uRos2021   @HeathrTurnr

Logistic regression

Suppose proportions $Y_i/n_i$ with $Y_i\sim \text{Binomial}(n_i,p_i)$. Then

$$\begin{array}{ccc}E\left(Y_i/n_i\right)& =p_i& \text{Var}\left(Y_i/n_i\right)=\frac{1}{n_i}{p}_i(1-p_i)\end{array}$$

The variance function is

$$V\left({\mu }_i\right)={\mu }_i(1-{\mu }_i)$$ Using a logit link maps the mean from $(0,1)$ to $(-\infty ,\infty )$ $$g\left({\mu }_i\right)=\text{logit}\left({\mu }_i\right)=\log \left(\frac{{\mu }_i}{1-{\mu }_i}\right)$$

heatherturner.net/talks/uRos2021   @HeathrTurnr

Log-linear model

Suppose counts $Y_i\sim \text{Poisson}\left({\lambda }_i\right)$. Then

$$\begin{array}{ccc}E\left(Y_i\right)& ={\lambda }_i& \text{Var}\left(Y_i\right)={\lambda }_i\end{array}$$

The variance function is

$$V\left({\mu }_i\right)={\mu }_i$$ Using a log link maps the mean from $(0,\infty )$ to $(-\infty ,\infty )$

$$g\left({\mu }_i\right)=\log \left({\mu }_i\right)$$

heatherturner.net/talks/uRos2021   @HeathrTurnr

Generalized Nonlinear Models

A generalized nonlinear model (GNM) is the same as a GLM except that we have

$$g\left(\mu \right)=\eta (x;\beta )$$

where $\eta (x;\beta )$ is nonlinear in the parameters $\beta$.

Equivalently an extension of nonlinear least squares model, where the variance of $Y$ is allowed to depend on the mean.

Using a nonlinear predictor can produce a more parsimonious and interpretable model.

heatherturner.net/talks/uRos2021   @HeathrTurnr

library(vcd)
mosaic(mentalHealth)

heatherturner.net/talks/uRos2021   @HeathrTurnr

Independence

A simple analysis of these data might be to test for independence of MHS and SES using a chi-squared test.

Equivalent to testing goodness-of-fit of the independence model $$\log \left({\mu }_{rc}\right)={\alpha }_r+{\beta }_c$$

Such a test compares the independence model to the saturated model

$$\log \left({\mu }_{rc}\right)={\alpha }_r+{\beta }_c+{\gamma }_{rc}$$ which may be over-complex.

heatherturner.net/talks/uRos2021   @HeathrTurnr

Row-column Association

One intermediate model is the Row-Column association model:

$$\log \left({\mu }_{rc}\right)={\alpha }_r+{\beta }_c+{\varphi }_r{\psi }_c$$ (Goodman, 1979), an example of a multiplicative interaction model.

For the mental health data:

Analysis of Deviance Table
Model 1: Freq ~ SES + MHS
Model 2: Freq ~ SES + MHS + Mult(SES, MHS)
Model 3: Freq ~ SES + MHS + SES:MHS
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1        15       47.4                     
2         8        3.6  7     43.8  2.3e-07
3         0        0.0  8      3.6     0.89

heatherturner.net/talks/uRos2021   @HeathrTurnr

gnm R package

• Provides a framework for specifying, fitting and criticizing GNMs in R
• Designed to be glm-like
  • common arguments, returned objects, methods, etc
  • equivalent results for linear predictors
• For family = gaussian equivalent to nls fit.

heatherturner.net/talks/uRos2021   @HeathrTurnr

Model Specification

• Linear terms in the model may be specified in the usual way, e.g.

y ∼ a + b + a:b

• Nonlinear terms must be specified using functions of class "nonlin"
  • specify structure of term, possible also labels & starting values
  • provided functions: Exp, Inv, Mult, MultHomog, Dref
  • custom functions

heatherturner.net/talks/uRos2021   @HeathrTurnr

RC model in gnm

RCmod <- gnm(Freq ~ SES + MHS + Mult(SES, MHS), family = poisson, 
             data = mentalHealth, verbose = FALSE, ofInterest = "Mult")
coef(RCmod)

Coefficients of interest:
       Mult(., MHS).SESA        Mult(., MHS).SESB        Mult(., MHS).SESC 
                 0.37314                  0.37649                  0.10062 
       Mult(., MHS).SESD        Mult(., MHS).SESE        Mult(., MHS).SESF 
                -0.04575                 -0.40728                 -0.70431 
    Mult(SES, .).MHSwell     Mult(SES, .).MHSmild Mult(SES, .).MHSmoderate 
                 0.76609                  0.07004                 -0.05557 
Mult(SES, .).MHSimpaired 
                -0.63370

heatherturner.net/talks/uRos2021   @HeathrTurnr

Parameterisation in (G)LMs

The independence model was defined in an over-parameterised form:

$$\begin{array}{cc}\log \left({\mu }_{rc}\right)& ={\alpha }_r+{\beta }_c\\ & =({\alpha }_r+1)+({\beta }_c-1)\\ & ={\alpha }_r^{\ast }+{\beta }_c^{\ast }\end{array}$$

glm will impose identifiability constraints, by default ${\alpha }_1={\beta }_1=0$

heatherturner.net/talks/uRos2021   @HeathrTurnr

Parameterisation in (G)NMs

In the Row-Column Association model

$$\begin{array}{cc}{\alpha }_r+{\beta }_r+{\varphi }_r{\psi }_c& ={\alpha }_r+({\beta }_r-{\psi }_c)+({\varphi }_r+1){\psi }_c\\ & ={\alpha }_r+{\beta }_r+\left(2{\varphi }_r\right)\left({\psi }_c/2\right)\end{array}$$

We need to constrain both location and scale for identifiability.

gnm works with overparameterised models

• will return one of infinite parameterisations at random
• by default, nonlinear parameters are not identifiable: standard error NA
• constrain parameters of interest, as required

heatherturner.net/talks/uRos2021   @HeathrTurnr

Background

Data from the ONS Longitudinal Survey, which links census and vital event data for 1% of the population of England and Wales.

Objective is to study the effect of social mobility on long-term limiting illness (LLTI)

a long-standing illness, health problem or handicap that limits a person's activities or the work they can do

Data include

• National Statistics socio-economic classification (7 classes, high man/prof to routine) in 1991 (origin) and 2001 (destination).
• Age, gender

heatherturner.net/talks/uRos2021   @HeathrTurnr

Bartley & Plewis Models

Bartley & Plewis (2007) modelled trajectories by combining social mobility effects with origin or destination effects.

For example, the origin + mobility model includes the term $${\theta }_{ij}=\{\begin{array}{cc}{\alpha }_i+{\delta }_0& i>j\\ {\alpha }_i+{\delta }_1& i=j\\ {\alpha }_i+{\delta }_2& i<j\end{array}$$

where ${\alpha }_i$ is the effect for origin $i$ and $i>j$ denotes moving to a more favourable destination class.

heatherturner.net/talks/uRos2021   @HeathrTurnr

GLM mobility models

The full models are logistic regression models for probability of LLTI in 2001.

Since the mobility models are GLMs, they can be fitted with glm, e.g.

dest_m <- glm(llti ~ age1991 + origin + mobility, 
              family = binomial, data = LS_m)

Age in 1991 was included as a covariate and men and women were modelled separately.

heatherturner.net/talks/uRos2021   @HeathrTurnr

Social Mobility Effects

The exponential of the mobility parameters gives the odds ratios of LLTI, e.g. for men

• given origin, odds of LLTI increased by downward mobility
• given destination, odds of LLTI decreased by downward mobility

heatherturner.net/talks/uRos2021   @HeathrTurnr

Weighted Residuals

The working residuals from the GLM fit can be averaged over each origin-destination combination to provide an indicator of fit for each trajectory

These average residuals can be standardized to be approximately N(0, 1) assuming the model is correct.

Passing a GLM fit to the mosaic function will produce a plot of the weighted residuals, with the p-value for a ${\chi }^2$ goodness-of-fit test.

heatherturner.net/talks/uRos2021   @HeathrTurnr

Although two sides of the same coin, the two models capture different features of the data.

heatherturner.net/talks/uRos2021   @HeathrTurnr

Diagonal Reference Model

Sobel (1991) proposed the diagonal reference model, which combines origin, destination and social mobility effects $$w_1{\gamma }_i+(1-w_1){\gamma }_j$$ The effect of moving from class $i$ to class $j$ is a weighted sum of the diagonal effects ${\gamma }_i$, where ${\gamma }_i$ is the effect for stable individuals in that class.

This is a alternative to adding main effects of origin and destination.

heatherturner.net/talks/uRos2021   @HeathrTurnr

Dref in gnm

A diagonal reference term may be specified via Dref(fac1, fac2, ...)

The weights are constrained to be non-negative and to sum to one by defining them as $$w_f=\frac{e^{{\delta }_f}}{{\sum }_i{e}^{{\delta }_i}}$$ and estimating unconstrained ${\delta }_f$ parameters.

heatherturner.net/talks/uRos2021   @HeathrTurnr

Following the previous analysis we fit the diagonal reference model with age as a covariate and fit models for men and women separately.

dref_m <- gnm(llti ~ age1991 + Dref(origin, dest),
              family = binomial, data = LS_m)
coef(summary(classMobility))

                                        Estimate Std. Error z value Pr(>|z|)
(Intercept)                             -1.30795         NA      NA       NA
age1991                                  0.06357   0.001221   52.08  < 2e-16
Dref(origin, dest)delta1                -0.05275         NA      NA       NA
Dref(origin, dest)delta2                -0.53378         NA      NA       NA
Dref(., .).origin|desthigh man/pro      -3.74444         NA      NA       NA
Dref(., .).origin|destlow man/pro       -3.29767         NA      NA       NA
Dref(., .).origin|destintermed          -2.74320         NA      NA       NA
Dref(., .).origin|destsmall empl/own ac -2.99956         NA      NA       NA
Dref(., .).origin|destlow sup/tech      -2.79887         NA      NA       NA
Dref(., .).origin|destsemi-routine      -2.47277         NA      NA       NA
Dref(., .).origin|destroutine           -2.27206         NA      NA       NA

heatherturner.net/talks/uRos2021   @HeathrTurnr

Diagonal Effects

The model for the stable individuals is an ordinary logistic regression

$$\text{logit}\left(p_{iik}\right)={\beta }_0+{\beta }_1{\text{age}}_k+{\gamma }_i$$ Setting ${\gamma }_1=0$, the diagonal effects are log odds ratios of LLTI in class $i$ against the first class for a given age

$$\begin{array}{cc}& \log \left(\frac{p_{iik}/(1-p_{iik})}{p_{11k}/(1-p_{11k})}\right)\\ =& ({\beta }_0+{\beta }_1{\text{age}}_k+{\gamma }_i)-({\beta }_0+{\beta }_1{\text{age}}_k)\\ =& {\gamma }_i\end{array}$$

heatherturner.net/talks/uRos2021   @HeathrTurnr

dref_m <- gnm(llti ~ age1991 + Dref(origin, dest),
              family = binomial, data = LS_m,  constrain = "(high man)")
coef(summary(classMobility))

                                        Estimate Std. Error z value Pr(>|z|)
(Intercept)                             -5.05239   0.063708 -79.305  < 2e-16
age1991                                  0.06357   0.001221  52.084  < 2e-16
Dref(origin, dest)delta1                      NA         NA      NA       NA
Dref(origin, dest)delta2                      NA         NA      NA       NA
Dref(., .).origin|desthigh man/pro       0.00000         NA      NA       NA
Dref(., .).origin|destlow man/pro        0.44677   0.053225   8.394  < 2e-16
Dref(., .).origin|destintermed           1.00124   0.063136  15.858  < 2e-16
Dref(., .).origin|destsmall empl/own ac  0.74488   0.052400  14.227  < 2e-16
Dref(., .).origin|destlow sup/tech       0.94557   0.053376  17.715  < 2e-16
Dref(., .).origin|destsemi-routine       1.27167   0.053763  23.653  < 2e-16
Dref(., .).origin|destroutine            1.47237   0.049871  29.524  < 2e-16

heatherturner.net/talks/uRos2021   @HeathrTurnr

Health Inequality

heatherturner.net/talks/uRos2021   @HeathrTurnr

DrefWeights

Setting ${\delta }_1=0$ is enough to constrain the weights, so that

$$\begin{array}{cccc}{w}_1& =\frac{1}{1+e^{{\delta }_2}}& w_2& =\frac{e^{{\delta }_2}}{1+e^{{\delta }_2}}\end{array}$$ DrefWeights() computes these weights and their standard errors, e.g.

 DrefWeights(dref_m)

$origin
 weight      se
0.61799 0.02649
$dest
 weight       se
0.38201 0.02649

heatherturner.net/talks/uRos2021   @HeathrTurnr

Diluting Effect of Social Mobility

The ratio of origin weight to destination weight quantifies the diluting effect of social mobility on health inequality

• 1:0 social mobility has no effect on individual
• 0:1 social mobility has no effect on inequality
• otherwise social mobility increases P(LLTI) in the upper classes and decreases P(LLTI) in the lower classes.

The larger the origin weight, the greater the diluting effect.

heatherturner.net/talks/uRos2021   @HeathrTurnr

Diagonal Weights

The diagonal weights for the LLTI models are

Since their destination class is given more weight, social mobility has a greater impact for women.

heatherturner.net/talks/uRos2021   @HeathrTurnr

Model Comparison

The models can be compared by the difference in deviance from the null model:

The diagonal reference model reduces the deviance the most despite requiring fewer degrees of freedom.

heatherturner.net/talks/uRos2021   @HeathrTurnr

Background

Data from the Human Mortality Database, which provides open, international access to population and mortality data from 41 countries.

Here we consider data on deaths for men aged 20 to 99 in Canada from 1921 to 2003.

The data comprise the variables

• Year: year death recorded
• Age: age at time of death
• Deaths: number of deaths
• Exposure: number at risk model

heatherturner.net/talks/uRos2021   @HeathrTurnr

Lee-Carter model

Lee & Carter (1992) proposed a model for age-specific population mortality rates that has been the basis of many subsequent analyses.

Suppose that death count $D_{ay}$ for individuals of age $a$ in year $y$ has mean ${\mu }_{ay}$ and quasi-Poisson variance $\varphi {\mu }_{ay}$.

Lee-Carter model

$$\begin{array}{cc}& \log \left({\mu }_{ay}/e_{ay}\right)={\alpha }_a+{\beta }_a{\gamma }_y,\\ \Rightarrow & \log \left({\mu }_{ay}\right)=\log \left(e_{ay}\right)+{\alpha }_a+{\beta }_a{\gamma }_y\end{array}$$

where $e_{ay}$ is the exposure (number of lives at risk).

heatherturner.net/talks/uRos2021   @HeathrTurnr

Lee-Carter in gnm

The Lee-Carter model can be fitted with gnm as follows

LCmodel <- gnm(Deaths ~ Mult(Exp(Age), Year),
               eliminate = Age, offset = log(Exposure),
               family = "quasipoisson")

• Exp(Age) is used to constrain the parameters representing the sensitivity of age group $a$ to have the same sign,
• Age is eliminated since it will typically have many levels,
• offset is used to add the log exposure with a coefficient of 1
• the quasipoisson family is used so that the dispersion parameter $\varphi$ is estimated rather than fixed to 1.

heatherturner.net/talks/uRos2021   @HeathrTurnr

Residuals are more dispersed for ages 25-35 and 50-65, as well as pre-1930

heatherturner.net/talks/uRos2021   @HeathrTurnr

Residuals by birth cohort (calendar year - age) show a clear trend

heatherturner.net/talks/uRos2021   @HeathrTurnr

Renshaw & Haberman extension

Renshaw & Haberman (2006) considered a cohort-based extension

$$\log \left({\mu }_{ay}/e_{ay}\right)={\alpha }_a+{\beta }_a^{\left(0\right)}{\gamma }_y+{\beta }_a^{\left(1\right)}{\lambda }_{y-a}$$

This is easy to specify in gnm

Deaths ~ Mult(Exp(Age), Year) + Mult(Exp(Age), Cohort)

But not so easy to fit!

• Year, Age and Cohort are directly related
• now have 405 multiplicative parameters, data for some cohorts is sparse

heatherturner.net/talks/uRos2021   @HeathrTurnr

Fitting the Renshaw & Haberman model

Step 1. Initially fit the model with age multipliers fixed to 1

$$\log \left({\mu }_{ay}/e_{ay}\right)={\alpha }_a+{\gamma }_y+{\lambda }_{y-a}$$

Step 2. Fit the full model with starting values from step 1, without constraining the age multipliers.

• Converges in 54 iterations
• Find the multiplier for age 99 in the first multiplicative term is negative

Step 3. Re-fit using Exp(Age) to constrain the age multipliers to be positive, setting $${\beta }_{99}^{\left(0\right)}=\exp (-1000)=0$$

heatherturner.net/talks/uRos2021   @HeathrTurnr

Residuals from Renshaw & Haberman model

heatherturner.net/talks/uRos2021   @HeathrTurnr

heatherturner.net/talks/uRos2021   @HeathrTurnr

Age multipliers

An advantage of re-parameterising the model such that $$\begin{array}{cccc}{\beta }_a^{\left(0\right)}& =\exp \left({\delta }_a^{\left(0\right)}\right)& {\beta }_a^{\left(1\right)}& =\exp \left({\delta }_a^{\left(1\right)}\right)\end{array}$$ is that simple contrasts (differences) of ${\delta }_a^{\left(n\right)}$ are identifiable.

We can obtain the contrasts and standard errors with getContrasts()

AgeMult2 <- pickCoef(APCmodel_GNM,  "Mult(Exp(.), Cohort)", fixed = TRUE)
getContrasts(APCmodel_GNM, AgeMult2, check = FALSE)$qvframe[1:3,]

                           estimate     SE quasiSE quasiVar
Mult(Exp(.), Cohort).Age20   0.0000 0.0000  0.0810 0.006561
Mult(Exp(.), Cohort).Age21  -0.2361 0.1294  0.1027 0.010543
Mult(Exp(.), Cohort).Age22  -0.1666 0.1290  0.1021 0.010420

heatherturner.net/talks/uRos2021   @HeathrTurnr

heatherturner.net/talks/uRos2021   @HeathrTurnr

Spline model

Define a piecewise linear spline over age:

Age_num <- as.numeric(as.character(Canada$Age))
Age_bs <- bs(Age_num, knots = c(seq(25, 70, by = 5), 80), 
             degree = 1)

Use instead of age in the Renshaw & Haberman model:

Deaths ~ Mult(Exp(Age_bs), Year) + Mult(Exp(Age_bs), Cohort)

heatherturner.net/talks/uRos2021   @HeathrTurnr

Model comparison

The spline model preserves the year and cohort multipliers which could be modelled using time series methods for forecasting as in Renshaw & Haberman (2006) - only depends on the differences between multipliers.

heatherturner.net/talks/uRos2021   @HeathrTurnr

Finding out more

Case studies with reproducibility materials: github.com/hturner/gnm-case-studies

• So far Lee-Carter and UNIDIFF
• Diagonal Reference and Rasch models to be added
• Working paper with D. Firth and I. Kosmidis in progress

See the gnm vignette for further examples

Also github.com/hturner/gnm-half-day-course and github.com/hturner/gnm-day-course

heatherturner.net/talks/uRos2021   @HeathrTurnr

References

Agresti, A. 2002. Categorical Data Analysis. 2nd ed. New York: Wiley.

Bartley, M, and I Plewis. 2007. 10.1111/j.1467-985X.2006.00464.x.

Goodman, L A. 1979. J. Amer. Statist. Assoc. 10.1080/01621459.1979.10481650.

Lee, R D, and L Carter. 1992. J. Amer. Statist. Assoc. 10.1080/01621459.1992.10475265.

Renshaw, A, and S Haberman. 2006. Insur Math Econ 10.1016/j.insmatheco.2005.12.001.

Sobel, M. E. 1981. Amer. Soc. Rev. 10.2307/2095086.