add test to see if "bestmodel" violates hereditary condition

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: BAS
 Version: 1.7.1.9000
-Date: 2023-12-6
+Date: 2023-12-7
 Title: Bayesian Variable Selection and Model Averaging using Bayesian Adaptive Sampling
 Authors@R: c(person("Merlise", "Clyde", email="[email protected]",
 	   		       role=c("aut","cre", "cph"),

R/bas_glm.R

@@ -430,11 +430,11 @@ bas.glm <- function(formula, family = binomial(link = "logit"),
 
   if (!inherits(betaprior, "prior")) stop("prior on coeeficients must be an object of type 'prior'")
 
-  loglik_null <- as.numeric(-0.5 * glm(Y ~ 1,
-                                       weights = weights,
-                                       offset = offset,
-                                       family = eval(call$family)
-                                       )$null.deviance)
+  null.deviance = glm(Y ~ 1,
+                      weights = weights,
+                      offset = offset,
+                      family = eval(call$family))$null.deviance
+  loglik_null <- as.numeric(-0.5 * null.deviance)
 
   betaprior$hyper.parameters$loglik_null <- loglik_null
   #  	browser()
@@ -542,7 +542,7 @@ bas.glm <- function(formula, family = binomial(link = "logit"),
 
   if (betaprior$class == "IC") df <- df - result$size + 1
   result$df <- df
-  result$R2 <- .R2.glm.bas(result$deviance, result$size, call)
+  result$R2 <- 1.0 - result$deviance/null.deviance
   result$n.vars <- p
   result$Y <- Yvec
   result$X <- X
@@ -637,19 +637,3 @@ bas.glm <- function(formula, family = binomial(link = "logit"),
   return(probs)
 }
 
-.R2.glm.bas <- function(deviance, size, call) {
-  n.models <- length(deviance)
-  null.model <- (1:n.models)[size == 1]
-  if (is.null(null.model)) {
-    null.deviance <- glm(eval(call$formula),
-      data = eval(call$data),
-      family = eval(call$family)
-    )$null.deviance
-  }
-  else {
-    null.deviance <- deviance[null.model]
-  }
-
-  R2 <- 1 - deviance / null.deviance
-  return(R2)
-}

tests/testthat/test-bas-glm.R

@@ -422,11 +422,19 @@ test_that("include always MCMC", {
   data("Pima.tr", package="MASS")
   pima_BAS = bas.glm(type ~ .,
                        data = Pima.tr, method = "MCMC",
-                       include.always = ~ bp,
+                       include.always = type ~ bp, 
                        betaprior = g.prior(g=100), family = binomial(),
                        modelprior = beta.binomial(1, 1))
   x = pima_BAS$probne0[match(c("Intercept", "bp") ,pima_BAS$namesx)]
   expect_equal(2, sum(x), tolerance=.002)
+
+  expect_equal(0, sum(pima_BAS$R2 < 0))
+
+ expect_warning(bas.glm(type ~ poly(bp,2),
+                     data = Pima.tr, method = "BAS",
+                     force.heredity = TRUE, bestmodel = c(1,0,1),
+                     betaprior = g.prior(g=100), family = binomial(),
+                     modelprior = beta.binomial(1, 1)))
 
 #  pima_BAS = bas.glm(type ~ .,
 #                     data = Pima.tr, method = "BAS",