text: #devtools::test("asremlPlus")
  context("model_selection3")
  asr3.lib <- "D:\\Analyses\\R oldpkg" 
  
  cat("#### Test for chooseModel.data.frame with asreml3\n")
  test_that("choose.model.data.frame_asreml3", {
    skip_if_not_installed("asreml")
    skip_on_cran()
    library(dae)
    library(asreml, lib.loc = asr3.lib)
    library(asremlPlus)
    print(packageVersion("asreml"))
    #'## Load data
    data("Ladybird.dat")
    
    #'## ANOVA of logits
    Ladybird.aov <- aov(logitP ~ Host*Cadavers*Ladybird + Error(Run/Plant), 
                        data=Ladybird.dat)
    summary(Ladybird.aov)
    
    #'## Mixed model analysis of logits 
    m <- asreml(logitP ~ Host*Cadavers*Ladybird, 
                random = ~ Run,
                rcov = ~ Run:Plant,
                data = Ladybird.dat)
    testthat::expect_true(all(summary(m)$varcomp$bound == c("B", "P"))) #shows bound Run component
    
    #'### Unconstrain Reps to make the analysis equivalent to ANOVA
    m <- setvarianceterms(m$call, terms = "Run", bounds = "U")
    summary(m)$varcomp #shows negative Run component
    testthat::expect_true(m$gammas["Run!Run.var"] < 0)
    
    #'### Use chooseModel.data.frame
    wald.tab <- wald.asreml(m, denDF = "numeric")$Wald
    testthat::expect_equal(nrow(wald.tab), 8)
    
    #'### Choose marginality-compliant model from wald.tab, obtaining marginality using pstructure
    Ladybird.pstr <- pstructure(formula = ~ Host*Cadavers*Ladybird, 
                                data = Ladybird.dat)
    HCL.marg <- marginality(Ladybird.pstr)
    testthat::expect_equal(nrow(HCL.marg), 7)
    sigmod <- chooseModel(wald.tab, terms.marginality = HCL.marg)
    testthat::expect_true(all(unlist(lapply(sigmod$sig.terms, function(term) term)) == 
                                c("Cadavers:Ladybird", "Host")))
    testthat::expect_equal(nrow(sigmod$choose.summary), 5)
    testthat::expect_true(all(names(sigmod$choose.summary) == c("terms", "DF", "denDF", "p", "action")))
    
    #'### Rechoose marginality-compliant model from wald.tab, but omit DF and denDF
    sigmod <- chooseModel(wald.tab, omit.DF = TRUE, terms.marginality = HCL.marg)
    testthat::expect_equal(nrow(sigmod$choose.summary), 5)
    testthat::expect_true(all(names(sigmod$choose.summary) == c("terms", "p", "action")))
    sigmod <- chooseModel(wald.tab, denDF = NA, terms.marginality = HCL.marg)
    testthat::expect_equal(nrow(sigmod$choose.summary), 5)
    testthat::expect_true(all(names(sigmod$choose.summary) == c("terms", "DF", "denDF", "p", "action")))
  
    #'### Specify the denDF argument in various ways, all of which result in the overwriting of denDF
    sigmod <- chooseModel(wald.tab, denDF = wald.tab$denDF, terms.marginality = HCL.marg)
    testthat::expect_equal(nrow(sigmod$choose.summary), 5)
    testthat::expect_true(all(names(sigmod$choose.summary) == c("terms", "DF", "denDF", "p", "action")))
    
    den.df <- wald.tab$denDF
    sigmod <- chooseModel(wald.tab, denDF = den.df, terms.marginality = HCL.marg)
    testthat::expect_equal(nrow(sigmod$choose.summary), 5)
    testthat::expect_true(all(names(sigmod$choose.summary) == c("terms", "DF", "denDF", "p", "action")))
    
    sigmod <- chooseModel(wald.tab, denDF = 59, terms.marginality = HCL.marg)
    testthat::expect_equal(nrow(sigmod$choose.summary), 5)
    testthat::expect_true(all(names(sigmod$choose.summary) == c("terms", "DF", "denDF", "p", "action")))
    
    new.tab <- wald.tab
    names(new.tab)[match("denDF", names(new.tab))] <- "den.df"
    sigmod <- chooseModel(wald.tab, denDF = 59, terms.marginality = HCL.marg)
    testthat::expect_equal(nrow(sigmod$choose.summary), 5)
    testthat::expect_true(all(names(sigmod$choose.summary) == c("terms", "DF", "denDF", "p", "action")))
    
  })
  
  cat("#### Test for chooseModel.asrtests with asreml3\n")
  test_that("choose.model.asrtests_asreml3", {
    skip_if_not_installed("asreml")
    skip_on_cran()
    library(dae)
    library(asreml, lib.loc = asr3.lib)
    library(asremlPlus)
    print(packageVersion("asreml"))
    data(WaterRunoff.dat)
    current.asr <- asreml(log.Turbidity ~ Benches + (Sources * (Type + Species)) * Date, 
                          random = ~Benches:MainPlots:SubPlots:spl(xDay), 
                          data = WaterRunoff.dat, keep.order = TRUE)
    current.asrt <- as.asrtests(current.asr, NULL, NULL)
   
    #some tests for validWaldTab
    testthat::expect_error(test.wald <- as.asrtests(current.asr, 
                                                    wald.tab = WaterRunoff.dat))
    asrt.wald <- testranfix(current.asrt, term = "Sources:Species", ssType = "conditional")
    testthat::expect_equal(ncol(asrt.wald$wald.tab), 6)
    testthat::expect_true("F.con" %in% colnames(asrt.wald$wald.tab))
    
    terms.treat <- c("Sources", "Type", "Species", 
                     "Sources:Type", "Sources:Species")
    terms <- sapply(terms.treat, 
                    FUN=function(term){paste("Date:",term,sep="")}, 
                    simplify=TRUE)
    terms <- c("Date", terms)
    terms <- unname(terms)
    marginality <-  matrix(c(1,0,0,0,0,0, 1,1,0,0,0,0,  1,0,1,0,0,0, 
                             1,0,1,1,0,0, 1,1,1,0,1,0, 1,1,1,1,1,1), nrow=6)
    rownames(marginality) <- terms
    colnames(marginality) <- terms
    choose <- chooseModel(current.asrt, marginality, denDF="algebraic")
    current.asrt <- choose$asrtests.obj
    sig.terms <- choose$sig.terms
    testthat::expect_equal(length(sig.terms), 2)
    testthat::expect_equal(sig.terms[[1]], "Date:Species")
    testthat::expect_equal(sig.terms[[2]], "Date:Sources")
  })
  
  cat("#### Test for testing at terms with asreml3\n")
  test_that("at_testing_asreml3", {
    skip_if_not_installed("asreml")
    skip_on_cran()
    library(dae)
    library(asreml, lib.loc = asr3.lib)
    library(asremlPlus)
    print(packageVersion("asreml"))
    #'## Load data
    data("Exp355.Control.dat")
    
    dat <- with(dat, dat[order(Genotype, Salt, NP_AMF, InTreat), ])
    
    #Fit model with quotes around AMF_plus 
    # NB must have quotes for character levels, 
    # but cannot have in testranfix term because not in wald.tab rownames
    current.asr <- asreml(fixed = TSP ~ Lane + xPosn + AMF*Genotype*NP + 
                            at(AMF, "AMF_plus"):per.col + (Genotype*NP):at(AMF, "AMF_plus"):per.col,
                          random = ~ spl(xPosn) + Position ,
                          rcov = ~ Genotype:idh(NP_AMF):InTreat,
                          keep.order=TRUE, data = dat, 
                          maxiter=50, na.method.X = "include")
    
    current.asrt <- asrtests(current.asr, NULL, NULL)
    current.asrt <- rmboundary(current.asrt)
    testthat::expect_equal(nrow(current.asrt$wald.tab),14)
    testthat::expect_equal(nrow(summary(current.asrt$asreml.obj)$varcomp),7)
    
    t.asrt <- testranfix(current.asrt, term = "Genotype:NP:at(AMF, AMF_plus):per.col", 
                         drop.fix.ns = TRUE)
    t.asrt$wald.tab
    testthat::expect_equal(nrow(t.asrt$wald.tab), 13)
    #Change position of at term in testranfix
    t.asrt <- testranfix(current.asrt, term = "at(AMF, AMF_plus):per.col:Genotype:NP", 
                         drop.fix.ns = TRUE)
    t.asrt$wald.tab
    testthat::expect_equal(nrow(t.asrt$wald.tab), 13)
    
    #Fit model with level index of 2
    current.asr <- asreml(fixed = TSP ~ Lane + xPosn + AMF*Genotype*NP + 
                            at(AMF, "AMF_plus"):per.col + (Genotype*NP):at(AMF, 2):per.col,
                          random = ~ spl(xPosn) + Position ,
                          rcov = ~ Genotype:idh(NP_AMF):InTreat,
                          keep.order=TRUE, data = dat, 
                          maxiter=50, na.method.X = "include")
    
    current.asrt <- asrtests(current.asr, NULL, NULL)
    current.asrt <- rmboundary(current.asrt)
    testthat::expect_equal(nrow(current.asrt$wald.tab),14)
    testthat::expect_equal(nrow(summary(current.asrt$asreml.obj)$varcomp),7)
    
    t.asrt <- testranfix(current.asrt, term = "at(AMF, AMF_plus):per.col:Genotype:NP", 
                         drop.fix.ns = TRUE)
    t.asrt$wald.tab
    testthat::expect_equal(nrow(t.asrt$wald.tab), 13)
    
    #Test for a numeric level that is not the same as the levels index (1:no.levels)
    current.asr <- asreml(fixed = TSP ~ at(Lane, 4) + xPosn + AMF*Genotype*NP + 
                            at(AMF, "AMF_plus"):per.col + (Genotype*NP):at(AMF, c(2)):per.col,
                          random = ~ spl(xPosn) + Position ,
                          rcov = ~ Genotype:idh(NP_AMF):InTreat,
                          keep.order=TRUE, data = dat, 
                          maxiter=50,  na.method.X = "include")
    current.asrt <- asrtests(current.asr, NULL, NULL)
    current.asrt <- rmboundary(current.asrt)
    current.asrt$wald.tab
    t.asrt <- testranfix(current.asrt, term = "at(Lane, 8)", 
                         drop.fix.ns = TRUE)
    t.asrt$wald.tab
    testthat::expect_equal(nrow(t.asrt$wald.tab), 14)
  })
  
  
  cat("#### Test for testing MET at terms with asreml3\n")
  test_that("at_multilevel_asreml3", {
    skip_if_not_installed("asreml")
    skip_on_cran()
    library(dae)
    library(asreml, lib.loc = asr3.lib)
    library(asremlPlus)
    print(packageVersion("asreml"))
  
    #'## Load data
    data(MET)
  
    #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, maxiter = 100, workspace = 1e08)
    
    summary(asreml.obj)$varcomp
    current.asrt <- asrtests(asreml.obj, NULL, NULL)
    testthat::expect_equal(nrow(current.asrt$wald.tab), 24)
    
    #Single term in at expresion with the level and drop.fix.ns = TRUE -- works
    t.asrt <- testranfix(current.asrt, 
                         term = "at(expt, mtnue10):vrow", 
                         drop.fix.ns = TRUE,
                         dDF.na = "residual", update = FALSE)
    testthat::expect_equal(nrow(t.asrt$wald.tab), 23)
    
    
    #Multiple fixed terms in an at expresion generates an error
    testthat::expect_error(t.asrt <- testranfix(current.asrt, 
                                                term = "at(expt, c(1:5)):rep", 
                                                drop.fix.ns = TRUE,
                                                dDF.na = "residual", update = FALSE))
    
    #Multiple random terms in an at expression  generates an error
    testthat::expect_error(t.asrt <- testranfix(current.asrt, 
                                                term = "at(expt, c(3,5,7)):dev(vcol)", 
                                                drop.ran.ns = TRUE,
                                                dDF.na = "residual", update = FALSE))
  
    #Single random term in an at expression - thinks absent
    t.asrt <- testranfix(current.asrt, 
                         term = "at(expt, tarlee13):dev(vcol)", 
                         drop.ran.ns = TRUE,
                         dDF.na = "residual", update = FALSE)
    testthat::expect_equal(t.asrt$test.summary$action[1], "Absent")
    
    #Test multiple at term with changeTerms
    t.asrt <- changeTerms(current.asrt, dropFixed = "at(expt, c(1:5)):rep", update = FALSE)
    testthat::expect_equal(nrow(t.asrt$wald.tab), 19)
  })
  
  
  cat("#### Test for at terms in testswapran with asreml3\n")
  test_that("at_testswapran_asreml3", {
    skip_if_not_installed("asreml")
    skip_on_cran()
    library(dae)
    library(asreml, lib.loc = asr3.lib)
    library(asremlPlus)
    print(packageVersion("asreml"))
    #'## Load data
    data(longit.dat)
    
    current.asr <- do.call(asreml, 
                           args=list(fixed = Area ~ Block + Treatments + Treatments:xDAP,
                                     random = ~ Block:Cart + at(Treatments):spl(xDAP, k = 10) + Treatments:DAP + 
                                       Block:Cart:spl(xDAP) + Block:Cart:xDAP,
                                     rcov = ~ Block:Cart:ar1h(DAP),
                                     keep.order=TRUE, data = longit.dat, maxiter=100))
    
    #'## Function to deal with bound variances - set to 1e-04
    fixBoundResidualVariances <-function(current.asr)
    {
      repeat
      {
        current.call <- current.asr$call
        vpR <- grepl("R!DAP", names(current.asr$gammas), fixed = TRUE)
        vpR <- vpR & names(current.asr$gammas.con) == "Boundary"
        vpR <- current.asr$gammas[vpR]
        (terms <- names(vpR))
        if (length(terms) == 0 || length(sum(vpR == 4)) > 5) break
        current.asr <- setvarianceterms(current.call, terms = terms, 
                                        bounds = "F", initial.values = 0.0001,
                                        ignore.suffices = FALSE, stepsize = 0.001)
      }
      invisible(current.asr)
    }
    current.asr <- fixBoundResidualVariances(current.asr)
    testthat::expect_true(all(table(summary(current.asr)$varcomp$bound) ==  c(2,46,1)))
    
    #'## Load starting model into an asrtests object
    current.asrt <- as.asrtests(current.asr, NULL, NULL, label = "Selected variance model")
    testthat::expect_false(current.asrt$asreml.obj$converge)
    
    #'### Test for Treatments:DAP deviations terms
    current.asrt <- testranfix(current.asrt, term = "Treatments:DAP",
                               positive.zero = TRUE)
    testthat::expect_equal(current.asrt$test.summary$action[2], "Dropped")
    testthat::expect_true(all(table(summary(current.asrt$asreml.obj)$varcomp$bound) ==  c(2,45,1)))
    
    #'### Test for different curvatures in splines
    current.asrt <- testswapran(current.asrt, oldterms = "at(Treatments):spl(xDAP, k = 10)",
                                newterms = "at(AMF):Zn:spl(xDAP, k = 10)",
                                simpler = TRUE,
                                label = "Heterogeneous Treatment splines")
    testthat::expect_true(getTestPvalue(current.asrt, label = "Heterogeneous Treatment splines") < 0.05)
    testthat::expect_true(names(current.asrt$asreml.obj$gammas[1]) == 
                            "at(Treatments, -,0):spl(xDAP, k = 10)")
    testthat::expect_true(names(current.asrt$asreml.obj$gammas[5]) == 
                            "at(Treatments, +,0):spl(xDAP, k = 10)")
    testthat::expect_true(all(abs(current.asrt$asreml.obj$vparameters[1:8] - 
                                    c(233.152932, 502.667930, 74.955973, 
                                      2.540186, 42.003197, 61.206138, 
                                      32.367734, 36.902978)) < 1e-02))
  })
  
  cat("#### Test for spline testing with asreml3\n")
  test_that("spl.asrtests_asreml3", {
    skip_if_not_installed("asreml")
    skip_on_cran()
    library(dae)
    library(asreml, lib.loc = asr3.lib)
    library(asremlPlus)
    print(packageVersion("asreml"))
    data(WaterRunoff.dat)
    current.asr <- do.call("asreml", 
                           args = list(fixed = log.Turbidity ~ Benches + 
                                         (Sources * (Type + Species)) * Date, 
                                       random = ~Benches:MainPlots:SubPlots:spl(xDay, k = 6), 
                                       keep.order = TRUE, data = WaterRunoff.dat))
    current.asrt <- as.asrtests(current.asr, NULL, NULL)
    
    #Test random splines
    current.asrt <- testranfix(current.asrt, term = "Benches:MainPlots:SubPlots:spl(xDay, k = 6)")
    current.asrt$test.summary
    
    testthat::expect_equal(nrow(current.asrt$test.summary), 1)
    testthat::expect_true(abs(current.asrt$test.summary$p - 0.08013755) < 1e-08)
    
    data(Wheat.dat)
    #'## Add cubic trend to Row so that spline is not bound
    Wheat.dat <- within(Wheat.dat, 
                        {
                          vRow <- as.numeric(Row)
                          vRow <- vRow - mean(unique(vRow))
                          yield <- yield + 10*vRow + 5 * (vRow^2) + 5 * (vRow^3)
                        })
    
    #'## Fit model using asreml
    asreml.obj <- asreml(fixed = yield ~ Rep + vRow + Variety, 
                         random = ~spl(vRow, k=6) + units, 
                         rcov = ~ar1(Row):ar1(Column), 
                         data = Wheat.dat, trace = FALSE)
    testthat::expect_true(summary(asreml.obj)$varcomp$constraint[1] == "Boundary")
    asreml.obj <- asreml(fixed = yield ~ Rep + vRow + Variety, 
                         random = ~spl(vRow) + units, 
                         rcov = ~ar1(Row):ar1(Column), 
                         data = Wheat.dat, trace = FALSE)
    testthat::expect_true(summary(asreml.obj)$varcomp$constraint[1] == "Boundary")
    
  })
  
  cat("#### Test for reparamSigDevn.asrtests with asreml3\n")
  test_that("reparamSigDevn.asrtests_asreml3", {
    skip_if_not_installed("asreml")
    skip_on_cran()
    library(dae)
    library(asreml, lib.loc = asr3.lib)
    library(asremlPlus)
    data(WaterRunoff.dat)
    current.asr <- asreml(fixed = log.Turbidity ~ Benches + Sources + Type + Species + 
                            Sources:Type + Sources:Species + Sources:Species:xDay + 
                            Sources:Species:Date, 
                          data = WaterRunoff.dat, keep.order = TRUE)
    current.asrt <- as.asrtests(current.asr, NULL, NULL)
    
    #Examine terms that describe just the interactions of Date and the treatment factors
    terms.treat <- c("Sources", "Type", "Species", "Sources:Type", "Sources:Species")
    date.terms <- sapply(terms.treat, 
                         FUN=function(term){paste("Date:",term,sep="")}, 
                         simplify=TRUE)
    date.terms <- c("Date", date.terms)
    date.terms <- unname(date.terms)
    treat.marginality <-  matrix(c(1,0,0,0,0,0, 1,1,0,0,0,0,  1,0,1,0,0,0, 
                                   1,0,1,1,0,0, 1,1,1,0,1,0, 1,1,1,1,1,1), nrow=6)
    rownames(treat.marginality) <- date.terms
    colnames(treat.marginality) <- date.terms
    choose <- chooseModel(current.asrt, treat.marginality, denDF="algebraic")
    current.asrt <- choose$asrtests.obj
    current.asr <- current.asrt$asreml.obj
    sig.date.terms <- choose$sig.terms
    
    #Remove all Date terms left in the fixed model
    terms <- "(Date/(Sources * (Type + Species)))"
    current.asrt <- changeTerms(current.asrt, dropFixed = terms)
    #if there are significant date terms, reparameterize to xDays + spl(xDays) + Date
    if (length(sig.date.terms) != 0)
    { #add lin + spl + devn for each to fixed and random models
      trend.date.terms <- sapply(sig.date.terms, 
                                 FUN=function(term){sub("Date","xDay",term)}, 
                                 simplify=TRUE)
      trend.date.terms <- paste(trend.date.terms,  collapse=" + ")
      current.asrt <- changeTerms(current.asrt, addFixed=trend.date.terms)
      trend.date.terms <- sapply(sig.date.terms, 
                                 FUN=function(term){sub("Date","spl(xDay)",term)}, 
                                 simplify=TRUE)
      trend.date.terms <- c(trend.date.terms, sig.date.terms)
      trend.date.terms <- paste(trend.date.terms,  collapse=" + ")
      current.asrt <- changeTerms(current.asrt, addRandom = trend.date.terms)
      current.asrt <- rmboundary.asrtests(current.asrt)
    }
    
    #Now test terms for sig date terms
    spl.terms <- sapply(terms.treat, 
                        FUN=function(term){paste("spl(xDay):",term,sep="")}, 
                        simplify=TRUE)
    spl.terms <- c("spl(xDay)",spl.terms)
    lin.terms <- sapply(terms.treat, 
                        FUN=function(term){paste(term,":xDay",sep="")}, 
                        simplify=TRUE)
    lin.terms <- c("xDay",lin.terms)
    systematic.terms <- c(terms.treat, lin.terms, spl.terms, date.terms)
    systematic.terms <- unname(systematic.terms)
    treat.marginality <-  matrix(c(1,0,0,0,0,0, 1,1,0,0,0,0,  1,0,1,0,0,0, 
                                   1,0,1,1,0,0, 1,1,1,1,1,0, 1,1,1,1,1,1), nrow=6)
    systematic.marginality <- kronecker(matrix(c(1,0,0,0, 1,1,0,0, 
                                                 1,1,1,0, 1,1,1,1), nrow=4), 
                                        treat.marginality)
    systematic.marginality <- systematic.marginality[-1, -1]
    rownames(systematic.marginality) <- systematic.terms
    colnames(systematic.marginality) <- systematic.terms
    choose <- chooseModel(current.asrt, systematic.marginality, 
                          denDF="algebraic", pos=TRUE)
    current.asrt <- choose$asrtests.obj
    
    #Check if any deviations are significant and, for those that are, go back to 
    #fixed dates
    current.asrt <- reparamSigDevn(current.asrt, choose$sig.terms, 
                                   trend.num = "xDay", devn.fac = "Date", 
                                   denDF = "algebraic")
    k <- match("Sources:Species:Date",rownames(current.asrt$wald.tab))
    testthat::expect_equal(nrow(current.asrt$wald.tab), 9)
    testthat::expect_equal(nrow(current.asrt$test.summary), 6)
    testthat::expect_true(!is.na(k))
  })
  
  cat("#### Test for changeModelOnIC with wheat94 using asreml3\n")
  test_that("changeModelOnIC_wheat94_asreml3", {
    skip_if_not_installed("asreml")
    skip_on_cran()
    library(dae)
    library(asreml)
    library(asremlPlus)
    ## use asremlPlus to analyse the 1994 wheat example from Gilmour et al. (1995)
    data(wheat94.dat)
    
    
    #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)
  
    #Use REML likelihood and BIC
    current.asrt <- as.asrtests(fm.max, NULL, label = "Maximal model", 
                                IClikelihood = "REML")
    current.asrt <- iterate(current.asrt)
    current.asrt <- rmboundary(current.asrt)
    testthat::expect_equal(nrow(current.asrt$test.summary), 2)
    
    #Drop random Row and Col terms
    current.asrt <- changeModelOnIC(current.asrt, dropRandom = "Row + Col", 
                                    label = "Drop Row + Col", 
                                    which.IC = "BIC", IClikelihood = "REML")
    testthat::expect_equal(current.asrt$test.summary$denDF[3], -2)
    testthat::expect_equal(current.asrt$test.summary$action[current.asrt$test.summary$terms == 
                                                              "Drop Row + Col"], "Swapped")
    
    #Drop random spl(Col) term
    current.asrt <- changeModelOnIC(current.asrt, dropRandom = "spl(Col)", 
                                    label = "Drop spl(Col)", 
                                    which.IC = "BIC", IClikelihood = "REML")
    testthat::expect_equal(current.asrt$test.summary$denDF[4], -1)
    testthat::expect_equal(current.asrt$test.summary$action[current.asrt$test.summary$terms == 
                                                              "Drop spl(Col)"], "Swapped")
    testthat::expect_true((abs(current.asrt$test.summary$BIC[4]) - 1.764507) < 1e-02)
    
    #Drop random units term
    current.asrt <- changeModelOnIC(current.asrt, dropRandom = "units", 
                                    label = "Drop units", 
                                    which.IC = "BIC", IClikelihood = "REML")
    testthat::expect_equal(current.asrt$test.summary$denDF[5], -1)
    testthat::expect_equal(current.asrt$test.summary$action[current.asrt$test.summary$terms == 
                                                              "Drop units"], "Unswapped")
    testthat::expect_true((abs(current.asrt$test.summary$AIC[5]) -60.520592) < 1e-02)
    
    mod <- printFormulae(current.asrt$asreml.obj)
    testthat::expect_equal(length(mod), 3)
    
  })
  
  
  cat("#### Test for changeModelOnIC example using asreml3\n")
  test_that("changeModelOnIC_Example_asreml3", {
    skip_if_not_installed("asreml")
    skip_on_cran()
    library(dae)
    library(asreml)
    library(asremlPlus)
    ## use asremlPlus to analyse the wheat (barley) example from section 8.6 of the asreml manual (Butler et al. 2010)
    data(Wheat.dat)
    
    #'## Fit maximal model
    current.asr <- asreml(yield ~ Rep + WithinColPairs + Variety, 
                          random = ~ Row + Column + units,
                          rcov = ~ ar1(Row):ar1(Column), 
                          data=Wheat.dat)
    current.asr <- update(current.asr)
    current.asrt <- as.asrtests(current.asr, NULL, NULL, 
                                label = "Maximal model", IClikelihood = "REML")
    current.asrt <- rmboundary(current.asrt)
    testthat::expect_true(current.asrt$asreml.obj$converge)
    testthat::expect_true(current.asrt$test.summary$action[1] == "Starting model")
    testthat::expect_equal(current.asrt$test.summary$DF[1], 0)
    testthat::expect_equal(current.asrt$test.summary$denDF[1], 5)
    testthat::expect_equal(nrow(summary(current.asrt$asreml.obj)$varcomp), 5)
    
    # Drop both Row and Column terms
    current.asrt <- changeModelOnIC(current.asrt, 
                                    dropRandom = "Row + Column", label = "Drop Row + Column",
                                    checkboundaryonly = TRUE,
                                    which.IC = "AIC", IClikelihood = "REML")
    testthat::expect_true(current.asrt$asreml.obj$converge)
    testthat::expect_equal(current.asrt$test.summary$denDF[3], -1)
    
    # Replace residual with model without Row autocorrelation
    current.asrt <- changeModelOnIC(current.asrt, 
                                    newResidual = "Row:ar1(Column)", 
                                    label="Row autocorrelation",
                                    IClikelihood = "REML")
    testthat::expect_true(current.asrt$asreml.obj$converge)
    testthat::expect_equal(current.asrt$test.summary$denDF[4], -2)
    testthat::expect_true((abs(current.asrt$test.summary$AIC[4]) - 24.158777) < 1e-03)
    
    mod <- printFormulae(current.asrt$asreml.obj)
    testthat::expect_equal(length(mod), 3)
    testthat::expect_true(grepl("units", mod[2], fixed = TRUE))
    
  })
  
  