Tutorial 9.3a - Randomized Complete Block ANOVA

20 Jul 2018

Overview

When single sampling units are selected amongst highly heterogeneous conditions, it is unlikely that these single units will adequately represent the populations and repeated sampling is likely to yield very different outcomes.

Tutorial 9.2a (Nested ANOVA), introduced the concept of employing sub-replicates that are nested within the main treatment levels as a means of absorbing some of the unexplained variability that would otherwise arise from designs in which sampling units are selected from amongst highly heterogeneous conditions. Such (nested) designs are useful in circumstances where the levels of the main treatment (such as burnt and un-burnt sites) occur at a much larger temporal or spatial scale than the experimental/sampling units (e.g. vegetation monitoring quadrats).

For circumstances in which the main treatments can be applied (or naturally occur) at the same scale as the sampling units (such as whether a stream rock is enclosed by a fish proof fence or not), an alternative design is available. In this design (randomized complete block design), each of the levels of the main treatment factor are grouped (blocked) together (in space and/or time) and therefore, whilst the conditions between the groups (referred to as `blocks') might vary substantially, the conditions under which each of the levels of the treatment are tested within any given block are far more homogeneous (see Figure below).

If any differences between blocks (due to the heterogeneity) can account for some of the total variability between the sampling units (thereby reducing the amount of variability that the main treatment(s) failed to explain), then the main test of treatment effects will be more powerful/sensitive.

As an simple example of a randomized complete block (RCB) design, consider an investigation into the roles of different organism scales (microbial, macro invertebrate and vertebrate) on the breakdown of leaf debris packs within streams. An experiment could consist of four treatment levels - leaf packs protected by fish-proof mesh, leaf packs protected by fine macro invertebrate exclusion mesh, leaf packs protected by dissolving antibacterial tablets, and leaf packs relatively unprotected as controls.

As an acknowledgement that there are many other unmeasured factors that could influence leaf pack breakdown (such as flow velocity, light levels, etc) and that these are likely to vary substantially throughout a stream, the treatments are to be arranged into groups or 'blocks' (each containing a single control, microbial, macro invertebrate and fish protected leaf pack). Blocks of treatment sets are then secured in locations haphazardly selected throughout a particular reach of stream. Importantly, the arrangement of treatments in each block must be randomized to prevent the introduction of some systematic bias - such as light angle, current direction etc.

Blocking does however come at a cost. The blocks absorb both unexplained variability as well as degrees of freedom from the residuals. Consequently, if the amount of the total unexplained variation that is absorbed by the blocks is not sufficiently large enough to offset the reduction in degrees of freedom (which may result from either less than expected heterogeneity, or due to the scale at which the blocks are established being inappropriate to explain much of the variation), for a given number of sampling units (leaf packs), the tests of main treatment effects will suffer power reductions.

Treatments can also be applied sequentially or repeatedly at the scale of the entire block, such that at any single time, only a single treatment level is being applied (see the lower two sub-figures above). Such designs are called repeated measures. A repeated measures ANOVA is to an single factor ANOVA as a paired t-test is to a independent samples t-test.

One example of a repeated measures analysis might be an investigation into the effects of a five different diet drugs (four doses and a placebo) on the food intake of lab rats. Each of the rats (`subjects') is subject to each of the four drugs (within subject effects) which are administered in a random order.

In another example, temporal recovery responses of sharks to bi-catch entanglement stresses might be simulated by analyzing blood samples collected from captive sharks (subjects) every half hour for three hours following a stress inducing restraint. This repeated measures design allows the anticipated variability in stress tolerances between individual sharks to be accounted for in the analysis (so as to permit more powerful test of the main treatments). Furthermore, by performing repeated measures on the same subjects, repeated measures designs reduce the number of subjects required for the investigation.

Essentially, this is a randomized complete block design except that the within subject (block) effect (e.g. time since stress exposure) cannot be randomized (the consequences of which are discussed in section on Sphericity).

To suppress contamination effects resulting from the proximity of treatment sampling units within a block, units should be adequately spaced in time and space. For example, the leaf packs should not be so close to one another that the control packs are effected by the antibacterial tablets and there should be sufficient recovery time between subsequent drug administrations.

In addition, the order or arrangement of treatments within the blocks must be randomized so as to prevent both confounding as well as computational complications (Sphericity). Whilst this is relatively straight forward for the classic randomized complete block design (such as the leaf packs in streams), it is logically not possible for repeated measures designs.

Blocking factors are typically random factors (see section~\ref{chpt:ANOVA.fixedVsRandomFactor}) that represent all the possible blocks that could be selected. As such, no individual block can truly be replicated. Randomized complete block and repeated measures designs can therefore also be thought of as un-replicated factorial designs in which there are two or more factors but that the interactions between the blocks and all the within block factors are not replicated.

Linear models

The linear models for two and three factor nested design are:
$$ \begin{align} y_{ij}&=\mu+\beta_{i}+\alpha_j + \varepsilon_{ij} &\hspace{2em} \varepsilon_{ij} &\sim\mathcal{N}(0,\sigma^2), \hspace{1em}\sum{}{\beta=0}\\ y_{ijk}&=\mu+\beta_{i} + \alpha_j + \gamma_{k} + \beta\alpha_{ij} + \beta\gamma_{ik} + \alpha\gamma_{jk} + \gamma\alpha\beta_{ijk} + \varepsilon_{ijk} \hspace{2em} (Model 1)\\ y_{ijk}&=\mu+\beta_{i} + \alpha_j + \gamma_{k} + \alpha\gamma_{jk} + \varepsilon_{ijk} \hspace{2em}(Model 2) \end{align} $$ where $\mu$ is the overall mean, $\beta$ is the effect of the Blocking Factor B, $\alpha$ and $\gamma$ are the effects of withing block Factor A and Factor C respectively and $\varepsilon$ is the random unexplained or residual component.

Tests for the effects of blocks as well as effects within blocks assume that there are no interactions between blocks and the within block effects. That is, it is assumed that any effects are of similar nature within each of the blocks. Whilst this assumption may well hold for experiments that are able to consciously set the scale over which the blocking units are arranged, when designs utilize arbitrary or naturally occurring blocking units, the magnitude and even polarity of the main effects are likely to vary substantially between the blocks.

The preferred (non-additive or `Model 1') approach to un-replicated factorial analysis of some bio-statisticians is to include the block by within subject effect interactions (e.g. $\beta\alpha$). Whilst these interaction effects cannot be formally tested, they can be used as the denominators in F-ratio calculations of their respective main effects tests (see the tables that follow).

Proponents argue that since these blocking interactions cannot be formally tested, there is no sound inferential basis for using these error terms separately. Alternatively, models can be fitted additively (`Model 2') whereby all the block by within subject effect interactions are pooled into a single residual term ($\varepsilon$). Although the latter approach is simpler, each of the within subject effects tests do assume that there are no interactions involving the blocks and that perhaps even more restrictively, that sphericity (see section Sphericity) holds across the entire design.

Null hypotheses

Separate null hypotheses are associated with each of the factors, however, blocking factors are typically only added to absorb some of the unexplained variability and therefore specific hypothesis tests associated with blocking factors are of lesser biological importance.

Factor A - the main within block treatment effect

Fixed (typical case)

H$_0(A)$: $\mu_1=\mu_2=...=\mu_i=\mu$ (the population group means of A (pooling B) are all equal)

The mean of population 1 (pooling blocks) is equal to that of population 2 and so on, and thus all population means are equal to an overall mean. No effect of A within each block (Model 2) or over and above the effect of blocks. If the effect of the $i^{th}$ group is the difference between the $i^{th}$ group mean and the overall mean ($\alpha_i = \mu_i - \mu$) then the H$_0$ can alternatively be written as:

H$_0(A)$: $\alpha_1 = \alpha_2 = ... = \alpha_i = 0$ (the effect of each group equals zero)
If one or more of the $\alpha_i$ are different from zero (the response mean for this treatment differs from the overall response mean), the null hypothesis is not true indicating that the treatment does affect the response variable.

Random

H$_0(A)$: $\sigma_\alpha^2=0$ (population variance equals zero)

There is no added variance due to all possible levels of A (pooling B).

Factor B - the blocking factor

Random (typical case)

H$_0(B)$: $\sigma_\beta^2=0$ (population variance equals zero)

There is no added variance due to all possible levels of B.

Fixed

H$_0(B)$: $\mu_{1}=\mu_{2}=...=\mu_{i}=\mu$ (the population group means of B are all equal)

H$_0(B)$: $\beta_{1} = \beta_{2}= ... = \beta_{i} = 0$ (the effect of each chosen B group equals zero)

The null hypotheses associated with additional within block factors, are treated similarly to Factor A above.

Analysis of Variance

Partitioning of the total variance sequentially into explained and unexplained components and F-ratio calculations predominantly follows the rules established in tutorials on Nested ANOVA and Factorial ANOVA. Randomized block and repeated measures designs can essentially be analysed as Model III ANOVAs. The appropriate unexplained residuals and therefore the appropriate F-ratios for each factor differ according to the different null hypotheses associated with different combinations of fixed and random factors and what analysis approach (Model 1 or 2) is adopted for the randomized block linear model (see the following tables).

In additively (Model 2) fitted models (in which block interactions are assumed not to exist and are thus pooled into a single residual term), hypothesis tests of the effect of B (blocking factor) are possible. However, since blocking designs are usually employed out of expectation for substantial variability between blocks, such tests are rarely of much biological interest.

         summary(aov(y~A+Error(B), data))

         library(nlme)
         summary(lme(y~A, random=~1|B, data))
         anova(lme(y~A, random=~1|B, data))
         #OR
         library(lme4)
         summary(lmer(y~(1|B)+A, data))

         library(nlme)
         summary(lme(y~1, random=~1|B/A, data))
         #OR
         library(lme4)
         summary(lmer(y~(1|B)+(1|A), data))

         summary(aov(y~Error(B/(A*C))+A*C, data))

         summary(aov(y~Error(B)+A*C, data))
         #or 
         library(nlme)
         summary(lme(y~A*C, random=~1|B, data))

         anova(lme(DV~A*C, random=~1|B), type='marginal')
         anova(lme(DV~A*C, random=~1|B, correlation=corAR1(form=~1|B)), type='marginal')

         library(lmer)
         summary(lmer(DV~(~1|A*C)+(1|B)), data) #for example
         # or
         summary(lm(y~A*B*C, data)) # and make manual F-ratio etc adjustments

         library(nlme)
         summary(lme(y~1, random=~1|B/(A*C), data))

Assumptions

As with other ANOVA designs, the reliability of hypothesis tests is dependent on the residuals being:

• normally distributed. Boxplots using the appropriate scale of replication (reflecting the appropriate residuals/F-ratio denominator (see Tables above) should be used to explore normality. Scale transformations are often useful.
• equally varied. Boxplots and plots of means against variance (using the appropriate scale of replication) should be used to explore the spread of values. Residual plots should reveal no patterns. Scale transformations are often useful.
• independent of one another. Although the observations within a block may not strictly be independent, provided the treatments are applied or ordered randomly within each block or subject, within block proximity effects on the residuals should be random across all blocks and thus the residuals should still be independent of one another. Nevertheless, it is important that experimental units within blocks are adequately spaced in space and time so as to suppress contamination or carryover effects.

Sphericity

Un-replicated factorial designs extend the usual equal variance (no relationship between mean and variance) assumption to further assume that the differences between each pair of within block treatments are equally varied across the blocks (see the following figure). To meet this assumption a matrix of variances (between pairs of observations within treatments) and covariances (between treatment pairs within each block) must display a pattern known as sphericity. Strickly, the variance-covariance matrix must display a very specific pattern of sphericity in which both variances and covariances are equal (compound symmetry), however an F-ratio will still reliably follow an F distribution provided basic sphericity holds.

Single factor ANOVA
		Variance-covariance structure T1 T2 T3 T4 T1 0.15 0.00 0.00 0.00 T2 0.00 0.15 0.00 0.00 T3 0.00 0.00 0.15 0.00 T4 0.00 0.00 0.00 0.15
Randomized Complete Block
		Variance-covariance structure T1 T2 T3 T4 T1 0.15 0.05 0.05 0.05 T2 0.05 0.15 0.05 0.05 T3 0.05 0.05 0.15 0.05 T4 0.05 0.05 0.05 0.15
Repeated Measures Design
		Variance-covariance structure T1 T2 T3 T4 T1 0.50 0.60 0.30 0.10 T2 0.60 0.50 0.60 0.30 T3 0.30 0.60 0.50 0.60 T4 0.10 0.30 0.60 0.50

Typically, un-replicated factorial designs in which the treatment levels have been randomly arranged (temporally and spatially) within each block (randomized complete block) should meet this sphericity assumption. Conversely, repeated measures designs that incorporate factors whose levels cannot be randomized within each block (such as distances from a source or time), are likely to violate this assumption. In such designs, the differences between treatments that are arranged closer together (in either space or time) are likely to be less variable (greater paired covariances) than the differences between treatments that are further apart.

Hypothesis tests are not very robust to substantial deviations from sphericity and consequently tend to have inflated type I errors. There are three broad techniques for compensating or tackling the issues of sphericity:

• reducing the degrees of freedom for F tests according to the degree of departure from sphericity (measured by epsilon ($\epsilon$)). The two main estimates of epsilon are Greenhouse-Geisser and Huynh-Feldt, the former of which is preferred (as it provides more liberal protection) unless its value is less than 0.5.
• perform a multivariate ANOVA (MANOVA). Although the sphericity assumption does not apply to such procedures, MANOVA's essentially test null hypotheses about the differences between multiple treatment pairs (and thus test whether an array of population means equals zero), and therefore assume multivariate normality - a difficult assumption to explore.
• fit a linear mixed effects (lme) model and incorporate a particular correlation structure into the modelling. The approximate form of the correlation structure can be specified up-front when fitting linear mixed effects models (via lme) and thus correlated data are more appropriately handled. A selection of variance-covariance structures appropriate for biological data are listed in Table~\ref{tab:randomizedBlockAnova-correlationStructures}. It is generally recommended that linear mixed effects models be fitted with a range of covariance structures. The "best" covariance structure is that the results in a better fit (as measured by either AIC, BIC or ANOVA) than a model fitted with a compound symmetry structure.

Block by treatment interactions

The presence of block by treatment interactions have important implications for models that incorporate a single within block factor as well as additive models involving two or more within block factors. In both cases, the blocking interactions and overall random errors are pooled into a residual term that is used as the denominator in F-ratio calculations (see the tables above).

Consequently, block by treatment interactions increase the denominator ($MS_{Resid}$) resulting in lower F-ratios (lower power). Moreover, the presence of strong blocking interactions would imply that any effects of the main factor are not consistent. Drawing conclusions from such an analysis (particularly in light of non-significant main effects) is difficult. Unless we assume that there are no block by within block interactions, non-significant within block effects could be due to either an absence of a treatment effect, or as a result of opposing effects within different blocks. As these block by within block interactions are unreplicated, they can neither be formally tested nor is it possible to perform main effects tests to diagnose non-significant within block effects.

Block by treatment interactions can be diagnosed by examining;

• interaction (cell means) plot. The mean ($n=1$) response for each level of the main factor is plotted against the block number. Parallel lines infer no block by treatment interaction.
• residual plot. A curvilinear pattern in which the residual values switch from positive to negative and back again (or visa versa) over the range of predicted values implies that the scale (magnitude but not polarity) of the main treatment effects differs substantially across the range of blocks. Scale transformations can be useful in removing such interactions.
• Tukey's test for non-additivity evaluated at $\alpha=0.10$ or even $\alpha=0.25$. This (curvilinear test) formally tests for the presence of a quadratic trend in the relationship between residuals and predicted values. As such, it too is only appropriate for simple interactions of scale.

$R^2$ approximations

Whilst $R^2$ is a popular goodness of fit metric in simple linear models, its use is rarely extended to (generalized) linear mixed effects models. The reasons for this include:

• there are numerous ways that $R^2$ could be defined for mixed effects models, some of which can result in values that are either difficult to interpret or illogical (for example negative $R^2$).
• perhaps as a consequence, software implementation is also largely lacking.

Nakagawa and Schielzeth (2013) discuss the issues associated with $R^2$ calculations and suggest a series of simple calculations to yield sensible $R^2$ values from mixed effects models.

An $R^2$ value quantifies the proportion of variance explained by a model (or by terms in a model) - the higher the value, the better the model (or term) fit. Nakagawa and Schielzeth (2013) offered up two $R^2$ for mixed effects models:

• Marginal $R^2$ - the proportion of total variance explained by the fixed effects. $$ \text{Marginal}~R^2 = \frac{\sigma^2_f}{\sigma^2_f + \sum^z_l{\sigma^2_l} + \sigma^2_d + \sigma^2_e} $$ where $\sigma^2_f$ is the variance of the fitted values (i.e. $\sigma^2_f = var(\mathbf{X\beta})$) on the link scale, $\sum^z_l{\sigma^2_l}$ is the sum of the $z$ random effects (including the residuals) and $\sigma^2_d$ and $\sigma^2_e$ are additional variance components appropriate when using non-Gaussian distributions.
• Conditional $R^2$ - the proportion of the total variance collectively explained by the fixed and random factors $$ \text{Conditional}~R^2 = \frac{\sigma^2_f + \sum^z_l{\sigma^2_l}}{\sigma^2_f + \sum^z_l{\sigma^2_l} + \sigma^2_d + \sigma^2_e} $$

ANOVA in R

Simple RCB

Scenario and Data

Imagine we has designed an experiment in which we intend to measure a response ($y$) to one of treatments (three levels; 'a1', 'a2' and 'a3'). Unfortunately, the system that we intend to sample is spatially heterogeneous and thus will add a great deal of noise to the data that will make it difficult to detect a signal (impact of treatment).

Thus in an attempt to constrain this variability you decide to apply a design (RCB) in which each of the treatments within each of 35 blocks dispersed randomly throughout the landscape. As this section is mainly about the generation of artificial data (and not specifically about what to do with the data), understanding the actual details are optional and can be safely skipped. Consequently, I have folded (toggled) this section away.

library(plyr)
set.seed(1)
nTreat <- 3
nBlock <- 35
sigma <- 5
sigma.block <- 12
n <- nBlock * nTreat
Block <- gl(nBlock, k = 1)
A <- gl(nTreat, k = 1, lab = LETTERS[1:nTreat])
dt <- expand.grid(A = A, Block = Block)
# Xmat <- model.matrix(~Block + A + Block:A, data=dt)
Xmat <- model.matrix(~-1 + Block + A, data = dt)
block.effects <- rnorm(n = nBlock, mean = 40, sd = sigma.block)
A.effects <- c(30, 40)
all.effects <- c(block.effects, A.effects)
lin.pred <- Xmat %*% all.effects

# OR
Xmat <- cbind(model.matrix(~-1 + Block, data = dt), model.matrix(~-1 + A, data = dt))
## Sum to zero block effects
block.effects <- rnorm(n = nBlock, mean = 0, sd = sigma.block)
A.effects <- c(40, 70, 80)
all.effects <- c(block.effects, A.effects)
lin.pred <- Xmat %*% all.effects



## the quadrat observations (within sites) are drawn from normal distributions with means according to
## the site means and standard deviations of 5
y <- rnorm(n, lin.pred, sigma)
data.rcb1 <- data.frame(y = y, expand.grid(A = A, Block = paste0("B", Block)))
head(data.rcb1)  #print out the first six rows of the data set

         y A Block
1 37.39761 A    B1
2 61.47033 B    B1
3 78.07370 C    B1
4 30.59803 A    B2
5 59.00035 B    B2
6 76.72575 C    B2

write.table(data.rcb1, file = "../downloads/data/data.rcb1.csv", quote = FALSE, row.names = FALSE, sep = ",")

Exploratory data analysis

Normality and Homogeneity of variance

boxplot(y ~ A, data.rcb1)

Conclusions:

• there is no evidence that the response variable is consistently non-normal across all populations - each boxplot is approximately symmetrical
• there is no evidence that variance (as estimated by the height of the boxplots) differs between the five populations. . More importantly, there is no evidence of a relationship between mean and variance - the height of boxplots does not increase with increasing position along the y-axis. Hence it there is no evidence of non-homogeneity

• transform the scale of the response variables (to address normality etc). Note transformations should be applied to the entire response variable (not just those populations that are skewed).

Block by within-Block interaction

library(car)
with(data.rcb1, interaction.plot(A, Block, y))

# OR with ggplot
library(ggplot2)
ggplot(data.rcb1, aes(y = y, x = A, group = Block, color = Block)) +
    geom_line() + guides(color = guide_legend(ncol = 3))

library(car)
residualPlots(lm(y ~ Block + A, data.rcb1))

           Test stat Pr(>|t|)
Block             NA       NA
A                 NA       NA
Tukey test    -0.885    0.376

# the Tukey's non-additivity test by itself can be obtained via an
# internal function within the car package
car:::tukeyNonaddTest(lm(y ~ Block + A, data.rcb1))

      Test     Pvalue 
-0.8854414  0.3759186 

# alternatively, there is also a Tukey's non-additivity test
# within the asbio package
library(asbio)
with(data.rcb1, tukey.add.test(y, A, Block))

Tukey's one df test for additivity 
F = 0.7840065   Denom df = 67    p-value = 0.3790855

Conclusions:

• there is no visual or inferential evidence of any major interactions between Block and the within-Block effect (A). Any trends appear to be reasonably consistent between Blocks.

Model fitting or statistical analysis

There are numerous ways of fitting a nested ANOVA in R.

library(nlme)
# random intercept
data.rcb1.lme <- lme(y ~ A, random = ~1 | Block, data.rcb1, method = "REML")
# random intercept/slope
data.rcb1.lme1 <- lme(y ~ A, random = ~A | Block, data.rcb1, method = "REML")
anova(data.rcb1.lme, data.rcb1.lme1)

               Model df      AIC      BIC    logLik   Test  L.Ratio p-value
data.rcb1.lme      1  5 722.1087 735.2336 -356.0544                        
data.rcb1.lme1     2 10 727.2001 753.4499 -353.6001 1 vs 2 4.908574  0.4271

library(lme4)
data.rcb1.lmer <- lmer(y ~ A + (1 | Block), data.rcb1, REML = TRUE)  #random intercept
data.rcb1.lmer1 <- lmer(y ~ A + (A | Block), data.rcb1, REML = TRUE,
    control = lmerControl(check.nobs.vs.nRE = "ignore"))  #random intercept/slope
anova(data.rcb1.lmer, data.rcb1.lmer1)

Data: data.rcb1
Models:
data.rcb1.lmer: y ~ A + (1 | Block)
data.rcb1.lmer1: y ~ A + (A | Block)
                Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
data.rcb1.lmer   5 729.03 742.30 -359.51   719.03                         
data.rcb1.lmer1 10 733.98 760.52 -356.99   713.98 5.0529      5     0.4095

library(glmmTMB)
data.rcb1.glmmTMB <- glmmTMB(y ~ A + (1 | Block), data.rcb1)  #random intercept
data.rcb1.glmmTMB1 <- glmmTMB(y ~ A + (A | Block), data.rcb1)  #random intercept/slope
anova(data.rcb1.glmmTMB, data.rcb1.glmmTMB1)

Data: data.rcb1
Models:
data.rcb1.glmmTMB: y ~ A + (1 | Block), zi=~0, disp=~1
data.rcb1.glmmTMB1: y ~ A + (A | Block), zi=~0, disp=~1
                   Df    AIC   BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
data.rcb1.glmmTMB   5 729.03 742.3 -359.51   719.03                        
data.rcb1.glmmTMB1 10                                          5           

data.rcb1.aov <- aov(y ~ A + Error(Block), data.rcb1)

Model evaluation

Residuals

As always, exploring the residuals can reveal issues of heteroscadacity, non-linearity and potential issues with autocorrelation. Note for lme() and lmer() residual plots use standardized (normalized) residuals rather than raw residuals as the former reflect changes to the variance-covariance matrix whereas the later do not.

The following function will be used for the production of some of the qqnormal plots.

qq.line = function(x) {
    # following four lines from base R's qqline()
    y <- quantile(x[!is.na(x)], c(0.25, 0.75))
    x <- qnorm(c(0.25, 0.75))
    slope <- diff(y)/diff(x)
    int <- y[1L] - slope * x[1L]
    return(c(int = int, slope = slope))
}

• lme
• lmer
• glmmTMB
• Traditional aov

plot(data.rcb1.lme)

qqnorm(resid(data.rcb1.lme))
qqline(resid(data.rcb1.lme))

library(sjPlot)
plot_grid(plot_model(data.rcb1.lme, type = "diag"))

plot(data.rcb1.lmer)

plot(fitted(data.rcb1.lmer), residuals(data.rcb1.lmer, type = "pearson",
    scaled = TRUE))

ggplot(fortify(data.rcb1.lmer), aes(y = .scresid, x = .fitted)) +
    geom_point()

QQline = qq.line(fortify(data.rcb1.lmer)$.scresid)
ggplot(fortify(data.rcb1.lmer), aes(sample = .scresid)) + stat_qq() +
    geom_abline(intercept = QQline[1], slope = QQline[2])

qqnorm(resid(data.rcb1.lmer))
qqline(resid(data.rcb1.lmer))

library(sjPlot)
plot_grid(plot_model(data.rcb1.lmer, type = "diag"))

ggplot(data = NULL, aes(y = resid(data.rcb1.glmmTMB, type = "pearson"),
    x = fitted(data.rcb1.glmmTMB))) + geom_point()

QQline = qq.line(resid(data.rcb1.glmmTMB, type = "pearson"))
ggplot(data = NULL, aes(sample = resid(data.rcb1.glmmTMB, type = "pearson"))) +
    stat_qq() + geom_abline(intercept = QQline[1], slope = QQline[2])

library(sjPlot)
plot_grid(plot_model(data.rcb1.glmmTMB, type = "diag"))  #not working yet - bug

Error in UseMethod("rstudent"): no applicable method for 'rstudent' applied to an object of class "glmmTMB"

par(mfrow = c(2, 2))
plot(lm(data.rcb1.aov))

Exploring model parameters

If there was any evidence that the assumptions had been violated, then we would need to reconsider the model and start the process again. In this case, there is no evidence that the test will be unreliable so we can proceed to explore the test statistics. As I had elected to illustrate multiple techniques for analysing this nested design, I will also deal with the summaries etc separately.

Partial effects plots

It is often useful to visualize partial effects plots while exploring the parameter estimates. Having a graphical representation of the partial effects typically makes it a lot easier to interpret the parameter estimates and inferences.

• lme
• lmer
• glmmTMB

library(effects)
plot(allEffects(data.rcb1.lme))

library(sjPlot)
plot_model(data.rcb1.lme, type = "eff", terms = "A")

library(effects)
plot(allEffects(data.rcb1.lmer))

library(sjPlot)
plot_model(data.rcb1.lmer, type = "eff", terms = "A")

These are not yet correct!

library(ggeffects)
# observation level effects averaged across margins
p = ggaverage(data.rcb1.glmmTMB, terms = "A")
p = cbind(p, A = levels(data.rcb1$A))
ggplot(p, aes(y = predicted, x = A)) + geom_pointrange(aes(ymin = conf.low,
    ymax = conf.high))

# marginal effects
p = ggpredict(data.rcb1.glmmTMB, terms = "A")
p = cbind(p, A = levels(data.rcb1$A))
ggplot(p, aes(y = predicted, x = A)) + geom_pointrange(aes(ymin = conf.low,
    ymax = conf.high))

Extractor functions

Parameter estimates

summary(data.rcb1.lme)

Linear mixed-effects model fit by REML
 Data: data.rcb1 
       AIC      BIC    logLik
  722.1087 735.2336 -356.0544

Random effects:
 Formula: ~1 | Block
        (Intercept) Residual
StdDev:    11.51409 4.572284

Fixed effects: y ~ A 
               Value Std.Error DF  t-value p-value
(Intercept) 43.03434  2.094074 68 20.55053       0
AB          28.45241  1.092985 68 26.03185       0
AC          40.15556  1.092985 68 36.73936       0
 Correlation: 
   (Intr) AB    
AB -0.261       
AC -0.261  0.500

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.78748258 -0.57867597 -0.07108159  0.49990644  2.33727672 

Number of Observations: 105
Number of Groups: 35 

intervals(data.rcb1.lme)

Approximate 95% confidence intervals

 Fixed effects:
               lower     est.    upper
(Intercept) 38.85568 43.03434 47.21300
AB          26.27140 28.45241 30.63343
AC          37.97455 40.15556 42.33658
attr(,"label")
[1] "Fixed effects:"

 Random Effects:
  Level: Block 
                   lower     est.    upper
sd((Intercept)) 8.964242 11.51409 14.78924

 Within-group standard error:
   lower     est.    upper 
3.864949 4.572284 5.409070 

anova(data.rcb1.lme)

            numDF denDF   F-value p-value
(Intercept)     1    68 1089.3799  <.0001
A               2    68  714.0295  <.0001

library(broom)
tidy(data.rcb1.lme, effects = "fixed", conf.int = TRUE)

# A tibble: 3 x 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)     43.0      2.09      20.6 6.99e-31
2 AB              28.5      1.09      26.0 4.40e-37
3 AC              40.2      1.09      36.7 1.40e-46

glance(data.rcb1.lme)

# A tibble: 1 x 5
  sigma logLik   AIC   BIC deviance
  <dbl>  <dbl> <dbl> <dbl> <lgl>   
1  4.57  -356.  722.  735. NA      

summary(data.rcb1.lmer)

Linear mixed model fit by REML ['lmerMod']
Formula: y ~ A + (1 | Block)
   Data: data.rcb1

REML criterion at convergence: 712.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.78748 -0.57868 -0.07108  0.49991  2.33728 

Random effects:
 Groups   Name        Variance Std.Dev.
 Block    (Intercept) 132.57   11.514  
 Residual              20.91    4.572  
Number of obs: 105, groups:  Block, 35

Fixed effects:
            Estimate Std. Error t value
(Intercept)   43.034      2.094   20.55
AB            28.452      1.093   26.03
AC            40.156      1.093   36.74

Correlation of Fixed Effects:
   (Intr) AB    
AB -0.261       
AC -0.261  0.500

confint(data.rcb1.lmer)

                2.5 %    97.5 %
.sig01       8.996205 14.784207
.sigma       3.851781  5.370067
(Intercept) 38.894075 47.174609
AB          26.311725 30.593100
AC          38.014876 42.296251

anova(data.rcb1.lmer)

Analysis of Variance Table
  Df Sum Sq Mean Sq F value
A  2  29855   14927  714.03

library(broom)
tidy(data.rcb1.lmer, effects = "fixed", conf.int = TRUE)

# A tibble: 3 x 6
  term        estimate std.error statistic conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
1 (Intercept)     43.0      2.09      20.6     38.9      47.1
2 AB              28.5      1.09      26.0     26.3      30.6
3 AC              40.2      1.09      36.7     38.0      42.3

glance(data.rcb1.lmer)

# A tibble: 1 x 6
  sigma logLik   AIC   BIC deviance df.residual
  <dbl>  <dbl> <dbl> <dbl>    <dbl>       <int>
1  4.57  -356.  722.  735.     719.         100

summary(data.rcb1.glmmTMB)

 Family: gaussian  ( identity )
Formula:          y ~ A + (1 | Block)
Data: data.rcb1

     AIC      BIC   logLik deviance df.resid 
   729.0    742.3   -359.5    719.0      100 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 Block    (Intercept) 128.79   11.348  
 Residual              20.31    4.506  
Number of obs: 105, groups:  Block, 35

Dispersion estimate for gaussian family (sigma^2): 20.3 

Conditional model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   43.034      2.064   20.85   <2e-16 ***
AB            28.452      1.077   26.41   <2e-16 ***
AC            40.156      1.077   37.28   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

confint(data.rcb1.glmmTMB)

                                   2.5 %    97.5 %  Estimate
cond.(Intercept)               38.989093 47.079596 43.034345
cond.AB                        26.341021 30.563794 28.452407
cond.AC                        38.044174 42.266947 40.155561
cond.Std.Dev.Block.(Intercept)  8.867121 14.524046 11.348413
sigma                           3.818554  5.318367  4.506492

summary(data.rcb1.aov)

Error: Block
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 34  14233   418.6               

Error: Within
          Df Sum Sq Mean Sq F value Pr(>F)    
A          2  29855   14927     714 <2e-16 ***
Residuals 68   1422      21                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Planned comparisons and pairwise post-hoc tests

As with non-heirarchical models, we can incorporate alternative contrasts for the fixed effects (other than the default treatment contrasts). The random factors must be sum-to-zero contrasts in order to ensure that the model is identifiable (possible to estimate true values of the parameters).

Likewise, post-hoc tests such as Tukey's tests can be performed.

library(multcomp)
summary(glht(data.rcb1.lme, linfct = mcp(A = "Tukey")))

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lme.formula(fixed = y ~ A, data = data.rcb1, random = ~1 | Block, 
    method = "REML")

Linear Hypotheses:
           Estimate Std. Error z value Pr(>|z|)    
B - A == 0   28.452      1.093   26.03   <2e-16 ***
C - A == 0   40.156      1.093   36.74   <2e-16 ***
C - B == 0   11.703      1.093   10.71   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

library(broom)
tidy(confint(glht(data.rcb1.lme, linfct = mcp(A = "Tukey"))))

# A tibble: 3 x 5
  lhs     rhs estimate conf.low conf.high
  <chr> <dbl>    <dbl>    <dbl>     <dbl>
1 B - A    0.     28.5    25.9       31.0
2 C - A    0.     40.2    37.6       42.7
3 C - B    0.     11.7     9.14      14.3

library(lsmeans)
lsmeans(data.rcb1.lme, pairwise ~ A)

$lsmeans
 A   lsmean       SE df lower.CL upper.CL
 A 43.03434 2.094074 34 38.77867 47.29001
 B 71.48675 2.094074 34 67.23108 75.74242
 C 83.18991 2.094074 34 78.93423 87.44558

Confidence level used: 0.95 

$contrasts
 contrast  estimate       SE df t.ratio p.value
 A - B    -28.45241 1.092985 68 -26.032  <.0001
 A - C    -40.15556 1.092985 68 -36.739  <.0001
 B - C    -11.70315 1.092985 68 -10.708  <.0001

P value adjustment: tukey method for comparing a family of 3 estimates 

summary(glht(data.rcb1.lme, linfct = lsm(pairwise ~ A)))

	 Simultaneous Tests for General Linear Hypotheses

Fit: lme.formula(fixed = y ~ A, data = data.rcb1, random = ~1 | Block, 
    method = "REML")

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)    
A - B == 0  -28.452      1.093  -26.03   <1e-10 ***
A - C == 0  -40.156      1.093  -36.74   <1e-10 ***
B - C == 0  -11.703      1.093  -10.71   <1e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

confint(glht(data.rcb1.lme, linfct = lsm(pairwise ~ A)))

	 Simultaneous Confidence Intervals

Fit: lme.formula(fixed = y ~ A, data = data.rcb1, random = ~1 | Block, 
    method = "REML")

Quantile = 2.3963
95% family-wise confidence level
 

Linear Hypotheses:
           Estimate lwr      upr     
A - B == 0 -28.4524 -31.0715 -25.8333
A - C == 0 -40.1556 -42.7746 -37.5365
B - C == 0 -11.7032 -14.3222  -9.0841

library(multcomp)
contr = rbind(`B-C` = c(0, 1, -1), `A-(B,C)` = c(1, -0.5, -0.5))
g = glht(data.rcb1.lme, linfct = mcp(A = contr))
summary(g, test = adjusted("none"))

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts


Fit: lme.formula(fixed = y ~ A, data = data.rcb1, random = ~1 | Block, 
    method = "REML")

Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)    
B-C == 0     -11.7032     1.0930  -10.71   <2e-16 ***
A-(B,C) == 0 -34.3040     0.9466  -36.24   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- none method)

