Nothing
tests/testthat/test-DI.R
In DImodels: Diversity-Interactions (DI) Models
## 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
# DImodel should be proper
expect_error(DI(prop = 3:5, data = sim0, DImodel = "AV1"),
regexp = "should be one of")
# 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_no <- DI(y = "yield", prop = 4:7, treat = "nitrogen", FG = c("G", "G", "H", "H"),
density = "density", DImodel = "FG", data = Switzerland,
estimate_theta = FALSE)
# DI_reference and DI_compare
expect_equal(DI_reference(mod_no, 4:7, Switzerland),
DI_compare(mod_no))
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(mod, 4:7, Switzerland, mod$coefficients["theta"]),
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))
})
# Test describe_model
test_that("describe_model works", {
data(sim2)
## Fit model
mod_FG <- DI(y = "response", FG = c("G", "G", "L", "L"),
prop = 3:6, data = sim2, DImodel = "FG")
## Describe model
expect_equal(describe_model(mod_FG),
"This model has 4 species with the Functional group interaction structure with FG values of (G, G, L, L).")
mod_FULL <- DI(y = "response", estimate_theta = TRUE,
prop = 3:6, data = sim2, DImodel = "FULL")
expect_equal(describe_model(mod_FULL),
"This model has 4 species with the Full pairwise interaction structure. The interaction terms have an associated theta (exponent on the interactions) value of '0.45'.")
mod_AV <- DI(y = "response", ID = c("ID1", "ID1", "ID2", "ID2"),
estimate_theta = TRUE,
prop = 3:6, data = sim2, DImodel = "AV")
expect_equal(describe_model(mod_AV),
"This model has 4 species with the Average interaction structure. The interaction terms have an associated theta (exponent on the interactions) value of '0.45'.The species identity effects are grouped as (ID1, ID1, ID2, ID2).")
mod_STR <- DI(y = "response", FG = c("G", "G", "L", "L"),
prop = 3:6, data = sim2, DImodel = "STR")
## Describe model
expect_equal(describe_model(mod_STR),
"This model doesn't have any species identity or interaction effects and only consists of experimental structures.")
mod_CUST <- DI(custom_formula = response ~ block, data = sim2)
## Describe model
expect_equal(describe_model(mod_CUST),
"This is a custom DI model.")
# Proper error is thrown
expect_error(describe_model(lm(1 ~ 1)),
regexp = "`model` should be regression object of class <DI>")
})
Try the DImodels package in your browser
Any scripts or data that you put into this service are public.
DImodels documentation built on May 29, 2024, 7:05 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## 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
# DImodel should be proper
expect_error(DI(prop = 3:5, data = sim0, DImodel = "AV1"),
regexp = "should be one of")
# 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_no <- DI(y = "yield", prop = 4:7, treat = "nitrogen", FG = c("G", "G", "H", "H"),
density = "density", DImodel = "FG", data = Switzerland,
estimate_theta = FALSE)
# DI_reference and DI_compare
expect_equal(DI_reference(mod_no, 4:7, Switzerland),
DI_compare(mod_no))
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(mod, 4:7, Switzerland, mod$coefficients["theta"]),
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))
})
# Test describe_model
test_that("describe_model works", {
data(sim2)
## Fit model
mod_FG <- DI(y = "response", FG = c("G", "G", "L", "L"),
prop = 3:6, data = sim2, DImodel = "FG")
## Describe model
expect_equal(describe_model(mod_FG),
"This model has 4 species with the Functional group interaction structure with FG values of (G, G, L, L).")
mod_FULL <- DI(y = "response", estimate_theta = TRUE,
prop = 3:6, data = sim2, DImodel = "FULL")
expect_equal(describe_model(mod_FULL),
"This model has 4 species with the Full pairwise interaction structure. The interaction terms have an associated theta (exponent on the interactions) value of '0.45'.")
mod_AV <- DI(y = "response", ID = c("ID1", "ID1", "ID2", "ID2"),
estimate_theta = TRUE,
prop = 3:6, data = sim2, DImodel = "AV")
expect_equal(describe_model(mod_AV),
"This model has 4 species with the Average interaction structure. The interaction terms have an associated theta (exponent on the interactions) value of '0.45'.The species identity effects are grouped as (ID1, ID1, ID2, ID2).")
mod_STR <- DI(y = "response", FG = c("G", "G", "L", "L"),
prop = 3:6, data = sim2, DImodel = "STR")
## Describe model
expect_equal(describe_model(mod_STR),
"This model doesn't have any species identity or interaction effects and only consists of experimental structures.")
mod_CUST <- DI(custom_formula = response ~ block, data = sim2)
## Describe model
expect_equal(describe_model(mod_CUST),
"This is a custom DI model.")
# Proper error is thrown
expect_error(describe_model(lm(1 ~ 1)),
regexp = "`model` should be regression object of class <DI>")
})
Try the DImodels package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.