Generalized linear mixed models

Authors

Douglas Bates

Phillip Alday

Published

2026-08-28

After working through this page you will be able to:

  • explain how GLMMs extend LMMs through a link function and a non-Gaussian conditional distribution;
  • contrast the role of the linear predictor in LMMs and GLMMs;
  • fit a logistic (Bernoulli) GLMM and summarize the results.
NoteBefore you start

Prerequisites: Analysis of the sleepstudy data; the Building Models section is helpful.

Datasets used: sleepstudy and contra (see the dataset catalog).

Load the packages to be used

Code 
using AlgebraOfGraphics
using CairoMakie
using DataFrameMacros
using DataFrames
using MixedModels
using MixedModelsMakie
using MixedModelsDatasets: dataset
using SMLP2026: fit_or_restore
using Statistics

const progress = isinteractive()

1 Matrix notation for the sleepstudy model

sleepstudy = DataFrame(dataset(:sleepstudy))
180×3 DataFrame
155 rows omitted
Row subj days reaction
String Int8 Float32
1 S308 0 249.56
2 S308 1 258.705
3 S308 2 250.801
4 S308 3 321.44
5 S308 4 356.852
6 S308 5 414.69
7 S308 6 382.204
8 S308 7 290.149
9 S308 8 430.585
10 S308 9 466.353
11 S309 0 222.734
12 S309 1 205.266
13 S309 2 202.978
⋮ ⋮ ⋮ ⋮
169 S371 8 350.781
170 S371 9 369.469
171 S372 0 269.412
172 S372 1 273.474
173 S372 2 297.597
174 S372 3 310.632
175 S372 4 287.173
176 S372 5 329.608
177 S372 6 334.482
178 S372 7 343.22
179 S372 8 369.142
180 S372 9 364.124 
m1 = let f = @formula reaction ~ 1 + days + (1 + days | subj)
  fit(MixedModel, f, sleepstudy; progress)
end
println(m1)

The response vector, y, has 180 elements. The fixed-effects coefficient vector, β, has 2 elements and the fixed-effects model matrix, X, is of size 180 × 2.

m1.y
180-element view(::Matrix{Float64}, :, 3) with eltype Float64:
 249.55999755859375
 258.7047119140625
 250.80059814453125
 321.4397888183594
 356.8518981933594
 414.6900939941406
 382.20379638671875
 290.1485900878906
 430.5852966308594
 466.3534851074219
   ⋮
 273.4739990234375
 297.5968017578125
 310.631591796875
 287.172607421875
 329.60760498046875
 334.4818115234375
 343.21990966796875
 369.1416931152344
 364.12359619140625
m1.β
2-element Vector{Float64}:
 251.4051060532072
  10.46728550560944
m1.X
180×2 Matrix{Float64}:
 1.0  0.0
 1.0  1.0
 1.0  2.0
 1.0  3.0
 1.0  4.0
 1.0  5.0
 1.0  6.0
 1.0  7.0
 1.0  8.0
 1.0  9.0
 ⋮    
 1.0  1.0
 1.0  2.0
 1.0  3.0
 1.0  4.0
 1.0  5.0
 1.0  6.0
 1.0  7.0
 1.0  8.0
 1.0  9.0

The second column of X is just the days vector and the first column is all 1’s.

There are 36 random effects, 2 for each of the 18 levels of subj. The “estimates” (technically, the conditional means or conditional modes) are returned as a vector of matrices, one matrix for each grouping factor. In this case there is only one grouping factor for the random effects so there is one one matrix which contains 18 intercept random effects and 18 slope random effects.

m1.b
1-element Vector{Matrix{Float64}}:
 [2.815658836479532 -40.04849255110683 … 0.723283769724057 12.118951000675557; 9.075536868353398 -8.644065257277084 … -0.9710555104626467 1.3106897770837675]
