## Test good scenarios
test_that("DI works with \"E\" model", {
  data("Switzerland")

  ## E model works
  mod_E <- DI(y = "yield", prop = 4:7, 
              density = "density", DImodel = "E",
              estimate_theta = TRUE, data = Switzerland)
  expect_equal(mod_E$coefficients,
               c(8.31251807959244, 8.36919642340669, 15.2626513599391, 11.7781500742191, 3.55873211379141, -0.137755987617646, 0.83571275568028),
               ignore_attr = TRUE)
  
  mod_E_theta <- DI(y = "yield", prop = 4:7, treat = "nitrogen",
              block = "density", DImodel = "E",
              estimate_theta = TRUE, data = Switzerland)
  expect_equal(mod_E_theta$coefficients,
               c(7.71767135989568, 7.77434970370992, 14.6678046402423, 11.1833033545223, 2.79449615220446, -0.137755987617647, 0.965375640775469, 0.760301181253625),
               ignore_attr = TRUE)
  
  mod_E_theta <- DI(y = "yield", prop = 4:7, treat = "nitrogen",
                     density = "density", DImodel = "E",
                     estimate_theta = TRUE, data = Switzerland)
  expect_equal(mod_E_theta$coef,
               c(8.5452910, 8.6019694, 15.4954243, 12.0109230, 2.794496, -0.9653756, 0.1377560, 0.7603012),
               ignore_attr = TRUE)

  mod_E_theta_custom <- DI(y = "yield", prop = 4:7, treat = "nitrogen",
                            block = "density", DImodel = "E",
                            theta = 0.5, data = Switzerland)
  expect_equal(mod_E_theta_custom$coef,
               c(7.71483343, 7.77151177, 14.66496671, 11.18046542, 1.276282, -0.13775599, 1.16877634, 0.500000),
               ignore_attr = TRUE)

  # Additional special conditions to cover all situations in code
  mod_E <- DI(y = "yield", prop = 4:7, treat = "nitrogen",
               density = "density", DImodel = "E",
               extra_formula = ~p1:nitrogen, data = Switzerland)
  expect_equal(mod_E$coef,
               c(8.5985566700200202, 8.52407309971615, 15.4175280362485392, 11.9330267505285228, 5.1554594166571261, -0.6406529039456219, 0.13775598761764632, -0.4003438816546424),
               ignore_attr = TRUE)

  mod_E <- DI(y = "yield", prop = 4:7, block = "nitrogen",
               density = "density", DImodel = "E",
               extra_formula = ~p1:nitrogen, data = Switzerland)
  expect_equal(mod_E$coef,
               c(8.5985566700200202, 8.52407309971615, 15.4175280362485392, 11.9330267505285228, 5.1554594166571261, -0.6406529039456219, 0.13775598761764632, -0.4003438816546424),
               ignore_attr = TRUE)

  mod_E <- DI(y = "yield", prop = 4:7, estimate_theta = T,
               DImodel = "E", extra_formula = ~p1:nitrogen, data = Switzerland)
  expect_equal(mod_E$coef,
               c(8.79816909765492, 8.32361781247475, 15.2170727490071, 11.7325714632871, 3.221081, -1.62146559961042, 0.804511561097006),
               ignore_attr = TRUE)

  mod_E <- DI(y = "yield", prop = 4:7, theta = 0.5,
               DImodel = "E", data = Switzerland)
  expect_equal(mod_E$coef,
               c(8.43218081507944, 8.48885915889368, 15.3823140954261, 11.8978128097061, 1.313556, 0.5),
               ignore_attr = TRUE)

})

