ergm.bridge.dindstart.llk() called via logLik.ergm() now inherits the…

… drop= control parameter from the original call, and ergm.bridge.llr() now handles the case of from= and to= having elements that are identical but unsubtractable (e.g., Inf). This fixes a bug that was causing NaN likelihoods for some models with terms dropped due to extreme sufficient statistics.

fixes #573

R/control.ergm.R

@@ -718,7 +718,7 @@ STATIC_MCMC_CONTROLS <- c("MCMC.samplesize", "MCMC.prop", "MCMC.prop.weights", "
 ADAPTIVE_MCMC_CONTROLS <- c("MCMC.effectiveSize", "MCMC.effectiveSize.damp", "MCMC.effectiveSize.maxruns", "MCMC.effectiveSize.burnin.pval", "MCMC.effectiveSize.burnin.min", "MCMC.effectiveSize.burnin.max", "MCMC.effectiveSize.burnin.nmin", "MCMC.effectiveSize.burnin.nmax", "MCMC.effectiveSize.burnin.PC", "MCMC.effectiveSize.burnin.scl", "obs.MCMC.effectiveSize")
 PARALLEL_MCMC_CONTROLS <- c("parallel","parallel.type","parallel.version.check")
 OBS_MCMC_CONTROLS <- c("MCMC.base.samplesize", "MCMC.base.effectiveSize", "MCMC.samplesize", "MCMC.effectiveSize", "MCMC.interval", "MCMC.burnin")
-MPLE_CONTROLS <- c("MPLE.samplesize","MPLE.type","MPLE.maxit")
+MPLE_CONTROLS <- c("MPLE.samplesize", "MPLE.type", "MPLE.maxit", "drop")
 
 remap_algorithm_MCMC_controls <- function(control, algorithm){
   CTRLS <- c(SCALABLE_MCMC_CONTROLS, STATIC_MCMC_CONTROLS, ADAPTIVE_MCMC_CONTROLS) %>% keep(startsWith,"MCMC.") %>% substr(6, 10000L)

R/control.ergm.bridge.R

@@ -24,6 +24,7 @@
 #' @param bridge.nsteps Number of geometric bridges to use.
 #' @param bridge.target.se If not `NULL`, if the estimated MCMC standard error of the likelihood estimate exceeds this, repeat the bridge sampling, accumulating samples.
 #' @param bridge.bidirectional Whether the bridge sampler first bridges from `from` to `to`, then from `to` to `from` (skipping the first burn-in), etc. if multiple attempts are required.
+#' @param drop See [control.ergm()].
 #' @param MCMC.burnin Number of proposals before any MCMC sampling is done. It
 #' typically is set to a fairly large number.
 #' @param MCMC.burnin.between Number of proposals between the bridges; typically, less and less is needed as the number of steps decreases.
@@ -50,6 +51,8 @@ control.ergm.bridge<-function(bridge.nsteps=16, # Number of geometric bridges to
                               bridge.target.se=NULL,
                               bridge.bidirectional = TRUE,
 
+                              drop = TRUE,
+
                               MCMC.burnin=MCMC.interval*128,
                               MCMC.burnin.between=max(ceiling(MCMC.burnin/sqrt(bridge.nsteps)), MCMC.interval*16),
                               MCMC.interval=128,

R/control.logLik.ergm.R

@@ -24,6 +24,8 @@ control.logLik.ergm<-function(bridge.nsteps=16,
                               bridge.target.se=NULL,
                               bridge.bidirectional = TRUE,
 
+                              drop = NULL,
+
                               MCMC.burnin=NULL,
                               MCMC.interval=NULL,
                               MCMC.samplesize=NULL,

R/ergm.bridge.R

@@ -153,7 +153,7 @@ ergm.bridge.llr<-function(object, response=NULL, reference=~Bernoulli, constrain
   }
 
   ## Miscellaneous settings
-  Dtheta.Du <- (to-from)[!state[[1]]$model$etamap$offsettheta]
+  Dtheta.Du <- ifelse(mapply(identical, to, from), 0, to - from)[!state[[1]]$model$etamap$offsettheta]
 
   ## Handle target statistics, if passed.
   if(!is.null(target.stats)){
@@ -339,15 +339,15 @@ ergm.bridge.dindstart.llk<-function(object, response=NULL, constraints=~., coef,
 
   message("Fitting the dyad-independent submodel...")
   if(is.null(coef.dind)){
-    ergm.dind<-suppressMessages(suppressWarnings(ergm(dind,basis=nw,estimate="MPLE",constraints=constraints,obs.constraints=obs.constraints,eval.loglik=FALSE,control=control.ergm(drop=FALSE, term.options=control$term.options, MPLE.max.dyad.types=control$MPLE.max.dyad.types), offset.coef = offset.dind)))
+    ergm.dind <- suppressMessages(suppressWarnings(ergm(dind, basis=nw, estimate="MPLE", constraints=constraints, obs.constraints=obs.constraints, eval.loglik=FALSE, control=control.ergm(drop=control$drop, term.options=control$term.options, MPLE.max.dyad.types=control$MPLE.max.dyad.types), offset.coef=offset.dind)))
     etamap.dind <- ergm.dind$etamap
     stats.dind <- ergm.dind$nw.stats
 
     eta.dind <- ergm.eta(coef(ergm.dind), ergm.dind$etamap)[!ergm.dind$etamap$offsetmap]
     eta.dind <- ifelse(is.na(eta.dind),0,eta.dind)
     llk.dind <- ergm.dind$mple.lik
   }else{
-    mple.dind <- suppressMessages(suppressWarnings(ergmMPLE(dind, output="matrix", constraints=constraints,obs.constraints=obs.constraints, control=control.ergm(drop=FALSE, term.options=control$term.options, MPLE.max.dyad.types=control$MPLE.max.dyad.types))))
+    mple.dind <- suppressMessages(suppressWarnings(ergmMPLE(dind, output="matrix", constraints=constraints,obs.constraints=obs.constraints, control=control.ergm(drop=control$drop, term.options=control$term.options, MPLE.max.dyad.types=control$MPLE.max.dyad.types))))
     etamap.dind <- attr(ergm.dind, "etamap")
     stats.dind <- summary(dind, basis=nw)
 

tests/testthat/test-drop.R

@@ -34,22 +34,25 @@ test_that("multiple covariates", {
   samplike.m <- as.matrix(samplike, matrix.type="adjacency")
   samplike.m[4:10,4:10] <- 0
 
-  truth <- c(logit((network.edgecount(samplike)-sum(samplike.m))/(network.dyadcount(samplike)-sum(samplike.m))),Inf)
+  truth <- c(edges = logit((network.edgecount(samplike)-sum(samplike.m))/(network.dyadcount(samplike)-sum(samplike.m))),
+             edgecov.samplike.m = Inf)
 
   maxed.mple <- ergm(samplike~edges+edgecov(samplike.m))
-  expect_true(all.equal(truth, coef(maxed.mple),check.attributes=FALSE))
+  expect_equal(coef(maxed.mple), truth)
 
   maxed.mcmc <- ergm(samplike~edges+edgecov(samplike.m), control=control.ergm(force.main=TRUE, MCMLE.maxit=10))
-  expect_true(all.equal(truth, coef(maxed.mcmc), check.attributes=FALSE,tolerance=0.1))
+  expect_equal(coef(maxed.mcmc), truth, tolerance = 0.05)
+  expect_equal(logLik(maxed.mcmc), logLik(maxed.mple), tolerance = 0.05, ignore_attr = TRUE)
 
-
-  truth <- c(logit((network.edgecount(samplike)-sum(samplike.m))/(network.dyadcount(samplike)-sum(samplike.m))),-Inf)
+  truth <- c(edges = logit((network.edgecount(samplike)-sum(samplike.m))/(network.dyadcount(samplike)-sum(samplike.m))),
+             `edgecov.-samplike.m` = -Inf)
 
   mined.mple <- ergm(samplike~edges+edgecov(-samplike.m))
-  expect_true(all.equal(truth, coef(mined.mple),check.attributes=FALSE))
+  expect_equal(coef(mined.mple), truth)
 
   mined.mcmc <- ergm(samplike~edges+edgecov(-samplike.m), control=control.ergm(force.main=TRUE, MCMLE.maxit=10))
-  expect_true(all.equal(truth, coef(mined.mcmc), check.attributes=FALSE, tolerance=0.1))
+  expect_equal(coef(mined.mcmc), truth, tolerance=0.05)
+  expect_equal(logLik(mined.mcmc), logLik(mined.mple), tolerance = 0.05, ignore_attr = TRUE)
 })
 
 # This is mainly to make sure it doesn't crash for dyad-dependent