Repeated Measures: Covariance Structures in ASReml

Dr. Samuel B Fernandes

2026-01-01

Like before, it starts with your questions.

Define your research questions/hypothesis.

We need a null hypothesis for each question, for example:

Does irrigation affect the yield of different cultivars differently?

\(H_0:i_1cv_1=i_2cv_2= ...=i_ncv_n=0\)
\(H_a:\) at least one combination of cultivar and irrication will yield differently.

Designing an experiment to answer your questions

First of all, I will load all packages I know will used downstream.

library(tidyverse)
library(ggpubr)
library(rstatix)
library(desplot)
library(asremlPlus)
library(asreml)
#> Online License checked out Mon Mar 30 10:33:42 2026

What design should I choose?

  • Do I have homogeneous experimental units within each block?
  • Can all factors in my experiment be just as easily randomized?

You can use FielD. Lets create a split plot design and save it as “Split-Plot_2025-04-04.csv”.

library(FielDHub)
run_app()

It is time for analysis

After running the experiment and collecting data with the hightest quality possible, Loading the data and creating a “fake” split-plot with the same response variable to illustrate the split-plot analysis.

Code
alfaRCBD <- read.table("ALFALFA.txt", header = TRUE, na.string = "*")
alfaSPLIT <- read_csv("Split-Plot_2025-04-04.csv") |>
    mutate(resp = alfaRCBD$Resp) |>
    rename(irrigation = "WHOLE_PLOT", cv = "SUB_PLOT", col = "COLUMN", block = "REP") |>
    select(-ID, -LOCATION, -TRT_COMB)
colnames(alfaSPLIT) <- tolower(colnames(alfaSPLIT))

Creating factors

alfaRCBD <- transform(alfaRCBD,
    cv = as.factor(Variety),
    b = factor(Block)
)

alfaSPLIT <- transform(alfaSPLIT,
    cv = factor(cv),
    block = factor(block),
    irrigation = factor(irrigation)
)

EDA - explatoraty data analysis

Understanding the experiment

desplot(
    alfaSPLIT,
    irrigation ~ row + col,
    col = irrigation,
    cex = .5,
    num = cv,
    out1 = block,
    out2 = cv
)

hist(alfaRCBD$Resp)

boxplot(alfaRCBD$Resp)

table(alfaRCBD$Variety, alfaRCBD$Block)
#>    
#>     1 2 3 4 5 6
#>   A 1 1 1 1 1 1
#>   B 1 1 1 1 1 1
#>   C 1 1 1 1 1 1
#>   D 1 1 1 1 1 1
#>   E 1 1 1 1 1 1
#>   F 1 1 1 1 1 1
#>   G 1 1 1 1 1 1
#>   H 1 1 1 1 1 1
#>   I 1 1 1 1 1 1
#>   J 1 1 1 1 1 1
#>   K 1 1 1 1 1 1
#>   L 1 1 1 1 1 1

Mixed effects model

Analysis using ASReml - RCBD

A mixed model with cv as a random effect by all the ANOVA assumptions still hold.

alfaRCBD <- alfaRCBD %>%
    arrange(cv, b)

rcbd1 <- asreml(
    fixed = Resp ~ b,
    random = ~cv,
    # residual=<error structure>, 
    data = alfaRCBD
)
#> ASReml Version 4.2 30/03/2026 10:33:43
#>           LogLik        Sigma2     DF     wall
#>  1      48.73447    0.06197422     66   10:33:43
#>  2      50.02113    0.05731910     66   10:33:43
#>  3      51.15010    0.05255230     66   10:33:43
#>  4      51.69752    0.04874910     66   10:33:43
#>  5      51.73661    0.04775121     66   10:33:43
#>  6      51.73696    0.04765360     66   10:33:43

Diagnostics

par(mfrow = c(2, 2))
plot(rcbd1)
par(mfrow = c(1, 1))

Dropping the assumtions

Dropping the assumtion of homogeneity of variances (“diag()” estructure) across cultivars.