# Most combinations are covered by the test cases in autoDI
# Checking for additional cases which weren't handled there
test_that("Additional test-cases for DI function", {
  data("Switzerland")
  available_models <- c("STR", "ID", "AV", "E",
                       "FG", "ADD", "FULL")
  
  # Model with treatment and density and extra formula
  coefficients <- list("STR" = c(14.7325931464327, 0.338889885147146, -0.238739909717057, -0.0357565936016199),
                       "ID" = c(12.7191573096262, 12.8599403616208, 19.8189815384772, 16.4000664930813, 0.0789205580693944, -0.368131033135327, -0.0480451360491652),
                       "AV" = c(8.37376681890422, 8.42410991369333, 15.3126245184171, 11.8231829008885, 13.8613027010167, -0.85707214528136, 0.17586229967076, 0.00361903522292134),
                       "E" = c(8.37376681890421, 8.42410991369333, 15.3126245184171, 11.8231829008885, 5.19798851288125, -0.857072145281361, 0.17586229967076, 0.00361903522292138),
                       "FG" = c(8.04275198699745, 8.09102486468054, 15.5573333454879, 12.066277338473, 18.2468914418737, 9.70477849123567, 0.697832803984883, -0.895087245635806, 0.188314586667893, 0.00480165465840326),
                       "ADD" = c(7.9211674672672, 8.16603686490083, 14.9115344207494, 12.6630338919457, 9.9818446117302, 8.4450637556713, 9.53731593345885, -0.127299299449124, -0.915704826082002, 0.195068114253214, 0.00544305113298685),
                       "FULL" = c(7.9211674672672, 8.16603686490083, 14.9115344207494, 12.6630338919457, 9.72985117122586, 24.0010539929582, 14.0697090606876, 22.1975434375367, 12.7996579039913, 0.712959437834094, -0.915704826082003, 0.195068114253215, 0.00544305113298689))
  models <- lapply(available_models, function(int_str){
    mod <- DI(y = "yield", prop = 4:7, 
              treat = "nitrogen", density = "density",
              extra_formula = ~ plot, FG = c("G", "G", "H", "H"),
              data = Switzerland, DImodel = int_str)
    expect_equal(mod$coefficients, coefficients[[int_str]], ignore_attr = TRUE)
  })
  
  # Model with treatment and extra formula
  coefficients <- list("STR" = c(14.3798837946354, 0.472226683397857, -0.0318349230648342),
                       "ID" = c(12.3938360952603, 12.5239439263854, 19.4746603897198, 16.0474206308018, -0.128420306719741, -0.0419468753200729),
                       "AV" = c(8.557347892295, 8.61333495504708, 15.5062508191916, 12.0212104610836, 13.7602667961902, -0.753432769542902, 0.000394896949268814),
                       "E" = c(8.55734789229499, 8.61333495504708, 15.5062508191915, 12.0212104610836, 5.16010004857132, -0.753432769542904, 0.000394896949268857),
                       "FG" = c(8.24122525582619, 8.29555168019982, 15.7631312061979, 12.276795854859, 18.1385234335331, 9.569598253824, 0.616277024715319, -0.783926834356985, 0.00134354297097426),
                       "ADD" = c(8.12851127157043, 8.37632785992546, 15.1254447318632, 12.8805635191488, 9.89883746064958, 8.38792815070444, 9.49008749067943, -0.164620580041107, -0.800463673245844, 0.00185799083219398),
                       "FULL" = c(8.12851127157042, 8.37632785992547, 15.1254447318632, 12.8805635191488, 9.58970841517837, 23.8708183990981, 13.949380629015, 22.0931793897904, 12.7052010184325, 0.62840971446269, -0.800463673245844, 0.00185799083219402))
  models <- lapply(available_models, function(int_str){
    mod <- DI(y = "yield", prop = 4:7, 
              treat = "nitrogen",
              extra_formula = ~ plot, FG = c("G", "G", "H", "H"),
              data = Switzerland, DImodel = int_str)
    expect_equal(mod$coefficients, coefficients[[int_str]], ignore_attr = TRUE)
  })

  # Model with just extra formula
  coefficients <- list("STR" = c(15.0069626872072, -0.0399466155307287),
                       "ID" = c(12.4310412685474, 12.5650242020063, 19.5187625341331, 16.0945446440076, -0.0441605423413487),
                       "AV" = c(8.8269711821793, 8.906469377341, 15.8177196107957, 12.351013621998, 13.5669928746415, -0.0130359273140788),
                       "E" = c(8.8269711821793, 8.906469377341, 15.8177196107957, 12.351013621998, 5.08762232799058, -0.0130359273140788),
                       "FG" = c(8.52915340159294, 8.60795970991156, 16.0806792043195, 12.6134336707359, 17.9368964812065, 9.25954556539232, 0.523075808493804, -0.0126406843101482),
                       "ADD" = c(8.42935873119389, 8.68891801506243, 15.4524555780312, 13.2219950563477, 9.68884597830234, 8.28101845613254, 9.42265158656329, -0.19258269370149, -0.0124262124953947),
                       "FULL" = c(8.42935873119389, 8.68891801506243, 15.4524555780312, 13.2219950563477, 9.27280723825924, 23.5933910126348, 13.7114270330073, 21.9188337911023, 12.5703292102002, 0.533011696686164, -0.0124262124953947))
  models <- lapply(available_models, function(int_str){
    mod <- DI(y = "yield", prop = 4:7, 
              extra_formula = ~ plot, FG = c("G", "G", "H", "H"),
              data = Switzerland, DImodel = int_str)
    expect_equal(mod$coefficients, coefficients[[int_str]], ignore_attr = TRUE)
  })
  
  # Model with density and extra formula
  coefficients <- list("STR" = c(15.2312041965283, -0.304786490592827, -0.0420291739082572),
                       "ID" = c(12.3324765061034, 12.4706928045875, 19.4277323846698, 16.0068157424997, 0.352692096189999, -0.0465788682945809),
                       "AV" = c(8.82706133670023, 8.90655213440679, 15.8177965992039, 12.3510848417485, 13.5672214100844, -0.000539836639348497, -0.013031701489559),
                       "E" = c(8.82706133670023, 8.90655213440679, 15.8177965992038, 12.3510848417485, 5.08770802878165, -0.00053983663934812, -0.013031701489559),
                       "FG" = c(8.52998003730568, 8.60871674423055, 16.0814216246409, 12.6141218147411, 17.939046732214, 9.26200409328415, 0.524917782616992, -0.0050756675955452, -0.0126009242479091),
                       "ADD" = c(8.43055794921086, 8.69006867899454, 15.4535466148142, 13.2230264659816, 9.6908851388531, 8.282631390706, 9.42410130360552, -0.191296194190481, -0.0075371760688523, -0.0123671497001649),
                       "FULL" = c(8.43055794921086, 8.69006867899454, 15.4535466148142, 13.2230264659816, 9.27645933338346, 23.5968798902277, 13.7147526930691, 21.921896442718, 12.5732286442846, 0.535747913239408, -0.0075371760688523, -0.0123671497001648))
  models <- lapply(available_models, function(int_str){
    mod <- DI(y = "yield", prop = 4:7, 
              density = "density", FG = c("G", "G", "H", "H"),
              extra_formula = ~ plot,
              data = Switzerland, DImodel = int_str)
    expect_equal(mod$coefficients, coefficients[[int_str]], ignore_attr = TRUE)
  })
})

