This is a transference of the ergm-private" "control_llik" branch to the public version ergm 4.0.

DESCRIPTION

@@ -49,6 +49,7 @@ Suggests:
     testthat (>= 3.0.2),
     tergm (>= 4.0.0),
     ergm.count (>= 4.0),
+    ergm.tapered (>= 1.0.0),
     networkDynamic (>= 0.10.1),
     withr,
     covr,

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.
@@ -275,6 +278,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?
@@ -420,6 +425,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
@@ -478,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(),
@@ -500,6 +510,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,
@@ -540,12 +551,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,
@@ -685,6 +697,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.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)

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,

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}}.
@@ -79,6 +81,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",
@@ -114,7 +117,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,
@@ -140,6 +143,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",

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)

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)
@@ -126,6 +126,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)

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,
                        verbose=verbose,
                        estimateonly=!calc.MCSE)
@@ -549,7 +552,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)

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"

R/ergm.bridge.R

@@ -154,6 +154,12 @@ ergm.bridge.llr<-function(object, response=NULL, reference=~Bernoulli, constrain
 
   ## Miscellaneous settings
   Dtheta.Du <- (to-from)[!state[[1]]$model$etamap$offsettheta]
+  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.
   if(!is.null(target.stats)){
@@ -167,8 +173,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=control$MCMC.esteq.exclude.statistics), `%*%`, Dtheta.Du)
+    }else{
+      sum(ergm.estfun(samp, theta, state[[1]]$model$etamap, exclude=control$MCMC.esteq.exclude.statistics) * Dtheta.Du)
+    }
   }