Combine fitted models with seemingly unrelated estimation

suest() combines two or more separately fitted models into one object with a joint robust covariance matrix, so you can test whether predictions and marginal effects differ across the models. Pass the result directly to marginaleffects::predictions(), marginaleffects::avg_comparisons(), marginaleffects::avg_slopes(), and marginaleffects::hypotheses().

Usage

suest(
  ...,
  model_names = NULL,
  observation_id = NULL,
  cluster = NULL,
  weight_type = NULL,
  survey_design = NULL
)

# S3 method for class 'suest_model'
coef(object, ...)

# S3 method for class 'suest_model'
vcov(object, ...)

# S3 method for class 'suest_model'
nobs(object, ...)

# S3 method for class 'suest_model'
print(x, ...)

Arguments

…: Two or more supported fitted model objects. For backward compatibility, a character vector supplied as the third unnamed argument is interpreted as model_names for a two-model system.
model_names: Optional character vector containing one display name per model. By default, the object names supplied in the call are used.
observation_id: Optional observation identifier used to align models fitted from different data objects. Supply one or more column names found in every model’s original data, such as "id" or c("id", "wave"), or a list containing an ID vector, matrix, or data frame already aligned to each model’s estimation sample. IDs must be complete and unique within each model. The default NULL uses data-source identity and model-frame row names.
cluster: Optional cluster identifier for a joint cluster-robust covariance matrix. Supply one or more column names found in every model’s original data, or a list containing a cluster vector, matrix, or data frame aligned to each model’s estimation sample. Cluster IDs must be complete. Shared observations must have the same cluster ID in every model; disjoint observations may still share clusters across models. Supported panel systems default to the panel identifier when cluster is omitted. Supplied clusters for these models must contain whole panels.
weight_type: Optional weight interpretation. The default NULL preserves the unweighted behavior and rejects nonunit estimation weights. Use "pweight" to treat model weights as sampling weights. In version 0.1.4, pweights are supported for linear, binary logit/probit, Poisson, negative-binomial, ordered logit/probit, and multinomial logit models.
survey_design: Full common one-stage survey::svydesign() object for supported survey::svyglm() models: Gaussian identity or binary quasibinomial() logit/probit. Requires observation_id column names. Retain the design before model-specific domain subsetting. Specify weights and clusters in this design, without weight_type or cluster. Ordinary models leave this argument NULL.
object, x: A "suest_model" object.

Value

An object of class "suest_model" containing the fitted models, their joint coefficient vector, a joint model-robust covariance matrix, and sample-alignment information.

Details

Two or more models are supported. Models can use identical, partially overlapping, or completely disjoint samples. When the model calls refer to the same data source, model-frame row names identify overlapping observations. Models fitted from different data objects are treated as disjoint by default because shared observations cannot be inferred safely. Use observation_id to identify their common observations explicitly.

All combinations of the supported scalar-response models can be combined, including models whose response variables or response scales differ. All combinations of the supported categorical-response models can be combined, and scalar and categorical models may appear in the same system. Results on different response scales are labeled separately for marginaleffects.

Ordinary linear models include an ancillary lnvar parameter, matching Stata’s regress/suest parameterization. Unweighted fits use log(RSS / df.residual); pweighted fits follow suest2’s reconstructed iweight-reference normalization. Negative-binomial models include log(theta) in the joint parameter vector. Ordered and multinomial models use analytic score and observed-information calculations for stable robust covariance estimation.

Aliased parameters are not supported. Nonunit weights are rejected unless weight_type = "pweight". Pweights must be finite and strictly positive. For observations included in any pair of models, evaluated weights must agree; weights may differ for observations unique to either model, including completely disjoint samples. Pweight support is available for linear, binary logit/probit, Poisson, negative-binomial, ordered logit/probit, and multinomial logit models.

Bias-reduced, adjusted-score, Firth, and penalized GLM fits are rejected because they do not use ordinary maximum-likelihood score equations. Ordinary glm and glm2 fits must have converged. Penalized nnet::multinom fits with nonzero decay are also rejected. Beta regression requires type = "ML"; bias-reduced and bias-corrected beta fits and penalized survreg fits are not supported. Conventional survreg fits with robust = TRUE use their model-based information for the joint sandwich covariance.

Survey models

Survey support combines coefficients from Gaussian identity-link and binary quasibinomial() logit/probit svyglm() fits under one common one-stage design. Gaussian and binary fits are not mixed in the same survey system. Strata and first-stage finite-population corrections are supported. Model-specific subsets and missing outcomes are aligned using observation IDs, with zero influence outside each model’s estimation sample. The full design retains PSUs outside all model samples. Each native coefficient covariance must be reproduced before the joint matrix is returned. Survey fits contain coefficients only and do not add ancillary parameters.

Replicate-weight, multistage, two-phase, calibrated, raked, post-stratified, and PPS designs are unsupported. Lonely-PSU options are restricted to fail, remove, or certainty, with survey.adjust.domain.lonely = FALSE. A singleton stratum with a certainty FPC is supported under fail. Extra fitting weights are unsupported; place offsets in the model formula.

Predictions and effects use the joint design-based coefficient covariance. Averaging treats the supplied covariate distribution as fixed; it does not add design uncertainty from estimating that distribution. Use suest_newdata() and wts = ".suest_weight" for model-specific weighted averages. Inference uses the usual asymptotic normal default. The full design’s degrees of freedom are recorded in $survey$design_df; no automatic survey t or F adjustment is applied.

Supported models

