Chapter 14 Linear mixed models

14.1 Random effects

Researchers often collect data in batches, for example

  1. Researchers interested in the effects of insectivorous birds on tree seedling performance in a forest stake out ten 1 m\(^2\) plots and use a wire-mesh cage to cover half of each plot.6 The cage allows insect herbivores into the seedlings inside but excludes insectivorous birds that eat the insects from the seedlings. In every plot, five seedlings are planted within the exclosure and five outside of the exclosure. At the end of the experiment, the total leaf mass is measured on all seedlings. Small, uncontrolled, environmental factors (including soil factors and density of insectivorous birds) will differ between plots but will be common to all seedlings within a plot and we would expect a common response to this uncontrolled variation on top of the differential response to each treatment. As a consequence, the ten measures of leaf mass within a plot are not independent.
  2. A nutrition researcher wants to compare the effect of glucose vs. fructose on glucose metabolism in humans. Ten individuals are recruited. Each individual has blood insulin measured 60 minutes after a noon meal over six successive days. The meal alternates between high glucose and high fructose on each day. Each individual has three measures under high glucose treatment and three measures under high fructose treatment. Small, uncontrolled, environmental factors (including metabolic variation, other meals, activity levels) will differ between the individuals but be common within an individual and we would expect a common response to this uncontrolled variation on top of the differential response to each treatment. As a consequence, the six measures of insulin within an individual are not independent.
  3. A researcher is using a mouse model to compare growth of a wildtype and engineered mutant strain of Staphylococcus. A small spot on both right and left forelimbs of ten mice is shaved and abraded. The two strains are randomly assigned to a side (so each mouse is infected with each strain). Small, uncontrolled, environmental factors (including immune responses) will differ between the mice but be common between the two limbs within a mouse and we would expect a common response to this uncontrolled variation on top of the differential response to each treatment. As a consequence, the two measures of growth within a mouse are not independent.
  4. An ecologist wants to measure the effect of an invasive plant on the reproduction of a native plant. They stake-out ten 2 m\(^2\) plots in a forest and divide each plot into four quadrants, with each quadrant assigned a different treatment: control, activated carbon (a procedural control), extract from the invasive plant’s leaves, and both activated carbon and extract from the invasive plant’s leaves. The response is seedling count. Small, uncontrolled, environmental factors (including soil, drainage, and light) will differ between plots but will will be common to all four quadrants within a plot and we would expect a common response to this uncontrolled variation on top of the differential response to each treatment. As a consequence, the four sets of counts within a plot are not independent.
  5. An ecologist wants to measure the effect of habitat on chick growth in a bird. Five individuals nest in artifical nest boxes built on the boundary between the forest and a large, agricultural field. Five other individuals nest in boxes built deep in the interior of the forest. Chicks in each nest are weighed 13 days after hatching. Small, uncontrolled, environmental factors (including parenting, food availablity, temperature, etc.) will differ between the nests but be common within the nests and we would expect a common response to this uncontrolled variation on top of the differential response to each treatment. As a consequence, the measures of weight within a nest are not independent.
  6. A physiologists has skeletal muscle cells growing in 5 control cultures, and 5 treated cultures. The \(Y\) variable is cell diameter, which is measured in 10 cells per culture. Small, uncontrolled, environmental factors (including chemical) will differ between cultures but will be common to all cells within a culture and we would expect a common response to this uncontrolled variation on top of the differential response to each treatment. As a consequence, the ten measures of diameter within a culture are not independent.
  7. A behavioral biologist wants to measure the effect of a predator fish on the preferred feeding location (open water or vegetation) of a prey fish. Ten tanks are set up with equal amounts of vegetated and unvegetated area. One-third of each tank is screened off to house a predator fish, which are added to five of the tanks. Ten prey fish are added to each tank. The response is minutes spent foraging in the open water as a fraction of total time foraging, which is measured in each fish in each tank. Small, uncontrolled, environmental factors (including temperature, water chemistry, light, and fish behavior) will differ between the tanks but be common within tanks and we would expect a common response to this uncontrolled variation on top of the differential response to each treatment. As a consequence, the ten measures of foraging of each fish within a tank are not independent.

