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 1.1509 0.679 1.6945 8 0.129 -0.415 2.717
#> 2 2 1 -0.3026 0.679 -0.4455 8 0.668 -1.869 1.264
#> 3 3 1 0.0118 0.679 0.0174 8 0.987 -1.554 1.578
#> 4 1 2 0.2756 0.679 0.4058 8 0.696 -1.291 1.842
#> 5 2 2 -0.8480 0.679 -1.2486 8 0.247 -2.414 0.718
#> 6 3 2 0.8651 0.679 1.2738 8 0.239 -0.701 2.431
#> 7 1 3 -0.1589 0.679 -0.2340 8 0.821 -1.725 1.407
#> 8 2 3 0.7936 0.679 1.1685 8 0.276 -0.773 2.360
#> 9 3 3 -0.4821 0.679 -0.7099 8 0.498 -2.048 1.084
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.4225 0.391 1.0813 10 0.305 -0.448 1.29
#> 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 -0.3026 0.677 -0.4470 10 0.664 -1.811 1.21
#> 6 3 2 0.0118 0.677 0.0175 10 0.986 -1.496 1.52
#> 7 1 3 NA NA NA NA NA NA NA
#> 8 2 3 -0.8480 0.677 -1.2530 10 0.239 -2.356 0.66
#> 9 3 3 0.8651 0.677 1.2783 10 0.230 -0.643 2.37
#> 10 1 4 NA NA NA NA NA NA NA
#> 11 2 4 0.7936 0.677 1.1726 10 0.268 -0.714 2.30
#> 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.423 0.391 1.08 10 0.305 -0.448 1.29
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