Extract Residual Variance-Covariance from ASReml-R

This function is specially useful for extracting residual variance covariance matrices from ASReml-R when running repeated measurements analysis.

Usage

extract_rcov(model = NULL, time = NULL, plot = NULL, vc_error = NULL)

Arguments

model

An asreml object.

time

An optional character string indicating the "Time". By default the function identifies this parameter.

plot

An optional character string indicating the "PlotID". By default the function identifies this parameter.

vc_error

An optional character string indicating the variance covariance. It can be "corv", "corh", "corgh", "us", "expv", "exph", "ar1v", "ar1h" or "ante". By using NULL the function tries to guess which was the variance-covariance used.

Value

An object with a list of:

corr_mat

A matrix with the residual correlation between time points.

vcov_mat

A matrix of the estimated residual variance-covariance between time points.

vc

A character string indicating the variance-covariance fitted.

Details

The expected residual variance covariance structure must be of the form: `~id(Plot):corv(Time)`, where `Plot` is a unique identifier of each experimental unit, and `Time`, represents the variable that contains the time when the experimental units were measured. This form also requires that the levels of the factor Time are nested in the levels of the factor Plot. If it is not in that form you can sort the dataset by using the following command `arrange(grassUV, Plant, Time)`.

Examples

if (FALSE) {
library(ggpubr)
library(agriutilities)
library(tidyverse)
library(asreml)

head(grassUV)
str(grassUV)

# Exploration -------------------------------------------------------------

grassUV %>%
  ggplot(
    aes(x = Time, y = y, group = Plant, color = Plant)
  ) +
  geom_point() +
  geom_line() +
  facet_wrap(~Tmt) +
  theme_minimal(base_size = 15)

tmp <- grassUV %>%
  group_by(Time, Plant) %>%
  summarise(mean = mean(y, na.rm = TRUE)) %>%
  spread(Time, mean) %>%
  column_to_rownames("Plant")

gg_cor(tmp, label_size = 5)

tmp %>%
  cor(use = "pairwise.complete.obs") %>%
  as.data.frame() %>%
  rownames_to_column(var = "Time") %>%
  gather("DAP2", "corr", -1) %>%
  type.convert(as.is = FALSE) %>%
  mutate(corr = ifelse(Time < DAP2, NA, corr)) %>%
  mutate(DAP2 = as.factor(DAP2)) %>%
  ggplot(
    aes(x = Time, y = corr, group = DAP2, color = DAP2)
  ) +
  geom_point() +
  geom_line() +
  theme_minimal(base_size = 15) +
  color_palette(palette = "jco") +
  labs(color = "Time", y = "Pearson Correlation")

# Modeling ----------------------------------------------------------------

# Identity variance model.
model_0 <- asreml(
  fixed = y ~ Time + Tmt + Tmt:Time,
  residual = ~ id(Plant):idv(Time),
  data = grassUV
)

# Simple correlation model; homogeneous variance form.
model_1 <- asreml(
  fixed = y ~ Time + Tmt + Tmt:Time,
  residual = ~ id(Plant):corv(Time),
  data = grassUV
)

# Exponential (or power) model; homogeneous variance form.
model_2 <- asreml(
  fixed = y ~ Time + Tmt + Tmt:Time,
  residual = ~ id(Plant):expv(Time),
  data = grassUV
)

# Exponential (or power) model; heterogeneous variance form.
model_3 <- asreml(
  fixed = y ~ Time + Tmt + Tmt:Time,
  residual = ~ id(Plant):exph(Time),
  data = grassUV
)

# Antedependence variance model of order 1
model_4 <- asreml(
  fixed = y ~ Time + Tmt + Tmt:Time,
  residual = ~ id(Plant):ante(Time),
  data = grassUV
)

# Autoregressive model of order 1; homogeneous variance form.
model_5 <- asreml(
  fixed = y ~ Time + Tmt + Tmt:Time,
  residual = ~ id(Plant):ar1v(Time),
  data = grassUV
)

# Autoregressive model of order 1; heterogeneous variance form.
model_6 <- asreml(
  fixed = y ~ Time + Tmt + Tmt:Time,
  residual = ~ id(Plant):ar1h(Time),
  data = grassUV
)

# Unstructured variance model.
model_7 <- asreml(
  fixed = y ~ Time + Tmt + Tmt:Time,
  residual = ~ id(Plant):us(Time),
  data = grassUV
)

# Model Comparison --------------------------------------------------------

models <- list(
  "id" = model_0,
  "cor" = model_1,
  "exp" = model_2,
  "exph" = model_3,
  "ante" = model_4,
  "ar1" = model_5,
  "ar1h" = model_6,
  "us" = model_7
)

summary_models <- data.frame(
  model = names(models),
  aic = unlist(lapply(models, \(x) summary(x)$aic)),
  bic = unlist(lapply(models, \(x) summary(x)$bic)),
  loglik = unlist(lapply(models, \(x) summary(x)$loglik)),
  nedf = unlist(lapply(models, \(x) summary(x)$nedf)),
  param = unlist(lapply(models, \(x) attr(summary(x)$aic, "param"))),
  row.names = NULL
)

summary_models %>%
  ggplot(
    aes(x = reorder(model, -bic), y = bic, group = 1)
  ) +
  geom_point(size = 2) +
  geom_text(aes(x = model, y = bic + 5, label = param)) +
  geom_line() +
  theme_minimal(base_size = 15) +
  labs(x = NULL)

# Extracting Variance Covariance Matrix -----------------------------------

extract_rcov(model_4)

covcor_heat(
  matrix = extract_rcov(model_1)$corr,
  legend = "none",
  size = 5
) + ggtitle(label = "Uniform Correlation (corv)")
covcor_heat(
  matrix = extract_rcov(model_2)$corr,
  legend = "none",
  size = 5
) + ggtitle(label = "Exponetial (expv)")
}