The batches – plots in experiment 1, individuals in experiment 2, mice in experiment 3, plots in experiment 4, nests in experiment 5, cultures in experiment 6, and tanks in experiment 7 – are often referred to as blocks or clusters. I’ll generally use the term “block” in this book. The blocks are the experimental units, meaning that it is at this level that the experimenter is controlling the conditions. There is systematic variation at two levels: among treatments due to treatment effects and among blocks due to block effects. This among-block variation is the random effect. An assumption of modeling random effects is that the blocks (plots/individuals/mice/nests/cultures/tanks) are a random sample of the blocks that could have been sampled. This is often not strictly true as blocks are often convenience samples.

The multiple measures within a block are often called repeated measures, especially if the block is an experimental animal such as a mouse or human. If multiple measures within a treatment level within a block (that is, within a \(block \times treatment\) combination) are taken over time, the data are longitudinal. Sometimes in cell biology, the subsamples within a treatment within a block are called “replicates”, as they are replicates of this \(block \times treatment\) combination, but this can be confusing because the treatments are replicated at the level of the block and not at the level of the subsamples within a treatment by block combination. The blocks are the independent experimental units. Instead the multiple measures of the response within a \(block \times treatment\) combination are subsamples.

Experiments 1 and 2 are examples of a complete randomized block with subsampling design. “Complete” means that each block has all treatment levels or combinations of levels if there is more than one factor. Experiments 3 and 4 are examples of a complete randomized block design. The blocks are complete but there is only one measure of the response per treatment. Experiments 5, 6, and 7 are examples of an incomplete randomized blocks design. The blocks are incomplete because they do not contain less than all treatment levels and combinations. In these examples, each block contains only one treatment level.

14.2 Random effects in statistical models

In all of the above examples, the researcher is interested in the treatment effect but not the variation due to differences among the blocks. The blocks are nuissance factors that add additional variance to the response, with the consequence that estimates of treatment effects are less precise, unless the varriance due to the blocks is explicitly modeled. Including block structure in the design and in the statistical model is known as blocking. A natural way to think about the block factor is as a random effect, meaning that plots in experiment 1 or the mice in experiment 3 are simply random samples from a population of plots or mice. Modeling this using the residual-error specification looks like

\[\begin{equation} y_{ij} = (\beta_{0} + \beta_{0j}) + (\beta_{1} + \beta_{1j}) x_i + \varepsilon_i \tag{14.1} \end{equation}\]

where \(i\) indexes the observation and \(j\) indexes the block (culture, plot, mouse, etc). The intercept parameter \(\beta_{0j}\) is a random intercept and the slope parameter \(\beta_{1j}\) is a random slope. The random intercept has a fixed component (\(\beta_0\)) that is common to all observations and a random component (\(\beta_{0j}\)) that is common within a block but differs among blocks (see table below). In the above equation, I’ve used parentheses to show how these combine into the random intercept that is unique for each block. Similarly, the random slope (treatment effect) has a fixed part (\(\beta_1\)) that is common to all observations and a random component (\(\beta_{1j}\)) that is common within a block but differs among blocks (see table below). Again, these are collected within a pair of parentheses in the equation above.

Table 14.1: The linear mixed model specified above estimates a fixed intercept and fixed slope (treatment effect) that are common to all observations and a random intercept and random slope for each block, each of which is common among observations within a block but differ among observations in different blocks.
block \(b_0\) \(b_{0j}\) \(b_1\) \(b_{1j}\)
1 \(b_0\) \(b_{0,j=1}\) \(b_1\) \(b_{1,j=1}\)
2 \(b_0\) \(b_{0,j=2}\) \(b_1\) \(b_{1,j=2}\)
3 \(b_0\) \(b_{0,j=3}\) \(b_1\) \(b_{1,j=3}\)
4 \(b_0\) \(b_{0,j=4}\) \(b_1\) \(b_{1,j=4}\)
5 \(b_0\) \(b_{0,j=5}\) \(b_1\) \(b_{1,j=5}\)
6 \(b_0\) \(b_{0,j=6}\) \(b_1\) \(b_{1,j=6}\)

