#devtools::test("asremlPlus")
context("model_selection")

cat("#### Test for wheat76 example with asreml42\n")
test_that("Wheat_asreml42", {
  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 initial model
  current.asr <- asreml(yield ~ Rep + WithinColPairs + Variety, 
                        random = ~ Row + Column + units,
                        residual = ~ ar1(Row):ar1(Column), 
                        maxit = 30, data=Wheat.dat)
  summary(current.asr)
  info <- infoCriteria(current.asr)
  testthat::expect_equal(info$varDF, 5)
  testthat::expect_lt(abs(info$AIC - 1346.768), 0.10)
  
  # Load current fit into an asrtests object
  current.asrt <- as.asrtests(current.asr, NULL, NULL, 
                              label = "Maximal model", IClikelihood = "full")
  testthat::expect_lt(abs(current.asrt$test.summary$AIC - 1653.098), 0.10)
  
  
  # Check for and remove any boundary terms
  current.asrt <- rmboundary(current.asrt, IClikelihood = "full")
  
  #Check term for within Column pairs
  current.asrt <- testranfix(current.asrt, term = "WithinColPairs", 
                             drop.fix.ns=TRUE, IClikelihood = "full")
  
  # Test nugget term
  current.asrt <- testranfix(current.asrt, "units", positive=TRUE, IClikelihood = "full")
  
  # Test Row autocorrelation
  current.asrt <- testresidual(current.asrt, "~ Row:ar1(Column)", 
                               label="Row autocorrelation", simpler=TRUE, 
                               IClikelihood = "full")
  
  # Test Col autocorrelation (depends on whether Row autocorrelation retained)
  p <- getTestPvalue(current.asrt, label = "Row autocorrelation")
  testthat::expect_true((abs(p - 2.314881e-06) < 1e-05))
  { 
    if (p <= 0.05)
    current.asrt <- testresidual(current.asrt, "~ ar1(Row):Column", 
                                 label="Col autocorrelation", simpler=TRUE,
                                 IClikelihood = "full", update=FALSE)
    else
      current.asrt <- testresidual(current.asrt, "~ Row:Column", 
                                   label="Col autocorrelation", simpler=TRUE,
                                   IClikelihood = "full", update=FALSE)
  }
  print(current.asrt)
  testthat::expect_equal(length(current.asrt), 3)
  testthat::expect_equal(nrow(current.asrt$wald.tab), 3)
  testthat::expect_equal(nrow(current.asrt$test.summary), 6)
  testthat::expect_lt(abs(current.asrt$test.summary$AIC[6] - 1651.329), 0.10)
  testthat::expect_lt(abs(current.asrt$test.summary$BIC[6] - 1756.701), 0.10)
  info <- infoCriteria(current.asrt$asreml.obj)
  testthat::expect_equal(info$varDF, 5)
  testthat::expect_lt(abs(info$AIC - 1353.762), 1e-03)
  
  # Get current fitted asreml object
  current.asr <- current.asrt$asreml.obj
  current.asr <- update(current.asr, aom=TRUE)
  
  
  # Do residuals-versus-fitted values plot
  plot(fitted(current.asr),residuals(current.asr))
  
  #Produce variogram and variogram faces plot (Stefanaova et al, 2009)
  if (requireNamespace("lattice"))
  {
    plot.varioGram(varioGram.asreml(current.asr))
  }
  V <- estimateV(current.asr)
  faces <- variofaces(current.asr, V=V, maxit=50, units="addtores",
                      ncores = parallel::detectCores())
  testthat::expect_equal(nrow(faces$face1), 10)
  testthat::expect_equal(nrow(faces$face2), 15)
  
  #Get Variety predictions and all pairwise prediction differences and p-values
  Var.diffs <- predictPlus(classify = "Variety", 
                           asreml.obj=current.asr, 
                           error.intervals="halfLeast",
                           wald.tab=current.asrt$wald.tab,
                           tables = "predictions", 
                           sortFactor = "Variety")
  testthat::expect_equal(nrow(Var.diffs$predictions), 25)
  testthat::expect_equal(ncol(Var.diffs$predictions), 6)
  testthat::expect_equal(nrow(Var.diffs$differences), 25)
  testthat::expect_equal(ncol(Var.diffs$differences), 25)
  testthat::expect_equal(nrow(Var.diffs$p.differences), 25)
  testthat::expect_equal(ncol(Var.diffs$p.differences), 25)
  testthat::expect_equal(length(Var.diffs$LSD), 8)
  testthat::expect_true("lower.halfLeastSignificant.limit" %in% names(Var.diffs$predictions))
  testthat::expect_equal(Var.diffs$backtransforms, NULL)
  testthat::expect_equal(as.character(Var.diffs$predictions$Variety[[1]]),"10")
  testthat::expect_silent(plotPvalues(Var.diffs))
  testthat::expect_silent(plotPvalues(Var.diffs, show.sig = TRUE, alpha = 0.05))
  
  #Test for single-value LSDs
  diffs <- predictPlus(classify = "Variety", 
                           asreml.obj=current.asr, 
                           error.intervals="halfLeast",
                           LSDtype = "fact", LSDby = "Variety",
                           wald.tab=current.asrt$wald.tab,
                           tables = "predictions", 
                           sortFactor = "Variety")
  testthat::expect_equal(nrow(diffs$LSD), 25)
  testthat::expect_true("lower.halfLeastSignificant.limit" %in% names(diffs$predictions))
  #Are the LSDs sorted?
  testthat::expect_true(all((diffs$predictions$upper.halfLeastSignificant.limit - 
                               diffs$predictions$lower.halfLeastSignificant.limit - 
                               diffs$LSD$meanLSD) < 1e-05))
  diffs$predictions$upper.halfLeastSignificant.limit - diffs$predictions$lower.halfLeastSignificant.limit
})


