lsmeans: Least-squares means

Description

Compute least-squares means for specified factors or factor combinations in a linear model, and optionally comparisons or contrasts among them.

Usage

lsmeans(object, specs, adjust=c("auto","tukey","sidak","bonferroni","none"),
  conf = .95, at, contr = list(), 
  cov.reduce = function(x, name) mean(x), 
  fac.reduce = function(coefs, lev) apply(coefs, 2, mean), 
  glhargs=NULL, ...)
  
## S3 method for class 'lsm':
print(x, omit=NULL, ...)

Arguments

object

A lm, aov (with no Error component), glm, lme, gls, lmer, or glmer object having at least one fixed factor among the predictors.

specs

A formula, or a list of formulas, specifying the desired families of least-squares means. The right-hand side of each formula specifies the desired factor levels. The optional left-hand side specifies what kind of comparisons or contrasts are desired. F

adjust

Adjustment method for the p values of tests of contrasts. "auto" uses the method returned in the "adjust" attribute of the contrast function; "tukey" computes p values using the Studentized range distribution with

conf

Desired confidence level for intervals. For robustness, you may specify either a fraction or a percentage; i.e., .975 and 97.5 yield the same results.

An optional named list or named vector of covariate values at which predictions are computed (give only one value for each covariate). If no value is found in at for a particular covariate, then cov.reduce is called.

contr

An optional named list. Each entry is itself a list or a data.frame specifying contrast coefficients. If the left-hand side of a formula in specs matches a name in contr, then those contrasts are estimated with the specified leas

cov.reduce

A function with arguments x and name that should return the value to use in prediction for the covariate with name name and values x. By default, the mean is used. If specified, the name arg

fac.reduce

A function of coefs and lev where lev is the level of a factor or factor combination at which a least-squares mean is calculated. The argument coefs is a matrix whose rows correspond to the combinations

glhargs

If this is a list, the object and specified contrasts are passed to the functionglht in the multcomp package, with the contents of glhargs as additional arguments

...

Additional argument(s) passed to the contrast function; see Details.

Object of class "lsm"

omit

Indexes of elements of x that you do not want printed.

Value

An object of class "lsm", which inherits from "list". Each element of the list is either a data.frame or an object of class "glht" (see the documentation for glht). (The latter occur only if glhargs is non-NULL.) Each element summarizes a family of least-squares means or contrasts among them. Each data.frame contains lsmeans or contrast estimates and associated quantities; in addition, there may be a mesg attribute with character string(s) providing information on multiplicity adjustments and such. The "lsm" class has only one method, print, which displays data.frame elements as-is along with any mesg attributes; and the summary of any glht elements.

Details

Least-squares means, popularized by SAS, are predictions from a linear model at combinations of specified factors. SAS's documentation describes them as ``predicted population margins---that is, they estimate the marginal means over a balanced population'' (SAS Institute 2012). In generalized linear models, least-squares means are marginal linear predictions that can be transformed back to the response scale via the inverse-link function. Unspecified factors and covariates are handled by summarizing the predictions over those factors and variables. For example, if the fitted model has formula response ~ x1 + x2 + treat where a1 and x2 are numeric and treat is a factor, the least-squares means will be the predicted response for each treatment, at some specified values of x1 and x2. By default, the means of the two covariates will be used, resulting in what ANOVA textbooks oftem call the adjusted means. We may use that at argument to instead make predictions at other values of x1 and x2. Now consider the model response ~ A + B + A:B, where A and B are both factors. If we ask for least-squares means for A, then at each level of A we are faced with a different prediction for each level of B. Blind (and default) use of least-squares means would result in these predictions being averaged together with equal weight, and this may be inappropriate, especially when the interaction effect is strong. Like most statistical calculations, it is possible to use least-squares means inappropriately. The fac.reduce argument at least expands one's options in producing meaningful results in multi-factor situations. One other note concerning covariates: One must be careful with covariates that depend on one another. For example, if a model contains covariates x and xsq where xsq = x^2, the default behavior will make predictions at x = mean(x) and xsq = mean(xsq), which probably isn't a valid combination (we need x = mean(x) and xsq = mean(x)^2). The inconsistency is avoided if the model specifis poly(x,2) (or even x + I(x^2)) instead of x + xsq, because then only x appears as a covariate and everything remains consistent when we substitute its mean value. The built-in contrast methods that can be used in specs formulas are pairwise, revpairwise, poly, trt.vs.ctrl, trt.vs.ctrl1, and trt.vs.ctrlk. They are implemented as functions pairwise.lsmc, etc. having the same names with .lsmc added. Users may write additional .lsmc functions that generate custom families of contrasts. See the documentation for pairwise.lsmc for an example. Degrees of freedom are currently not provided for lme or glme objects, or for mer objects arising from generalized linear models; in those cases, asymptotic results are printed, and this is emphasized by displaying NA for the defrees of freedom. For linear mer objects, degrees of freedom are computed using the Kenward and Roger (1997) method, provided the pbkrtest package is installed (the package is loaded if needed.) Moreover, in that case, the adjusted covariance matrix from the vcovAdj() function in the pbkrtest package is used to calculate standard errors. See Halekoh and H�jsgaard (2012) and the documentation for KRmodcomp for more details. Degrees of freedom are not passed to glht except in the case of lm objects. If the model contains a matrix among its predictors, each column is averaged using the function specified in cov.reduce. There is no provision for matrices in the at argument.

