Without interaction
at a special value
# Compute least-squares means for specified factors or factor combinations in a linear model
require(lsmeans)
## Loading required package: lsmeans
## Loading required package: emmeans
## The 'lsmeans' package is now basically a front end for 'emmeans'.
## Users are encouraged to switch the rest of the way.
## See help('transition') for more information, including how to
## convert old 'lsmeans' objects and scripts to work with 'emmeans'.
# Use dataset 'fiber'
### model 1 without interaction, Covariance example
fiber.lm <- lm(strength ~ diameter + machine, data = fiber)
# adjusted means and comparison
lsmeans(fiber.lm, list (pairwise ~ machine ) )
## $`lsmeans of machine`
## machine lsmean SE df lower.CL upper.CL
## A 40.4 0.724 11 38.8 42.0
## B 41.4 0.744 11 39.8 43.1
## C 38.8 0.788 11 37.1 40.5
##
## Confidence level used: 0.95
##
## $`pairwise differences of machine`
## 1 estimate SE df t.ratio p.value
## A - B -1.04 1.01 11 -1.024 0.5781
## A - C 1.58 1.11 11 1.431 0.3596
## B - C 2.62 1.15 11 2.283 0.1005
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# adjusted means and contrast, treating machine C as control
lsmeans (fiber.lm, "machine", contr = "trt.vs.ctrlk")
## $lsmeans
## machine lsmean SE df lower.CL upper.CL
## A 40.4 0.724 11 38.8 42.0
## B 41.4 0.744 11 39.8 43.1
## C 38.8 0.788 11 37.1 40.5
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## A - C 1.58 1.11 11 1.431 0.3062
## B - C 2.62 1.15 11 2.283 0.0796
##
## P value adjustment: dunnettx method for 2 tests
# adjusted means and comparison, after setting "diameter" as a specific value
lsmeans(fiber.lm, list (pairwise ~ machine ), at = list(diameter =20) )
## $`lsmeans of machine`
## machine lsmean SE df lower.CL upper.CL
## A 36.4 0.928 11 34.4 38.5
## B 37.5 0.988 11 35.3 39.7
## C 34.9 0.726 11 33.3 36.5
##
## Confidence level used: 0.95
##
## $`pairwise differences of machine`
## 1 estimate SE df t.ratio p.value
## A - B -1.04 1.01 11 -1.024 0.5781
## A - C 1.58 1.11 11 1.431 0.3596
## B - C 2.62 1.15 11 2.283 0.1005
##
## P value adjustment: tukey method for comparing a family of 3 estimates
With interaction
### model 2 with interaction
fiber.lm2 <- lm(strength ~ diameter*machine, data=fiber)
# means when diameter is its mean
(lsmeans(fiber.lm2, ~ machine | diameter))
## diameter = 24.1:
## machine lsmean SE df lower.CL upper.CL
## A 40.2 0.777 9 38.5 42.0
## B 41.6 0.858 9 39.7 43.5
## C 38.5 0.966 9 36.3 40.7
##
## Confidence level used: 0.95
# means when setting diameter as its mean
lsmeans(fiber.lm2, list (pairwise ~ machine ), at = list(diameter =24.1) )
## NOTE: Results may be misleading due to involvement in interactions
## $`lsmeans of machine`
## machine lsmean SE df lower.CL upper.CL
## A 40.2 0.779 9 38.4 41.9
## B 41.6 0.861 9 39.6 43.5
## C 38.5 0.962 9 36.3 40.7
##
## Confidence level used: 0.95
##
## $`pairwise differences of machine`
## 1 estimate SE df t.ratio p.value
## A - B -1.39 1.16 9 -1.194 0.4858
## A - C 1.68 1.24 9 1.357 0.4017
## B - C 3.07 1.29 9 2.374 0.0956
##
## P value adjustment: tukey method for comparing a family of 3 estimates
### model 3 for category covariates
warp.lm <- lm(breaks ~ wool * tension, data = warpbreaks)
# means after adjusting the interaction effect
(lsmeans(warp.lm, ~ wool | tension))
## tension = L:
## wool lsmean SE df lower.CL upper.CL
## A 44.6 3.65 48 37.2 51.9
## B 28.2 3.65 48 20.9 35.6
##
## tension = M:
## wool lsmean SE df lower.CL upper.CL
## A 24.0 3.65 48 16.7 31.3
## B 28.8 3.65 48 21.4 36.1
##
## tension = H:
## wool lsmean SE df lower.CL upper.CL
## A 24.6 3.65 48 17.2 31.9
## B 18.8 3.65 48 11.4 26.1
##
## Confidence level used: 0.95
# set tension ="M"
lsmeans(warp.lm, list (pairwise ~ wool ), at = list(tension ="M"))
## NOTE: Results may be misleading due to involvement in interactions
## $`lsmeans of wool`
## wool lsmean SE df lower.CL upper.CL
## A 24.0 3.65 48 16.7 31.3
## B 28.8 3.65 48 21.4 36.1
##
## Confidence level used: 0.95
##
## $`pairwise differences of wool`
## 1 estimate SE df t.ratio p.value
## A - B -4.78 5.16 48 -0.926 0.3589
lsmeans(warp.lm, list (pairwise ~ wool ), at = list(tension ="H"))
## NOTE: Results may be misleading due to involvement in interactions
## $`lsmeans of wool`
## wool lsmean SE df lower.CL upper.CL
## A 24.6 3.65 48 17.2 31.9
## B 18.8 3.65 48 11.4 26.1
##
## Confidence level used: 0.95
##
## $`pairwise differences of wool`
## 1 estimate SE df t.ratio p.value
## A - B 5.78 5.16 48 1.120 0.2682
Compare slopes and test
### Trends, obtain slopes
fiber.lst=( lstrends(fiber.lm, "machine", var="diameter") )
fiber.lst
## machine diameter.trend SE df lower.CL upper.CL
## A 0.954 0.114 11 0.703 1.21
## B 0.954 0.114 11 0.703 1.21
## C 0.954 0.114 11 0.703 1.21
##
## Confidence level used: 0.95
# pairwise comparisons thereof slopes
pairs(fiber.lst)
## contrast estimate SE df t.ratio p.value
## A - B 0 0 11 NaN NaN
## A - C 0 0 11 NaN NaN
## B - C 0 0 11 NaN NaN
##
## P value adjustment: tukey method for comparing a family of 3 estimates
fiber.lst2=( lstrends(fiber.lm2, "machine", var="diameter") )
fiber.lst2
## machine diameter.trend SE df lower.CL upper.CL
## A 1.104 0.194 9 0.666 1.54
## B 0.857 0.224 9 0.351 1.36
## C 0.864 0.208 9 0.394 1.33
##
## Confidence level used: 0.95
Least-squares means with interaction for LMM
library(lme4)
## Loading required package: Matrix
##############
m1 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy)
lsmeans(m1, ~ Days , at = list(Days =20) )
## Days lsmean SE df lower.CL upper.CL
## 20 461 15.4 103 430 491
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
lsmeans(m1, ~ Days )
## Days lsmean SE df lower.CL upper.CL
## 4.5 299 9.05 17 279 318
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95