library(tidyverse)
library(gtsummary)
library(SASmixed)
library(emmeans)
library(dfmtbx)SAS for Mixed Models but in R
Chapter 1 - Randomized Blocks Designs
A randomized complete block design is an experimental design where a treatment is randomized within blocks, and each block receives all possible treatments of interest.
Suppose we wanted to study the effect of a treatment on measures of depression. During recruiting, we may wish to block on sex. That is we would randomize each male to one of the treatments (e.g. active treatment or placebo) so that roughly equal numbers of males are represented in each condition. Similarly, we would randomize each female into one of the conditions for the same reason. The aim is to control for potential differences between males and females. This approach ensures that each treatment condition is represented within each block, which helps isolate the treatment effect while accounting for variability due to sex. This setup exemplifies a randomized complete block design (RCBD).
Bond data set description
In chapter 1 of Littell et al. 1996, the authors present the Bond data set. This data set can be obtained from the SASmixed package in R and is also in the sasuser library in SAS version 9.4.
The data come from a study that investigated the pressure required to break a bond between two metal components made of three different elements (nickel, iron, and copper). The metal components also vary in their physical shape, referred to as ingot type, which serves as the blocking factor in the design. The study includes seven distinct ingot types. Within each ingot type, all three metals are represented. This results in a total of 21 observations (three per ingot) where the pressure to break the bond between the two bonded components is recorded as the outcome variable. The outcome is measured in thousands of pounds per square inch. The objective of the study is to determine which metal forms the strongest bond, while controlling for variation due to ingot shape. Since all three treatments are represented in each ingot, the study exemplifies a randomized complete block design.
# Set all column names to lower case
names(Bond) <- tolower(names(Bond))
# as_tibble converts values to factor
bond <- Bond %>%
mutate(metal = factor(case_match(metal,
"c" ~ "copper",
"i" ~ "iron",
"n" ~ "nickel"))) %>%
as_tibble()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> nickel, iron, copper, nickel, iron, copper, nickel, iron, cop…
$ 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 | nickel | 1 |
| 71.90 | iron | 1 |
| 72.20 | copper | 1 |
| 67.50 | nickel | 2 |
| 68.80 | iron | 2 |
| 66.40 | copper | 2 |
bond %>%
group_by(metal) %>%
summarize_variable(pressure) %>%
format_gt_tbl()| metal | n | min | max | median | mean | sd | se |
|---|---|---|---|---|---|---|---|
| copper | 7 | 66.40 | 74.50 | 69.00 | 70.19 | 3.11 | 1.18 |
| iron | 7 | 68.80 | 84.90 | 74.20 | 75.90 | 6.14 | 2.32 |
| nickel | 7 | 65.80 | 76.00 | 72.70 | 71.10 | 4.26 | 1.61 |
bond %>%
group_by(ingot) %>%
summarize_variable(pressure) %>%
format_gt_tbl()| ingot | n | min | max | median | mean | sd | se |
|---|---|---|---|---|---|---|---|
| 1 | 3 | 67.00 | 72.20 | 71.90 | 70.37 | 2.92 | 1.69 |
| 2 | 3 | 66.40 | 68.80 | 67.50 | 67.57 | 1.20 | 0.69 |
| 3 | 3 | 74.50 | 82.60 | 76.00 | 77.70 | 4.31 | 2.49 |
| 4 | 3 | 67.30 | 78.10 | 72.70 | 72.70 | 5.40 | 3.12 |
| 5 | 3 | 73.10 | 74.20 | 73.20 | 73.50 | 0.61 | 0.35 |
| 6 | 3 | 65.80 | 70.80 | 68.70 | 68.43 | 2.51 | 1.45 |
| 7 | 3 | 69.00 | 84.90 | 75.60 | 76.50 | 7.99 | 4.61 |
Bond data analyzed as a two-way ANOVA
One way to analyze these data is through a two-way ANOVA, specifying both Metal and Ingot as fixed effects. The results of this analysis revealed a significant effect of Metal (𝑝=0.013), indicating evidence that the mean pressure values differ across metal types. Similarly, there is a significant effect of Ingot (𝑝=0.015), suggesting that the mean pressure values vary across ingot types. These findings support the conclusion that both the type of metal and the ingot shape influence the pressure required to break the bonds between the metal components tested.
mod_1.2 <- aov(pressure ~ metal + ingot, data = bond)
anova(mod_1.2) %>%
broom::tidy() %>%
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 |
Bond data set analyzed as mixed effects regression model
One drawback of analyzing these data using the ANOVA framework is that the inferences are limited to the specific metals and ingot types included in the study. If we wished to generalize our findings to a broader population of ingot shapes beyond the seven tested, we can instead use a mixed effects model. In this approach, ingot is specified as a random effect, allowing us to treat the seven ingots as a representative sample from a larger population. This broadens the scope of the inference and accounts for variability due to ingot type without estimating a fixed effect for each ingot. When we don’t have to estimate a parameter for each ingot, we have more degrees of freedom left over for estimating error variance and greater power to detect treatment effects, which results in greater statistical efficiency.
# A mixed effects regression model can be specified with the lmerTest() package. It's formula syntax is similar to that of the base R lm() function. The "predictor" in the formulat (1 | Ingot) specifies that we wish to model a random intercept for each Ingot.# Modeled as a mixed model
mod_1.3 <- lmerTest::lmer(pressure ~ metal + (1 | ingot), data = bond)Model fitting information
To display a summary of measure like goodness of fit or model convergence information, the glance() function can be used.
mod_1.3 %>%
broom.mixed::glance() %>%
format_gt_tbl()| nobs | sigma | logLik | AIC | BIC | REMLcrit | df.residual |
|---|---|---|---|---|---|---|
| 21 | 3.22 | −53.90 | 117.79 | 123.01 | 107.79 | 16 |
Type 3 Tests of Fixed Effects
The type 3 tests of fixed effects can be displayed by passing the model object to the anova() function. The output is in turn piped into the tidy() function from the broom package to return a tibble. The type 3 test of fixed effect reveals whether or not all metal means are equal. The fixed effect for metal 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 mirror those obtained in the ANOVA approach because the data are balanced for the main effect of metal, but this may not always be the case.
anova(mod_1.3) %>%
broom::tidy() %>%
format_gt_tbl()Warning: The column names NumDF and DenDF in ANOVA output were not recognized or
transformed.| term | sumsq | meansq | NumDF | DenDF | statistic | p.value |
|---|---|---|---|---|---|---|
| metal | 131.90 | 65.95 | 2 | 12.00 | 6.36 | 0.013 |
Model output
Model output including the parameter estimates can be displayed with the broom.mixed::tidy() function. Since metal is coded as a categorical variable, the interpretation of the fixed effects is slightly different than that of a continuous variable. 1) The estimate for the fixed effect (Intercept) represents the average predicted pressure for the reference level metal since there is only one fixed effect. In this case the reference level is copper. Thus, 70.2 represents the mean predicted pressure to break a copper bond. 2) The estimate for the fixed effect Metali represents the difference between the reference (copper) and iron. The difference is statistically significant (p = 0.00612) which means that on average the pressure required to break the iron ingot is 5.71 units higher than that of copper. 3) The estimate for the fixed effect Metaln is the difference between nickel and copper. The finding is not significant (p = 0.605). There is no evidence to suggest that copper and nickel have different breaking pressures. 4) The random effect sd__(Intercept) represents the estimated standard deviation of the ingot-specific intercepts in the mixed model for the reference level. It captures how much the average pressure for copper might vary from ingot to ingot due to random ingot-level effects. 5) The Random effect sd__Observation is the standard deviation of the residuals.
mod_1.3 %>%
broom.mixed::tidy() %>%
format_gt_tbl()| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | (Intercept) | 70.19 | 1.77 | 39.75 | 11.61 | < 0.001 | |
| fixed | metaliron | 5.71 | 1.72 | 3.32 | 12.00 | 0.006 | |
| fixed | metalnickel | 0.91 | 1.72 | 0.53 | 12.00 | 0.605 | |
| ran_pars | ingot | sd__(Intercept) | 3.38 | ||||
| ran_pars | Residual | sd__Observation | 3.22 |
The model output above already displays the pairwise differences between copper and iron, and between copper and nickel, because copper is the reference. To display all pairwise differences, we can use the emmeans package.
# Set the degrees of freedom to Satterthwaite otherwise emmeans() defaults to Kenward-Rogers
emm_options(lmer.df = "satterthwaite")
# All emmeans and pairwise diffs
emm_grid <- emmeans(mod_1.3, ~ metal) Here, the output of the pairs() function can be used as input to the tidy() function from the broom package. The first row displays the pairwise difference between copper (c) and iron (i) which is the same estimate, standard error, statistic, and p-value displayed in the model output. Similarly, the second rows displays the pairwise differences between copper and nickel and displays the same values from the model output. The pairwise difference that is not displayed in the model output is that between iron and nickel.
# Display the pdiffs
broom::tidy(pairs(emm_grid, adjust = "none")) %>%
select(-null.value) %>%
format_gt_tbl()| term | contrast | estimate | std.error | df | statistic | p.value |
|---|---|---|---|---|---|---|
| metal | copper - iron | −5.71 | 1.72 | 12.00 | −3.32 | 0.006 |
| metal | copper - nickel | −0.91 | 1.72 | 12.00 | −0.53 | 0.605 |
| metal | iron - nickel | 4.80 | 1.72 | 12.00 | 2.79 | 0.016 |
In addition to displaying the pairwise differences, specific contrasts can be specified to test whether the mean of nickel is significantly different than 0. Instead of displaying all of the pairwise differences, one could also display the difference between copper and iron only to supplement the model output.
# Test if Nickel's mean is different from 0 like SAS estimate
contrast(emm_grid,
method = list("Nickel vs 0" = c(0, 0, 1)),
adjust = "none") %>%
broom::tidy() %>%
select(-null.value) %>%
format_gt_tbl()| term | contrast | estimate | std.error | df | statistic | p.value |
|---|---|---|---|---|---|---|
| metal | Nickel vs 0 | 71.10 | 1.77 | 11.61 | 40.27 | < 0.001 |
# Test if Copper's mean is different from Iron like SAS estimate
contrast(emm_grid,
method = list("Copper vs Iron" = c(1, -1, 0)),
adjust = "none") %>%
broom::tidy() %>%
select(-null.value) %>%
format_gt_tbl()| term | contrast | estimate | std.error | df | statistic | p.value |
|---|---|---|---|---|---|---|
| metal | Copper vs Iron | −5.71 | 1.72 | 12.00 | −3.32 | 0.006 |
Partially balanced incomplete block design
There are two broad categories of randomized block designs, complete and incomplete. In the randomized complete block design all treatment levels are represented in each block. However, there may be situations where it is not feasible to represent all treatments in each block. In these cases then these designs would be considered incomplete block designs.
There are two types of incomplete block designss, balanced (BIBD) and partially balanced (PBID). Suppose a situation where there are 4 treatments, A, B, C, and D and each block will have two treatments assigned. In the BIBD, each block will have a unique pair of treatments (e.g. AB, AC, AD, BC, and CD) and there will be an equal number of blocks for a given treatment combination. In the PBID, not all treatment combinations are represented equally.
The following data set comes from a PBID where the blocking factor contains 15 levels (15 block types). Each block contains 4 out of 15 treatments that are randomly assigned. Each treatment is represented 4 times. The design is icomplete because not all treatments are represented in each block, and the design is partially balanced because not all treatment combinations within a block are represented equally, although each treatment is represented the same number of times.
To make it more concrete the blocking variable is a plot of land. Each plot of land is divided into 4 sub-areas where one treatment is applied. The outcome (response) is pounds of seed cotton per plot.
# As_tibble automatically converts to factors
pbib <- as_tibble(PBIB)
names(pbib) <- tolower(names(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)Model fitting information
mod_1.5.2 %>%
broom.mixed::glance() %>%
format_gt_tbl()| nobs | sigma | logLik | AIC | BIC | REMLcrit | df.residual |
|---|---|---|---|---|---|---|
| 60 | 0.293 | −25.99 | 85.98 | 121.59 | 51.98 | 43 |
Tests of Fixed Effects
The results from the Type III tests of fixed effects revealed no statistically significant effect of treatment. This suggests that the mean pounds of seed cotton do not differ significantly across treatment groups.
# Solution for fixed effects, will match SAS when ddfm is set to satterwaithe
anova(mod_1.5.2) %>%
broom::tidy() %>%
format_gt_tbl()Warning: The column names NumDF and DenDF in ANOVA output were not recognized or
transformed.| term | sumsq | meansq | NumDF | DenDF | statistic | p.value |
|---|---|---|---|---|---|---|
| treatment | 1.83 | 0.131 | 14 | 36.23 | 1.53 | 0.149 |
Model output
Similar to the RCBD model above the fixed effect estimate for the (Intercept) represents the predicted value for treatment1. Each subsequent fixed effect estimate for a given treatment is the difference between treatment1 and the displayed treatment. For example the row displaying treatment7 indicates that the differences between treatment1 and treatment7 is -0.03 units and suggests that the predicted value for treatment7 is 2.82 - 0.03, which equals 2.79.
mod_1.5.2 %>%
broom.mixed::tidy() %>%
format_gt_tbl()| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | (Intercept) | 2.82 | 0.166 | 16.93 | 44.14 | < 0.001 | |
| fixed | treatment2 | −0.41 | 0.222 | −1.86 | 36.25 | 0.072 | |
| fixed | treatment3 | −0.36 | 0.222 | −1.63 | 36.25 | 0.111 | |
| fixed | treatment4 | −0.03 | 0.222 | −0.15 | 36.25 | 0.880 | |
| fixed | treatment5 | −0.01 | 0.222 | −0.06 | 36.25 | 0.955 | |
| fixed | treatment6 | 0.09 | 0.227 | 0.41 | 37.70 | 0.684 | |
| fixed | treatment7 | −0.03 | 0.222 | −0.13 | 36.25 | 0.898 | |
| fixed | treatment8 | −0.04 | 0.222 | −0.16 | 36.25 | 0.872 | |
| fixed | treatment9 | 0.07 | 0.222 | 0.33 | 36.25 | 0.742 | |
| fixed | treatment10 | −0.33 | 0.222 | −1.47 | 36.25 | 0.150 | |
| fixed | treatment11 | 0.08 | 0.227 | 0.36 | 37.70 | 0.723 | |
| fixed | treatment12 | 0.24 | 0.222 | 1.06 | 36.25 | 0.296 | |
| fixed | treatment13 | −0.20 | 0.222 | −0.90 | 36.25 | 0.374 | |
| fixed | treatment14 | −0.33 | 0.222 | −1.47 | 36.25 | 0.150 | |
| fixed | treatment15 | 0.04 | 0.222 | 0.19 | 36.25 | 0.852 | |
| ran_pars | block | sd__(Intercept) | 0.22 | ||||
| ran_pars | Residual | sd__Observation | 0.29 |
Estimates & Contrasts
# All estimated marginal means and pairwise differences by treatment
emm <- emmeans(mod_1.5.2, ~ treatment)
# Display the lsmeans
broom::tidy(emm) %>%
format_gt_tbl()| treatment | estimate | std.error | df | statistic | p.value |
|---|---|---|---|---|---|
| 1 | 2.82 | 0.166 | 44.14 | 16.93 | < 0.001 |
| 2 | 2.41 | 0.166 | 44.14 | 14.45 | < 0.001 |
| 3 | 2.45 | 0.166 | 44.14 | 14.75 | < 0.001 |
| 4 | 2.78 | 0.166 | 44.14 | 16.73 | < 0.001 |
| 5 | 2.80 | 0.166 | 44.14 | 16.86 | < 0.001 |
| 6 | 2.91 | 0.166 | 44.14 | 17.49 | < 0.001 |
| 7 | 2.79 | 0.166 | 44.14 | 16.76 | < 0.001 |
| 8 | 2.78 | 0.166 | 44.14 | 16.72 | < 0.001 |
| 9 | 2.89 | 0.166 | 44.14 | 17.37 | < 0.001 |
| 10 | 2.49 | 0.166 | 44.14 | 14.97 | < 0.001 |
| 11 | 2.90 | 0.166 | 44.14 | 17.42 | < 0.001 |
| 12 | 3.05 | 0.166 | 44.14 | 18.34 | < 0.001 |
| 13 | 2.62 | 0.166 | 44.14 | 15.73 | < 0.001 |
| 14 | 2.49 | 0.166 | 44.14 | 14.97 | < 0.001 |
| 15 | 2.86 | 0.166 | 44.14 | 17.18 | < 0.001 |
# Display the first 5 and last 5 rows of pdiffs
broom::tidy(pairs(emm, adjust = "none")) %>%
as.data.frame() %>%
psych::headTail() %>%
mutate(contrast = ifelse(is.na(contrast), "...", contrast)) %>%
select(-term, -null.value) %>%
format_gt_tbl()| contrast | estimate | std.error | df | statistic | p.value |
|---|---|---|---|---|---|
| treatment1 - treatment2 | 0.41 | 0.22 | 36.25 | 1.86 | 0.07 |
| treatment1 - treatment3 | 0.36 | 0.22 | 36.25 | 1.63 | 0.11 |
| treatment1 - treatment4 | 0.03 | 0.22 | 36.25 | 0.15 | 0.88 |
| treatment1 - treatment5 | 0.01 | 0.22 | 36.25 | 0.06 | 0.95 |
| ... | ... | ... | ... | ... | ... |
| treatment12 - treatment15 | 0.19 | 0.22 | 36.25 | 0.87 | 0.39 |
| treatment13 - treatment14 | 0.13 | 0.22 | 36.25 | 0.57 | 0.57 |
| treatment13 - treatment15 | -0.24 | 0.22 | 36.25 | -1.09 | 0.28 |
| treatment14 - treatment15 | -0.37 | 0.22 | 36.25 | -1.66 | 0.11 |
# 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) %>%
format_gt_tbl()| treatment | emmean | SE | df | lower.CL | upper.CL |
|---|---|---|---|---|---|
| 1 | 2.82 | 0.166 | 44.14 | 2.48 | 3.15 |
# 2) trt 1 mean, including the random effect
# This is the mean of predicted treatment 1 values from all blocks
# 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
preds <- predict(mod_1.5.2, newdata = newdata, re.form = NULL)
# Average across blocks
mean(preds)[1] 2.817523# 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") %>%
as.data.frame() %>%
format_gt_tbl()| contrast | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|
| Treat 1 vs Treat 2 | 0.412 | 0.222 | 36.25 | 1.86 | 0.072 |
Chapter 2 - Common Mixed Models
Balanced Split-Plot Experiment
A split-plot experiment is a type of factorial design where some factors are applied to large experimental units (whole plots) while other factors are applied to smaller units (subplots) nested within them. This design is often implemented due to practical constraints as some treatments are difficult or costly to randomize across all units. An example from agricultural research is when the irrigation method to a large plot of land is difficult to change, whereas the fertilizer type is much easier to manipulate in smaller plots within the larger one. In this case the irrigation method is the whole-plot factor and the fertilizer type is the sub-plot factor. This type of design introduces two sources of variation that must be accounted for in the statistical analysis.
The semiconductor data set is an example of a balanced split-plot design from an industrial setting. In this experiment, we want to know if different manufacturing processes lead to differences in the electrical resistance of wafers used to make computer chips. Twelve wafers were drawn randomly from a lot and three wafers were randomly assigned to one of four different manufacturing processes (et). Within each wafer, the electrical resistance was measured at four different positions. Position represents the subplot as it is nested within wafer.
semiconductor <- as_tibble(Semiconductor)
names(semiconductor) <- tolower(names(semiconductor))Upon inspecting the data, one may notice that although there are 12 unique wafers, the numerical assignment for wafer is repeated within each treatment. To align the data in a framework that may be more frequently encountered in social and behavioral science research, and to make the model specification a bit more explicit, I decided to create a unique identifier.
table(semiconductor$wafer)
1 2 3
16 16 16 A unique identifier (waferID) can be created by concatenating the et and wafer variables. This results in an identifier where the first digit represents the manufacturing process and the second digit represents the wafer within that group. Another option would be to assign each wafer a unique values from 1 to 12, but the data need to be arranged by et and wafer before doing so. Both approaches work equally well.
# Create a unique wafer ID
semiconductor <- semiconductor %>%
mutate(waferID = as_factor(str_c(et, wafer)))
# Alternative way of creating waferID
# semiconductor %>%
# arrange(et, wafer) %>%
# mutate(waferID = as_factor(rep(1:12, each = 4)))
table(semiconductor$waferID)
11 12 13 21 22 23 31 32 33 41 42 43
4 4 4 4 4 4 4 4 4 4 4 4 Now that we have our unique identifier created, we can proceed to specify a linear mixed effects model. Since wafers are randomly selected and it would be wise to generalize beyond those 12 that were selected, we can specify waferID as our random effect and include a 1 to sepcify that each wafer gets their own intercept. For this example we will treat et and position as fixed effects, since there or only 4 treatments we are interested in and since position is consistent across wafers. Finally we add an interaction term of et by position to ensure that the manufacturing process does not lead to inconsistent results by position which would not be ideal in this type of setting.
mod_2.3.1 <- lmerTest::lmer(resistance ~ et + position + et:position + (1 | waferID), data = semiconductor)Model fitting information
broom.mixed::glance(mod_2.3.1) %>%
format_gt_tbl()| nobs | sigma | logLik | AIC | BIC | REMLcrit | df.residual |
|---|---|---|---|---|---|---|
| 48 | 0.333 | −25.33 | 86.65 | 120.33 | 50.65 | 30 |
Tests of Fixed Effects
Since both et and position are categorical variables, we will skip inspecting the model output and focus on the tests of fixed effects. The interaction is not significant which leads us to conclude that the effect of position does not depend on et. However, there is a statistically significant main effect for position (p = 0.0345). There is evidence that the mean resistance varies by position. Finally, the effect of et is not significant. There is no evidence that the mean resistance differs across et levels.
anova(mod_2.3.1) %>%
broom::tidy() %>%
format_gt_tbl()| term | sumsq | meansq | NumDF | DenDF | statistic | p.value |
|---|---|---|---|---|---|---|
| et | 0.65 | 0.216 | 3 | 8 | 1.94 | 0.202 |
| position | 1.13 | 0.376 | 3 | 24 | 3.39 | 0.035 |
| et:position | 0.81 | 0.090 | 9 | 24 | 0.81 | 0.613 |
Least Squares Means by et
We can display the least squares means by et using the emmeans() function from the emmeans package.
# LS means for et
emmeans(mod_2.3.1, ~ et) %>%
broom::tidy() %>%
format_gt_tbl()| et | estimate | std.error | df | statistic | p.value |
|---|---|---|---|---|---|
| 1 | 5.63 | 0.211 | 8 | 26.66 | < 0.001 |
| 2 | 5.97 | 0.211 | 8 | 28.27 | < 0.001 |
| 3 | 6.09 | 0.211 | 8 | 28.85 | < 0.001 |
| 4 | 6.33 | 0.211 | 8 | 30.01 | < 0.001 |
Least Squares Means by position
Similarly, we can display the least squares means by position. Inspection suggests that the mean resistance at position 3 is lower compared to the other positions, which is consistent with the significant main effect of position.
# LS means for position
emmeans(mod_2.3.1, ~ position) %>%
broom::tidy() %>%
format_gt_tbl()| position | estimate | std.error | df | statistic | p.value |
|---|---|---|---|---|---|
| 1 | 6.02 | 0.134 | 18.68 | 44.78 | < 0.001 |
| 2 | 6.13 | 0.134 | 18.68 | 45.62 | < 0.001 |
| 3 | 5.75 | 0.134 | 18.68 | 42.75 | < 0.001 |
| 4 | 6.11 | 0.134 | 18.68 | 45.44 | < 0.001 |
To list all of the pairwise comparisons by position we can pipe the emmeans object to the pairs() function. The pairs() function can take any number of adjustment methods for multiple comparisons. Here I’m showing the results from unadjusted tests as they would be comparable to those obtained in SAS with the default settings. Two pairwise comparisons are statistically significant: position 2 vs. position 3, where position 2 has significantly higher mean resistance than position 3, and position 3 vs. position 4, where position 3 has significantly lower mean resistance than position 4. In both cases, the general interpretation is that position 3 shows lower mean resistance relative to other positions.
emmeans(mod_2.3.1, ~ position) %>%
pairs(adjust = "none") %>%
broom::tidy() %>%
select(-term, -null.value) %>%
format_gt_tbl()| contrast | estimate | std.error | df | statistic | p.value |
|---|---|---|---|---|---|
| position1 - position2 | −0.11 | 0.136 | 24 | −0.83 | 0.413 |
| position1 - position3 | 0.27 | 0.136 | 24 | 2.01 | 0.056 |
| position1 - position4 | −0.09 | 0.136 | 24 | −0.65 | 0.522 |
| position2 - position3 | 0.39 | 0.136 | 24 | 2.84 | 0.009 |
| position2 - position4 | 0.03 | 0.136 | 24 | 0.18 | 0.856 |
| position3 - position4 | −0.36 | 0.136 | 24 | −2.66 | 0.014 |
Least Squares Means et * position
To display all cell means, we modify the formula in the emmeans() function.
# Cell means
emmeans(mod_2.3.1, ~ position:et) %>%
broom::tidy() %>%
format_gt_tbl()| position | et | estimate | std.error | df | statistic | p.value |
|---|---|---|---|---|---|---|
| 1 | 1 | 5.61 | 0.269 | 18.68 | 20.87 | < 0.001 |
| 2 | 1 | 5.45 | 0.269 | 18.68 | 20.27 | < 0.001 |
| 3 | 1 | 5.55 | 0.269 | 18.68 | 20.65 | < 0.001 |
| 4 | 1 | 5.89 | 0.269 | 18.68 | 21.89 | < 0.001 |
| 1 | 2 | 5.99 | 0.269 | 18.68 | 22.29 | < 0.001 |
| 2 | 2 | 6.19 | 0.269 | 18.68 | 23.01 | < 0.001 |
| 3 | 2 | 5.77 | 0.269 | 18.68 | 21.44 | < 0.001 |
| 4 | 2 | 5.92 | 0.269 | 18.68 | 22.00 | < 0.001 |
| 1 | 3 | 6.14 | 0.269 | 18.68 | 22.82 | < 0.001 |
| 2 | 3 | 6.35 | 0.269 | 18.68 | 23.60 | < 0.001 |
| 3 | 3 | 5.77 | 0.269 | 18.68 | 21.47 | < 0.001 |
| 4 | 3 | 6.09 | 0.269 | 18.68 | 22.66 | < 0.001 |
| 1 | 4 | 6.34 | 0.269 | 18.68 | 23.58 | < 0.001 |
| 2 | 4 | 6.55 | 0.269 | 18.68 | 24.37 | < 0.001 |
| 3 | 4 | 5.90 | 0.269 | 18.68 | 21.93 | < 0.001 |
| 4 | 4 | 6.54 | 0.269 | 18.68 | 24.32 | < 0.001 |
Estimates & Contrasts
# 1) et1 vs et2
# Get marginal means for et
emm <- emmeans(mod_2.3.1, ~ et)
# Contrast et1 vs et2
contrast(emm, method = list("et1 vs et2" = c(1, -1, 0, 0))) %>%
broom::tidy() %>%
select(-term, -null.value) %>%
format_gt_tbl()| contrast | estimate | std.error | df | statistic | p.value |
|---|---|---|---|---|---|
| et1 vs et2 | −0.34 | 0.298 | 8 | −1.14 | 0.288 |
# 2) Position1 vs Position 2
emm <- emmeans(mod_2.3.1, ~ position)
# Contrast pos1 vs pos2
contrast(emm, method = list("pos1 vs pos2" = c(1, -1, 0, 0))) %>%
broom::tidy() %>%
select(-term, -null.value) %>%
format_gt_tbl()| contrast | estimate | std.error | df | statistic | p.value |
|---|---|---|---|---|---|
| pos1 vs pos2 | −0.11 | 0.136 | 24 | −0.83 | 0.413 |
# 3) Simple effect contrasts pos1 vs pos2 at each level of et
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") %>%
as_tibble() %>%
filter(et == 1) %>%
format_gt_tbl()| contrast | et | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| pos1 vs pos2 by et | 1 | 0.163 | 0.272 | 24 | 0.600 | 0.554 |
#4 ) Simple effect contrasts
emm <- emmeans(mod_2.3.1, ~ et | position)
# Contrast et1 vs et2 within each level of position
contrast(emm,
method = list("et1 vs et2 by pos" = c(1, -1, 0, 0)),
by = "position",
infer = TRUE) %>%
as_tibble() %>%
filter(position == 1) %>%
format_gt_tbl()| contrast | position | estimate | SE | df | lower.CL | upper.CL | t.ratio | p.value |
|---|---|---|---|---|---|---|---|---|
| et1 vs et2 by pos | 1 | −0.38 | 0.380 | 18.68 | −1.18 | 0.417 | −1.00 | 0.330 |
# 5) The mean of et1 and et2 vs the mean of et3 and et4 at position 2
contrast(emm,
method = list("et1,et2 vs et3,4 by pos" = c(.5, .5, -.5, -.5)),
by = "position") %>%
as_tibble() %>%
filter(position == 2) %>%
format_gt_tbl()| contrast | position | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| et1,et2 vs et3,4 by pos | 2 | −0.63 | 0.269 | 18.68 | −2.35 | 0.030 |
Simple Effect Estimation and Testing
## Simple effect estimation and testing
# Equivalent to the slice option in SAS
# Instead of typing out all possible pairwise differences at each level of pos,
# the pairs() function can print all pairwise differences of of et within each
# level of position. Since this table can be long, the full output is not
# displayed
# pairs(emm, adjust = "none")
# To follow along in the textbook, the results here are filtered to
# position == 2.
pairs(emm, adjust = "none") %>%
as_tibble() %>%
filter(position == 2) %>%
format_gt_tbl()| contrast | position | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| et1 - et2 | 2 | −0.74 | 0.380 | 18.68 | −1.94 | 0.068 |
| et1 - et3 | 2 | −0.90 | 0.380 | 18.68 | −2.36 | 0.029 |
| et1 - et4 | 2 | −1.10 | 0.380 | 18.68 | −2.90 | 0.009 |
| et2 - et3 | 2 | −0.16 | 0.380 | 18.68 | −0.42 | 0.679 |
| et2 - et4 | 2 | −0.37 | 0.380 | 18.68 | −0.96 | 0.347 |
| et3 - et4 | 2 | −0.21 | 0.380 | 18.68 | −0.54 | 0.593 |
Split-Plot Experiments with Whole Plots in Randomized Blocks
cult <- Cultivation %>% as_tibble()
names(cult) <- names(cult) %>% tolower()
cult# A tibble: 24 × 4
block cult inoc drywt
<fct> <fct> <fct> <dbl>
1 1 a con 27.4
2 1 a dea 29.7
3 1 a liv 34.5
4 1 b con 29.4
5 1 b dea 32.5
6 1 b liv 34.4
7 2 a con 28.9
8 2 a dea 28.7
9 2 a liv 33.4
10 2 b con 28.7
# ℹ 14 more rows# Because cultivar is nested within block
cult_mod_1 <- lmerTest::lmer(drywt ~ cult + inoc + cult:inoc + (1 | block + block:cult), data = cult)
cult_mod_1 %>% broom::glance()# A tibble: 1 × 7
nobs sigma logLik AIC BIC REMLcrit df.residual
<int> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
1 24 0.840 -32.5 83.1 93.7 65.1 15cult_mod_1 %>%
anova()Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
cult 0.537 0.537 1 3 0.7616 0.4471
inoc 118.176 59.088 2 12 83.7631 8.919e-08 ***
cult:inoc 1.826 0.913 2 12 1.2942 0.3098
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cult_mod_2 <- lmerTest::lmer(drywt ~ cult + inoc + cult:inoc + (1 | block), data = cult)
cult_mod_2 %>% broom::glance()# A tibble: 1 × 7
nobs sigma logLik AIC BIC REMLcrit df.residual
<int> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
1 24 1.09 -34.2 84.5 93.9 68.5 16anova(cult_mod_1, cult_mod_2)refitting model(s) with ML (instead of REML)Data: cult
Models:
cult_mod_2: drywt ~ cult + inoc + cult:inoc + (1 | block)
cult_mod_1: drywt ~ cult + inoc + cult:inoc + (1 | block + block:cult)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
cult_mod_2 8 89.322 98.746 -36.661 73.322
cult_mod_1 9 86.755 97.358 -34.378 68.755 4.5664 1 0.0326 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cult_mod_3 <- lm(drywt ~ cult + inoc + cult:inoc, data = cult)
anova(cult_mod_2, cult_mod_3)refitting model(s) with ML (instead of REML)Data: cult
Models:
cult_mod_3: drywt ~ cult + inoc + cult:inoc
cult_mod_2: drywt ~ cult + inoc + cult:inoc + (1 | block)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
cult_mod_3 7 96.252 104.498 -41.126 82.252
cult_mod_2 8 89.322 98.746 -36.661 73.322 8.9302 1 0.002805 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Unbalanced Split-Plot Experiments
unb_cult <- cult %>% filter(!(cult == "a" & block == 1))
lmerTest::lmer(drywt ~ cult + inoc + cult:inoc + (1 | block + block:cult), data = unb_cult) %>% anova()Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
cult 0.260 0.260 1 2.2825 0.4211 0.5756
inoc 90.282 45.141 2 10.0000 73.0110 1.082e-06 ***
cult:inoc 2.177 1.089 2 10.0000 1.7607 0.2213
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Multilocation Trial
multi <- Multilocation %>% as_tibble()
names(multi) <- names(multi) %>% tolower()
# Block is nested within location
lmerTest::lmer(adj ~ location + trt + location:trt + (1 | block/location), data = multi) %>% anova()boundary (singular) fit: see help('isSingular')Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
location 6.9475 0.86844 8 18 25.1148 3.001e-08 ***
trt 1.2217 0.40725 3 54 11.7774 4.803e-06 ***
location:trt 0.9966 0.04152 24 54 1.2008 0.2828
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Same model specifies a random intercept for block, and one for the interaction between block and location
lmerTest::lmer(adj ~ location + trt + location:trt + (1 | block) + (1 | block:location), data = multi) %>% anova() boundary (singular) fit: see help('isSingular')Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
location 6.9475 0.86844 8 18 25.1148 3.001e-08 ***
trt 1.2217 0.40725 3 54 11.7774 4.803e-06 ***
location:trt 0.9966 0.04152 24 54 1.2008 0.2828
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Chapter 3 - Analysis of Repeated Measures
weights <- Weights %>% as_tibble()
names(weights) <- names(weights) %>% tolower()
weights <- weights %>% mutate(time = factor(time))# Matches SASfmm compound symmetry approach, which seems to be the default
lmerTest::lmer(
strength ~ program + time + program:time + (1 | subject/program),
data = weights ) %>%
anova()boundary (singular) fit: see help('isSingular')Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
program 7.337 3.6686 2 54 3.0651 0.054842 .
time 53.354 8.8924 6 324 7.4297 1.821e-07 ***
program:time 43.000 3.5834 12 324 2.9939 0.000544 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Since subj is the unique subjectID use that to specify the model to avoid using the nested/interaction syntax
lmerTest::lmer(
strength ~ program + time + program:time + (1 | subj),
data = weights ) %>%
anova()Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
program 7.337 3.6686 2 54 3.0651 0.054842 .
time 53.354 8.8924 6 324 7.4297 1.821e-07 ***
program:time 43.000 3.5834 12 324 2.9939 0.000544 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# To specify compound symmetry explicitly
nlme::lme(
strength ~ program * time,
random = ~ 1 | subj,
correlation = nlme::corCompSymm(form = ~ 1 | subj),
data = weights) %>%
anova() numDF denDF F-value p-value
(Intercept) 1 324 38242.27 <.0001
program 2 54 3.07 0.0548
time 6 324 7.37 <.0001
program:time 12 324 2.99 0.0005# To specify AR(1)
nlme::lme(
strength ~ program * time,
random = ~ 1 | subj,
correlation = nlme::corAR1(form = ~ 1 | subj),
data = weights) %>%
anova() numDF denDF F-value p-value
(Intercept) 1 324 38920.32 <.0001
program 2 54 3.18 0.0493
time 6 324 4.57 0.0002
program:time 12 324 1.32 0.2067# Specify unstructured covariance
weights <- weights %>%
mutate(time_fct = time,
time_num = as.numeric(time))
nlme::lme(
strength ~ program * time,
random = ~ 1 | subj,
correlation = nlme::corSymm(form = ~ time_num | subj),
weights = nlme::varIdent(form = ~ 1 | time_fct),
data = weights,
control = nlme::lmeControl(opt = "nlminb", msMaxIter = 200, msMaxEval = 200)) %>%
anova() numDF denDF F-value p-value
(Intercept) 1 324 47711.26 <.0001
program 2 54 2.77 0.0715
time 6 324 6.99 <.0001
program:time 12 324 1.57 0.0998Chapter 4 - Random Effects Models
Chapter 5 - ANCOVA
Chapter 6 - BLUPs
Chapter 7 - Random Coefficient Models
Chapter 8 - Heterogeneous Variance Models
Chapter 9 - Spatial Variability
Chapter 10 - Case Studies
References
Littell, R. C., Milliken, G. A., Stroup, W. W., & Wolfinger, R. D. (1996). SAS system for mixed models. SAS Institute Inc.