Fitting Generalized Linear Models

Darren Glosemeyer

Wolfram Research

What is a Generalized Linear Model?

First, let's start with the Linear Model.

The Linear Model

Given a set of predictor variables (or functions) (Subscript[x, 1]), …  (Subscript[x, n]) and a vector of observed responses (Overscript[y,  ⇀])

(Overscript[y,  ⇀]) =X.(Overscript[β,  ⇀]) +(Overscript[ϵ,  ⇀])

(Subscript[ϵ, i])~N(0,(σ^2))

Where X is a design matrix, (Overscript[β, ⇀]) is a parameter vector to be determined, and the elements of (Overscript[ϵ, ⇀]) are independent identically distributed normal errors.

Equivalently, for each fitted value \(TraditionalForm\`\(y\_i\)\&^\), we can write

\(TraditionalForm\`\(y\_i\)\&^ = ∑\+\(j = 1\)\%n\(\( β\&^\)\_j\) x\_\(i,j\)\)

with measured y values coming a normal distribution centered at the true value (Subscript[μ, i]) given the measured values of (Subscript[x, 1]), …  (Subscript[x, n]).

(Subscript[Y, i])~N((Subscript[μ, i]),(σ^2))

The linear model is fitted by minimizing the sum or squared errors, which is equivalent to maximizing the likelihood function under the assumed normal distribution.

The Generalization

The Generalized Linear Model (GLM) takes the following form.

\(TraditionalForm\`\(y\_i\)\&^ = \(g\^\(-1\)\)(∑\+\(j = 1\)\%n\(\( β\&^\)\_j\) x\_\(i,j\))\)

(Subscript[Y, i])~𝒟((Subscript[μ, i]),ϕ var((Subscript[μ, i])))

Where g is a smooth function called the link function, and the i th random variable (Subscript[Y, i]) is assumed to be from a distribution 𝒟 with mean (Subscript[μ, i]) and variance ϕ var((Subscript[μ, i])), which is a function of the mean (Subscript[μ, i]) and a dispersion constant ϕ, and the (Subscript[Y, i])s are assumed to be independent.

Some Common Cases

Some of the more common GLMs include linear regression, logistic regression, probit regression, loglinear models for count data, and exponential dispersion models.

({{"Distribution", "Canonical Link", "Alternate Links"}, {"Normal(μ,\!\(σ\^2\))", "μ", "log(μ), \!\(μ\^λ\), \!\(\(μ\^λ - 1\)\/λ\) assuming x>0"}, {"\!\(\(Binomial \((n, p)\)\)\/n\)", "log(\!\(p\/\(1 - p\)\))", "\!\(Φ\^\(-1\)\)(p), log(-log(1-p)), -log(-log(p))"}, {"Poisson(μ)", "log(μ)", "\!\(μ\^λ\), \!\(\(μ\^λ - 1\)\/λ\)"}, {"Gamma(α,β)", "\!\(1\/μ\)", "log(μ), \!\(μ\^λ\), \!\(\(μ\^λ - 1\)\/λ\)"}, {"InverseGaussian(μ,λ)", "\!\(1\/μ\^2\)", "log(μ), (μ^λ), (((μ^λ) - 1)/λ)"}})

Fitting a GLM

To fit a generalized linear model, we want to maximize the likelihood function

(∏_(i = 1)^m) f(
        (Subscript[y, i]);  (g^(-1))((∑_(j = 1)^n) (Subscript[β, j]) 
              (Subscript[x, i,j])))

or equivalently, maximize the loglikelihood function

(∑_(i = 1)^m) 
        log (f(
            (Subscript[y, i]); (g^(-1))((∑_(j = 1)^n) (Subscript[β, j]) 
                  (Subscript[x, i,j]))) )

with respect to the parameter vector  (Overscript[β, ⇀])={(Subscript[β, 1]), … (Subscript[β, n])} for the assumed response distribution 𝒟.

If the distribution 𝒟 is not the normal distribution, minimizing the sum of squared errors no longer gives the optimum value for (Overscript[β,  ⇀]) .

Instead, an optimum (Overscript[β, ⇀]) can be obtained via direct maximum likelihood or iteratively reweighted least squares.

An Example

First consider the loglinear model for count data.  In this model, we have

g ((Subscript[μ, i])) = log((Subscript[μ, i]))

(Subscript[Y, i])~Poisson((Subscript[μ, i]))

var((Subscript[μ, i])) = (Subscript[μ, i])

For a single response modeled a linear combination of basis functions, the loglikelihood can be constructed as follows.

<< Statistics`

PoissonLogLike[beta_, basis_, resp_] := 
      PowerExpand[
        Log[PDF[PoissonDistribution[Exp[beta . basis]], resp]]]

PoissonLogLike[{(Subscript[β, 1]), (Subscript[β, 2]), (Subscript[β, 3])}, {(Subscript[x, 1]), (Subscript[x, 2]), 
        (Subscript[x, 3])}, y]

Likelihood Maximization via Newton's Method

Generalized linear models are well suited for maximum likelihood estimation via Newton's method.

poissondata = 
        Block[{$RandomState}, SeedRandom[1]; 
          Table[{x, 
              Random[PoissonDistribution[
                  Exp[{2.3, 1.4, - .6} . {1, x, x^2}]]]}, {x,  .1, 
              3,  .1}]];

ListPlot[poissondata, PlotStyle -> {PointSize[ .02]}]

The loglikelihood for the data set can be obtained as a sum of the loglikelihoods for the individual data points.

loglikefun[b1_, b2_, b3_] = 
        Total[MapThread[PoissonLogLike[{b1, b2, b3}, {1, #1, #1^2}, #2] &, 
            Transpose[poissondata]]];

Good initial estimates for the parameter values can often be obtained from a least squares fitting.

linres = 
      LinearAlgebra`LeastSquares[DesignMatrix[poissondata, {x, x^2}, x], 
        Log[poissondata[[All, -1]]]]

The Poisson loglikelihood can then be maximized using FindMaximum.

poissonmaxlike = 
      FindMaximum[loglikefun[b1, b2, b3], 
        Evaluate[Transpose[{{b1, b2, b3}, linres}]], 
        Method -> Newton]

The maximum likelihood estimate gives a slightly higher loglikelihood value than the least squares fitting.

Apply[loglikefun, linres]

Iteratively Reweighted Least Squares (IRLS)

An alternative to directly maximizing the loglikelihood function is iteratively reweighted least squares (see, for instance, McCullagh and Nelder).

The basic idea is to iteratively fit the linear model

E((Subscript[((Overscript[z,  ⇀]) ), (j)]))=X.(Subscript[((Overscript[β,  ⇀]) ), (j)])  with regression weights(Subscript[W, (j)])

until (Subscript[(Overscript[β, ⇀]), (j)]) (or the loglikelihood for the original model) converges to within a specified tolerance, where the values of (Subscript[(Overscript[z, ⇀]), (j)]) and (Subscript[(Overscript[W, ⇀]), (j)]) at the j th iteration are

(Subscript[((Overscript[z,  ⇀]) ), (j)])=X.(Subscript[((Overscript[β,  ⇀]) ), (j)])+((Overscript[y,  ⇀]) -(g^(-1))(X.(Subscript[((Overscript[β,  ⇀]) ), (j)])))(Subscript[\((\( \)\(∂\(g(μ)\)\/∂μ\))\), μ=(g^(-1))(X.(Subscript[((Overscript[β,  ⇀]) ), (j)]))])

(Subscript[((Overscript[W,  ⇀]) ), (j)])=((((Subscript[(∂(μ/∂)g(μ)), μ=(g^(-1))(X.(Subscript[((Overscript[β,  ⇀]) ), (j)]))])))^2)(1/(var((g^(-1))(X.(Subscript[((Overscript[β,  ⇀]) ), (j)])))))

This method is similar to Newton's method with expected values of the response instead of observed values.  Initial values can be computed from the observed response or a minor adjustment of the observed response.

The following function fits a GLM using the iteratively reweighted least squares.  In addition to the data, basis functions, and variables, the link function for the GLM and the variance function var(μ) for the assumed response distribution are also needed.

glmFit[data_, basis_, vars_, link_, varfun_, wts_:  1, 
        max_:  100]:=
Block[{xmat, response, start, step, beta, invlink, partder, 
            wtfun},
invlink = 
              Evaluate[
                  If[MatchQ[
                      link[#], # | #^_], (#^(1/
                            Exponent[link[#], #])), 
                    Simplify[
                      mu /. Internal`DeactivateMessages[Solve[
                              link[mu] == #, {mu}][[1]]]]]] &;
xmat = 
              DesignMatrix[data, basis, vars, 
                IncludeConstant -> False];
response = data[[All, -1]];
start = 
              LinearAlgebra`LeastSquares[Sqrt[wts]*xmat, 
                Sqrt[wts]*
                  Map[link[
                        Which[# == 1,  .99, # == 0,  .01, 
                          True, #]] &, response]];
partder = Evaluate[D[link[#], #]] &;
wtfun = 
              Evaluate[
                  Simplify[(D[invlink[x], x] /. x -> link[#])^2/
                      varfun[#]]] &;
step[beta_]:=Block[{linearResponse, invlinkLinResponse, z, 
                    weightvals},
linearResponse = xmat . beta;
invlinkLinResponse=Map[invlink,linearResponse];
z = (linearResponse + (response - invlinkLinResponse)*Map[partder, invlinkLinResponse]);
weightvals=wts*Map[wtfun,invlinkLinResponse];
LinearAlgebra`LeastSquares[
                      Sqrt[weightvals]*xmat, Sqrt[weightvals]*z]
];
NestWhile[step, 
              start, (Norm[#1 - #2] > 
                    (10^(-8)) + Norm[#2]*(10^(-8))) &, 2, max]
]

In the loglinear model, the link is Log, and the variance function is var(μ)=μ.

glmfitres=glmFit[poissondata, {1, x, x^2}, {x}, Log, Identity]

{b1, b2, b3} /. poissonmaxlike[[2]]

{Apply[loglikefun, glmfitres], poissonmaxlike[[1]]} // 
      InputForm

The fitted parameters and the loglikelihood agree with the direct maximum likelihood results to well within the tolerances.

The iteratively reweighted least squares approach will tend to be faster in general because it does not deal directly with the gradient of the loglikelihood function.

Do[FindMaximum[loglikefun[b1, b2, b3], 
          Evaluate[Transpose[{{b1, b2, b3}, linres}]], 
          Method -> Newton], {10}] // Timing

Do[glmFit[poissondata, {1, x, x^2}, {x}, Log, 
            Identity],{10}]//Timing

A Binomial Example

Logistic and probit regression are common GLMs for probabilities based on the binomial distribution.

(Subscript[Y, i])~b((Subscript[n, i]), (Subscript[p, i]))

var((Subscript[p, i])) = (Subscript[p, i]) (1 - (Subscript[p, i]))

For the logistic regression model,

g ((Subscript[p, i])) = 
      log(((Subscript[p, i])/(1 - (Subscript[p, i]))))

For the probit model

g ((Subscript[p, i])) = (Φ^(-1))(
        (Subscript[p, i]))

where (Φ^(-1))is the inverse cdf for the standard normal distribution.

Here we generate some random data that could be modeled by a binomial GLM.

binomdata = 
        Block[{$RandomState}, SeedRandom[1234]; 
          Table[With[{n = 
                  Random[Integer, {1, 20}]}, {x, {Random[
                    BinomialDistribution[
                      n, (Exp[#]/(1 + Exp[#])) &[{-1, 
                            2} . {1, x}]]], n}}], {x,  .1, 
              3,  .1}]];
Take[binomdata, 5]

The responses are pairs of number of successes and total number of trials. It will also be useful to express the data in terms of success probabilities.

probdata = binomdata;
probdata[[All, -1]] = 
        Apply[Divide, Transpose[probdata[[All, -1]]]];

ListPlot[probdata, PlotStyle -> {PointSize[ .02]}]

First, define the link functions of interest.

logit[x_] = Log[x/(1 - x)]

probit[x_] = Quantile[NormalDistribution[0, 1], x]

The loglikelihood for a single binomial response is

BinomialLogLike[beta_, basis_, invlink_, resp_] := 
      PowerExpand[
        Log[PDF[BinomialDistribution[resp[[2]], invlink[beta . basis]], 
            resp[[1]]]]]

where resp is a list of the number of successes and the total number of trials.

The loglikelihood function assuming a model of the form invlink[b1+b2*x] is

binloglikefun[b1_, b2_, invlink_] = 
        Total[MapThread[BinomialLogLike[{b1, b2}, {1, #1}, invlink, #2] &, 
            Transpose[binomdata]]];

Likelihood Maximization via Newton's Method

The direct maximum likelihood approach is basically the same as in the Poisson case, however, the probabilities are weighted by the associated number of trials.  Observed probabilities of 1 and 0 much also be adjusted slightly when obtaining a starting point.

wtvals = binomdata[[All, -1, -1]];

logitstart = 
      LinearAlgebra`LeastSquares[
        Sqrt[wtvals]*DesignMatrix[probdata, {1, x}, x], 
        Sqrt[wtvals]*
          Map[Which[# == 0, logit[ .01], # == 1, logit[ .99], 
                True, logit[#]] &, probdata[[All, -1]]]]

This gives the logistic regression fitting.

logitmaxlike = 
      FindMaximum[binloglikefun[b1, b2, (Exp[#]/(1 + Exp[#])) &], 
        Transpose[{{b1, b2}, logitstart}], Method -> Newton]

probitstart = 
      LinearAlgebra`LeastSquares[
        Sqrt[wtvals]*DesignMatrix[probdata, {1, x}, x], 
        Sqrt[wtvals]*
          Map[Which[# == 0, Log[ .01/ .99], # == 1, 
                Log[ .99/ .01], True, probit[#]] &, 
            probdata[[All, -1]]]]

And here is the probit fitting.

probitmaxlike = 
      FindMaximum[binloglikefun[b1, b2, (1 + Erf[#/Sqrt[2]])/2 &], 
        Transpose[{{b1, b2}, probitstart}], Method -> Newton]

glmFit can be used to same the same models with probdata as the data to be fit and the numbers of trials for each observation as weights.

glmlogit=glmFit[probdata, {1, x}, {x}, 
        logit, (# (1 - #)) &, wtvals]

glmprobit=glmFit[probdata, {1, x}, {x}, 
        probit, (# (1 - #)) &, wtvals]

Again, the parameters vectors are very close

Norm[glmlogit - {b1, b2} /. logitmaxlike[[2]]]

Norm[glmprobit - {b1, b2} /. probitmaxlike[[2]]]

and the loglikelihoods are nearly identical.

{logitmaxlike[[1]], 
        binloglikefun[b1, b2, (Exp[#]/(1 + Exp[#])) &] /. 
          Thread[Rule[{b1, b2}, glmlogit]]} // InputForm

{probitmaxlike[[1]], 
        binloglikefun[b1, b2, (1 + Erf[#/Sqrt[2]])/2 &] /. 
          Thread[Rule[{b1, b2}, glmprobit]]} // InputForm

Conclusion

We have seen examples of fitting three common generalized linear models via two methods. The same strategies can be used to fit other types of generalized linear models.

In the examples used, convergence of the parameter vectors in the IRLS method also gave convergence of the loglikelihood to within similar tolerances.  While this is often the case, it is not always the case. Additional checks on the loglikelihood can be added in the iterative fitting process to require convergence in the loglikelihood. Likewise, checking for and removal of small imaginary parts may be needed within the IRLS process clean numerical errors that may arise from some weighting functions.

References

McCullagh, P. and Nelder, J. A. (1989), Generalized Linear Models, Second Edition, Chapman & Hall/CRC, London.