## Special conditions for the ADD model along with theta_CI
test_that("ADD model works with theta", {
  data("Switzerland")
  
  mod_ADD_treat <- DI(y = "yield", treat = "nitrogen",
                      DImodel = "ADD", prop = 4:7, 
                      data = Switzerland,
                      estimate_theta = TRUE)
  expect_equal(mod_ADD_treat$coef,
              c(8.20728328146857, 8.44732385981388, 15.2043119670558, 13.0134282595101, 6.04716805325336, 5.09190694619538, 5.80285734859481, -0.936494203249058, -0.941661834080826, 0.78595066963828),
              ignore_attr = TRUE)
  #8.18542119397663, 8.4347651913239, 15.1857578074747, 12.9427523389733, 9.85581823286135, 8.35831709222467, 9.46561091111296, -0.183962680694284, -0.740738874359283
  mod_ADD <- DI(y = "yield", 
                DImodel = "ADD", prop = 4:7, 
                data = Switzerland,
                estimate_theta = TRUE)
  expect_equal(mod_ADD$coef,
              c(7.8513433527586, 8.09323716891178, 14.8497165412014, 12.6422322301889, 7.32032604290266, 6.22099797598873, 7.03400285269713, -0.545604749698444, 0.856025879261015),
              ignore_attr = TRUE)
  #7.88878936443855, 8.13813336178581, 14.8891259779366, 12.6461205094352, 10.1178646560904, 8.6203635154537, 9.72765733434199, 0.0780837425347466
  
  # Confidence interval for theta
  theta_CI(mod_ADD)
})