confint(g)

	 Simultaneous Confidence Intervals

Multiple Comparisons of Means: User-defined Contrasts


Fit: lme.formula(fixed = y ~ A, data = data.rcb1, random = ~1 | Block, 
    method = "REML")

Quantile = 2.2364
95% family-wise confidence level
 

Linear Hypotheses:
             Estimate lwr      upr     
B-C == 0     -11.7032 -14.1475  -9.2588
A-(B,C) == 0 -34.3040 -36.4209 -32.1871

library(broom)
tidy(confint(g))

# A tibble: 2 x 5
  lhs       rhs estimate conf.low conf.high
  <chr>   <dbl>    <dbl>    <dbl>     <dbl>
1 B-C        0.    -11.7    -14.1     -9.26
2 A-(B,C)    0.    -34.3    -36.4    -32.2 

# OR manually
cmat = cbind(`B-C` = c(0, 1, -1), `A-(B,C)` = c(1, -0.5, -0.5))
crossprod(cmat)

        B-C A-(B,C)
B-C       2     0.0
A-(B,C)   0     1.5

newdata = data.frame(A = levels(data.rcb1$A))
Xmat = model.matrix(~A, data = newdata)
Xmat = t(cmat) %*% Xmat
coefs = fixef(data.rcb1.lme)
fit = as.vector(coefs %*% t(Xmat))
se = sqrt(diag(Xmat %*% vcov(data.rcb1.lme) %*% t(Xmat)))
q = qt(0.975, df = nrow(data.rcb1.lme$data) - length(coefs) - 2)
newdata = data.frame(Contrast = rownames(Xmat), fit = fit, lower = fit - q * se, upper = fit + q * se)
newdata

        Contrast       fit     lower     upper
