We propose a new joint mean and correlation regression model for correlated multivariate discrete responses, that simultaneously regresses the mean of each response against a set of covariates, and the correlations between responses against a set of similarity/distance measures. A set of joint estimating equations are formulated to construct an estimator of both the mean regression coefficients and the correlation regression parameters. Under a general setting where the number of responses can tend to infinity, the joint estimator is demonstrated to be consistent and asymptotically normally distributed, with differing rates of convergence due to the mean regression coefficients being heterogeneous across responses. An iterative estimation procedure is developed to obtain parameter estimates in the required, constrained parameter space. We apply the proposed model to a multivariate abundance dataset comprising overdispersed counts of 38 Carabidae ground beetle species sampled throughout Scotland, along with information about the environmental conditions of each site and the traits of each species. Results show in particular that the relationships between the mean abundances of various beetle species and environmental covariates are different and that beetle total length has statistically important effect in driving the correlations between the species. Simulations demonstrate the strong finite sample performance of the proposed estimator in terms of point estimation and inference.

Multivariate binary data are widely collected in many disciplines, including in finance, psychometrics, and ecology. Often, many binary responses are driven by only a subset of predictors, and groups of binary responses may exhibit similar effects to a predictor. Motivated by a survey containing presence-absence records of 22 demersal fish species recorded across the U.S. Northeast shelf, we propose a novel method for simultaneous coefficient clustering and variable selection in multivariate binary data using a penalized Ising regression model. The Ising regression model formulates an explicit joint distribution for a binary response vector via main and pairwise interaction coefficients, where the former is modeled as a function of covariates and the latter captures conditional dependence relationships between responses. To cluster coefficients within each covariate, and encourage sparsity across both covariates and pairwise interaction coefficients, we propose to augment the Ising regression model with adaptive fused lasso and adaptive lasso penalties. Such a structured penalty to encourage simultaneous sparsity and groupings is aligned with goals of achieving sparse species-covariate relationships, and homogeneity of environmental responses across species, in the motivating demersal fish survey. Through a reparametrization, we show that the proposed estimator can be efficiently obtained by fitting a single, adaptive lasso logistic regression model. Simulation studies and an application to the demersal fish survey demonstrate competitive performance of the proposed method relative to several existing Ising regression models for multivariate binary data, and leads to an interpretable and parsimonious set of response-covariate relationships.

Additive models offer a flexible framework for modeling nonlinear relationships between predictors and a continuous response variable, where nonparametric terms are commonly modeled via penalized splines. However, such models can be sensitive to outliers, leading to potentially unreliable estimation and inferential results. Existing methods and software for robust additive models tend to be computationally burdensome, and can come with restrictions such as permitting only one nonparametric term or not offering uncertainty quantification. In this article, we propose a fast, robust approach to fitting additive models based on the gamma-divergence, which we refer to as gamma-divergence additive models or GDAMs. Specifically, we apply gamma-divergence to the restricted maximum likelihood function of the additive model, based on treating the smoothing coefficients as random effects. This leads to an efficient minorization-maximization algorithm that adaptively downweights the impact of outlying observations via a set of normalized power density weights. Simulation studies and an application to data concerning the distribution of federal grants across the United States confirm that GDAMs perform similarly to or better than many existing (non-)robust additive modeling methods under varying degrees of contamination, while also being computationally faster and more scalable.

Generalized Linear Mixed Models (GLMMs) are widely used for analysing clustered data. One well-established method of overcoming the integral in the marginal likelihood function for GLMMs is penalized quasi-likelihood (PQL) estimation, although to date there are few asymptotic distribution results relating to PQL estimation for GLMMs in the literature. In this paper, we establish large sample results for PQL estimators of the parameters and random effects in independent-cluster GLMMs, when both the number of clusters and the cluster sizes go to infinity. This is done under two distinct regimes: conditional on the random effects (essentially treating them as fixed effects) and unconditionally (treating the random effects as random). Under the conditional regime, we show the PQL estimators are asymptotically normal around the true fixed and random effects. Unconditionally, we prove that while the estimator of the fixed effects is asymptotically normally distributed, the correct asymptotic distribution of the so-called prediction gap of the random effects may in fact be a normal scale-mixture distribution under certain relative rates of growth. A simulation study is used to verify the finite sample performance of our theoretical results.

