extract_rcov {agriutilities}

Extract Residual Variance-Covariance from ASReml-R

Description

This function is specially useful when running repeated measurements models.

Usage

extract_rcov(model = NULL, time = "Time", plot = "Plot", vc_error = NULL)

Arguments

model

An asreml object

time

A character string indicating the "Time"

plot

A character string indicating the "PlotID"

vc_error

An optional character string indicating the variance covariance. Can be "corv", "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.

Examples

## Not run: 
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, function(x) summary(x)$aic)),
  bic = unlist(lapply(models, function(x) summary(x)$bic)),
  loglik = unlist(lapply(models, function(x) summary(x)$loglik)),
  nedf = unlist(lapply(models, function(x) summary(x)$nedf)),
  param = unlist(lapply(models, function(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 -----------------------------------

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

## End(Not run)