cat("#### Test for wheat76 example using AIC with asreml42\n")
test_that("Wheat_asreml42", {
  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)
  ## The AIC is used in asremlPlus-package.RD
  data(Wheat.dat)
  
  # Fit initial model
  current.asr <- asreml(yield ~ Rep + WithinColPairs + Variety, 
                        random = ~ Row + Column + units,
                        residual = ~ ar1(Row):ar1(Column), 
                        maxit = 30, data=Wheat.dat)
  summary(current.asr)
  info <- infoCriteria(current.asr)
  testthat::expect_equal(info$varDF, 5)
  testthat::expect_lt(abs(info$AIC - 1346.768), 0.10)
  
  # Load current fit into an asrtests object
  current.asrt <- as.asrtests(current.asr, NULL, NULL, 
                              label = "Maximal model", IClikelihood = "full")
  testthat::expect_lt(abs(current.asrt$test.summary$AIC - 1653.098), 0.10)
  
  
  # Check for and remove any boundary terms
  current.asrt <- rmboundary(current.asrt, IClikelihood = "full")
  
  #Check term for within Column pairs
  current.asrt <- testranfix(current.asrt, term = "WithinColPairs", 
                             drop.fix.ns=TRUE, IClikelihood = "full")
  
  # Test nugget term
  current.asrt <- changeModelOnIC(current.asrt, addRandom = "units", 
                                  label = "units", IClikelihood = "full")
  
  # Test Row autocorrelation
  current.asrt <- changeModelOnIC(current.asrt, newResidual = "Row:ar1(Column)", 
                                  label="Row autocorrelation", 
                                  IClikelihood = "full")
  
  # Test Col autocorrelation (depends on whether Row autocorrelation retained)
  result <- getTestEntry(current.asrt, label = "Row autocorrelation")$action
  testthat::expect_true(grepl("Unswapped", result))
  { 
    if (grepl("Unswapped", result) || grepl("Unchanged", result))
      current.asrt <- changeModelOnIC(current.asrt, newResidual = "ar1(Row):Column", 
                                      label="Col autocorrelation", 
                                      IClikelihood = "full", update=FALSE)
    else
      current.asrt <- changeModelOnIC(current.asrt, newResidual = "Row:Column", 
                                      label="Col autocorrelation", 
                                      IClikelihood = "full", update=FALSE)
  }
  print(current.asrt)
  testthat::expect_equal(length(current.asrt), 3)
  testthat::expect_equal(nrow(current.asrt$wald.tab), 3)
  testthat::expect_equal(nrow(current.asrt$test.summary), 8)
  testthat::expect_lt(abs(current.asrt$test.summary$AIC[8] - 54.18178706), 0.10)
  testthat::expect_lt(abs(current.asrt$test.summary$BIC[8] - 51.17115177), 0.10)
  info <- infoCriteria(current.asrt$asreml.obj)
  testthat::expect_equal(info$varDF, 5)
  testthat::expect_lt(abs(info$AIC - 1353.762), 1e-03)
})