## Special conditions for the FG model
test_that("FG model works with theta", {
  data("Switzerland")
  
  mod_FG_treat <- DI(y = "yield", treat = "nitrogen",
                      DImodel = "FG", prop = 4:7,
                      FG = c("G", "G", "H", "H"),
                      data = Switzerland,
                      estimate_theta = TRUE)
  expect_equal(mod_FG_treat$coef,
               c(8.31620467951036, 8.3728830233246, 15.8716577347265, 12.3871564490065, 10.4737628188265, 3.88086134464133, -2.04718179494063, -0.972602380260069, 0.752389619705663),
               ignore_attr = TRUE)
  #8.18542119397663, 8.4347651913239, 15.1857578074747, 12.9427523389733, 9.85581823286135, 8.35831709222467, 9.46561091111296, -0.183962680694284, -0.740738874359283
  mod_FG <- DI(y = "yield", 
                DImodel = "FG", prop = 4:7, 
                FG = c("G", "G", "H", "H"),
                data = Switzerland,
                estimate_theta = TRUE)
  expect_equal(mod_FG$coef,
               c(7.94604929097243, 8.00272763478668, 15.4934048896756, 12.0089036039556, 12.8506780250522, 5.7074593892259, -1.04515673145053, 0.829932131271245),
               ignore_attr = TRUE)
  #7.88878936443855, 8.13813336178581, 14.8891259779366, 12.6461205094352, 10.1178646560904, 8.6203635154537, 9.72765733434199, 0.0780837425347466
  
})

## Ensure specific models can't be fit for models with less than four species
test_that("Correct number of species are present", {
  data("sim1")
  
  sim1$p5 <- 1 - sim1$p1
  sim1$p6 <- 1 - sim1$p1 - sim1$p2
  
  expect_error(DI(data = sim1, y = "response", 
                  prop = c("p1", "p5"), 
                  DImodel= "E"),
               "you must have > 2 species to fit model E")

  expect_error(DI(data = sim1, y = "response", 
                  prop = c("p1", "p5"), 
                  DImodel= "AV"),
               "you must have > 2 species to fit model AV")
  
  expect_error(DI(data = sim1, y = "response", 
                  prop = c("p1", "p5"), 
                  DImodel= "FG",
                  FG = c("G", "G")),
               "you must have > 2 species to fit model FG")
  
  expect_error(DI(data = sim1, y = "response", 
                  prop = c("p1", "p2", "p6"), 
                  DImodel= "ADD"),
               "you must have > 3 species to fit model ADD")

})

## Richness vs DI works
test_that("Richness vs DI works", {
  data("Switzerland")
  
  ## compare the richness model with DI alternatives
  t1 <- richness_vs_DI(y = "yield", prop = 4:7, data = Switzerland)
  m1 <- DI(y = "yield", prop = 4:7, data = Switzerland,
           DImodel = "AV", estimate_theta = TRUE)
  expect_equal(t1$coefficients, m1$coefficients)
    
  ## include the density effects in the linear predictors of the four models
  t2 <- richness_vs_DI(y = "yield", prop = 4:7, data = Switzerland, extra_formula = ~ density)
  m2 <- DI(y = "yield", prop = 4:7, data = Switzerland, DImodel = "AV",
           extra_formula = ~ density, estimate_theta = TRUE)
  expect_equal(t2$coefficients, m2$coefficients)
})

## S3 methods for class DI (extract, AIC, BIC)
test_that("S3 methods work", {
  data("Switzerland")
  
  mod <- DI(y = "yield", prop = 4:7, treat = "nitrogen", FG = c("G", "G", "H", "H"), 
            density = "density", DImodel = "FG", data = Switzerland)
  
  # Extract function
  int_terms <- DI_data(prop = 4:7, data = Switzerland, FG = c("G", "G", "H","H"))
  int_terms[3:5] <- lapply(int_terms[3:5], as.data.frame)
  expect_equal(extract(mod),
               int_terms,
               ignore_attr = TRUE)
  
  # AIC 
  expect_equal(AIC(mod), 260.57313407399)
  
  # AICc 
  expect_equal(AICc(mod), 264.43278)
  
  # BIC
  expect_equal(BIC(mod), 282.768211125751)
  
  # BICc
  expect_equal(BICc(mod), 290.911120732231)
  
})

## Custom formula
test_that("Custom formula example", {
  data("Switzerland")
  # Basic custom formula
  mod <- DI(custom_formula = "yield ~ p1 + p2 + p3 + p4",
            data = Switzerland)
  lm_ver <- stats::glm(yield ~ p1 + p2 + p3 + p4, 
                       data = Switzerland)
  expect_equal(mod$coef, 
               lm_ver$coef)
})

