A General Framework for Random Effects Models for Binary, Ordinal, Count Type and Continuous Dependent Variables Including Variable Selection

Random effects models are a strong tool to model the heterogeneity of clustered responses. By postulating the existence of unobserved latent variables, the so-called random effects, which are shared by the measurement within a cluster, correlation between the measurements within clusters is introduced. The clusters or units can refer to persons in repeated measurement trials or to larger units as, for example, schools with measurements referring to performance scores of students.

Detailed expositions of linear mixed models, which are typically used for continuous dependent variables, are found in Hsiao (1986), Lindsey (1993), and Jones (1993). Models for binary variables and counts are often discusses within the framework of generalized mixed models, see, for example, McCulloch and Searle (2001). Random effects models for ordinal dependent variables were considered by Harville and Mee (1984), Jansen (1990), Tutz and Hennevogl (1996) and Hartzel et al. (2001). Mixed model versions for continuous bounded data in the form of rates and proportions that take values in the interval $(0,1)$, have been considered by Qiu et al. (2008) based on the simplex model and by Bonat et al. (2015) who propagate beta distribution models. Several R packages are available to fit generalized mixed linaer models, for example, glmmTMB (Brooks et al., 2017) for various continuous and discrete distributions, lme4 for continuous and binary data, ordinal and MultOrdRS (Schauberger, 2024) for ordinal data.

Generalized mixed models within the generalized linear model framework as well as extended approaches as the models for bounded continuous data postulate familiar fixed distributions for the responses. This is different in the approach propagated here although there is some overlap with generalized linear models. The mixed thresholds model used here gains its flexibility concerning distributional assumptions by using two components, a response function, which is a distribution function, and a thresholds function that modifies the distribution. In the simplest case, by assuming a linear thresholds function, the distribution of the responses follows the response function. Thus, familiar linear Gaussian response models are obtained but also linear models with quite different, possibly skewed distribution functions are available. Distributions with a restricted support, for example if responses are observed in an interval or are positive only, are obtained by using non-linear thresholds functions. They are also useful when modeling discrete data, which typically are restricted to a specific range, for example, count data take only values $0,1,\dots$ and ordered categorical responses can be coded by $0,1,\dots,k$, where numbers only represent the order of values. In the case of binary and ordered categorical responses cumulative generalized mixed model is a special case of the proposed thresholds model.

More concrete, let $y_{i1,\dots,}y_{im}$ denote the observations on unit $i$ ($i=1,\dots,n$), which can be continuous or discrete. In addition, let $\boldsymbol{x}_{ij},\boldsymbol{z}_{ij}$ denote covariates associated with response $y_{ij}$. Then, the Mixed Thresholds Model (MTM) has the form

where $F(.)$ is a strictly increasing distribution function and $\delta_{j}(.))$ is a non-decreasing measurement-specific function defined on the support $S$ of the dependent variables, referred to as thresholds function. The predictor $\eta_{ij}=\boldsymbol{z}_{ij}^{T}\boldsymbol{b}_{i}+\boldsymbol{x}_{ij}^{T}{\boldsymbol{\beta}}_{j}-\delta_{j}(y))$ contains two components linked to explanatory variables.

• The term $\boldsymbol{z}_{ij}^{T}\boldsymbol{b}_{i}$ contains the cluster-specific effects $\boldsymbol{b}_{i}$. They are assumed to vary independently across clusters and are assumed to follow a specific distribution, typically the normal distribution, $\boldsymbol{b}_{i}\thicksim N(\boldsymbol{0},\boldsymbol{\Sigma})$.

• The term $\boldsymbol{x}_{ij}^{T}{\boldsymbol{\beta}}_{j}$ contains the effects of $\boldsymbol{x}_{ij}$ on the dependent variable. The parameters ${\boldsymbol{\beta}}_{j}$ are fixed population-specific parameters.

The distribution of the dependent variables is crucially determined by the choice of the distribution function $F(.)$, also referred to as response function, and the thresholds function $\delta_{j}(.)$. Specific choices yield models that are in common use in random effects modeling. Other choices widen the toolbox yielding models that show better fits than classical approaches. Thresholds models have been considered before in the form of item thresholds model (Tutz, 2022), which are latent trait models that aim at measurement and do not contain any covariates. The objective of the models considered here is quite different. The focus is on the effect of covariates on responses, with covariates that can take quite dfferent forms. They can be measurement-specific, unit-specific or both, which raises problems, in particular with the possible number of parameters involved.

