Quick reference for supported objects and options

Here is an alphabetical list of model classes that are supported, andthe arguments that apply. Detailed documentation follows, with objectsgrouped by the code in the “Group” column. Scroll down or follow thelinks to those groups for more information.

Group A – “Standard” or minimally supported models

Models in this group, such aslm, do not have unusualfeatures that need special support; hence no extra arguments are needed.Some may requiredata in the call.

Group B – Beta regression

The additionalmode argument forbetaregobjects has possible values of"response","link","precision","phi.link","variance", and"quantile", which have thesame meaning as thetype argument inpredict.betareg – with the addition that"phi.link" is like"link", but for theprecision portion of the model. Whenmode = "quantile" isspecified, the additional argumentquantile (a numericscalar or vector) specifies which quantile(s) to compute; the default is0.5 (the median). Also in"quantile" mode, an additionalvariablequantile is added to the reference grid, and itslevels are the values supplied. Modes"variance" and"quantile" cause an implied re-grid.

Back to quick reference

Group C – Count models

Two optional arguments –mode andlin.pred– are provided. Themode argument has possible values"response" (the default),"count","zero", or"prob0".lin.pred islogical and defaults toFALSE.

Withlin.pred = FALSE, the results are comparable tothose returned bypredict(..., type = "response"),predict(..., type = "count"),predict(..., type = "zero"), orpredict(..., type = "prob")[, 1]. See the documentation forpredict.hurdle andpredict.zeroinfl. Note thatspecifyinglin.pred = FALSE causes re-gridding to takeplace.

The optionlin.pred = TRUE only applies tomode = "count" andmode = "zero". The resultsreturned are on the linear-predictor scale, with the same transformationas the link function in that part of the model. The predictions for areference grid withmode = "count",lin.pred = TRUE, andtype = "response" will bethe same as those obtained withlin.pred = FALSE andmode = "count"; however, any EMMs derived from these gridswill be different, because the averaging is done on the log-count scaleand the actual count scale, respectively – thereby producing geometricmeans versus arithmetic means of the predictions.

If thevcov. argument is used (see details in thedocumentation forref_grid), it must yield a matrix of thesame size as would be obtained usingvcov.hurdle orvcov.zeroinfl with itsmodel argument set to("full", "count", "zero") in respective correspondence withmode of("mean", "count", "zero"). Ifvcov. is a function, it must support themodelargument.

Back to quick reference

Group E – GEE models