rcbd2 <- asreml(
    fixed = Resp ~ b,
    random = ~cv,
    residual = ~ diag(cv):id(b),
    data = alfaRCBD
)
#> ASReml Version 4.2 30/03/2026 10:33:43
#>           LogLik        Sigma2     DF     wall
#>  1      48.02611           1.0     66   10:33:43  (  1 restrained)
#>  2      52.37703           1.0     66   10:33:43  (  1 restrained)
#>  3      52.85416           1.0     66   10:33:43  (  1 restrained)
#>  4      57.36453           1.0     66   10:33:43  (  1 restrained)
#>  5      58.08428           1.0     66   10:33:43  (  1 restrained)
#>  6      58.15785           1.0     66   10:33:43
#>  7      58.16005           1.0     66   10:33:43
plot(rcbd2)

Dropping the assumtions

Dropping the assumtion of independent residuals (“ar1()” estructure) across cultivars.

rcbd3 <- asreml(
    fixed = Resp ~ b,
    random = ~cv,
    residual = ~ ar1(cv):id(b),
    data = alfaRCBD
)
#> ASReml Version 4.2 30/03/2026 10:33:43
#>           LogLik        Sigma2     DF     wall
#>  1      49.93931    0.06124234     66   10:33:43
#>  2      53.07193    0.05568745     66   10:33:43
#>  3      55.48466    0.05307777     66   10:33:43
#>  4      56.29486    0.05281485     66   10:33:43
#>  5      56.43731    0.05330333     66   10:33:43
#>  6      56.44170    0.05381629     66   10:33:43
#>  7      56.44196    0.05399564     66   10:33:43
plot(rcbd3)

Dropping the assumtions

Dropping the assumtions of independence and homogeneity of variances (“ar1()” estructure) across cultivars.

rcbd4 <- asreml(
    fixed = Resp ~ b,
    random = ~cv,
    residual = ~ ar1h(cv):id(b),
    data = alfaRCBD
)
#> ASReml Version 4.2 30/03/2026 10:33:44
#>           LogLik        Sigma2     DF     wall
#>  1      49.59856           1.0     66   10:33:44  (  1 restrained)
#>  2      57.09737           1.0     66   10:33:44  (  1 restrained)
#>  3      61.70343           1.0     66   10:33:44  (  1 restrained)
#>  4      63.34096           1.0     66   10:33:44  (  1 restrained)
#>  5      63.35612           1.0     66   10:33:44  (  1 restrained)
#>  6      63.77299           1.0     66   10:33:44  (  1 restrained)
#>  7      63.44874           1.0     66   10:33:44  (  1 restrained)
#>  8      63.97328           1.0     66   10:33:44
#>  9      63.96703           1.0     66   10:33:44
#> 10      63.63862           1.0     66   10:33:44  (  1 restrained)
#> 11      62.72847           1.0     66   10:33:44
#> 12      63.54572           1.0     66   10:33:44
#> 13      63.95277           1.0     66   10:33:44
plot(rcbd4)
rbind(
    infoCriteria(rcbd1),
    infoCriteria(rcbd2),
    infoCriteria(rcbd3),
    infoCriteria(rcbd4)
)
#>   fixedDF varDF NBound        AIC        BIC   loglik
#> 1       0     2      0  -99.47393  -95.09462 51.73696
#> 2       0    12      2  -92.32010  -66.04425 58.16005
#> 3       0     3      0 -106.88391 -100.31495 56.44196
#> 4       0    14      1  -99.90553  -69.25037 63.95277

Analysis of a Split-Plot - ASReml

alfaSPLIT <- alfaSPLIT %>%
    arrange(irrigation, cv, block, row, col)

split1 <- asreml(
    fixed = resp ~ irrigation * cv,
    random = ~ block + block:irrigation,
    residual = ~ id(units),
    data = alfaSPLIT
)
#> ASReml Version 4.2 30/03/2026 10:33:44
#>           LogLik        Sigma2     DF     wall
#>  1      11.59703     0.1189617     48   10:33:44  (  1 restrained)
#>  2      12.45475     0.1160740     48   10:33:44  (  1 restrained)
#>  3      12.72741     0.1127543     48   10:33:44  (  1 restrained)
#>  4      12.83207     0.1099238     48   10:33:44  (  1 restrained)
#>  5      12.83719     0.1092731     48   10:33:44