References

Halekoh, U. and H�jsgaard, S. (2012), A Kenward-Roger Approximation and parametric bootsrap methods for tests lin linear mixed models -- the R package pbkrtest, submitted. Kenward, M.G. and Roger, J.H. (1997), Small sample inference for fixed effects from restricted maximum likelihood, Biometrics, 53, 983--997. SAS Institute Inc. (2012), Online documentation; Shared concepts; LSMEANS statement, http://support.sas.com/documentation/cdl/en/statug/63962/HTML/default/viewer.htm#statug_introcom_a0000003362.htm, accessed August 15, 2012.

Examples

Run this code

require(lsmeans)

### Covariance example (from Montgomery Design (7th ed.), p.591)
fiber = data.frame(
  machine = rep(c("A","B","C"), each=5),
  strength = c(36,41,39,42,49, 40,48,39,45,44, 35,37,42,34,32),
  diameter = c(20,25,24,25,32, 22,28,22,30,28, 21,23,26,21,15))
fiber.lm = lm(strength ~ diameter + machine, data = fiber)

# adjusted means and comparisons, treating machine C as control
lsmeans (fiber.lm, trt.vs.ctrlk ~ machine)


### Factorial experiment
warp.lm = lm(breaks ~ wool * tension, data = warpbreaks)
#-- We only need to see the wool*tension means listed once ...
print(lsmeans (warp.lm,  list(pairwise ~ wool | tension,  poly ~ tension | wool)),
    omit=3)


### Unbalanced split-plot example ###
#-- The imbalance biases the variance estimates somewhat
require(nlme)
Oats.lme = lme(yield ~ factor(nitro) + Variety, random = ~1 | Block/Variety, 
    subset = -c(1,2,3,5,8,13,21,34,55), data=Oats)
lsmeans(Oats.lme, list(poly ~ nitro, pairwise ~ Variety))

# Compare with lmer result (lsmeans provides df, adjusted SEs)
require(lme4)
Oats.lmer = lmer(yield ~ factor(nitro) + Variety + (1 | Block/Variety), 
    subset = -c(1,2,3,5,8,13,21,34,55), data=Oats)
#-- require(pbkrtest) #-- (loaded as needed by lsmeans)
lsmeans(Oats.lmer, list(poly ~ nitro, pairwise ~ Variety))

# Use glht (multcomp) to do comparisons (but does not use adjusted vcov)
#-- require(multcomp) #-- (loaded as needed by lsmeans)
lsmeans(Oats.lmer, pairwise ~ Variety, glhargs=list(df=9.5))

# Custom contrasts
lsmeans(Oats.lmer, my.own ~ Variety, 
  contr = list(my.own = list(G.vs.M = c(1,-1,0), GM.vs.V = c(.5,.5,-1))))

Run the code above in your browser using DataLab