cat("#### Test for wheat76 addtoSummary with asreml42\n")
test_that("Wheat_addtoSummary_asreml42", {
  skip_if_not_installed("asreml")
  skip_on_cran()
  library(dae)
  library(asreml)
  library(asremlPlus)
  data(Wheat.dat)
  
  ## Fit an autocorrelation model
  ar1.asr <- do.call(asreml, 
                     list(yield ~ Rep + WithinColPairs + Variety, 
                          random = ~ Row + Column + units,
                          residual = ~ ar1(Row):ar1(Column), 
                          data=Wheat.dat))
  ar1.asrt <- as.asrtests(ar1.asr, NULL, NULL, 
                          label = "Autocorrelation model")
  ar1.asrt <- rmboundary.asrtests(ar1.asrt)
  testthat::expect_equal(nrow(ar1.asrt$test.summary),2)
  
  ## Fit a tensor spline
  Wheat.dat <- within(Wheat.dat, 
                      {
                        cRow <- dae::as.numfac(Row)
                        cRow <- cRow - mean(unique(cRow))
                        cColumn <- dae::as.numfac(Column)
                        cColumn <- cColumn - mean(unique(cColumn))
                      })
  ts.asr <- do.call(asreml,
                    list(yield ~ Rep + cRow + cColumn + WithinColPairs + 
                           Variety, 
                         random = ~ spl(cRow) + spl(cColumn) + 
                           dev(cRow) + dev(cColumn) + 
                           spl(cRow):cColumn + cRow:spl(cColumn) + 
                           spl(cRow):spl(cColumn),
                         residual = ~ Row:Column, 
                         data=Wheat.dat))
  ts.asrt <- as.asrtests(ts.asr, NULL, NULL, 
                         label = "Tensor spline model")
  ts.asrt <- rmboundary.asrtests(ts.asrt)
  testthat::expect_equal(nrow(ts.asrt$test.summary),3)

  ar1.ic <- infoCriteria(ar1.asrt$asreml.obj)
  ts.ic <- infoCriteria(ts.asrt$asreml.obj)
  if (ar1.ic$AIC < ts.ic$AIC)
  {
    ic.diff <- ar1.ic - ts.ic
    new.asrt <- ar1.asrt 
    new.asrt$test.summary <- addto.test.summary(ar1.asrt$test.summary, 
                                                terms = "Compare ar1 to ts", 
                                                DF = ic.diff$varDF, 
                                                AIC = ic.diff$AIC, BIC = ic.diff$BIC, 
                                                action = "Chose ar1")
  } else
  {
    ic.diff <- ts.ic - ar1.ic
    new.asrt <- ts.asrt
    new.asrt$test.summary <- addto.test.summary(ts.asrt$test.summary, 
                                                terms = "Compare ar1 to ts", 
                                                DF = ic.diff$varDF, 
                                                AIC = ic.diff$AIC, BIC = ic.diff$BIC, 
                                                action = "Chose ts")
  }
  testthat::expect_equal(nrow(new.asrt$test.summary),3)
})