Source-linked AI summary
mgm: Estimating Time-Varying Mixed Graphical Models in High-Dimensional Data
Jonas M. B. Haslbeck, Lourens J. Waldorp
TL;DR
mgm addresses the need to model systems containing mixed variable types and evolving relationships over time. It provides estimation, prediction, stability assessment, and reproducible examples for mixed graphical and vector autoregressive models, including time-varying forms. The package extends related implementations with k-order and mixed time-varying models, while its interpretation is constrained by certain penalty and sampling limitations.
Problem
Existing graphical-model tools are limited for datasets combining variable types and for systems whose relationships change over time.
Method
mgm uses penalized generalized linear-model neighborhood selection to estimate stationary and kernel-weighted time-varying k-order MGMs and mixed VAR models.
Results
mgm provides estimation, prediction, resampling-based stability assessment, flexible sampling, qgraph-compatible outputs, and fully reproducible usage examples.
Takeaways & Limitations
The package supports analysis of mixed-variable systems across time and can represent temporally evolving dependencies relevant to organizational, biological, medical, and psychiatric processes.
Takeaways & Limitations
Different penalties for different edge types are not currently supported, and declining local sample size can reduce sensitivity in later time points of time-varying analyses.
Abstract
from arXiv · showhide
We present the R-package mgm for the estimation of k-order Mixed Graphical Models (MGMs) and mixed Vector Autoregressive (mVAR) models in high-dimensional data. These are a useful extensions of graphical models for only one variable type, since data sets consisting of mixed types of variables (continuous, count, categorical) are ubiquitous. In addition, we allow to relax the stationarity assumption of both models by introducing time-varying versions MGMs and mVAR models based on a kernel weighting approach. Time-varying models offer a rich description of temporally evolving systems and allow to identify external influences on the model structure such as the impact of interventions. We provide the background of all implemented methods and provide fully reproducible examples that illustrate how to use the package.
1. Introduction
mgm is an R package for estimating mixed graphical and vector autoregressive models, including time-varying versions, in data containing multiple variable types. It extends existing implementations with broader model coverage, prediction, stability assessment, visualization compatibility, and reproducible examples.
- Motivation and contribution: mgm estimates k-order Mixed Graphical Models and mixed Vector Autoregressive models for mixed-variable data, including stationary and time-varying versions.The package also supports mVAR models with specified lags.
- Functionality and implementation: The package provides estimation algorithms, prediction, nodewise errors, resampling-based stability assessment, sampling functions, and qgraph-compatible outputs.It supplies reproducible examples and uses qgraph for visualization rather than providing its own plotting functions.
- Motivation and contribution: Mixed-variable datasets combine continuous, count, categorical, binary, ordinal, and skewed continuous variables in applications such as marketing and medicine.These data structures motivate extensions beyond graphical models designed for a single variable type.
- Motivation and contribution: Time-varying MGMs describe relations at the same time point, while mVAR models describe relations across specified time lags.Allowing parameters to vary over time supports analysis of evolving systems and potential intervention impacts.
- Functionality and implementation: The package is implemented in R and uses glmnet to fit penalized generalized linear models for neighborhood selection.glmnet is implemented in Fortran and optimized for computational efficiency.
- Related implementations: mgm extends related packages by supporting k-order MGMs, mVAR models with arbitrary lag sets, prediction, resampling-based stability assessment, and time-varying mixed models.The paper describes mgm as the first package supporting time-varying MGMs and mVAR models.
2. Background
The paper develops graphical-model foundations and extends them to mixed-variable, higher-order, and time-varying models estimated through the mgm framework. It also describes regularization choices and computational considerations for these models.
- Graphical-model foundations: Undirected graphical models encode conditional independence statements through graph separation and the Global Markov property.For strictly positive exponential-family distributions, the Global Markov property is equivalent to Markov factorization into clique functions.
- Mixed Graphical Models: Mixed Graphical Models combine arbitrary univariate exponential-family conditional distributions, including variables with different measurement types.Their canonical parameters depend on products of sufficient-statistic functions of neighboring variables, with interaction order bounded by k.
- Mixed Graphical Models: For continuous variables with pairwise interactions, an MGM reduces to a multivariate Gaussian model parameterized by intercepts and a matrix of partial correlations.Higher-order interactions add further parameters, while categorical variables can require several parameters per interaction.
- Scope and limitations: mgm does not implement normalizability constraints for mixed densities because its underlying glmnet implementation cannot specify them.The paper notes that constrained estimation algorithms exist for pairwise MGMs.
- Estimation and regularization: The EBIC selects the regularization parameter using fit, model size, dimensionality, and a tuning parameter γ that trades sensitivity against precision.The paper adopts γ = 0.25 as a default, while recommending simulation-based selection when the application-specific trade-offs are important.
- Time-varying models: For each estimation point te, time-varying estimation produces a parameter vector θ_te and graph G_te, with complexity O(|E|p log(p^2k−1)).The computational cost scales with the number of estimation points as well as the number of variables and interaction order.
3. Usage and Examples
The examples show how mgm supports installation, simulation, estimation, prediction, visualization, and stability assessment for mixed graphical models. A pairwise MGM example demonstrates specifying variable types, interaction order, regularization, and recovered dependencies.
- Usage workflow: The examples cover sampling, estimation, prediction, visualization, and stability assessment for the four supported model types.Sampling functions enable simulation-based evaluation, and the analyses are fully reproducible.
- Estimating a pairwise MGM: A pairwise MGM is estimated from mixed data by specifying variable types, category levels, interaction order, and regularization selection.The example uses k = 2 and selects the penalty parameter by 10-fold cross-validation.
- Estimating a pairwise MGM: mgm() returns weighted adjacency matrices, interaction parameters, intercepts, and nodewise regression models.The pairwise output also stores signs when they are defined.
- Results: The example correctly recovers pairwise dependencies 1-4, 2-3, and 1-2.The interaction 1-4 has weight 0.4676953 and positive sign.
- Interpreting interactions: Signs are defined by default for interactions between non-categorical variables, whereas interactions involving categorical variables with more than two categories have no sign.Categorical interactions can involve multiple parameters, so their overall sign is undefined.
- Regularization: The penalty can be changed from the default ℓ1 penalty to elastic-net or ℓ2 regularization through alphaSeq.alphaSeq = 1 gives ℓ1 regularization, while alphaSeq = 0 gives ℓ2 regularization.
Making Predictions from Mixed Graphical Models
The predict() function generates predictions and nodewise errors for mixed graphical models using user-selected error functions. The example reports separate measures for continuous and categorical variables and illustrates their visualization with the estimated network.
- Prediction outputs: predict() takes an estimated model and same-format data, then returns predictions, categorical probabilities, and nodewise prediction errors.Users can select error functions separately for continuous and categorical variables.
- Prediction errors: Normalized accuracy is defined as 1 − max{p0, p1, . . . , pm}, where the p values are marginal category probabilities.For binary variables, p1 = 1−p0.
- Prediction errors: Continuous variables are evaluated with RMSE and R2, while categorical variables are evaluated with accuracy and normalized accuracy.Additional customary error functions can also be supplied.
- Out-of-sample evaluation: Using the estimation data for prediction produces within-sample errors, whereas generalization requires fresh test data in the same format.New data can be passed directly to predict().
- Visualization: The estimated network can display interaction parameters together with nodewise errors using qgraph.Edge colors encode interaction signs, edge widths encode parameter magnitude, and node rings display selected error measures.
- Visualization: In the illustrated network, the green edge between variables 1 and 2 denotes a positive Gaussian relationship, while grey edges mark categorical interactions without defined signs.Edge width is proportional to the corresponding edge-parameter size.
Bootstrap Sampling Distributions
The package uses nonparametric bootstrap resampling to characterize edge-weight stability, while illustrating reproducible data generation for mixed-variable MGMs.
- Bootstrap procedure: The resample() function obtains empirical sampling distributions for MGM parameter estimates using B nonparametric bootstrap samples.It accepts a fitted MGM, data, the number of bootstrap samples, and requested lower and upper quantiles.
- Bootstrap procedure: The bootstrap output stores edge-weight distributions in a B×p×p array alongside quantiles, fitted models, function calls, and per-sample runtimes.For example, the distribution for edge 3-4 is stored in res_obj$bootParameters[, 3, 4].
- Results: Edges 1-2 and 1-4 were far from zero, had small standard deviations, and were nonzero in 100% of bootstrap estimates.The plot summarizes each edge using the arithmetic mean, the proportion of nonzero estimates, and the 0.05 and 0.95 quantiles.
- Results: Edges 3-4, 2-3, and 1-3, absent from the true graph, had sampling distributions near zero and substantially smaller nonzero-estimate proportions.This contrasts their bootstrap behavior with the recovered edges 1-2 and 1-4.
- Reproducible example: The simulated MGM example uses two Gaussian nodes, two categorical nodes, and three pairwise interactions, with mgmsampler() generating 500 observations.The example specifies variable types, levels, thresholds, Gaussian standard deviations, interaction factors, and reproducible sampling settings.
Application: Autism and Well-being
The package estimates and visualizes a pairwise MGM for 3,521 individuals with autism-spectrum diagnoses, linking variables across demographic, psychological, social, and medical domains.
- Dataset: The autism dataset contains responses from 3,521 individuals to 28 questions spanning demographic, psychological, social-environment, and medical measurements.The data include mixed variable types.
- Model estimation: The example estimates a pairwise MGM with k = 2 and selects regularization parameters using EBIC with γ = 0.25.The fitted model is visualized through its weighted adjacency matrix, edge colors, and domain groupings.
- Visualization: The visualization uses green for positive relationships, red for negative relationships, grey for categorical interactions without a defined sign, and edge width for absolute edge weight.Node colors identify Demographics, Psychological, Social Environment, and Medical domains.
- Findings: Age and age of diagnosis show a strong positive relationship, while unfinished educations and work satisfaction show a plausible negative relationship.These examples illustrate how the estimated graph can expose interpretable associations.
- Findings: Well-being is strongly connected to satisfaction with social contacts and integration in society.All three variables are categorical, so showInteraction() can inspect their multi-parameter interactions in detail.
Estimating higher-order Mixed Graphical Models
The package supports k-order MGMs that retain higher-order interaction structure and provides factor-graph and conditional-table visualizations for interpreting those interactions.
- Model specification: Setting k = 3 estimates all pairwise and three-way interactions in the PTSD symptom dataset, with categorical states modeled explicitly under overparameterization.The standard and overparameterized representations are statistically equivalent, but may support different interpretations.
- Estimated structure: The fitted model contained two nonzero pairwise interactions and five nonzero three-way interactions, including Dreams–Flashbacks–Upset.The interaction list stores indicators, strengths, and the parameters defining each interaction.
- Factor-graph interpretation: A factor graph distinguishes variable nodes from interaction factors, preserving whether dependencies arise from pairwise or three-way interactions.Circle nodes represent variables, while factor nodes connect to two or three variables; edge width represents interaction strength.
- Three-way interaction: The Dreams–Upset relationship considerably depends on whether Flashbacks is present or absent.Conditioning the cross-table on Flashbacks reveals moderation of the lower-order relationship by the third variable.
3.2. Stationary mixed VAR models
The package fits mixed VAR models for heterogeneous time series, recovers simulated lagged dependencies, computes nodewise predictions, and visualizes effects with prediction errors.
- Model estimation: The mixed VAR example contains four categorical and two Gaussian variables, with lag-1 effects from 6 to 5, 5 to 1, and 3 to 1.The model uses EBIC with γ = .25 to select regularization parameters.
- Model structure: The mVAR output stores lagged edge weights, signs when defined, and the full parameterization across specified lags.Diagonal entries represent autoregressive effects, while off-diagonal entries represent cross-lagged effects; columns predict rows.
- Recovery results: All three true cross-lagged effects were recovered, and all other effects were correctly set to zero.This result is reported for the simulated mixed VAR dataset.
- Prediction: The predict() function returns predictions, categorical probabilities, and nodewise errors using RMSE, explained variance, and normalized accuracy measures.Users can also supply customary error functions through errorCon and errorCat.
- Prediction results: Node 1 had the highest normalized accuracy, nodes 2–4 had normalized accuracy of 0, and node 5 had nonzero explained variance because node 6 predicted it.Node 6 had zero explained variance because no other node predicted it.
- Visualization: Figure 7 combines lagged interaction parameters with nodewise errors, using green for positive effects, grey for undefined signs, and width for absolute edge weight.The visualization shows explained variance for Gaussian variables and normalized accuracy for categorical variables.
Sampling from mixed VAR model
The mVAR sampling example constructs a mixed model with six variables, one lag, categorical and Gaussian nodes, and explicitly specified lagged effects. The fitted model is visualized separately by lag, with many nonzero lag-one coefficients in the fMRI example.
- Model specification: The sampled mVAR model contains six variables: four categorical variables with levels 2, 2, 4, and 4, plus two Gaussian variables.The model uses one lag of order one and samples N = 200 observations.
- Model specification: Gaussian nodes receive one threshold and a standard deviation, whereas categorical nodes receive separate thresholds for their categories.Thresholds correspond to intercept terms, and standard-deviation entries for non-Gaussian variables are ignored.
- Lagged effects: The coefficient array stores lagged effects across source nodes, target nodes, category levels, and the specified lag.The example includes continuous-to-continuous, continuous-to-binary, and categorical-to-categorical effects.
- Fitted example: The fitted resting-state fMRI mVAR model uses lags 1, 2, and 3 and represents 68 × 68 × 3 interaction parameters in separate lag-specific network plots.All 68 variables are modeled as conditional Gaussians.
- Fitted example: Many coefficients are nonzero at lag one, whereas only a few are nonzero at lags two and three.The example notes that fitting only lag one could reduce variance in further analyses.
3.3. Time-varying Mixed Graphical Model
Time-varying MGMs estimate local networks by combining observations according to a Gaussian-kernel bandwidth. The bandwidth is selected by time-stratified cross-validation and then used to estimate models at chosen time points.
- Data and motivation: The time-varying MGM example analyzes 150 gene expressions measured at 67 time points across the fruit fly life span in a high-dimensional setting.The number of variables exceeds the number of observations.
- Bandwidth selection: The bandwidth σ controls how many nearby time points influence each local model through Gaussian weighting on the normalized interval [0, 1].Smaller bandwidths emphasize observations close to the estimation point, whereas larger bandwidths distribute weight more broadly.
- Bandwidth selection: Time-stratified cross-validation leaves out each estimation-time observation, predicts it with candidate models, and evaluates RMSE or 0/1-loss.The procedure can use specified fold and sequence sizes to reduce computational cost.
- Bandwidth selection: For the fruit-fly dataset, σ = 0.3 minimizes the prediction error among the candidate bandwidths.If an endpoint candidate wins, the candidate sequence should be extended in that direction.
- Final estimation: Using the selected bandwidth, the final model is estimated at 20 equally spaced points on [0, 1], with pairwise parameters stored as time-indexed arrays.The example uses bandwidth = 0.3 and stores the third estimation point separately from the others.
Making Predictions from time-varying Mixed Graphical Model
Predictions from time-varying MGMs can be computed without fitting a separate model at every measurement time. The example also shows that changing edge counts may reflect declining local sample size and sensitivity rather than true network changes.
- Prediction strategies: Estimating a model at every time point would require n fitted models, so mgm provides two lower-cost prediction strategies.These strategies operate from models estimated at selected time points.
- Prediction strategies: The weighted method averages predictions from all estimation-point models using their kernel weights at each time point.The closestModel method instead uses the model estimated at the nearest estimation point.
- Prediction output: The prediction output includes predicted values, categorical probabilities, the true data matrix, and local nodewise errors indexed by estimation point.The example invokes predict() with tvMethod = "weighted".
- Fruit-fly application: In the fruit-fly analysis, the estimated edge count decreases across time as the local sample size becomes very small near the series end.The figure compares edge count, local sample size, measurement availability, and networks at three estimation points.
- Fruit-fly application: Unequal sampling and sparse late-time data reduce sensitivity, making it difficult to determine whether late edges are truly absent or undetected.The paper highlights roughly constant sampling frequency as important for interpreting temporal changes.
Making Predictions from time-varying mixed VAR model
Time-varying mixed VAR predictions use either kernel-weighted models or the closest estimation-point model, while the symptom-data example visualizes changing lagged effects over time. The analysis uses only measurements with sufficient preceding observations.
- Prediction: Predictions from a time-varying mVAR model require choosing between weighted averaging across estimation-point models and the closestModel strategy.The prediction call also supplies day and notification information describing measurement consecutiveness.
- Prediction: The prediction output contains predicted values, categorical probabilities, the true data matrix, and local nodewise errors indexed by estimation point.This output structure parallels the time-varying MGM prediction output.
- Visualization: The symptom-data networks change considerably over time, including a Tired autoregressive effect that is strong initially and decreases almost monotonously.The authors interpret this pattern as suggesting that a stationary model is not appropriate for these data.
4. Concluding Comments
The mgm package supports stationary and time-varying mixed graphical and vector autoregressive models, with tools for prediction, stability assessment, simulation, and reproducible examples. Planned extensions target higher-order interactions, tuning-parameter selection, and alternative estimators.
- mgm fits stationary and time-varying Mixed Graphical Models and mixed Vector Autoregressive Models.
- Prediction, nodewise-error, and resampling methods support model assessment and stability evaluation.
- Flexible sampling functions enable simulation-based evaluation of estimation algorithms in specified situations.
- Planned extensions include accessible higher-order interaction inspection, additional tuning-parameter selection methods, and estimators beyond ℓα-penalized regression.
Journal of Statistical Software http://www.jstatsoft.org/
The supplied passages identify the Journal of Statistical Software and provide publication metadata fields for volume, issue, submission, and acceptance dates.
- The publication is associated with the Journal of Statistical Software website.
- The journal metadata includes fields for volume and issue numbers.
- The metadata includes submitted and accepted date fields.