Multivariate correlated outcomes occur across disciplines, including ecology, social sciences, and psychometrics. This paper focuses on clustering these outcomes across observational units, specifically, finding groups of units with the same “outcome profile". Our motivation comes from bioregionalization in ecology, which aims to cluster sites into bioregions with the same species profiles, where site membership can depend on environmental or habitat covariates. To accomplish this, we propose finite mixtures of generalized estimating equations (MixGEE). Unlike existing approaches to model-based bioregionalization, MixGEE partitions sites into regions while accounting for between-species correlations through a region-specific working correlation structure. Thus, each region is characterized by a marginal species mean vector and a between-species correlation matrix. Unlike likelihood-based finite mixture models, MixGEE does not require a full joint distribution of the multivariate outcomes. Instead, we construct a pseudo-posterior probability for region membership motivated by the large-sample distribution of the estimating equation. This leads to an iterative algorithm alternating between updating these probabilities and solving weighted estimating equations. We determine the number of regions using cross-validation based on predictive performance for held-out sites and species components, and use a clustered Dirichlet random-weight bootstrap for uncertainty quantification. Simulations demonstrate reliable estimation and inference under various correlation structures and more stable selection of the number of groups than methods that ignore dependence. Applying MixGEE to presence–absence records of fish species around the Kerguelen Plateau reveals three distinct fish assemblage profiles with heterogeneous occurrence and within-site correlation patterns.

In insurance markets, claim costs are highly variable, heavy-tailed, and difficult to predict. At the same time, policyholder retention and lapse behavior (customer churn) are critical determinants of long-term profitability and solvency. Most existing models in the literature treat claim costs and lapses as independent, overlooking potential latent associations that arise from adverse selection and unobserved heterogeneity. In this article, we introduce a joint modeling framework that simultaneously captures individual-level claim costs and churn behavior, using multivariate Tweedie regression with shared random effects. This framework integrates claim cost dynamics with lapse risk, allowing insurers to more accurately predict costs, classify policyholder profitability, and design retention or pricing strategies. Applying our approach to data from the Wisconsin Local Government Property Insurance Fund, we demonstrate that accounting for dependence between claim risk and lapse risk improves out-of-sample prediction and yields actionable insights for customer valuation and management.

When fitting generalized linear mixed models, choosing the random effects distribution is an important decision. As random effects are unobserved, misspecification of their distribution is a real possibility. Thus, the consequences of random effects misspecification for point prediction and prediction inference of random effects in generalized linear mixed models need to be investigated. A combination of theory, simulation, and a real application is used to explore the effect of using the common normality assumption for the random effects distribution when the correct specification is a mixture of normal distributions, focusing on the impacts on point prediction, mean squared prediction errors, and prediction intervals. Results show that the level of shrinkage for the predicted random effects can differ greatly under the two random effect distributions, and so is susceptible to misspecification. Also, the unconditional mean squared prediction errors for the random effects are almost always larger under the misspecified normal random effects distribution, while results for the mean squared prediction errors conditional on the random effects are more complicated but remain generally larger under the misspecified distribution (especially when the true random effect is close to the mean of one of the component distributions in the true mixture distribution). Results for prediction intervals indicate that the overall coverage probability is, in contrast, not greatly impacted by misspecification. It is concluded that misspecifying the random effects distribution can affect prediction of random effects, and greater caution is recommended when adopting the normality assumption in generalized linear mixed models.

Linear mixed models (LMMs) are a popular class of methods for analyzing longitudinal and clustered data. However, such models can be sensitive to outliers, and this can lead to biased inference on model parameters and inaccurate prediction of random effects if the data are contaminated. We propose a new approach to robust estimation and inference for LMMs using a hierarchical gamma-divergence, which offers an automated, data-driven approach to downweight the effects of outliers occurring in both the error and the random effects, using normalized powered density weights. For estimation and inference, we develop a computationally scalable minorization-maximization algorithm for the resulting objective function, along with a clustered bootstrap method for uncertainty quantification and a Hyvarinen score criterion for selecting a tuning parameter controlling the degree of robustness. Under suitable regularity conditions, we show the resulting robust estimates can be asymptotically controlled even under a heavy level of (covariate-dependent) contamination. Simulation studies demonstrate hierarchical gamma-divergence consistently outperforms several currently available methods for robustifying LMMs. We also illustrate the proposed method using data from a multi-center AIDS cohort study. Supplementary materials for this article are available online.