B-C          B-C -11.70315 -13.87160  -9.53470
A-(B,C)  A-(B,C) -34.30399 -36.18192 -32.42605

library(multcomp)
summary(glht(data.rcb1.lmer, linfct = mcp(A = "Tukey")))

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lme4::lmer(formula = y ~ A + (1 | Block), data = data.rcb1, REML = TRUE)

Linear Hypotheses:
           Estimate Std. Error z value Pr(>|z|)    
B - A == 0   28.452      1.093   26.03   <2e-16 ***
C - A == 0   40.156      1.093   36.74   <2e-16 ***
C - B == 0   11.703      1.093   10.71   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

library(broom)
tidy(confint(glht(data.rcb1.lmer, linfct = mcp(A = "Tukey"))))

# A tibble: 3 x 5
  lhs     rhs estimate conf.low conf.high
  <chr> <dbl>    <dbl>    <dbl>     <dbl>
1 B - A    0.     28.5    25.9       31.0
2 C - A    0.     40.2    37.6       42.7
3 C - B    0.     11.7     9.14      14.3

library(lsmeans)
lsmeans(data.rcb1.lmer, pairwise ~ A)

$lsmeans
 A   lsmean       SE    df lower.CL upper.CL
 A 43.03434 2.094074 40.93 38.80504 47.26364
 B 71.48675 2.094074 40.93 67.25746 75.71605
 C 83.18991 2.094074 40.93 78.96061 87.41920

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 