These models all have more than one covariance estimate available,and it may be selected by supplying a string as thevcov.method argument. It is partially matched with theavailable choices shown in the quick reference. Ingeeseandgeeglm, the aliases"robust" (for"vbeta") and"naive" (for"vbeta.naiv" are also accepted.

If a matrix or function is supplied asvcov.method, itis interpreted as avcov. specification as described for... in the documentation forref_grid.

Group G – Generalized linear models and relatives

Most models in this group receive only standard support as inGroup A, but typically the tests and confidence intervalsare asymptotic. Thus thedf column for tabular results willbeInf.

Some objects in this group may require that the original or referencedataset be provided when callingref_grid() oremmeans().

Forcoxph objects, the estimates we obtain arecomparable to runningpredict.coxph() withreference = "zero"; that is, no covariate centering isdone. The user may useat to specify adjusted covariatevalues. For example, the default reference grid sets each covariate toits mean, so estimates comparable topredict.coxph(..., reference = "sample") would be obtainedby specifyingat = list(x = 0).

In the case ofmgcv::gam objects, there are optionalfreq andunconditional arguments as isdetailed in the documentation formgcv::vcov.gam(). Bothdefault toFALSE. The value ofunconditionalmatters only iffreq = FALSE andobject$Vc isnon-null.

Formgcv::gamm objects,emmeans() resultsare based on theobject$gam part. Unfortunately, that ismissing itscall component, so the user must supply it inthecall argument (e.g.,call = quote(gamm(y ~ s(x), data = dat))) or give thedataset in thedata argument. Alternatively (andrecommended), you may first setobject$gam$call to thequoted call ahead of time. Thewhat arguments are used toselect which model formula to use:"location", "scale"apply togaulss andgevlss families,"shape" applies only togevlss, and"rate", "prob.gt.0" apply toziplss.

Withgam::Gam objects, standard errors are estimatedusing a bootstrap method when there are any smoothers involved.Accordingly, there is an optionalnboot argument that setsthe number of bootstrap replications used to estimate the variances andcovariances of the smoothing portions of the model. Generally, it isbetter to use models fitted viamgcv::gam() rather thangam::gam().

Back to quick reference

Group H –gamlss models

Thewhat argument has possible values of"mu" (default),"sigma","nu", or"tau" depending on which part of the model you want resultsfor. Currently, there is no support when the selected part of the modelcontains a smoothing method likepb().

Group I – Multiple models (via imputation or averaging)

These objects are the results of fitting several models withdifferent predictor subsets or imputed values. Thebhat andV slots are obtained via averaging and, in the case ofmultiple imputation, adding a multiple of the between-imputationcovariance per Rubin’s rules, along with an associateddegrees-of-freedom adjustment (Barnard & Rubin 1999). In the case ofmira models with model classes not supported byemmeans,GitHub issue#446 includes a functionpool_estimates_for_qdrg() thatmay be useful for obtaining results viaqdrg(). Anotheruseful link may be a page that showshowto pool severalemmeans results; that is, instead ofpooling the models and then runningemmeans(), we do justthe reverse.

Support forMuMIn::averaging objects may be somewhatdodgy, as it is not clear that all supported model classes will work.The objectmust have a"modelList" attribute(obtained by constructing the object explicitly from a model list or byincludingfit = TRUE in the call). And each model should befitted withdata as anamed argument inthe call; or else provide adata argument in the call toemmeans() orref_grid(). Only ``full’’averaging is supported; conditional averaging can result innon-positive-definite covariance matrices, and so cannot be considered.No estimability checking is done at present (not clear what we even meanby it).

Also, be aware that support foraveraging objects doesnot pay attention to the class of the models being averaged(nor anyemmeans options associated with that class, suchas alternative modes or d.f. methods), and that you can only obtaindirect results of linear predictions or back-transformations thereof(type = "response"). It doesnot take apartmultivariate models, nor multiple-intercept models (e.g. ordinal ones).(But keep reading…)

Finally, note that special care is needed with models having multiplecomponents (e.g. hurdle models), where there is essentially more thanone set of coefficients. We can handle only one at a time. Asubset argument must be provided to specify whichcoefficients to use; that can be anamed vector of integers(where the names are the names of the actual model terms), or thespecial character values"prefix:pfx","pfx",orwrap:wrp which picks out all coefficients whose namesare prefixed bypfx (e.g.,pfxtrtB) or wrappedbywrp (e.g.,wrp(trtB)). Aformula option is also available for specifying anappropriate formula other thanobject$formula when that isnot suitable.

In many cases (especially with multivariate or ordinal models), youare better off making a copy of one of the averaged models that has allrequired terms, and hacking that object by replacing the coefficients bycoef(object, full = TRUE) and the covariance matrix byvcov(object, full = TRUE). You need to be careful to matchthe names of the coefficients correctly, and to use the same indexes topermute the rows and columns of the covariance matrix. You then use theemmeans support, including any options available, for theclass of the model object, rather than for classaveraging.An example is given towards the end ofIssue 442.

Group K –gls andlme models

ThesigmaAdjust argument is a logical value thatdefaults toTRUE. It is comparable to theadjustSigma option innlme::summary.lme (thename-mangling is to avoid conflicts with the often-usedadjust argument), and determines whether or not adegrees-of-freedom adjustment is performed with models fitted using theML method.

The optionalmode argument affects the degrees offreedom. Themode = "satterthwaite" option determinesdegrees of freedom via the Satterthwaite method: Ifs^2 isthe estimate of some variance, then its Satterthwaite d.f. is2*s^4 / Var(s^2). In case our numerical methods for thisfail, we also offermode = "appx-satterthwaite" as abackup, by which quantities related toVar(s^2) areobtained by randomly perturbing the response values. Currently, only"appx-satterthwaite" is available forlmeobjects, and it is used if"satterthwaite" is requested.Becauseappx-satterthwaite is simulation-based, results mayvary if the same analysis is repeated. Anextra.iterargument may be added to request additional simulation runs (at[possibly considerable] cost of repeating the model-fitting that manymore times). (Note: Previously,"appx-satterthwaite" wastermed"boot-satterthwaite"; this is still supported forbackward compatibility. The “boot” was abandoned because it is really anapproximation method, not a bootstrap method in the sense as manystatistical methods.)

An alternative method is"df.error" (forgls) and"containment" (forlme).df.error is just the error degrees of freedom for themodel, minus the number of extra random effects estimated; it generallyover-estimates the degrees of freedom. Theasymptotic modesimply sets the degrees of freedom to infinity."containment" mode (forlme models) determinesthe degrees of freedom for the coarsest grouping involved in thecontrast or linear function involved, so it tends to under-estimate thedegrees of freedom. The default ismode = "auto", whichuses Satterthwaite if there are estimated random effects and thenon-Satterthwaite option otherwise.

User reports indicate that models with special terms likepoly() are not adequately supported bygls inthat the needed basis is not recoverable from itstermscomponent. This is not a problem withlme.

Theextra.iter argument is ignored unless the d.f.method is (or defaults to)appx-satterthwaite.

Back to quick reference

Group L –lmerMod models

There is an optionallmer.df argument that defaults toget_EMM_option("lmer.df") (which in turn defaults to"kenward-roger"). The possible values are"kenward-roger","satterthwaite", and"asymptotic" (these are partially matched andcase-insensitive). With"kenward-roger", d.f. are obtainedusing code from thepbkrtest package, if installed.With"satterthwaite", d.f. are obtained using code from thelmerTest package, if installed. With"asymptotic", or if the needed package is not installed,d.f. are set toInf. (For backward compatibility, the usermay specifymode in lieu oflmer.df.)

A by-product of the Kenward-Roger method is that the covariancematrix is adjusted usingpbkrtest::vcovAdj(). This canrequire considerable computation; so to avoid that overhead, the usershould opt for the Satterthwaite or asymptotic method; or, for backwardcompatibility, may disable the use ofpbkrtest viaemm_options(disable.pbkrtest = TRUE) (this does not disablethepbkrtest package entirely, just its use inemmeans). The computation time required depends roughlyon the number of observations,N, in the design matrix (becausea major part of the computation involves inverting anN xN matrix). Thus,pbkrtest is automaticallydisabled ifN exceeds the value ofget_emm_option("pbkrtest.limit"), for which the factorydefault is 3000. (The user may also specifypbkrtest.limitordisable.pbkrtest as an argument in the call toemmeans() orref_grid())

Similarly to the above, thedisable.lmerTest andlmerTest.limit options or arguments affect whetherSatterthwaite methods can be implemented.

Thedf argument may be used to specify some otherdegrees of freedom. Note that ifdf andmethod = "kenward-roger" are both specified, the covariancematrix is adjusted but the K-R degrees of freedom are not used.

Finally, note that a user-specified covariance matrix (via thevcov. argument) will also disable the Kenward-Roger method;in that case, the Satterthwaite method is used in place ofKenward-Roger.

Back to quick reference

Group M – Multivariate models

When there is a multivariate response, the different responses aretreated as if they were levels of a factor – namedrep.measby default. Themult.name argument may be used to changethis name. Themult.levs argument may specify a named listof one or more sets of levels. If this has more than one element, thenthe multivariate levels are expressed as combinations of the namedfactor levels via the functionbase::expand.grid.

N - Multinomial responses

The reference grid includes a pseudo-factor with the same name andlevels as the multinomial response. (If the response is an expression,the name of that pseudo-factor will be the first name in the expression;e.g., if the model formula iscbind(col1, col2) ~ trt, thegrid factors will becbind andtrt and thelevels ofcbind will be1 and2.)(You can change the assigned name for the multinomial response via themult.resp argument.)

There is an optionalmode argument which should match"prob" or"latent". Withmode = "prob", the reference-grid predictions consist ofthe estimated multinomial probabilities – and this implies a re-griddingso no link functions are passed on. The"latent" modereturns the linear predictor, re-centered so that it averages to zeroover the levels of the response variable (similar to sum-to-zerocontrasts). Thus each latent variable can be regarded as the logprobability at that level minus the average log probability over alllevels.

There are two optional arguments:mode andrescale (which defaults toc(0, 1)).

Please note that, because the probabilities sum to 1 (and the latentvalues sum to 0) over the multivariate-response levels, all sensibleresults fromemmeans() must involve that response as one ofthe factors. For example, ifresp is a response withk levels,emmeans(model, ~ resp | trt) will yieldthe estimated multinomial distribution for eachtrt; butemmeans(model, ~ trt) will just yield the averageprobability of 1/k for eachtrt.

Back to quick reference

Group O - Ordinal responses

The reference grid for ordinal models will include all variables thatappear in the main model as well as those in thescale ornominal models (if provided). There are two optionalarguments:mode (a character string) andrescale (which defaults toc(0, 1)).mode should match one of"latent" (thedefault),"linear.predictor","cum.prob","exc.prob","prob","mean.class",or"scale" – see the quick reference and note which aresupported. With the exception of"linear.predictor", all ofthese modes do an implied regrid.

Withmode = "latent", the reference-grid predictions aremade on the scale of the latent variable implied by the model. The scaleand location of this latent variable are arbitrary, and may be alteredviarescale. The predictions are multiplied byrescale[2], then added torescale[1]. Keep inmind that the scaling is related to the link function used in the model;for example, changing from a probit link to a logistic link will inflatethe latent values by around\(\pi/\sqrt{3}\), all other things beingequal.rescale has no effect for other values ofmode. Even though the latent means comprise a re-scaling ofthe linear predictor, we regard this as a re-gridding, and thecumulative link function is not included in the reference grid.

Withmode = "linear.predictor",mode = "cum.prob", andmode = "exc.prob", theboundaries between categories (i.e., thresholds) in the ordinal responseare included in the reference grid as a pseudo-factor namedcut. The reference-grid predictions are then of thecumulative probabilities at each threshold (formode = "cum.prob"), exceedance probabilities (one minuscumulative probabilities, formode = "exc.prob"), or thelink function thereof (formode = "linear.predictor").

Withmode = "prob", a pseudo-factor with the same nameas the model’s response variable is created, and the grid predictionsare of the probabilities of each class of the ordinal response. With"mean.class", the returned results are means of the ordinalresponse, interpreted as a numeric value from 1 to the number ofclasses, using the"prob" results as the estimatedprobability distribution for each case.

Withmode = "scale", and the fitted object incorporatesa scale model, EMMs are obtained for the factors in the scale model(with a log response) instead of the response model. The grid isconstructed using only the factors in the scale model.

Any grid point that is non-estimable by either the location or thescale model (if present) is set toNA, and any EMMsinvolving such a grid point will also be non-estimable. A consequence ofthis is that if there is a rank-deficientscale model, thenall latent responses become non-estimable because thepredictions are made using the average log-scale estimate.

rms models have an additionalmode. Withmode = "middle" (this is the default), the middle interceptis used, comparable to the default forrms::Predict(). Thisis quite similar in concept tomode = "latent", where allintercepts are averaged together.

Back to quick reference

Group P – Other packages

Models in this group have theiremmeans supportprovided by another package (usually the one that implements themodel-fitting procedure). Users should refer to the packagedocumentation for details onemmeans support. In somecases, a package’s models may have been supported here inemmeans; if so, the other package’s support overridesit.

Group Q – Quantile regression

The elements oftau are included in the reference gridas a pseudo-factor namedtau. In these models, thecovariance matrix is obtained via the model’ssummary()method withcovariance = TRUE. The user may specify one ormore of the other arguments forsummary (e.g.,se = "boot") or to be passed to the...argument.

A caveat is that when there is more than onetau value,we do not have estimates of the covariances between regressioncoefficients associated with differenttaus. Thus, acontrast involving differenttaus can be estimated but itsSE will beNA. Also, due toNAs in thecovariance matrix, the"mvt" adjustment is unavailable.

Note: Older versions of thisrq andrqs supportrequiredtau; now it isoptional. Onlytau values included inobject$tau are allowed; others are silently ignored. It ismore efficient (and less memory-greedy) to specifytau thanto subset them using theat argument toref_grid.

Group S – Sampling (MCMC) methods

Models fitted using MCMC methods contain a sample from the posteriordistribution of fixed-effect coefficients. In some cases (e.g., resultsofMCMCpack::MCMCregress() andMCMCpack::MCMCpoisson()), the object may include a"call" attribute thatemmeans() can use toreconstruct the data and obtain a basis for the EMMs. If not, aformula anddata argument are provided thatmay help produce the right results. In addition, thecontrasts specifications are not necessarily recoverablefrom the object, so the system default must match what was actually usedin fitting the model.

Thesummary.emmGrid() method provides credibilityintervals (HPD intervals) of the results, and ignores thefrequentist-oriented arguments (infer,adjust,etc.) Anas.mcmc() method is provided that creates anmcmc object that can be summarized or plotted using thecoda package (or others that support those objects). Itprovides a posterior sample of EMMs, or contrasts thereof, for the givenreference grid, based on the posterior sample of the fixed effects fromthe model object.

InMCMCglmm objects, thedata argument isrequired; however, if you save it as a member of the model object (e.g.,object$data = quote(mydata)), that removes the necessity ofspecifying it in each call. The special keywordtrait isused in some models. When the response is multivariate and numeric,trait is generated automatically as a factor in thereference grid, and the argumentsmult.levels can be usedto name its levels. In other models such as a multinomial model, use themode argument to specify the type of model, andtrait = <factor name> to specify the name of the datacolumn that contains the levels of the factor response.

Thebrms package version 2.13 and later, has its ownemmeans support. Refer to the documentation in thatpackage.

Back to quick reference

Group V –aovList objects (also used withafex_aov objects)

Support for these objects is limited. To avoid strong biases in thepredictions, it is strongly recommended that when fitting the model, thecontrasts attribute of all factors should be of a type thatsums to zero – for example,"contr.sum","contr.poly", or"contr.helmert" butnot"contr.treatment". If that is found not to bethe case, the model is re-fitted using sum-to-zero contrasts (thusrequiring additional computation). Doing so doesnot remove allbias in the EMMs unless the design is perfectly balanced, and anannotation is added to warn of that. This bias cancels out when doingcomparisons and contrasts.

Only intra-block estimates of covariances are used. That is, if afactor appears in more than one error stratum, only the covariancestructure from its lowest stratum is used in estimating standard errors.Degrees of freedom are obtained using the Satterthwaite method. Ingeneral,aovList support is best with balanced designs,with due caution in the use of contrasts. If avcov.argument is supplied, it must yield a single covariance matrix for theunique fixed effects (not a set of them for each error stratum). In thatcase, the degrees of freedom are set toNA.

Back to quick reference

Index of all vignette topics