LS-means (least squares means, also known as population means and as marginal means) for a range of model types.
LSmeans(object, effect = NULL, at = NULL, level = 0.95, ...)
# Default S3 method
LSmeans(object, effect = NULL, at = NULL, level = 0.95, ...)
# S3 method for class 'lmerMod'
LSmeans(object, effect = NULL, at = NULL, level = 0.95, adjust.df = TRUE, ...)
popMeans(object, effect = NULL, at = NULL, level = 0.95, ...)
# Default S3 method
popMeans(object, effect = NULL, at = NULL, level = 0.95, ...)
# S3 method for class 'lmerMod'
popMeans(object, effect = NULL, at = NULL, level = 0.95, adjust.df = TRUE, ...)Model object
A vector of variables. For each configuration of these the estimate will be calculated.
A list of values of covariates (including levels of some factors) to be used in the calculations
The level of the (asymptotic) confidence interval.
Additional arguments; currently not used.
Should denominator degrees of freedom be adjusted?
A dataframe with results from computing the contrasts.
There are restrictions on the formulas allowed in the model object.
For example having y ~ log(x) will cause an error. Instead one
must define the variable logx = log(x) and do y ~ logx.
LSmeans and popMeans are synonymous.
Notice that LSmeans and LE_matrix
fails if the model formula contains an offset (as one would
have in connection with e.g. Poisson regression.
## Two way anova:
data(warpbreaks)
m0 <- lm(breaks ~ wool + tension, data=warpbreaks)
m1 <- lm(breaks ~ wool * tension, data=warpbreaks)
LSmeans(m0)
#> estimate std.error statistic df p.value lwr upr
#> 1 28.1 1.58 17.8 50 2.74e-23 25 31.3
LSmeans(m1)
#> estimate std.error statistic df p.value lwr upr
#> 1 28.1 1.49 18.9 48 6.98e-24 25.2 31.1
## same as:
K <- LE_matrix(m0);K
#> (Intercept) woolB tensionM tensionH
#> [1,] 1 0.5 0.3333333 0.3333333
linest(m0, K)
#> estimate std.error statistic df p.value lwr upr
#> 1 28.1 1.58 17.8 50 2.74e-23 25 31.3
K <- LE_matrix(m1);K
#> (Intercept) woolB tensionM tensionH woolB:tensionM woolB:tensionH
#> [1,] 1 0.5 0.3333333 0.3333333 0.1666667 0.1666667
linest(m1, K)
#> estimate std.error statistic df p.value lwr upr
#> 1 28.1 1.49 18.9 48 6.98e-24 25.2 31.1
LE_matrix(m0, effect="wool")
#> List of 2
#> $ new.fact.lev:List of 1
#> ..$ wool: chr [1:2] "A" "B"
#> $ grid.data :'data.frame': 2 obs. of 1 variable:
#> ..$ wool: chr [1:2] "A" "B"
#> (Intercept) woolB tensionM tensionH
#> [1,] 1 0 0.3333333 0.3333333
#> [2,] 1 1 0.3333333 0.3333333
LSmeans(m0, effect="wool")
#> List of 2
#> $ new.fact.lev:List of 1
#> ..$ wool: chr [1:2] "A" "B"
#> $ grid.data :'data.frame': 2 obs. of 1 variable:
#> ..$ wool: chr [1:2] "A" "B"
#> wool estimate std.error statistic df p.value lwr upr
#> 1 A 31.0 2.24 13.9 50 8.81e-19 26.5 35.5
#> 2 B 25.3 2.24 11.3 50 2.26e-15 20.8 29.7
LE_matrix(m1, effect="wool")
#> List of 2
#> $ new.fact.lev:List of 1
#> ..$ wool: chr [1:2] "A" "B"
#> $ grid.data :'data.frame': 2 obs. of 1 variable:
#> ..$ wool: chr [1:2] "A" "B"
#> (Intercept) woolB tensionM tensionH woolB:tensionM woolB:tensionH
#> [1,] 1 0 0.3333333 0.3333333 0.0000000 0.0000000
#> [2,] 1 1 0.3333333 0.3333333 0.3333333 0.3333333
LSmeans(m1, effect="wool")
#> List of 2
#> $ new.fact.lev:List of 1
#> ..$ wool: chr [1:2] "A" "B"
#> $ grid.data :'data.frame': 2 obs. of 1 variable:
#> ..$ wool: chr [1:2] "A" "B"
#> wool estimate std.error statistic df p.value lwr upr
#> 1 A 31.0 2.11 14.7 48 1.91e-19 26.8 35.3
#> 2 B 25.3 2.11 12.0 48 4.71e-16 21.0 29.5
LE_matrix(m0, effect=c("wool", "tension"))
#> List of 2
#> $ new.fact.lev:List of 2
#> ..$ wool : chr [1:2] "A" "B"
#> ..$ tension: chr [1:3] "L" "M" "H"
#> $ grid.data :'data.frame': 6 obs. of 2 variables:
#> ..$ wool : chr [1:6] "A" "B" "A" "B" ...
#> ..$ tension: chr [1:6] "L" "L" "M" "M" ...
#> (Intercept) woolB tensionM tensionH
#> [1,] 1 0 0 0
#> [2,] 1 1 0 0
#> [3,] 1 0 1 0
#> [4,] 1 1 1 0
#> [5,] 1 0 0 1
#> [6,] 1 1 0 1
LSmeans(m0, effect=c("wool", "tension"))
#> List of 2
#> $ new.fact.lev:List of 2
#> ..$ wool : chr [1:2] "A" "B"
#> ..$ tension: chr [1:3] "L" "M" "H"
#> $ grid.data :'data.frame': 6 obs. of 2 variables:
#> ..$ wool : chr [1:6] "A" "B" "A" "B" ...
#> ..$ tension: chr [1:6] "L" "L" "M" "M" ...
#> wool tension estimate std.error statistic df p.value lwr upr
#> 1 A L 39.3 3.16 12.42 50 6.68e-17 32.9 45.6
#> 2 B L 33.5 3.16 10.60 50 2.22e-14 27.1 39.9
#> 3 A M 29.3 3.16 9.26 50 2.00e-12 22.9 35.6
#> 4 B M 23.5 3.16 7.43 50 1.27e-09 17.1 29.9
#> 5 A H 24.6 3.16 7.77 50 3.83e-10 18.2 30.9
#> 6 B H 18.8 3.16 5.94 50 2.72e-07 12.4 25.1
LE_matrix(m1, effect=c("wool", "tension"))
#> List of 2
#> $ new.fact.lev:List of 2
#> ..$ wool : chr [1:2] "A" "B"
#> ..$ tension: chr [1:3] "L" "M" "H"
#> $ grid.data :'data.frame': 6 obs. of 2 variables:
#> ..$ wool : chr [1:6] "A" "B" "A" "B" ...
#> ..$ tension: chr [1:6] "L" "L" "M" "M" ...
#> (Intercept) woolB tensionM tensionH woolB:tensionM woolB:tensionH
#> [1,] 1 0 0 0 0 0
#> [2,] 1 1 0 0 0 0
#> [3,] 1 0 1 0 0 0
#> [4,] 1 1 1 0 1 0
#> [5,] 1 0 0 1 0 0
#> [6,] 1 1 0 1 0 1
LSmeans(m1, effect=c("wool", "tension"))
#> List of 2
#> $ new.fact.lev:List of 2
#> ..$ wool : chr [1:2] "A" "B"
#> ..$ tension: chr [1:3] "L" "M" "H"
#> $ grid.data :'data.frame': 6 obs. of 2 variables:
#> ..$ wool : chr [1:6] "A" "B" "A" "B" ...
#> ..$ tension: chr [1:6] "L" "L" "M" "M" ...
#> wool tension estimate std.error statistic df p.value lwr upr
#> 1 A L 44.6 3.65 12.22 48 2.43e-16 37.2 51.9
#> 2 B L 28.2 3.65 7.74 48 5.47e-10 20.9 35.6
#> 3 A M 24.0 3.65 6.58 48 3.23e-08 16.7 31.3
#> 4 B M 28.8 3.65 7.89 48 3.22e-10 21.4 36.1
#> 5 A H 24.6 3.65 6.73 48 1.88e-08 17.2 31.9
#> 6 B H 18.8 3.65 5.15 48 4.84e-06 11.4 26.1
## Regression; two parallel regression lines:
data(Puromycin)
m0 <- lm(rate ~ state + log(conc), data=Puromycin)
## Can not use LSmeans / LE_matrix here because of
## the log-transformation. Instead we must do:
Puromycin$lconc <- log( Puromycin$conc )
m1 <- lm(rate ~ state + lconc, data=Puromycin)
LE_matrix(m1)
#> (Intercept) stateuntreated lconc
#> [1,] 1 0.5 -1.905247
LSmeans(m1)
#> estimate std.error statistic df p.value lwr upr
#> 1 126 2.37 53.2 20 5.12e-23 121 131
LE_matrix(m1, effect="state")
#> List of 2
#> $ new.fact.lev:List of 1
#> ..$ state: chr [1:2] "treated" "untreated"
#> $ grid.data :'data.frame': 2 obs. of 1 variable:
#> ..$ state: chr [1:2] "treated" "untreated"
#> (Intercept) stateuntreated lconc
#> [1,] 1 0 -1.905247
#> [2,] 1 1 -1.905247
LSmeans(m1, effect="state")
#> List of 2
#> $ new.fact.lev:List of 1
#> ..$ state: chr [1:2] "treated" "untreated"
#> $ grid.data :'data.frame': 2 obs. of 1 variable:
#> ..$ state: chr [1:2] "treated" "untreated"
#> state lconc estimate std.error statistic df p.value lwr upr
#> 1 treated -1.91 139 3.29 42.3 20 4.94e-21 132 146
#> 2 untreated -1.91 114 3.43 33.1 20 6.04e-19 107 121
LE_matrix(m1, effect="state", at=list(lconc=3))
#> List of 2
#> $ new.fact.lev:List of 1
#> ..$ state: chr [1:2] "treated" "untreated"
#> $ grid.data :'data.frame': 2 obs. of 1 variable:
#> ..$ state: chr [1:2] "treated" "untreated"
#> (Intercept) stateuntreated lconc
#> [1,] 1 0 3
#> [2,] 1 1 3
LSmeans(m1, effect="state", at=list(lconc=3))
#> List of 2
#> $ new.fact.lev:List of 1
#> ..$ state: chr [1:2] "treated" "untreated"
#> $ grid.data :'data.frame': 2 obs. of 1 variable:
#> ..$ state: chr [1:2] "treated" "untreated"
#> state lconc estimate std.error statistic df p.value lwr upr
#> 1 treated 3 299 9.35 31.9 20 1.23e-18 279 318
#> 2 untreated 3 273 9.70 28.2 20 1.41e-17 253 294
## Non estimable contrasts
## ## Make balanced dataset
dat.bal <- expand.grid(list(AA=factor(1:2), BB=factor(1:3),
CC=factor(1:3)))
dat.bal$y <- rnorm(nrow(dat.bal))
## ## Make unbalanced dataset
# 'BB' is nested within 'CC' so BB=1 is only found when CC=1
# and BB=2,3 are found in each CC=2,3,4
dat.nst <- dat.bal
dat.nst$CC <-factor(c(1, 1, 2, 2, 2, 2, 1, 1, 3, 3,
3, 3, 1, 1, 4, 4, 4, 4))
mod.bal <- lm(y ~ AA + BB * CC, data=dat.bal)
mod.nst <- lm(y ~ AA + BB : CC, data=dat.nst)
LSmeans(mod.bal, effect=c("BB", "CC"))
#> List of 2
#> $ new.fact.lev:List of 2
#> ..$ BB: chr [1:3] "1" "2" "3"
#> ..$ CC: chr [1:3] "1" "2" "3"
#> $ grid.data :'data.frame': 9 obs. of 2 variables:
#> ..$ BB: chr [1:9] "1" "2" "3" "1" ...
#> ..$ CC: chr [1:9] "1" "1" "1" "2" ...
#> BB CC estimate std.error statistic df p.value lwr upr
#> 1 1 1 0.7356 0.66 1.114 8 0.2976 -0.787 2.2583
#> 2 2 1 -1.5583 0.66 -2.360 8 0.0460 -3.081 -0.0357
#> 3 3 1 -0.6908 0.66 -1.046 8 0.3260 -2.213 0.8318
#> 4 1 2 -0.8061 0.66 -1.221 8 0.2569 -2.329 0.7165
#> 5 2 2 0.1523 0.66 0.231 8 0.8233 -1.370 1.6750
#> 6 3 2 -1.1967 0.66 -1.812 8 0.1075 -2.719 0.3259
#> 7 1 3 1.2685 0.66 1.921 8 0.0909 -0.254 2.7912
#> 8 2 3 -1.4100 0.66 -2.135 8 0.0652 -2.933 0.1126
#> 9 3 3 -0.0876 0.66 -0.133 8 0.8977 -1.610 1.4350
LSmeans(mod.nst, effect=c("BB", "CC"))
#> List of 2
#> $ new.fact.lev:List of 2
#> ..$ BB: chr [1:3] "1" "2" "3"
#> ..$ CC: chr [1:4] "1" "2" "3" "4"
#> $ grid.data :'data.frame': 12 obs. of 2 variables:
#> ..$ BB: chr [1:12] "1" "2" "3" "1" ...
#> ..$ CC: chr [1:12] "1" "1" "1" "2" ...
#> BB CC estimate std.error statistic df p.value lwr upr
#> 1 1 1 0.399 0.440 0.907 10 0.3855 -0.581 1.380
#> 2 2 1 NA NA NA NA NA NA NA
#> 3 3 1 NA NA NA NA NA NA NA
#> 4 1 2 NA NA NA NA NA NA NA
#> 5 2 2 -1.558 0.762 -2.044 10 0.0681 -3.257 0.140
#> 6 3 2 -0.691 0.762 -0.906 10 0.3861 -2.389 1.007
#> 7 1 3 NA NA NA NA NA NA NA
#> 8 2 3 0.152 0.762 0.200 10 0.8456 -1.546 1.851
#> 9 3 3 -1.197 0.762 -1.570 10 0.1475 -2.895 0.502
#> 10 1 4 NA NA NA NA NA NA NA
#> 11 2 4 -1.410 0.762 -1.850 10 0.0941 -3.108 0.288
#> 12 3 4 NA NA NA NA NA NA NA
LSmeans(mod.nst, at=list(BB=1, CC=1))
#> List of 2
#> $ new.fact.lev:List of 2
#> ..$ BB: num 1
#> ..$ CC: num 1
#> $ grid.data :'data.frame': 1 obs. of 2 variables:
#> ..$ BB: chr "1"
#> ..$ CC: chr "1"
#> BB CC estimate std.error statistic df p.value lwr upr
#> 1 1 1 0.399 0.44 0.907 10 0.385 -0.581 1.38
LSmeans(mod.nst, at=list(BB=1, CC=2))
#> List of 2
#> $ new.fact.lev:List of 2
#> ..$ BB: num 1
#> ..$ CC: num 2
#> $ grid.data :'data.frame': 1 obs. of 2 variables:
#> ..$ BB: chr "1"
#> ..$ CC: chr "2"
#> BB CC estimate std.error statistic df p.value lwr upr
#> 1 1 2 NA NA NA NA NA NA NA
## Above: NA's are correct; not an estimable function
if( require( lme4 )){
warp.mm <- lmer(breaks ~ -1 + tension + (1|wool), data=warpbreaks)
LSmeans(warp.mm, effect="tension")
class(warp.mm)
fixef(warp.mm)
coef(summary(warp.mm))
vcov(warp.mm)
if (require(pbkrtest))
vcovAdj(warp.mm)
}
#> Loading required package: lme4
#> Loading required package: Matrix
#> List of 2
#> $ new.fact.lev:List of 1
#> ..$ tension: chr [1:3] "L" "M" "H"
#> $ grid.data :'data.frame': 3 obs. of 1 variable:
#> ..$ tension: chr [1:3] "L" "M" "H"
#> Loading required package: pbkrtest
#> 3 x 3 Matrix of class "dgeMatrix"
#> tensionL tensionM tensionH
#> tensionL 13.344136 5.846482 5.846482
#> tensionM 5.846482 13.344136 5.846482
#> tensionH 5.846482 5.846482 13.344136
LSmeans(warp.mm, effect="tension")
#> List of 2
#> $ new.fact.lev:List of 1
#> ..$ tension: chr [1:3] "L" "M" "H"
#> $ grid.data :'data.frame': 3 obs. of 1 variable:
#> ..$ tension: chr [1:3] "L" "M" "H"
#> tension estimate std.error statistic df p.value lwr upr
#> 1 L 36.4 3.65 9.96 2.54 0.00423 23.47 49.3
#> 2 M 26.4 3.65 7.22 2.54 0.00935 13.47 39.3
#> 3 H 21.7 3.65 5.93 2.54 0.01509 8.75 34.6