The role the response and the thresholds function play in modeling the response distribution will become obvious when considering specific choices. In Section 2 the case of continuous responses is considered with linear and non-linear thresholds functions. Section 3 is devoted to discrete data with infinite and finite support. Joint modeling, which allows for different types of responses, in particular a mixture of continuous and discrete responses, are considered in Section 4. Marginal likelihood estimation methods are given in Section 5. In Section 6 the problem of obtaining sparser representations is addressed. Several small examples are used to demonstrate the versatility of the approach. They are meant for illustration, no in-depth investigation of effects is given.

We start with models that contain linear thresholds functions. The model class comprises the classical normal response model but allows for alternative distributions. Then we consider models for restricted support, which call for non-linear thresholds functions.

Discrete data come in two forms, with infinite support and finite support. We will consider first the case of infinite support than the case where the response is in categories. In the latter typically only ordinal scale level is assumed for the dependent variable.

A strength of the thresholds model is that dependent variables can be of various types. It allows for some of the measurements to be continous while others are binary, ordinal or given as counts. Also the combination of continuous measurements with differing support can be modeled in a joint random effects model. There have been some approaches to joint modeling of different types of responses in specificsettings, see, for example, Ivanova et al. (2016) with a focus on ordinal variables or Loeys et al. (2011), where a joint modeling approach for reaction time and accuracy in psycholinguistic experiments has been proposed. However, no general random effects model that allows to combine different types of responses seems available

The flexibility of thresholds models to account for various types of measurement is due to the general form of the model. Since the same model form, which specifies $P(Y_{ij}>y)$, applies to different types of measurement it is straightforward to obtain a joint model simply by allowing for different distributions (continuous or discrete) and specifying the thresholds function accordingly. For example, if measurement 1 is continuous and measurement 2 ordered categorical, one can choose for the first measurement the linear or logarithmic thresholds function (depending on the support) and for the second measurement the logistic thresholds function. The common random effects will account for the correlation between measurements without the need for an explicit new concept for the correlation between a continuous and an ordered categorical variable.

A general form of the threshold function term is given by $\delta_{j}(y))={\boldsymbol{\Phi}}_{j}(y)^{T}\boldsymbol{\delta}_{j}$, where ${\boldsymbol{\Phi}}_{j}(y)=(\Phi_{j0}(y),\dots,\Phi_{jM}(y))$ is a vector that contains functions of $y$. If the same threshold function is used for all measurements it can be specified as ${\boldsymbol{\Phi}}_{j}(y)^{T}=(1,g(y))$. Measurement-specific thresholds functions can be obtained by ${\boldsymbol{\Phi}}_{j}(y)^{T}=(1,g_{j}(y))$. But also more general vectors of functions can be useful. Using this parameterization, the model has the form

For continuous measurement $j$ the density is

where ${\boldsymbol{\Phi}}^{{}^{\prime}}_{j}(y)=(\Phi^{{}^{\prime}}_{j0}(y),\dots,\Phi^{{}^{\prime}}_{jM}(y))$ contains the derivatives of the components and $f(.)$ is the density linked to $F(.)$. If ${\boldsymbol{\Phi}}_{j}(y)^{T}=(1,g(y))$ one obtains the simpler form

For discrete data with $Y_{ij}\in\{0,1,\dots\}$. the discrete density is given by

When assuming that random effects are normally distributed with mean zero and covariance matrix $\boldsymbol{\mathit{\Sigma}}$ the marginal log-likelihood has the form

where $f_{0,\boldsymbol{\mathit{\Sigma}}}(.)$ is the density of the normal distribution with mean zero and covariance matrix $\boldsymbol{\mathit{\Sigma}}$. The corresponding log-likelihood is given by

Maximization of the marginal log-likelihood can be obtained by using Gauss Hermite quadrature, which has been used in generalized mixed models, see, for example, Anderson and Aitkin (1985), Liu and Pierce (1994), Hedeker and Gibbons (1994), Rodríguez (2008), Gueorguieva (2001). An alternatively is the EM, which was considered for mixed models among others by Bock and Aitkin (1981) Anderson and Aitkin (1985). Overviews on inference tools for generalized moxed models are found in Jiang and Nguyen (2007), McCulloch and Searle (2001).

As is seen from the fears data the effects of covariates can vary strongly across measurements. For example, the variable EastWest has positive values for some items but negative values for other items. People living in the Eastern part of Germany tend to show stronger fears concerning globalization and refugees but lesser fears concerning the climate crisis than people living in the Western part of Germany. Thus, the cultural imprint from living in different countries with differing experiences seems to still be present.

The inclusion of measurement-specific fixed effects ${\boldsymbol{\beta}}_{j}$ make the general model considered here very flexible. A possible downside is that the model contains many parameters. Each variable comes with as many parameters as measurements are used. However, not each variable necessarily has effects that vary across measurements. Selecting those variables that have global effects, that is effects that do not vary across measurements, would yield sparser models with an easier interpretation.

