Improved code for valued ERGM

Improved code for valued ERGM

Minor documentation improvements

DESCRIPTION

@@ -4,8 +4,8 @@ Title: Tapered Exponential-Family Models for Networks
 Version: 1.1-0
 Date: 2022-08-31
 Authors@R: c(
-	  person("Mark S.", "Handcock", role=c("aut","cre"), email="[email protected]"),
-	  person("Pavel N.", "Krivitsky", role=c("aut"), email="pavel@uow.edu.au"),
+	  person("Mark S.", "Handcock", role=c("aut","cre"), email="[email protected]", comment=c(ORCID="0000-0002-9985-2785")),
+	  person("Pavel N.", "Krivitsky", role=c("aut"), email="pavel@statnet.org"),
 	  person("Ian", "Fellows", role=c("aut"), email="[email protected]"))
 Description: A set of terms and functions implementing Tapered exponential-family random graph models (ERGMs). Tapered ERGMs are a modification of ERGMs that reduce the effects of phase transitions,
     and with properly chosen hyper-parameters, provably removes all multiphase

R/InitWtErgmTerm.Kpenalty.R

@@ -0,0 +1,31 @@
+################################################################################
+
+#' @templateVar name WeightedKpenalty
+#' @title Weighted Kpenalty 
+#' @description 
+#'    Adds the terms specified in \code{formula} to the model \emph{and}
+#'    imposes the quadratic penalty of Fellows and Handcock
+#'    (2017). 
+#' @usage
+#' # binary: Kpenalty(formula=NULL, coef=NULL, m=NULL)
+#'
+#' @param formula a valid formula for a  standard ERGM
+#' @param coef a numeric vector of coefficients giving the penalty coefficients \eqn{\beta} for the tapering of the terms.
+#'   If \code{NULL} is passed, the tapering coefficients are set to \code{1/(4*summary(formula))}, the default
+#'   in Fellows and Handcock (2017).
+#'   If a numeric vector is given, there are interpreted as the tapering coefficients of the terms in the
+#'   model, including the terms enclosed in \code{offset()}.
+#'   If a numeric scalar is given, it is interpreted as a multiplier of the default tapering coefficients , that is,
+#' @param m numeric vector. If given, is
+#'    the value of the network statistic relative to which the model is
+#'    tapered. If omitted, it defaults to that of the model's LHS
+#'    network if \code{formula} is one-sided and that of the network on
+#'    the LHS of \code{formula} if it is two-sided.
+#'
+#' @concept tapered
+#' @import ergm statnet.common network stats
+InitWtErgmTerm.Kpenalty <- function(nw, arglist, response=NULL, ...){
+  out <- InitErgmTerm.Kpenalty(nw, arglist, response=response, ...)
+  out$name <- "wttaper_term"
+  out
+}

R/ergm.tapered-package.R

@@ -4,6 +4,9 @@
 #' diagnose, and simulate from Tapered exponential-family random graph models (ERGMs).
 #' For a list of functions type: \code{help(package='ergm.tapered')}
 #' 
+#' A good place to start is the vignette at \url{https://github.com/statnet/ergm.tapered}
+#' and the first two referenced papers below.
+#
 #' For a complete list of the functions, use \code{library(help="ergm.tapered")} or
 #' read the rest of the manual. 
 #' 
@@ -23,14 +26,14 @@
 #' \code{\link{mcmc.diagnostics}} and \code{\link{gof.ergm.tapered}}).
 #' 
 #' For detailed information on how to download and install the software, go to
-#' the \code{\link[=ergm-package]{ergm}} website: \url{https://statnet.org}. A
+#' the \code{\link[=ergm.tapered-package]{ergm.tapered}} website: \url{https://statnet.org}. A
 #' tutorial, support newsgroup, references and links to further resources are
 #' provided there.
 #' 
 #' @name ergm.tapered-package
 #' @docType package
-#' @author Mark S. Handcock \email{handcock@@stat.ucla.edu},\cr Pavel N. Krivitsky
-#' \email{krivitsky@@stat.psu.edu}, and\cr Ian E. Fellows
+#' @author Mark S. Handcock \email{[email protected]},\cr Pavel N. Krivitsky
+#' \email{[email protected]}, and\cr Ian E. Fellows
 #' \email{[email protected]}
 #' 
 #' Maintainer: Mark S. Handcock \email{handcock@@stat.ucla.edu}