only(m1.b)   # only one grouping factor
2×18 Matrix{Float64}:
 2.81566  -40.0485   -38.4332   22.8323   …  -24.7104    0.723284  12.119
 9.07554   -8.64407   -5.51337  -4.65876       4.65974  -0.971056   1.31069

There is a model matrix, Z, for the random effects. In general it has one chunk of columns for the first grouping factor, a chunk of columns for the second grouping factor, etc.

In this case there is only one grouping factor.

Int.(first(m1.reterms))
180×36 Matrix{Int64}:
 1  0  0  0  0  0  0  0  0  0  0  0  0  …  0  0  0  0  0  0  0  0  0  0  0  0
 1  1  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 1  2  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 1  3  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 1  4  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 1  5  0  0  0  0  0  0  0  0  0  0  0  …  0  0  0  0  0  0  0  0  0  0  0  0
 1  6  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 1  7  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 1  8  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 1  9  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮        ⋱     ⋮              ⋮              ⋮
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  1  1
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  1  2
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  1  3
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  1  4
 0  0  0  0  0  0  0  0  0  0  0  0  0  …  0  0  0  0  0  0  0  0  0  0  1  5
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  1  6
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  1  7
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  1  8
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  1  9

The defining property of a linear model or linear mixed model is that the fitted values are linear combinations of the fixed-effects parameters and the random effects. We can write the fitted values as

m1.X * m1.β + only(m1.reterms) * vec(only(m1.b))
180-element Vector{Float64}:
 254.22076488968673
 273.7635872636495
 293.3064096376124
 312.8492320115752
 332.3920543855381
 351.93487675950087
 371.4776991334638
 391.02052150742657
 410.5633438813894
 430.10616625535226
   ⋮
 275.30203233657596
 287.0800076192691
 298.8579829019624
 310.6359581846556
 322.4139334673488
 334.191908750042
 345.9698840327352
 357.7478593154284
 369.5258345981216
fitted(m1)   # just to check that these are indeed the same as calculated above
180-element Vector{Float64}:
 254.22076488968673
 273.7635872636496
 293.3064096376124
 312.8492320115753
 332.3920543855381
 351.9348767595009
 371.4776991334638
 391.0205215074266
 410.5633438813894
 430.10616625535226
   ⋮
 275.30203233657596
 287.08000761926917
 298.8579829019624
 310.6359581846556
 322.4139334673488
 334.19190875004205
 345.9698840327352
 357.7478593154284
 369.5258345981216

In symbols we would write the linear predictor expression as \[ \boldsymbol{\eta} = \mathbf{X}\boldsymbol{\beta} +\mathbf{Z b} \] where \(\boldsymbol{\eta}\) has 180 elements, \(\boldsymbol{\beta}\) has 2 elements, \(\bf b\) has 36 elements, \(\bf X\) is of size 180 × 2 and \(\bf Z\) is of size 180 × 36.

For a linear model or linear mixed model the linear predictor is the mean response, \(\boldsymbol\mu\). That is, we can write the probability model in terms of a 180-dimensional random variable, \(\mathcal Y\), for the response and a 36-dimensional random variable, \(\mathcal B\), for the random effects as \[ \begin{aligned} (\mathcal{Y} | \mathcal{B}=\bf{b}) &\sim\mathcal{N}(\bf{ X\boldsymbol\beta + Z b},\sigma^2\bf{I})\\\\ \mathcal{B}&\sim\mathcal{N}(\bf{0},\boldsymbol{\Sigma}_{\boldsymbol\theta}) . \end{aligned} \] where \(\boldsymbol{\Sigma}_\boldsymbol{\theta}\) is a 36 × 36 symmetric covariance matrix that has a special form - it consists of 18 diagonal blocks, each of size 2 × 2 and all the same.

Recall that this symmetric matrix can be constructed from the parameters \(\boldsymbol\theta\), which generate the lower triangular matrix \(\boldsymbol\lambda\), and the estimate \(\widehat{\sigma^2}\).

m1.θ
3-element Vector{Float64}:
 0.9292259065867788
 0.018165149307614335
 0.2226454658393477