$contrasts
 contrast  estimate       SE df t.ratio p.value
 A - B    -28.45241 1.092985 68 -26.032  <.0001
 A - C    -40.15556 1.092985 68 -36.739  <.0001
 B - C    -11.70315 1.092985 68 -10.708  <.0001

P value adjustment: tukey method for comparing a family of 3 estimates 

summary(glht(data.rcb1.lmer, linfct = lsm(pairwise ~ A)))

	 Simultaneous Tests for General Linear Hypotheses

Fit: lme4::lmer(formula = y ~ A + (1 | Block), data = data.rcb1, REML = TRUE)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)    
A - B == 0  -28.452      1.093  -26.03   <2e-16 ***
A - C == 0  -40.156      1.093  -36.74   <2e-16 ***
B - C == 0  -11.703      1.093  -10.71   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

confint(glht(data.rcb1.lmer, linfct = lsm(pairwise ~ A)))

	 Simultaneous Confidence Intervals

Fit: lme4::lmer(formula = y ~ A + (1 | Block), data = data.rcb1, REML = TRUE)

Quantile = 2.3959
95% family-wise confidence level
 

Linear Hypotheses:
           Estimate lwr      upr     
A - B == 0 -28.4524 -31.0711 -25.8337
A - C == 0 -40.1556 -42.7743 -37.5369
B - C == 0 -11.7032 -14.3218  -9.0845

library(multcomp)
contr = rbind(`B-C` = c(0, 1, -1), `A-(B,C)` = c(1, -0.5, -0.5))
g = glht(data.rcb1.lmer, linfct = mcp(A = contr))
summary(g, test = adjusted("none"))

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts


Fit: lme4::lmer(formula = y ~ A + (1 | Block), data = data.rcb1, REML = TRUE)

Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)    
B-C == 0     -11.7032     1.0930  -10.71   <2e-16 ***
A-(B,C) == 0 -34.3040     0.9466  -36.24   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- none method)

confint(g)

	 Simultaneous Confidence Intervals

Multiple Comparisons of Means: User-defined Contrasts


Fit: lme4::lmer(formula = y ~ A + (1 | Block), data = data.rcb1, REML = TRUE)

Quantile = 2.2364
95% family-wise confidence level
 

Linear Hypotheses:
             Estimate lwr      upr     
B-C == 0     -11.7032 -14.1475  -9.2588
A-(B,C) == 0 -34.3040 -36.4209 -32.1871

# OR manually
cmat = cbind(`B-C` = c(0, 1, -1), `A-(B,C)` = c(1, -0.5, -0.5))
crossprod(cmat)

        B-C A-(B,C)
B-C       2     0.0
A-(B,C)   0     1.5

newdata = data.frame(A = levels(data.rcb1$A))
Xmat = model.matrix(~A, data = newdata)
Xmat = t(cmat) %*% Xmat
coefs = fixef(data.rcb1.lmer)
fit = as.vector(coefs %*% t(Xmat))
se = sqrt(diag(Xmat %*% vcov(data.rcb1.lmer) %*% t(Xmat)))
q = qt(0.975, df = df.residual(data.rcb1.lmer))
newdata = data.frame(Contrast = rownames(Xmat), fit = fit, lower = fit - q * se, upper = fit + q * se)
newdata

  Contrast       fit     lower     upper
1      B-C -11.70315 -13.87160  -9.53470
2  A-(B,C) -34.30399 -36.18192 -32.42605

library(multcomp)
tuk.mat <- contrMat(n = table(data.rcb1$A), type = "Tukey")
newdata = data.frame(A = levels(data.rcb1$A))
Xmat = model.matrix(~A, data = newdata)
Xmat = tuk.mat %*% Xmat
coefs = fixef(data.rcb1.glmmTMB)[[1]]
fit = as.vector(coefs %*% t(Xmat))
se = sqrt(diag(Xmat %*% vcov(data.rcb1.glmmTMB)[[1]] %*% t(Xmat)))
q = qt(0.975, df = df.residual(data.rcb1.glmmTMB))
newdata = data.frame(Contrast = rownames(Xmat), fit = fit, lower = fit - q * se, upper = fit + q * se)
newdata

      Contrast      fit     lower    upper
B - A    B - A 28.45241 26.315159 30.58966
C - A    C - A 40.15556 38.018312 42.29281
C - B    C - B 11.70315  9.565905 13.84040

library(lsmeans)
recover.data.glmmTMB <- function(object, ...) {
    fcall <- getCall(object)
    recover.data(fcall, delete.response(terms(object)), attr(model.frame(object),
        "na.action"), ...)
}
lsm.basis.glmmTMB <- function(object, trms, xlev, grid, vcov., mode = "asymptotic",
    component = "cond", ...) {
    if (mode != "asymptotic")
        stop("only asymptotic mode is available")
    if (component != "cond")
        stop("only tested for conditional component")
    if (missing(vcov.))
        V <- as.matrix(vcov(object)[[component]]) else V <- as.matrix(.my.vcov(object, vcov.))
    dfargs = misc = list()
    if (mode == "asymptotic") {
        dffun = function(k, dfargs) NA
    }
    ## use this? misc = .std.link.labels(family(object), misc)
    contrasts = attr(model.matrix(object), "contrasts")
    m = model.frame(trms, grid, na.action = na.pass, xlev = xlev)
    X = model.matrix(trms, m, contrasts.arg = contrasts)
    bhat = fixef(object)[[component]]
    if (length(bhat) < ncol(X)) {
        kept = match(names(bhat), dimnames(X)[[2]])
        bhat = NA * X[1, ]
        bhat[kept] = fixef(object)[[component]]
        modmat = model.matrix(trms, model.frame(object), contrasts.arg = contrasts)
        nbasis = estimability::nonest.basis(modmat)
    } else nbasis = estimability::all.estble
    list(X = X, bhat = bhat, nbasis = nbasis, V = V, dffun = dffun,
        dfargs = dfargs, misc = misc)
}
lsmeans(data.rcb1.glmmTMB, pairwise ~ A)

$lsmeans
 A   lsmean       SE df asymp.LCL asymp.UCL
 A 43.03434 2.063942 NA  38.98909  47.07960
 B 71.48675 2.063942 NA  67.44150  75.53200
 C 83.18991 2.063942 NA  79.14465  87.23516

Confidence level used: 0.95 

$contrasts
 contrast  estimate       SE df z.ratio p.value
 A - B    -28.45241 1.077258 NA -26.412  <.0001
 A - C    -40.15556 1.077258 NA -37.276  <.0001
 B - C    -11.70315 1.077258 NA -10.864  <.0001

P value adjustment: tukey method for comparing a family of 3 estimates 

# summary(glht(data.rcb1.lmer, linfct=lsm(pairwise ~ A)))
# confint(glht(data.rcb1.lmer, linfct=lsm(pairwise ~ A)))

cmat = cbind(`B-C` = c(0, 1, -1), `A-(B,C)` = c(1, -0.5, -0.5))
crossprod(cmat)

        B-C A-(B,C)
B-C       2     0.0
A-(B,C)   0     1.5

newdata = data.frame(A = levels(data.rcb1$A))
Xmat = model.matrix(~A, data = newdata)
Xmat = t(cmat) %*% Xmat
coefs = fixef(data.rcb1.glmmTMB)[[1]]
fit = as.vector(coefs %*% t(Xmat))
se = sqrt(diag(Xmat %*% vcov(data.rcb1.glmmTMB)[[1]] %*% t(Xmat)))
q = qt(0.975, df = df.residual(data.rcb1.glmmTMB))
newdata = data.frame(Contrast = rownames(Xmat), fit = fit, lower = fit - q * se, upper = fit + q * se)
newdata

        Contrast       fit    lower      upper
B-C          B-C -11.70315 -13.8404  -9.565905
A-(B,C)  A-(B,C) -34.30398 -36.1549 -32.453072

Predictions

As with other linear models, it is possible to generate predicted values from the fitted model. Since the linear mixed effects model (with random intercepts) captures information on the levels of the random effects, we can indicate multiple hierarchy from which predictions could be generated. For example, do we wish to predict a new value of $Y$ from a specific level of $A$ regardless of the level of the random effect(s) - this is like predicting a new value at a random level of $A$ and is the typical case. Alternatively we could be interested in predicting the value of $Y$ at a specific level of $A$ and for a specific level of the random factor (in this case $Block$).

Note, the predict() function does not provide confidence or prediction intervals for mixed effects models. It they are wanted then they need to be calculated manually.

predict(data.rcb1.lme, newdata = data.frame(A = levels(data.rcb1$A)), level = 0)

[1] 43.03434 71.48675 83.18991
attr(,"label")
[1] "Predicted values"

library(ggeffects)
ggpredict(data.rcb1.lme, terms = "A", x.as.factor = TRUE)

# A tibble: 3 x 5
  x     predicted conf.low conf.high group
  <fct>     <dbl>    <dbl>     <dbl> <fct>
1 A          43.0     38.9      47.1 1    
2 B          71.5     67.4      75.6 1    
3 C          83.2     79.1      87.3 1    

library(effects)
as.data.frame(Effect(focal = "A", mod = data.rcb1.lme))

  A      fit       se    lower    upper
1 A 43.03434 2.094074 38.88076 47.18793
2 B 71.48675 2.094074 67.33317 75.64034
3 C 83.18991 2.094074 79.03632 87.34349

library(lsmeans)
lsmeans(data.rcb1.lme, eff ~ A)$lsmeans

 A   lsmean       SE df lower.CL upper.CL
 A 43.03434 2.094074 34 38.77867 47.29001
 B 71.48675 2.094074 34 67.23108 75.74242
 C 83.18991 2.094074 34 78.93423 87.44558

Confidence level used: 0.95 

newdata = expand.grid(A = levels(data.rcb1$A), Block = levels(data.rcb1$Block)[c(3, 5)])
predict(data.rcb1.lme, newdata = newdata, level = 1)

      B3       B3       B3       B5       B5       B5 
39.33369 67.78610 79.48925 49.68656 78.13898 89.84213 
attr(,"label")
[1] "Predicted values"

# OR
newdata1 <- expand.grid(A = levels(data.rcb1$A), Block = levels(data.rcb1$Block)[c(3, 5)])
augment(data.rcb1.lme, newdata = newdata1)

# A tibble: 6 x 3
  A     Block .fitted
  <fct> <fct>   <dbl>
1 A     B3       39.3
2 B     B3       67.8
3 C     B3       79.5
4 A     B5       49.7
5 B     B5       78.1
6 C     B5       89.8

# Manual confidence intervals
newdata1 <- expand.grid(A = levels(data.rcb1$A), Block = levels(data.rcb1$Block)[c(3, 5)])
coefs <- as.matrix(coef(data.rcb1.lme))[c(3, 5), ]
Xmat <- model.matrix(~A, newdata1)
## Split the Xmat and Coefs up by Blocks
Xmat.list = split(as.data.frame(Xmat), newdata1$Block)
coefs.list = split(coefs, rownames(coefs))
## Perform matrix multiplication listwise
fit = unlist(Map("%*%", coefs.list, lapply(Xmat.list, function(x) t(as.matrix(x)))))
se <- sqrt(diag(Xmat %*% vcov(data.rcb1.lme) %*% t(Xmat)))
q = qt(0.975, df = nrow(data.rcb1.lme$data) - length(coefs) - 1)
(newdata1 = cbind(newdata1, fit = fit, lower = fit - q * se, upper = fit + q * se))

    A Block      fit    lower    upper