## Special scenarios which can throw errors
test_that("DI fails with appropriate message", {
  data("sim0")
  # To suppress rounding error warning
  sim0$p3 <- 1- sim0$p1 - sim0$p2
  
  # y can't be missing
  expect_error(DI(prop = 3:5, data = sim0),
               regexp = "You must supply a response variable name or column index through the argument 'y'")
  
  # Theta can't be negative
  expect_error(DI(y = "response", prop = 3:5, DImodel = "E",
                  data = sim0, theta = -1),
               regexp = "Please choose a positive value for theta")
  
  # Can't specify both custom and extra formula
  expect_error(DI(custom_formula = "response ~ p1 + p2 + p3",
                  extra_formula = "~ + p1:p2", data = sim0),
               regexp = "Please provide either custom_formula or extra_formula; not the two at the same time.")
  
  # Can't fit ADD model with 3 species
  expect_error(DI(y = "response", prop = 3:5, DImodel = "ADD",
                  data = sim0),
               regexp = "> 3 species")
  
  # Can't fit AV and E model with 2 species
  sim0$p4 <- 1 - sim0$p1
  expect_error(DI(y = "response", prop = c("p1", "p4"), 
                  DImodel = "AV", data = sim0),
               regexp = "> 2 species")
  expect_error(DI(y = "response", prop = c("p1", "p4"), 
                  DImodel = "E", data = sim0),
               regexp = "> 2 species")
  
  # FG should be specified when fitting FG model
  expect_error(DI(y = "response", prop = 3:5, DImodel = "FG",
                  data = sim0),
               regexp = "The argument FG must be specified alongside DImodel = 'FG'")
  
  data("Switzerland")
  # Can't estimate theta when using custom_formula
  expect_error(DI(custom_formula = "yield ~ p1 + p2 + p3 + p4",
                  data = Switzerland, estimate_theta = TRUE),
               regexp = "theta estimation not available when custom_formula is supplied")
  
  ## Warnings
  # Don't specify theta and estimate_theta, estimate_theta takes precedence 
  expect_warning(DI(y = "yield", prop = 4:7, DImodel = AV,
                  data = Switzerland, theta = 0.5, estimate_theta = TRUE),
                regexp = "By specifying estimate_theta as TRUE, DI is overriding the specified theta value")
  
  # If custom_formula and DImodel both specified custom_formula takes precendence
  expect_warning(DI(y = "yield", prop = 4:7, DImodel = "FULL",
                    data = Switzerland, custom_formula = "yield ~ p1 + p2"),
                 regexp = "fitting custom DI model using supplied custom_formula instead of model")
  
})  

## Testing for reference model and model comparison along with anovaglm
test_that("Interior functions in DI_internal work", {
  data("Switzerland")
  
  mod <- DI(y = "yield", prop = 4:7, treat = "nitrogen", FG = c("G", "G", "H", "H"), 
            density = "density", DImodel = "FG", data = Switzerland,
            estimate_theta = TRUE)
  
  # DI_reference and DI_compare
  expect_equal(DI_reference,
               DI_compare(mod))
  
  # anovaDIglm
  obj <- anovaDIglm(mod)
  expect_equal(obj$`Resid. Dev`,
               c(13310.9015415132, 9103.40986879718, 6326.62608653603, 2402.17197250788, 407.144041330928, 145.443529829053, 140.640221756244, 139.663529974278, 127.40531115179, 127.082707045673),
               ignore_attr = TRUE)
  
  # Rao test
  obj <- anovaDIglm(mod, test = "Rao")
  expect_equal(obj$`Rao`,
               c(NA, 4207.491672716, 2776.78378226115, 3924.45411402815, 1995.02793117695, 261.700511501874, 4.80330807280922, 0.976691781966224, 12.258218822488, 0.322604106116728),
               ignore_attr = TRUE)
  
  ## anova_glmList
  mod1 <- DI(y = "yield", prop = 4:7, treat = "nitrogen",  
             density = "density", DImodel = "AV", 
             data = Switzerland,
             estimate_theta = TRUE)
  mod2 <- DI(y = "yield", prop = 4:7, treat = "nitrogen",  
             density = "density", DImodel = "FULL", 
             data = Switzerland,
             estimate_theta = TRUE)
  
  expect_equal(anova_glmlist(list(mod1, mod, mod2))$Deviance,
               c(NA, 37.8055122016968, 11.4973620431856))
  
  expect_equal(anova_glmlist(list(mod1, mod, mod2), test = "Rao")$Rao,
               c(NA, 37.8141225527972, 11.5427716510587))
})