We introduce Gradient Boosted Mixed Models (GBMixed), a framework which extends boosting to clustered data by jointly modeling the mean and variance components in a linear mixed model via likelihood-based gradients. GBMixed estimates a nonparametric fixed effects function characterizing the overall mean of the response, while also allowing the random effects covariance matrix along with the residual variance to depend on covariates in a flexible manner. We demonstrate how GBMixed facilitates covariate-dependent random effect predictions, and subsequently point predictions and prediction intervals for individual treatment effects, that can adapt between population-level and cluster-level information. Simulations and applications to two real-world datasets demonstrate that GBMixed can accurately recover complex nonlinear fixed effect functions and covariate-dependent covariances in a linear mixed model, while also improving point and probabilistic predictive performance compared with several existing approaches such as parametric linear mixed models, Natural Gradient Boosting, and Gaussian Process Boosting.

Sufficient dimension reduction (SDR) is a popular class of regression methods which aim to find a small number of linear combinations of covariates that capture all the information of the responses i.e., a central subspace. The majority of current methods for SDR focus on the setting of independent observations, while the few techniques that have been developed for clustered data assume the linear transformation is identical across clusters. In this article, we introduce random effects SDR, where cluster-specific random effect central subspaces are assumed to follow a distribution on the Grassmann manifold, and the random effects distribution is characterized by a covariance matrix that captures the heterogeneity between clusters in the SDR process itself. We incorporate random effect SDR within a model-based inverse regression framework. Specifically, we propose a random effects principal fitted components model, where a two-stage algorithm is used to estimate the overall fixed effect central subspace, and predict the cluster-specific random effect central subspaces. We demonstrate the consistency of the proposed estimators, while simulation studies demonstrate the superior performance of the proposed approach compared to global and cluster-specific SDR approaches. We also present extensions of the above model to handle mixed predictors, demonstrating how random effects SDR can be achieved in the case of mixed continuous and binary covariates. Applying the proposed methods to study the longitudinal association between the life expectancy of women and socioeconomic variables across 117 countries, we find log income per capita, infant mortality, and income inequality are the main drivers of a two-dimensional fixed effect central subspace, although there is considerable heterogeneity in how the country-specific central subspaces are driven by the predictors.

Background Over the past decade, joint species distribution models (JSDMs) and model-based ordination have emerged as powerful tools for the analysis of community ecology data. Generalized linear latent variable models (GLLVMs) offer a flexible framework for multivariate analysis of a wide range of data types, based on including a small number of latent variables to perform dimension reduction while accounting for residual correlation between species. Fast estimation methods The R package gllvm implements a wide range of GLLVMs, with estimation performed via fast approximate likelihood-based techniques; including the recently proposed extended variational approximation, which is applicable to almost any combination of response type and link function. Since its original development and accompanying software paper, the gllvm package has undergone a significant overhaul, consolidating its place as a general framework for joint modeling of community ecology datasets. Expanded functionalities Some of the key new features of gllvm include model-based constrained and concurrent ordination methods, capacity to account for nested/hierarchical sampling designs, and (phylogenetic) random effects. On top of this, other notable improvements include a great expansion of the response types that it can handle, enhanced capabilities of GLLVM inference, selection and prediction, and an easier-to-use interface for model fitting.