B31 A    B3 39.33369 35.17807 43.48931
B32 B    B3 67.78610 63.63048 71.94173
B33 C    B3 79.48925 75.33363 83.64488
B51 A    B5 49.68656 45.53094 53.84218
B52 B    B5 78.13898 73.98335 82.29460
B53 C    B5 89.84213 85.68650 93.99775

## Manual prediction invervals
sigma = sigma(data.rcb1.lme)
(newdata1 = cbind(newdata1, lowerP = fit - q * (se * sigma), upperP = fit + q * (se * sigma)))

    A Block      fit    lower    upper   lowerP    upperP
B31 A    B3 39.33369 35.17807 43.48931 20.33301  58.33437
B32 B    B3 67.78610 63.63048 71.94173 48.78542  86.78679
B33 C    B3 79.48925 75.33363 83.64488 60.48857  98.48994
B51 A    B5 49.68656 45.53094 53.84218 30.68588  68.68724
B52 B    B5 78.13898 73.98335 82.29460 59.13829  97.13966
B53 C    B5 89.84213 85.68650 93.99775 70.84144 108.84281

predict(data.rcb1.lmer, newdata = data.frame(A = levels(data.rcb1$A)), re.form = NA)

       1        2        3 
43.03434 71.48675 83.18991 

library(ggeffects)
ggpredict(data.rcb1.lmer, terms = "A", x.as.factor = TRUE)

# A tibble: 3 x 5
  x     predicted conf.low conf.high group
  <fct>     <dbl>    <dbl>     <dbl> <fct>
1 A          43.0     38.9      47.1 1    
2 B          71.5     67.4      75.6 1    
3 C          83.2     79.1      87.3 1    

library(effects)
as.data.frame(Effect(focal = "A", mod = data.rcb1.lmer))

  A      fit       se    lower    upper
1 A 43.03434 2.094074 38.88076 47.18793
2 B 71.48675 2.094074 67.33317 75.64034
3 C 83.18991 2.094074 79.03632 87.34349

library(lsmeans)
lsmeans(data.rcb1.lmer, eff ~ A)$lsmeans

 A   lsmean       SE    df lower.CL upper.CL
 A 43.03434 2.094074 40.93 38.80504 47.26364
 B 71.48675 2.094074 40.93 67.25746 75.71605
 C 83.18991 2.094074 40.93 78.96061 87.41920

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 

newdata = expand.grid(A = levels(data.rcb1$A), Block = levels(data.rcb1$Block)[c(3, 5)])
predict(data.rcb1.lmer, newdata = newdata, re.form = ~1 | Block)

       1        2        3        4        5        6 
39.33369 67.78610 79.48925 49.68656 78.13898 89.84213 

# OR
newdata1 <- expand.grid(A = levels(data.rcb1$A), Block = levels(data.rcb1$Block)[c(3, 5)])
augment(data.rcb1.lmer, newdata = newdata1)

  A Block  .fitted      .mu .offset .sqrtXwt .sqrtrwt .weights     .wtres
1 A    B3 39.33369 36.45695       0        1        1        1  0.9406589
2 B    B3 67.78610 64.90937       0        1        1        1 -3.4390334
3 C    B3 79.48925 76.61252       0        1        1        1  1.4611797
4 A    B5 49.68656 33.09453       0        1        1        1 -2.4964946
5 B    B5 78.13898 61.54694       0        1        1        1 -2.5465860
6 C    B5 89.84213 73.25009       0        1        1        1  3.4756614

# Manual confidence intervals
newdata1 <- expand.grid(A = levels(data.rcb1$A), Block = levels(data.rcb1$Block)[c(3, 5)])
coefs <- as.matrix(coef(data.rcb1.lmer)$Block)[c(3, 5), ]
Xmat <- model.matrix(~A, newdata1)
## Split the Xmat and Coefs up by Blocks
Xmat.list = split(as.data.frame(Xmat), newdata1$Block)
coefs.list = split(coefs, rownames(coefs))
## Perform matrix multiplication listwise
fit = unlist(Map("%*%", coefs.list, lapply(Xmat.list, function(x) t(as.matrix(x)))))
se <- sqrt(diag(Xmat %*% vcov(data.rcb1.lmer) %*% t(Xmat)))
q = qt(0.975, df = df.residual(data.rcb1.lmer))
(newdata1 = cbind(newdata1, fit = fit, lower = fit - q * se, upper = fit + q * se))

    A Block      fit    lower    upper
B31 A    B3 39.33369 35.17911 43.48827
B32 B    B3 67.78610 63.63152 71.94069
B33 C    B3 79.48925 75.33467 83.64384
B51 A    B5 49.68656 45.53198 53.84115
B52 B    B5 78.13898 73.98439 82.29356
B53 C    B5 89.84213 85.68754 93.99671

## Manual prediction invervals
sigma = sigma(data.rcb1.lmer)
(newdata1 = cbind(newdata1, lowerP = fit - q * (se * sigma), upperP = fit + q * (se * sigma)))

    A Block      fit    lower    upper   lowerP    upperP
B31 A    B3 39.33369 35.17911 43.48827 20.33776  58.32962
B32 B    B3 67.78610 63.63152 71.94069 48.79017  86.78204
B33 C    B3 79.48925 75.33467 83.64384 60.49332  98.48519
B51 A    B5 49.68656 45.53198 53.84115 30.69063  68.68250
B52 B    B5 78.13898 73.98439 82.29356 59.14304  97.13491
B53 C    B5 89.84213 85.68754 93.99671 70.84619 108.83806

predict(data.rcb1.glmmTMB, newdata = data.frame(A = levels(data.rcb1$A)), re.form = NA)

