From 2d0a718ac1226605e071df630a05d1534bfeec85 Mon Sep 17 00:00:00 2001
From: "Mark S. Handcock" <handcock@ucla.edu>
Date: Wed, 30 Jun 2021 17:56:29 -0700
Subject: [PATCH 1/8] This is a transference of the ergm-private"
 "control_llik" branch to the public version ergm 4.0.

The branch consolidates control arguments to the likelihood functions used in ergm.estimate into a single list variable called "control.llik".
This simplifies the calls to the functions.

The branch also adds a simple hook in "ergm.estimate" for the Kpenalty likelihood functions. These are required for the "ergm'tapered" package to interoperate with "ergm".

Also fixed a few documentation typos in "ergm".
---
 DESCRIPTION                  |  1 +
 R/control.ergm.R             | 13 ++++++++-
 R/control.simulate.R         |  6 +++-
 R/ergm.CD.fixed.R            |  9 ++++--
 R/ergm.MCMLE.R               | 13 +++++----
 R/ergm.estimate.R            | 47 ++++++++++++++----------------
 R/ergm.getMCMCsample.R       |  4 +--
 R/ergm.llik.R                | 55 ++++++++++++++++++------------------
 R/ergm.llik.obs.R            | 36 +++++++++++------------
 R/ergm.robmon.R              |  6 ++++
 R/ergm.stocapprox.R          |  7 +++++
 R/ergm_estfun.R              | 28 +++++++++++++-----
 R/logLik.ergm.R              |  2 +-
 R/parallel.utils.R           |  2 +-
 R/summary.ergm.R             |  2 +-
 R/zzz.R                      |  2 +-
 man/control.ergm.Rd          | 14 +++++++--
 man/control.simulate.ergm.Rd |  9 +++++-
 man/ergm-terms.Rd            |  2 +-
 man/ergm.estfun.Rd           | 10 +++----
 20 files changed, 165 insertions(+), 103 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index f907ac4a8..f068c4076 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -48,6 +48,7 @@ Suggests:
     tergm (>= 4.0.0),
     ergm.count (>= 4.0),
     ergm.userterms (>= 3.10.0),
+    ergm.tapered (>= 1.0.0),
     networkDynamic (>= 0.10.1),
     covr
 SystemRequirements: OpenMPI
diff --git a/R/control.ergm.R b/R/control.ergm.R
index f45b0824e..82092f3ad 100644
--- a/R/control.ergm.R
+++ b/R/control.ergm.R
@@ -277,6 +277,8 @@
 #' may increase the target ESS to reduce the MCMC standard error.
 #' @param MCMLE.metric Method to calculate the loglikelihood approximation.
 #' See Hummel et al (2010) for an explanation of "lognormal" and "naive".
+#' "Kpenalty" is used to specify the use of a tapering model with a penalty for
+#' kurtosis. It is used by the \code{ergm.tapered} package and is not usually user specified.
 #' @param MCMLE.method Deprecated. By default, ergm uses \code{trust}, and
 #' falls back to \code{optim} with Nelder-Mead method when \code{trust} fails.
 #' @param MCMLE.dampening (logical) Should likelihood dampening be used?
@@ -427,6 +429,10 @@
 #'   for CD.
 #' 
 #' @param loglik See \code{\link{control.ergm.bridge}}
+#' @param MCMC.esteq.exclude.statistics character string pattern. statistics with names that contain this string will be excluded from 
+#' the estimating equations. For example, this supports tapering models which may have terms with non-specific means.
+#' @param MCMLE.varweight numeric multiplier for covariance term in loglikelihood
+#' where 0.5 is the "true" value but this can be increased or decreased.
 #' @template term_options
 #' @template control_MCMC_parallel
 #' @template seed