Linear mixed models are called “mixed models” because they are a mix of fixed and random components. Another useful way to specify this model is to think about it hierarchically, using

\[\begin{align} y_{ij} &= \beta_{0j} + \beta_{1j}x_i + \varepsilon_i \\ \varepsilon_i &\sim N(0, \sigma) \\ \beta_{0j} &= \beta_{0} + N(0, \sigma_{0}) \\ \beta_{1j} &= \beta_{1} + N(0, \sigma_{1}) \tag{14.2} \end{align}\]

The first line states that the response is a function of a block-specific intercept and a block specific slope plus some error that is unique to each observation. The third and fourth lines state that these block-specific effects are themselves a function of a common effect and a random term that is unique to each block. That is, we have a hierarchical or multi-level structure to the model. Line 1 is the top level and the effects that are specified in line 1 are a function of effects at a second, lower level, which are specified in lines 3 and 4. Because of this structure, linear mixed models are sometimes called hierarchical or multi-level models.

Finally, it’s useful to think how to specify a linear mixed model using the random-draw specification, as this leads naturally to generalized linear mixed models, or GLMMs.

\[\begin{align} y_{ij} &\sim N(\mu_{ij}, \sigma) \\ \mu_{ij} &=\beta_{0j} + \beta_{1j}x_i \\ \beta_{0j} &\sim N(\beta_0, \sigma_0) \\ \beta_{1j} &\sim N(\beta_1, \sigma_1) \tag{14.3} \end{align}\]

14.3 Linear mixed models are flexible

The linear mixed model in Equation (14.1) specifies both a random intercept and a random slope but a researcher might limit the random effect to the intercept only, or less commonly, the slope only. Excluding the random slope from Equation (14.1) results in the model

\[\begin{equation} y_{ij} = (\beta_{0} + \beta_{0j}) + \beta_{1}x_i + \varepsilon_i \tag{14.4} \end{equation}\]

We might use a random-intercept-only model if we think that features of the block would effect the mean response among blocks but not effect the difference in treatment level (or treatment effect) among blocks. For example, differences in the immune systems among the individual mice in experiment 3 might effect growth in both the wild-type and engineered strains of staph but won’t effect the difference in growth between wild-type and engineered strains from one mouse to another.

Not more than you should know – For more complex mixed models, matrix algebra makes the specification of the model much more manageable than the scalar algebra in ??.

\[\begin{equation} \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{Zu} + \boldsymbol{\varepsilon} \end{equation}\]

where \(\mathbf{y}\) is the vector of the response, \(\mathbf{X}\boldsymbol{\beta}\) is the linear predictor of fixed effects and \(\mathbf{Zu}\) is the linear predictor of random effects. \(\mathbf{X}\) is the model matrix for the fixed effects and \(\boldsymbol{\beta}\) is the vector of fixed-effect terms (the fixed part of the intercept (\(\beta_0\)), including the fixed-effect coefficients for each of the

14.4 Visualizing block effects

To visualize random effects due to block, Let’s create fake data that look something like experiments 1 or 2, with a single factor with two treatment levels, \(k=10\) blocks, and \(n=3\) measures for each treatment level within each block. This is a randomized complete block design with subsampling and has a total of \(N=2 \times k \times n\) measures of \(Y\) (and rows of the data.table).

Visualizing random effects. A) The response in the two treatment levels. B) The same data but separated by block. The blue line is at the control mean and the yellow line is at the treated mean. The black dots are the mean response within a block.

Figure 14.1: Visualizing random effects. A) The response in the two treatment levels. B) The same data but separated by block. The blue line is at the control mean and the yellow line is at the treated mean. The black dots are the mean response within a block.

