Nothing
cat("#### Test estimateV at, fa, rr functions for MET data with asreml42\n")
test_that("MET_estimateV_asreml42", {
skip_if_not_installed("asreml")
skip_on_cran()
library(dae)
library(asreml)
library(asremlPlus)
#Function too replace design matrix colnames when computing V
replaceColnames <- function(des)
{
colnam <- colnames(des)
colnam <- strsplit(colnam, split = ":")
colnam <- lapply(colnam, strsplit, split = "_")
colnam <- lapply(colnam,
function(nam)
paste(sapply(nam, function(fac) fac[1]), collapse = ":"))
colnames(des) <- colnam
return(des)
}
data(MET)
asreml.options(design = TRUE, keep.order=TRUE)
#Test at
asreml.obj <-asreml(fixed = GY.tha ~ + at(expt, c(1:5)):rep + at(expt, c(1)):vrow +
at(expt, c(2,3,6,7)):colblocks +
at(expt, c(1:5,7)):vcol + Genotype*Condition*expt,
random = ~ at(expt, c(1)):dev(vrow) + at(expt, c(2)):spl(vcol) +
at(expt, c(3,5,7)):dev(vcol) + at(expt, c(7)):units,
data=comb.dat, maxit = 100, workspace = "1Gb")
summary(asreml.obj)$varcomp
ranterms <- names(asreml.obj$G.param)
n <- nrow(comb.dat)
design <- asreml.obj$design
design <- replaceColnames(design)
V.g <- matrix(0, nrow = n, ncol = n)
for (term in ranterms)
{
cols <- grep(term, colnames(design), fixed = TRUE)
V.g <- V.g + asreml.obj$vparameters[term] * (design[, cols] %*%
t(as.matrix(design[, cols])))
}
V.g <- asreml.obj$sigma2 * (V.g + mat.I(n))
Vat <- estimateV(asreml.obj)
testthat::expect_true(all(abs(Vat - V.g) < 1e-06))
#Test fa
asreml.obj <-asreml(fixed = GY.tha ~ + at(expt, c(1:5)):rep + at(expt, c(1)):vrow +
at(expt, c(2,3,6,7)):colblocks +
at(expt, c(1:5,7)):vcol + Condition*expt,
random = ~ fa(exptCond, k = 2):Genotype +
at(expt, c(1)):dev(vrow) + at(expt, c(2)):spl(vcol) +
at(expt, c(3,5,7)):dev(vcol) + at(expt, c(7)):units,
data=comb.dat, maxit = 100, workspace = "1Gb")
summary(asreml.obj)$varcomp
ranterms <- names(asreml.obj$G.param)
n <- nrow(comb.dat)
V.g <- matrix(0, nrow = n, ncol = n)
# for (term in ranterms[2:7])
# {
# cols <- grep(term, colnames(design), fixed = TRUE)
# print(length(cols))
# V.g <- V.g + asreml.obj$vparameters[term] * (asreml.obj$design[, cols] %*%
# t(as.matrix(asreml.obj$design[, cols])))
# }
# term <- ranterms[1]
# term <- "fa(exptCond, k = 2)"
# cols <- grep(term, colnames(asreml.obj$design), fixed = TRUE)[1:1364]
# vp <- asreml.obj$vparameters[names(asreml.obj$vparameters)[grep(term,
# names(asreml.obj$vparameters), fixed = TRUE)]]
# spec.var <- diag(vp[grepl("!var", names(vp), fixed = TRUE)])
# loads <- matrix(vp[grepl("!fa", names(vp), fixed = TRUE)], ncol = 2)
# Gfa <- loads %*% t(loads) + spec.var
# V.g <- V.g + (asreml.obj$design[, cols] %*% kronecker(Gfa, mat.I(62)) %*%
# t(as.matrix(asreml.obj$design[, cols])))
# V.g <- asreml.obj$sigma2 * (V.g + mat.I(n))
design <- asreml.obj$design
design <- replaceColnames(design)
for (term in ranterms[2:7])
{
cols <- grep(term, colnames(design), fixed = TRUE)
V.g <- V.g + asreml.obj$vparameters[term] * (design[, cols] %*% t(as.matrix(design[, cols])))
}
term <- ranterms[1]
cols <- grep(term, colnames(design), fixed = TRUE)[1:1364]
vp <- asreml.obj$vparameters[names(asreml.obj$vparameters)[grep(term,
names(asreml.obj$vparameters), fixed = TRUE)]]
spec.var <- diag(vp[grepl("!var", names(vp), fixed = TRUE)])
loads <- matrix(vp[grepl("!fa", names(vp), fixed = TRUE)], ncol = 2)
Gfa <- loads %*% t(loads) + spec.var
V.g <- V.g + (design[, cols] %*% kronecker(Gfa, mat.I(62)) %*%
t(as.matrix(design[, cols])))
V.g <- asreml.obj$sigma2 * (V.g + mat.I(n))
Vfa <- estimateV(asreml.obj)
testthat::expect_true(all(abs(Vfa - V.g) < 1e-06))
#Test rr
asreml.obj <-asreml(fixed = GY.tha ~ + at(expt, c(1:5)):rep + at(expt, c(1)):vrow +
at(expt, c(2,3,6,7)):colblocks +
at(expt, c(1:5,7)):vcol + Condition*expt,
random = ~ rr(exptCond, k = 2):Genotype +
at(expt, c(1)):dev(vrow) + at(expt, c(2)):spl(vcol) +
at(expt, c(3,5,7)):dev(vcol) + at(expt, c(7)):units,
data=comb.dat, maxit = 100, workspace = "1Gb")
summary(asreml.obj)$varcomp
ranterms <- names(asreml.obj$G.param)
n <- nrow(comb.dat)
V.g <- matrix(0, nrow = n, ncol = n)
design <- asreml.obj$design
design <- replaceColnames(design)
for (term in ranterms[2:7])
{
cols <- grep(term, colnames(design), fixed = TRUE)
V.g <- V.g + asreml.obj$vparameters[term] * (design[, cols] %*% t(as.matrix(design[, cols])))
}
term <- ranterms[1]
cols <- grep(term, colnames(design), fixed = TRUE)[1:1364]
vp <- asreml.obj$vparameters[names(asreml.obj$vparameters)[grep(term,
names(asreml.obj$vparameters), fixed = TRUE)]]
loads <- matrix(vp[grepl("!fa", names(vp), fixed = TRUE)], ncol = 2)
Gfa <- loads %*% t(loads)
V.g <- V.g + (design[, cols] %*% kronecker(Gfa, mat.I(62)) %*%
t(as.matrix(design[, cols])))
V.g <- asreml.obj$sigma2 * (V.g + mat.I(n))
Vrr <- estimateV(asreml.obj)
testthat::expect_true(all(abs(Vrr - V.g) < 1e-06))
asreml.options(design = FALSE)
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.