In fisheries ecology, species abundance data are often collected by multiple surveys, each with unique characteristics. This article focuses on Atlantic sea scallop abundance data along the northeast coast of the United States, collected from two bottom trawl surveys which cover a larger spatial domain but have low catch efficiency, and a dredge survey which is more efficient but limited to domains where the species are believed to be present. To model such data, integrated species distribution models (ISDMs) have been proposed to incorporate information from multiple surveys, by including common environmental effects along with correlated survey-specific spatial fields. However, while flexible, these ISDMs can be susceptible to overfitting, which can complicate interpretability of the shared environmental effects and potentially lead to poor predictive performance. To overcome these drawbacks, we introduce a novel single index ISDM, built from a single index (with spatial random effects) that represents a latent measure of the true species distribution, and survey-specific catch efficiency functions which map the single index to the survey-specific expected catch. Our results show that the single index ISDM offers more meaningful interpretations of the environmental effects and survey catch efficiency differences, while potentially achieving better predictive performance than existing ISDMs.

Modeling directional time series data such as wind or ocean current direction presents several interesting challenges. Standard linear time series techniques do not account for the circularity of the observations, while existing circular modeling approaches typically work best when the entirety of the data only spans a small arc. Motivated by a series of fire burn experiments collecting wind direction in south eastern Australia, we propose a new method for directional time series data by combining hidden Markov models with a conditional circular distribution given the latent state. The resulting circular hidden Markov model, or cHMM, can allow for multimodality and/or varying amounts of circular dispersion over time. Furthermore, by utilizing a von Mises distribution whose mean direction depends on previous observations, we can accommodate strong serial correlations within a specific hidden state. We employ direct maximum likelihood estimation to fit the cHMM, making use of recursive probability formulas to efficiently compute the marginal log-likelihood function, and examine three approaches to perform forecasting based on extrapolating the latent state sequence and then direction observations conditional on this sequence. Simulation studies demonstrate the validity of our estimation procedure, while an application to three motivating wind direction datasets from fire burn experiments reveals that cHMMs produces similar or better point/probabilistic forecasting performance compared with several established time series methods.

Generalized estimating equations (GEEs) are a popular statistical method for longitudinal data analysis, requiring specification of the first 2 marginal moments of the response along with a working correlation matrix to capture temporal correlations within a cluster. When it comes to prediction at future/new time points using GEEs, a standard approach adopted by practitioners and software is to base it simply on the marginal mean model. In this article, we propose an alternative approach to prediction for independent cluster GEEs. By viewing the GEE as solving an iterative working linear model, we borrow ideas from universal kriging to construct an adjusted predictor that exploits working cross-correlations between the current and new observations within the same cluster. We establish theoretical conditions for the adjusted GEE predictor to outperform the standard GEE predictor. Simulations and an application to longitudinal data on the growth of sitka spruces demonstrate that, even when we misspecify the working correlation structure, adjusted GEE predictors can achieve better performance relative to standard GEE predictors, the so-called "oracle" GEE predictor using all time points, and potentially even cluster-specific predictions from a generalized linear mixed model.

We study three commonly applied measures of uncertainty for random effects prediction in generalized linear mixed models (GLMMs), namely the unconditional and conditional mean squared errors of prediction (UMSEP and CMSEP, respectively), and the unconditional variance of the prediction gap used by the popular R package for glmmTMB. We demonstrate that, although the three theoretical measures differ in how they quantify uncertainty, the resulting estimators all turn out to be very similar in form. We derive asymptotic results regarding the consistency of the three measures of uncertainty, and in doing so resolve a contradiction between theoretical and empirical results for the glmmTMB variance estimator by re-interpreting it conditionally on a finite subset of the random effects. Our results have important implications for predictive inference in GLMMs, particularly around the legitimacy and implications of coupling these measures with a normality assumption to construct prediction intervals for the random effects.

1. Joint species distribution models (JSDMs) have gained considerable traction among ecologists over the past decade, due to their capacity to answer a wide range of questions at both the species- and the community-level. The family of generalized linear latent variable models in particular has proven popular for building JSDMs, being able to handle many response types including presence-absence data, biomass, overdispersed and/or zero-inflated counts. 2. We extend latent variable models to handle percent cover data, with vegetation, sessile invertebrate, and macroalgal cover data representing the prime examples of such data arising in community ecology. 3. Sparsity is a commonly encountered challenge with percent cover data. Responses are typically recorded as percentages covered per plot, though some species may be completely absent or present, i.e., have 0% or 100% cover respectively, rendering the use of beta distribution inadequate. 4. We propose two JSDMs suitable for percent cover data, namely a hurdle beta model and an ordered beta model. We compare the two proposed approaches to a beta distribution for shifted responses, transformed presence-absence data, and an ordinal model for percent cover classes. Results demonstrate the hurdle beta JSDM was generally the most accurate at retrieving the latent variables and predicting ecological percent cover data.