One strategy to identify variables with global effect is testing. Likelihood ratio test for hypotheses of the form

can be used to test if single variables are global. Table 3 shows the results of testing for the fears data. It is seen that all of the variables should be considered as having effects that vary over measurements.

An alternative way to obtain sparser representations is to use penalization methods that have been widely used since the advent of the Lasso (Tibshirani, 1996). Instead of maximizing the usual marginal log-likelihood one maximizes the penalized log-likelihood

where $l(.)$ is the usual log-likelihood and $P_{\boldsymbol{\lambda}}(\{{\boldsymbol{\beta}}_{j}\})$ is a penalty term depending on a tuning parameter $\lambda$.

A penalty term that enforces the building of global variables and simultaneously select variables is given by

where $||.||$ is the $L_{2}$-norm, ${\boldsymbol{\beta}}_{s,\text{dif}}^{T}=(\beta_{1s}-\beta_{2s},\beta_{2s}-\beta_{3s},\dots)$ contains all differences of effects between pairs of measurement for the $s$th variable, and ${\boldsymbol{\beta}}_{.s}^{T}=(\beta_{1s}\dots,\beta_{ms})$ is the vector that collects all parameters linked to the $s$-th variable. It is a penalty term built from two terms. The first term enforces the building of global variables. It is a grouped lasso type penalty that aims at the fusion of parameters. It enforces that whole groups of differences are set to zero with the group referring to differences of effects of single variables. In common grouped lasso penalties as considered by Yuan and Lin (2006) and Meier et al. (2008) groups refer to parameters linked to variables variables, not differences of parameters linked to variables. Fusion penalties have been used in particular to combine effects of adjacent categories when modeling the effects of categorical predictors, see Tutz and Gertheiss (2016), Gertheiss and Tutz (2023) and to enforce selection and shrinkage to global effects, for example, by Oelker and Tutz (2017). The second term is a selection term. It enforces that all parameters linked to specific variables are set to zero.

Figure 5 shows the resulting coefficient paths plotted against $\log(1+\lambda)$. For growing $\lambda$ coefficients become more similar (global) while the shrinkage toward zero (elimination of variables) sets in rather late. Zero estimates are obtained for very large values of $\lambda$, well beyond the values that are shown. This is sensible since the variables certainly have an effect on the levels of fear. The last panel shows the negative cross-validation log-likelihood (5-fold cross-validation). It is seen that the smallest levels are obtained if $\log(1+\lambda)$ is slightly above 2. In this range coefficients are still measurement-specific, which supports the testing results.

It has been demonstrated that the thresholds model is very versatile and can be adapted to quite different distribution functions. Moreover, it can be used in joint modeling of different responses, which has not been considered in the literature since appropriate models have not been available. A rare exception is the modeling of survival data where joint modeling of survival times and longitudinal data has been considered, see, for example, Hsieh et al. (2006).

We used marginal estimation based on integration methods, but alternative estimation from the toolbox of mixed models could be used. Also alternative regularization methods could be useful when trying to find sparser representations, for an overview on basic regularization method including boosting, see Hastie et al. (2009).

The class of mixed thresholds models could be made even more flexible by letting the data decide which thresholds function is appropriate. In particular, B-splines as considered extensively by Eilers and Marx (2021) could be useful, which has been demonstrated by Tutz (2022) in the item-response setting without covariates.

For simplicity, let the model be given in the form

where $\tilde{\eta}_{ij}=\boldsymbol{z}_{ij}^{T}\boldsymbol{b}_{i}+\boldsymbol{x}_{ij}^{T}{\boldsymbol{\beta}}_{j}$ and $F(.)$ is a strictly increasing distribution function.

Proof: For linear item function the thresholds model has the form $P(Y_{ij}>y|\boldsymbol{b}_{i},\boldsymbol{z}_{ij},\boldsymbol{x}_{ij})=F(\tilde{\eta}_{ij}-\delta_{j}(y))$. The corresponding distribution function is

which has the density

Thus, the mean is given by

With $\eta=\tilde{\eta}_{ij}-\delta_{0j}-\delta_{j}y$ and $d\eta/dy=-\delta_{j}$ one obtains

where $\mu_{F}=\int yf(y)dy$ is a parameter that depends on $F$ only.

The variance is given by

where $\operatorname{var}_{F}=\int({\eta-\mu})^{2}f(\eta)d\eta$.

Proof: The mean is given by

With $\eta=\tilde{\eta}_{ij}-\delta_{0j}-\delta_{j}g(y)$ one obtains

since $\beta_{0j}=\delta_{0j}-\mu_{F}$.

The variance is given by