What if it doesn’t converge?

If the statistical model did not reach convergence, you can try the update function.

split1 <- update.asreml(split1)
#> ASReml Version 4.2 30/03/2026 10:33:44
#>           LogLik        Sigma2     DF     wall
#>  1      12.83720     0.1092318     48   10:33:44
#>  2      12.83720     0.1092318     48   10:33:44

Variance components

The function vc()from the package lucid can extract the variance components from asreml, lme4, mmer, nlme and mcmc.list objects.

lucid::vc(split1)
#>            effect  component std.error z.ratio bound %ch
#>             block 0.00000017        NA      NA     B  NA
#>  block:irrigation 0.03426      0.03073     1.1     P   0
#>           units!R 0.1092       0.02329     4.7     P   0

Predicted values (you can think of adjusted means)

Fixed effects: estimates

  • Best linear unbiased estimates (BLUEs)

Random effects: Predictions

  • Best linear unbiased predictions (BLUPs)

How do we get predicted values (you can think of adjusted means).

preds <- predict(split1, classify = "cv")$pvals
#> ASReml Version 4.2 30/03/2026 10:33:44
#>           LogLik        Sigma2     DF     wall
#>  1      12.83720     0.1092317     48   10:33:44
#>  2      12.83720     0.1092317     48   10:33:44
#>  3      12.83720     0.1092317     48   10:33:44

If you want a tukey test, you can use the library biometryassist

pred.out <- biometryassist::multiple_comparisons(split1, classify = "cv")

Plotting the results

autoplot(pred.out, label_height = 0.5)

Repeated Measures

grassUV <- read.csv("grass.csv")

grassUV <- transform(
    grassUV,
    Tmt = factor(Tmt),
    Plant = factor(Plant),
    Time = factor(Time),
    HeightID = factor(HeightID)
)

ggplot(grassUV, aes(x = Time, y = y, group = Plant, colour = Plant)) +
    geom_line() +
    facet_grid(. ~ Tmt)

identity

id_asreml <- asreml(y ~ Tmt * Time, data = grassUV)
#> ASReml Version 4.2 30/03/2026 10:33:45
#>           LogLik        Sigma2     DF     wall
#>  1     -209.4418      286.3099     60   10:33:45

compound symetry

cs_asreml <- asreml(y ~ Tmt * Time,
    residual = ~ Plant:cor(Time),
    data = grassUV,
    trace = FALSE
)

residual variance from asreml

summary(cs_asreml)$varcomp[1, 1]
#> [1] 286.3089

correlation (between measures taken on the same plants in multiple times) from asreml

summary(cs_asreml)$varcomp[2, 1]
#> [1] 0.5581911

AR1

ar1_asreml <- asreml(y ~ Tmt * Time,
    residual = ~ Plant:ar1(Time),
    data = grassUV
)
#> ASReml Version 4.2 30/03/2026 10:33:45
#>           LogLik        Sigma2     DF     wall
#>  1     -205.4312      252.5043     60   10:33:45  (  1 restrained)
#>  2     -193.2995      192.1892     60   10:33:45  (  1 restrained)
#>  3     -182.9788      203.0901     60   10:33:45
#>  4     -181.3071      350.9370     60   10:33:45
#>  5     -181.0678      280.9927     60   10:33:45
#>  6     -181.0631      289.0975     60   10:33:45
#>  7     -181.0630      287.7454     60   10:33:45

Diagonal

diag_asreml <- asreml(y ~ Tmt * Time,
    residual = ~ Plant:diag(Time),
    data = grassUV
)
#> ASReml Version 4.2 30/03/2026 10:33:45
#>           LogLik        Sigma2     DF     wall
#>  1     -194.8709           1.0     60   10:33:45
#>  2     -193.9222           1.0     60   10:33:45
#>  3     -193.1067           1.0     60   10:33:45
#>  4     -192.7882           1.0     60   10:33:45
#>  5     -192.7819           1.0     60   10:33:45

Unstructured