Figure 14.1A shows the response as a function of treatment. The responses are nicely symmetric around the treatment means (the blue and yellow lines). A linear model (and generalized linear models, more generally) assumes that a response, conditional on the \(X\), are independent. Figure 14.1B shows how this assumption is violated for the simulated data. That pattern of residuals within a block around the treatment means does not look at all random. Instead, there is a distinct pattern within a block for the points to cluster either below the treatment means or above the treatment means. In blocks a, b, c, g, and h, all or most of the responses are below their treatment mean (for example in block a, all the blue points are below the blue line and 2 of 3 yellow points are below the yellow line). In blocks, e, f, i, and j, all or most of the responses are above their treatment mean (for example, in block i, all three yellow points are above the yellow line and 2 of 3 blue points are above the blue line). In other words, the responses within a block covary together. For a linear model, this is known as correlated error.

14.5 Linear mixed models can increase precision of point estimates

Block effects are differences in expected mean response among blocks due to unmeasured factors that are shared within blocks but not among blocks. A classical linear model fails to model this component of the total variance in the response, and as a consquence, this block-specific variance becomes part of the error variance. One way to visualize this is by moving the random intercept and random slope components of equation (14.1) to the right and combining it with the observation-specific error

\[\begin{equation} y_{ij} = \beta_{0} + \beta_{1} x_i + (\beta_{0j} + \beta_{1j} + \varepsilon_i) \tag{14.4} \end{equation}\]

% which shows that the random effects \(\beta_{0j}\) and \(\beta_{1j}\) are modeled as error in a linear model. As a consequence, the residual variance is larger and, therefore, the standard errors of point estimates, including means, coefficients of the model, and contrasts from the model, are larger. Here is the table of model coefficients of the data in Figure 14.1 fitted using a classical linear model

##               Estimate Std. Error  t value     Pr(>|t|)
## (Intercept) 10.0235749  0.2439544 41.08791 1.393185e-44
## TreatmentT+  0.6884787  0.3450036  1.99557 5.068466e-02

and a linear mixed model

##               Estimate Std. Error       df   t value     Pr(>|t|)
## (Intercept) 10.0235749  0.3954692 9.000000 25.346031 1.114196e-09
## TreatmentT+  0.6884787  0.2700725 8.999999  2.549237 3.123383e-02

The linear mixed model has increased precision (a smaller SE of the estimates) because it estimates the value of \(\beta_{0j}\) and \(\beta_{1j}\) for each block. The linear model does not estimate these parameters and the variance in these parameters is swept into the residual variance.

NHST blues – A paired t-test is equivalent to the special case of a linear mixed model with a single factor with two treatment levels, \(k\) blocks, and a single measure of each treatment level within each block. A good example is the wild type vs. engineered staph count in mice in experiment 3 above. A linear mixed model is much more flexible than a paired t-test because it allows a researcher to add treatment levels, additional factors, and covariates to the model. In addition, a linear mixed model can handle missing data.

Here is fake data similar in design to experiment 3 with a single factor with two treatment levels and both levels applied to the same experimental unit.

set.seed(2)
n <- 10 # number of mice (blocks)
x <- rep(c("WT","T+"), each=n) # treatments
id <- rep(letters[1:n], 2) # block id
y <- c(rnorm(n, mean=10), rnorm(n, mean=11))
fake_data <- data.table(Y=y, X=x, ID=id)

The t-test p-value is

t.test(Y~X, data=fake_data, paired=TRUE)$p.value
## [1] 0.05336815

and the coefficient table of the fixed effect in the linear mixed model is

coef(summary(lme(Y~X, random = ~1|ID, correlation=corCompSymm(form=~1|ID), data=fake_data)))
##                  Value Std.Error DF   t-value      p-value
## (Intercept) 11.1797704 0.3438775  9 32.510914 1.212113e-10
## XWT         -0.9686188 0.4358740  9 -2.222245 5.336815e-02

14.6 Linear mixed models are used to avoid pseudoreplication

14.7 Linear mixed models shrink coefficients by partial pooling

In experiment 1 above, there are 10 sites (maybe different woodlots). In each plot, five seedlings are planted inside a cage and five outside the cage. The cage excludes insectivorous birds but not herbivorous insects. The researchers are investigating how birds affect plant growth indirectly – by eating insects that feed on the plants. The response is total leaf area in each seedling.

