statnet · handcock · Jul 1, 2021 · Jul 1, 2021 · Jul 9, 2021 · Jul 17, 2021
diff --git 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
@@ -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
@@ -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
@@ -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.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.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.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.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,20 +153,23 @@ 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,
                          lognormal=llik.grad.obs.IS,
                          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,
                         lognormal=llik.hessian.obs.IS,
                         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,20 +178,23 @@ 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,
                          lognormal=llik.grad.IS,
                          logtaylor=llik.grad.IS,
                          Median.Likelihood=llik.grad.IS,
                          EF.Likelihood=llik.grad.IS,
+                         Kpenalty=ergm.tapered::llik.grad.Kpenalty,
                          llik.grad.IS)
     Hessianfn <- switch(metric,
                         Likelihood=llik.hessian.IS,
                         lognormal=llik.hessian.IS,
                         logtaylor=llik.hessian.IS,
                         Median.Likelihood=llik.hessian.IS,
                         EF.Likelihood=llik.hessian.IS,
+                        Kpenalty=ergm.tapered::llik.hessian.Kpenalty,
                         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
     #
@@ -335,12 +330,11 @@ 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,
-                        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,