library(tidyverse)
library(gtsummary)
library(dfmtbx)
library(SASmixed)
library(emmeans)SAS for Mixed Models
# Use sum contrasts for unordered factors
options(contrasts = c("contr.sum", "contr.poly"))
# Default contrasts
options(contrasts = c("contr.treatment", "contr.poly"))Chapter 1 - Randomized Blocks Designs
Bond data set description
The Bond data set can be obtained from the SASmixed package. The SASmixed package contains the data set used int Litell et. al’s 1996 version of SAS System for Mixed Models. The Bond data set is also a pre-packaged dataset in SAS 9.4 found in sasuser.bond.
# as_tibble converts values to factor
bond <- as_tibble(Bond)The data describe an experiment where the breaking pressure of different metals and ingots are tested. There are a total of 7 different ingot types and three different metal, resulting in 21 observations since each metal type is represented in each ingot type. This a Randomized complete block design.
glimpse(bond)Rows: 21
Columns: 3
$ pressure <dbl> 67.0, 71.9, 72.2, 67.5, 68.8, 66.4, 76.0, 82.6, 74.5, 72.7, 7…
$ Metal <fct> n, i, c, n, i, c, n, i, c, n, i, c, n, i, c, n, i, c, n, i, c
$ Ingot <fct> 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7# Display the first 6 rows of the data
bond %>%
head() %>%
format_gt_tbl()| pressure | Metal | Ingot |
|---|---|---|
| 67.00 | n | 1 |
| 71.90 | i | 1 |
| 72.20 | c | 1 |
| 67.50 | n | 2 |
| 68.80 | i | 2 |
| 66.40 | c | 2 |
bond %>%
group_by(Metal) %>%
summarize_variable(pressure)# A tibble: 3 × 8
Metal n min max median mean sd se
<fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 c 7 66.4 74.5 69 70.2 3.11 1.18
2 i 7 68.8 84.9 74.2 75.9 6.14 2.32
3 n 7 65.8 76 72.7 71.1 4.26 1.61bond %>%
group_by(Ingot) %>%
summarize_variable(pressure)# A tibble: 7 × 8
Ingot n min max median mean sd se
<fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 3 67 72.2 71.9 70.4 2.92 1.69
2 2 3 66.4 68.8 67.5 67.6 1.20 0.694
3 3 3 74.5 82.6 76 77.7 4.31 2.49
4 4 3 67.3 78.1 72.7 72.7 5.4 3.12
5 5 3 73.1 74.2 73.2 73.5 0.608 0.351
6 6 3 65.8 70.8 68.7 68.4 2.51 1.45
7 7 3 69 84.9 75.6 76.5 7.99 4.61 Bond data analyzed as a two-way ANOVA
One way to analyze these data is through a two-way ANOVA. Modeled as a two-way anova, but it’s not two between subjects and not repeated measures, because each observation is independent. The main effect of Metal is significant and indicates that the mean breaking pressure of each metal are not equal to one another. The main effect for ingot is also significant which indicates that the mean breaking pressure across ingots are not equal.
mod_1.2 <- aov(pressure ~ Metal + Ingot, data = bond)
broom::tidy(anova(mod_1.2)) %>%
format_gt_tbl()| term | df | sumsq | meansq | statistic | p.value |
|---|---|---|---|---|---|
| Metal | 2 | 131.90 | 65.95 | 6.36 | 0.013 |
| Ingot | 6 | 268.29 | 44.71 | 4.31 | 0.015 |
| Residuals | 12 | 124.46 | 10.37 |
The output of the fixed effects from the ANOVA reveals a statistically significant effect of Metal which indicates that there strong evidence that the breaking pressure means across metals are not equal. Similarly, the main effect of Ingot reveals a statistically significan effect that suggest the mean breaking pressure between Ingots are not equal.
Bond data set analyzed as mixed effects regression model
# Modeled as a mixed model
mod_1.3 <- lmerTest::lmer(pressure ~ Metal + (1 | Ingot), data = bond)
broom.mixed::tidy(mod_1.3)# A tibble: 5 × 8
effect group term estimate std.error statistic df p.value
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) 70.2 1.77 39.8 11.6 9.05e-14
2 fixed <NA> Metali 5.71 1.72 3.32 12.0 6.12e- 3
3 fixed <NA> Metaln 0.914 1.72 0.531 12.0 6.05e- 1
4 ran_pars Ingot sd__(Intercept) 3.38 NA NA NA NA
5 ran_pars Residual sd__Observation 3.22 NA NA NA NA # Both Metal and Ingot are coded as factors
# 1) Fixed effect (Intercept) represents the average predicted pressure for the
# reference level (copper).
# 2) Fixed effect Metali represents the difference between copper and iron. The
# difference is statistically significant which means that pressure required to
# break the iron ingot is 5.71 units higher than copper.
# 3) Fixed effect Metaln is the difference between nickel and copper. The finding
# is not significant. There is no evidence to suggest that copper and nickel are
# different.
# 4) Random effect sd__(Intercept) the estimated standard deviation of a distribution of
# intercepts centered on 70.2. Which is essentially the distribution of the means
# of iron from all ingots.
# 5) Random effect sd__Observation is the sd of there residuals.
broom::tidy(anova(mod_1.3))Warning: The column names NumDF and DenDF in ANOVA output were not recognized or
transformed.# A tibble: 1 × 7
term sumsq meansq NumDF DenDF statistic p.value
<chr> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
1 Metal 132. 66.0 2 12.0 6.36 0.0131# Solution to the fixed effect of Metal which tests if there is a significant
# differencs between the means of Metali, Metalc, and Metaln. The fixed effect
# is significant p = 0.0131, which indicates that there is at least one
# signficant pairwise difference between the different metals. There is strong
# evidence of differences between metal means.
# These results are the same as the aov() approach, because the data are
# balanced for the main effect of metal.emm_options(lmer.df = "satterthwaite")
# All emmeans and pairwise diffs
emm <- emmeans(mod_1.3, list(pairwise ~ Metal), adjust = "none")
summary(emm)$`emmeans of Metal`
Metal emmean SE df lower.CL upper.CL
c 70.2 1.77 11.6 66.3 74.0
i 75.9 1.77 11.6 72.0 79.8
n 71.1 1.77 11.6 67.2 75.0
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
$`pairwise differences of Metal`
1 estimate SE df t.ratio p.value
c - i -5.714 1.72 12 -3.320 0.0061
c - n -0.914 1.72 12 -0.531 0.6050
i - n 4.800 1.72 12 2.788 0.0164
Degrees-of-freedom method: satterthwaite # Test if Nickel's mean is different from 0 like SAS estimate
contrast(emm, method = list("Nickel vs 0" = c(0, 0, 1)), adjust = "none") contrast estimate SE df t.ratio p.value
Nickel vs 0 71.1 1.77 11.6 40.271 <.0001
Degrees-of-freedom method: satterthwaite # Test if Copper's mean is different from Iron like SAS estimate
contrast(emm, method = list("Copper vs Iron" = c(1, -1, 0)), adjust = "none") contrast estimate SE df t.ratio p.value
Copper vs Iron -5.71 1.72 12 -3.320 0.0061
Degrees-of-freedom method: satterthwaite # SAS contrast is a companion to SAS estimate# As_tibble automatically converts to factors
pbib <- as_tibble(PBIB)
# Order the levels of Treatment and Block
pbib <- pbib %>%
mutate(across(Treatment:Block, ~ factor(.x, levels = c(1:15))))mod_1.5.2 <- lmerTest::lmer(response ~ Treatment + (1 | Block), data = pbib)
# broom.mixed::tidy(mod_1.5.2)
# Solution for fixed effects, will match SAS when ddfm is set to satterwaithe
broom::tidy(anova(mod_1.5.2))Warning: The column names NumDF and DenDF in ANOVA output were not recognized or
transformed.# A tibble: 1 × 7
term sumsq meansq NumDF DenDF statistic p.value
<chr> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
1 Treatment 1.83 0.131 14 36.2 1.53 0.149# All estimated marginal means and pairwise differences by Treatment
emm <- emmeans(mod_1.5.2, list(pairwise ~ Treatment), adjust = "none")
summary(emm)$`emmeans of Treatment`
Treatment emmean SE df lower.CL upper.CL
1 2.82 0.166 44.1 2.48 3.15
2 2.41 0.166 44.1 2.07 2.74
3 2.45 0.166 44.1 2.12 2.79
4 2.78 0.166 44.1 2.45 3.12
5 2.80 0.166 44.1 2.47 3.14
6 2.91 0.166 44.1 2.58 3.25
7 2.79 0.166 44.1 2.45 3.12
8 2.78 0.166 44.1 2.45 3.12
9 2.89 0.166 44.1 2.56 3.23
10 2.49 0.166 44.1 2.16 2.83
11 2.90 0.166 44.1 2.56 3.23
12 3.05 0.166 44.1 2.72 3.39
13 2.62 0.166 44.1 2.28 2.95
14 2.49 0.166 44.1 2.16 2.83
15 2.86 0.166 44.1 2.52 3.19
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
$`pairwise differences of Treatment`
1 estimate SE df t.ratio p.value
Treatment1 - Treatment2 0.412208 0.222 36.2 1.856 0.0716
Treatment1 - Treatment3 0.362579 0.222 36.2 1.633 0.1112
Treatment1 - Treatment4 0.033692 0.222 36.2 0.152 0.8802
Treatment1 - Treatment5 0.012625 0.222 36.2 0.057 0.9550
Treatment1 - Treatment6 -0.093171 0.227 37.7 -0.410 0.6841
Treatment1 - Treatment7 0.028538 0.222 36.2 0.129 0.8985
Treatment1 - Treatment8 0.035917 0.222 36.2 0.162 0.8724
Treatment1 - Treatment9 -0.073789 0.222 36.2 -0.332 0.7416
Treatment1 - Treatment10 0.326461 0.222 36.2 1.470 0.1501
Treatment1 - Treatment11 -0.081176 0.227 37.7 -0.357 0.7229
Treatment1 - Treatment12 -0.235299 0.222 36.2 -1.060 0.2963
Treatment1 - Treatment13 0.199753 0.222 36.2 0.900 0.3743
Treatment1 - Treatment14 0.326211 0.222 36.2 1.469 0.1505
Treatment1 - Treatment15 -0.041710 0.222 36.2 -0.188 0.8521
Treatment2 - Treatment3 -0.049628 0.222 36.2 -0.223 0.8244
Treatment2 - Treatment4 -0.378515 0.222 36.2 -1.705 0.0968
Treatment2 - Treatment5 -0.399583 0.222 36.2 -1.799 0.0803
Treatment2 - Treatment6 -0.505379 0.222 36.2 -2.276 0.0289
Treatment2 - Treatment7 -0.383670 0.227 37.7 -1.689 0.0995
Treatment2 - Treatment8 -0.376291 0.222 36.2 -1.695 0.0987
Treatment2 - Treatment9 -0.485996 0.222 36.2 -2.189 0.0352
Treatment2 - Treatment10 -0.085747 0.222 36.2 -0.386 0.7016
Treatment2 - Treatment11 -0.493384 0.222 36.2 -2.222 0.0326
Treatment2 - Treatment12 -0.647506 0.227 37.7 -2.850 0.0070
Treatment2 - Treatment13 -0.212454 0.222 36.2 -0.957 0.3450
Treatment2 - Treatment14 -0.085996 0.222 36.2 -0.387 0.7008
Treatment2 - Treatment15 -0.453918 0.222 36.2 -2.044 0.0483
Treatment3 - Treatment4 -0.328887 0.222 36.2 -1.481 0.1472
Treatment3 - Treatment5 -0.349955 0.222 36.2 -1.576 0.1237
Treatment3 - Treatment6 -0.455751 0.222 36.2 -2.052 0.0474
Treatment3 - Treatment7 -0.334042 0.222 36.2 -1.504 0.1412
Treatment3 - Treatment8 -0.326662 0.227 37.7 -1.438 0.1587
Treatment3 - Treatment9 -0.436368 0.222 36.2 -1.965 0.0571
Treatment3 - Treatment10 -0.036118 0.222 36.2 -0.163 0.8717
Treatment3 - Treatment11 -0.443756 0.222 36.2 -1.998 0.0532
Treatment3 - Treatment12 -0.597878 0.222 36.2 -2.692 0.0107
Treatment3 - Treatment13 -0.162826 0.227 37.7 -0.717 0.4780
Treatment3 - Treatment14 -0.036368 0.222 36.2 -0.164 0.8708
Treatment3 - Treatment15 -0.404290 0.222 36.2 -1.821 0.0769
Treatment4 - Treatment5 -0.021068 0.222 36.2 -0.095 0.9249
Treatment4 - Treatment6 -0.126863 0.222 36.2 -0.571 0.5713
Treatment4 - Treatment7 -0.005155 0.222 36.2 -0.023 0.9816
Treatment4 - Treatment8 0.002225 0.222 36.2 0.010 0.9921
Treatment4 - Treatment9 -0.107481 0.227 37.7 -0.473 0.6389
Treatment4 - Treatment10 0.292769 0.222 36.2 1.318 0.1956
Treatment4 - Treatment11 -0.114869 0.222 36.2 -0.517 0.6081
Treatment4 - Treatment12 -0.268991 0.222 36.2 -1.211 0.2336
Treatment4 - Treatment13 0.166061 0.222 36.2 0.748 0.4594
Treatment4 - Treatment14 0.292519 0.227 37.7 1.287 0.2058
Treatment4 - Treatment15 -0.075403 0.222 36.2 -0.340 0.7361
Treatment5 - Treatment6 -0.105796 0.222 36.2 -0.476 0.6366
Treatment5 - Treatment7 0.015913 0.222 36.2 0.072 0.9433
Treatment5 - Treatment8 0.023292 0.222 36.2 0.105 0.9170
Treatment5 - Treatment9 -0.086413 0.222 36.2 -0.389 0.6994
Treatment5 - Treatment10 0.313836 0.227 37.7 1.381 0.1753
Treatment5 - Treatment11 -0.093801 0.222 36.2 -0.422 0.6752
Treatment5 - Treatment12 -0.247923 0.222 36.2 -1.116 0.2716
Treatment5 - Treatment13 0.187129 0.222 36.2 0.843 0.4049
Treatment5 - Treatment14 0.313587 0.222 36.2 1.412 0.1664
Treatment5 - Treatment15 -0.054335 0.227 37.7 -0.239 0.8123
Treatment6 - Treatment7 0.121709 0.222 36.2 0.548 0.5870
Treatment6 - Treatment8 0.129088 0.222 36.2 0.581 0.5646
Treatment6 - Treatment9 0.019383 0.222 36.2 0.087 0.9309
Treatment6 - Treatment10 0.419632 0.222 36.2 1.890 0.0668
Treatment6 - Treatment11 0.011995 0.227 37.7 0.053 0.9582
Treatment6 - Treatment12 -0.142127 0.222 36.2 -0.640 0.5262
Treatment6 - Treatment13 0.292925 0.222 36.2 1.319 0.1954
Treatment6 - Treatment14 0.419383 0.222 36.2 1.889 0.0670
Treatment6 - Treatment15 0.051461 0.222 36.2 0.232 0.8180
Treatment7 - Treatment8 0.007379 0.222 36.2 0.033 0.9737
Treatment7 - Treatment9 -0.102326 0.222 36.2 -0.461 0.6477
Treatment7 - Treatment10 0.297923 0.222 36.2 1.342 0.1881
Treatment7 - Treatment11 -0.109714 0.222 36.2 -0.494 0.6242
Treatment7 - Treatment12 -0.263836 0.227 37.7 -1.161 0.2528
Treatment7 - Treatment13 0.171216 0.222 36.2 0.771 0.4457
Treatment7 - Treatment14 0.297674 0.222 36.2 1.341 0.1884
Treatment7 - Treatment15 -0.070248 0.222 36.2 -0.316 0.7536
Treatment8 - Treatment9 -0.109706 0.222 36.2 -0.494 0.6243
Treatment8 - Treatment10 0.290544 0.222 36.2 1.308 0.1990
Treatment8 - Treatment11 -0.117094 0.222 36.2 -0.527 0.6012
Treatment8 - Treatment12 -0.271216 0.222 36.2 -1.221 0.2298
Treatment8 - Treatment13 0.163836 0.227 37.7 0.721 0.4753
Treatment8 - Treatment14 0.290294 0.222 36.2 1.307 0.1994
Treatment8 - Treatment15 -0.077628 0.222 36.2 -0.350 0.7287
Treatment9 - Treatment10 0.400249 0.222 36.2 1.802 0.0798
Treatment9 - Treatment11 -0.007388 0.222 36.2 -0.033 0.9736
Treatment9 - Treatment12 -0.161510 0.222 36.2 -0.727 0.4717
Treatment9 - Treatment13 0.273542 0.222 36.2 1.232 0.2259
Treatment9 - Treatment14 0.400000 0.227 37.7 1.761 0.0864
Treatment9 - Treatment15 0.032078 0.222 36.2 0.144 0.8859
Treatment10 - Treatment11 -0.407637 0.222 36.2 -1.836 0.0746
Treatment10 - Treatment12 -0.561760 0.222 36.2 -2.530 0.0159
Treatment10 - Treatment13 -0.126708 0.222 36.2 -0.571 0.5718
Treatment10 - Treatment14 -0.000249 0.222 36.2 -0.001 0.9991
Treatment10 - Treatment15 -0.368171 0.227 37.7 -1.620 0.1135
Treatment11 - Treatment12 -0.154122 0.222 36.2 -0.694 0.4921
Treatment11 - Treatment13 0.280930 0.222 36.2 1.265 0.2139
Treatment11 - Treatment14 0.407388 0.222 36.2 1.835 0.0748
Treatment11 - Treatment15 0.039466 0.222 36.2 0.178 0.8599
Treatment12 - Treatment13 0.435052 0.222 36.2 1.959 0.0578
Treatment12 - Treatment14 0.561510 0.222 36.2 2.529 0.0159
Treatment12 - Treatment15 0.193588 0.222 36.2 0.872 0.3891
Treatment13 - Treatment14 0.126458 0.222 36.2 0.569 0.5725
Treatment13 - Treatment15 -0.241464 0.222 36.2 -1.087 0.2840
Treatment14 - Treatment15 -0.367922 0.222 36.2 -1.657 0.1062
Degrees-of-freedom method: satterthwaite # Estimate statements
# 1) treat 1 mean. Matches SAS w ddfm satterwaithe, but does not give a t-value for the test of whether or not the mean differs from 0
emmeans(mod_1.5.2, ~ Treatment) %>%
as_tibble() %>%
filter(Treatment == 1)# A tibble: 1 × 6
Treatment emmean SE df lower.CL upper.CL
<fct> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 2.82 0.166 44.1 2.48 3.15# In R to get the test of treat 1 mean we need a contrast function
contrast(emm,
method = list("Treat 1 vs 0" = c(1, rep(0, length(levels(pbib$Treatment)) - 1))),
adjust = "none") contrast estimate SE df t.ratio p.value
Treat 1 vs 0 2.82 0.166 44.1 16.931 <.0001
Degrees-of-freedom method: satterthwaite # 2) trt 1 mean, including the random effect
# The following produces the same estimate as #1, but no se.
# Create newdata for Treatment 1 across all blocks
newdata <- expand.grid(Treatment = "1", Block = unique(pbib$Block))
# Predictions including random effects (BLUPs)
preds <- predict(mod_1.5.2, newdata = newdata, re.form = NULL)
# Average across blocks
mean(preds)[1] 2.817523f <- function(fit) {
mean(predict(fit, newdata = newdata, re.form = NULL))
}
set.seed(123)
bootres <- lme4::bootMer(mod_1.5.2, f, nsim = 1000)
mean(bootres$t) # estimate[1] 2.812753sd(bootres$t) # bootstrap SE (closest to SAS)[1] 0.1638689# 3) trt 1 blk 1
# Estimate the mean of treat 1 blk 1
# The following produced the correct estimate, but no se or p value
# Create a new data frame with the combination you want
newdata <- data.frame(Treatment = "1", Block = "1")
# Predict the estimated mean trt 1 blk 1
predict(mod_1.5.2, newdata = newdata, re.form = NULL) # includes random effect for Block1 1
2.528808 # 4) trt 1 vs trt 2
contrast(emm,
method = list("Treat 1 vs Treat 2" = c(1, -1, rep(0, length(levels(pbib$Treatment)) - 2))),
adjust = "none") contrast estimate SE df t.ratio p.value
Treat 1 vs Treat 2 0.412 0.222 36.2 1.856 0.0716
Degrees-of-freedom method: satterthwaite Chapter 2 - Common Mixed Models
semiconductor <- as_tibble(Semiconductor)
mod_2.3.1 <- lmerTest::lmer(resistance ~ ET * position + (1 | ET:Wafer), data = semiconductor)
# summary(mod_2.3.1)
# broom.mixed::tidy(mod_2.3.1)
broom::tidy(anova(mod_2.3.1))Warning: The column names NumDF and DenDF in ANOVA output were not recognized or
transformed.# A tibble: 3 × 7
term sumsq meansq NumDF DenDF statistic p.value
<chr> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
1 ET 0.647 0.216 3 8.00 1.94 0.202
2 position 1.13 0.376 3 24.0 3.39 0.0345
3 ET:position 0.809 0.0899 9 24.0 0.809 0.613 broom.mixed::glance(mod_2.3.1)# A tibble: 1 × 7
nobs sigma logLik AIC BIC REMLcrit df.residual
<int> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
1 48 0.333 -25.3 86.7 120. 50.7 30# LS means for ET
emmeans(mod_2.3.1, ~ ET)NOTE: Results may be misleading due to involvement in interactions ET emmean SE df lower.CL upper.CL
1 5.63 0.211 8 5.14 6.11
2 5.97 0.211 8 5.48 6.45
3 6.09 0.211 8 5.60 6.57
4 6.33 0.211 8 5.85 6.82
Results are averaged over the levels of: position
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95 # LS means for position
emmeans(mod_2.3.1, ~ position)NOTE: Results may be misleading due to involvement in interactions position emmean SE df lower.CL upper.CL
1 6.02 0.134 18.7 5.74 6.30
2 6.13 0.134 18.7 5.85 6.42
3 5.75 0.134 18.7 5.47 6.03
4 6.11 0.134 18.7 5.83 6.39
Results are averaged over the levels of: ET
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95 # LS means for ET * Position
# Swap the terms to get the output to mimic SAS
emmeans(mod_2.3.1, ~ position:ET) position ET emmean SE df lower.CL upper.CL
1 1 5.61 0.269 18.7 5.05 6.18
2 1 5.45 0.269 18.7 4.89 6.01
3 1 5.55 0.269 18.7 4.99 6.12
4 1 5.89 0.269 18.7 5.32 6.45
1 2 5.99 0.269 18.7 5.43 6.56
2 2 6.19 0.269 18.7 5.62 6.75
3 2 5.77 0.269 18.7 5.20 6.33
4 2 5.92 0.269 18.7 5.35 6.48
1 3 6.14 0.269 18.7 5.57 6.70
2 3 6.35 0.269 18.7 5.78 6.91
3 3 5.77 0.269 18.7 5.21 6.34
4 3 6.09 0.269 18.7 5.53 6.66
1 4 6.34 0.269 18.7 5.78 6.90
2 4 6.55 0.269 18.7 5.99 7.12
3 4 5.90 0.269 18.7 5.33 6.46
4 4 6.54 0.269 18.7 5.98 7.10
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95 # Estimates -------------------------------------------------------------------
# 1)
# Get marginal means for ET
emm <- emmeans(mod_2.3.1, ~ ET)NOTE: Results may be misleading due to involvement in interactions# Contrast ET1 vs ET2
contrast(emm, method = list("et1 vs et2" = c(1, -1, 0, 0))) contrast estimate SE df t.ratio p.value
et1 vs et2 -0.34 0.298 8 -1.139 0.2875
Results are averaged over the levels of: position
Degrees-of-freedom method: satterthwaite # 2)
emm <- emmeans(mod_2.3.1, ~ position)NOTE: Results may be misleading due to involvement in interactions# Contrast ET1 vs ET2
contrast(emm, method = list("pos1 vs pos2" = c(1, -1, 0, 0))) contrast estimate SE df t.ratio p.value
pos1 vs pos2 -0.113 0.136 24 -0.833 0.4132
Results are averaged over the levels of: ET
Degrees-of-freedom method: satterthwaite # 3) Simple effect contrast pos1 vs pos 2 in et1
emm <- emmeans(mod_2.3.1, ~ position | ET)
# Contrast pos1 vs pos2 within ET = 1
contrast(emm, method = list("pos1 vs pos2 by et" = c(1, -1, 0, 0)), by = "ET")ET = 1:
contrast estimate SE df t.ratio p.value
pos1 vs pos2 by et 0.163 0.272 24 0.600 0.5541
ET = 2:
contrast estimate SE df t.ratio p.value
pos1 vs pos2 by et -0.193 0.272 24 -0.710 0.4844
ET = 3:
contrast estimate SE df t.ratio p.value
pos1 vs pos2 by et -0.210 0.272 24 -0.771 0.4480
ET = 4:
contrast estimate SE df t.ratio p.value
pos1 vs pos2 by et -0.213 0.272 24 -0.784 0.4409
Degrees-of-freedom method: satterthwaite #4 )
emm <- emmeans(mod_2.3.1, ~ ET | position)
# Contrast pos1 vs pos2 within ET = 1
contrast(emm, method = list("et1 vs et2 by pos" = c(1, -1, 0, 0)),
by = "position") %>%
as_tibble() %>%
filter(position == 1)# A tibble: 1 × 7
contrast position estimate SE df t.ratio p.value
<fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
1 et1 vs et2 by pos 1 -0.380 0.380 18.7 -0.999 0.330# To get the confidence interval use the formula: [estimate +/ t_crit * SE]
# round 18.7 down to 18 to be more conservative
t_crit <- qt(.975, 18)
upr <- -0.380 + t_crit * 0.38
lwr <- -0.380 - t_crit * 0.38
# 5)
contrast(emm, method = list("et1,et2 vs et3,4 by pos" = c(.5, .5, -.5, -.5)),
by = "position") %>%
as_tibble() %>%
filter(position == 2)# A tibble: 1 × 7
contrast position estimate SE df t.ratio p.value
<fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
1 et1,et2 vs et3,4 by pos 2 -0.632 0.269 18.7 -2.35 0.0300# For the main effect of ET, the 0.05 critical value can be determined for the following, which can then be used in the formula: [estimate +/ t_crit * SE] to create confidence intervals. These are also displayed by emmeans.
qt(0.975, df = 8)[1] 2.306004## Simple effect estimation and testing
pairs(emm, adjust = "none")position = 1:
contrast estimate SE df t.ratio p.value
ET1 - ET2 -0.38000 0.38 18.7 -0.999 0.3305
ET1 - ET3 -0.52333 0.38 18.7 -1.376 0.1851
ET1 - ET4 -0.72667 0.38 18.7 -1.911 0.0715
ET2 - ET3 -0.14333 0.38 18.7 -0.377 0.7105
ET2 - ET4 -0.34667 0.38 18.7 -0.912 0.3736
ET3 - ET4 -0.20333 0.38 18.7 -0.535 0.5992
position = 2:
contrast estimate SE df t.ratio p.value
ET1 - ET2 -0.73667 0.38 18.7 -1.937 0.0680
ET1 - ET3 -0.89667 0.38 18.7 -2.358 0.0295
ET1 - ET4 -1.10333 0.38 18.7 -2.901 0.0093
ET2 - ET3 -0.16000 0.38 18.7 -0.421 0.6788
ET2 - ET4 -0.36667 0.38 18.7 -0.964 0.3473
ET3 - ET4 -0.20667 0.38 18.7 -0.543 0.5933
position = 3:
contrast estimate SE df t.ratio p.value
ET1 - ET2 -0.21333 0.38 18.7 -0.561 0.5815
ET1 - ET3 -0.22000 0.38 18.7 -0.578 0.5698
ET1 - ET4 -0.34333 0.38 18.7 -0.903 0.3781
ET2 - ET3 -0.00667 0.38 18.7 -0.018 0.9862
ET2 - ET4 -0.13000 0.38 18.7 -0.342 0.7363
ET3 - ET4 -0.12333 0.38 18.7 -0.324 0.7493
position = 4:
contrast estimate SE df t.ratio p.value
ET1 - ET2 -0.03000 0.38 18.7 -0.079 0.9380
ET1 - ET3 -0.20667 0.38 18.7 -0.543 0.5933
ET1 - ET4 -0.65333 0.38 18.7 -1.718 0.1023
ET2 - ET3 -0.17667 0.38 18.7 -0.465 0.6476
ET2 - ET4 -0.62333 0.38 18.7 -1.639 0.1179
ET3 - ET4 -0.44667 0.38 18.7 -1.175 0.2549
Degrees-of-freedom method: satterthwaite