Let’s say we want to know the treatment effect in each of these sites. There are several ways of estimating this.

  1. Fit \(k\) separate models, one for each site. The intercept (control mean) and slope (treatment effect) parameters for each site are estimated independently from all other sites. Conequently, the model parameters are computed using no pooling. For the estimation of the \(\beta\) terms, this is equivalent to a single, factorial linear model with \(Site\) modeled as a fixed effect (this is not true for the estimate of the standard errors of these terms since these are computed from the residual sum of squares of the model. For balanced data, all of the “intercept” or “slope” terms will have the same SE in the factorial analysis but differ among the \(k\) independent analyses).

  2. Fit a linear model to all the data combined as if these were from a single site, and assign the intercept and treatment effect paramters to all sites. The model parameters are computed using complete pooling.

  3. Fit a linear mixed model to all the data, using site as a random factor to estimate both random intercepts and slopes. Similar to the no-pooling analysis, there will be different intercept (control mean) and slope (treatment effect) estimates for each site, but unlike the no-pooling analysis, these estimates are computed by combining information from the other sites. The information used to estimate parameters in a linear mixed model is somewhere in between no pooling and complete pooling and is sometimes called partial pooling.

The consequence of partial pooling in a linear mixed model is that site intercepts (control means) are pulled toward the single intercept in the complete-pooling analysis and the site slopes (treatment effects) are pulled toward the single slope in the complete-pooling analysis. This has the consequence that the differences in parameter estimates among sites are shrunk toward zero. A consequence of this shrinkage is that the variance of the intercept estimates or of the slope estimates is smaller than that in the no-pooling analysis. Figure ?? shows this shrinkage effect using fake data simulating the seedling experiment.

Shrinkage estimates of the treatment effect in a linear mixed model. The grey line is the estimate using complete pooling (so there is only one estimate which is assigned to each site). In general, the partial-pooling (linear mixed model) estimates (yellow) are generally closer to the complete pooling estimate than the no-pooling (separate linear models) estimates (blue). More specifically, if the no-pooling estimate is far from the complete pooling estimate, the partial pooling estimate is *much* closer to the complete pooling estimate. The consequence of partial pooling is that the differences among the estimates are shrunk toward zero.

Figure 14.2: Shrinkage estimates of the treatment effect in a linear mixed model. The grey line is the estimate using complete pooling (so there is only one estimate which is assigned to each site). In general, the partial-pooling (linear mixed model) estimates (yellow) are generally closer to the complete pooling estimate than the no-pooling (separate linear models) estimates (blue). More specifically, if the no-pooling estimate is far from the complete pooling estimate, the partial pooling estimate is much closer to the complete pooling estimate. The consequence of partial pooling is that the differences among the estimates are shrunk toward zero.

The linear mixed model estimates of the treatment effects for each site are a type of shrinkage estimate and a linear mixed model is one kind of shrinkage estimator. Shrinkage estimates have fascinating properties:

  1. the variance of shrinkage estimates is less than that of ordinary least squares estimates (no-pooling, or using the block as a fixed factor)
  2. shrinkage estimates are biased but the OLS estimates are not. This means that the expected value of a coefficient from the linear mixed model does not equal the true (parameter) value! Or, more formally, \(\mathrm{E}(b_j) \ne \beta_j\).
  3. the mean square error of shrinkage estimates will be smaller than that for OLS estimates.

The first property was discussed above and shown in Figure 14.2. The second property raises the question, if we want to estimate the treatment effects within each site, why would we ever want to use \(Site\) as a random instead of fixed effect? The answer is the third property, which can be summarized as, “if we were to replicate the experiment many times, the shrinkage estimates will be, on average, less wrong (or closer to the true value) than the OLS estimates, where”wrong" is the absolute deviation from the true value."

