Estimation Techniques for Generalized Linear Mixed Models with Binary Outcomes: Application in Medicine

2.1 Linear models

Linear models are the most common statistical models in regression analysis. It is commonly applied for its simplicity in many different statistical regression frameworks. A one-dimensional type of the model is defined:

$y=\beta_o+\beta_{i 1} x+e_i \quad i=1,2....$          (1)

where, x and y represent the predictor and response variable respectively. β_o and β₁ are model parameters. The e_i values are error terms, which follows normal distribution.

Models involving multiple predictors can be represented using a matrix format where p indicates number of parameters β₁, ..., β_p. The model can be represented as follows:

$y_i=\beta_1+\beta_2 x_{i 2}+\beta_3 x_{i 3}+\ldots \ldots \ldots \ldots+\beta_p x_{i p}+e_i$          (2)

We can rewrite Eq. (2) in vector as follows:

$y=X \beta+e$          (3)

Assuming an intercept in X such that x_i1=1, the matrix format of the model for y and x can be written in matrix form as:

$\begin{aligned} Y=\left[\begin{array}{c}y_1 \\ \vdots \\ y_n\end{array}\right], X & =\left[\begin{array}{rrrrr}1 & x_{12} & x_{13} & & x_p \\ 1 & x_{22} & x_{23} & \cdots & x_{1 p} \\ & \vdots & & \ddots & \vdots \\ 1 & x_{n 2} & x_{n 3} & \cdots & x_{n p}\end{array}\right] \\ \beta & =\left[\begin{array}{c}\beta_1 \\ \vdots \\ \beta_p\end{array}\right], e=\left[\begin{array}{c}e_1 \\ \vdots \\ e_n\end{array}\right]\end{aligned}$          (4)

Following Eq. (1), e_i~N(0, σ), E(e_i)=0; Var(e_i)=σ². In this case, a maximum likelihood method is engaged for parameter estimation. For classical linear models (LMs), the estimated parameters using maximum likelihood estimation technique is equivalent to minimized least squares method.

2.2 Maximum likelihood estimation techniques

The point in the parameter space that does the maximization of the likelihood function is referred to as maximum likelihood estimate. The maximum likelihood is flexible, and it has become a dominant means in statistical inference.

Since y~N(xβ,σ²l), the function of the likelihood is given by:

$L\left(\beta, \sigma^2\right)=\left(\frac{1}{2 \pi^2}\right)^{\frac{n}{2}} \exp \left(\frac{1}{2 \sigma^2}(y-x \beta)^{\prime}(y-x \beta)\right)$          (5)

Using logarithms, we have:

$\begin{aligned} L\left(\beta, \sigma^2\right)= & -\frac{n}{2} \log 2 \pi-\frac{n}{2} \log \sigma^2-\frac{1}{2 \sigma^2}(y- x \beta)^{\prime}(y-x \beta)= \\& -\frac{n}{2} \log 2 \pi-\frac{n}{2} \log \sigma^2-\frac{1}{2 \sigma^2}(' y-\left.2 y^{\prime} x \beta+\beta^{\prime} x^{\prime} x \beta\right)\end{aligned}$          (6)

To maximize the likelihood function, we set Eq. (6) to zero and take partial derivatives with respect to β and σ²:

$\frac{\vartheta}{\vartheta \beta} \log L\left(\beta, \sigma^2\right)=-\frac{n}{2 \sigma^2}+\frac{1}{2 \sigma^4}(y-x \beta)^{\prime}(y-x \beta)$          (7)

Equating to zero, we have:

$\hat{\beta}\left(x^{\prime} x\right)^{-\prime} x^{\prime} y$ and $\hat{\sigma}^2=\frac{1}{n}(y-x \beta)^{\prime}(y-x \beta)$          (8)

Assuming a variable Y relates to k predictor variables x₁, x₂, x₃, ..., x_k and a disturbance term e. If we have a sample of $n$-observation on y and the X's, as we have in Eq. (2) we can write:

$y_i=\beta_a+\beta_1 X_{1 i}+\cdots+\beta_k X_{k i}+e_i \quad i=1,2, \ldots, n$          (9)

The beta coefficients and the parameter of the $e$ distribution is unknown. To estimate the unknown parameters, where β_o=intercept and β₁, …, β_k=partial slope coefficients.

Unknown parameters can be obtain using matrix:

$\beta=\left(X^{\prime} X\right)^{-\prime} X^{\prime} Y$          (10)

The sum of squared residuals of n-residuals of (y-xβ) can be expressed as:

$\begin{aligned} & \sum_{i=1}^n e_i^2=e^{\prime} e=(y-x \beta)^{\prime}(y-x \beta)\quad=Y^{\prime} Y-2 \beta^{\prime} X^{\prime} X^{\prime} Y+\beta^{\prime} X^{\prime} X \beta\end{aligned}$          (11)

Noting that β'X'Y is a scalar that equals to its transpose Y'Xβ. Minimizing the sum of squared residuals, and differentiating Eq. (11), we have:

$\begin{gathered}\frac{\delta\left(e^{\prime} e\right)}{\delta \beta}=-2 X^{\prime} Y+2 X^{\prime} X \beta \\ \beta=\left(X^{\prime} X\right)^{-\prime} X^{\prime} Y \\ \beta=\left(\begin{array}{lllll}\beta_0 & \beta_1 & \ldots & \beta_k\end{array}\right)^{\prime}\end{gathered}$          (12)

where,

$\begin{aligned}

& X^{\prime} X=\left[\begin{array}{cccc}

1 & 1 & \cdots & 1 \\

X_{11} & X_{12} & \cdots & X_{1 n} \\

\vdots & \vdots & \vdots & \vdots \\

X_{k 1} & X_{k 2} & \cdots & X_{k n}

\end{array}\right]\left[\begin{array}{cccc}

1 & X_{11} & \cdots & X_{k 1} \\

1 & X_{12} & \cdots & X_{k 2} \\

\vdots & \vdots & \vdots & \vdots \\

1 & X_{1 n} & \cdots & X_{k n}

\end{array}\right] \\

& {\left[\begin{array}{cccc}

1 & 1 & \ldots & 1 \\

X_{11} & X_{12} & \ldots & X_{1 n} \\

\vdots & \vdots & \vdots & \vdots \\

X_{k 1} & X_{k 2} & \ldots & X_{k n}

\end{array}\right]\left[\begin{array}{c}

Y_1 \\

Y_2 \\

\vdots \\

Y_n

\end{array}\right]=\left[\begin{array}{c}

\sum Y_i \\

\sum X_{1 i} Y_{1 i} \\

\vdots \\

\sum X_{k i} Y_{k i}

\end{array}\right]}

\end{aligned}$          (13)

2.3 Generalized linear models

Using central limit theorem assumptions, many real events are not normally distributed. Linear models (LMs) commonly used sufficient models for tackling such assumptions. Generalized linear models (GLMs) tackles the non-normality challenges by providing solutions to such statistical analyses. Nelder and Wedderburn [20] established the theory of generalized linear models (GLMs). The generalized linear models are extensions of classical linear models that allows the mean of a population to depend on a linear predictor through a nonlinear link function. GLM uses probability distributions that belongs to the exponential class of family. GLM consists of the following components, it is of the form:

$f\left(y_i ; \theta_i\right)=\exp \left(\frac{y_i \theta_i-b\left(\theta_i\right)}{a_i(\varphi)}+c\left(y_i \varphi\right)\right)$          (14)

For a random sample Y₁, ..., Y_n, the linear component is defined as ɳ_i=X_iβ where i=1, 2, ..., n for some vector parameters β=(β₁, ..., β_p and covariate X=(x_i1, ..., x_ip) associated with observations Y_i. For a careful choice of α_i≠0, b and c, φ represent the scale parameters  representing the natural parameter. A monotonic differentiable link function g describes how the expected response μ_i=E(Y_i) is related to the linear predictor g(μ_i)=η_i. Link functions express a connection between a function E(Y) to a linear predictor ɳ. Classical linear regression models express the identity link g(μ_i)=ɳ_i. Gamma or Poisson distributions use links that are strictly positive like in claim counts and severity. The summary of link functions for exponential family members of distributions are expressed in Table 1.

Table 1. Link functions of some exponential distributions

2.4 Parameter estimation of GLMMs

The condition where various methods of estimating the parameters of GLMMs remains unclear in literature which this study aim to shed more light on. Estimating parameters of models is fundamental in statistical analyses. GLMMs parameters are described as fixed and random-effect parameters. Normal response variables help in fixing models using maximum likelihood (ML).

Treatments with equal sample sizes (i.e. balanced design with all nested random effects like classical ANOVA methods based on computing differences of sums of squares with ML approaches). However, this equivalence breaks down for more complex LMMs or GLMMs. To fix the problem of finding the ML estimates, one must integrate likelihoods over all possible values of the random effects. To describe this, a computation of the likelihood can be expressed as:

$L=\int f_y / u(y / u) f u(u) d u$          (15)

Or individually

$L=\int \prod_{i-1}^n \frac{f y_i}{u\left(\frac{y_i}{u}\right)} f u(u) d u$          (16)

Evaluating exponential family functions, we have a likelihood equation of the form:

$L=\int \prod_i^n \exp \left[\frac{y_i \theta_i-b\left(\theta_i\right)}{\phi}+c\left(y_i ; \phi\right)\right] \frac{1}{\sigma_u \sqrt{2 \pi}} \exp \left[\frac{-(x-\mu)^2}{2 \sigma_u^2}\right] d u$        (17)

With corresponding log-likelihood:

$L=\log \left[\int \prod_i^n \exp \left[\frac{y_i \theta_i-b\left(\theta_i\right)}{\phi}+c\left(y_i ; \phi\right)\right] \frac{1}{\sigma_u \sqrt{2 \pi}} \exp \left[\frac{-(x-\mu)^2}{2 \sigma_u^2}\right] d u\right]$         (18)

From Eq. (18), resolving to analytical solution is not in view, therefore numerical approximation procedures can be adopted in estimating the likelihood. Different methods of approximation are now considered. We now examine three different methods.

2.5 Penalized quasi-likelihood estimation procedures

If the exact distribution of data is unknown, then estimation of variance component is sorted out using quasi-likelihood methods by introducing a penalty term on random effects. The purpose of including penalty is to avoid some arbitrary values of random effects to force the random effect approximate to zero. Approximation for vector response data y_i is given by:

$y \approx \mu_i+\varepsilon_i=E\left(y_i / b\right)+\epsilon_i$          (19)

where, $e=h\left(x_i^{\prime} \beta+Z_i^{\prime} b\right)+\epsilon_i$ and $\beta=$ fixed vector parameter. Using Taylor expansion in Eq. (19), we have:

$\begin{gathered}y_i \approx h\left(x_i^{\prime}+Z_i^{\prime} b\right)+h^{\prime}\left(x_i^{\prime} \hat{\beta}+Z_i^{\prime} \hat{b}\right) x_i^{\prime}(\beta-\hat{\beta})+ h^{\prime}\left(x_i^{\prime} \hat{\beta}+Z_i^{\prime} \hat{b}\right) Z_i^{\prime}(\beta-\hat{\beta})+\varepsilon_i\end{gathered}$          (20)

$=\mu_i+V\left(\mu_i\right) x_i^{\prime}(\beta-\hat{\beta})+V(\hat{\mu}) Z_i^{\prime}(\beta-\hat{\beta})+\varepsilon_i$          (21)

Eq. (21) is expressed as $y_i^*=V_i^{-\prime}\left(y_i-\mu_i\right)+x_i \hat{\beta}+Z_i \hat{b}_i$, where, $y_i^*$=pseudo - response.

The mode fitting is done iteratively using algorithms that give conditional variance for as:

${var}\left(y_i / b\right)=a_i^{\prime}(\emptyset) V\left(\mu_i\right)$          (22)

where, ∅ reprents the dispersed Parameter.

2.6 Laplace with approximation methods

This engages Taylor expression in an exponential format for approximates integrals:

$\int e^{h(u)} d u$          (23)

where, u=q-dimensinal vector and h(u) sufficiency smooth function.

Taylor expansion second order for h in u_o can be described by:

$h(u) \approx h\left(u_o\right)+\frac{1}{2}\left(u-u_o\right)^{\prime} h^{\prime \prime}\left(u_o\right)\left(u-u_o\right)$          (24)

From Eq. (24) using a Laplace approximate function:

$\int e^{h(u)} d u \approx \exp \left[h\left(u_o\right)\right](2 \pi)^{\frac{q}{2}}\left|-h^{\prime \prime}\left(u_o\right)\right|^{\frac{-1}{2}}$          (25)

Approximating the likelihood function with Laplace, we have:

$\begin{gathered}L=\log \int f y / u(y / u) f u(u) d u=\log \int \exp [\log f y / u(y / u)]+\log f u(u) d u\end{gathered}$          (26)

$\frac{\partial^2 \log f u}{\partial u \partial u|}=-D^{-1}$          (27)

Exponential family with chain rules, we have the following:

$\frac{\frac{\partial \operatorname{logfy}}{u\left(\frac{y}{u}\right)}}{\partial u}=\frac{1}{\varphi} \sum_i\left(y_i \frac{\partial \theta_i}{\partial u}-\frac{\partial b\left(\theta_i\right) \partial \theta_i}{\partial \theta_i \partial u}\right)=\frac{1}{\varphi} \sum_i\left(y_i-\mu_i\right) \frac{1}{v\left(\mu_i\right)} \frac{1}{g^{\prime \prime}\left(u_i\right) Z_i^{\prime}}$          (28)

$\begin{aligned} & \frac{\partial^2 h(u)}{\partial u \partial u^{\prime}}=\frac{\partial}{\partial u^{\prime}}\left(\frac{1}{\varphi} Z^{\prime} W \Delta(y-\mu)-D^{-1} u\right)= \frac{1}{\varphi}\left(-Z^{\prime} W \Delta \frac{\partial \mu}{\partial u^{\prime}}+Z^{\prime} \frac{\partial W}{\partial u^{\prime}}(y-\mu)-D^{-1}\right)\end{aligned}$          (29)

$\frac{\partial^2 h(u)}{\partial u \partial u^{\prime}}=-\frac{1}{\varphi}\left(Z^{\prime} W Z D+I\right) D^{-1}$          (30)

$\begin{aligned} L \approx & \log f y / u\left(y / u_o\right)-\frac{1}{2} u_o^{\prime} D^{-1} u_o- \frac{1}{2} \log \left|\left(\frac{1}{\varphi} Z^{\prime} W Z D+I\right) D^{-1}\right|\end{aligned}$          (31)

$\begin{gathered}\frac{\partial l}{\partial \beta}=\frac{\partial \log f y / u\left(y / u_o\right)}{\partial \beta}+\frac{\partial}{\partial \beta} \frac{1}{2} \log \left|\frac{Z^{\prime} W Z D}{\varphi}+I\right| \approx \frac{1}{\varphi} X^{\prime} W \Delta(y-\mu)\end{gathered}$          (32)

Changing in β provides the estimate of β and u and by resolving the Eq. (32), we have:

$\frac{1}{\varphi} X^{\prime} W \Delta(y-\mu)=0$ and $\frac{1}{\varphi} Z^{\prime} W \Delta(y-\mu)=D^{-1} u$

2.7 Gauss-Hermite quadrature steps

The steps involving Gauss-Hermite quadrature (AGQ) apply a Gaussian technique, where gaussian functions replace the factor exp(-Z²) with suitable shift in the weights and approximation points by following the outline and subjecting it to GLMMs. Let ϕ(t:μ,σ) represent a probability density function of the normal distribution with μ and standard deviation .Let’s define a function g(t) in such a way that g(t)>0 is unimodal (i.e. with unique mode) and is sufficiently smooth. The aim is to approximate  by transformation. To establish this, we replace the gauss-Hermite quadrature for the integral:

$\int f(t) \phi(t, \mu, \sigma) d t$          (33)

It approximates the likelihood by taking optimal subdivisions at which to solve the integrand. Updates information from the initial fit to increase precision. This involves transforming sampling nodes x_i to t_i from exp[Z_i] to ϕ(t:μ,σ) and $t=\mu+\sqrt{2 \sigma Z_i}$.

However, we can sample the integral in region of g(t) with μ as the mode of g(t) and $\sigma=\frac{1}{\sqrt{j}}, j=-\frac{\partial^2}{\partial t^2} \log (g(t))$ and defining $h(t)$ as $\frac{g(t)}{\phi(t ; \mu, \sigma)}$.

Rewriting the integral for g(t), we have:

$\left.\int g(t) d t=\int h(t) \phi(t, \mu, \sigma)\right) d t$          (34)

Using the transformed Gauss-Hermite quadrature, we have:

$\int g(t) d t=\sqrt{2 \sigma} \sum_{i=1}^Q w_i^* g\left(\mu+\sqrt{2 \sigma Z_i}\right)$          (35)

Using a GLMM form, we present a procedure of a single random effect. This effect is seen as being clustered into different clusters. All clusters with random effects which are distributed as u_i~N(0, σ²). Hence we determine the posterior mode of u_i which depends on the parameters β, ϕ and σ and u_i. We replace this by the current estimate (and in the first step a well-chosen value) β^*, ϕ^* and σ^*. Applying these estimates on u_i, we maximize Using u_i as the mode for ui together with the gauss-Hermite quadrature to approximate.

$\int f Y / u\left(y i / u_i\right) f u\left(u_i\right) d u_i \approx \sum_{i=1}^Q w_i^*\left(\prod_{j=1}^n f Y / u\left(y_{i j} Z_i^*\right)\right)$          (36)

u_i represent the size of the cluster i, y_ij shows j-th element of cluster i with have the adaptive weights given by:

$W_i^*=\sqrt{2 \sigma_i} w_i \exp \left(Z_i^2\right) \phi\left(Z_i^* ; 0,1\right)$          (37)

where, σ_i represent approximation for $\sigma^{-1} u_i \sim N(0,1)$. Hence, we have linear explanatory variables $x_i^{\prime}+\sigma Z_i^*$ for $f Y / u\left(y_{i j} Z_i^*\right)$.

2.8 Data description

National Health Insurance Scheme (NHIS) data of three health facilities in Ota, Ogun State, Nigeria were collected. Medicare Private Hospital and Ogun State Government Hospital NHIS records of patients were collated during the visitation. Records of follow-up status of individual patients were picked as dependent variables while sex, age, number of diagnoses and blood group of concerned patients were categorized as predictor variables.

This consists of clinical data of 1500 patients visiting the facilities between July 2016 and July 2017. The data set was collated to examine the impact of predictors on a binary response variable (follow-up). Follow-up, Sex (gender), Age is biological age, Ndiagnosis stands for the number of diagnoses, Bgroup stands for blood group, coded as (A=1, B=2, AB=2 and O=4), Gnotyp stands for genotype, coded as (AA=1, AS=2, SS=3), Sstatus stands for smoking status, coded as (smoke=1, not smoke=0). For this research, the binary response was based on whether a patient was on follow-up or not following advice given by the physician. The descriptive statistics can be found in the Table A1.

The area of application centers around doctor's follow-up on Patient’s based on some variables. We illustrate applications of different GLMMs estimation techniques in analyzing the clinical binary response data with an inclusion of random effects. Leveraging on the model coefficients, factors associated with medical follow ups may influence frequent medical diagnoses. Four GLMM techniques were fit to the data and allowed the intercept to vary across the clusters. The model parameters were estimated using penalized quasi-likelihood (PQL), Adaptive Gaussian Hermite quadrature (QUAD), and Laplace approximation techniques.

Appropriate modeling for medical data with binary response has some obvious and important implications for patient diagnoses in Health care system. In this paper, the observation of multiple correlation across covariates is an indication of the presence of random effects because of repeated measures on predictor variables. This made it extremely difficult to model such data using traditional generalized linear model, which does not accommodate random effects. The choice of GLMMs was to take care of multi-collinearity effects across covariates to avoid model misspecification.

The statistical inference drawn from this study shows the fitness and superiority of a gauss-Hermite quadrature estimation technique over other estimation techniques based on some selected model performance metrics for the real data set.

In addition, the study was not only meant to establish the fitness or superiority of the model (glmmGHQ) over others but to also validate effects of clinical follow-up on Patient's medical diagnoses over a period of twelve months leveraging on the fitness of the model (glmmGHQ). Table 4 shows increase in medical diagnoses for more follow-up status of individuals. Using the statistics, for every unit increase in follow-up status, there is a corresponding increase of 0.05471 of number of diagnoses. It is also an observation from Table 4 that with a unit increase in the follow-up status, there is a decrease of 0.01473 in smoking habit of patients. The follow-up effects also trigger an increase of 0.03497 factor in older individual patients. Older patients tend to have more diagnoses than other age categories according to the study results.

Leveraging research results, patients with prompt response to clinical follow-up had a corresponding increasing number of diagnoses, which helped in early detection of ailments that could have been deadly. The result analysis showed patients with follow-up status with a decrease in smoking attributes, which might be a result of clinical counselling during medical appointments.

Based on the research results, the gauss-hermite Quadrature technique might be preferred in estimating model parameters for generalized linear mixed models with binary responses better model results. Investigating effects of risk factors in health care systems may also be carried out using a gauss-Hermite Quadrature technique.