λ = only(m1.λ)  # with multiple grouping factors there will be multiple λ's
2×2 LinearAlgebra.LowerTriangular{Float64, Matrix{Float64}}:
 0.929226    ⋅ 
 0.0181651  0.222645
Σ = varest(m1) * (λ * λ')
2×2 Matrix{Float64}:
 565.516   11.0551
  11.0551  32.6822

Compare the diagonal elements to the Variance column of

VarCorr(m1)
Column Variance Std.Dev Corr.
subj (Intercept) 565.51589 23.78058
days 32.68220 5.71683 +0.08
Residual 654.94103 25.59182

2 Linear predictors in LMMs and GLMMs

Writing the model for \(\mathcal Y\) as \[ (\mathcal{Y} | \mathcal{B}=\bf{b})\sim\mathcal{N}(\bf{ X\boldsymbol\beta + Z b},\sigma^2\bf{I}) \] may seem like over-mathematization (or “overkill”, if you prefer) relative to expressions like \[ y_i = \beta_1 x_{i,1} + \beta_2 x_{i,2}+ b_1 z_{i,1} +\dots+b_{36} z_{i,36}+\epsilon_i \] but this more abstract form is necessary for generalizations.

The way that I read the first form is

The conditional distribution of the response vector, \(\mathcal Y\), given that the random effects vector, \(\mathcal B =\bf b\), is a multivariate normal (or Gaussian) distribution whose mean, \(\boldsymbol\mu\), is the linear predictor, \(\boldsymbol\eta=\bf{X\boldsymbol\beta+Zb}\), and whose covariance matrix is \(\sigma^2\bf I\). That is, conditional on \(\bf b\), the elements of \(\mathcal Y\) are independent normal random variables with constant variance, \(\sigma^2\), and means of the form \(\boldsymbol\mu = \boldsymbol\eta = \bf{X\boldsymbol\beta+Zb}\).

So the only things that differ in the distributions of the \(y_i\)’s are the means and they are determined by this linear predictor, \(\boldsymbol\eta = \bf{X\boldsymbol\beta+Zb}\).

3 Generalized Linear Mixed Models

Consider first a GLMM for a vector, \(\bf y\), of binary (i.e. yes/no) responses. The probability model for the conditional distribution \(\mathcal Y|\mathcal B=\bf b\) consists of independent Bernoulli distributions where the mean, \(\mu_i\), for the i’th response is again determined by the i’th element of a linear predictor, \(\boldsymbol\eta = \mathbf{X}\boldsymbol\beta+\mathbf{Z b}\).

However, in this case we will run into trouble if we try to make \(\boldsymbol\mu=\boldsymbol\eta\) because \(\mu_i\) is the probability of “success” for the i’th response and must be between 0 and 1. We can’t guarantee that the i’th component of \(\boldsymbol\eta\) will be between 0 and 1. To get around this problem we apply a transformation to take \(\eta_i\) to \(\mu_i\). For historical reasons this transformation is called the inverse link, written \(g^{-1}\), and the opposite transformation - from the probability scale to an unbounded scale - is called the link, g.

Each probability distribution in the exponential family (which is most of the important ones), has a canonical link which comes from the form of the distribution itself. The details aren’t as important as recognizing that the distribution itself determines a preferred link function.

For the Bernoulli distribution, the canonical link is the logit or log-odds function, \[ \eta = g(\mu) = \log\left(\frac{\mu}{1-\mu}\right), \] (it’s called log-odds because it is the logarithm of the odds ratio, \(p/(1-p)\)) and the canonical inverse link is the logistic \[ \mu=g^{-1}(\eta)=\frac{1}{1+\exp(-\eta)}. \] This is why fitting a binary response is sometimes called logistic regression.

For later use we define a Julia logistic function. See this presentation for more information than you could possibly want to know on how Julia converts code like this to run on the processor.

increment(x) = x + one(x)
logistic(η) = inv(increment(exp(-η)))
logistic (generic function with 1 method)

To reiterate, the probability model for a Generalized Linear Mixed Model (GLMM) is \[ \begin{aligned} (\mathcal{Y} | \mathcal{B}=\bf{b}) &\sim\mathcal{D}(\bf{g^{-1}(X\boldsymbol\beta + Z b)},\phi)\\\\ \mathcal{B}&\sim\mathcal{N}(\bf{0},\Sigma_{\boldsymbol\theta}) . \end{aligned} \] where \(\mathcal{D}\) is the distribution family (such as Bernoulli or Poisson), \(g^{-1}\) is the inverse link and \(\phi\) is a scale parameter for \(\mathcal{D}\) if it has one. The important cases of the Bernoulli and Poisson distributions don’t have a scale parameter - once you know the mean you know everything you need to know about the distribution. (For those following the presentation, this poem by John Keats is the one with the couplet “Beauty is truth, truth beauty - that is all ye know on earth and all ye need to know.”)

3.1 An example of a Bernoulli GLMM

The contra dataset in the MixedModels package is from a survey on the use of artificial contraception by women in Bangladesh.

contra = DataFrame(dataset(:contra))
1934×6 DataFrame
1909 rows omitted
Row dist urban urbdist livch age use
String String String String Float32 String
1 D01 Y U01 3+ 18.44 N
2 D01 Y U01 0 -5.56 N
3 D01 Y U01 2 1.44 N
4 D01 Y U01 3+ 8.44 N
5 D01 Y U01 0 -13.56 N
6 D01 Y U01 0 -11.56 N
7 D01 Y U01 3+ 18.44 N
8 D01 Y U01 3+ -3.56 N
9 D01 Y U01 1 -5.56 N
10 D01 Y U01 3+ 1.44 N
11 D01 Y U01 0 -11.56 Y
12 D01 Y U01 0 -2.56 N
13 D01 Y U01 1 -4.56 N
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮
1923 D61 N R61 0 -11.56 Y
1924 D61 N R61 3+ 1.44 N
1925 D61 N R61 1 -5.56 N
1926 D61 N R61 3+ 14.44 N
1927 D61 N R61 3+ 19.44 N
1928 D61 N R61 2 -9.56 Y
1929 D61 N R61 2 -2.56 N
1930 D61 N R61 3+ 14.44 N
1931 D61 N R61 2 -4.56 N
1932 D61 N R61 3+ 14.44 N
1933 D61 N R61 0 -13.56 N
1934 D61 N R61 3+ 10.44 N 
combine(groupby(contra, :dist), nrow)
60×2 DataFrame
35 rows omitted
Row dist nrow
String Int64
1 D01 117
2 D02 20
3 D03 2
4 D04 30
5 D05 39
6 D06 65
7 D07 18
8 D08 37
9 D09 23
10 D10 13
11 D11 21
12 D12 29
13 D13 24
⋮ ⋮ ⋮
49 D49 4
50 D50 19
51 D51 37
52 D52 61
53 D53 19
54 D55 6
55 D56 45
56 D57 27
57 D58 33
58 D59 10
59 D60 32
60 D61 42

The information recorded included woman’s age, the number of live children she has, whether she lives in an urban or rural setting, and the political district in which she lives.

The age was centered. Unfortunately, the version of the data to which I had access did not record what the centering value was.

A data plot, Figure 1, shows that the probability of contraception use is not linear in age - it is low for younger women, higher for women in the middle of the range (assumed to be women in late 20’s to early 30’s) and low again for older women (late 30’s to early 40’s in this survey).

If we fit a model with only the age term in the fixed effects, that term will not be significant. This doesn’t mean that there is no “age effect”, it only means that there is no significant linear effect for age.

Code 
draw(
  data(
    @transform(
      contra,
      :numuse = Int(:use == "Y"),
      :urb = ifelse(:urban == "Y", "Urban", "Rural"),
      :age = Float64(:age)
    )
  ) *
  mapping(
    :age => "Centered age (yr)",
    :numuse => "Frequency of contraception use";
    col=:urb,
    color=:livch,
  ) *
  smooth();
  figure=(; size=(800, 450)),
)
Figure 1: Smoothed relative frequency of contraception use versus centered age for women in the 1989 Bangladesh Fertility Survey 
contrasts = Dict(
  :urban => HelmertCoding(),
  :livch => DummyCoding(), # default, but no harm in being explicit
)
dist = Bernoulli()
gm1 = let
  form = @formula(
    use ~ 1 + age + abs2(age) + urban + livch + (1 | dist)
  )
  fit(MixedModel, form, contra, dist;  contrasts, progress)
end
Est. SE z p σ_dist
(Intercept) -0.6863 0.1686 -4.07 <1e-04 0.4786
age 0.0036 0.0092 0.39 0.6994
abs2(age) -0.0046 0.0007 -6.29 <1e-09
urban: Y 0.3483 0.0600 5.81 <1e-08
livch: 1 0.8146 0.1622 5.02 <1e-06
livch: 2 0.9157 0.1851 4.95 <1e-06
livch: 3+ 0.9143 0.1858 4.92 <1e-06

Notice that the linear term for age is not significant but the quadratic term for age is highly significant. We usually retain the lower order term, even if it is not significant, if the higher order term is significant.

Notice also that the parameter estimates for the treatment contrasts for livch are similar. Thus the distinction of 1, 2, or 3+ children is not as important as the contrast between having any children and not having any. Those women who already have children are more likely to use artificial contraception.

Furthermore, the women without children have a different probability vs age profile than the women with children. To allow for this we define a binary children factor and incorporate an age&children interaction.

VarCorr(gm1)
Column Variance Std.Dev
dist (Intercept) 0.2291 0.4786

Notice that there is no “residual” variance being estimated. This is because the Bernoulli distribution doesn’t have a scale parameter.

3.2 Convert livch to a binary factor

@transform!(contra, :children = ifelse(:livch ≠ "0", "Y", "N"))
# add the associated contrast specifier
contrasts[:children] = EffectsCoding()
EffectsCoding(nothing, nothing)
gm2 = let
  form = @formula(
    use ~
      1 +
      age * children +
      abs2(age) +
      children +
      urban +
      (1 | dist)
  )
  fit(MixedModel, form, contra, dist; contrasts, progress)
end
Est. SE z p σ_dist
(Intercept) -0.3614 0.1275 -2.84 0.0046 0.4756
age -0.0131 0.0110 -1.19 0.2350
children: Y 0.6055 0.1035 5.85 <1e-08
abs2(age) -0.0058 0.0008 -6.89 <1e-11
urban: Y 0.3568 0.0602 5.93 <1e-08
age & children: Y 0.0342 0.0127 2.69 0.0072
Code 
let
  mods = [gm2, gm1]
  DataFrame(;
    model=[:gm2, :gm1],
    npar=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    BIC=bic.(mods),
    AICc=aicc.(mods),
  )
end
2×6 DataFrame
Row model npar deviance AIC BIC AICc
Symbol Int64 Float64 Float64 Float64 Float64
1 gm2 7 2364.92 2379.18 2418.15 2379.24
2 gm1 8 2372.46 2388.73 2433.27 2388.81

Because these models are not nested, we cannot do a likelihood ratio test. Nevertheless we see that the deviance is much lower in the model with age & children even though the 3 levels of livch have been collapsed into a single level of children. There is a substantial decrease in the deviance even though there are fewer parameters in model gm2 than in gm1. This decrease is because the flexibility of the model - its ability to model the behavior of the response - is being put to better use in gm2 than in gm1.

At present the calculation of the geomdof as sum(influence(m)) is not correctly defined in our code for a GLMM so we need to do some more work before we can examine those values.

3.3 Using urban&dist as a grouping factor

It turns out that there can be more difference between urban and rural settings within the same political district than there is between districts.

dist_mean = combine(groupby(contra, :dist), 
                           :use => (x -> mean(x .== "Y")) => "dist_mean")
plt = data(dist_mean) * mapping(:dist_mean => "Distribution of district means") * AlgebraOfGraphics.density()
draw(plt)
dist_urban_mean = combine(groupby(contra, [:dist, :urban]), 
                           :use => (x -> mean(x .== "Y")) => "dist_urban_mean")
plt = data(dist_urban_mean) * mapping(:dist_urban_mean => "Distribution of district × urban means"; color=:urban) * AlgebraOfGraphics.density()
draw(plt)
dum_sorter = combine(groupby(dist_urban_mean, :dist),
                     :dist_urban_mean => diff => :urban_rural_diff)    
transform!(dum_sorter, :urban_rural_diff => ByRow(abs); renamecols=false)
all_dists = DataFrame(; dist=unique(dist_urban_mean.dist))
dum_sorter = leftjoin(all_dists, dum_sorter; on=:dist)
transform!(dum_sorter, 
          :urban_rural_diff => ByRow(x -> coalesce(x, 0));
          renamecols=false)
sort!(dum_sorter, :urban_rural_diff)

plt = data(dist_urban_mean) * 
    mapping(:dist_urban_mean => "Proportion contraception use",
            :dist => sorter(dum_sorter.dist) => "District") * 
    (mapping(; color=:urban) * visual(Scatter) + 
     mapping(; group=:dist) * visual(Lines))
draw(plt; figure=(;size=(500, 950), title="Distribution of district × urban means"), legend=(; position=:top))

To model this difference we build a model with urban&dist as a grouping factor.

gm3 = let
  form = @formula(
    use ~
      1 +
      age * children +
      abs2(age) +
      children +
      urban +
      (1 | urban & dist)
  )
  fit(MixedModel, form, contra, dist; contrasts, progress)
end
Est. SE z p σ_urban & dist
(Intercept) -0.3421 0.1269 -2.70 0.0070 0.5761
age -0.0129 0.0112 -1.16 0.2463
children: Y 0.6067 0.1045 5.80 <1e-08
abs2(age) -0.0056 0.0008 -6.66 <1e-10
urban: Y 0.3935 0.0859 4.58 <1e-05
age & children: Y 0.0332 0.0128 2.59 0.0096
Code 
let
  mods = [gm3, gm2, gm1]
  DataFrame(;
    model=[:gm3, :gm2, :gm1],
    npar=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    BIC=bic.(mods),
    AICc=aicc.(mods),
  )
end
3×6 DataFrame
Row model npar deviance AIC BIC AICc
Symbol Int64 Float64 Float64 Float64 Float64
1 gm3 7 2353.82 2368.48 2407.45 2368.54
2 gm2 7 2364.92 2379.18 2418.15 2379.24
3 gm1 8 2372.46 2388.73 2433.27 2388.81

Notice that the parameter count in gm3 is the same as that of gm2 - the thing that has changed is the number of levels of the grouping factor- resulting in a much lower deviance for gm3. This reinforces the idea that a simple count of the number of parameters to be estimated does not always reflect the complexity of the model.

gm2
Est. SE z p σ_dist
(Intercept) -0.3614 0.1275 -2.84 0.0046 0.4756
age -0.0131 0.0110 -1.19 0.2350
children: Y 0.6055 0.1035 5.85 <1e-08
abs2(age) -0.0058 0.0008 -6.89 <1e-11
urban: Y 0.3568 0.0602 5.93 <1e-08
age & children: Y 0.0342 0.0127 2.69 0.0072 
gm3
Est. SE z p σ_urban & dist
(Intercept) -0.3421 0.1269 -2.70 0.0070 0.5761
age -0.0129 0.0112 -1.16 0.2463
children: Y 0.6067 0.1045 5.80 <1e-08
abs2(age) -0.0056 0.0008 -6.66 <1e-10
urban: Y 0.3935 0.0859 4.58 <1e-05
age & children: Y 0.0332 0.0128 2.59 0.0096

The coefficient for age may be regarded as insignificant but we retain it for two reasons: we have a term of age² (written abs2(age)) in the model and we have a significant interaction age & children in the model.

3.4 Predictions for some subgroups

For a “typical” district (random effect near zero) the predictions on the linear predictor scale for a woman whose age is near the centering value (i.e. centered age of zero) are:

using Effects
design = Dict(
  :children => ["Y", "N"], :urban => ["Y", "N"], :age => [0.0]
)
preds = effects(design, gm3)
4×7 DataFrame
Row children age urban use: Y err lower upper
String Float64 String Float64 Float64 Float64 Float64
1 Y 0.0 Y 0.658034 0.150523 0.507511 0.808558
2 N 0.0 Y -0.555368 0.2305 -0.785868 -0.324868
3 Y 0.0 N -0.128908 0.113012 -0.24192 -0.0158953
4 N 0.0 N -1.34231 0.221603 -1.56391 -1.12071

We can plot this with a few more values for age:

design = Dict(
  :children => ["Y", "N"], 
  :urban => ["Y", "N"], 
  :age => -10:10
)
preds = effects(design, gm3; level=0.95)
base = data(preds) * mapping(:age; 
                            color=:children, 
                            col=:urban => renamer(["N" => "rural", "Y" => "urban"])) 
                             
                             
lines = mapping("use: Y") * visual(Lines)
bands = mapping(:lower, :upper) * visual(Band; alpha=0.3)


draw(base * (lines + bands),
     legend = (; position = :top, 
               framevisible=false),
     axis=(; ylabel="Log odds of contraception use",
           xlabel="Centered age"))

We can also plot this on the response scale, i.e. the probability scale:

preds = effects(design, gm3; invlink=AutoInvLink(), level=0.95)
base = data(preds) * mapping(:age; 
                            color=:children, 
                            col=:urban => renamer(["N" => "rural", "Y" => "urban"])) 
                             
                             
lines = mapping("use: Y") * visual(Lines)
bands = mapping(:lower, :upper) * visual(Band; alpha=0.3)


draw(base * (lines + bands);
     legend = (; position = :top, 
               framevisible=false),
     axis=(; ylabel="Probability of contraception use",
           xlabel="Centered age",
           limits=(nothing, (0, 1))))

4 Summarizing the results

  • From the data plot we can see a quadratic trend in the probability by age.
  • The patterns for women with children are similar and we do not need to distinguish between 1, 2, and 3+ children.
  • We do distinguish between those women who do not have children and those with children. This shows up in a significant age & children interaction term.

5 See Also

  • The binary-response GLMM chapter of Embrace Uncertainty — a full treatment of the same contraception data, including link functions and interpretation on the odds scale.

6 Exercises

  1. Interpret on the right scale. For the logistic GLMM fit to the contra data, a coefficient is reported on the logit scale. How do you turn it into a statement about probability, and why can’t you read it directly as a probability change?

Coefficients are additive on the log-odds (logit) scale; exponentiating gives an odds ratio. Because the logistic link is nonlinear, the same log-odds change corresponds to a different probability change depending on the baseline probability, so there is no single “probability per unit” — you evaluate predicted probabilities at specific covariate values (e.g. with an effects/marginal-means calculation).

  1. Link and distribution. What two choices distinguish a GLMM from an LMM, and which part of the fit(MixedModel, ...) call encodes them?

A GLMM adds (i) a conditional distribution for the response (e.g. Bernoulli(), Poisson()) and (ii) a link function relating the linear predictor to the conditional mean. Both are passed to fit/GeneralizedLinearMixedModel — the distribution as a positional argument and the link via the link keyword (each distribution has a canonical default link).

This page was rendered from git revision 29e8d33 using Quarto 1.10.18 and Julia 1.12.7.

Back to top