When shrinkage estimators were first discovered, the third property surprised stasticians. The third property has profound consequences. Consider a scenario where researchers want to compare the performance of a new expression vector to that of an existing expression vector on protein production using E. coli. The researchers have ten different E. coli strains and are interested in strain-specific effects because they will choose the three strains with the largest effects for further testing. The researchers measure the response of each strain five times.

Table 14.2: Effect of new expression vector on protein production in ten strains of E. coli using a fixed effect factorial model and linear mixed model.
Strain \(\beta_{1j}\) fixed \(b_{1j}\) random \(b_{1j}\)
a 0.91 1.07 0.98
b 0.87 0.94 0.85
c 0.90 -1.03 0.30
d 0.81 0.64 0.63
e 1.09 1.00 1.07
f 0.62 0.91 1.14
g 1.33 2.26 1.36
h 1.27 1.48 0.96
i 1.61 0.57 1.13
j 0.89 1.50 0.93

The table above shows the true strain-specific effect and both the fixed (OLS) and random (LMM) effect estimates. The largest OLS estimate is 70% larger than the true effect and the strain with the largest true effect is not among the top three biggest OLS estimates (its ranking is 9/10). By contrast, the LMM estimates are closer to the true effects and the top strain is among the three largest LMM estimates.

These results are specific to these fake data but more generally, 1) the largest OLS estimates are inflated (larger error from the true effect), relative to the largest LMM estimates 2) overall, the LMM estimates will be closer than the OLS estimates to the true effects

To understand this, rank order the treatment effects for each strain. An individual strain’s position in this rank is the sum of the true effect for that strain and some random error. Because OLS, relative to shrinkage estimates, have greater variance in the estimate (that is, the random error component is bigger), the biggest effects estimated by OLS are more likely to be big because of the error component, compared to shrinkage estimates.

Not more than you want to know – Shrinkage estimators are not only useful when we are interested in block-specific effects but are also useful for estimating effects when there are multiple responses. For example, consider a researcher interested in measuring the effects of some exercise treatment on gene expression in adipose cells. The researcher measures expression levels in 10,000 genes. Given the typical content in undergraduate biostatics courses, a researcher would probably model these responses using 10,000 t-tests, or equivalently, 10,000 separate linear models. If the tests were ranked by \(p\)-value or absolute effect size, many of the genes with largest absolute effect would be there because of a large error component and many of the largest effects would be massively overinflated. Re-imagining the design as a single, linear mixed model with each gene modeled as a block would lead to a rank order in which the biggest measured effects more closely approximate the true effects.

14.8 Working in R

The major function for working with linear mixed models is lmer() from the lme4 package. An older, and still sometimes used and useful function is lme() from the nlme package. The authors of the lme4 package argue that the df in a linear mixed model are too approximate for a useful \(p\)-value and, consequently, the lme function does not return a \(p\)-value. Many biological researchers want a \(p\)-value and typically use the lmerTest package to get this.

14.8.1 coral data

Source Zill, J. A., Gil, M. A., & Osenberg, C. W. (2017). When environmental factors become stressors: interactive effects of vermetid gastropods and sedimentation on corals. Biology letters, 13(3), 20160957.

Dryad source https://datadryad.org/resource/doi:10.5061/dryad.p59n8

file name “VermetidSedimentData_ZillGilOsenberg_DRYAD.xlsx”

folder <- "Data from When environmental factors become stressors- interactive effects of vermetid gastropods and sedimentation on corals"
fn <- "VermetidSedimentData_ZillGilOsenberg_DRYAD.xlsx"
sheet_i <- "Coral Growth Rate Data"
file_path <- paste(data_path, folder, fn, sep="/")
coral <- data.table(read_excel(file_path, sheet=sheet_i))
setnames(coral, old=colnames(coral), new=clean_label(colnames(coral)))
coral[, Vermetids:=factor(Vermetids)]
coral[, Sediment:=factor(Sediment)]

lmer adds the random component to the formula. lme adds the random component as a separate argument

# to reproduce the results
# observation 2 should be excluded from the analysis
inc <- c(1, 3:nrow(coral))