us_asreml <- asreml(y ~ Tmt * Time,
    residual = ~ Plant:us(Time),
    data = grassUV,
    trace = FALSE
)
rbind(
    infoCriteria(id_asreml),
    infoCriteria(cs_asreml),
    infoCriteria(ar1_asreml),
    infoCriteria(diag_asreml),
    infoCriteria(us_asreml)
)
#>   fixedDF varDF NBound      AIC      BIC    loglik
#> 1       0     1      0 420.8836 422.9779 -209.4418
#> 2       0     2      0 397.7535 401.9422 -196.8768
#> 3       0     2      0 366.1259 370.3146 -181.0630
#> 4       0     5      1 395.5637 406.0354 -192.7819
#> 5       0    15      1 346.0712 377.4863 -158.0356

Wald test for fixed effects

wald(id_asreml)
#> Wald tests for fixed effects.
#> Response: y
#> Terms added sequentially; adjusted for those above.
#> 
#>               Df Sum of Sq Wald statistic Pr(Chisq)    
#> (Intercept)    1    305303        1066.34 < 2.2e-16 ***
#> Tmt            1      8707          30.41 3.495e-08 ***
#> Time           4     27554          96.24 < 2.2e-16 ***
#> Tmt:Time       4      2581           9.02   0.06072 .  
#> residual (MS)          286                             
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(aov(y ~ Tmt * Time, data = grassUV)) # this option uses a Chisq2 test
#> Anova Table (Type II tests)
#> 
#> Response: y
#>            Sum Sq Df F value    Pr(>F)    
#> Tmt        8707.0  1 30.4112 7.810e-07 ***
#> Time      27554.1  4 24.0597 6.614e-12 ***
#> Tmt:Time   2581.1  4  2.2538    0.0738 .  
#> Residuals 17178.6 60                      
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Spatial adjustment

The model below won’t work because of the order the data was arranged.

split2 <- asreml(
    fixed = resp ~ irrigation * cv,
    random = ~ block + block:irrigation,
    residual = ~ row:col,
    data = alfaSPLIT
)

Arranging the data by row and column

alfaSPLIT <- alfaSPLIT %>%
    arrange(row, col)

split2 <- asreml(
    fixed = resp ~ irrigation * cv,
    random = ~ block + block:irrigation,
    residual = ~ row:col,
    data = alfaSPLIT,
    trace = F
)

Using an AR1 structure for rows and columns

alfaSPLIT <- alfaSPLIT %>%
    arrange(row, col) |>
    mutate(row = factor(row), col = factor(col))

split3 <- asreml(
    fixed = resp ~ irrigation * cv,
    random = ~ block + block:irrigation,
    residual = ~ ar1(row):ar1(col),
    data = alfaSPLIT,
    trace = F
)

Comparing models with and without spatial adjustment

infoCriteria(split1)
#>   fixedDF varDF NBound       AIC       BIC  loglik
#> 1       0     2      1 -21.67441 -17.93201 12.8372
infoCriteria(split3)
#>   fixedDF varDF NBound      AIC       BIC   loglik
#> 1       0     4      1 -27.7667 -20.28189 17.88335

Pvalues for fixed effects

wald(split1)
#> Wald tests for fixed effects.
#> Response: resp
#> Terms added sequentially; adjusted for those above.
#> 
#>               Df Sum of Sq Wald statistic Pr(Chisq)    
#> (Intercept)    1    38.534         352.78    <2e-16 ***
#> irrigation     1     0.192           1.76    0.1852    
#> cv            11     0.562           5.15    0.9237    
#> irrigation:cv 11     0.758           6.94    0.8040    
#> residual (MS)        0.109                             
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
wald(split3)
#> Wald tests for fixed effects.
#> Response: resp
#> Terms added sequentially; adjusted for those above.
#> 
#>               Df Sum of Sq Wald statistic Pr(Chisq)    
#> (Intercept)    1   24.9437        218.621    <2e-16 ***
#> irrigation     1    0.1667          1.461    0.2268    
#> cv            11    0.8835          7.743    0.7361    
#> irrigation:cv 11    1.2929         11.332    0.4159    
#> residual (MS)       0.1141                             
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1