Restricted maximum likelihood (REML) estimation is a widely accepted and frequently used method for fitting linear mixed models, with its principal advantage being that it produces less biased estimates of the variance components. However, the concept of REML does not immediately generalize to the setting of non-normally distributed responses, and it is not always clear the extent to which, either asymptotically or in finite samples, such generalizations reduce the bias of variance component estimates compared to standard unrestricted maximum likelihood estimation. In this article, we review various attempts that have been made over the past four decades to extend REML estimation in generalized linear mixed models. We establish four major classes of approaches, namely approximate linearization, integrated likelihood, modified profile likelihoods, and direct bias correction of the score function, and show that while these four classes may have differing motivations and derivations, they often arrive at a similar if not the same REML estimate. We compare the finite sample performance of these four classes, along with methods for REML estimation in hierarchical generalized linear models, through a numerical study involving binary and count data, with results demonstrating that all approaches perform similarly well reducing the finite sample size bias of variance components. Overall, we believe REML estimation should more widely adopted by practitioners using generalized linear mixed models, and that the exact choice of which REML approach to use should, at this point in time, be driven by software availability and ease of implementation.

In many applications of multivariate longitudinal mixed models, it is reasonable to assume that each response is informed by only a subset of covariates. Moreover, one or more responses may exhibit the same relationship to a particular covariate for example, if they are capturing the same underlying aspect of an individual physical, mental, and emotional health. To address the above challenges, we propose a method for simultaneous clustering and variable selection of fixed effect coefficients in multivariate mixed models. We achieve this in a computationally scalable manner via a composite likelihood approach: separate mixed models are first fitted to each response, after which the model estimates are combined into a single quadratic form resembling a multivariate Wald statistic. We then augment this with fusion- and sparsity-inducing penalties based on broken adaptive ridge regression. Simulation studies demonstrate that the proposed composite quadratic estimator is similar to or better than several existing techniques for fixed effects selection in (univariate) mixed models while being computationally much more efficient. We apply the proposed method to longitudinal panel data from Australia to quantify how an individual's overall health, assessed via a set of eight composite scores, evolves as a function of various demographic and lifestyle variables.

The ecological niche is a fundamental concept in ecology that can be used in order better understand species relationships. The overlap in species niches provides a measure of the likelihood for species to co-occur. Most approaches that quantify niche overlap have been based on distance and similarity indices, for pairwise combinations of species. In this paper, we suggest that niche overlap can be calculated from the predictions of a model. Using a statistical model to predict niche overlap provides various benefits, includes the possibility to adjust the model to properties of the data. We demonstrate this using an example dataset of an ecological community of Foraminifera species, to which we fit a generalized linear latent variable model (GLLVM). GLLVMs are a flexible class of models that allow to estimate the distribution of species using both measured environmental predictors and residual covariation between species. We demonstrate how to calculate niche overlap from GLLVMs for any combination of species, and separately for different environments. Predicting niche overlap from a model further expands the toolset available to ecologists for the exploration of species co-occurrence patterns.

Structural equation models (SEMs) are commonly used to study the structural relationship between observed variables and latent constructs. Recently, Bayesian fitting procedures for SEMs have received more attention thanks to their potential to facilitate the adoption of more flexible model structures, and variational approximations have been shown to provide fast and accurate inference for Bayesian analysis of SEMs. However, the application of variational approximations is currently limited to very simple, elemental SEMs. We develop mean-field variational Bayes algorithms for two SEM formulations for data that present non-Gaussian features such as skewness and multimodality. The proposed models exploit the use of mixtures of Gaussians, include covariates for the analysis of latent traits and consider missing data. We also examine two variational information criteria for model selection that are straightforward to compute in our variational inference framework. The performance of the MFVB algorithms and information criteria is investigated in a simulated data study and a real data application.