# specification using lmer
# random intercept only
fit.lmer1 <- lmer(GrowthRate ~ Vermetids*Sediment + (1|Block), data=coral[inc])
# random intercept and slope
fit.lmer2 <- lmer(GrowthRate ~ Vermetids*Sediment + (Vermetids|Block) + (Sediment|Block), data=coral[inc])
# random intercept and slope
fit.lmer3 <- lmer(GrowthRate ~ Vermetids*Sediment + (Vermetids|Block) + (Sediment|Block), data=coral[inc])
# to include the interaction as a random effect we'd need subsampling within each factorial treatment combination

# specification using lme
fit.lme <- lme(GrowthRate ~ Vermetids*Sediment, random= ~1|Block, data=coral[inc])

# results using lmer fit
coefficients(fit.lmer1)
## $Block
##    (Intercept)  Vermetids1 Sediment1 Vermetids1:Sediment1
## 1     1.205030 0.004655556 0.2852037           -0.7735482
## 2     1.336057 0.004655556 0.2852037           -0.7735482
## 3     1.213975 0.004655556 0.2852037           -0.7735482
## 4     1.262806 0.004655556 0.2852037           -0.7735482
## 6     1.320778 0.004655556 0.2852037           -0.7735482
## 7     1.201253 0.004655556 0.2852037           -0.7735482
## 8     1.314433 0.004655556 0.2852037           -0.7735482
## 9     1.199842 0.004655556 0.2852037           -0.7735482
## 10    1.361526 0.004655556 0.2852037           -0.7735482
## 
## attr(,"class")
## [1] "coef.mer"
coefficients(summary(fit.lmer1))
##                          Estimate Std. Error       df     t value
## (Intercept)           1.268411111  0.1541678 30.42810  8.22747196
## Vermetids1            0.004655556  0.2091414 22.94176  0.02226033
## Sediment1             0.285203739  0.2160141 23.53066  1.32030164
## Vermetids1:Sediment1 -0.773548184  0.3006696 23.24572 -2.57275179
##                          Pr(>|t|)
## (Intercept)          3.129452e-09
## Vermetids1           9.824328e-01
## Sediment1            1.994339e-01
## Vermetids1:Sediment1 1.693478e-02
fit.emm <- emmeans(fit.lmer1, specs=c("Vermetids", "Sediment"))
summary(contrast(fit.emm, method="revpairwise", adjust="none"), infer=c(TRUE, TRUE))
##  contrast      estimate        SE    df   lower.CL    upper.CL t.ratio
##  1,0 - 0,0  0.004655556 0.2091414 23.04 -0.4279489  0.43725999   0.022
##  0,1 - 0,0  0.285203739 0.2167702 23.62 -0.1625690  0.73297646   1.316
##  0,1 - 1,0  0.280548184 0.2167702 23.62 -0.1672245  0.72832090   1.294
##  1,1 - 0,0 -0.483688889 0.2091414 23.04 -0.9162933 -0.05108446  -2.313
##  1,1 - 1,0 -0.488344444 0.2091414 23.04 -0.9209489 -0.05574001  -2.335
##  1,1 - 0,1 -0.768892628 0.2167702 23.62 -1.2166653 -0.32111991  -3.547
##  p.value
##   0.9824
##   0.2009
##   0.2081
##   0.0300
##   0.0286
##   0.0017
## 
## Confidence level used: 0.95
fit.lmer1.ml <- lmer(GrowthRate ~ Vermetids*Sediment + (1|Block), data=coral[inc], REML=FALSE)
# random intercept and slope
fit.lmer3.ml <- lmer(GrowthRate ~ Vermetids*Sediment + (Vermetids|Block) + (Sediment|Block), data=coral[inc], REML=FALSE)

AIC(fit.lmer1.ml)
## [1] 52.80349
AIC(fit.lmer3.ml)
## [1] 61.04197

The formula for lmer


  1. Giffard, B., Corcket, E., Barbaro, L., & Jactel, H. (2012). Bird predation enhances tree seedling resistance to insect herbivores in contrasting forest habitats. Oecologia, 168(2), 415-424