Nothing
#devtools::test("asremlPlus")
context("model_selection")
asr3.lib <- "D:\\Analyses\\R oldpkg"
cat("#### Test for REMLRT with asreml3\n")
test_that("REMLRT_asreml3", {
skip_if_not_installed("asreml")
skip_on_cran()
library(dae)
library(asremlPlus)
loadASRemlVersion(3, lib.loc = asr3.lib)
## use asremlPlus to analyse the wheat (barley) example from section 8.6 of the asreml manual (Butler et al. 2010)
data(Wheat.dat)
# Fit initial model
m1.asr <- asreml(yield ~ Rep + WithinColPairs + Variety,
random = ~ Row + Column + units,
rcov = ~ ar1(Row):ar1(Column),
data=Wheat.dat)
summary(m1.asr)$varcomp
info <- infoCriteria(m1.asr)
testthat::expect_equal(info$varDF, 5)
testthat::expect_lt(abs(info$AIC - 1346.766), 1e-03)
#Fit model without the units term
m2.asr <- asreml(yield ~ Rep + WithinColPairs + Variety,
random = ~ Row + Column,
rcov = ~ ar1(Row):ar1(Column),
data=Wheat.dat)
summary(m2.asr)$varcomp
info <- infoCriteria(m2.asr)
testthat::expect_equal(info$varDF, 4)
testthat::expect_lt(abs(info$AIC - 1352.941), 1e-03)
test <- REMLRT(m2.asr, m1.asr)
testthat::expect_lt(abs(test$p - 0.004232946), 1e-03)
testthat::expect_equal(test$DF, 1)
test <- REMLRT(m2.asr, m1.asr, DF = 1)
testthat::expect_lt(abs(test$p - 0.004232946), 1e-03)
testthat::expect_equal(test$DF, 1)
m3.asr <- asreml(yield ~ Rep + WithinColPairs + Variety,
random = ~ Row + Column,
rcov = ~ Row:Column,
data=Wheat.dat)
summary(m3.asr)$varcomp
test3 <- REMLRT(m3.asr, m1.asr)
testthat::expect_lt(abs(test3$p - 2.596812e-13), 1e-03)
testthat::expect_equal(test3$DF, 2)
test3 <- REMLRT(m3.asr, m1.asr, DF = 3)
testthat::expect_lt(abs(test3$p - 1.603828e-12), 1e-03)
testthat::expect_equal(test3$DF, 3)
info <- infoCriteria(m3.asr, IClikelihood = "REML")
testthat::expect_equal(info$fixedDF, 0)
testthat::expect_equal(info$varDF, 3)
testthat::expect_lt(abs(info$AIC - 1400.719), 5e-03)
testthat::expect_lt(abs(info$BIC - 1409.056), 5e-03)
testthat::expect_lt(abs(info$loglik - m3.asr$loglik), 1e-06)
})
cat("#### Test for wheat76 example with asreml3\n")
test_that("Wheat_asreml3", {
skip_if_not_installed("asreml")
skip_on_cran()
library(asremlPlus)
loadASRemlVersion(3, lib.loc = asr3.lib)
## Fit several models to the wheat data and caclulate their ICs
data(Wheat.dat)
# Fit initial model
m.max <- asreml(yield ~ Rep + WithinColPairs + Variety,
random = ~ Row + Column + units,
rcov = ~ ar1(Row):ar1(Column),
data=Wheat.dat)
#Drop term for within Column pairs
m1 <- asreml(yield ~ Rep + Variety,
random = ~ Row + Column + units,
rcov = ~ ar1(Row):ar1(Column),
data=Wheat.dat)
#Drop nugget term
m2 <- asreml(yield ~ Rep + WithinColPairs + Variety,
random = ~ Row + Column,
rcov = ~ ar1(Row):ar1(Column),
data=Wheat.dat)
#Drop Row autocorrelation
m3 <- asreml(yield ~ Rep + WithinColPairs + Variety,
random = ~ Row + Column + units,
rcov = ~ Row:ar1(Column),
data=Wheat.dat)
#Drop Col autocorrelation
m4 <- asreml(yield ~ Rep + WithinColPairs + Variety,
random = ~ Row + Column + units,
rcov = ~ ar1(Row):Column,
data=Wheat.dat)
mods.asr <- list(m.max, m1, m2, m3, m4)
ic <- infoCriteria(mods.asr, IClikelihood = "REML")
testthat::expect_equal(nrow(ic), 5)
testthat::expect_true(all(ic$fixedDF == 0))
testthat::expect_true(all(ic$varDF == c(5, 5, 4, 4, 5)))
testthat::expect_true(all(abs(ic$AIC - c(1346.766,1353.762,1352.941,1363.901,1393.475)) < 1e-01))
testthat::expect_true(abs(ic$BIC[1] - 1360.662) < 1)
testthat::expect_true(abs(ic$loglik[1] - (-668.3832)) < 1e-01)
})
cat("#### Test for IC with wheat94 using asreml3\n")
test_that("IC_wheat94_asreml3", {
skip_if_not_installed("asreml")
skip_on_cran()
library(dae)
library(asremlPlus)
loadASRemlVersion(3, lib.loc = asr3.lib)
## use asremlPlus to analyse the wheat (barley) example from section 8.6 of the asreml manual (Butler et al. 2010)
data(wheat94.dat)
### Start with a simple model
fm0 <- asreml(yield ~ 1,
random = ~ Variety + Block,
data = wheat94.dat)
fm0 <- update(fm0)
current.asrt <- as.asrtests(fm0, NULL, NULL,
label = "Simple model", IClikelihood = "REML")
testthat::expect_equal(nrow(current.asrt$wald.tab), 1)
#Add autocorrelation
current.asrt <- changeTerms(current.asrt, newResidual = "ar1(Col):ar1(Row)",
label = "Add autocorrelation", IClikelihood = "REML")
testthat::expect_true(abs(diff(current.asrt$test.summary$AIC)) - 312.1018 < 1e-03)
#Add units term
current.asrt <- changeTerms(current.asrt, addRandom = "units",
label = "Add units", IClikelihood = "REML")
vpar3 <- current.asrt$asreml.obj$gammas[1:3]
current.asrt <- iterate(current.asrt)
testthat::expect_true(all(abs(current.asrt$asreml.obj$gammas[1:3] - vpar3) < 1e-04))
testthat::expect_true(abs(current.asrt$asreml.obj$loglik - -1563.459) < 1e-03)
testthat::expect_equal(nrow(current.asrt$wald.tab), 1)
#Add random Row and Col terms
current.asrt <- changeTerms(current.asrt, addRandom = "Row + Col",
label = "Add Row + Col", IClikelihood = "REML")
current.asrt <- iterate(current.asrt)
#check that denDf for current model is the same as the number of variance parameters
testthat::expect_true(
nrow(summary(current.asrt$asreml.obj)$varcomp) ==
current.asrt$test.summary$denDF[current.asrt$test.summary$terms == "Add Row + Col"])
#Add fixed lin(Row) and lin(Col) terms
current.asrt <- changeTerms(current.asrt, addFixed = "lin(Row) + lin(Col)",
label = "Add lin(Row) + lin(Col)", IClikelihood = "REML")
#three fixed parameters?
testthat::expect_equal(nrow(current.asrt$wald.tab), 3)
#Add random spl(Col) term
current.asrt <- changeTerms(current.asrt, addRandom = "spl(Col)",
label = "Add spl(Col)",
IClikelihood = "REML")
testthat::expect_equal(
nrow(summary(current.asrt$asreml.obj)$varcomp),
current.asrt$test.summary$denDF[current.asrt$test.summary$terms == "Add spl(Col)"][1])
#Restart with fixed Rowcode and Colcode covariates, units and autocorrelation
fm6 <- asreml(yield ~ Rowcode + Colcode,
random = ~ Variety + Block + units,
rcov = ~ ar1(Col):ar1(Row),
data = wheat94.dat)
fm6 <- update(fm6)
current.asrt <- as.asrtests(fm6, wald.tab = NULL,
test.summary = current.asrt$test.summary,
label = "Basic + Row/Col covariates", IClikelihood = "REML")
testthat::expect_true(tail(current.asrt$test.summary$action,1) == "Starting model")
#Add random Row and Col terms
current.asrt <- changeTerms(current.asrt, addRandom = "Row + Col",
label = "Add Row + Col", IClikelihood = "REML")
current.asrt <- iterate(current.asrt)
#Add fixed lin(Row) and lin(Col) terms
current.asrt <- changeTerms(current.asrt, addFixed = "lin(Row) + lin(Col)",
label = "Add lin(Row) + lin(Col)", IClikelihood = "REML")
#Add random spl(Col) term
current.asrt <- changeTerms(current.asrt, addRandom = "spl(Col)",
label = "Add spl(Col)", IClikelihood = "REML")
current.asrt <- iterate(current.asrt)
testthat::expect_equal(nrow(current.asrt$test.summary), 14)
print(current.asrt$test.summary, omit.columns = "p")
#Start with Maximal model
fm.max <- asreml(yield ~ lin(Row) + lin(Col) + Rowcode + Colcode,
random = ~ Variety + Block + Row + spl(Col) + Col + units,
rcov = ~ ar1(Col):ar1(Row),
data = wheat94.dat)
current.asrt <- as.asrtests(fm.max, NULL, NULL,
label = "Maximal model", IClikelihood = "REML")
current.asrt <- iterate(current.asrt)
testthat::expect_true(tail(current.asrt$test.summary$action,1) == "Starting model")
testthat::expect_equal(current.asrt$test.summary$DF, 0)
testthat::expect_equal(current.asrt$test.summary$denDF, 8)
testthat::expect_equal(nrow(summary(current.asrt$asreml.obj)$varcomp), 9) #includes bound Block
current.asrt <- rmboundary(current.asrt)
testthat::expect_equal(nrow(summary(current.asrt$asreml.obj)$varcomp),
current.asrt$test.summary$denDF[1])
#Drop random Row and Col terms
current.asrt <- changeTerms(current.asrt, dropRandom = "Row + Col",
label = "Drop Row + Col", IClikelihood = "REML")
testthat::expect_equal(nrow(summary(current.asrt$asreml.obj)$varcomp),
current.asrt$test.summary$denDF[3])
#Drop random spl(Col) term
current.asrt <- changeTerms(current.asrt, dropRandom = "spl(Col)",
label = "Drop spl(Col)", IClikelihood = "REML")
testthat::expect_equal(nrow(summary(current.asrt$asreml.obj)$varcomp),
current.asrt$test.summary$denDF[4])
testthat::expect_true((abs(diff(current.asrt$test.summary$BIC[3:4])) - 4.062308) < 1e-05)
#Use hypothesis testing with the maximal model
current.asrt <- as.asrtests(fm.max, NULL, test.summary = current.asrt$test.summary,
label = "Maximal model", IClikelihood = "REML")
current.asrt <- iterate(current.asrt)
current.asrt <- rmboundary(current.asrt)
testthat::expect_equal(nrow(current.asrt$test.summary), 6)
#Test random Row term
current.asrt <- testranfix(current.asrt, term = "Row", alpha = 0.20)
#Test random Col term
current.asrt <- testranfix(current.asrt, term = "Col", alpha = 0.20)
current.asrt <- iterate(current.asrt)
#test random spl(Col) term
if (getTestPvalue(current.asrt, label = "Col") > 0.05)
current.asrt <- testranfix(current.asrt, term = "spl(Col)", alpha = 0.20)
#Test units term
current.asrt <- testranfix(current.asrt, term = "units", alpha = 0.20)
testthat::expect_equal(nrow(current.asrt$test.summary), 10)
})
cat("#### Test for getFormulae with wheat94 using asreml\n")
test_that("Formulae_wheat94_asreml3", {
skip_if_not_installed("asreml")
skip_on_cran()
library(dae)
library(asremlPlus)
loadASRemlVersion(3, lib.loc = asr3.lib)
## use asremlPlus to analyse the wheat (barley) example from section 8.6 of the asreml manual (Butler et al. 2010)
data(wheat94.dat)
fm.max <- asreml(yield ~ lin(Row) + lin(Col) + Rowcode + Colcode,
random = ~ Variety + Block + Row + spl(Col) + Col + units,
rcov = ~ ar1(Col):ar1(Row),
data = wheat94.dat)
mod <- getFormulae(fm.max, which = "all", envir = sys.frame(sys.nframe()))
testthat::expect_true(all(unlist(lapply(mod, function(form) is.null(form) |
inherits(form, what = "formula")))))
testthat::expect_true(all(names(mod) == c("fixed", "random", "rcov", "sparse")))
#Print fitted model
testthat::expect_equal(length(mod), 4)
testthat::expect_true(is.null(mod$sparse))
p <- printFormulae(fm.max, which = "all")
testthat::expect_true(all(nchar(p) > 11))
testthat::expect_equal(length(p), 4)
p <- printFormulae(fm.max, expanded = TRUE)
testthat::expect_equal(length(p), 3)
p <- printFormulae(fm.max, which = c("fixed", "random"))
testthat::expect_equal(length(p), 2)
p <- printFormulae(fm.max, which = "fixed")
testthat::expect_equal(length(p), 1)
#Test when have formulae are in a character or list
fix.mod <- mod$fixed
fm.max <- asreml(fixed = fix.mod,
random = mod$random,
rcov = mod$rcov,
data = wheat94.dat)
mod <- getFormulae(fm.max, which = "all", envir = sys.frame(sys.nframe()))
testthat::expect_true(all(unlist(lapply(mod, function(form) is.null(form) |
inherits(form, what = "formula")))))
testthat::expect_equal(length(mod), 4)
testthat::expect_true(is.null(mod$sparse))
p <- printFormulae(fm.max, which = "all")
testthat::expect_true(all(nchar(p) > 11))
})
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.