Generalized Linear Mixed Model Estimation Methods in Sample Survey Statistics
I. Introduction to Generalized Linear Mixed Models in the Survey Context
Generalized Linear Mixed Models (GLMMs) have emerged as a sophisticated and flexible class of statistical models, finding increasing application in the analysis of complex data structures prevalent in sample survey statistics. Their utility stems from an ability to handle non-normally distributed response variables and accommodate correlations among observations, features that are common in survey research. Understanding GLMMs requires appreciating their evolutionary path from simpler linear model frameworks and recognizing the unique analytical demands posed by sample survey data.
A. Defining Generalized Linear Mixed Models (GLMMs): An Evolution of Linear Models
The core concept of a Generalized Linear Mixed Model (GLMM) is an extension of both Generalized Linear Models (GLMs) and Linear Mixed Models (LMMs).1 GLMMs are specifically designed for analytical situations where the outcome variable does not follow a normal distribution and/or where the observations exhibit dependencies, often due to hierarchical or clustered sampling designs inherent in many surveys.3
The progression from GLMs to GLMMs involves two key generalizations. Firstly, GLMs extend ordinary linear regression by permitting the response variable to follow a distribution from the exponential family (e.g., binomial for binary outcomes, Poisson for count data) rather than being restricted to a normal distribution. Furthermore, GLMs introduce a link function, g(.), which relates the mean of the response variable, E(Yi​)=μi​, to a linear predictor of fixed effects, ηi​=Xi​β, such that g(μi​)=ηi​.5 GLMMs build upon this foundation by incorporating random effects, say Zi​ui​, into the linear predictor, so that g(μi​)=Xi​β+Zi​ui​. These random effects are instrumental in modeling correlation and unobserved heterogeneity among observations.1 As such, GLMMs are considered a general modeling framework, with linear models (LM), LMMs, and GLMs being special cases.2
Secondly, the transition from LMMs to GLMMs involves broadening the scope of response variable distributions. LMMs are tailored for normally distributed (Gaussian) response data where both fixed and random effects are present, allowing for the modeling of various correlation structures.2 GLMMs extend this capability to accommodate non-normal response distributions, making them more versatile for the diverse types of outcomes encountered in survey research.4
The power of GLMMs lies in their capacity to simultaneously consider fixed effects, which represent systematic influences such as treatment effects or demographic characteristics, and random effects, which capture sources of variability such as differences between survey clusters (e.g., schools, hospitals) or the correlation inherent in repeated measurements on the same individuals over time.4 This dual capacity makes GLMMs exceptionally well-suited for analyzing the complex, often hierarchical, data structures that arise from sample surveys.3 The development of GLMMs represents a significant convergence of methodologies: the flexibility of GLMs in handling various data types is combined with the LMM framework's ability to model complex dependency structures. When applied to survey data, this signifies a crucial advancement, enabling these sophisticated models to be adapted to the practical realities and inherent complexities of large-scale population data collection methods.
B. The Unique Challenges of Sample Survey Data
Sample survey data frequently present unique challenges that distinguish them from data collected under simpler experimental designs. "Complex survey data" are typically characterized by one or more of the following features: stratification, clustering, unequal probabilities of selection (resulting in the need for survey weights), and finite population corrections.9 These design elements are deliberately employed to enhance statistical efficiency, control costs, or target specific subpopulations, but they result in deviations from the simple random sampling (SRS) assumption that underpins many classical statistical techniques.9
The impact of these features on statistical analysis is profound:
●	Clustering: Observations within the same cluster (e.g., individuals in the same household, students in the same school, or enumeration areas known as Primary Sampling Units - PSUs) often exhibit greater similarity to each other than to observations from different clusters. This intra-cluster correlation violates the assumption of independence that is fundamental to many standard statistical procedures.10
●	Stratification: The division of the population into homogeneous subgroups (strata) prior to sampling can improve the precision of estimates for the total population and for these subgroups. However, this stratification must be accounted for in variance estimation to reflect the design properly.
●	Unequal Weights: Survey weights are assigned to sample units to compensate for unequal probabilities of selection, to adjust for nonresponse, and to calibrate sample estimates to known population totals. If these weights are ignored in the analysis, the resulting estimates of population parameters can be severely biased.10
Failure to appropriately account for these complex design features during analysis can lead to erroneous conclusions. Specifically, it can result in biased parameter estimates, incorrect standard errors (often underestimated), and consequently, invalid statistical inferences regarding p-values and confidence intervals.10 Standard statistical formulas and software procedures that assume SRS are generally not appropriate for complex survey data.9 The move towards more sophisticated models like GLMMs for analyzing survey data is therefore not merely an academic pursuit of statistical elegance, but a practical necessity driven by the imperative to obtain valid and nuanced inferences that accurately reflect the population structure and the intricacies of the sampling process. The alternative—ignoring these critical design features—can lead to demonstrably misleading or incorrect conclusions, as evidenced by studies showing altered inferences when survey weights are omitted.10
C. Why GLMMs are Suited for Analyzing Complex Survey Data
GLMMs offer a powerful and appropriate framework for analyzing complex survey data due to several inherent characteristics that align well with the challenges posed by such data:
●	Modeling Hierarchical Structures: Survey data frequently possess a hierarchical or multilevel nature. For instance, individuals may be nested within households, households within PSUs, and PSUs within regions. GLMMs, as a class of multilevel models, are naturally equipped to handle such nested structures through the incorporation of random effects at different levels of the hierarchy.7 The term "multilevel model" describes a category of statistical models, and GLMMs are a specific type particularly well-suited for these hierarchical datasets encountered in surveys.13
●	Handling Correlated Data: The random effects component in GLMMs provides an explicit mechanism for modeling the correlation that arises between observations within the same cluster (e.g., due to shared environments or characteristics) or from repeated measurements on the same survey unit over time.1 This directly addresses the violation of independence often found in survey data.
●	Flexibility for Non-Normal Outcomes: Many outcomes measured in surveys are not normally distributed. Examples include binary variables (e.g., employment status, presence of a health condition), count variables (e.g., number of hospital visits, frequency of a behavior), or ordinal responses. GLMMs accommodate a wide range of response distributions from the exponential family, allowing the analyst to choose a distribution appropriate for the specific outcome variable.2 Data encountered in social sciences and health research often deviate significantly from normality, making the flexibility of GLMMs particularly valuable.16
●	Model-Based Approach to Design Features: GLMMs provide a model-based approach to analysis. In this paradigm, the dependencies and complexities introduced by the survey sampling design are often accounted for through the careful specification of the model, particularly the structure of the random effects.10 For example, clustering can be directly modeled by including random effects for the clusters.
The following table provides a concise comparison of GLMs, LMMs, and GLMMs, highlighting the features that make GLMMs a comprehensive tool.
Table 1: Comparison of GLM, LMM, and GLMM