@@ -502,6 +508,7 @@ control.ergm<-function(drop=TRUE,
                        MCMC.maxedges=Inf,
                        MCMC.addto.se=TRUE,
                        MCMC.packagenames=c(),
+                       MCMC.esteq.exclude.statistics=NULL,
 
                        SAN.maxit=4,
                        SAN.nsteps.times=8,
@@ -542,12 +549,13 @@ control.ergm<-function(drop=TRUE,
 
                        MCMLE.min.depfac=2,
                        MCMLE.sampsize.boost.pow=0.5,
+                       MCMLE.varweight=0.5,
 
                        MCMLE.MCMC.precision=if(startsWith("confidence", MCMLE.termination[1])) 0.1 else 0.005,
                        MCMLE.MCMC.max.ESS.frac=0.1,
                        MCMLE.metric=c("lognormal", "logtaylor",
                          "Median.Likelihood",
-                         "EF.Likelihood", "naive"),
+                         "EF.Likelihood", "Kpenalty", "naive"),
                        MCMLE.method=c("BFGS","Nelder-Mead"),
                        MCMLE.dampening=FALSE,
                        MCMLE.dampening.min.ess=20,
@@ -691,6 +699,9 @@ ADAPTIVE_MCMC_CONTROLS <- c("MCMC.effectiveSize", "MCMC.effectiveSize.damp", "MC
 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.max.dyad.types","MPLE.samplesize","MPLE.type","MPLE.maxit")
+STATIC_MCMLE_CONTROLS <- c("MCMC.esteq.exclude.statistics", "MCMLE.varweight", 
+                           "MCMLE.dampening", "MCMLE.dampening.min.ess", "MCMLE.dampening.level")
+STATIC_TAPERING_CONTROLS <- c("MCMLE.kurtosis.prior", "MCMLE.kurtosis.location", "MCMLE.kurtosis.scale", "MCMLE.kurtosis.penalty")
 
 remap_algorithm_MCMC_controls <- function(control, algorithm){
   CTRLS <- c(SCALABLE_MCMC_CONTROLS, STATIC_MCMC_CONTROLS, ADAPTIVE_MCMC_CONTROLS) %>% keep(startsWith,"MCMC.") %>% substr(6, 10000L)
diff --git a/R/control.simulate.R b/R/control.simulate.R
index 668bc5bca..bf99e45aa 100644
--- a/R/control.simulate.R
+++ b/R/control.simulate.R
@@ -47,6 +47,8 @@
 #' @template term_options
 #' @template control_MCMC_parallel
 #' @template control_MCMC_packagenames
+#' @param MCMC.esteq.exclude.statistics character string pattern. statistics with names that contain this string will be excluded from 
+#' the estimating equations. For example, this supports tapering models which may have terms with non-specific means.
 #' @template control_dots
 #' @return A list with arguments as components.
 #' @seealso \code{\link{simulate.ergm}}, \code{\link{simulate.formula}}.
@@ -73,6 +75,7 @@ control.simulate.formula.ergm<-function(MCMC.burnin=MCMC.interval*16,
                                         
                                         MCMC.maxedges=Inf,
                                         MCMC.packagenames=c(),
+                                        MCMC.esteq.exclude.statistics=NULL,
                                         
                                         MCMC.runtime.traceplot=FALSE,  
                                         network.output="network",
@@ -108,7 +111,7 @@ control.simulate.formula<-control.simulate.formula.ergm
 #' @description While the others supply a full set of simulation
 #'   settings, `control.simulate.ergm` when passed as a control
 #'   parameter to [simulate.ergm()] allows some settings to be
-#'   inherited from the ERGM stimation while overriding others.
+#'   inherited from the ERGM estimation while overriding others.
 #' 
 #' @export control.simulate.ergm
 control.simulate.ergm<-function(MCMC.burnin=NULL,
@@ -128,6 +131,7 @@ control.simulate.ergm<-function(MCMC.burnin=NULL,
                                 
                                 MCMC.maxedges=Inf,
                                 MCMC.packagenames=NULL,
+                                MCMC.esteq.exclude.statistics=NULL,
                                 
                                 MCMC.runtime.traceplot=FALSE,
                                 network.output="network",
diff --git a/R/ergm.CD.fixed.R b/R/ergm.CD.fixed.R
index c364cac77..99fd35379 100644
--- a/R/ergm.CD.fixed.R
+++ b/R/ergm.CD.fixed.R
@@ -122,6 +122,11 @@ ergm.CD.fixed <- function(init, nw, model,
   mcmc.init <- init
   finished <- FALSE
 
+  control.llik <- list(MCMLE.metric=control$CD.metric,
+                       MCMLE.varweight=control$MCMLE.varweight, MCMLE.dampening=control$CD.dampening,
+                       MCMLE.dampening.min.ess=control$CD.dampening.min.ess,
+                       MCMLE.dampening.level=control$CD.dampening.level)
+
   for(iteration in 1:control$CD.maxit){
     if(iteration == control$CD.maxit) finished <- TRUE
     if(verbose){
@@ -233,10 +238,8 @@ ergm.CD.fixed <- function(init, nw, model,
                        calc.mcmc.se=FALSE,
                        hessianflag=control$main.hessian,
                        method=control$CD.method,
-                       dampening=control$CD.dampening,
-                       dampening.min.ess=control$CD.dampening.min.ess,
-                       dampening.level=control$CD.dampening.level,
                        metric=control$CD.metric,
+                       control.llik=control.llik,
                        steplen=steplen,
                        verbose=verbose,
                        estimateonly=!finished)
diff --git a/R/ergm.MCMLE.R b/R/ergm.MCMLE.R
index 2f10003ff..3c385a796 100644
--- a/R/ergm.MCMLE.R
+++ b/R/ergm.MCMLE.R
@@ -211,6 +211,11 @@ ergm.MCMLE <- function(init, nw, model,
     rm(state)
   }
 
+  control.llik <- list()
+  control.transfer <- c(STATIC_MCMLE_CONTROLS, STATIC_TAPERING_CONTROLS)
+  for(arg in control.transfer)
+    if(!is.null(control[[arg]]))
+      control.llik[arg] <- list(control[[arg]])
   
   for(iteration in 1:control$MCMLE.maxit){
     if(verbose){
@@ -292,7 +297,7 @@ ergm.MCMLE <- function(init, nw, model,
     }
 
     # Compute the sample estimating equations and the convergence p-value. 
-    esteqs <- ergm.estfun(statsmatrices, theta=mcmc.init, model=model)
+    esteqs <- ergm.estfun(statsmatrices, theta=mcmc.init, model=model, exclude=control$MCMC.esteq.exclude.statistics)
     esteq <- as.matrix(esteqs)
     if(isTRUE(all.equal(apply(esteq,2,stats::sd), rep(0,ncol(esteq)), check.names=FALSE))&&!all(esteq==0))
       stop("Unconstrained MCMC sampling did not mix at all. Optimization cannot continue.")
@@ -302,7 +307,7 @@ ergm.MCMLE <- function(init, nw, model,
                              nonident_action = control$MCMLE.nonident,
                              nonvar_action = control$MCMLE.nonvar)
 
-    esteqs.obs <- if(obs) ergm.estfun(statsmatrices.obs, theta=mcmc.init, model=model) else NULL
+    esteqs.obs <- if(obs) ergm.estfun(statsmatrices.obs, theta=mcmc.init, model=model, exclude=control$MCMC.esteq.exclude.statistics) else NULL
     esteq.obs <- if(obs) as.matrix(esteqs.obs) else NULL
 
     # Update the interval to be used.
@@ -411,10 +416,8 @@ ergm.MCMLE <- function(init, nw, model,
                        calc.mcmc.se=control$MCMLE.termination == "precision" || (control$MCMC.addto.se && last.adequate) || iteration == control$MCMLE.maxit,
                        hessianflag=control$main.hessian,
                        method=control$MCMLE.method,
-                       dampening=control$MCMLE.dampening,
-                       dampening.min.ess=control$MCMLE.dampening.min.ess,
-                       dampening.level=control$MCMLE.dampening.level,
                        metric=control$MCMLE.metric,
+                       control.llik=control.llik,
                        steplen=steplen, steplen.point.exp=control$MCMLE.steplength.point.exp,
                        verbose=verbose,
                        estimateonly=!calc.MCSE)
diff --git a/R/ergm.estimate.R b/R/ergm.estimate.R
index 751d8aac9..df626d9f1 100644
--- a/R/ergm.estimate.R
+++ b/R/ergm.estimate.R
@@ -37,10 +37,10 @@
 #   trace           : a non-negative interger specifying how much tracing
 #                     information should be printed by the <optim> routine;
 #                     default=6*'verbose'
-#   dampening       : (logical) should likelihood dampening be used?
-#  dampening.min.ess: effective sample size below which dampening is used
-#   dampening.level : proportional distance from boundary of the convex hull
-#                     move
+#   control.llik    : A list of control values for the log-likelihood functions. These include:
+#     MCMLE.dampening       : (logical) should likelihood dampening be used?
+#    MCMLE.dampening.min.ess: effective sample size below which dampening is used
+#     MCMLE.dampening.level : proportional distance from boundary of the convex hull
 #   estimateonly    : whether only the estimates (vs. the estimates and the
 #                     standard errors) should be calculated; default=FALSE
 #
@@ -61,9 +61,7 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
                         method="Nelder-Mead",
                         calc.mcmc.se=TRUE, hessianflag=TRUE,
                         verbose=FALSE, trace=6*verbose,
-                        dampening=FALSE,
-                        dampening.min.ess=100,
-                        dampening.level=0.1,
+                        control.llik=control.logLik.ergm(),
                         steplen=1, steplen.point.exp=1,
                         cov.type="normal",# cov.type="robust", 
                         estimateonly=FALSE, ...) {
@@ -147,9 +145,6 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
     
   # Choose appropriate loglikelihood, gradient, and Hessian functions
   # depending on metric chosen and also whether obsprocess==TRUE
-  # Also, choose varweight multiplier for covariance term in loglikelihood
-  # where 0.5 is the "true" value but this can be increased or decreased
-  varweight <- 0.5
   if (verbose) { message("Using ", metric, " metric (see control.ergm function).") }
   if (obsprocess) {
     loglikelihoodfn <- switch(metric,
@@ -158,6 +153,7 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
                               logtaylor=llik.fun.obs.lognormal,
                               Median.Likelihood=llik.fun.obs.robust,
                               EF.Likelihood=llik.fun.obs.lognormal,
+                              Kpenalty=ergm.tapered::llik.fun.obs.Kpenalty,
                               llik.fun.obs.IS)
     gradientfn <- switch(metric,
                          Likelihood=llik.grad.obs.IS,
@@ -165,6 +161,7 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
                          logtaylor=llik.grad.obs.IS,
                          Median.Likelihood=llik.grad.obs.IS,
                          EF.Likelihood=llik.grad.obs.IS,
+                         Kpenalty=ergm.tapered::llik.grad.obs.Kpenalty.numDeriv,
                          llik.grad.obs.IS)
     Hessianfn <- switch(metric,
                         Likelihood=llik.hessian.obs.IS,
@@ -172,6 +169,7 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
                         logtaylor=llik.hessian.obs.IS,
                         Median.Likelihood=llik.hessian.obs.IS,
                         EF.Likelihood=llik.hessian.obs.IS,
+                        Kpenalty=ergm.tapered::llik.hessian.obs.Kpenalty.numDeriv,
                         llik.hessian.obs.IS)
   } else {
     loglikelihoodfn <- switch(metric,
@@ -180,6 +178,7 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
                               logtaylor=llik.fun.logtaylor,
                               Median.Likelihood=llik.fun.median,
                               EF.Likelihood=llik.fun.EF,
+                              Kpenalty=ergm.tapered::llik.fun.Kpenalty,
                               llik.fun.IS)
     gradientfn <- switch(metric,
                          Likelihood=llik.grad.IS,
@@ -187,6 +186,7 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
                          logtaylor=llik.grad.IS,
                          Median.Likelihood=llik.grad.IS,
                          EF.Likelihood=llik.grad.IS,
+                         Kpenalty=ergm.tapered::llik.grad.Kpenalty.numDeriv,
                          llik.grad.IS)
     Hessianfn <- switch(metric,
                         Likelihood=llik.hessian.IS,
@@ -194,6 +194,7 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
                         logtaylor=llik.hessian.IS,
                         Median.Likelihood=llik.hessian.IS,
                         EF.Likelihood=llik.hessian.IS,
+                        Kpenalty=ergm.tapered::llik.hessian.Kpenalty.numDeriv,
                         llik.hessian.IS)
   }
   
@@ -275,11 +276,8 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
                       parscale=rep(1,length(guess)), minimize=FALSE,
                       xsim=xsim,
                       xsim.obs=xsim.obs,
-                      varweight=varweight,
-                      dampening=dampening,
-                      dampening.min.ess=dampening.min.ess,
-                      dampening.level=dampening.level,
-                      eta0=eta0, etamap=etamap.no),
+                      eta0=eta0, etamap=etamap.no,
+                      control.llik=control.llik),
             silent=FALSE)
     Lout$par<-Lout$argument
     if(inherits(Lout,"try-error")) {
@@ -292,11 +290,8 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
                                      reltol=nr.reltol),
                         xsim=xsim,
                         xsim.obs=xsim.obs,
-                        varweight=varweight,
-                        dampening=dampening,
-                        dampening.min.ess=dampening.min.ess,
-                        dampening.level=dampening.level,
-                        eta0=eta0, etamap=etamap.no),
+                        eta0=eta0, etamap=etamap.no,
+                        control.llik=control.llik),
               silent=FALSE)
       if(inherits(Lout,"try-error")){
         message(paste("No direct MLE exists!"))
@@ -323,8 +318,8 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
     gradienttheta <- llik.grad.IS(theta=Lout$par,
                         xsim=xsim,
                         xsim.obs=xsim.obs,
-                        varweight=varweight,
-                        eta0=eta0, etamap=etamap.no)
+                        eta0=eta0, etamap=etamap.no,
+                        control.llik=control.llik)
     gradient <- rep(NA, length=length(init))
     gradient[!model$etamap$offsettheta] <- gradienttheta
     #
@@ -338,9 +333,8 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
       Lout$hessian <- Hessianfn(theta=Lout$par,
                         xsim=xsim.orig,
                         xsim.obs=xsim.orig.obs,
-                        varweight=varweight,
-                        eta0=eta0, etamap=etamap.no
-                        )
+                        eta0=eta0, etamap=etamap.no,
+                        control.llik=control.llik)
     }
     
     covar <- matrix(NA, ncol=length(theta), nrow=length(theta))
@@ -354,7 +348,8 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
     if(calc.mcmc.se){
       if (verbose) message("Starting MCMC s.e. computation.")
       mcse.metric <-
-        if ((metric == "lognormal" || metric == "Likelihood") && length(model$etamap$curved) == 0) "lognormal"
+        if ((metric %in% c("lognormal", "Likelihood", "Kpenalty")) 
+          && length(model$etamap$curved) == 0) "lognormal"
         else "IS"
       mc.cov <- ergm.MCMCse(model = model, theta = theta, init = init,
                             statsmatrices = statsmatrices,
diff --git a/R/ergm.getMCMCsample.R b/R/ergm.getMCMCsample.R
index ba1882c1a..38212560f 100644
--- a/R/ergm.getMCMCsample.R
+++ b/R/ergm.getMCMCsample.R
@@ -207,7 +207,7 @@ ergm_MCMC_sample <- function(state, control, theta=NULL,
         if(verbose) message("Increasing thinning to ",interval,".")
       }
       
-      esteq <- lapply(sms, function(sm) NVL3(theta, ergm.estfun(sm, ., as.ergm_model(state[[1]])), sm[,!as.ergm_model(state[[1]])$etamap$offsetmap,drop=FALSE])) %>%
+      esteq <- lapply(sms, function(sm) NVL3(theta, ergm.estfun(sm, ., as.ergm_model(state[[1]]), exclude=control$MCMC.esteq.exclude.statistics), sm[,!as.ergm_model(state[[1]])$etamap$offsetmap,drop=FALSE])) %>%
         lapply.mcmc.list(mcmc, start=1, thin=interval) %>% lapply.mcmc.list(`-`)
 
       if(control.parallel$MCMC.runtime.traceplot){
@@ -260,7 +260,7 @@ ergm_MCMC_sample <- function(state, control, theta=NULL,
     if(!is.null(nws)) nws <- map(outl, "saved")
     
     if(control.parallel$MCMC.runtime.traceplot){
-      lapply(sms, function(sm) NVL3(theta, ergm.estfun(sm, ., as.ergm_model(state[[1]])), sm[,!as.ergm_model(state[[1]])$etamap$offsetmap,drop=FALSE])) %>% lapply(mcmc, start=control.parallel$MCMC.burnin+1, thin=control.parallel$MCMC.interval) %>% as.mcmc.list() %>% window(., thin=thin(.)*max(1,floor(niter(.)/1000))) %>% plot(ask=FALSE,smooth=TRUE,density=FALSE)
+      lapply(sms, function(sm) NVL3(theta, ergm.estfun(sm, ., as.ergm_model(state[[1]]), exclude=control$MCMC.esteq.exclude.statistics), sm[,!as.ergm_model(state[[1]])$etamap$offsetmap,drop=FALSE])) %>% lapply(mcmc, start=control.parallel$MCMC.burnin+1, thin=control.parallel$MCMC.interval) %>% as.mcmc.list() %>% window(., thin=thin(.)*max(1,floor(niter(.)/1000))) %>% plot(ask=FALSE,smooth=TRUE,density=FALSE)
     }
   }
 
diff --git a/R/ergm.llik.R b/R/ergm.llik.R
index 2feabae49..4a4b23977 100644
--- a/R/ergm.llik.R
+++ b/R/ergm.llik.R
@@ -32,7 +32,7 @@
 #                of the 'etamap'
 #   xsim       : the matrix of simulated statistics 
 #   probs      : the probability weight for each row of the stats matrix
-#   varweight  : the weight by which the variance of the base predictions will be 
+#   control.llik$MCMLE.varweight  : the weight by which the variance of the base predictions will be 
 #                scaled; the name of this param was changed from 'penalty' to better 
 #                reflect what this parameter actually is; default=0.5, which is the 
 #                "true"  weight, in the sense that the lognormal approximation is
@@ -60,9 +60,10 @@
 #####################################################################################
 
 llik.fun.lognormal <- function(theta, xsim, xsim.obs=NULL,
-                     varweight=0.5,
-                     dampening=FALSE,dampening.min.ess=100, dampening.level=0.1,
-                     eta0, etamap){
+                               eta0, etamap, 
+                               control.llik=control.logLik.ergm(),
+                               ...){
+
   # Convert theta to eta
   eta <- ergm.eta(theta, etamap)
 
@@ -72,7 +73,7 @@ llik.fun.lognormal <- function(theta, xsim, xsim.obs=NULL,
   basepred <- xsim %*% etaparam
   mb <- lweighted.mean(basepred,lrowweights(xsim))
   vb <- lweighted.var(basepred,lrowweights(xsim))
-  - mb - varweight*vb
+  - mb - control.llik$MCMLE.varweight*vb
 }
 
 
@@ -82,9 +83,9 @@ llik.fun.lognormal <- function(theta, xsim, xsim.obs=NULL,
 #####################################################################################
 
 llik.grad.IS <- function(theta, xsim,  xsim.obs=NULL,
-                     varweight=0.5,
-                     dampening=FALSE,dampening.min.ess=100, dampening.level=0.1,
-                     eta0, etamap){
+                         eta0, etamap, 
+                         control.llik=control.logLik.ergm(),
+                         ...){
 
   # Obtain canonical parameters incl. offsets and difference from sampled-from
   eta <- ergm.eta(theta, etamap)
@@ -114,9 +115,9 @@ llik.grad.IS <- function(theta, xsim,  xsim.obs=NULL,
 #       this is equation (3.5) of Hunter & Handcock (2006)
 #####################################################################################
 llik.hessian.IS <- function(theta, xsim, xsim.obs=NULL,
-                     varweight=0.5,
-                     dampening=FALSE,dampening.min.ess=100, dampening.level=0.1,
-                     eta0, etamap){
+                            eta0, etamap, 
+                            control.llik=control.logLik.ergm(),
+                            ...){
   # Obtain canonical parameters incl. offsets and difference from sampled-from
   eta <- ergm.eta(theta, etamap)
   etaparam <- eta-eta0
@@ -141,9 +142,9 @@ llik.hessian.IS <- function(theta, xsim, xsim.obs=NULL,
 #####################################################################################
 
 llik.fun.EF <- function(theta, xsim, xsim.obs=NULL,
-                     varweight=0.5,
-                     dampening=FALSE,dampening.min.ess=100, dampening.level=0.1,
-                     eta0, etamap){
+                        eta0, etamap, 
+                        control.llik=control.logLik.ergm(),
+                        ...){
   eta <- ergm.eta(theta, etamap)
   etaparam <- eta-eta0
   basepred <- xsim %*% etaparam
@@ -161,9 +162,9 @@ llik.fun.EF <- function(theta, xsim, xsim.obs=NULL,
 #####################################################################################
 
 llik.fun.IS <- function(theta, xsim, xsim.obs=NULL, 
-                     varweight=0.5,
-                     dampening=FALSE,dampening.min.ess=100, dampening.level=0.1,
-                     eta0, etamap){
+                        eta0, etamap, 
+                        control.llik=control.logLik.ergm(),
+                        ...){
   # Obtain canonical parameters incl. offsets and difference from sampled-from
   eta <- ergm.eta(theta, etamap)
   etaparam <- eta-eta0
@@ -179,9 +180,9 @@ llik.fun.IS <- function(theta, xsim, xsim.obs=NULL,
 #####################################################################################
 
 llik.fun.median <- function(theta, xsim, xsim.obs=NULL,
-                     varweight=0.5,
-                     dampening=FALSE,dampening.min.ess=100, dampening.level=0.1,
-                     eta0, etamap){
+                            eta0, etamap, 
+                            control.llik=control.logLik.ergm(),
+                            ...){
   # Convert theta to eta
   eta <- ergm.eta(theta, etamap)
 
@@ -201,27 +202,27 @@ llik.fun.median <- function(theta, xsim, xsim.obs=NULL,
 # 
 # This is the log-likelihood ratio (and not its negative)
 #
-  - (mb + varweight*sdb*sdb)
+  - (mb + control.llik$MCMLE.varweight*sdb*sdb)
 }
 
 llik.fun.logtaylor <- function(theta, xsim, xsim.obs=NULL, 
-		                     varweight=0.5,
-	 	                     dampening=FALSE,dampening.min.ess=100, dampening.level=0.1, 
-	 	                     eta0, etamap){ 
+                               eta0, etamap, 
+                               control.llik=control.logLik.ergm(),
+                               ...){
 	 	  # Convert theta to eta 
 	 	  eta <- ergm.eta(theta, etamap) 
  	 
  	  # Calculate approximation to l(eta) - l(eta0) using a lognormal approximation 
 	 	  etaparam <- eta-eta0 
 	 	  # 
-	 	  if (dampening) { 
+	 	  if (control.llik$MCMLE.dampening) { 
 	 	    #if theta_extended is (almost) outside convex hull don't trust theta 
-	 	    eta_extended <- eta + etaparam*dampening.level 
+	 	    eta_extended <- eta + etaparam*control.llik$MCMLE.dampening.level 
 	 	    expon_extended <- xsim %*% (eta_extended - eta0) 
 	 	    wts <- exp(expon_extended) 
 	 	    ess <- ceiling(sum(wts)^2/sum(wts^2)) 
 	 	#   https://xianblog.wordpress.com/2010/09/24/effective-sample-size/ 
-	 	    if(!is.na(ess) && {ess<dampening.min.ess}){ return(-Inf) } #.005*length(wts)) 
+	 	    if(!is.na(ess) && {ess<control.llik$MCMLE.dampening.min.ess}){ return(-Inf) } #.005*length(wts)) 
 	 	  } 
 	 	 
 	 	  basepred <- xsim %*% etaparam 
diff --git a/R/ergm.llik.obs.R b/R/ergm.llik.obs.R
index b915ee3d2..c9a2f73c3 100644
--- a/R/ergm.llik.obs.R
+++ b/R/ergm.llik.obs.R
@@ -33,7 +33,7 @@
 #   probs      : the probability weight for each row of the stats matrix
 #   xsim.obs  : the 'xsim' counterpart for observation process
 #   probs.obs : the 'probs' counterpart for observation process
-#   varweight  : the weight by which the variance of the base predictions will be
+#   control.llik$MCMLE.varweight  : the weight by which the variance of the base predictions will be
 #                scaled; the name of this param was changed from 'penalty' to better
 #                reflect what this parameter actually is; default=0.5, which is the
 #                "true"  weight, in the sense that the lognormal approximation is
@@ -58,9 +58,9 @@
 #        normally  distributed so that exp(eta * stats) is lognormal
 #####################################################################################                           
 llik.fun.obs.lognormal <- function(theta, xsim, xsim.obs=NULL,
-                     varweight=0.5,
-                     dampening=FALSE,dampening.min.ess=100, dampening.level=0.1,
-                     eta0, etamap){
+                                   eta0, etamap,
+                                   control.llik=control.logLik.ergm(),
+                                   ...){
   eta <- ergm.eta(theta, etamap)
   etaparam <- eta-eta0
 # These lines standardize:
@@ -78,7 +78,7 @@ llik.fun.obs.lognormal <- function(theta, xsim, xsim.obs=NULL,
 # 
 # This is the log-likelihood ratio (and not its negative)
 #
-  (mm + varweight*vm) - (mb + varweight*vb)
+  (mm + control.llik$MCMLE.varweight*vm) - (mb + control.llik$MCMLE.varweight*vb)
 }
 
 
@@ -88,9 +88,9 @@ llik.fun.obs.lognormal <- function(theta, xsim, xsim.obs=NULL,
 #   llg: the gradient of the not-offset eta parameters with ??
 #####################################################################################
 llik.grad.obs.IS <- function(theta, xsim,  xsim.obs=NULL,
-                      varweight=0.5,
-                      dampening=FALSE,dampening.min.ess=100, dampening.level=0.1,
-                      eta0, etamap){
+                             eta0, etamap,
+                             control.llik=control.logLik.ergm(),
+                             ...){
   # Obtain canonical parameters incl. offsets and difference from sampled-from
   eta <- ergm.eta(theta, etamap)
   etaparam <- eta-eta0
@@ -117,9 +117,9 @@ llik.grad.obs.IS <- function(theta, xsim,  xsim.obs=NULL,
 #####################################################################################
 
 llik.hessian.obs.IS <- function(theta, xsim, xsim.obs=NULL,
-                     varweight=0.5,
-                     dampening=FALSE,dampening.min.ess=100, dampening.level=0.1,
-                     eta0, etamap){
+                                eta0, etamap,
+                                control.llik=control.logLik.ergm(),
+                                ...){
   # Obtain canonical parameters incl. offsets and difference from sampled-from
   eta <- ergm.eta(theta, etamap)
   etaparam <- eta-eta0
@@ -149,9 +149,9 @@ llik.hessian.obs.IS <- function(theta, xsim, xsim.obs=NULL,
 #####################################################################################
 
 llik.fun.obs.IS <- function(theta, xsim, xsim.obs=NULL, 
-                     varweight=0.5,
-                     dampening=FALSE,dampening.min.ess=100, dampening.level=0.1,
-                     eta0, etamap){
+                            eta0, etamap,
+                            control.llik=control.logLik.ergm(),
+                            ...){
   # Obtain canonical parameters incl. offsets and difference from sampled-from
   eta <- ergm.eta(theta, etamap)
   etaparam <- eta-eta0
@@ -170,9 +170,9 @@ llik.fun.obs.IS <- function(theta, xsim, xsim.obs=NULL,
 #####################################################################################
 
 llik.fun.obs.robust<- function(theta, xsim, xsim.obs=NULL,
-                     varweight=0.5,
-                     dampening=FALSE,dampening.min.ess=100, dampening.level=0.1,
-                     eta0, etamap){
+                               eta0, etamap,
+                               control.llik=control.logLik.ergm(),
+                               ...){
   eta <- ergm.eta(theta, etamap)
   etaparam <- eta-eta0
 
@@ -191,7 +191,7 @@ llik.fun.obs.robust<- function(theta, xsim, xsim.obs=NULL,
 # 
 # This is the log-likelihood ratio (and not its negative)
 #
-  (mm + varweight*vm*vm) - (mb + varweight*vb*vb)
+  (mm + control.llik$MCMLE.varweight*vm*vm) - (mb + control.llik$MCMLE.varweight*vb*vb)
 }
 
 llik.fun.obs.robust<- llik.fun.obs.IS 
diff --git a/R/ergm.robmon.R b/R/ergm.robmon.R
index 40b4ff58c..f116f0304 100644
--- a/R/ergm.robmon.R
+++ b/R/ergm.robmon.R
@@ -54,6 +54,11 @@ ergm.robmon <- function(init, nw, model,
 
   control <- remap_algorithm_MCMC_controls(control, "RM")
 
+  control.llik <- list()
+  for(arg in STATIC_MCMLE_CONTROLS)
+    if(!is.null(control[[arg]]))
+      control.llik[arg] <- list(control[[arg]])
+
   #phase 1:  Estimate diagonal elements of D matrix (covariance matrix for init)
   n1 <- control$SA.phase1_n
   if(is.null(n1)) {n1 <- 7 + 3 * model$etamap$etalength} #default value
@@ -164,6 +169,7 @@ message(paste("theta new:",theta,""))
                    hessianflag=control$main.hessian,
                    method=control$MCMLE.method,
                    metric=control$MCMLE.metric,
+                   control.llik=control.llik,
                    verbose=verbose)
 
   ve$sample <- ergm.sample.tomcmc(ve$sample, control)
diff --git a/R/ergm.stocapprox.R b/R/ergm.stocapprox.R
index 880e1e1c2..a2ee352f1 100644
--- a/R/ergm.stocapprox.R
+++ b/R/ergm.stocapprox.R
@@ -57,6 +57,12 @@ ergm.stocapprox <- function(init, nw, model,
     stats <- model$nw.stats - NVL(model$target.stats,model$nw.stats)
     target.stats <- model$target.stats
   })
+
+  control.llik <- list(MCMLE.metric=control$MCMLE.metric, MCMLE.trustregion=control$SA.trustregion,
+                       MCMLE.varweight=control$MCMLE.varweight, MCMLE.dampening=control$MCMLE.dampening,
+                       MCMLE.dampening.min.ess=control$MCMLE.dampening.min.ess,
+                       MCMLE.dampening.level=control$MCMLE.dampening.level)
+
 # message(paste("Phase 1: ",n1,"iterations"))
 # message(paste(" (interval=",control$MCMC.interval,")",sep=""))
   nw.orig <- nw
@@ -127,6 +133,7 @@ ergm.stocapprox <- function(init, nw, model,
                    hessianflag=control$main.hessian,
                    method=control$MCMLE.method,
                    metric=control$MCMLE.metric,
+                   control.llik=control.llik,
                    verbose=verbose)
 #
 # Important: Keep R-M (pre-NR) theta
diff --git a/R/ergm_estfun.R b/R/ergm_estfun.R
index 7599c8162..6c4807535 100644
--- a/R/ergm_estfun.R
+++ b/R/ergm_estfun.R
@@ -30,36 +30,50 @@
 #'   \eqn{q}-vectors of estimating function values.
 #' @keywords internal
 #' @export
-ergm.estfun <- function(stats, theta, model, ...){
+ergm.estfun <- function(stats, theta, model, exclude=NULL, ...){
   UseMethod("ergm.estfun")
 }
 
 #' @describeIn ergm.estfun Method for numeric vectors of length \eqn{p}.
 #' @export
-ergm.estfun.numeric <- function(stats, theta, model, ...){
+ergm.estfun.numeric <- function(stats, theta, model, exclude=NULL, ...){
   etamap <- if(is(model, "ergm_model")) model$etamap else model
   estf <- c(ergm.etagradmult(theta,stats,etamap))[!etamap$offsettheta]
   names(estf) <- (if(is(model, "ergm_model")) param_names(model, FALSE) else names(theta))[!etamap$offsettheta]
+  # Exclude statistics that are not expected to have a specific mean
+  if(!is.null(exclude)){
+    x.exclude <- match(exclude,names(estf))
+    if(!is.na(x.exclude)){
+     estf <- estf[-x.exclude]
+    }
+  }
   -estf
 }
 
 #' @describeIn ergm.estfun Method for matrices with \eqn{p} columns.
 #' @export
-ergm.estfun.matrix <- function(stats, theta, model, ...){
+ergm.estfun.matrix <- function(stats, theta, model, exclude=NULL, ...){
   etamap <- if(is(model, "ergm_model")) model$etamap else model
   estf <- t(ergm.etagradmult(theta,t(as.matrix(stats)),etamap))[,!etamap$offsettheta,drop=FALSE]
   colnames(estf) <- (if(is(model, "ergm_model")) param_names(model, FALSE) else names(theta))[!etamap$offsettheta]
+# Exclude statistics that are not expected to have a specific mean
+  if(!is.null(exclude)){
+    x.exclude <- match(exclude,colnames(as.matrix(estf)))
+    if(!is.na(x.exclude)){
+     estf <- estf[,-x.exclude,drop=FALSE]
+    }
+  }
   -estf
 }
 
 #' @describeIn ergm.estfun Method for [`mcmc`] objects with \eqn{p} variables.
 #' @export
-ergm.estfun.mcmc <- function(stats, theta, model, ...){
-  mcmc(ergm.estfun(as.matrix(stats), theta, model), start=start(stats), end=end(stats), thin=thin(stats))
+ergm.estfun.mcmc <- function(stats, theta, model, exclude=NULL, ...){
+  mcmc(ergm.estfun(as.matrix(stats), theta, model, exclude=exclude), start=start(stats), end=end(stats), thin=thin(stats))
 }
 
 #' @describeIn ergm.estfun Method for  [`mcmc.list`] objects with \eqn{p} variables.
 #' @export
-ergm.estfun.mcmc.list <- function(stats, theta, model, ...){
-  lapply.mcmc.list(stats, ergm.estfun, theta, model)
+ergm.estfun.mcmc.list <- function(stats, theta, model, exclude=NULL, ...){
+  lapply.mcmc.list(stats, ergm.estfun, theta, model, exclude=exclude)
 }
diff --git a/R/logLik.ergm.R b/R/logLik.ergm.R
index e4d935281..cc99e2a16 100644
--- a/R/logLik.ergm.R
+++ b/R/logLik.ergm.R
@@ -88,7 +88,7 @@ logLik.ergm<-function(object, add=FALSE, force.reeval=FALSE, eval.loglik=add ||
 
   bridge.names <- names(formals(control.ergm.bridge))
   control.bridge <- control[intersect(bridge.names, names(control))]
-  control.transfer <- intersect(c("MCMC.samplesize", SCALABLE_MCMC_CONTROLS, STATIC_MCMC_CONTROLS, PARALLEL_MCMC_CONTROLS, MPLE_CONTROLS), bridge.names)
+  control.transfer <- intersect(c("MCMC.samplesize", SCALABLE_MCMC_CONTROLS, STATIC_MCMC_CONTROLS, PARALLEL_MCMC_CONTROLS, MPLE_CONTROLS, STATIC_MCMLE_CONTROLS), bridge.names)
   for(arg in control.transfer)
     if(is.null(control.bridge[[arg]]))
       control.bridge[[arg]] <- object$control[[arg]]
diff --git a/R/parallel.utils.R b/R/parallel.utils.R
index 9e410f38c..aeaf3fb75 100644
--- a/R/parallel.utils.R
+++ b/R/parallel.utils.R
@@ -86,7 +86,7 @@ ergm.MCMC.packagenames <- local({
 #' multiple threads run in parallel. This is done using the
 #' [OpenMP](https://www.openmp.org/) API.
 #'
-#' However, this introduces a nontrivial amont of computational
+#' However, this introduces a nontrivial amount of computational
 #' overhead. See below for a list of the major factors affecting
 #' whether it is worthwhile.}}
 #'
diff --git a/R/summary.ergm.R b/R/summary.ergm.R
index 42c8f3092..9e3f1f7ef 100644
--- a/R/summary.ergm.R
+++ b/R/summary.ergm.R
@@ -87,7 +87,7 @@ summary.ergm <- function (object, ...,
     warn(paste0("This object was fit with ", sQuote("ergm"), " version ", objver, " or earlier. Summarizing it with version ", nextver, " or later may return incorrect results or fail."))
   }
 
-  if("digits" %in% names(list(...))) warn("summary.ergm() no lnger takes a digits= argument.")
+  if("digits" %in% names(list(...))) warn("summary.ergm() no longer takes a digits= argument.")
   control <- object$control
   pseudolikelihood <- object$estimate=="MPLE"
   independence <- NVL(object$MPLE_is_MLE, is.dyad.independent(object))
diff --git a/R/zzz.R b/R/zzz.R
index dd1eb825c..464b4dcaf 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -17,7 +17,7 @@
 }
 
 .onLoad <- function(libname, pkgname){
-  # . is used as a placeholder by stantet.common::NVL3().
+  # . is used as a placeholder by statnet.common::NVL3().
   utils::globalVariables(".")
 
   default_options(ergm.eval.loglik=TRUE,
diff --git a/man/control.ergm.Rd b/man/control.ergm.Rd
index 49652bebb..e5bcadcc1 100644
--- a/man/control.ergm.Rd
+++ b/man/control.ergm.Rd
@@ -38,6 +38,7 @@ control.ergm(
   MCMC.maxedges = Inf,
   MCMC.addto.se = TRUE,
   MCMC.packagenames = c(),
+  MCMC.esteq.exclude.statistics = NULL,
   SAN.maxit = 4,
   SAN.nsteps.times = 8,
   SAN = control.san(term.options = term.options, SAN.maxit = SAN.maxit,
@@ -69,11 +70,12 @@ control.ergm(
   obs.MCMC.impute.default_density = function(nw) 2/network.size(nw),
   MCMLE.min.depfac = 2,
   MCMLE.sampsize.boost.pow = 0.5,
+  MCMLE.varweight = 0.5,
   MCMLE.MCMC.precision = if (startsWith("confidence", MCMLE.termination[1])) 0.1 else
     0.005,
   MCMLE.MCMC.max.ESS.frac = 0.1,
   MCMLE.metric = c("lognormal", "logtaylor", "Median.Likelihood", "EF.Likelihood",
-    "naive"),
+    "Kpenalty", "naive"),
   MCMLE.method = c("BFGS", "Nelder-Mead"),
   MCMLE.dampening = FALSE,
   MCMLE.dampening.min.ess = 20,
@@ -371,6 +373,9 @@ algorithm to the estimates' standard errors.}
 statistic functions in addition to those autodetected. This argument should
 not be needed outside of very strange setups.}
 
+\item{MCMC.esteq.exclude.statistics}{character string pattern. statistics with names that contain this string will be excluded from
+the estimating equations. For example, this supports tapering models which may have terms with non-specific means.}
+
 \item{SAN.maxit}{When \code{target.stats} argument is passed to
 \code{\link[=ergm]{ergm()}}, the maximum number of attempts to use \code{\link{san}}
 to obtain a network with statistics close to those specified.}
@@ -472,6 +477,9 @@ imputation density when number of non-missing dyads is too low.}
 
 \item{MCMLE.min.depfac, MCMLE.sampsize.boost.pow}{When using adaptive MCMC effective size, and methods that increase the MCMC sample size, use \code{MCMLE.sampsize.boost.pow} as the power of the boost amount (relative to the boost of the target effective size), but ensure that sample size is no less than \code{MCMLE.min.depfac} times the target effective size.}
 
+\item{MCMLE.varweight}{numeric multiplier for covariance term in loglikelihood
+where 0.5 is the "true" value but this can be increased or decreased.}
+
 \item{MCMLE.MCMC.precision, MCMLE.MCMC.max.ESS.frac}{\code{MCMLE.MCMC.precision} is a vector of upper bounds on the standard
 errors induced by the MCMC algorithm, expressed as a percentage of the total
 standard error. The MCMLE algorithm will terminate when the MCMC standard
@@ -482,7 +490,9 @@ If effective sample size is used (see \code{MCMC.effectiveSize}), then ergm
 may increase the target ESS to reduce the MCMC standard error.}
 
 \item{MCMLE.metric}{Method to calculate the loglikelihood approximation.
-See Hummel et al (2010) for an explanation of "lognormal" and "naive".}
+See Hummel et al (2010) for an explanation of "lognormal" and "naive".
+"Kpenalty" is used to specify the use of a tapering model with a penalty for
+kurtosis. It is used by the \code{ergm.tapered} package and is not usually user specified.}
 
 \item{MCMLE.method}{Deprecated. By default, ergm uses \code{trust}, and
 falls back to \code{optim} with Nelder-Mead method when \code{trust} fails.}
diff --git a/man/control.simulate.ergm.Rd b/man/control.simulate.ergm.Rd
index 4a708dc89..b7d8fbf6c 100644
--- a/man/control.simulate.ergm.Rd
+++ b/man/control.simulate.ergm.Rd
@@ -21,6 +21,7 @@ control.simulate.formula.ergm(
   MCMC.effectiveSize.order.max = NULL,
   MCMC.maxedges = Inf,
   MCMC.packagenames = c(),
+  MCMC.esteq.exclude.statistics = NULL,
   MCMC.runtime.traceplot = FALSE,
   network.output = "network",
   term.options = NULL,
@@ -45,6 +46,7 @@ control.simulate(
   MCMC.effectiveSize.order.max = NULL,
   MCMC.maxedges = Inf,
   MCMC.packagenames = c(),
+  MCMC.esteq.exclude.statistics = NULL,
   MCMC.runtime.traceplot = FALSE,
   network.output = "network",
   term.options = NULL,
@@ -69,6 +71,7 @@ control.simulate.formula(
   MCMC.effectiveSize.order.max = NULL,
   MCMC.maxedges = Inf,
   MCMC.packagenames = c(),
+  MCMC.esteq.exclude.statistics = NULL,
   MCMC.runtime.traceplot = FALSE,
   network.output = "network",
   term.options = NULL,
@@ -94,6 +97,7 @@ control.simulate.ergm(
   MCMC.effectiveSize.order.max = NULL,
   MCMC.maxedges = Inf,
   MCMC.packagenames = NULL,
+  MCMC.esteq.exclude.statistics = NULL,
   MCMC.runtime.traceplot = FALSE,
   network.output = "network",
   term.options = NULL,
@@ -166,6 +170,9 @@ sample size and the variance for the Geweke diagnostic.}
 statistic functions in addition to those autodetected. This argument should
 not be needed outside of very strange setups.}
 
+\item{MCMC.esteq.exclude.statistics}{character string pattern. statistics with names that contain this string will be excluded from
+the estimating equations. For example, this supports tapering models which may have terms with non-specific means.}
+
 \item{MCMC.runtime.traceplot}{Logical: If \code{TRUE}, plot traceplots of the MCMC
 sample.}
 
@@ -212,7 +219,7 @@ function.
 While the others supply a full set of simulation
 settings, \code{control.simulate.ergm} when passed as a control
 parameter to \code{\link[=simulate.ergm]{simulate.ergm()}} allows some settings to be
-inherited from the ERGM stimation while overriding others.
+inherited from the ERGM estimation while overriding others.
 }
 \details{
 This function is only used within a call to the \code{\link{simulate}}
diff --git a/man/ergm-terms.Rd b/man/ergm-terms.Rd
index 2cb4d631b..be6e6b5b1 100644
--- a/man/ergm-terms.Rd
+++ b/man/ergm-terms.Rd
@@ -506,7 +506,7 @@
       taking \eqn{y} and removing any edges for which
       \eqn{f_{i,j}(y_{i,j}) = 0}{f[i,j] (y[i,j])=0}.
 
-      For convenience, the term in \code{filter} can be a part of a simple logical or comparison operation: (e.g., \code{~!nodematch("A")} or \code{~absdiff("X")>3}), which filters on \eqn{f_{i,j}(y_{i,j}) \square 0}{f[i,j] (y[i,j]) \%OP\% 0} instead.
+      For convenience, the term in \code{filter} can be a part of a simple logical or comparison operation: (e.g., \code{~!nodematch("A")} or \code{~absdiff("X")>3}), which filters on \eqn{f_{i,j}(y_{i,j}) \fbox{$\phantom{5}$} 0}{f[i,j] (y[i,j]) \%OP\% 0} instead.
     }
 
     \item{\code{Log(formula, log0=-1/sqrt(.Machine$double.eps))} (binary) (operator), \code{Log(formula, log0=-1/sqrt(.Machine$double.eps))} (valued) (operator)}{\emph{Take a natural logarithm of a network's statistic}:
diff --git a/man/ergm.estfun.Rd b/man/ergm.estfun.Rd
index 624e6d6d0..d5ab901e2 100644
--- a/man/ergm.estfun.Rd
+++ b/man/ergm.estfun.Rd
@@ -8,15 +8,15 @@
 \alias{ergm.estfun.mcmc.list}
 \title{Compute the Sample Estimating Function Values of an ERGM.}
 \usage{
-ergm.estfun(stats, theta, model, ...)
+ergm.estfun(stats, theta, model, exclude = NULL, ...)
 
-\method{ergm.estfun}{numeric}(stats, theta, model, ...)
+\method{ergm.estfun}{numeric}(stats, theta, model, exclude = NULL, ...)
 
-\method{ergm.estfun}{matrix}(stats, theta, model, ...)
+\method{ergm.estfun}{matrix}(stats, theta, model, exclude = NULL, ...)
 
-\method{ergm.estfun}{mcmc}(stats, theta, model, ...)
+\method{ergm.estfun}{mcmc}(stats, theta, model, exclude = NULL, ...)
 
-\method{ergm.estfun}{mcmc.list}(stats, theta, model, ...)
+\method{ergm.estfun}{mcmc.list}(stats, theta, model, exclude = NULL, ...)
 }
 \arguments{
 \item{stats}{An object representing sample statistics with observed values subtracted out.}

From 2dc98f431bb8a5301933edb6d0ceeef71359528b Mon Sep 17 00:00:00 2001
From: "Mark S. Handcock" <handcock@ucla.edu>
Date: Thu, 1 Jul 2021 11:51:21 -0700
Subject: [PATCH 2/8] Added in exclusion for "Taper_Penalty" term from
 mcmc.diagnostics.

---
 R/mcmc.diagnostics.ergm.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/mcmc.diagnostics.ergm.R b/R/mcmc.diagnostics.ergm.R
index b8c0d6cdb..03582fd97 100644
--- a/R/mcmc.diagnostics.ergm.R
+++ b/R/mcmc.diagnostics.ergm.R
@@ -207,8 +207,8 @@ mcmc.diagnostics.ergm <- function(object,
 
   if(esteq){
     if (!is.null(coef(object)) && !is.null(object$etamap)) {
-      sm <- ergm.estfun(sm, theta=coef(object), model=object$etamap) %>% lapply.mcmc.list(`-`)
-      if(!is.null(sm.obs)) sm.obs <- ergm.estfun(sm.obs, theta=coef(object), model=object$etamap) %>% lapply.mcmc.list(`-`)
+      sm <- ergm.estfun(sm, theta=coef(object), model=object$etamap, exclude=object$control$MCMC.esteq.exclude.statistics) %>% lapply.mcmc.list(`-`)
+      if(!is.null(sm.obs)) sm.obs <- ergm.estfun(sm.obs, theta=coef(object), model=object$etamap, exclude=object$control$MCMC.esteq.exclude.statistics) %>% lapply.mcmc.list(`-`)
     }
   }
 

From e709ad3dec3b2af8c7cc46a8dd948f7c2cdc69dd Mon Sep 17 00:00:00 2001
From: "Mark S. Handcock" <handcock@ucla.edu>
Date: Fri, 16 Jul 2021 17:20:16 -0700
Subject: [PATCH 3/8] Extend the exclude code and robustify it.

---
 R/ergm.MCMCse.R | 14 +++++++++++---
 R/ergm.bridge.R | 11 +++++++++--
 R/ergm_estfun.R | 21 ++++++++++++++++++++-
 3 files changed, 40 insertions(+), 6 deletions(-)

diff --git a/R/ergm.MCMCse.R b/R/ergm.MCMCse.R
index fc7d82d00..4deaa8684 100644
--- a/R/ergm.MCMCse.R
+++ b/R/ergm.MCMCse.R
@@ -24,7 +24,7 @@
 #'   latter are always 1 for lognormal metric.
 #' @noRd
 ergm.MCMCse <- function(model, theta, init, statsmatrices, statsmatrices.obs,
-                        H, H.obs, metric = c("IS", "lognormal")) {
+                        H, H.obs, metric = c("IS", "lognormal"), exclude=NULL) {
   metric <- match.arg(metric)
 
   # Transform theta to eta
@@ -39,7 +39,7 @@ ergm.MCMCse <- function(model, theta, init, statsmatrices, statsmatrices.obs,
   av <- colMeans.mcmc.list(statsmatrices)
 # av <- apply(statsmatrices,2,median)
   xsims <- sweep.mcmc.list(statsmatrices, av, "-")
-  gsims <- ergm.estfun(xsims, theta=theta, model=model)
+  gsims <- ergm.estfun(xsims, theta=theta, model=model, exclude=exclude)
   xobs <- -av
   xsims <- xsims[,!offsetmap, drop=FALSE]
   xsim <- as.matrix(xsims)
@@ -49,7 +49,7 @@ ergm.MCMCse <- function(model, theta, init, statsmatrices, statsmatrices.obs,
    av.obs <- colMeans.mcmc.list(statsmatrices.obs)
 #  av.obs <- apply(statsmatrices.obs, 2, median)
    xsims.obs <- sweep.mcmc.list(statsmatrices.obs, av.obs,"-")
-   gsims.obs <- ergm.estfun(xsims.obs, theta=theta, model=model)
+   gsims.obs <- ergm.estfun(xsims.obs, theta=theta, model=model, exclude=exclude)
    xsims.obs <- xsims.obs[,!offsetmap, drop=FALSE]
    xsim.obs <- as.matrix(xsims.obs)
    gsim.obs <- as.matrix(gsims.obs)
@@ -127,6 +127,14 @@ ergm.MCMCse <- function(model, theta, init, statsmatrices, statsmatrices.obs,
     mc.cov[!novar,!novar] <- mc.cov0
   }
 
+  if(!is.null(exclude)){
+    x.exclude <- match(exclude,names(theta))
+    if(!is.na(x.exclude)){
+     offsettheta[x.exclude] <- TRUE
+    }
+  }
+
+
   mc.cov.offset[!offsettheta,!offsettheta] <- mc.cov
 
   rownames(mc.cov.offset) <- colnames(mc.cov.offset) <- param_names(model)
diff --git a/R/ergm.bridge.R b/R/ergm.bridge.R
index 6c657df2c..410a497c3 100644
--- a/R/ergm.bridge.R
+++ b/R/ergm.bridge.R
@@ -128,6 +128,10 @@ ergm.bridge.llr<-function(object, response=NULL, reference=~Bernoulli, constrain
 
   ## Miscellaneous settings
   Dtheta.Du <- (to-from)[!state[[1]]$model$etamap$offsettheta] / control$bridge.nsteps
+  x.exclude <- match("Taper_Penalty",names(Dtheta.Du))
+  if(!is.na(x.exclude)){
+    Dtheta.Du <- Dtheta.Du[-x.exclude]
+  }
 
   ## Handle target statistics, if passed.
   if(!is.null(target.stats)){
@@ -141,8 +145,11 @@ ergm.bridge.llr<-function(object, response=NULL, reference=~Bernoulli, constrain
 
   ## Helper function to calculate Dtheta.Du %*% Deta.Dtheta %*% g(y)
   llrsamp <- function(samp, theta){
-    if(is.mcmc.list(samp)) lapply.mcmc.list(ergm.estfun(samp, theta, state[[1]]$model$etamap), `%*%`, Dtheta.Du)
-    else sum(ergm.estfun(samp, theta, state[[1]]$model$etamap) * Dtheta.Du)
+    if(is.mcmc.list(samp)){
+      lapply.mcmc.list(ergm.estfun(samp, theta, state[[1]]$model$etamap, exclude="Taper_Penalty"), `%*%`, Dtheta.Du)
+    }else{
+      sum(ergm.estfun(samp, theta, state[[1]]$model$etamap, exclude="Taper_Penalty") * Dtheta.Du)
+    }
   }
 
 
diff --git a/R/ergm_estfun.R b/R/ergm_estfun.R
index 6c4807535..641606d04 100644
--- a/R/ergm_estfun.R
+++ b/R/ergm_estfun.R
@@ -40,6 +40,16 @@ ergm.estfun.numeric <- function(stats, theta, model, exclude=NULL, ...){
   etamap <- if(is(model, "ergm_model")) model$etamap else model
   estf <- c(ergm.etagradmult(theta,stats,etamap))[!etamap$offsettheta]
   names(estf) <- (if(is(model, "ergm_model")) param_names(model, FALSE) else names(theta))[!etamap$offsettheta]
+  if(is(model, "ergm_model")){
+   names(estf) <- param_names(model, FALSE)[!etamap$offsettheta]
+  }else{
+    if(!is.null(names(theta))){
+      names(estf) <- names(theta)[!etamap$offsettheta]
+    }else{
+      a <- names(stats)[!etamap$offsettheta]
+      if(length(a)==length(estf)) names(estf) <- a
+    }
+  }
   # Exclude statistics that are not expected to have a specific mean
   if(!is.null(exclude)){
     x.exclude <- match(exclude,names(estf))
@@ -55,7 +65,16 @@ ergm.estfun.numeric <- function(stats, theta, model, exclude=NULL, ...){
 ergm.estfun.matrix <- function(stats, theta, model, exclude=NULL, ...){
   etamap <- if(is(model, "ergm_model")) model$etamap else model
   estf <- t(ergm.etagradmult(theta,t(as.matrix(stats)),etamap))[,!etamap$offsettheta,drop=FALSE]
-  colnames(estf) <- (if(is(model, "ergm_model")) param_names(model, FALSE) else names(theta))[!etamap$offsettheta]
+  if(is(model, "ergm_model")){
+   colnames(estf) <- param_names(model, FALSE)[!etamap$offsettheta]
+  }else{
+    if(!is.null(names(theta))){
+      colnames(estf) <- names(theta)[!etamap$offsettheta]
+    }else{
+      a <- colnames(stats)[!etamap$offsettheta]
+      if(length(a)==ncol(estf)) colnames(estf) <- a
+    }
+  }
 # Exclude statistics that are not expected to have a specific mean
   if(!is.null(exclude)){
     x.exclude <- match(exclude,colnames(as.matrix(estf)))

From da7e4c1bef5034f586ed8fe956e5825ba23b075a Mon Sep 17 00:00:00 2001
From: "Mark S. Handcock" <handcock@ucla.edu>
Date: Fri, 16 Jul 2021 17:22:49 -0700
Subject: [PATCH 4/8] Added direct deterministic gradients and hessian for
 Kpenalty.

These required extensive algebra and computation. However, they dramatically improve the computational speed and accuracy of the Kpenalty method.
---
 R/ergm.estimate.R | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/R/ergm.estimate.R b/R/ergm.estimate.R
index df626d9f1..9ea9c76fd 100644
--- a/R/ergm.estimate.R
+++ b/R/ergm.estimate.R
@@ -186,7 +186,7 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
                          logtaylor=llik.grad.IS,
                          Median.Likelihood=llik.grad.IS,
                          EF.Likelihood=llik.grad.IS,
-                         Kpenalty=ergm.tapered::llik.grad.Kpenalty.numDeriv,
+                         Kpenalty=ergm.tapered::llik.grad.Kpenalty,
                          llik.grad.IS)
     Hessianfn <- switch(metric,
                         Likelihood=llik.hessian.IS,
@@ -194,7 +194,7 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
                         logtaylor=llik.hessian.IS,
                         Median.Likelihood=llik.hessian.IS,
                         EF.Likelihood=llik.hessian.IS,
-                        Kpenalty=ergm.tapered::llik.hessian.Kpenalty.numDeriv,
+                        Kpenalty=ergm.tapered::llik.hessian.Kpenalty,
                         llik.hessian.IS)
   }
   
@@ -330,7 +330,7 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
     mc.cov <- matrix(NA, length(theta), length(theta))
     covar <- NA
     if(!hessianflag || steplen!=1){
-      Lout$hessian <- Hessianfn(theta=Lout$par,
+      Lout$hessian <- llik.hessian.IS(theta=Lout$par,
                         xsim=xsim.orig,
                         xsim.obs=xsim.orig.obs,
                         eta0=eta0, etamap=etamap.no,

From 236c9ba7621c7085b2fda9fd3d53a561343a531b Mon Sep 17 00:00:00 2001
From: handcock <handcock@ucla.edu>
Date: Fri, 26 Aug 2022 17:34:59 -0700
Subject: [PATCH 5/8] Added a MPLE.save.xmat control argument

Added a MPLE.save.xmat control argument to return the design matrix from the MPLE fit (if any). The current and retained default is to remove it.
---
 R/control.ergm.R     | 4 ++++
 R/ergm.R             | 2 +-
 R/ergm.mple.R        | 2 ++
 man/control.ergm.Rd  | 4 ++++
 man/ergm-parallel.Rd | 2 +-
 5 files changed, 12 insertions(+), 2 deletions(-)

diff --git a/R/control.ergm.R b/R/control.ergm.R
index 6f1f6dc42..2e473ca71 100644
--- a/R/control.ergm.R
+++ b/R/control.ergm.R
@@ -128,6 +128,9 @@
 #'   necessary. Note that this can be very dangerous unless you know
 #'   what you are doing.
 #'
+#' @param MPLE.save.xmat If `TRUE`, the design matrix of change statistics
+#'   from the pseudo-likelihood computation will be returned.
+#'
 #' @template control_MCMC_prop
 #'
 #' @param MCMC.interval Number of proposals between sampled statistics.
@@ -484,6 +487,7 @@ control.ergm<-function(drop=TRUE,
                        MPLE.nonident=c("warning","message","error"),
                        MPLE.nonident.tol=1e-10,
                        MPLE.constraints.ignore=FALSE,
+                       MPLE.save.xmat=FALSE,
 
                        MCMC.prop=trim_env(~sparse),
                        MCMC.prop.weights="default", MCMC.prop.args=list(),
diff --git a/R/ergm.R b/R/ergm.R
index 64d37617c..33f79ce3e 100644
--- a/R/ergm.R
+++ b/R/ergm.R
@@ -754,7 +754,7 @@ ergm <- function(formula, response=NULL,
                                                                nonvar_action = control$MPLE.nonvar))
          )
 
-  initialfit$xmat.full <- NULL # No longer needed but takes up space.
+  if(!control$MPLE.save.xmat) initialfit$xmat.full <- NULL # No longer needed but takes up space.
 
   estimate.desc <- switch(estimate,
                           MPLE = if(MPLE.is.MLE) "Maximum Likelihood"
diff --git a/R/ergm.mple.R b/R/ergm.mple.R
index 59a29f696..255f85f32 100644
--- a/R/ergm.mple.R
+++ b/R/ergm.mple.R
@@ -95,6 +95,7 @@ ergm.mple<-function(nw, fd, m, init=NULL,
    mplefit.summary <- mplefit
   }else{
    if(MPLEtype=="logitreg"){
+    glm.result <- NULL
     mplefit <- model.matrix(terms(pl$zy ~ .-1,data=data.frame(pl$xmat)),
                            data=data.frame(pl$xmat))
     mplefit <- ergm.logitreg(x=mplefit, y=pl$zy, m=m, wt=pl$wend,
@@ -197,6 +198,7 @@ ergm.mple<-function(nw, fd, m, init=NULL,
       mple.lik.null = structure(
         ERRVL(try(logLik(mplefit.null), silent=TRUE), -mplefit.null$deviance/2),
         nobs = nobs, df = df, class="logLik"),
+      glm.result = if(save.xmat) glm.result,
       xmat.full = if(save.xmat) pl$xmat.full
       ),
       class="ergm")
diff --git a/man/control.ergm.Rd b/man/control.ergm.Rd
index 0be36a73b..26fe689b2 100644
--- a/man/control.ergm.Rd
+++ b/man/control.ergm.Rd
@@ -21,6 +21,7 @@ control.ergm(
   MPLE.nonident = c("warning", "message", "error"),
   MPLE.nonident.tol = 1e-10,
   MPLE.constraints.ignore = FALSE,
+  MPLE.save.xmat = FALSE,
   MCMC.prop = trim_env(~sparse),
   MCMC.prop.weights = "default",
   MCMC.prop.args = list(),
@@ -296,6 +297,9 @@ missingness. This can be used to avert evaluating and storing the
 necessary. Note that this can be very dangerous unless you know
 what you are doing.}
 
+\item{MPLE.save.xmat}{If \code{TRUE}, the design matrix of change statistics
+from the pseudo-likelihood computation will be returned.}
+
 \item{MCMC.prop}{Specifies the proposal (directly) and/or
 a series of "hints" about the structure of the model being
 sampled. The specification is in the form of a one-sided formula
diff --git a/man/ergm-parallel.Rd b/man/ergm-parallel.Rd
index cbaef6d46..03f758bd0 100644
--- a/man/ergm-parallel.Rd
+++ b/man/ergm-parallel.Rd
@@ -114,7 +114,7 @@ sampling by evaluating qualified terms' change statistics in
 multiple threads run in parallel. This is done using the
 \href{https://www.openmp.org/}{OpenMP} API.
 
-However, this introduces a nontrivial amont of computational
+However, this introduces a nontrivial amount of computational
 overhead. See below for a list of the major factors affecting
 whether it is worthwhile.}}
 

From fdd727d88ed71c6f6407c9ef9d0f9d4a0dc4545a Mon Sep 17 00:00:00 2001
From: handcock <handcock@ucla.edu>
Date: Fri, 26 Aug 2022 17:39:51 -0700
Subject: [PATCH 6/8] Added a MCMC.esteq.exclude.statistics control argument

control$MCMC.esteq.exclude.statistics allows network statistics to not be included in the estimating equations. There equations determine most convergence diagnostics. For some models there are some terms that are not zero at convergence and should be removed.  For example, this supports tapering models which may have terms with non-specific means. The estimating equations are used in a number of places and the code follows those places.

Specifically, it is a vector of  strings being the names of the statistics to be excluded from
the estimating equations.
---
 R/control.ergm.bridge.R    |  3 +++
 R/ergm.MCMLE.R             |  9 +++++++--
 R/ergm.bridge.R            | 12 +++++++-----
 R/ergm.estimate.R          |  3 +--
 R/ergm.getMCMCsample.R     | 10 ++++++----
 man/control.ergm.bridge.Rd |  6 +++++-
 6 files changed, 29 insertions(+), 14 deletions(-)

diff --git a/R/control.ergm.bridge.R b/R/control.ergm.bridge.R
index 6e5d65102..9e5ece980 100644
--- a/R/control.ergm.bridge.R
+++ b/R/control.ergm.bridge.R
@@ -41,6 +41,8 @@
 #' @template control_MCMC_parallel
 #' @template seed
 #' @template control_MCMC_packagenames
+#' @param MCMC.esteq.exclude.statistics character string pattern. statistics with names that contain this string will be excluded from 
+#' the estimating equations. For example, this supports tapering models which may have terms with non-specific means.
 #' @template control_dots
 #' @return A list with arguments as components.
 #' @seealso [ergm.bridge.llr()]
@@ -69,6 +71,7 @@ control.ergm.bridge<-function(bridge.nsteps=16, # Number of geometric bridges to
 
                               MCMC.maxedges=Inf,
                               MCMC.packagenames=c(),
+                              MCMC.esteq.exclude.statistics=NULL,
                               
                               term.options=list(),
                               seed=NULL,
diff --git a/R/ergm.MCMLE.R b/R/ergm.MCMLE.R
index 98b697532..27acb00bf 100644
--- a/R/ergm.MCMLE.R
+++ b/R/ergm.MCMLE.R
@@ -297,7 +297,12 @@ ergm.MCMLE <- function(init, nw, model,
     }
 
     # Compute the sample estimating equations and the convergence p-value. 
-    esteqs <- ergm.estfun(statsmatrices, theta=mcmc.init, model=model, exclude=control$MCMC.esteq.exclude.statistics)
+    #esteqs <- ergm.estfun(statsmatrices, theta=mcmc.init, model=model, exclude=control$MCMC.esteq.exclude.statistics)
+    esteqs <- ergm.estfun(statsmatrices, theta=mcmc.init, model=model)
+    if(!is.null(control$MCMC.esteq.exclude.statistics)){
+      x.exclude <- match(control$MCMC.esteq.exclude.statistics,colnames(esteqs[[1]]))
+      if(!is.na(x.exclude)){ esteqs[[1]] <- esteqs[[1]][,-x.exclude,drop=FALSE] }
+    }
     esteq <- as.matrix(esteqs)
     if(isTRUE(all.equal(apply(esteq,2,stats::sd), rep(0,ncol(esteq)), check.names=FALSE))&&!all(esteq==0))
       stop("Unconstrained MCMC sampling did not mix at all. Optimization cannot continue.")
@@ -552,7 +557,7 @@ ergm.MCMLE <- function(init, nw, model,
           control$MCMC.effectiveSize <- round(control$MCMC.effectiveSize * prec.scl)
           if(control$MCMC.effectiveSize/control$MCMC.samplesize>control$MCMLE.MCMC.max.ESS.frac) control$MCMC.samplesize <- control$MCMC.effectiveSize/control$MCMLE.MCMC.max.ESS.frac
           # control$MCMC.samplesize <- round(control$MCMC.samplesize * prec.scl)
-          message("Increasing target MCMC sample size to ", control$MCMC.samplesize, ", ESS to",control$MCMC.effectiveSize,".")
+          message("Increasing target MCMC sample size to ", control$MCMC.samplesize, ", ESS to ",control$MCMC.effectiveSize,".")
         } else { # Fixed-interval sampling
           control$MCMC.samplesize <- round(control$MCMC.samplesize * prec.scl)
           control$MCMC.burnin <- round(control$MCMC.burnin * prec.scl)
diff --git a/R/ergm.bridge.R b/R/ergm.bridge.R
index 7ef172e36..eb120862d 100644
--- a/R/ergm.bridge.R
+++ b/R/ergm.bridge.R
@@ -154,9 +154,11 @@ ergm.bridge.llr<-function(object, response=NULL, reference=~Bernoulli, constrain
 
   ## Miscellaneous settings
   Dtheta.Du <- (to-from)[!state[[1]]$model$etamap$offsettheta]
-  x.exclude <- match("Taper_Penalty",names(Dtheta.Du))
-  if(!is.na(x.exclude)){
-    Dtheta.Du <- Dtheta.Du[-x.exclude]
+  if(!is.null(control$MCMC.esteq.exclude.statistics)){
+    x.exclude <- match(control$MCMC.esteq.exclude.statistics,names(Dtheta.Du))
+    if(!is.na(x.exclude)){
+      Dtheta.Du <- Dtheta.Du[-x.exclude]
+    }
   }
 
   ## Handle target statistics, if passed.
@@ -172,9 +174,9 @@ ergm.bridge.llr<-function(object, response=NULL, reference=~Bernoulli, constrain
   ## Helper function to calculate Dtheta.Du %*% Deta.Dtheta %*% g(y)
   llrsamp <- function(samp, theta){
     if(is.mcmc.list(samp)){
-      lapply.mcmc.list(ergm.estfun(samp, theta, state[[1]]$model$etamap, exclude="Taper_Penalty"), `%*%`, Dtheta.Du)
+      lapply.mcmc.list(ergm.estfun(samp, theta, state[[1]]$model$etamap, exclude=control$MCMC.esteq.exclude.statistics), `%*%`, Dtheta.Du)
     }else{
-      sum(ergm.estfun(samp, theta, state[[1]]$model$etamap, exclude="Taper_Penalty") * Dtheta.Du)
+      sum(ergm.estfun(samp, theta, state[[1]]$model$etamap, exclude=control$MCMC.esteq.exclude.statistics) * Dtheta.Du)
     }
   }
 
diff --git a/R/ergm.estimate.R b/R/ergm.estimate.R
index 51ff8d99d..24959cf4b 100644
--- a/R/ergm.estimate.R
+++ b/R/ergm.estimate.R
@@ -194,7 +194,7 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
                         logtaylor=llik.hessian.IS,
                         Median.Likelihood=llik.hessian.IS,
                         EF.Likelihood=llik.hessian.IS,
-                        Kpenalty=ergm.tapered::llik.hessian.Kpenalty,
+                        Kpenalty=llik.hessian.IS,
                         llik.hessian.IS)
   }
   
@@ -276,7 +276,6 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
                       eta0=eta0, etamap=etamap.no,
                       control.llik=control.llik),
             silent=FALSE)
-    Lout$par<-Lout$argument
 
     if(inherits(Lout,"try-error")) {
       message("MLE could not be found. Trying Nelder-Mead...")
diff --git a/R/ergm.getMCMCsample.R b/R/ergm.getMCMCsample.R
index 77998c09b..a4a4556a9 100644
--- a/R/ergm.getMCMCsample.R
+++ b/R/ergm.getMCMCsample.R
@@ -239,10 +239,12 @@ ergm_MCMC_sample <- function(state, control, theta=NULL,
         if(verbose) message("Increasing thinning to ",interval,".")
       }
       
-<
-      esteq <- lapply(sms, function(sm) NVL3(theta, ergm.estfun(sm, ., as.ergm_model(state[[1]]), exclude=control$MCMC.esteq.exclude.statistics), sm[,!as.ergm_model(state[[1]])$etamap$offsetmap,drop=FALSE])) %>%
+      esteq <- lapply(sms, function(sm) NVL3(theta, ergm.estfun(sm, ., as.ergm_model(state[[1]]),
+        exclude=control$MCMC.esteq.exclude.statistics), sm[,!as.ergm_model(state[[1]])$etamap$offsetmap,drop=FALSE])) %>%
         lapply.mcmc.list(mcmc, start=1, thin=interval) %>% lapply.mcmc.list(`-`)
 
+      if("Taper_Penalty" %in% dimnames(esteq[[1]])[[2]]) {browser()}
+
       if(control.parallel$MCMC.runtime.traceplot){
         plot(window(esteq, thin=thin(esteq)*max(1,floor(niter(esteq)/1000)))
             ,ask=FALSE,smooth=TRUE,density=FALSE)
@@ -293,8 +295,8 @@ ergm_MCMC_sample <- function(state, control, theta=NULL,
     if(!is.null(nws)) nws <- map(outl, "saved")
     
     if(control.parallel$MCMC.runtime.traceplot){
-      lapply(sms, function(sm) NVL3(theta, ergm.estfun(sm, ., as.ergm_model(state0[[1]])
-exclude=control$MCMC.esteq.exclude.statistics), sm[,!as.ergm_model(state0[[1]])$etamap$offsetmap,drop=FALSE])) %>% lapply(mcmc, start=control.parallel$MCMC.burnin+1, thin=control.parallel$MCMC.interval) %>% as.mcmc.list() %>% window(., thin=thin(.)*max(1,floor(niter(.)/1000))) %>% plot(ask=FALSE,smooth=TRUE,density=FALSE)
+      lapply(sms, function(sm) NVL3(theta, ergm.estfun(sm, ., as.ergm_model(state0[[1]]),
+        exclude=control$MCMC.esteq.exclude.statistics), sm[,!as.ergm_model(state0[[1]])$etamap$offsetmap,drop=FALSE])) %>% lapply(mcmc, start=control.parallel$MCMC.burnin+1, thin=control.parallel$MCMC.interval) %>% as.mcmc.list() %>% window(., thin=thin(.)*max(1,floor(niter(.)/1000))) %>% plot(ask=FALSE,smooth=TRUE,density=FALSE)
     }
   }
 
diff --git a/man/control.ergm.bridge.Rd b/man/control.ergm.bridge.Rd
index bce01c524..b5e05a773 100644
--- a/man/control.ergm.bridge.Rd
+++ b/man/control.ergm.bridge.Rd
@@ -26,6 +26,7 @@ control.ergm.bridge(
   obs.MCMC.prop.args = MCMC.prop.args,
   MCMC.maxedges = Inf,
   MCMC.packagenames = c(),
+  MCMC.esteq.exclude.statistics = NULL,
   term.options = list(),
   seed = NULL,
   parallel = 0,
@@ -51,7 +52,7 @@ control.logLik.ergm(
   obs.MCMC.prop = MCMC.prop,
   obs.MCMC.prop.weights = MCMC.prop.weights,
   obs.MCMC.prop.args = MCMC.prop.args,
-  MCMC.maxedges = NULL,
+  MCMC.maxedges = Inf,
   MCMC.packagenames = NULL,
   term.options = NULL,
   seed = NULL,
@@ -113,6 +114,9 @@ specifying additional arguments to proposal.}
 statistic functions in addition to those autodetected. This argument should
 not be needed outside of very strange setups.}
 
+\item{MCMC.esteq.exclude.statistics}{character string pattern. statistics with names that contain this string will be excluded from
+the estimating equations. For example, this supports tapering models which may have terms with non-specific means.}
+
 \item{term.options}{A list of additional arguments to be passed to term initializers. See \code{\link[=term.options]{? term.options}}.}
 
 \item{seed}{Seed value (integer) for the random number generator.  See

From 546408c1d225d0ada6859abb4fc8819eff1f13f1 Mon Sep 17 00:00:00 2001
From: handcock <handcock@ucla.edu>
Date: Fri, 26 Aug 2022 19:21:10 -0700
Subject: [PATCH 7/8] Added tests and fixed important typo

It now passes R CMD check
---
 DESCRIPTION            | 1 -
 R/ergm.getMCMCsample.R | 4 ++--
 2 files changed, 2 insertions(+), 3 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 2f7c29d46..109b33bf4 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -49,7 +49,6 @@ Suggests:
     testthat (>= 3.0.2),
     tergm (>= 4.0.0),
     ergm.count (>= 4.0),
-    ergm.userterms (>= 3.10.0),
     ergm.tapered (>= 1.0.0),
     networkDynamic (>= 0.10.1),
     withr,
diff --git a/R/ergm.getMCMCsample.R b/R/ergm.getMCMCsample.R
index a4a4556a9..781945356 100644
--- a/R/ergm.getMCMCsample.R
+++ b/R/ergm.getMCMCsample.R
@@ -239,8 +239,8 @@ ergm_MCMC_sample <- function(state, control, theta=NULL,
         if(verbose) message("Increasing thinning to ",interval,".")
       }
       
-      esteq <- lapply(sms, function(sm) NVL3(theta, ergm.estfun(sm, ., as.ergm_model(state[[1]]),
-        exclude=control$MCMC.esteq.exclude.statistics), sm[,!as.ergm_model(state[[1]])$etamap$offsetmap,drop=FALSE])) %>%
+      esteq <- lapply(sms, function(sm) NVL3(theta, ergm.estfun(sm, ., as.ergm_model(state0[[1]]),
+        exclude=control$MCMC.esteq.exclude.statistics), sm[,!as.ergm_model(state0[[1]])$etamap$offsetmap,drop=FALSE])) %>%
         lapply.mcmc.list(mcmc, start=1, thin=interval) %>% lapply.mcmc.list(`-`)
 
       if("Taper_Penalty" %in% dimnames(esteq[[1]])[[2]]) {browser()}

From b1ad98175e580ae828561badce1a049651e38f23 Mon Sep 17 00:00:00 2001
From: handcock <handcock@ucla.edu>
Date: Sun, 28 Aug 2022 13:08:59 -0700
Subject: [PATCH 8/8] Removed simplification that roils parallel

---
 R/ergm.MCMLE.R | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/R/ergm.MCMLE.R b/R/ergm.MCMLE.R
index 27acb00bf..25870384a 100644
--- a/R/ergm.MCMLE.R
+++ b/R/ergm.MCMLE.R
@@ -297,12 +297,7 @@ ergm.MCMLE <- function(init, nw, model,
     }
 
     # Compute the sample estimating equations and the convergence p-value. 
-    #esteqs <- ergm.estfun(statsmatrices, theta=mcmc.init, model=model, exclude=control$MCMC.esteq.exclude.statistics)
-    esteqs <- ergm.estfun(statsmatrices, theta=mcmc.init, model=model)
-    if(!is.null(control$MCMC.esteq.exclude.statistics)){
-      x.exclude <- match(control$MCMC.esteq.exclude.statistics,colnames(esteqs[[1]]))
-      if(!is.na(x.exclude)){ esteqs[[1]] <- esteqs[[1]][,-x.exclude,drop=FALSE] }
-    }
+    esteqs <- ergm.estfun(statsmatrices, theta=mcmc.init, model=model, exclude=control$MCMC.esteq.exclude.statistics)
     esteq <- as.matrix(esteqs)
     if(isTRUE(all.equal(apply(esteq,2,stats::sd), rep(0,ncol(esteq)), check.names=FALSE))&&!all(esteq==0))
       stop("Unconstrained MCMC sampling did not mix at all. Optimization cannot continue.")