R/ergm.tapered.R

@@ -85,7 +85,11 @@ ergm.tapered <- function(formula, r=2, beta=NULL, tau=NULL, tapering.centers=NUL
 
   # Determine the dyadic independence terms
   nw <- ergm.getnetwork(formula)
-  m<-ergm_model(formula, nw, response=response, ...)
+  primary.response <- response
+  ergm_preprocess_response(nw,primary.response)
+
+  proposal <- list(auxiliaries=NULL)
+  m<-ergm_model(formula, nw, x=NVL3(proposal$auxiliaries,list(proposal=.)), term.options=control$term.options, ...)
 
   if(is.null(tapering.centers)) tapering.centers <- target.stats
 
@@ -219,7 +223,7 @@ ergm.tapered <- function(formula, r=2, beta=NULL, tau=NULL, tapering.centers=NUL
     message("\n")
   }
 
-  re.names <- names(summary(newformula))
+  re.names <- names(summary(newformula, response=response))
   if(!fixed){
     control$init <- c(control$init,1)
     names(control$init)[length(control$init)] <- "Taper_Penalty"
@@ -233,12 +237,14 @@ ergm.tapered <- function(formula, r=2, beta=NULL, tau=NULL, tapering.centers=NUL
     }
   }
 
-  # fit ergm
-  fit.MPLE.control <- control
-  fit.MPLE.control$init <- NULL
-  fit.MPLE.control$MPLE.save.xmat <- TRUE
-  fit.MPLE <- ergm(reformula, control=fit.MPLE.control, estimate="MPLE",
-                   response=response, constraints=constraints, reference=reference, eval.loglik=eval.loglik, verbose=verbose, ...)
+  if(!is.valued(nw)){
+  # fit binary ergm
+    fit.MPLE.control <- control
+    fit.MPLE.control$init <- NULL
+    fit.MPLE.control$MPLE.save.xmat <- TRUE
+    fit.MPLE <- ergm(reformula, control=fit.MPLE.control, estimate="MPLE",
+                     response=response, constraints=constraints, reference=reference, eval.loglik=eval.loglik, verbose=verbose, ...)
+  }
   if(is.null(target.stats)){
     fit <- ergm(newformula, control=control,
                 response=response, constraints=constraints, reference=reference, eval.loglik=eval.loglik, verbose=verbose, ...)
@@ -287,8 +293,10 @@ ergm.tapered <- function(formula, r=2, beta=NULL, tau=NULL, tapering.centers=NUL
   fcoef[seq_along(fulltau)[!is.na(nm)]] <- fcoef[nm[!is.na(nm)]]
   fit$tapering.coefficients <- fulltau
   fit$taudelta.offset <- 2*fulltau*as.vector(apply(sapply(fit$sample,function(x){apply(x[,-ncol(x)],2,sd)}),1,mean))
-  fit$taudelta.mean <- apply((2*fit.MPLE$glm.result$value$model[,1]-1)*sweep(fit.MPLE$xmat.full,2,fulltau,"*"),2,weighted.mean,weight=fit.MPLE$glm.result$value$prior.weights)
-  fit$taudelta.mad <- apply((2*fit.MPLE$glm.result$value$model[,1]-1)*sweep(abs(fit.MPLE$xmat.full),2,fulltau,"*"),2,weighted.mean,weight=fit.MPLE$value$glm.result$prior.weights)
+  if(!is.valued(nw)){
+    fit$taudelta.mean <- apply((2*fit.MPLE$glm.result$value$model[,1]-1)*sweep(fit.MPLE$xmat.full,2,fulltau,"*"),2,weighted.mean,weight=fit.MPLE$glm.result$value$prior.weights)
+    fit$taudelta.mad <- apply((2*fit.MPLE$glm.result$value$model[,1]-1)*sweep(abs(fit.MPLE$xmat.full),2,fulltau,"*"),2,weighted.mean,weight=fit.MPLE$value$glm.result$prior.weights)
+  }
 
   # post process fit to alter Hessian etc
   if(is.null(tapering.centers)){