Feature	General Linear Model (GLM)	Linear Mixed Model (LMM)	Generalized Linear Mixed Model (GLMM)	
Response Variable Distribution	Exponential family (e.g., Normal, Binomial, Poisson) 5	Normal (Gaussian) only 2	Exponential family (e.g., Normal, Binomial, Poisson, Gamma) 1	
Link Function	Various (e.g., Identity, Logit, Log) 5	Identity only	Various (e.g., Identity, Logit, Log) 1	
Nature of Effects	Fixed effects only	Fixed and Random effects 17	Fixed and Random effects 1	
Handles Correlated Data	No (assumes independence)	Yes, via random effects and/or correlated error structures 2	Yes, primarily via random effects 1	
Typical Application	Regression for non-normal independent data	Regression for normal, correlated/hierarchical data	Regression for non-normal, correlated/hierarchical data 1	
This comparative overview underscores how GLMMs integrate the distributional flexibility of GLMs with the capacity of LMMs to model dependencies through random effects, establishing them as a versatile and powerful analytical tool for the multifaceted data typically generated by sample surveys.
II. Core Components of Generalized Linear Mixed Models
Generalized Linear Mixed Models are composed of several key elements that work in concert to define the relationship between predictors and a (potentially non-normal and correlated) response variable. These components include the linear predictor incorporating fixed and random effects, a link function, and a conditional distribution for the response.
A. The Linear Predictor: Fixed Effects and Covariates
At the heart of a GLMM is the linear predictor, denoted as η. This component is a linear combination of fixed effects parameters and their corresponding explanatory variables (covariates).2 If β represents the vector of unknown fixed-effects parameters and X is the matrix of covariates, the fixed-effects part of the linear predictor is Xβ. This systematic component reflects the influence of observed covariates on the response.5
Fixed effects represent influences that are considered constant across all individuals or groups in the population, or they are the primary factors whose effects are of direct research interest, such as the impact of a specific policy, an intervention (treatment effect), or demographic characteristics.7 In the context of GLMMs, these fixed effects describe population-average relationships, assuming these relationships are consistent across the different levels of random factors.7
B. Incorporating Random Effects: Modeling Hierarchies and Dependencies
A defining feature of GLMMs, distinguishing them from GLMs, is the inclusion of random effects in the linear predictor. If u is a vector of random effects and Z is the design matrix for these random effects, the linear predictor becomes η=Xβ+Zu. These random effects are unobserved variables that are assumed to be drawn from a probability distribution, typically a normal distribution with a mean of zero and a variance-covariance matrix G (whose parameters, the variance components, are estimated from the data).1
Random effects serve several crucial purposes in GLMMs:
1.	Modeling Correlation: They explicitly account for the correlation among observations that are grouped or clustered, such as repeated measurements on the same individual over time, or individuals sampled within the same household or geographic area.1 The shared random effect for units within the same cluster induces this correlation.
2.	Accounting for Unobserved Heterogeneity: Random effects capture variability between clusters or individuals that is not explained by the observed fixed effects covariates.2 This is particularly important in survey data where clusters (e.g., PSUs) can vary substantially in unmeasured ways.
3.	Addressing Overdispersion: In the analysis of count or proportion data, the observed variance sometimes exceeds what is predicted by simpler distributions like the Poisson or binomial. This phenomenon, known as overdispersion, can be modeled by including random effects (e.g., an observation-level random effect).1
4.	Enabling Subject-Specific Inference: While fixed effects provide population-averaged estimates, random effects allow for predictions or inferences specific to particular subjects or clusters (e.g., estimating a particular school's effect).20
The decision to treat an effect as fixed or random is a critical modeling choice. Generally, effects are considered random if their levels are thought of as a sample from a larger population of levels, and the primary interest is in the overall variability attributable to these levels rather than the effect of each specific level.1 For instance, if survey PSUs are randomly selected to represent a larger set of PSUs, their effects might be modeled as random. In contrast, if the specific levels of a factor are of direct interest and are not considered a sample (e.g., specific treatment categories), they are modeled as fixed. While statistical tests like the Hausman test can sometimes guide this decision by comparing fixed and random effects models, practical considerations and the research question often play a more significant role.19 The specification and estimation of these random effects are central to the utility of mixed models for the correlated, hierarchical data structures frequently encountered in surveys, as they provide the mechanism to model both within-cluster homogeneity and between-cluster heterogeneity.
C. Link Functions: Connecting the Linear Predictor to the Response Mean
The link function, g(.), provides the crucial connection between the linear predictor, η, and the conditional mean of the response variable, μ=E(Y∣u).2 The linear predictor η=Xβ+Zu can take any real value, but the mean of the response variable is often constrained (e.g., probabilities μ for binomial data must be between 0 and 1; counts μ for Poisson data must be non-negative). The link function transforms the mean to ensure it maps to the full range of the linear predictor: g(μ)=η. Conversely, the inverse link function, g−1(.), maps the linear predictor back to the original scale of the mean: μ=g−1(η).
Common link functions used in GLMMs include 2:
●	Identity link: μ=η. Used for normally distributed responses where the mean itself is modeled linearly. This is the link function in standard LMMs.
●	Logit link: log(1−μμ​)=η. Used for binomial (binary or proportional) data. The quantity 1−μμ​ is the odds.
●	Probit link: Φ−1(μ)=η, where Φ−1 is the inverse cumulative distribution function (CDF) of the standard normal distribution. Also used for binary data, it assumes an underlying continuous latent variable that is normally distributed.6
●	Log link: log(μ)=η. Used for Poisson or Negative Binomial distributed count data.
●	Cumulative Logit link: Used for ordered categorical (multinomial) responses.2
A canonical link function is one that transforms the mean of the response to the natural parameter of its assumed exponential family distribution.5 For example, the logit link is canonical for the binomial distribution, and the log link is canonical for the Poisson distribution. Canonical links often have desirable statistical properties and can simplify estimation and interpretation.6 The interplay between the chosen link function and the assumed distribution for the response variable is fundamental; an appropriate pairing ensures that the model's predictions are valid for the scale of the response variable and that parameters are interpretable.
D. Conditional Distributions for Various Response Types
In a GLMM, the response variable Y, conditional on the random effects u, is assumed to follow a probability distribution from the exponential family.1 This family is a broad class of distributions that includes the normal, binomial, Poisson, gamma, and negative binomial distributions, among others. The choice of the specific distribution is dictated by the nature of the response variable being analyzed.2
Common conditional distributions and their typical applications in survey contexts include:
●	Normal (Gaussian) Distribution: For continuous response variables that are, or can be transformed to be, approximately normally distributed and are unbounded.2 If random effects are included, this leads to an LMM, a special case of GLMM.
●	Binomial Distribution: For response variables that represent the number of "successes" in a fixed number of trials (e.g., number of household members with health insurance out of total members) or for binary (0/1) outcomes (e.g., an individual is employed or not).2 The Bernoulli distribution is a special case of the binomial where the number of trials is one.2
●	Poisson Distribution: For response variables that are counts of events occurring within a fixed period or region (e.g., number of doctor visits in a year, number of crimes in a PSU).2
●	Negative Binomial Distribution: Also for count data, particularly when there is overdispersion, meaning the variance of the counts is greater than the mean. This is common in survey data where unobserved heterogeneity can inflate variance.2
●	Gamma Distribution: For continuous, positively skewed response variables, such as costs or durations.5
●	Beta Distribution: For continuous response variables that are proportions or percentages constrained between 0 and 1 (e.g., proportion of a budget spent on housing).2
A key assumption is that observations are independent conditional on the random effects.7 The random effects structure is what explicitly models the dependence among observations within the same group or cluster. Unobserved heterogeneity between individuals or clusters in survey data often leads to overdispersion in count or binomial outcomes. Random effects in GLMMs provide a natural mechanism to capture this extra-Poisson or extra-Binomial variation, leading to more accurate standard errors and more reliable inferences than would be obtained from a simpler GLM that ignores this feature.1
The following table summarizes common distributions, their typical survey data applications, and associated link functions.
Table 2: Common Distributions, Typical Data Types, and Canonical/Common Link Functions in GLMMs

Distribution	Typical Survey Data Type	Data Range	Canonical Link Function	Common Alternative Link(s)	
Normal (Gaussian)	Continuous measurements (e.g., income, blood pressure)	(−∞,∞)	Identity 5	None	
Binomial	Number of successes in n trials (e.g., individuals with attribute in household)	0,1,...,n	Logit 5	Probit, Cloglog 6	
Bernoulli	Binary outcome (e.g., yes/no, employed/unemployed)	0,1	Logit 5	Probit, Cloglog 6	
Poisson	Counts of events (e.g., number of health visits, defects per item)	0,1,2,...	Log 5	Identity (less common)	
Negative Binomial	Counts of events with overdispersion (e.g., accident counts, disease incidence)	0,1,2,...	Log 2	Identity (less common)	
Gamma	Continuous, positive, skewed data (e.g., insurance claims, duration data)	(0,∞)	Inverse (1/μ)	Log, Identity	
Beta	Continuous proportions/percentages (e.g., test scores as proportion, recovery rates)	(0,1)	Logit 2	Probit, Cloglog	
III. Estimation Methods for GLMM Parameters
Estimating the parameters of a GLMM—both the fixed effects (β) and the variance components of the random effects (elements of G and the residual variance/dispersion parameter)—is considerably more complex than for LMMs or GLMs. This complexity arises primarily because the likelihood function for a GLMM involves an integral over the distribution of the random effects, which typically does not have a closed-form analytical solution.1 Consequently, various numerical approximation and computational techniques have been developed.
A. Likelihood-Based Approaches
Likelihood-based methods aim to find parameter estimates that maximize the probability of observing the sample data. However, the direct evaluation of the GLMM likelihood, L(β,G∣y)=∫∏i​f(yi​∣Xi​,Zi​,β,ui​)p(ui​∣G)dui​, is challenging due to this integration over the random effects ui​.18
1. Maximum Likelihood (ML) Estimation
ML estimation seeks parameter values that maximize the (log-)likelihood function.18 For GLMMs, this maximization requires numerical techniques to approximate the integral and then find the maximum of this approximated likelihood.1 While ML estimates for fixed effects are generally consistent, ML estimates for variance components can be biased downwards, particularly in samples with few clusters or groups, as they do not account for the degrees of freedom lost in estimating the fixed effects.
2. Restricted Maximum Likelihood (REML) Estimation
REML estimation is designed to reduce the bias observed in ML estimates of variance components.1 It achieves this by maximizing a modified likelihood function that is integrated over the fixed effects, effectively focusing on the part of the data that provides information about the variance components, independent of the fixed effects estimates.21 REML is widely preferred for estimating variance components in LMMs.21 Its application to GLMMs is more intricate due to the non-linear link function and non-normal data. Adaptations, often involving pseudo-likelihood concepts (such as Residual Pseudo-Likelihood, RSPL), are employed to extend the REML principle to GLMMs.2 For instance, the SAS PROC GLIMMIX procedure uses RSPL by default when fitting GLMMs to non-normal data.2
3. Numerical Integration Techniques (Approximating the Likelihood)
Given the intractability of the GLMM likelihood integral, several numerical methods are used for its approximation:
●	Adaptive Gauss-Hermite Quadrature (AGHQ): This is a highly accurate numerical integration technique, particularly effective when the number of random effects per cluster is small (e.g., random intercepts or random intercepts and one random slope).8 AGHQ approximates the integral by a weighted sum of the integrand evaluated at specific points (nodes). The "adaptive" nature refers to the centering and scaling of these nodes based on the mode and curvature of the integrand for each cluster, which improves accuracy.23 The number of quadrature points used (e.g., nAGQ parameter in R's glmer) affects accuracy and computational time; more points yield greater accuracy but increase computation.23
●	Laplace Approximation: This method provides a faster but generally less accurate approximation to the likelihood integral compared to AGHQ, especially if the true conditional distribution of random effects given data is far from normal or if variance components are large.8 It approximates the integrand with a Gaussian function centered at its mode. The Laplace approximation is often the default or only feasible option for models with multiple correlated random effects or crossed random effects, where AGHQ becomes computationally prohibitive.24
●	Penalized Quasi-Likelihood (PQL) and Marginal Quasi-Likelihood (MQL): These are older, computationally faster methods based on iterative linearization of the model and fitting a working LMM at each step.22 PQL tends to perform better than MQL but can still produce biased estimates, especially for binary responses or when cluster sizes or event counts per cluster are small.22 They are generally considered less accurate than AGHQ or Laplace for many GLMM scenarios.
●	Integrated Nested Laplace Approximations (INLA): While often discussed in a Bayesian context (see below), INLA is a sophisticated numerical approximation method that can also be viewed as an advanced likelihood approximation technique for latent Gaussian models, a class that includes GLMMs.8 It directly approximates posterior marginal distributions, avoiding the need for MCMC simulations.
The estimation of GLMM parameters is thus a cascade of approximations. Nearly all methods grapple with the intractable likelihood, employing numerical strategies (quadrature, Laplace), methodological simplifications (PQL, pseudo-likelihood), or sampling-based approaches (MCMC). A clear understanding of these approximations' nature and limitations is vital for selecting an appropriate method and correctly interpreting the results.
B. Pseudo-Likelihood Approaches for Complex Surveys
When GLMMs are applied to data from complex sample surveys, standard likelihood-based estimation methods face an additional challenge: they typically assume simple random sampling (conditional on covariates). Complex survey designs, with features like stratification, clustering, and unequal probabilities of selection (requiring survey weights), violate this assumption. Pseudo-likelihood (PL) methods have been developed to adapt likelihood-based estimation to the survey context by modifying the likelihood function to incorporate survey weights.12 The goal is to obtain estimators of model parameters that are approximately unbiased with respect to the survey design.
The PL approach for survey data often involves maximizing a weighted sum of individual or cluster-level contributions to a likelihood-like function, where the weights are the survey weights.20 This effectively re-weights the sample to represent the target finite population. For instance, in SAS PROC GLIMMIX, observation-level weights can be specified using the OBSWEIGHT option in the MODEL statement, and cluster-level weights via the WEIGHT option in the RANDOM statement.29
Many PL estimation procedures are iterative. For example, some methods involve iteratively fitting a weighted Gaussian LMM to a linearized or "working" dependent variable, where the survey weights are incorporated into the weighting matrix at each iteration.20 Residual Pseudo-Likelihood (RSPL) is a REML-like variant of PL often used for estimating variance components in this context.2
A critical component of PL estimation with survey data is the subsequent estimation of variances for the parameter estimates. Standard errors derived from the inverse of the observed or expected Fisher information matrix of the pseudo-likelihood are generally incorrect because they do not account for the full impact of the survey design (especially clustering and stratification) or the fact that the PL is an approximation. Therefore, robust variance estimation techniques, such as Taylor series linearization (sandwich estimators) or replication methods (e.g., bootstrap, jackknife), are commonly employed in conjunction with PL to obtain valid standard errors and confidence intervals.11 Pseudo-likelihood thus emerges as a pragmatic and vital adaptation for survey data, directly addressing the critical issue of unequal sampling weights by embedding them within the estimation objective function. While other methods like ML or Bayesian approaches might find it more challenging to incorporate weights or may necessitate intricate modeling of the sampling process itself, PL provides a more direct, albeit approximate, pathway to design-adjusted estimates.
C. Bayesian Estimation Methods
Bayesian methods offer an alternative paradigm for GLMM parameter estimation. This approach combines the likelihood of the data (given the parameters) with prior distributions specified for all model parameters (fixed effects, variance components) to derive a joint posterior distribution of the parameters.8 Inferences are then based on this posterior distribution.
Markov Chain Monte Carlo (MCMC) Methods:
Due to the complexity of the posterior distribution in GLMMs (which, like the likelihood, involves intractable integrals), MCMC methods are the traditional computational workhorse for Bayesian GLMMs.8 Algorithms like Gibbs sampling and Metropolis-Hastings are used to generate a sequence of samples from the posterior distribution. After a burn-in period and assessment of convergence, these samples can be used to estimate posterior means, medians, credible intervals, and other summaries for the parameters.
●	Advantages: MCMC is highly flexible, capable of handling very complex models, various prior specifications, and an arbitrary number of random effects. It is generally considered to produce accurate characterizations of the posterior distribution given sufficient run time.22
●	Disadvantages: MCMC methods are computationally intensive and can be very slow to converge, especially for large datasets or complex models.22 Proper implementation requires careful attention to convergence diagnostics, choice of priors (which can be influential, especially for variance components), and thinning of chains.
Integrated Nested Laplace Approximations (INLA):
INLA has gained prominence as a faster alternative to MCMC for Bayesian inference in latent Gaussian models, a class that encompasses GLMMs.8 INLA uses a series of Laplace approximations to directly compute highly accurate approximations to the posterior marginal distributions of the parameters, thereby avoiding stochastic MCMC simulations.
●	Advantages: Significantly faster than MCMC for many GLMMs, making Bayesian analysis more feasible for larger datasets or more complex random effect structures where MCMC would be prohibitive.26
●	Limitations: While powerful, INLA may be less flexible than MCMC for extremely complex or non-standard prior distributions or for models that do not fit neatly into the latent Gaussian framework.
Challenges in Bayesian GLMMs for Surveys:
Applying Bayesian GLMMs to complex survey data introduces specific challenges:
●	Prior Specification: Defining appropriate and non-informative (or weakly informative) prior distributions, particularly for variance components, can be difficult and may substantially influence posterior inferences, especially with limited data at higher levels of a hierarchy.26
●	Incorporating Survey Weights: The integration of survey weights into a fully Bayesian framework is not as straightforward as in PL methods. Common approaches include using a weighted likelihood in the Bayesian update (leading to a pseudo-Bayesian analysis), attempting to model the sampling process explicitly (which is often very complex and requires information about the design not always available), or using methods like Bayesian pseudo-empirical likelihood.27
●	Computational Cost: Even with INLA, fitting Bayesian GLMMs to very large survey datasets can be computationally demanding. MCMC methods are often impractical for such datasets without significant computational resources or algorithmic advances.30
The landscape of Bayesian GLMMs is continually evolving. While MCMC remains a foundational technique, the speed and practicality of INLA have made Bayesian approaches more accessible. For survey data specifically, the development of robust and computationally efficient methods that fully integrate complex design features within a coherent Bayesian framework is an active area of research, with hybrid approaches sometimes borrowing concepts from pseudo-likelihood to bridge the gap.27
The estimation process in GLMMs invariably involves tackling both fixed effects parameters and variance components. Methods like REML (and its pseudo-likelihood counterparts such as RSPL) are specifically tailored to enhance the estimation of variance components.1 The accuracy of these variance component estimates is paramount, as they directly influence the standard errors of the fixed effects estimates and the empirical Bayes predictions of the random effects, thereby impacting the entire inferential structure of the mixed model.
The following table provides a comparative overview of the main estimation methods for GLMMs, particularly highlighting their relevance and challenges in the context of sample survey data.
Table 3: Overview of Estimation Methods for GLMMs in Survey Context

Estimation Method	Brief Description	Primary Likelihood Approximated/Targeted	How Survey Weights are Typically Handled	Common Numerical Techniques Used	Key Advantages	Key Limitations/Challenges in Survey Context	
Maximum Likelihood (ML)	Finds parameters maximizing the likelihood of observed data.	Full marginal likelihood (approximated).	Not directly; assumes SRS. May be used within pseudo-likelihood frameworks if weights are part of the model structure (e.g., inverse probability weighting as a covariate, rarely optimal).	Numerical integration (AGHQ, Laplace).8	Standard likelihood theory applies (asymptotically).	Biased variance components. Computationally intensive. Does not inherently account for complex survey design.1	
Restricted Maximum Likelihood (REML)	Maximizes a likelihood function independent of fixed effects to reduce bias in variance components.	Restricted marginal likelihood (approximated).	Similar to ML; primarily for variance components. Adaptations (e.g., RSPL) used in PL.	Numerical integration (AGHQ, Laplace); iterative algorithms.	Less biased variance components than ML.21	More complex for GLMMs than LMMs. Does not inherently account for survey design beyond what PL adaptations offer.2	
Pseudo-Likelihood (PL) / Penalized Quasi-Likelihood (PQL)	Maximizes a modified likelihood or quasi-likelihood function incorporating survey weights.	Approximates marginal likelihood or uses quasi-likelihood principles.	Directly incorporated into the objective function (e.g., weighted sum of log-likelihood contributions).12	Iteratively Reweighted Least Squares (IRLS); fitting weighted Gaussian LMMs iteratively.20	Directly accounts for survey weights. Computationally less intensive than full Bayesian MCMC for large surveys. Robust variance estimation available.28	Estimates are approximate. Properties depend on approximation quality. Potential bias for some data types (e.g., binary with few events per cluster for PQL).22	
Bayesian MCMC	Samples from the posterior distribution of parameters using Markov Chain Monte Carlo algorithms.	Full posterior distribution (derived from likelihood and priors).	Complex; may use weighted likelihood (pseudo-Bayesian), model the sampling process, or use Bayesian pseudo-empirical likelihood.27	Gibbs sampling, Metropolis-Hastings.8	Highly flexible for complex models/priors. Provides full posterior distribution for inference. Accurate with sufficient computation.22	Computationally very intensive, slow convergence. Prior specification can be challenging. Incorporating survey weights correctly is difficult.22	
Bayesian INLA	Deterministic approximation to posterior marginals for latent Gaussian models.	Approximates posterior marginals.	Similar to MCMC; depends on how likelihood is formulated with weights.	Integrated Nested Laplace Approximations.8	Much faster than MCMC for many models. Provides accurate approximations.26	Less flexible than MCMC for some priors/models. Correct handling of survey weights still an area of development.	
IV. Addressing Complex Survey Design Features in GLMM Estimation
The validity of inferences drawn from GLMMs applied to sample survey data hinges critically on the appropriate incorporation of key design features: survey weights, clustering, and stratification. These elements are not mere nuisances but reflect the intricate process by which the sample was drawn from the finite population and must be addressed to ensure that estimates are representative and standard errors are correctly calculated.
A. Incorporating Survey Weights
Survey weights are fundamental to producing estimates that accurately reflect the target population from which the sample was drawn. These weights are typically constructed to account for differential probabilities of selection, to adjust for nonresponse at various stages, and to calibrate sample totals to known population benchmarks.9 Omitting survey weights when they are informative—that is, when the weights are correlated with the outcome variable even after conditioning on model covariates—will generally lead to biased estimates of population parameters.10
Methods of Incorporation in GLMMs:
The predominant method for incorporating survey weights into GLMMs is through pseudo-likelihood estimation. In this approach, the standard log-likelihood function is modified so that the contribution of each sampled unit (or cluster) is weighted by its survey weight.12 This effectively creates a "pseudo-population" that the model is fitted to, aiming to achieve design-consistency for the estimated parameters. For example, software like SAS PROC GLIMMIX allows specification of weights in the MODEL statement (e.g., OBSWEIGHT=wgt_var) for observation-level effects or, more commonly in mixed models, in the RANDOM statement (e.g., WEIGHT=cluster_wgt_var) to apply weights at a specific level of the hierarchy.29
In iterative fitting procedures common to many pseudo-likelihood algorithms (such as those based on iteratively reweighted least squares or PQL), the survey weights are typically incorporated into the weighting matrix used at each step of the algorithm.20
Incorporating weights into Bayesian GLMMs is more complex. While some approaches involve using a weighted likelihood within the Bayesian framework (a form of pseudo-Bayesian inference) 27, or attempting to model the sampling mechanism itself (often highly challenging), this remains an area of active research.
Weighting at Different Levels:
Complex surveys may provide weights applicable at different levels of a data hierarchy (e.g., person weights, household weights, PSU weights). When fitting multilevel GLMMs to such data, it is often necessary to incorporate these weights at the corresponding levels of the model.10 For instance, if a two-level model has students nested within schools, and separate weights are available for students and schools, the estimation procedure should ideally accommodate both sets of weights. Software implementations vary in their capacity to handle multilevel weights directly within the GLMM framework.
Scaling of Weights:
Raw survey weights can sometimes be highly variable or take on very large values, which can lead to numerical instability during model estimation or produce inefficient estimates, particularly for variance components. Scaling the weights—for example, by normalizing them to sum to the effective sample size within each cluster, or to have an average of one across the sample—is a common practice to mitigate these issues.10 However, the method of scaling can itself influence the estimates, especially of variance components, and requires careful consideration.28 There is no universally optimal scaling strategy, and the choice may depend on the specific survey design and analytical goals. The decision of at which level to apply weights in a multilevel GLMM and whether/how to scale them are critical choices that demand careful thought from the analyst, as they can significantly affect both parameter estimates and computational performance.
Impact on Estimates:
Properly incorporated survey weights aim to yield approximately unbiased and consistent estimates of finite population parameters. The use of weights, along with other design features like clustering and stratification, also necessitates robust methods for variance estimation to ensure that standard errors, confidence intervals, and hypothesis tests are valid.11
B. Accounting for Clustering
Clustering is a common feature of survey designs, where ultimate sampling units are selected in groups (e.g., individuals within households, households within PSUs) primarily for reasons of cost-effectiveness and operational feasibility.9 This design choice typically leads to correlation among observations within the same cluster, as units within a cluster tend to share common characteristics or environmental influences.
The primary mechanism by which GLMMs account for such clustering is through the inclusion of random effects specific to the clusters.1 For example, in a survey of students clustered within schools, a random intercept for schools would capture unobserved school-level factors that affect all students in that school similarly, thereby inducing a correlation among their outcomes. The variance of this random effect quantifies the extent of between-cluster heterogeneity.
In statistical software, specifying clustering typically involves identifying the variable that defines the clusters in the model syntax. For instance, in SAS PROC GLIMMIX, the SUBJECT= option within the RANDOM statement is used to define the grouping structure for random effects.2 In R's lme4 package, this is done using the (1 | cluster_variable) syntax in the model formula.
Failure to account for clustering when it is present and influential generally leads to an underestimation of standard errors for fixed effects estimates. This, in turn, can result in overly narrow confidence intervals and an inflated probability of Type I errors (false positives) in hypothesis testing.
C. Handling Stratification
Stratification involves partitioning the population into mutually exclusive and exhaustive subgroups (strata) based on known characteristics (e.g., geographic region, urban/rural status) and then drawing independent samples from each stratum.9 This technique is often used to improve sampling efficiency, ensure representation of key subpopulations, and allow for the production of stratum-specific estimates.
The approach to handling stratification in GLMMs can be less direct than for clustering:
1.	Implicitly through PSU Specification: In many multistage designs, PSUs are sampled within strata. If the GLMM includes random effects for PSUs (often as the highest level of clustering in the model), and variance estimation correctly identifies these PSUs as nested within strata (or if strata are the highest level units), then stratification is implicitly accounted for in the variance estimation process, particularly when using Taylor series linearization or replication methods.11
2.	Fixed Effects for Strata: If the number of strata is small and stratum-specific effects are of substantive interest or need to be controlled for, stratum indicator variables can be included as fixed effects in the model. However, this approach can become unwieldy and lead to numerical instability if the number of strata is large.32
3.	Robust Variance Estimation: Design-based variance estimators (e.g., sandwich estimators) used in conjunction with pseudo-likelihood methods can be specifically formulated to incorporate stratification information in the calculation of standard errors for the GLMM parameters.28 This is often the most common and practical approach.
While GLMMs inherently address clustering via their random effects structure, the handling of stratification and the impact of weighting on variance often relies on robust variance estimation techniques applied post-hoc or integrated with pseudo-likelihood methods. This suggests a practical division of labor: the model's random effects structure captures the correlation from clustering, while specialized variance estimation methods address the broader survey design complexities, including stratification and weighting.
It is also important to recognize the interdependence of these design features. PSUs are often defined within strata, and survey weights are calculated based on selection probabilities at each stage of a potentially stratified and clustered sampling process. Analytical strategies must therefore consider these interconnections, for example, by correctly specifying the nesting of PSUs within strata for variance estimation or by using weights that accurately reflect all stages of the complex design.9
V. Software Implementation for GLMMs with Survey Data
The practical application of GLMMs to complex survey data relies heavily on the capabilities of statistical software packages. Several major packages offer functionalities for fitting GLMMs, though their ability to seamlessly integrate all aspects of complex survey designs (particularly multilevel weights and robust variance estimation for all design features) can vary.
A. Overview of Common Statistical Packages
●	SAS: A widely adopted software in official statistics, research institutions, and industry. PROC GLIMMIX is the primary SAS procedure for fitting GLMMs, offering extensive control over model specification, including distributions, link functions, random effects structures, and estimation methods.32 While SAS has dedicated survey analysis procedures like PROC SURVEYLOGISTIC and PROC SURVEYREG for GLMs, these do not directly handle mixed models. Thus, fitting GLMMs to survey data in SAS typically involves using PROC GLIMMIX with careful specification of weights and requests for empirical (sandwich) variance estimators.12 PROC NLMIXED can also fit some GLMMs but requires more manual specification.
●	R: An open-source statistical programming language and environment, R boasts a vast collection of packages contributed by its user community. For GLMMs, prominent packages include lme4 (with the glmer function), glmmTMB (offering extended families and zero-inflation), and MCMCglmm (for Bayesian GLMMs using MCMC).8 The survey package (with svyglm) is the standard for analyzing complex survey data with GLMs but does not directly fit mixed models.12 Integrating lme4 or glmmTMB with full survey design considerations (especially for variance estimation reflecting all design features) often requires additional programming, specialized packages (e.g., WeMix, BIFIEsurvey), or careful manual adjustments.
●	Stata: A popular package in econometrics, sociology, and public health research. Stata provides robust capabilities for both GLMMs and complex survey data analysis. The meglm command (part of the me suite for multilevel mixed-effects models) is Stata’s official tool for GLMMs.24 The powerful user-written command gllamm has historically been a go-to for complex GLMMs and can incorporate survey weights.28 Stata's svy prefix can be used with many estimation commands to account for survey design features in variance estimation, though its direct application to all aspects of meglm or gllamm for parameter estimation in GLMMs requires careful attention to documentation.12
●	Other Packages: SUDAAN is a specialized software package designed explicitly for the analysis of data from complex sample surveys and can handle various models.35 SPSS also offers procedures for fitting GLMMs, often relying on methods like Penalized Quasi-Likelihood (PQL) for estimation.25
The choice of software often depends on the analyst's familiarity, the specific complexities of the model and survey design, and the availability of particular estimation algorithms or variance estimation options. The features available within these packages can significantly influence the analytical strategies that are practically feasible.
B. Specific Procedures and Functions
SAS: PROC GLIMMIX
●	The MODEL statement is used to specify the fixed effects, the response variable, its distribution (via DIST=option), and the link function (via LINK=option).
●	The RANDOM statement defines the random effects. The SUBJECT= option specifies the clustering (grouping) variable. The WEIGHT= option within the RANDOM statement can be used to incorporate survey weights at the level of the random effect (e.g., cluster-level weights).29 An OBSWEIGHT= option in the MODEL statement can handle observation-level weights if appropriate for the model.
●	The CLASS statement lists all categorical variables, including those defining clusters or strata if they are part of the model structure.
●	The EMPIRICAL=CLASSICAL (or EMPIRICAL=MBB etc.) option in the PROC GLIMMIX statement requests empirical sandwich estimators for the covariance matrix of fixed effects, which are crucial for robust variance estimation with complex survey data.29
●	The METHOD= option (e.g., METHOD=QUAD for adaptive quadrature, METHOD=LAPLACE for Laplace approximation, METHOD=RSPL for residual pseudo-likelihood) specifies the core estimation algorithm.2
R: lme4::glmer, glmmTMB, survey::svyglm
●	glmer() (from lme4): Uses a formula-based syntax, e.g., response ~ fixed_predictor1 + (1 | cluster_id). A weights argument is available for specifying prior weights. However, using this argument with survey weights requires careful consideration, as glmer itself does not automatically produce standard errors adjusted for the full complex survey design (stratification, multi-stage clustering beyond the specified random effect).
●	glmmTMB: Offers similar formula syntax but supports a wider range of distributions (e.g., beta, negative binomial, zero-inflated models) and more complex covariance structures for random effects. It also has a weights argument. Similar caveats about full survey variance estimation apply as with glmer.
●	survey::svyglm(): This function is designed for fitting GLMs (not mixed models) to complex survey data. It requires a svydesign object (created using svydesign()) which encapsulates all survey design information (PSUs, strata, weights). It produces design-adjusted parameter estimates and standard errors for GLMs.
●	For full GLMM analysis with complex survey features in R, researchers might need to combine outputs from mixed-model packages with variance estimation techniques from the survey package or use specialized packages designed for this purpose.
Stata: meglm, gllamm, svy prefix
●	meglm: This is Stata's official suite for multilevel GLMs. The syntax is typically meglm dependent_var fixed_predictors | | level2_id: random_effect_vars | | level3_id:..., family(distribution) link(linkfunction). It supports various distributions and links and can incorporate weights (pweight, iweight, fweight) at different levels of the hierarchy using options within the random-effects specification.24 Robust standard errors are often default or can be requested.
●	gllamm: A highly flexible user-written program. It uses a more general syntax and can fit a wide array of multilevel and latent variable models, including GLMMs. It has specific options for incorporating survey weights at different levels (e.g., pweight(stubname) where stubname1, stubname2 are weight variables for level 1, level 2, etc.) and can be combined with robust standard error estimation to account for clustering and stratification.28
●	svy: prefix: This prefix can be applied to many Stata estimation commands to adjust for complex survey design features specified via svyset. While svy: can be used with some simpler regression models, its direct and full integration with all features of meglm or gllamm for the parameter estimation of GLMMs (especially regarding how weights are handled within the mixed model likelihood versus in the variance estimation) needs careful verification from Stata's documentation.24 Often, weights are specified directly in the mixed-model command (meglm, gllamm), and robust variance options are used.
The complexity of GLMMs, particularly when combined with the intricacies of survey data, means that users must possess a strong understanding of both the statistical theory and the specific software implementation. Default settings may not always be appropriate, and simply obtaining output without convergence warnings does not guarantee a valid or meaningful model. There is also no perfectly standardized approach across all packages for handling every nuance of survey design features within GLMMs, particularly for multilevel weights and complex stratification schemes. This necessitates that analysts become proficient with the specific capabilities and limitations of their chosen software.
C. Syntax Examples for Specifying Weights, Random Effects, and Survey Design Features
Illustrative syntax snippets demonstrate how these components are typically specified.
SAS (PROC GLIMMIX):

SAS

PROC GLIMMIX DATA=my_survey_data METHOD=QUAD EMPIRICAL=CLASSICAL;
  CLASS cluster_id stratum_id; /* Define clustering and stratification variables */
  MODEL binary_outcome = fixed_cov1 fixed_cov2 / DIST=BINOMIAL LINK=LOGIT SOLUTION;
  RANDOM INTERCEPT / SUBJECT=cluster_id WEIGHT=cluster_weight_var TYPE=VC;
  /* Optional: OBSWEIGHT obs_level_weight_var; if observation-level weights used in MODEL */
  /* Stratification often handled by EMPIRICAL estimator if PSUs are defined within strata */
RUN;

This example fits a logistic GLMM with a random intercept for cluster_id, incorporates cluster_weight_var, and uses an empirical sandwich estimator for robust standard errors.29
R (lme4 style, simplified for weights; full survey SEs need more):

R

library(lme4)
# svy_design_object <- svydesign(...) # From survey package
# my_survey_data <- svy_design_object$variables
# my_survey_data$final_weight <- weights(svy_design_object)

glmm_model <- glmer(binary_outcome ~ fixed_cov1 + fixed_cov2 + (1 | cluster_id),
                    data = my_survey_data,
                    family = binomial(link = "logit"),
                    weights = final_weight)
# Note: Standard errors from summary(glmm_model) are NOT fully survey-adjusted.
# Further steps with packages like 'survey' or 'clubSandwich' would be needed for robust SEs.

This shows basic glmer syntax with a weights argument. Achieving fully survey-robust inference requires more advanced techniques beyond standard lme4 output.23
Stata (meglm):

Stata

svyset psu_id [pweight=final_weight], strata(stratum_id) // Set survey design first (for SEs if compatible)

meglm binary_outcome fixed_cov1 fixed_cov2 |
| cluster_id: , family(binomial) link(logit) vce(robust) pw(obs_level_weight_var_if_different_or_handled_by_meglm)
/*
  Stata's meglm can also handle pweights directly at different levels, e.g.,
  meglm... |
| cluster_id: , pweight(cluster_weight_var)...
  The interaction of svy prefix with meglm for weights in likelihood vs. SEs needs care.
  Often, weights are passed to meglm, and vce(robust) or vce(cluster cluster_id) is used.
*/

This meglm example indicates a logistic model. Stata's meglm commands allow pweight options directly for different levels.24 The vce(robust) or vce(cluster clustvar) option provides robust standard errors.
Stata (gllamm):

Stata

/* Define weight variables: pw1 for level 1, pw2 for level 2 (cluster) */
gllamm binary_outcome fixed_cov1 fixed_cov2, i(cluster_id) link(logit) family(binomial) pweight(pw) adapt robust
/* Assumes pw1 and pw2 exist in the dataset. pw is the stubname. */

This gllamm example demonstrates its ability to handle multilevel weights via the pweight(stubname) option and requests robust standard errors.28
The following table summarizes key software options for fitting GLMMs with complex survey data.
Table 4: Software Options for GLMMs with Complex Survey Data

Feature	SAS (PROC GLIMMIX)	R (lme4::glmer, glmmTMB)	Stata (meglm suite, gllamm)	
Key Procedure(s)/Function(s)	PROC GLIMMIX 34	glmer() (lme4), glmmTMB() (glmmTMB) 8	meglm commands, gllamm (user-written) 24	
Support for Distributions/Links	Extensive (via DIST=, LINK=) 34	Good (family argument in glmer); Very extensive in glmmTMB (incl. zero-inflation, beta) 8	Extensive in meglm (family(), link() options); Very extensive in gllamm 24	
Specification of Random Effects	Flexible via RANDOM statement (SUBJECT=, TYPE=) 34	Formula-based `(random_terms	group_var)inglmer/glmmTMB` 23	
Handling of Survey Weights (Obs & Cluster Level)	WEIGHT= in RANDOM for cluster-level, OBSWEIGHT= in MODEL for obs-level.29	weights argument in glmer/glmmTMB for obs-level. Multilevel/cluster weights often require custom handling or specialized packages (e.g., WeMix).	pweight() option at different levels in meglm and gllamm.24	
Robust Variance Estimation for Survey Design	EMPIRICAL=CLASSICAL (sandwich estimator) accounts for clustering/weighting.29	Requires external packages (e.g., survey, clubSandwich) or manual calculation to combine with glmer/glmmTMB output for full survey design.	vce(robust) or vce(cluster clustvar) in meglm. robust option in gllamm. svy prefix can be used with some commands if compatible.24	
Key Strengths for Survey GLMMs	Well-established; detailed control over estimation; explicit multilevel weight handling in RANDOM.	Open-source; vast array of packages; cutting-edge methods often appear first. glmmTMB very flexible.	Strong survey data capabilities (svy); meglm is official and powerful; gllamm is extremely flexible for complex models and weights.	
Potential Limitations/Considerations	Can be computationally intensive; syntax can be verbose; licensing costs.	Full survey variance estimation with GLMMs not seamlessly integrated into core lme4/glmmTMB; requires more user expertise.	gllamm is user-written (though widely respected). Ensuring full compatibility of svy with all meglm features needs care.	
VI. Applications in Sample Survey Statistics
GLMMs have proven to be invaluable tools across a range of applications in sample survey statistics, primarily due to their ability to handle hierarchically structured data, correlated observations, and non-normal response variables while attempting to account for complex survey designs. Their versatility makes them applicable to diverse fields such as public health, sociology, economics, and education research that rely heavily on survey data.
A. Small Area Estimation (SAE)
Small Area Estimation (SAE) refers to the collection of statistical techniques used to produce reliable estimates of parameters of interest (e.g., means, totals, proportions) for small geographic areas or demographic subgroups for which the direct survey estimates are too unreliable due to small sample sizes.21 GLMMs play a crucial role in SAE by "borrowing strength" across related small areas or from auxiliary information (such as census data or administrative records) to improve the precision of estimates.21
In the SAE context, GLMMs can be formulated as unit-level models (where the model is specified for individual survey respondents) or area-level models (where the model relates direct estimates for areas to area-level covariates). Unit-level GLMMs typically include fixed effects for individual and area-level covariates and random effects for the small areas to capture area-specific deviations not explained by the covariates.21 These random area effects are key to borrowing strength. For instance, in estimating the prevalence of a health condition in small counties, a logistic GLMM might use individual demographic data and county-level characteristics as fixed effects, with random county effects to account for unexplained variation between counties.30
Once the GLMM is fitted using the available survey data, the estimated parameters (fixed effects and variance components) are used to predict the desired characteristic for each small area. Because these are model-based estimates, assessing their uncertainty is critical. Bayesian methods, often implemented via MCMC, are commonly used to obtain credible intervals for small area estimates.30 Frequentist approaches also employ techniques like parametric or non-parametric bootstrapping and Monte Carlo simulations based on the fitted GLMM to construct confidence intervals.21 The need for reliable estimates for small domains, where direct survey data is sparse, has been a significant impetus for the development and application of mixed models, including GLMMs, in survey statistics, as these models provide a formal framework for leveraging broader data patterns and auxiliary information.
B. Analysis of Longitudinal Survey Data
Longitudinal surveys, which collect data from the same sample units (e.g., individuals, households, or firms) at multiple points in time, are essential for studying change, developmental trajectories, and causal processes.7 A key characteristic of longitudinal data is the inherent correlation between repeated measurements on the same unit.
GLMMs are a dominant and highly suitable methodology for analyzing longitudinal survey data.7 In this framework, repeated observations are conceptualized as being nested within individuals (or other units of observation). Random effects are used to model this within-unit correlation and individual-specific heterogeneity in trajectories over time. For example, a random intercept can account for baseline differences between individuals, and a random slope for time can model individual variation in the rate of change.7 This allows for a nuanced understanding of how effects vary across the population.
When longitudinal data are collected via complex survey designs, it is imperative to account for features such as survey weights, clustering (beyond the individual level, e.g., PSUs), and stratification in the GLMM analysis. This often involves the pseudo-likelihood estimation methods and robust variance estimation techniques discussed in previous sections.32
C. Modeling Hierarchical Health or Social Survey Data
Many health, social, and educational surveys possess an intrinsically hierarchical or multilevel structure. Examples include surveys of students nested within classrooms and schools, patients nested within clinics or hospitals, or individuals nested within households and geographic areas (e.g., PSUs or counties).7 GLMMs are naturally suited to analyze such data.
By specifying random effects at different levels of the hierarchy, GLMMs can appropriately account for the correlations among observations within the same higher-level units and estimate the variance attributable to each level. For instance, a study examining factors influencing health service utilization might include individual-level characteristics as fixed effects, with random effects for households and communities to capture unobserved contextual influences and correlations.7 The response variables in such surveys are often non-normal (e.g., binary health status, count of services used, ordinal satisfaction levels), further underscoring the utility of GLMMs. The applications of GLMMs across SAE, longitudinal analysis, and general hierarchical modeling highlight their remarkable versatility. They serve not as a specialized tool for niche problems but as a fundamental methodology for a broad spectrum of complex survey analysis challenges where data exhibit correlation, non-normality, or hierarchical structuring.
D. Illustrative Case Studies
Two case studies exemplify the application of GLMMs to complex survey data.
1. PISA Data Analysis (Reading Proficiency)
The Programme for International Student Assessment (PISA) provides data on student proficiency in various subjects, collected through a complex multistage survey design involving stratification and unequal probabilities of selection, with students nested within schools.28 A common research interest is to model factors associated with student proficiency.
●	GLMM Specification: A typical GLMM for dichotomous reading proficiency (e.g., proficient vs. not proficient) would use a logit link function. Fixed effects could include student-level variables (e.g., socio-economic status, gender) and school-level variables (e.g., school type, resources). Random effects for schools (e.g., random intercepts) are included to account for correlation among students within the same school and to estimate between-school variability.28
●	Handling Survey Features: The analysis described by Rabe-Hesketh and Skrondal (2006) used a pseudo-likelihood estimation approach implemented in the Stata program gllamm. This involved incorporating inverse probability weights to make inferences to the student population. A sandwich estimator was used for standard errors to account for the complex sampling design (stratification and clustering). The study also highlighted the importance of appropriate scaling for level 1 (student-level) weights, noting that regression coefficients for dichotomous responses can be severely biased if weight scaling is not handled carefully. It was conjectured that deriving approximate marginal effects from the GLMM might yield estimates less affected by weight scaling issues.28
2. MMP HIV Provider Survey (Adequacy of Risk Reduction Services)
The Medical Monitoring Project (MMP) HIV Provider Survey utilized a two-stage probability-proportionate-to-size sample design, first sampling states/territories and then HIV facilities, with all providers within selected facilities included. The data included provider-level and facility-level survey weights.10 The research aimed to identify factors associated with providers delivering adequate drug use and sexual risk reduction services.
●	GLMM Specification: Multilevel logistic regression models were fitted for the binary outcomes (adequate vs. inadequate services). Fixed effects included various provider- and facility-level covariates. Random intercepts for facilities were included to capture between-facility variation in service provision.10
●	Handling Survey Features: The models were fitted using SAS PROC GLIMMIX (or equivalently, Stata's svy: melogit). Sampling weights were incorporated at the appropriate levels (provider weights, scaled to sum to facility sample sizes to avoid overstating precision).
●	Key Findings: A crucial finding was that ignoring the sampling weights led to noticeably different inferences for several fixed-effects coefficients and a clear overstatement of the random facility intercept variability. The weighted models also explained more of the variance in the outcomes at each level compared to unweighted models. This case study strongly underscored the necessity of incorporating survey weights in multilevel models when analyzing complex survey data to avoid erroneous conclusions.10
These case studies powerfully demonstrate that the successful application of GLMMs to survey data is not merely a matter of choosing a sophisticated statistical model. It critically depends on the meticulous integration of the model's statistical components (fixed and random effects, distributions, links) with a thorough understanding and incorporation of the survey design features (weights, clusters, strata). Overlooking either aspect can significantly compromise the validity and interpretability of the research findings.
VII. Challenges, Best Practices, and Advanced Topics
While GLMMs offer a powerful framework for analyzing complex survey data, their application is not without challenges. Practitioners must navigate issues related to computational demands, model selection intricacies, interpretation nuances, and specific data features like overdispersion or zero-inflation. Adherence to best practices and awareness of advanced topics are crucial for robust and meaningful analyses.
A. Computational Intensity and Convergence Issues
Fitting GLMMs can be computationally intensive, particularly with large survey datasets, models involving numerous random effects (e.g., crossed random effects or multiple random slopes), or complex random effects covariance structures.1 The numerical integration techniques required to approximate the likelihood (such as adaptive Gauss-Hermite quadrature or MCMC sampling) are inherently demanding. For instance, the use of ML in GLMMs can be limited by the need to evaluate high-dimensional integrals, and MCMC methods carry a substantial computational overhead.1 SAS PROC GLIMMIX, for example, can become very slow with many strata or clusters.32
Furthermore, the optimization algorithms used to find parameter estimates may encounter convergence problems. They might fail to converge to a stable solution, converge very slowly, or converge to a local maximum of the (pseudo-)likelihood surface rather than the desired global maximum.22 Convergence issues can be exacerbated by model misspecification, sparse data in certain cells defined by covariates and random effects, high correlations among parameters, poorly scaled predictor variables, or inappropriate starting values for the optimization algorithm. Problems with survey weight scaling can also contribute to instability.
Several strategies can be employed to address these computational challenges:
●	Model Simplification: If a complex model fails to converge, consider simplifying the random effects structure (e.g., removing random slopes, using simpler covariance structures like variance components instead of unstructured).
●	Starting Values: Providing good starting values for fixed effects and variance components can significantly aid convergence. Some software allows using estimates from a simpler model (e.g., a GLM) as starting values.34
●	Optimization Algorithm Tuning: Software packages often provide options to control the optimization process, such as changing the algorithm, increasing the maximum number of iterations, or adjusting convergence tolerances.34
●	Alternative Estimation Methods: If a highly accurate method like AGHQ with many quadrature points fails or is too slow, one might switch to a faster but potentially less accurate method like Laplace approximation or AGHQ with fewer points.23 For Bayesian models, INLA might be a faster alternative to MCMC.26
●	Data Scaling: Standardizing or centering continuous predictor variables can sometimes improve numerical stability.
●	Software Choice: Different software packages implement algorithms differently, and one package might succeed where another struggles for a particular model or dataset.
Beyond the mathematical formulations, successfully fitting GLMMs, especially with complex and potentially problematic survey data, often involves a degree of practical skill and iterative troubleshooting. This "art" of model fitting—diagnosing convergence issues, selecting appropriate starting values, and judiciously simplifying models when necessary—is as vital as a firm grasp of the underlying statistical theory.
B. Model Selection and Diagnostics in the Survey Context
Choosing the "best" GLMM from a set of candidate models and diagnosing its adequacy are more complex tasks than for simpler models, especially when dealing with survey data.
Model Selection: Standard information criteria like AIC (Akaike Information Criterion) and BIC (Bayesian Information Criterion) are based on the maximized likelihood value. When pseudo-likelihood estimation is used, as is common with survey weights, the theoretical justification for AIC and BIC may not hold directly, or they may require adjustments to account for the survey design and the nature of the pseudo-likelihood. Similarly, for Bayesian models, criteria like DIC (Deviance Information Criterion) can be problematic in random-effects models if not used cautiously.26 Model selection in the survey GLMM context often relies more on a combination of theoretical considerations (which variables should be in the model based on subject-matter knowledge), the significance of parameters (using robust standard errors), changes in variance components, and, where possible, comparisons of predictive performance on hold-out samples (though this is also complicated by survey designs).
Model Diagnostics: Assessing the fit of a GLMM and checking its assumptions is challenging. Standard residuals from GLMs are not directly applicable due to the presence of random effects and non-normal responses. Various types of residuals can be defined for GLMMs (e.g., Pearson, deviance, raw residuals on the linear predictor scale), but their interpretation requires care. Diagnostic plots based on these residuals can help check for patterns, outliers, or misspecification of the link function or distribution. Specialized R packages like DHARMa and HLMdiag provide tools for GLMM diagnostics by simulating data from the fitted model and comparing to observed data.8 When analyzing survey data, diagnostics should ideally also consider the influence of the survey design features, though this is an area where methodology is still developing.
C. Interpretation of Fixed and Random Effects with Survey Weights
The interpretation of parameters from GLMMs fitted with survey data requires careful consideration:
●	Fixed Effects: When survey weights are appropriately incorporated and the model is correctly specified, the fixed-effects coefficients (β) are typically interpreted as estimates of population-average effects of the covariates on the transformed mean of the response (via the link function). However, if the sampling design is informative (i.e., selection probabilities are related to the outcome variable even after conditioning on model covariates), the interpretation can be more nuanced, and simply weighting may not fully recover the true population parameters without explicitly modeling the informativeness.
●	Random Effects (Variance Components): The estimated variances of the random effects quantify the magnitude of unexplained variability between clusters (e.g., PSUs, schools) or subjects at different levels of the hierarchy, after accounting for fixed effects. They indicate the degree of correlation or dependency among observations within the same group.
●	Conditional vs. Marginal Interpretations: GLMM coefficients inherently have a subject-specific or cluster-specific (conditional) interpretation. That is, they describe the effect of a covariate on an individual or cluster with a given set of random effects. For many policy or public health applications, population-averaged (marginal) effects are of greater interest. These marginal effects average over the distribution of random effects. While for linear mixed models with an identity link, conditional and marginal fixed effects are the same, for GLMMs with non-linear link functions, they differ. Marginal effects can be derived from a fitted GLMM, often through numerical integration or simulation, but this requires additional computational steps.4 It has been suggested that these derived marginal effects might be less sensitive to issues like the scaling of level 1 survey weights in certain models.28
D. Handling Overdispersion and Zero-Inflation
Survey data, particularly counts or proportions, often exhibit features that violate the assumptions of simpler distributions.
●	Overdispersion: This occurs when the observed variance in the data is greater than what is predicted by the nominal distribution (e.g., variance > mean for Poisson, or more variance than np(1−p) for binomial).2 Overdispersion can arise from unobserved heterogeneity or additional clustering not captured by the model. Ignoring it leads to underestimated standard errors and potentially spurious findings. GLMMs can address overdispersion in several ways:
○	By including an observation-level random effect (an extra random term for each observation).
○	By choosing a distribution that inherently allows for overdispersion, such as the Negative Binomial distribution for count data (instead of Poisson) or the Beta-Binomial distribution for proportions (instead of Binomial).2 Many software packages for GLMMs support these distributions.
●	Zero-Inflation: Count data in surveys frequently exhibit an excess of zero values compared to what would be expected under a standard Poisson or Negative Binomial distribution (e.g., number of cigarettes smoked, number of chronic conditions).12 This can occur if there are two underlying processes: one generating zeros (e.g., structural non-participation) and another generating counts for those who do participate. Standard GLMMs may not fit such data well. Specialized zero-inflated models, such as Zero-Inflated Poisson (ZIP) or Zero-Inflated Negative Binomial (ZINB) GLMMs, are designed for this situation. These models typically combine a logistic (or probit) model for the probability of an excess zero with a count model for the non-excess zeros.8
E. Limitations of Current Methods and Areas for Future Research
Despite their power, current methods for GLMM estimation with survey data have limitations and point to areas needing further research:
●	Bias in Pseudo-Likelihood Estimators: PL methods, while pragmatic for incorporating survey weights, are based on approximations to the true likelihood. These approximations can introduce bias into parameter estimates, particularly for variance components or with certain data types like sparse binary data or when clusters are small.22 The properties of PL estimators under various complex survey designs are not always fully understood.
●	Computational Scalability: Applying GLMMs with sophisticated random effects structures or advanced estimation techniques (like AGHQ with many points or long MCMC runs) to very large-scale national surveys with tens or hundreds of thousands of records remains a significant computational hurdle.30
●	Fully Bayesian GLMMs for Complex Surveys: Developing computationally efficient and theoretically sound methods for fully Bayesian GLMMs that can rigorously incorporate all features of complex survey designs (especially multilevel weights, intricate stratification, and PSUs from different selection stages) is an ongoing challenge.26 The interplay between prior specification and survey weighting needs more exploration.
●	Model Diagnostics for Survey GLMMs: There is a need for more comprehensive and accessible diagnostic tools specifically tailored for GLMMs fitted to complex survey data, which can help assess model adequacy while accounting for the design.
●	Software Integration and Usability: Continued development of statistical software that seamlessly integrates robust GLMM fitting procedures with comprehensive handling of complex survey design features (including user-friendly specification of multilevel weights and appropriate variance estimation) is crucial for wider and more accurate application.
The application of GLMMs involves a near-constant reliance on approximations—be it numerical integration for the likelihood, methodological adjustments in pseudo-likelihood, or sampling-based approximations in MCMC. This implies that reported results (p-values, confidence intervals, parameter estimates) carry not only sampling error but also error due to these approximations. This is a fundamental consideration for practitioners, especially when interpreting results from highly complex models or when findings are marginal. Furthermore, a persistent tension exists between the desire to specify highly complex random effect structures in GLMMs to mirror intricate data dependencies and the challenge of simultaneously and fully accounting for all complexities of a survey design. Often, pragmatic compromises are made, such as using robust standard errors to cover some design aspects rather than explicitly modeling every feature within the likelihood framework itself.
VIII. Concluding Remarks and Recommendations
Generalized Linear Mixed Models represent a significant advancement in the statistical analysis of sample survey data, offering the flexibility to model diverse response types and the power to account for the hierarchical and correlated nature inherent in many survey designs. Their proper application, however, demands a sophisticated understanding of both the model's theoretical underpinnings and the specific complexities introduced by the survey sampling methodology.
A. Summary of Key Estimation Principles for GLMMs in Survey Statistics
The analysis of survey data using GLMMs is guided by several core principles. Firstly, the choice of a GLMM is motivated by the need to handle response variables that may not be normally distributed and, crucially, to address the correlations among observations that arise from clustered sampling designs or longitudinal data collection. Random effects are the primary mechanism within GLMMs for modeling these dependencies and capturing unobserved heterogeneity between sampling units or over time.
Secondly, it is paramount that the complex features of the survey design—particularly survey weights, clustering, and stratification—are appropriately accounted for in the estimation process. Failure to do so can lead to biased parameter estimates and invalid standard errors, thereby undermining statistical inference. Survey weights are typically incorporated using pseudo-likelihood estimation methods, which modify the likelihood function to reflect the sample's representativeness of the target population. Clustering is naturally handled by the random effects structure corresponding to the survey clusters (e.g., PSUs). Stratification is often addressed through robust variance estimation techniques or by including strata indicators as fixed effects if feasible.
Thirdly, due to the analytical intractability of the GLMM likelihood function, parameter estimation relies on numerical approximation techniques. These include likelihood-based methods employing numerical integration (such as Adaptive Gauss-Hermite Quadrature or Laplace approximation), pseudo-likelihood approaches (often involving iterative fitting of weighted linear mixed models), and Bayesian methods (typically using MCMC or INLA). Each approach has its own set of advantages, limitations, and computational demands, particularly when adapted for survey data.
Finally, robust variance estimation (e.g., sandwich estimators, Taylor series linearization, or replication methods like bootstrapping) is essential when working with GLMMs in the survey context. These methods provide standard errors for parameter estimates that more accurately reflect the variability due to the complex sampling design and the use of pseudo-likelihood or other approximate estimation techniques.
Successfully navigating the application of GLMMs in survey statistics necessitates a holistic approach. This involves an integrated understanding of statistical theory pertaining to GLMMs, deep knowledge of survey methodology and sampling design principles, and proficiency in computational practice using appropriate statistical software. These elements are not isolated domains but are intrinsically linked facets of a comprehensive analytical endeavor.
B. Recommendations for Practitioners
Given the complexities involved, practitioners aiming to use GLMMs for analyzing complex survey data should consider the following recommendations:
1.	Thoroughly Understand the Data and Survey Design: Before model fitting, gain a comprehensive understanding of the survey methodology, including how clusters and strata were defined, how weights were calculated (and if multilevel weights exist), and the nature of the variables (distribution, potential for overdispersion or zero-inflation).
2.	Select Appropriate Software and Estimation Methods Carefully: Choose a statistical package that offers the necessary flexibility for specifying the desired GLMM (distributions, links, random effects structures) and can adequately handle the survey design features, particularly weights and robust variance estimation. Be aware of the specific algorithms implemented (e.g., PQL, Laplace, AGHQ) and their known properties and limitations.
3.	Start with Simpler Models: Begin with a more parsimonious model and gradually increase complexity (e.g., by adding random slopes or more complex covariance structures for random effects). At each step, carefully check for model convergence, the plausibility of parameter estimates (especially variance components), and changes in model fit (where appropriate design-adjusted criteria can be used).
4.	Pay Meticulous Attention to Survey Weights: Critically evaluate how survey weights are to be incorporated. If multilevel weights are present, ensure the software and model can accommodate them at the correct levels. Consider the potential need for weight scaling, understand its implications (especially for variance components), and document the approach taken.
5.	Prioritize Robust Variance Estimation: Always employ variance estimation methods that are appropriate for the complex survey design and the GLMM estimation technique used. This typically involves using sandwich estimators or replication methods that account for weighting, clustering, and stratification.
6.	Conduct Thorough Model Diagnostics: Despite the challenges, make every effort to assess model fit and check assumptions. Examine residuals (appropriately defined for GLMMs), assess the impact of influential observations or clusters, and check the plausibility of the random effects distribution.
7.	Acknowledge Model Limitations and Approximation: Be transparent about the estimation methods used and any approximations involved. Interpret results cautiously, especially if the model is very complex, data are sparse in some areas, or convergence was difficult to achieve.
8.	Consult with Experts: Given the advanced nature of these methods, seeking collaboration or consultation with statisticians who have specific expertise in fitting GLMMs to complex survey data can be invaluable.
9.	Stay Informed of Methodological Developments: The field of GLMMs for survey data is an active area of research. New estimation techniques, diagnostic tools, and software capabilities continue to emerge. Staying abreast of these developments is important for applying best practices.
While the ideal is a perfectly specified model that flawlessly accounts for every nuance of the data and survey design, practitioners often operate within constraints of data availability, computational resources, and methodological limitations. Therefore, a degree of pragmatism is often required, balancing theoretical rigor with practical feasibility to achieve the most reliable and meaningful inferences possible. The careful and informed application of GLMMs, guided by these principles and recommendations, can yield powerful insights from complex sample survey data.
Works cited
1.	ecommons.cornell.edu, accessed May 15, 2025, https://ecommons.cornell.edu/bitstream/handle/1813/31937/BU-1340-MA.pdf;sequence=1
2.	griffin.uga.edu, accessed May 15, 2025, https://griffin.uga.edu/files/2023/05/CountsGLMMfin1.pdf
3.	1. Introduction Generalized linear mixed models (GLMMs) are a powerful class of - McMaster Math, accessed May 15, 2025, https://ms.mcmaster.ca/~bolker/classes/s4c03/notes/GLMM_Bolker_draft5.pdf
4.	The value of generalized linear mixed models for data analysis in the plant sciences, accessed May 15, 2025, https://www.frontiersin.org/journals/horticulture/articles/10.3389/fhort.2024.1423462/full
5.	Generalized linear model - Wikipedia, accessed May 15, 2025, https://en.wikipedia.org/wiki/Generalized_linear_model
6.	Generalized Linear Models Link Function The logistic equation is stated in terms of the probability that Y = 1, which is π, and, accessed May 15, 2025, https://web.pdx.edu/~newsomj/cdaclass/ho_glm.pdf
7.	Generalized Linear Mixed Model in SPSS - GLMM Explained, accessed May 15, 2025, https://spssanalysis.com/generalized-linear-mixed-model-in-spss/
8.	CRAN Task View: Mixed, Multilevel, and Hierarchical Models in R - R Project, accessed May 15, 2025, https://cran.r-project.org/view=MixedModels
9.	support.sas.com, accessed May 15, 2025, https://support.sas.com/resources/papers/proceedings15/2600-2015.pdf
10.	Weighted Multilevel Models: A Case Study - PMC, accessed May 15, 2025, https://pmc.ncbi.nlm.nih.gov/articles/PMC4605185/
11.	Accounting for Clustering in the Analysis of MEPS Data: Frequently Asked Questions (FAQ), accessed May 15, 2025, https://meps.ahrq.gov/survey_comp/hc_clustering_faq.pdf
12.	%SURVEYGENMOD Macro: An Alternative to Deal with Complex Survey Design for the GENMOD Procedure - SAS Support, accessed May 15, 2025, https://support.sas.com/resources/papers/proceedings17/0268-2017.pdf
13.	Multilevel models vs GLMMs for correlated clustered data - Cross Validated, accessed May 15, 2025, https://stats.stackexchange.com/questions/184012/multilevel-models-vs-glmms-for-correlated-clustered-data
14.	Multilevel Modeling Approach for Hierarchical Data an Empirical Investigation, accessed May 15, 2025, https://www.scirp.org/journal/paperinformation?paperid=137832
15.	Exploring Generalized Estimating Equations (GEE) in Applied Models - Number Analytics, accessed May 15, 2025, https://www.numberanalytics.com/blog/exploring-generalized-estimating-equations-gee-applied-models
16.	Report Quality of Generalized Linear Mixed Models in Psychology: A Systematic Review, accessed May 15, 2025, https://www.frontiersin.org/journals/psychology/articles/10.3389/fpsyg.2021.666182/full
17.	documentation.sas.com, accessed May 15, 2025, https://documentation.sas.com/doc/en/statug/15.1/statug_intromod_sect010.htm#:~:text=Models%20in%20which%20all%20effects,effects%20and%20random%2Deffects%20terms.
18.	Generalized linear mixed models: a practical guide for ecology and evolution - UCF College of Sciences, accessed May 15, 2025, https://sciences.ucf.edu/biology/d4lab/wp-content/uploads/sites/23/2020/02/Bolker-et-al-2009-TREE.pdf
19.	When to use fixed effects vs random effects vs mixed-effects models? : r/AskEconomics, accessed May 15, 2025, https://www.reddit.com/r/AskEconomics/comments/qkj1gc/when_to_use_fixed_effects_vs_random_effects_vs/
20.	Generalized Linear Mixed Models: A Pseudo-Likelihood Approach - ResearchGate, accessed May 15, 2025, https://www.researchgate.net/publication/265367587_Generalized_Linear_Mixed_Models_A_Pseudo-Likelihood_Approach
21.	cimat.repositorioinstitucional.mx, accessed May 15, 2025, https://cimat.repositorioinstitucional.mx/jspui/bitstream/1008/351/2/TE%20517.pdf
22.	GLMM FAQ - GitHub Pages, accessed May 15, 2025, https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html
23.	glmer: Fitting Generalized Linear Mixed-Effects Models in lme4 - rdrr.io, accessed May 15, 2025, https://rdrr.io/cran/lme4/man/glmer.html
24.	mecloglog — Multilevel mixed-effects complementary log-log regression - Title Description Quick start Menu, accessed May 15, 2025, https://www.stata.com/manuals14/memecloglog.pdf
25.	Does GLMM use penalized quasi-likelihood for non-normal data? | SPSS Statistics, accessed May 15, 2025, https://community.ibm.com/community/user/discussion/does-glmm-use-penalized-quasi-likelihood-for-non-normal-data
26.	Bayesian inference for generalized linear mixed models - PMC, accessed May 15, 2025, https://pmc.ncbi.nlm.nih.gov/articles/PMC2883299/
27.	Bayesian pseudo-empirical-likelihood intervals for complex surveys - SciSpace, accessed May 15, 2025, https://scispace.com/pdf/bayesian-pseudo-empirical-likelihood-intervals-for-complex-5ejiqlrhv1.pdf
28.	www.gllamm.org, accessed May 15, 2025, http://www.gllamm.org/JRSSAsurvey_06.pdf
29.	SAS Help Center: Pseudo-likelihood Estimation for Weighted Multilevel Models, accessed May 15, 2025, https://documentation.sas.com/doc/en/amocdc/8.2/pgmsascdc/statug/statug_glimmix_details25.htm
30.	Constructing Statistical Intervals for Small Area Estimates Based on ..., accessed May 15, 2025, https://pmc.ncbi.nlm.nih.gov/articles/PMC9336217/
31.	Matching with Sampling Weights - CRAN, accessed May 15, 2025, https://cran.r-project.org/web/packages/MatchIt/vignettes/sampling-weights.html
32.	support.sas.com, accessed May 15, 2025, https://support.sas.com/resources/papers/proceedings20/4937-2020.pdf
33.	FAQs: Specifying survey weights in gllamm, accessed May 15, 2025, http://www.gllamm.org/faqs/estimation/pfound.html
34.	PROC GLIMMIX Statement - SAS Help Center, accessed May 15, 2025, https://documentation.sas.com/doc/ja/pgmsascdc/v_041/statug/statug_glimmix_syntax01.htm
35.	Choosing Statistical Software - Medical College of Wisconsin, accessed May 15, 2025, https://www.mcw.edu/-/media/MCW/Departments/Biostatistics/choosingstatisticalsoftware61512.pdf?la=en&hash=343ED91E2B4D7A07E371640FE408A7C1
36.	Incorporating survey weights with weightthem in R - Stack Overflow, accessed May 15, 2025, https://stackoverflow.com/questions/78211930/incorporating-survey-weights-with-weightthem-in-r
37.	Modern methods for longitudinal data analysis, capabilities, caveats and cautions - PMC, accessed May 15, 2025, https://pmc.ncbi.nlm.nih.gov/articles/PMC5434286/
38.	Generalized linear mixed-effects models for studies using different sets of stimuli across conditions - Frontiers, accessed May 15, 2025, https://www.frontiersin.org/journals/psychology/articles/10.3389/fpsyg.2022.955722/full
39.	Pseudo Maximum Likelihood Approach for the Analysis of Multivariate Left-Censored Longitudinal Data - PMC, accessed May 15, 2025, https://pmc.ncbi.nlm.nih.gov/articles/PMC5138145/