Error in predict.glmmTMB(data.rcb1.glmmTMB, newdata = data.frame(A = levels(data.rcb1$A)), : re.form not yet implemented

library(ggeffects)
ggpredict(data.rcb1.glmmTMB, terms = "A", x.as.factor = TRUE)

# A tibble: 3 x 5
  x     predicted conf.low conf.high group
  <fct>     <dbl>    <dbl>     <dbl> <fct>
1 A          36.5     31.3      41.6 1    
2 B          64.9     59.8      70.0 1    
3 C          76.6     71.5      81.7 1    

library(effects)
as.data.frame(Effect(focal = "A", mod = data.rcb1.glmmTMB))

  A      fit       se    lower    upper
1 A 43.03434 2.094074 38.88076 47.18793
2 B 71.48675 2.094074 67.33317 75.64034
3 C 83.18991 2.094074 79.03632 87.34349

library(lsmeans)
lsmeans(data.rcb1.glmmTMB, eff ~ A)$lsmeans

 A   lsmean       SE df asymp.LCL asymp.UCL
 A 43.03434 2.063942 NA  38.98909  47.07960
 B 71.48675 2.063942 NA  67.44150  75.53200
 C 83.18991 2.063942 NA  79.14465  87.23516

Confidence level used: 0.95 

newdata = expand.grid(A = levels(data.rcb1$A), Block = levels(data.rcb1$Block)[c(3, 5)])
predict(data.rcb1.glmmTMB, newdata = newdata, re.form = ~1 | Block)

Error in predict.glmmTMB(data.rcb1.glmmTMB, newdata = newdata, re.form = ~1 | : re.form not yet implemented

# OR
newdata1 <- expand.grid(A = levels(data.rcb1$A), Block = levels(data.rcb1$Block)[c(3, 5)])
augment(data.rcb1.glmmTMB, newdata = newdata1)

Error: No augment method for objects of class glmmTMB

# Manual confidence intervals
newdata1 <- expand.grid(A = levels(data.rcb1$A), Block = levels(data.rcb1$Block)[c(3, 5)])
coefs <- as.matrix(coef(data.rcb1.glmmTMB)$Block)[c(3, 5), ]

Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x), : 'data' must be of a vector type, was 'NULL'

Xmat <- model.matrix(~A, newdata1)
## Split the Xmat and Coefs up by Blocks
Xmat.list = split(as.data.frame(Xmat), newdata1$Block)
coefs.list = split(coefs, rownames(coefs))
## Perform matrix multiplication listwise
fit = unlist(Map("%*%", coefs.list, lapply(Xmat.list, function(x) t(as.matrix(x)))))
se <- sqrt(diag(Xmat %*% vcov(data.rcb1.glmmTMB) %*% t(Xmat)))

Error in Xmat %*% vcov(data.rcb1.glmmTMB): requires numeric/complex matrix/vector arguments

q = qt(0.975, df = df.residual(data.rcb1.glmmTMB))
(newdata1 = cbind(newdata1, fit = fit, lower = fit - q * se, upper = fit + q * se))

    A Block      fit    lower    upper
B31 A    B3 39.33369 35.17911 43.48827
B32 B    B3 67.78610 63.63152 71.94069
B33 C    B3 79.48925 75.33467 83.64384
B51 A    B5 49.68656 45.53198 53.84115
B52 B    B5 78.13898 73.98439 82.29356
B53 C    B5 89.84213 85.68754 93.99671

## Manual prediction invervals
sigma = sigma(data.rcb1.glmmTMB)
(newdata1 = cbind(newdata1, lowerP = fit - q * (se * sigma), upperP = fit + q * (se * sigma)))

    A Block      fit    lower    upper   lowerP    upperP
B31 A    B3 39.33369 35.17911 43.48827 20.61109  58.05629
B32 B    B3 67.78610 63.63152 71.94069 49.06351  86.50870
B33 C    B3 79.48925 75.33467 83.64384 60.76666  98.21185
B51 A    B5 49.68656 45.53198 53.84115 30.96396  68.40916
B52 B    B5 78.13898 73.98439 82.29356 59.41638  96.86157
B53 C    B5 89.84213 85.68754 93.99671 71.11953 108.56472

$R^2$ approximations

library(MuMIn)
r.squaredGLMM(data.rcb1.lme)

      R2m       R2c 
0.6516126 0.9525456 

library(MuMIn)
r.squaredGLMM(data.rcb1.lmer)

      R2m       R2c 
0.6516126 0.9525456 

source(system.file("misc/rsqglmm.R", package = "glmmTMB"))
my_rsq(data.rcb1.glmmTMB)

$family
[1] "gaussian"

$link
[1] "identity"

$Marginal
[1] 0.6581639

$Conditional
[1] 0.9534379

The fixed effect of A (within Block) accounts for approximately 65.16% of the total variation in Y. The random effect of Block accounts for approximately 30.09% of the total variation in Y and collectively, the hierarchical level of Block (containing the fixed effect) explains approximately 95.25% of the total variation in Y.

Graphical summary

It is relatively trivial to produce a summary figure based on the raw data. Arguably a more satisfying figure would be one based on the modelled data.

• lme
• lmer
• glmmTMB

The new effects package seems to have a bug when calculating effects

library(effects)
data.rcb1.eff = as.data.frame(allEffects(data.rcb1.lme)[[1]])
# OR
data.rcb1.eff = as.data.frame(Effect("A", data.rcb1.lme))

ggplot(data.rcb1.eff, aes(y = fit, x = A)) + geom_pointrange(aes(ymin = lower,
    ymax = upper)) + scale_y_continuous("Y") + theme_classic()

# OR using lsmeans
fit = summary(ref.grid(data.rcb1.lme), infer = TRUE)
ggplot(fit, aes(y = prediction, x = A)) + geom_pointrange(aes(ymin = lower.CL,
    ymax = upper.CL)) + scale_y_continuous("Y") + theme_classic()

# OR
fit = summary(lsmeans(data.rcb1.lme, eff ~ A)$lsmean)
ggplot(fit, aes(y = lsmean, x = A)) + geom_pointrange(aes(ymin = lower.CL,
    ymax = upper.CL)) + scale_y_continuous("Y") + theme_classic()

Or manually

newdata = data.frame(A = levels(data.rcb1$A))
Xmat = model.matrix(~A, data = newdata)
coefs = fixef(data.rcb1.lme)
fit = as.vector(coefs %*% t(Xmat))
se = sqrt(diag(Xmat %*% vcov(data.rcb1.lme) %*% t(Xmat)))
q = qt(0.975, df = nrow(data.rcb1.lme$data) - length(coefs) -
    2)
newdata = cbind(newdata, fit = fit, lower = fit - q * se, upper = fit +
    q * se)
ggplot(newdata, aes(y = fit, x = A)) + geom_pointrange(aes(ymin = lower,
    ymax = upper)) + scale_y_continuous("Y") + theme_classic()

The new effects package seems to have a bug when calculating effects

library(effects)
data.rcb1.eff = as.data.frame(allEffects(data.rcb1.lmer)[[1]])
# OR
data.rcb1.eff = as.data.frame(Effect("A", data.rcb1.lmer))

ggplot(data.rcb1.eff, aes(y = fit, x = A)) + geom_pointrange(aes(ymin = lower,
    ymax = upper)) + scale_y_continuous("Y") + theme_classic()

# OR using lsmeans
fit = summary(ref.grid(data.rcb1.lmer), infer = TRUE)
ggplot(fit, aes(y = prediction, x = A)) + geom_pointrange(aes(ymin = lower.CL,
    ymax = upper.CL)) + scale_y_continuous("Y") + theme_classic()

# OR
fit = summary(lsmeans(data.rcb1.lmer, eff ~ A)$lsmean)
ggplot(fit, aes(y = lsmean, x = A)) + geom_pointrange(aes(ymin = lower.CL,
    ymax = upper.CL)) + scale_y_continuous("Y") + theme_classic()

Or manually

newdata = data.frame(A = levels(data.rcb1$A))
Xmat = model.matrix(~A, data = newdata)
coefs = fixef(data.rcb1.lmer)
fit = as.vector(coefs %*% t(Xmat))
se = sqrt(diag(Xmat %*% vcov(data.rcb1.lmer) %*% t(Xmat)))
q = qt(0.975, df = df.residual(data.rcb1.lmer))
newdata = cbind(newdata, fit = fit, lower = fit - q * se, upper = fit +
    q * se)
ggplot(newdata, aes(y = fit, x = A)) + geom_pointrange(aes(ymin = lower,
    ymax = upper)) + scale_y_continuous("Y") + theme_classic()

The new effects package seems to have a bug when calculating effects

fit = summary(ref.grid(data.rcb1.glmmTMB), infer = TRUE)
ggplot(fit, aes(y = prediction, x = A)) + geom_pointrange(aes(ymin = asymp.LCL,
    ymax = asymp.UCL)) + scale_y_continuous("Y") + theme_classic()

# OR
fit = summary(lsmeans(data.rcb1.glmmTMB, eff ~ A)$lsmean)
ggplot(fit, aes(y = lsmean, x = A)) + geom_pointrange(aes(ymin = asymp.LCL,
    ymax = asymp.UCL)) + scale_y_continuous("Y") + theme_classic()

Or manually

newdata = data.frame(A = levels(data.rcb1$A))
Xmat = model.matrix(~A, data = newdata)
coefs = fixef(data.rcb1.glmmTMB)[[1]]
fit = as.vector(coefs %*% t(Xmat))
se = sqrt(diag(Xmat %*% vcov(data.rcb1.glmmTMB)[[1]] %*% t(Xmat)))
q = qt(0.975, df = df.residual(data.rcb1.glmmTMB))
newdata = cbind(newdata, fit = fit, lower = fit - q * se, upper = fit +
    q * se)
ggplot(newdata, aes(y = fit, x = A)) + geom_pointrange(aes(ymin = lower,
    ymax = upper)) + scale_y_continuous("Y") + theme_classic()

RCB (repeated measures) ANOVA in R - continuous within and AR1 temporal autocorrelation

Scenario and Data

Imagine now that we has designed an experiment to investigate the effects of a continuous predictor ($x$, for example time) on a response ($y$). Again, the system that we intend to sample is spatially heterogeneous and thus will add a great deal of noise to the data that will make it difficult to detect a signal (impact of treatment).

Thus in an attempt to constrain this variability, we again decide to apply a design (RCB) in which each of the levels of X (such as time) treatments within each of 35 blocks dispersed randomly throughout the landscape. As this section is mainly about the generation of artificial data (and not specifically about what to do with the data), understanding the actual details are optional and can be safely skipped. Consequently, I have folded (toggled) this section away.

Random data incorporating the following properties

• the number of times = 10
• the number of blocks containing treatments = 34
• the number of treatments = 2
• mean slope (rate of change in response over time) = 30
• mean difference in slope between the first and second treatments = 20
• mean intercept (value of response at time 0 = 200
• the variability (standard deviation) between blocks of the same treatment = 30
• the variability (standard deviation) in slope = 50

library(tidyverse)
set.seed(1)
slope <- 30
intercept <- 200
nBlock <- 34
nTime <- 10
sd.res = 50
sd.block = 10
sd.slope = 0.1
cor_int_slope = 0.4
rho = 0.8
## define the random effects (intercepts and slopes)
rand.effects = MASS::mvrnorm(n = nBlock, mu = c(0, 0), Sigma = matrix(c(1, cor_int_slope, cor_int_slope,
    1), nrow = 2), empirical = TRUE) %>% as.data.frame %>% mutate(V1 = V1 * sd.block + intercept, V2 = V2 *
    sd.slope + slope)
apply(rand.effects, 2, mean)

 V1  V2 
200  30 

apply(rand.effects, 2, sd)

  V1   V2 
10.0  0.1 

cor(rand.effects[, 1], rand.effects[, 2])

[1] 0.4

## define the residuals (in this case with an AR1 dependency structure)
cmat <- (sd.res^2) * rho^abs(outer(0:(nTime - 1), 0:(nTime - 1), "-"))
eps <- c(t(MASS::mvrnorm(nBlock, mu = rep(0, nTime), Sigma = cmat, empirical = TRUE)))
## create the predictor data
data.rm = data.frame(Block = factor(paste0("B", rep(1:nBlock, each = nTime))), Time = rep(1:nTime, times = nBlock))
## calculate the linear predictor
Xmat = model.matrix(~-1 + Block + Block:Time, data = data.rm)
block.effects = rand.effects[, 1]
slope.effects = rand.effects[, 2]
lin.pred = c(block.effects, slope.effects)
data.rm$y = as.vector(lin.pred %*% t(Xmat)) + eps

ggplot(data.rm, aes(y = y, x = Time)) + geom_point() + geom_smooth(method = "lm") + facet_wrap(~Block)

## cmat <- (sd.res^2)*rho^abs(outer(0:(nTime-1),0:(nTime-1),'-')) errs <-
## MASS::mvrnorm(nBlock,mu=rep(0,nTime),Sigma=cmat, empirical=TRUE)

## library(tidyverse) set.seed(123) slope <- 30 intercept <- 200 treatment <- 20 nBlock <- 34 nTreat
## <- 2 nTime <- 10 sigma <- 9 sigma.block <- 3# 30 n <- nBlock*nTime Block <- gl(nBlock, k=1) Time <-
## 1:10 rho <- 0.8

## dt <- expand.grid(Time=Time,Block=Block) %>% mutate(Treat =
## gl(2,nTime*(nBlock/2),labels=c('A','B'))) Xmat <- model.matrix(~-1+Block + Treat*Time, data=dt)
## block.effects <- rnorm(n = nBlock, mean = intercept, sd = sigma.block) #A.effects <- c(30,40)
## all.effects <- c(block.effects,treatment,slope,0) lin.pred <- Xmat %*% all.effects cmat <-
## sigma*rho^abs(outer(0:(nTime-1),0:(nTime-1),'-')) errs <-
## MASS::mvrnorm(nBlock,mu=rep(0,nTime),Sigma=cmat^2) c(t(errs)) y = lin.pred + c(t(errs)) data.rm <-
## data.frame(y=y, Block=paste0('B',dt$Block), Treat=dt$Treat, Time=dt$Time)

## ggplot(data.rm, aes(y=y, x=Time, color=Treat)) + geom_smooth(method='lm') + geom_point() +
## facet_wrap(~Block)


## anova(lme(y ~ 1, random = ~1 | Block, data.rm, method = 'REML'),lme(y ~ 1, random = ~Time | Block,
## data.rm, method = 'REML'))
data.rm.lme <- lme(y ~ Time, random = ~1 | Block, data.rm, method = "REML")
# random intercept/slope
data.rm.lme1 <- lme(y ~ Time, random = ~Time | Block, data.rm, method = "REML")
anova(data.rm.lme, data.rm.lme1)

             Model df      AIC      BIC    logLik   Test  L.Ratio p-value
data.rm.lme      1  4 3468.512 3483.804 -1730.256                        
data.rm.lme1     2  6 3393.773 3416.711 -1690.886 1 vs 2 78.73915  <.0001

summary(data.rm.lme)

Linear mixed-effects model fit by REML
 Data: data.rm 
       AIC      BIC    logLik
  3468.512 3483.804 -1730.256

Random effects:
 Formula: ~1 | Block
        (Intercept) Residual
StdDev:    36.62552 35.12878

Fixed effects: y ~ Time 
            Value Std.Error  DF t-value p-value
(Intercept)   200  7.509425 305 26.6332       0
Time           30  0.663280 305 45.2298       0
 Correlation: 
     (Intr)
Time -0.486

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-3.13197800 -0.57164745 -0.03286899  0.63997697  2.71992377 

Number of Observations: 340
Number of Groups: 34 

summary(data.rm.lme1)

Linear mixed-effects model fit by REML
 Data: data.rm 
       AIC      BIC    logLik
  3393.773 3416.711 -1690.886

Random effects:
 Formula: ~Time | Block
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 54.452435 (Intr)
Time         7.131168 -0.73 
Residual    27.930278       

Fixed effects: y ~ Time 
            Value Std.Error  DF  t-value p-value
(Intercept)   200  9.895208 305 20.21180       0
Time           30  1.331842 305 22.52519       0
 Correlation: 
     (Intr)
Time -0.748

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.97489572 -0.60719333  0.09430909  0.63735615  2.59035786 

Number of Observations: 340
Number of Groups: 34 

plot(ACF(data.rm.lme1, resType = "normalized"), alpha = 0.05)

data.rm.lme3 = update(data.rm.lme1, correlation = corAR1(form = ~Time | Block))
plot(ACF(data.rm.lme3, resType = "normalized"), alpha = 0.05)

summary(data.rm.lme3)

Linear mixed-effects model fit by REML
 Data: data.rm 
       AIC      BIC    logLik
  3312.423 3339.184 -1649.211

Random effects:
 Formula: ~Time | Block
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev     Corr  
(Intercept) 15.1253205 (Intr)
Time         0.5839737 -0.096
Residual    48.8427118       

Correlation Structure: AR(1)
 Formula: ~Time | Block 
 Parameter estimate(s):
      Phi 
0.7959455 
Fixed effects: y ~ Time 
            Value Std.Error  DF  t-value p-value
(Intercept)   200  9.311507 305 21.47880       0
Time           30  1.227616 305 24.43762       0
 Correlation: 
     (Intr)
Time -0.722

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.63891387 -0.69032054  0.08288653  0.68842633  2.39784979 

Number of Observations: 340
Number of Groups: 34 

data.rm.lme4 <- lme(y ~ Time, random = ~1 | Block, data.rm, method = "REML", correlation = corAR1())
plot(ACF(data.rm.lme4, resType = "normalized"), alpha = 0.05)

summary(data.rm.lme4)

Linear mixed-effects model fit by REML
 Data: data.rm 
       AIC      BIC    logLik
  3308.424 3327.539 -1649.212

Random effects:
 Formula: ~1 | Block
        (Intercept) Residual
StdDev:    14.67827 49.02625

Correlation Structure: AR(1)
 Formula: ~1 | Block 
 Parameter estimate(s):
     Phi 
0.797383 
Fixed effects: y ~ Time 
            Value Std.Error  DF  t-value p-value
(Intercept)   200  9.324885 305 21.44798       0
Time           30  1.226774 305 24.45438       0
 Correlation: 
     (Intr)
Time -0.724

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.65722844 -0.69471076  0.08010752  0.68511193  2.39738678 

Number of Observations: 340
Number of Groups: 34 

## cmat <- sigma*rho^abs(outer(0:(nTime-1),0:(nTime-1),'-')) errs <-
## MASS::mvrnorm(nBlock,mu=rep(0,nTime),Sigma=cmat) ranef <- rnorm(nBlock,mean=0,sd=sigma.block) d <-
## data.frame(Block=paste0('B',rep(1:nBlock,each=nTime))) eta <- ranef[as.numeric(d$Block)] +
## c(t(errs)) ## unpack errors by row eta = eta + mu <- identity(eta) d$y <- identity(mu) d$Time <-
## factor(rep(1:nTime,nBlock)) data.rm = d %>% mutate(Treat =
## gl(nTreat,nTime*(nBlock/nTreat),labels=c('A','B')), Time=as.numeric(as.character(Time)))


## data.rm.lme <- lme(y ~ Treat * Time, random = ~1 | Block, data.rm, method = 'REML') # random
## intercept/slope data.rm.lme1 <- lme(y ~ Treat * Time, random = ~Time | Block, data.rm, method =
## 'REML') data.rm.lme2 <- lme(y ~ Treat * Time, random = ~Treat * Time | Block, data.rm, method =
## 'REML') anova(data.rm.lme, data.rm.lme1, data.rm.lme2)



## dt <- expand.grid(Time=Time,Block=Block) %>% mutate(Treat =
## gl(2,nTime*(nBlock/2),labels=c('A','B'))) Xmat <- model.matrix(~-1+Block + Treat*Time, data=dt)
## block.effects <- rnorm(n = nBlock, mean = intercept, sd = sigma.block) #A.effects <- c(30,40)
## all.effects <- c(block.effects,0,slope,treatment) lin.pred <- Xmat %*% all.effects

## y <- rnorm(n,lin.pred,sigma) data.rm <- data.frame(y=y, Block=paste0('B',dt$Block), Treat=dt$Treat,
## Time=dt$Time) data.rm = data.rm %>% group_by(Block,Treat) %>% do({ x=.  eps=NULL eps=rho* eps[1]=0

## })

## ## the quadrat observations (within sites) are drawn from ## normal distributions with means
## according to the site means ## and standard deviations of 5 eps <- NULL eps[1] <- 0 for (j in 2:n)
## { eps[j] <- rho*eps[j-1] #residuals } y <- rnorm(n,lin.pred,sigma)+eps

## #OR eps <- NULL # first value cant be autocorrelated eps[1] <- rnorm(1,0,sigma) for (j in 2:n) {
## eps[j] <- rho*eps[j-1] + rnorm(1, mean = 0, sd = sigma) #residuals } y <- lin.pred + eps data.rm <-
## data.frame(y=y, Block=paste0('B',dt$Block), Treat=dt$Treat, Time=dt$Time)

head(data.rm)  #print out the first six rows of the data set

  Block Time        y
1    B1    1 167.8633
2    B1    2 165.8682
3    B1    3 178.6709
4    B1    4 239.3880
5    B1    5 279.6270
6    B1    6 285.1919

write.table(data.rm, file = "../downloads/data/data.rm.csv", quote = FALSE, row.names = FALSE, sep = ",")
library(ggplot2)
ggplot(data.rm, aes(y = y, x = Time)) + geom_smooth(method = "lm") + geom_point() + facet_wrap(~Block)

Exploratory data analysis

Normality and Homogeneity of variance

# ggplot(data.rm) + geom_boxplot(aes(y=y, x=Treat)) ggplot(data.rm) + geom_boxplot(aes(y=y,
# x=as.factor(Time), color=Treat))
ggplot(data.rm) + geom_boxplot(aes(y = y, x = as.factor(Time)))

Conclusions:

• there is no evidence that the response variable is consistently non-normal across all populations - each boxplot is approximately symmetrical
• there is no evidence that variance (as estimated by the height of the boxplots) differs substantially between all populations. . More importantly, there is no evidence of a relationship between mean and variance - the height of boxplots does not increase with increasing position along the y-axis. Hence it there is no evidence of non-homogeneity

• transform the scale of the response variables (to address normality etc). Note transformations should be applied to the entire response variable (not just those populations that are skewed).

Block by within-Block interaction

Explore whether random intercept is adequate, or whether it is necessary to employ random intercept and random slope. If the slopes differ dramatically between different blocks ( i.e. there are block by within-block interactions), we might want to consider random slopes..

library(car)
with(data.rm, interaction.plot(Time, Block, y))

ggplot(data.rm, aes(y = y, x = Time, color = Block, group = Block)) + geom_line() + guides(color = guide_legend(ncol = 3))

library(car)
# residualPlots(lm(y~Block+Treat*Time, data.rm))
residualPlots(lm(y ~ Block + Time, data.rm))

           Test stat Pr(>|t|)
Block             NA       NA
Time           0.000    1.000
Tukey test    -0.099    0.921

# the Tukey's non-additivity test by itself can be obtained via an internal function within the car
# package car:::tukeyNonaddTest(lm(y~Block+Treat*Time, data.rm))
car:::tukeyNonaddTest(lm(y ~ Block + Time, data.rm))

       Test      Pvalue 
-0.09912413  0.92103972 

Conclusions:

• there is no visual or inferential evidence of any major interactions between Block and the within-Block effect (Time). Any trends appear to be reasonably consistent between Blocks.
• it is therefore unlikely that we are going to need to model in random slopes.

Model fitting or statistical analysis

One of the predictor variables is Time. The order of Time treatment levels cannot be randomized and therefore there may be dependency issues. Observations that are closer together in time are likely to be less variable than observations spaced further appart temporally. That is, the data might be temporally autocorrelated.

As with the previous example, there are numerous ways of fitting a repeated measures model. However, if we need to account for temporal autocorrelation, our options are more limited.

Importantly, random intercept/slope models and random intercept models that also incorporate some form of serial dependency structure (such as a first order autoregressive) both allow the within Block dependency structure to vary over time. Once random noise is added ontop of this, it can be very difficult to pick up the subtle differences between these two temporal influences. It is even more difficult to tease the two appart.

library(nlme)
# random intercept data.rm.lme <- lme(y~Treat*Time, random=~1|Block, data.rm, method='REML')
data.rm.lme <- lme(y ~ Time, random = ~1 | Block, data.rm, method = "REML")
# random intercept/slope data.rm.lme1 <- lme(y~Treat*Time, random=~Time|Block, data.rm,
# method='REML') data.rm.lme2 <- lme(y~Treat*Time, random=~Treat*Time|Block, data.rm, method='REML')
data.rm.lme1 <- lme(y ~ Time, random = ~Time | Block, data.rm, method = "REML")
# anova(data.rm.lme, data.rm.lme1, data.rm.lme2)
anova(data.rm.lme, data.rm.lme1)

             Model df      AIC      BIC    logLik   Test  L.Ratio p-value
data.rm.lme      1  4 3468.512 3483.804 -1730.256                        
data.rm.lme1     2  6 3393.773 3416.711 -1690.886 1 vs 2 78.73915  <.0001

library(lme4)
data.rm.lmer <- lmer(y ~ Time + (1 | Block), data.rm, REML = TRUE)  #random intercept
data.rm.lmer1 <- lmer(y ~ Time + (Time | Block), data.rm, REML = TRUE,
    control = lmerControl(check.nobs.vs.nRE = "ignore"))  #random intercept/slope
anova(data.rm.lmer, data.rm.lmer1)

Data: data.rm
Models:
object: y ~ Time + (1 | Block)
..1: y ~ Time + (Time | Block)
       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
object  4 3475.1 3490.4 -1733.6   3467.1                             
..1     6 3401.8 3424.7 -1694.9   3389.8 77.358      2  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

library(glmmTMB)
data.rm.glmmTMB <- glmmTMB(y ~ Time + (1 | Block), data.rm)  #random intercept
data.rm.glmmTMB1 <- glmmTMB(y ~ Time + (Time | Block), data.rm)  #random intercept/slope
anova(data.rm.glmmTMB, data.rm.glmmTMB1)

Data: data.rm
Models:
data.rm.glmmTMB: y ~ Time + (1 | Block), zi=~0, disp=~1
data.rm.glmmTMB1: y ~ Time + (Time | Block), zi=~0, disp=~1
                 Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
data.rm.glmmTMB   4 3475.1 3490.4 -1733.6   3467.1                             
data.rm.glmmTMB1  6 3401.8 3424.7 -1694.9   3389.8 77.358      2  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

data.rm.aov <- aov(y ~ Treat * Time + Error(Block), data.rm)

Error in eval(expr, envir, enclos): object 'Treat' not found

Model evaluation

Temporal autocorrelation

Before proceeding any further we should really explore whether we are likely to have an issue with temporal autocorrelation. The models assume that there are no temporal dependency issues. A good way to explore this is to examine the autocorrelation function. Essentially, this involves looking at the degree of correlation between residuals associated with times of incrementally greater temporal lags.

We have already observed that a model with random intercepts and random slopes fits better than a model with just random intercepts. It is possible that this is due to temporal autocorrelation. The random intercept/slope model might have fit the temporally autocorrelated data better, but if this is due to autocorrelation, then the random intercepts/slope model does not actually adress the underlying issue. Consequently, it is important to explore autocorrelation for both models and if there is any evidence of temporal autocorrelation, refit both models.

We can visualize the issue via linear mixed model formulation: $$ \begin{align} y_i &\sim{} N(\mu_i, \sigma^2)\\ \mu_i &= \mathbf{X}\boldsymbol{\beta} + \mathbf{Z}\mathbf{b}\\ \mathbf{b} &\sim{} MVN(0, \Sigma)\\ \end{align} $$ With a bit or rearranging, $\mathbf{b}$ and $\Sigma$ can represent a combination of random intercepts ($\mathbf{b}_0$) and autocorrelated residuals ($\Sigma_{AR}$): $$ \begin{align} b_{ij} &= \mathbf{b}_{0,ij} + \varepsilon_{ij}\\ \varepsilon_i &\sim{} \Sigma_{AR}\\ \end{align} $$ where $i$ are the blocks and $j$ are the observations.

Residuals

As always, exploring the residuals can reveal issues of heteroscadacity, non-linearity and potential issues with autocorrelation. Note for lme() and lmer() residual plots use standardized (normalized) residuals rather than raw residuals as the former reflect changes to the variance-covariance matrix whereas the later do not.

The following function will be used for the production of some of the qqnormal plots.

qq.line = function(x) {
    # following four lines from base R's qqline()
    y <- quantile(x[!is.na(x)], c(0.25, 0.75))
    x <- qnorm(c(0.25, 0.75))
    slope <- diff(y)/diff(x)
    int <- y[1L] - slope * x[1L]
    return(c(int = int, slope = slope))
}

• lme
• lmer
• glmmTMB

plot(data.rm.lme3)

qqnorm(resid(data.rm.lme3))
qqline(resid(data.rm.lme3))

library(sjPlot)
plot_grid(plot_model(data.rm.lme3, type = "diag"))

plot(data.rm.lmer3)

plot(fitted(data.rm.lmer3), residuals(data.rm.lmer3, type = "pearson",
    scaled = TRUE))

ggplot(data.rm %>% cbind(do.call("cbind", fortify(data.rm.lmer3))),
    aes(y = .scresid, x = .fitted)) + geom_point()

QQline = qq.line(as.vector(fortify(data.rm.lmer3)$.scresid))
ggplot(data.rm %>% cbind(do.call("cbind", fortify(data.rm.lmer3))),
    aes(sample = .scresid)) + stat_qq() + geom_abline(intercept = QQline[1],
    slope = QQline[2])

qqnorm(resid(data.rm.lmer3))
qqline(resid(data.rm.lmer3))

ggplot(data = NULL, aes(y = resid(data.rm.glmmTMB3, type = "pearson"),
    x = fitted(data.rm.glmmTMB3))) + geom_point()

QQline = qq.line(resid(data.rm.glmmTMB3, type = "pearson"))
ggplot(data = NULL, aes(sample = resid(data.rm.glmmTMB3, type = "pearson"))) +
    stat_qq() + geom_abline(intercept = QQline[1], slope = QQline[2])

par(mfrow = c(2, 2))
plot(lm(data.rm.aov))

Error in stats::model.frame(formula = data.rm.aov, drop.unused.levels = TRUE): object 'data.rm.aov' not found

Conclusions: there are no issues obvious from the residuals.

Exploring model parameters

If there was any evidence that the assumptions had been violated, then we would need to reconsider the model and start the process again. In this case, there is no evidence that the test will be unreliable so we can proceed to explore the test statistics. As I had elected to illustrate multiple techniques for analysing this nested design, I will also deal with the summaries etc separately.

Partial effects plots

It is often useful to visualize partial effects plots while exploring the parameter estimates. Having a graphical representation of the partial effects typically makes it a lot easier to interpret the parameter estimates and inferences.

• lme
• lmer
• glmmTMB

library(effects)
plot(allEffects(data.rm.lme3, residuals = TRUE))

library(sjPlot)
plot_model(data.rm.lme3, type = "eff", terms = "Time")

Unfortunately, these don't work for the hack AR1 structure..

library(effects)
plot(allEffects(data.rm.lmer3))

Error in Effect(predictors, mod, vcov. = vcov., ...): model formula should not contain calls to
  factor(), as.factor(), ordered(), as.ordered(), as.numeric(), or as.integer();
  see 'Warnings and Limitations' in ?Effect

library(sjPlot)
plot_model(data.rm.lmer3, type = "eff", terms = "A")

Error in effects::Effect(focal.predictors = terms, mod = model, xlevels = xl, : model formula should not contain calls to
  factor(), as.factor(), ordered(), as.ordered(), as.numeric(), or as.integer();
  see 'Warnings and Limitations' in ?Effect

## but we can do it manually
newdata = with(data.rm, data.frame(Time = seq(min(Time), max(Time),
    len = 100)))
Xmat = model.matrix(~Time, data = newdata)
coefs = fixef(data.rm.glmmTMB3)[[1]]
fit = as.vector(coefs %*% t(Xmat))
se = sqrt(diag(Xmat %*% vcov(data.rm.glmmTMB3)[[1]] %*% t(Xmat)))
q = qt(0.975, df = df.residual(data.rm.glmmTMB3))
newdata = cbind(newdata, data.frame(fit = fit, lower = fit - q * se,
    upper = fit + q * se))
ggplot(newdata, aes(y = fit, x = Time)) + geom_line() + geom_ribbon(aes(ymin = lower,
    ymax = upper), fill = "blue", alpha = 0.3, color = NA)

Parameter estimates

summary(data.rm.lme3)

Linear mixed-effects model fit by REML
 Data: data.rm 
       AIC      BIC    logLik
  3308.424 3327.539 -1649.212

Random effects:
 Formula: ~1 | Block
        (Intercept) Residual
StdDev:    14.67827 49.02625

Correlation Structure: AR(1)
 Formula: ~Time | Block 
 Parameter estimate(s):
     Phi 
0.797383 
Fixed effects: y ~ Time 
            Value Std.Error  DF  t-value p-value
(Intercept)   200  9.324885 305 21.44798       0
Time           30  1.226774 305 24.45438       0
 Correlation: 
     (Intr)
Time -0.724

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.65722844 -0.69471076  0.08010752  0.68511193  2.39738678 

Number of Observations: 340
Number of Groups: 34 

intervals(data.rm.lme3)

Approximate 95% confidence intervals

 Fixed effects:
                lower est.     upper
(Intercept) 181.65075  200 218.34925
Time         27.58599   30  32.41401
attr(,"label")
[1] "Fixed effects:"

 Random Effects:
  Level: Block 
                    lower     est.    upper
sd((Intercept)) 0.8372059 14.67827 257.3459

 Correlation structure:
        lower     est.     upper
Phi 0.6484983 0.797383 0.8875083
attr(,"label")
[1] "Correlation structure:"

 Within-group standard error:
   lower     est.    upper 
36.93680 49.02625 65.07260 

anova(data.rm.lme3)

            numDF denDF  F-value p-value
(Intercept)     1   305 3215.819  <.0001
Time            1   305  598.017  <.0001

library(broom)
tidy(data.rm.lme3, effects = "fixed", conf.int = TRUE)

# A tibble: 2 x 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)    200.       9.32      21.4 7.33e-63
2 Time            30.0      1.23      24.5 7.25e-74

glance(data.rm.lme3)

# A tibble: 1 x 5
  sigma logLik   AIC   BIC deviance
  <dbl>  <dbl> <dbl> <dbl> <lgl>   
1  49.0 -1649. 3308. 3328. NA      

summary(data.rm.lmer3)

Linear mixed model fit by REML ['lmerMod']

REML criterion at convergence: 3298.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-0.39534 -0.08754  0.00350  0.07771  0.33821 

Random effects:
 Groups   Name        Variance Std.Dev. Corr                                        
 Block    (Intercept)  185.051 13.603                                               
 Block.1  Time1       2423.592 49.230                                               
          Time2       2423.592 49.230   0.80                                        
          Time3       2423.592 49.230   0.64 0.80                                   
          Time4       2423.592 49.230   0.52 0.64 0.80                              
          Time5       2423.592 49.230   0.42 0.52 0.64 0.80                         
          Time6       2423.592 49.230   0.33 0.42 0.52 0.64 0.80                    
          Time7       2423.592 49.230   0.27 0.33 0.42 0.52 0.64 0.80               
          Time8       2423.592 49.230   0.21 0.27 0.33 0.42 0.52 0.64 0.80          
          Time9       2423.592 49.230   0.17 0.21 0.27 0.33 0.42 0.52 0.64 0.80     
          Time10      2423.592 49.230   0.14 0.17 0.21 0.27 0.33 0.42 0.52 0.64 0.80
 Residual                9.209  3.035                                               
Number of obs: 340, groups:  Block, 34

Fixed effects:
                 Estimate Std. Error t value
(Intercept)       200.000      9.332   21.43
as.numeric(Time)   30.000      1.229   24.41

Correlation of Fixed Effects:
            (Intr)
as.nmrc(Tm) -0.724

confint(data.rm.lmer3)

Error in getOptfun(optimizer): couldn't find optimizer function 

anova(data.rm.lmer3)

Analysis of Variance Table
                 Df Sum Sq Mean Sq F value
as.numeric(Time)  1 5489.2  5489.2  596.09

library(broom)
tidy(data.rm.lmer3, effects = "fixed", conf.int = TRUE)

# A tibble: 2 x 6
  term             estimate std.error statistic conf.low conf.high
  <chr>               <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
1 (Intercept)         200.       9.33      21.4    182.      218. 
2 as.numeric(Time)     30.0      1.23      24.4     27.6      32.4

glance(data.rm.lmer3)

# A tibble: 1 x 6
  sigma logLik   AIC   BIC deviance df.residual
  <dbl>  <dbl> <dbl> <dbl>    <dbl>       <int>
1  3.03 -1649. 3416. 3642.    3306.         281

summary(data.rm.glmmTMB3)

 Family: gaussian  ( identity )
Formula:          y ~ Time + (1 | Block) + ar1(as.factor(Time) - 1 | Block)
Data: data.rm

     AIC      BIC   logLik deviance df.resid 
      NA       NA       NA       NA      334 

Random effects:

Conditional model:
 Groups   Name             Variance  Std.Dev. Corr      
 Block    (Intercept)      2.254e+02 15.015             
 Block.1  as.factor(Time)1 2.334e+03 48.314   0.79 (ar1)
 Residual                  3.240e-04  0.018             
Number of obs: 340, groups:  Block, 34

Dispersion estimate for gaussian family (sigma^2): 0.000324 

Conditional model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  199.999      9.207   21.72   <2e-16 ***
Time          30.000      1.214   24.71   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

confint(data.rm.glmmTMB3)

Error in is.nan(x): default method not implemented for type 'list'

$R^2$ approximations

library(MuMIn)
r.squaredGLMM(data.rm.lme3)

      R2m       R2c 
0.7398128 0.7612169 

library(MuMIn)
r.squaredGLMM(data.rm.lmer3)

Error in r.squaredGLMM.merMod(data.rm.lmer3): random term names do not match those in model matrix. 
Have 'options(contrasts)' changed since the model was fitted?

source(system.file("misc/rsqglmm.R", package = "glmmTMB"))
my_rsq(data.rm.lmer3)

$family
[1] "gaussian"

$link
[1] "identity"

$Marginal
[1] 0.9745773

$Conditional
[1] 0.9987948

source(system.file("misc/rsqglmm.R", package = "glmmTMB"))
my_rsq(data.rm.glmmTMB3)

$family
[1] "gaussian"

$link
[1] "identity"

$Marginal
[1] 0.970617

$Conditional
[1] 1

Graphical summary

It is relatively trivial to produce a summary figure based on the raw data. Arguably a more satisfying figure would be one based on the modelled data.

• lme
• lmer
• glmmTMB

The new effects package seems to have a bug when calculating effects

library(effects)
data.rm.eff = as.data.frame(allEffects(data.rm.lme3, xlevels = 100)[[1]])
# OR
data.rm.eff = as.data.frame(Effect("Time", data.rm.lme3, xlevels = 100))

ggplot(data.rm.eff, aes(y = fit, x = Time)) + geom_line() + geom_ribbon(aes(ymin = lower,
    ymax = upper), fill = "blue", color = NA, alpha = 0.3) +
    scale_y_continuous("Y") + theme_classic()

# OR using lsmeans
at = list(Time = seq(min(data.rm$Time), max(data.rm$Time), len = 100))
fit = summary(ref.grid(data.rm.lme3, at = at), infer = TRUE)
ggplot(fit, aes(y = prediction, x = Time)) + geom_line() + geom_ribbon(aes(ymin = lower.CL,
    ymax = upper.CL), fill = "blue", color = NA, alpha = 0.3) +
    scale_y_continuous("Y") + theme_classic()

# OR
at = list(Time = seq(min(data.rm$Time), max(data.rm$Time), len = 100))
fit = summary(lsmeans(data.rm.lme3, eff ~ Time, at = at)$lsmean)
ggplot(fit, aes(y = lsmean, x = Time)) + geom_line() + geom_ribbon(aes(ymin = lower.CL,
    ymax = upper.CL), fill = "blue", color = NA, alpha = 0.3) +
    scale_y_continuous("Y") + theme_classic()

Or manually

newdata = with(data.rm, data.frame(Time = seq(min(Time), max(Time),
    len = 100)))
Xmat = model.matrix(~Time, data = newdata)
coefs = fixef(data.rm.lme3)
fit = as.vector(coefs %*% t(Xmat))
se = sqrt(diag(Xmat %*% vcov(data.rm.lme3) %*% t(Xmat)))
q = qt(0.975, df = nrow(data.rm.lme3$data) - length(coefs) -
    2)
newdata = cbind(newdata, fit = fit, lower = fit - q * se, upper = fit +
    q * se)
ggplot(newdata, aes(y = fit, x = Time)) + geom_line() + geom_ribbon(aes(ymin = lower,
    ymax = upper), fill = "blue", color = NA, alpha = 0.3) +
    scale_y_continuous("Y") + theme_classic()

It seems that manually is our only option.

newdata = with(data.rm, data.frame(Time = seq(min(Time), max(Time),
    len = 100)))
Xmat = model.matrix(~Time, data = newdata)
coefs = fixef(data.rm.lmer3)
fit = as.vector(coefs %*% t(Xmat))
se = sqrt(diag(Xmat %*% vcov(data.rm.lmer3) %*% t(Xmat)))
q = qt(0.975, df = df.residual(data.rm.lmer3))
newdata = cbind(newdata, fit = fit, lower = fit - q * se, upper = fit +
    q * se)
ggplot(newdata, aes(y = fit, x = Time)) + geom_line() + geom_ribbon(aes(ymin = lower,
    ymax = upper), fill = "blue", color = NA, alpha = 0.3) +
    scale_y_continuous("Y") + theme_classic()

Manually

newdata = with(data.rm, data.frame(Time = seq(min(Time), max(Time),
    len = 100)))
Xmat = model.matrix(~Time, data = newdata)
coefs = fixef(data.rm.glmmTMB3)[[1]]
fit = as.vector(coefs %*% t(Xmat))
se = sqrt(diag(Xmat %*% vcov(data.rm.glmmTMB3)[[1]] %*% t(Xmat)))
q = qt(0.975, df = df.residual(data.rm.glmmTMB3))
newdata = cbind(newdata, fit = fit, lower = fit - q * se, upper = fit +
    q * se)
ggplot(newdata, aes(y = fit, x = Time)) + geom_line() + geom_ribbon(aes(ymin = lower,
    ymax = upper), fill = "blue", color = NA, alpha = 0.3) +
    scale_y_continuous("Y") + theme_classic()

RCB (repeated measures) ANOVA in R - categoical and continuous within and AR1 temporal autocorrelation

Scenario and Data

Imagine now that we has designed an experiment to investigate the effects of a continuous predictor ($x$, for example time) on a response ($y$) and as with the previous example, the data are collected repeatedly from within a number of blocks. We can extend this example to include a treatment with three levels (A, B and C), each of which occur within each block. Thus there are now continuous and categorical predictors (and their interaction) within each block.

Thus in an attempt to constrain this variability, we again decide to apply a design (RCB) in which each of the levels of X (such as time) treatments within each of 35 blocks dispersed randomly throughout the landscape. As this section is mainly about the generation of artificial data (and not specifically about what to do with the data), understanding the actual details are optional and can be safely skipped. Consequently, I have folded (toggled) this section away.

library(tidyverse)
set.seed(1)
slope <- 30
intercept <- 200
treatmentB = 90
treatmentC = 100
nTreat = 3
nBlock <- 34
nTime <- 10
sd.res = 50
sd.block = 10
sd.slope = 0.1
cor_int_slope = 0.4
rho = 0.8
## define the random effects (intercepts and slopes)
rand.effects = MASS::mvrnorm(n = nBlock, mu = c(0, 0), Sigma = matrix(c(1, cor_int_slope, cor_int_slope,
    1), nrow = 2), empirical = TRUE) %>% as.data.frame %>% mutate(V1 = V1 * sd.block + intercept, V2 = V2 *
    sd.slope + slope)
apply(rand.effects, 2, mean)