Chapter 6 Generalized linear models
When the variable of interest is not continuous or the assumptions of the classical linear model are not satisfied, generalized linear models provide an appropriate alternative for modeling the relationships of interest. This approach was introduced by Nelder & Wedderburn (1972), who showed that several widely used statistical methods, such as logistic regression, log-linear models for count data, and certain advanced regression models for continuous variables, could be integrated within a single unified theoretical framework. This formulation made it possible to extend regression tools beyond the settings in which the response variable follows a normal distribution.
The need for this generalization arises because classical linear models assume, among other things, normality and homogeneity of variance in the response variable, conditions that are rarely met when categorical variables, proportions, or counts are analyzed. In these situations, generalized linear models offer a flexible framework for appropriately modeling the relationship between variables by selecting a probability distribution suitable for the response and a link function that connects the mean of that distribution with the linear predictor.
To illustrate the models presented in this chapter, the survey_data database and the BigCity population are used, following the same structure used in previous chapters. After loading tidyverse (Wickham, Averick, et al., 2026), survey (Lumley, 2024), srvyr (Freedman Ellis & Schneider, 2024), broom (Robinson et al., 2026), and jtools (Long, 2026), the following code block creates two new categorical variables (poverty and unemployment) from the original variables Poverty and Employment, respectively.
options(digits = 4)
library(tidyverse)
library(survey)
library(srvyr)
library(broom)
library(jtools)
survey_data <- readRDS("Data/encuesta.rds")
data("BigCity", package = "TeachingSampling")
survey_design <- survey_data %>%
as_survey_design(
strata = Stratum,
ids = PSU,
weights = wk,
nest = TRUE
)
survey_design <- survey_design %>%
mutate(
poverty = factor(
ifelse(Poverty != "NotPoor", 1, 0),
levels = c(0, 1),
labels = c("=Poor", "=Not poor")
),
unemployment = factor(ifelse(Employment == "Unemployed", 1, 0),
levels = c(0, 1),
labels = c("=Unemployed", "=Not unemployed")
)
)6.1 Logistic regression model
When the variable of interest is dichotomous, taking values of zero or one, the classical linear model is not appropriate because it can generate predictions outside the interval \((0,1)\) and because the variance of the response depends on its mean. To address these limitations, the logistic regression model is used; it is one of the most widely used generalized linear models for analyzing binary variables. In the context of household surveys, logistic regression can be used to analyze the association between a condition of interest, such as poverty, and sociodemographic characteristics of individuals or households.
Consider a dichotomous variable \(Y_k \in \{0,1\}\), whose expectation is \(\theta_k=Pr(Y_k=1)\). The logistic regression model relates this probability to a set of explanatory variables through the logit link function, defined as
\[ \log\left(\frac{\theta_k}{1-\theta_k}\right) = \beta_0+\beta_1x_{k1}+\cdots+\beta_px_{kp} =\mathbf{x}_{k}\boldsymbol{\beta} \]
Where \(\boldsymbol{\beta}\) is the vector of regression coefficients associated with the covariates. This logit transformation maps the restricted probability interval \((0,1)\) onto the entire real line, allowing the relationship between the response and the auxiliary variables to be modeled through a linear predictor. Equivalently, the probability of success can be expressed by the following relationship:
\[ \theta_k = Pr(Y_k = 1 \mid \mathbf{x}_k) = \frac{\exp(\mathbf{x}_k\boldsymbol{\beta})} {1+\exp(\mathbf{x}_k\boldsymbol{\beta})} \]
Equivalently, the probability that the unit does not have the characteristic of interest is \(Pr(Y_k = 0 \mid \mathbf{x}_k) = 1-\theta_k\). Under this parameterization, each coefficient in \(\boldsymbol{\beta}\) represents the effect of a covariate on the logarithm of the odds that the event of interest occurs, holding the other variables in the model constant. The parameters can be estimated using the pseudo-maximum likelihood approach. This method incorporates the expansion factors into the conventional likelihood function in order to obtain estimators that are consistent with respect to the target population. Following Heeringa et al. (2017), the pseudo-likelihood function for the logistic model is defined as
\[ PL(\boldsymbol{\beta}) = \prod_{k} \left[ \theta_{k}^{y_{k}} \left\{ 1-\theta_{k} \right\}^{1-y_{k}} \right]^{w_{k}} \]
Where the subscript \(k\) represents an observed unit in the sample. Likewise, \(\theta_{k}\) denotes the probability that unit \(k\) in PSU \(i\) in stratum \(h\) has the characteristic of interest, \(y_{k}\) corresponds to the observed value of the dichotomous variable for that unit, and \(w_{k}\) represents the expansion factor associated with the observation.
The contribution of each observation to the estimation process is determined by its respective expansion factor, allowing the function to appropriately reflect the characteristics of the complex survey design. Because maximizing the pseudo-likelihood is equivalent to maximizing its logarithm, the problem is therefore reduced to optimizing the following function:
\[ \ell(\boldsymbol{\beta}) = \sum_{k} w_{k} \left[ y_{k}\log(\theta_{k}) + (1-y_{k}) \log(1-\theta_{k}) \right] \]
This expression is the objective function used by numerical optimization algorithms to obtain the estimator of the model parameters. It is important to note that the preceding function is not a likelihood in the strict sense under the complex survey design. For this reason, the term pseudo-likelihood is used, and the estimators obtained by maximizing it are known as pseudo-maximum likelihood estimators (Pseudo Maximum Likelihood Estimators, PMLE).
According to Molina & Skinner (1992), under regularity conditions, these estimators are consistent and asymptotically normal, forming the basis of the inference procedures implemented in most statistical packages for the analysis of complex survey data.
Once the pseudo-maximum likelihood estimators have been obtained, their variances are estimated using Taylor linearization techniques (Binder, 1983), which explicitly incorporate the characteristics of the complex sampling design. The variance-covariance matrix of the estimated parameters, denoted by \(\text{Var}(\hat{\boldsymbol{\beta}})=\boldsymbol{\Sigma}\), is obtained from a linear approximation of the estimating equations around the true value of the parameters.
This matrix describes the precision of the estimators, since its diagonal elements correspond to the variances of each coefficient, while the off-diagonal elements represent the covariances between pairs of coefficients. Its estimation is the basis for calculating standard errors, constructing confidence intervals, and carrying out hypothesis tests on the model parameters.
Because the data come from a complex sample design, it is necessary to appropriately incorporate the weights, stratification, and clustering in order to obtain consistent estimators of the parameters and their standard errors. In R, this task can be carried out using the svyglm() function from the survey package, which extends generalized linear models to the context of complex surveys.
Continuing with the example survey, and in order to analyze the factors associated with poverty status, a logistic regression model is fitted using expenditure, employment status, and sex as explanatory variables. Because the data come from a survey with a complex sample design, estimation is carried out using the svyglm() function from the survey package, specifying family = binomial to indicate that the response variable is dichotomous and that the logit link function will be used. The results of the fit are presented in Table 6.1.
logit_model <- svyglm(
formula = poverty ~ Expenditure + Employment + Sex,
family = binomial,
design = survey_design
)In the model specification, the argument formula = poverty ~ Expenditure + Employment + Sex defines the relationship between the response variable poverty and the covariates included in the linear predictor. For categorical variables, the function automatically generates the corresponding indicator variables and estimates their effects relative to a reference category. The argument design = survey_design, in turn, incorporates the sample design information. The resulting object, logit_model, contains the estimated parameters, their precision measures, and the statistics needed to make inferences about the relationship between poverty status and the explanatory variables considered in the model.
The results in Table 6.1 show the estimates of the regression coefficients, together with their standard errors, confidence intervals, and p-values. The model indicates that a higher level of expenditure is associated with a lower probability of being in poverty. Likewise, compared with the reference category (unemployed), both inactive and employed persons have a lower propensity to be poor, with this effect being more pronounced among the latter. In contrast, sex has a very small influence on the probability of poverty once the effect of the other covariates included in the model is controlled for.
## 2.5 % 97.5 %
## (Intercept) 1.448618 2.926867
## Expenditure -0.007044 -0.003752
## EmploymentInactive -1.618341 -0.108547
## EmploymentEmployed -2.076705 -0.737132
## SexMale -0.309026 0.314060| Parameter | term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|---|
| 1 | (Intercept) | 2.188 | 0.373 | 5.863 | 0.000 |
| 2 | Expenditure | -0.005 | 0.001 | -6.497 | 0.000 |
| 3 | EmploymentInactive | -0.863 | 0.381 | -2.266 | 0.025 |
| 4 | EmploymentEmployed | -1.407 | 0.338 | -4.161 | 0.000 |
| 5 | SexMale | 0.003 | 0.157 | 0.016 | 0.987 |
Figure 6.1 presents the distribution of the estimators of the regression coefficients. It can be seen that, for all covariates except sex, the confidence intervals do not include the value zero, providing evidence of a statistically significant association. In contrast, the confidence interval associated with the sex variable contains the value zero, so there is not enough evidence to conclude that this variable has a significant effect once the other covariates included in the model are controlled for.
To visualize the distribution of the coefficients, plot_summs() from the jtools package is used, together with ggstance (Henry et al., 2024) for the horizontal geometry of the intervals.
Figure 6.1: Distribution of the parameters of the baseline logistic model
The model can be extended by incorporating interaction terms to assess whether the effect of one explanatory variable depends on the values taken by another. In this case, an interaction between sex and employment status is included through the term Sex:Employment, with the aim of analyzing whether the relationship between employment status and the probability of poverty is the same for men and women. Without this term, the model assumes that the effect of employment status on poverty is identical for both sexes. By incorporating the interaction, however, this effect is allowed to vary between men and women, capturing possible differences in the way labor-market status is associated with poverty in each group. The results of the fitted model are presented in Table 6.2.
interaction_logit_model <- svyglm(
formula = poverty ~ Expenditure + Employment + Sex + Sex:Employment,
family = binomial,
design = survey_design
)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 1.770 | 0.610 | 2.900 | 0.004 |
| Expenditure | -0.005 | 0.001 | -6.473 | 0.000 |
| EmploymentInactive | -0.395 | 0.594 | -0.665 | 0.507 |
| EmploymentEmployed | -1.059 | 0.570 | -1.860 | 0.066 |
| SexMale | 0.587 | 0.748 | 0.785 | 0.434 |
| EmploymentInactive:SexMale | -0.846 | 0.860 | -0.984 | 0.327 |
| EmploymentEmployed:SexMale | -0.481 | 0.760 | -0.633 | 0.528 |
When the interaction between sex and employment status is incorporated, expenditure continues to show a negative and statistically significant association with the probability of poverty, remaining the covariate with the strongest explanatory evidence in the model. In contrast, the main effects of employment status lose statistical importance once the interaction terms are included, suggesting greater uncertainty in the estimation of their effects. Neither the main effect of sex nor the coefficients associated with the interactions between sex and employment status are statistically significant at the 5% level. Consequently, the inclusion of the interaction terms does not appear to provide additional relevant information compared with the model without interaction.
Figure 6.2 presents a comparison of the coefficient estimates and their corresponding confidence intervals for the models with and without interaction. It can be seen that incorporating the interaction terms increases the uncertainty associated with several of the parameters, as reflected in wider confidence intervals. Likewise, the interaction coefficients have estimates close to zero and wide confidence intervals that include that value, which is consistent with the lack of statistical evidence for concluding that the effect of employment status on the probability of poverty differs between men and women.
Figure 6.2: Comparison of the parameters of the logistic model with and without interactions
6.2 Multinomial regression model
In household surveys, it is common to find response variables with more than two categories, such as employment status, whose possible states include employed, unemployed, and inactive. When these categories are mutually exclusive and do not have a natural order, an appropriate tool for modeling their relationship with a set of covariates is multinomial logistic regression, which is an extension of the binary logistic regression model.
This model makes it possible to estimate the probability of belonging to each category of the response variable as a function of continuous or categorical explanatory variables. Its application requires the response categories to be exhaustive and mutually exclusive and requires that there be no severe multicollinearity among the covariates. In addition, for continuous covariates, a linear relationship is assumed between them and the logarithm of the odds ratios (log-odds) of each category relative to a reference category.
The multinomial logistic model specifies the probability \(\theta_{k(j)}\) that unit \(k\) belongs to category \(j\), with \(j=1,\ldots,J\), given the covariate vector \(\mathbf{x}_k\), through the following expression:
\[ Pr(Y_{k} = j \mid \mathbf{x}_{k}) = \theta_{k(j)} = \frac{\exp(\mathbf{x}_{k}\boldsymbol{\beta}_j)} {\sum_{l=1}^{J} \exp(\mathbf{x}_{k}\boldsymbol{\beta}_l)} \]
Where \(\boldsymbol{\beta}_j\) is the vector of regression coefficients for the \(j\)-th category. This formulation directly expresses the probabilities of belonging to each category; in practice, it is more convenient to rewrite the model using a reference category. Alternative parameterizations facilitate the interpretation of the parameters, since each set of coefficients describes the effect of the covariates on the logarithm of the odds ratios (log-odds) of belonging to a given category compared with the reference category.
The model parameters are estimated using the pseudo-maximum likelihood approach. Following Heeringa et al. (2017), the multinomial pseudo-likelihood function is defined as
\[ PL(\boldsymbol{\beta} \mid \mathbf{X}) = \prod_{k} \left\{ \prod_{j=1}^{J} \theta_{k(j)}^{y_{k(j)}} \right\}^{w_{k}} \]
Where \(k\) represents the observation units belonging to the complex sample. Likewise, \(y_{k(j)}\) is an indicator variable that takes the value one when that unit belongs to category \(j\) and zero otherwise, and \(w_{k}\) corresponds to the expansion factor associated with the observed unit. For mathematical and computational convenience, the parameters are estimated from the pseudo-log-likelihood, whose maximization leads to the same estimators as the original pseudo-likelihood.
\[ \ell(\boldsymbol{\beta}) = \sum_{k} w_k \sum_{j=1}^{J} y_{k(j)} \log(\theta_{k(j)}) \]
As mentioned above, the estimators obtained through this procedure are consistent and asymptotically normal under regular conditions (Molina & Skinner, 1992). Statistical inference for the parameters of interest is based on variance estimates obtained through Taylor linearization techniques that explicitly incorporate the characteristics of the complex sample design, making it possible to calculate standard errors, confidence intervals, and hypothesis tests appropriate for survey data (Binder, 1983).
For fitting the model using R, as an initial reference, the proportions of persons by employment status are estimated, restricting the units of interest to those older than 15 years, as presented in Table 6.3.
survey_design %>%
filter(Age >= 15) %>%
group_by(Employment) %>%
summarise(proportion = survey_mean(vartype = c("se", "ci")))| Employment | proportion | proportion_se | proportion_low | proportion_upp |
|---|---|---|---|---|
| Unemployed | 0.043 | 0.007 | 0.029 | 0.057 |
| Inactive | 0.384 | 0.015 | 0.354 | 0.414 |
| Employed | 0.573 | 0.014 | 0.545 | 0.601 |
To model a multinomial response variable, the svy_vglm() function from the svyVGAM package (Lumley, 2026) is used; it allows multinomial regression models to be fitted while explicitly incorporating the characteristics of the complex survey design. In this example, a model is specified in which activity status is explained as a function of age, sex, and area of residence. Unlike conventional multinomial models, svy_vglm() calculates the variances of the estimators while accounting for the complexity of the sample design, including stratification, clustering, and expansion weights. In this way, both the point estimates and their standard errors and associated inferences are consistent with the survey design.
library(svyVGAM)
survey_design_15 <- survey_design %>%
filter(Age >= 15)
multinomial_model <- svy_vglm(
formula = Employment ~ Age + Sex + Zone,
design = survey_design_15,
crit = "coef",
family = multinomial(refLevel = "Unemployed")
)The argument family = multinomial(refLevel = "Unemployed") defines a multinomial logistic model using the Unemployed category as the reference. Under this specification, the logarithms of the odds ratios (log-odds) between each of the Employment categories and the reference category are estimated simultaneously. The parameters are obtained through pseudo-maximum likelihood, incorporating the survey expansion factors into the objective function.
Because broom::tidy() cannot be applied directly to objects of class svyVGAM, the following helper function is defined to structure the model results in a standardized tabular format5:
tidy.svyVGAM <- function(x, conf.int = FALSE, conf.level = 0.95,
exponentiate = FALSE, ...) {
ret <- as_tibble(summary(x)$coeftable, rownames = "term")
colnames(ret) <- c("term", "estimate", "std.error", "statistic", "p.value")
coefs <- tibble::enframe(stats::coef(x), name = "term", value = "estimate")
ret <- left_join(coefs, ret, by = c("term", "estimate"))
if (conf.int) {
ci <- broom:::broom_confint_terms(x, level = conf.level, ...)
ret <- dplyr::left_join(ret, ci, by = "term")
}
if (exponentiate) ret <- broom:::exponentiate(ret)
ret %>%
tidyr::separate(term, into = c("term", "y.level"), sep = ":") %>%
arrange(y.level) %>%
relocate(y.level, .before = term)
}Table 6.4 presents the estimated coefficients of the multinomial logistic model, using the Unemployed category as the reference. The results show that age has a positive and statistically significant effect on the relative probabilities of belonging to categories 1 and 2 compared with unemployment. Likewise, the variable SexMale has negative and significant coefficients, indicating that men have lower relative probabilities than women of being in those categories. By contrast, the area-of-residence variable does not show statistically significant effects at the conventional 5% level.
| y.level | term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|---|
| 1 | (Intercept) | 2.556 | 0.550 | 4.648 | 0.000 |
| 1 | Age | 0.022 | 0.010 | 2.117 | 0.034 |
| 1 | SexMale | -2.185 | 0.317 | -6.901 | 0.000 |
| 1 | ZoneUrban | -0.254 | 0.404 | -0.627 | 0.531 |
| 2 | (Intercept) | 2.306 | 0.461 | 4.998 | 0.000 |
| 2 | Age | 0.018 | 0.009 | 2.056 | 0.040 |
| 2 | SexMale | -0.578 | 0.277 | -2.091 | 0.037 |
| 2 | ZoneUrban | 0.039 | 0.361 | 0.109 | 0.913 |
Figure 6.3 shows the confidence intervals of the estimated coefficients for each equation of the multinomial model. The visualization is built with dotwhisker::dwplot() (Solt & Hu, 2025), which makes it possible to compare coefficients and intervals on a common scale, and gtools::stars.pval() (Warnes et al., 2023) is used to encode statistical significance. Covariates whose intervals do not include the value zero show evidence of a statistically significant association with the response variable, while those whose intervals contain zero do not show significant effects at the confidence level considered.
multinomial_results %>%
mutate(
model = if_else(y.level == 1, "Inactive", "Employed"),
sig = gtools::stars.pval(p.value)
) %>%
dotwhisker::dwplot(
dodge_size = 0.3,
vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)
) +
guides(color = guide_legend(reverse = TRUE)) +
theme_bw() +
theme(legend.position = "top")
Figure 6.3: Confidence intervals for the coefficients of the multinomial model
6.3 Gamma regression model
Gamma regression is an extension of generalized linear models suitable for continuous, strictly positive response variables whose variability increases with the mean. The Gamma model is an appropriate alternative for analyzing strictly positive continuous variables whose variability increases as their expected value grows. This behavior is common in economic and social variables, such as household income or consumption expenditure.
If \(Y_{k}\), the variable of interest observed for unit \(k\), follows a Gamma distribution with mean \(\mu_{k}\) and dispersion parameter \(\phi\), the model relates the expected value of the response variable to a set of covariates through a link function. The most commonly used link function is the logarithmic link, which ensures that predictions are always positive. In this case, the model is expressed as
\[ \log(\mu_{k}) = \mathbf{x}_{k}\boldsymbol{\beta} \]
Where \(\mathbf{x}_{k}\) corresponds to the vector of explanatory variables associated with the observed unit and \(\boldsymbol{\beta}\) is the vector of unknown model parameters. Equivalently, we have:
\[ \mu_{k} = E(Y_{k}\mid \mathbf{x}_{k}) = \exp(\mathbf{x}_{k}\boldsymbol{\beta}), \]
With the logarithmic link, each parameter can be interpreted as the effect of a covariate on the logarithm of the expected mean, so that \(\exp(\beta_j)\) represents the multiplicative change associated with a one-unit increase in the corresponding covariate, holding the other variables in the model constant. Under the Gamma distribution, the conditional density function of \(Y_{k}\) can be written as
\[ f(y_{k}\mid \mu_{k},\phi) = \frac{1} {\Gamma(1/\phi) (\phi\mu_{k})^{1/\phi}} y_{k}^{\frac{1}{\phi}-1} \exp\left( -\frac{y_{k}} {\phi\mu_{k}} \right), \qquad y_{k}>0, \]
where \(\Gamma(\cdot)\) denotes the Gamma function and \(\phi\) represents the dispersion parameter. Under the inference approach based on pseudo-maximum likelihood, the objective function optimized to obtain the estimators of the model parameters is given by
\[ \ell(\boldsymbol{\beta},\phi) = \sum_{k} w_{k} \log \left[ f(y_{hik}\mid\mu_{k},\phi) \right] \]
Where \(k\) represents the observation units in the sample and \(w_k\) corresponds to the expansion factor associated with each of them. This function defines the optimization criterion used to estimate the model parameters under the pseudo-maximum likelihood approach. In R, it is first necessary to define the sampling design.
gamma_design <- survey_data %>%
as_survey_design(
strata = Stratum,
ids = PSU,
weights = wk,
nest = TRUE
)The following code fits a Gamma regression model for the household income variable (Income), incorporating the complex survey design defined in the gamma_design object. The model uses household expenditure (Expenditure), age (Age), sex (Sex), and area of residence (Zone) as explanatory variables. A logarithmic link function is also specified, ensuring positive predictions and allowing the effects of the explanatory variables to be interpreted in multiplicative terms on expected income. The estimated coefficients are presented in Table 6.5:
gamma_model <- svyglm(
formula = Income ~ Expenditure + Age + Sex + Zone,
design = gamma_design,
family = Gamma(link = "log")
)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 5.409 | 0.098 | 55.184 | 0.000 |
| Expenditure | 0.002 | 0.000 | 16.891 | 0.000 |
| Age | 0.002 | 0.001 | 2.430 | 0.017 |
| SexMale | 0.048 | 0.038 | 1.260 | 0.210 |
| ZoneUrban | 0.151 | 0.089 | 1.701 | 0.092 |
Because the model uses a logarithmic link, the coefficients can be interpreted as multiplicative effects on expected income. The coefficient associated with household expenditure is positive and statistically significant (\(p<0.001\)), indicating that, holding the other variables in the model constant, a one-unit increase in expenditure is associated with an approximate \(0.2\%\) increase in expected income, since \(\exp(0.002)\approx 1.002\). Age also has a positive and significant effect (\(p=0.017\)). On average, each additional year of age is associated with an increase of around \(0.2\%\) in expected income, holding the other covariates constant. By contrast, the sex and area-of-residence variables do not show sufficient statistical evidence to conclude that they influence income.