stats::lm()
restricted survey-weighted Gaussian identity and binary quasibinomial() logit/probit models from survey::svyglm()
binary logit, probit, and complementary-log-log models from stats::glm() or glm2::glm2()
Poisson log-link models from stats::glm() or glm2::glm2()
other GLMs using identity, log, logit, probit, complementary-log-log, or log-log links
fractional-response GLMs using quasibinomial() with logit, probit, complementary-log-log, or a user-supplied log-log link
negative-binomial log-link models from MASS::glm.nb()
parametric survival and censored-regression models from survival::survreg() with a common scale, including Gaussian interval regression
beta regressions from betareg::betareg() using logit, probit, complementary-log-log, or log-log mean links
Poisson and negative-binomial zero-inflated models from pscl::zeroinfl()
truncated Gaussian regressions from truncreg::truncreg()
left-, right-, and two-limit censored Gaussian regressions from censReg::censReg(); response-scale predictions are the latent mean
unweighted two-stage least squares from fixest::feols() without absorbed fixed effects; the original data object must remain available
heteroskedastic binary probit and logit from Rchoice::hetprob()
maximum-likelihood instrumental-variable probit from Rchoice::ivpml(); response predictions use the average structural probability
bivariate probit from mvProbit::mvProbit() when both equations use the same regressors; fit with intGrad = TRUE and finalHessian = TRUE
unweighted individual fixed-effects, between-effects, and Swamy-Arora random-effects linear panel models from plm::plm(); balanced random- effects panels have the closest Stata parity, while unbalanced panels can retain small engine-specific differences. Within-model prediction uncertainty conditions on estimation-sample means and propagates slope uncertainty. At matching evaluation means, the covariance is structurally degenerate: average-level confidence intervals and hypothesis tests are unsupported, and numerical SEs can be missing or nearly zero. Slope and finite-change comparisons remain supported.
unweighted single-level random-intercept Gaussian panel models from nlme::lme() fitted with method = "ML"
unweighted individual random-intercept binary logit and probit models from pglm::pglm() fitted with model = "random", effect = "individual", and R = 12; response predictions integrate over the random effect
unweighted gamma random-effects Poisson log models from pglm::pglm() fitted with model = "random", effect = "individual", and other = "sd"; the final parameter is exposed as gamma variance alpha
unweighted binomial-logit, Poisson-log, and negative-binomial NB2 log models from glmmTMB::glmmTMB() with one grouping variable and one conditional random intercept; the final parameter is the log random-intercept standard deviation, and response predictions integrate over the Gaussian random effect. NB2 models must use the default constant dispersion model (dispformula = ~1); its estimated log size parameter log_phi precedes log_sigma, with conditional variance mu + mu^2/exp(log_phi). Weights, offsets, and zero inflation are unsupported; NB2 additionally excludes mapped or constrained parameters. Its Laplace log likelihood must exceed the zero-random-effect NB2 log likelihood at the same fixed effects and dispersion by more than sqrt(.Machine$double.eps) * max(1, abs(logLik(model))). This numerical boundary check is not a significance test and does not guarantee an interior global maximum
unweighted binomial-logit, Poisson-log, and NB2-log glmmTMB models with one correlated random intercept and numeric slope, (1 + x | id), using an unstructured covariance matrix. The slope must be a single untransformed numeric column with a syntactically valid name. The nuisance parameters are log_sd_intercept, log_sd_slope, and atanh_rho. NB2 requires the default estimated constant dispersion (dispformula = ~1), includes log_phi before these three parameters, and uses the NB2 boundary check above. Count response predictions are exp(X beta + (var_intercept + 2*x*cov_intercept_slope + x^2*var_slope)/2). Binomial response predictions integrate the logistic probability over this Gaussian variance; responses must be Bernoulli. Link predictions are X beta. Near-singular random covariance is rejected in a centered, standardized predictor basis. Weights, offsets, zero inflation, constraints, diagonal covariance, multiple slopes, slope-only terms, and other random-slope families are unsupported. Random-slope systems must contain models of the same family; the slope variable must be supplied for response predictions even when it is absent from the fixed-effects formula. Native model covariance is used as supplied. Poorly scaled predictors can produce inaccurate native numerical curvature even with convergence and a positive-definite Hessian; center and rescale continuous predictors before fitting, then compare predictions and effects after converting them to the same original units. For logit random-slope avg_slopes(), use numderiv = list("fdcenter", eps = 1e-4) and compare results at nearby steps (for example, 5e-5 and 2e-4). An unstable standard error requires further investigation; changing the finite-difference step does not repair an inaccurate native model covariance
unweighted GEE from geepack::geeglm(): Gaussian identity, binary logit/probit/cloglog, and Poisson log, with independence or exchangeable correlation and numeric outcomes
ordered logit and probit models from MASS::polr()
ordered logit and probit models from ordinal::clm() with flexible thresholds, proportional effects, and no scale model
multinomial logit models from nnet::multinom()

Examples

dat <- mtcars
dat$am <- factor(dat$am)

model1 <- glm(am ~ wt, family = binomial(), data = dat)
model2 <- glm(am ~ wt + hp, family = binomial(), data = dat)

fit <- suest(model1, model2, model_names = c("Base", "Adjusted"))
fit

Seemingly Unrelated Estimation
Models: Base + Adjusted 
Model types: Base=logit, Adjusted=logit 
Model engines: Base=stats::glm, Adjusted=stats::glm 
Comparison scale: predicted probabilities 
Observations: Base=32, Adjusted=32 
Overlapping observations: 32 
Union observations: 32 
Parameters: 5

effects <- marginaleffects::avg_comparisons(fit, variables = "wt", newdata = dat)
marginaleffects::hypotheses(effects, hypothesis = difference ~ revpairwise)

          Hypothesis Estimate Std. Error    z Pr(>|z|)    S  2.5 % 97.5 %
 (Base) - (Adjusted)   0.0779     0.0141 5.54   <0.001 25.0 0.0503  0.105