diff --git a/.gitignore b/.gitignore index 31b6ba1..580e89e 100644 --- a/.gitignore +++ b/.gitignore @@ -95,3 +95,4 @@ dkms.conf .Rproj.user *.print R/ergm.tapered.R.withcomments +R/ergm.tapered.R.nosancall diff --git a/DESCRIPTION b/DESCRIPTION index 4d23d7a..550157f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -23,7 +23,7 @@ Suggests: testthat License: GPL-3 Roxygen: list(markdown = TRUE) Encoding: UTF-8 -RoxygenNote: 7.2.1 +RoxygenNote: 7.2.3 Remotes: github::statnet/statnet.common@master, github::statnet/network@master, github::statnet/ergm@tapered diff --git a/NAMESPACE b/NAMESPACE index 5b8e969..40187aa 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -6,6 +6,7 @@ S3method(summary,ergm.tapered) export(control.ergm.tapered) export(ergm.mvtapered) export(ergm.tapered) +export(ergm.tapered.formula) export(llik.fun.Kpenalty) export(llik.fun.obs.Kpenalty) export(llik.grad.Kpenalty) @@ -14,6 +15,7 @@ export(llik.grad.obs.Kpenalty.numDeriv) export(llik.hessian.Kpenalty) export(llik.hessian.Kpenalty.numDeriv) export(llik.hessian.obs.Kpenalty.numDeriv) +export(simulate_ergm.tapered) import(ergm) import(network) import(statnet.common) diff --git a/R/control.ergm.tapered.R b/R/control.ergm.tapered.R index d58a629..4df3fb0 100644 --- a/R/control.ergm.tapered.R +++ b/R/control.ergm.tapered.R @@ -40,6 +40,8 @@ #' errors are below the precision bound, and the Hummel step length is 1 for #' two consecutive iterations. This is an experimental feature [ergm()]. #' The default value in [ergm.tapered()] is 0.15, higher than in [ergm()]. +#' @param estimate.tapered.bias logical. It \code{TRUE} the bias of the estimated estimates due to tapering is estimated. +#' If this is \code{FALSE} the MPLE of the untapered model is not computed, saving computational time. #' @param loglik list List of additional control arguments for the loglik. See \code{\link{control.logLik.ergm.tapered}}. #' @param \dots Additional arguments, passed to other functions This argument #' is helpful because it collects any control parameters that have been @@ -71,6 +73,7 @@ control.ergm.tapered<-function( MCMLE.kurtosis.penalty=2.0, MCMLE.termination=c("precision","confidence", "Hummel", "Hotelling", "none"), MCMLE.MCMC.precision=0.005, + estimate.tapered.bias=TRUE, # MCMC.effectiveSize.maxruns=8, loglik=control.logLik.ergm.tapered(), ... @@ -88,4 +91,4 @@ control.ergm.tapered<-function( set.control.class(c("control.ergm", "control.ergm.tapered")) } -STATIC_TAPERING_CONTROLS <- c("MCMLE.kurtosis.location", "MCMLE.kurtosis.scale", "MCMLE.kurtosis.penalty") +STATIC_TAPERING_CONTROLS <- c("MCMLE.kurtosis.location", "MCMLE.kurtosis.scale", "MCMLE.kurtosis.penalty", "estimate.tapered.bias") diff --git a/R/ergm.tapered.R b/R/ergm.tapered.R old mode 100755 new mode 100644 index 48c65e0..df99eb5 --- a/R/ergm.tapered.R +++ b/R/ergm.tapered.R @@ -22,16 +22,23 @@ #' @param response {Name of the edge attribute whose value is to be #' modeled in the valued ERGM framework. Defaults to \code{NULL} for #' simple presence or absence, modeled via a binary ERGM.} -#' @param reference {A one-sided formula specifying -#' the reference measure (\eqn{h(y)}) to be used. -#' See help for \link[=ergm-references]{ERGM reference measures} implemented in the -#' \strong{\link[=ergm-package]{ergm}} package.} #' @param constraints {A formula specifying one or more constraints #' on the support of the distribution of the networks being modeled, #' using syntax similar to the \code{formula} argument, on the -#' right-hand side. See \link[=ergm-constraints]{ERGM constraints} -#' \code{\link{ergm}} for details.} -#' @param control An object of class control.ergm. Passed to the ergm function. +#' right-hand side. See \link[=ergm-constraints]{ERGM constraints}.} +#' @param reference {A one-sided formula specifying +#' the reference measure (\eqn{h(y)}) to be used. +#' See help for \link[=ergm-references]{ERGM reference measures} implemented in the +#' \strong{\link[=ergm-package]{ergm}} package for details.} +#' @param estimate {If "MPLE," then the maximum pseudolikelihood estimator +#' is returned. If "MLE" (the default), then an approximate maximum likelihood +#' estimator is returned. For certain models, the MPLE and MLE are equivalent, +#' in which case this argument is ignored. (To force MCMC-based approximate +#' likelihood calculation even when the MLE and MPLE are the same, see the +#' \code{force.main} argument of \code{\link{control.ergm}}. If "CD" (\emph{EXPERIMENTAL}), +#' the Monte-Carlo contrastive divergence estimate is returned. ) +#' } +#' @param control An object of class control.erg.tapered. Passed to the ergm function. #' @param fixed A `logical`: if this is \code{TRUE}, the tapering is fixed at the passed values. #' If this is \code{FALSE}, the tapering is estimated using a kurtosis penalized likelihood. #' See Blackburn and Handcock (2022) for an explanation. @@ -66,6 +73,7 @@ ergm.tapered <- function(formula, r=2, beta=NULL, tau=NULL, tapering.centers=NULL, target.stats=NULL, family="taper", taper.terms="all", response=NULL, constraints=~., reference=~Bernoulli, + estimate=c("MLE", "MPLE", "CD"), control = control.ergm.tapered(), fixed=TRUE, verbose=FALSE, eval.loglik=TRUE, ...){ if(!methods::is(control,"control.ergm.tapered")){ @@ -93,6 +101,9 @@ ergm.tapered <- function(formula, r=2, beta=NULL, tau=NULL, tapering.centers=NUL if(is.null(tapering.centers)) tapering.centers <- target.stats + estimate <- match.arg(estimate) + arg_list <- list(...) + arg_list$estimate <- NULL # set tapering terms if(is.character(taper.terms) & length(taper.terms)==1){ @@ -126,7 +137,11 @@ ergm.tapered <- function(formula, r=2, beta=NULL, tau=NULL, tapering.centers=NUL reformula <- append_rhs.formula(trimmed_formula,taper_formula) } attr(taper.terms,"env") <- NULL - taper.stats <- summary(append_rhs.formula(nw ~.,taper.terms), response=response, ...) + if(is.null(target.stats)){ + taper.stats <- summary(append_rhs.formula(nw ~.,taper.terms), response=response, ...) + }else{ + taper.stats <- target.stats + } # set tapering coefficient tau <- switch(family, @@ -230,20 +245,58 @@ ergm.tapered <- function(formula, r=2, beta=NULL, tau=NULL, tapering.centers=NUL } } - if(!is.valued(nw)){ + if(!is.valued(nw) & control$estimate.tapered.bias){ # 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(length(arg_list)>0){ + fit.MPLE <- ergm(reformula, control=fit.MPLE.control, estimate="MPLE", + response=response, constraints=constraints, reference=reference, eval.loglik=eval.loglik, verbose=verbose, arg_list) + }else{ + 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)){ + if(!is.null(target.stats) & estimate != "MPLE"){ + target.stats.aug <- c(target.stats,NA) + names(target.stats.aug) <- re.names + names(target.stats) <- re.names[1:length(target.stats)] + for( i in 1:20){ + sanformula <- statnet.common::nonsimp_update.formula(newformula,nw ~ .) + env$.taper.center <- target.stats + env$.taper.coef <- tau + env$nw <- nw + environment(sanformula) <- env + if(verbose){ + message(sprintf("Computing SAN network\n")) + } + nw <- san(sanformula,target.stats=target.stats, + control=control.san(SAN.maxit=10, SAN.exclude.statistics="Taper_Penalty", parallel=control$parallel)) + } + env$.taper.center <- target.stats + env$.taper.coef <- tau + env$nw <- nw + environment(sanformula) <- env + #if(verbose){ + message(sprintf("SAN network compared to target statistics:\n")) + print(cbind(target.stats, summary(sanformula)[-length(target.stats)-1])) + #} + newformula <- sanformula + #tostats <- ostats + #names(tostats) <- re.names + #control$SAN$SAN.exclude.statistics <- "Taper_Penalty" + #fit <- ergm(newformula, control=control, #target.stats=tostats, #offset.coef=tau, + # response=response, constraints=constraints, reference=reference, eval.loglik=eval.loglik, verbose=verbose, ...) + } + if(length(arg_list)>0){ fit <- ergm(newformula, control=control, - response=response, constraints=constraints, reference=reference, eval.loglik=eval.loglik, verbose=verbose, ...) + response=response, constraints=constraints, reference=reference, estimate=estimate, + eval.loglik=eval.loglik, verbose=verbose, arg_list) }else{ - fit <- ergm(newformula, control=control, target.stats=ostats, offset.coef=tau, - response=response, constraints=constraints, reference=reference, eval.loglik=eval.loglik, verbose=verbose, ...) + fit <- ergm(newformula, control=control, + response=response, constraints=constraints, reference=reference, estimate=estimate, + eval.loglik=eval.loglik, verbose=verbose) } fit$tapering.centers <- taper.stats @@ -269,6 +322,7 @@ ergm.tapered <- function(formula, r=2, beta=NULL, tau=NULL, tapering.centers=NUL } fit$orig.formula <- formula + if(estimate == "MLE"){ if(fixed){ sample <- as.matrix(fit$sample)[,1:npar,drop=FALSE] fit$hessian <- fit$hessian[1:npar,1:npar] @@ -286,9 +340,12 @@ 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)) - if(!is.valued(nw)){ + if(!is.valued(nw) & control$estimate.tapered.bias){ 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) + }else{ + fit$taudelta.mean <- rep(NA, length(fulltau)) + fit$taudelta.mad <- rep(NA, length(fulltau)) } # post process fit to alter Hessian etc @@ -310,6 +367,7 @@ ergm.tapered <- function(formula, r=2, beta=NULL, tau=NULL, tapering.centers=NUL fit$hessian <- -fit$hessian } } + } class(fit) <- c("ergm.tapered",family,class(fit)) diff --git a/R/ergm.tapered.formula.R b/R/ergm.tapered.formula.R new file mode 100644 index 0000000..aa8e2fe --- /dev/null +++ b/R/ergm.tapered.formula.R @@ -0,0 +1,181 @@ +#' Create an Tapered ERGM formula specifying an Tapered ERGM model as a standard ERGM +#' @param formula An ergm formula to fit +#' @param r The scaling factor to use for the heuristic of setting beta equal to r standard deviations of the observed statistics +#' @param beta The tapering parameters, expressed as in Fellows and Handcock (2017). If not NULL, these override the heuristics (r). +#' @param tau The tapering parameters, expressed as natural parameters. If not NULL, these override the beta and the heuristics (r). +#' @param tapering.centers The centers of the tapering terms. If null, these are taken to be the mean value parameters. +#' @param target.stats {vector of "observed network statistics," +#' if these statistics are for some reason different than the +#' actual statistics of the network on the left-hand side of +#' \code{formula}. +#' Equivalently, this vector is the mean-value parameter values for the +#' model. If this is given, the algorithm finds the natural +#' parameter values corresponding to these mean-value parameters. +#' If \code{NULL}, the mean-value parameters used are the observed +#' statistics of the network in the formula. +#' } +#' @param family The type of tapering used. This should either be the \code{stereo} or \code{taper}, the +#' tapering model of Fellows and Handcock (2017). +#' @param taper.terms Specification of the tapering used. If the character variable "dependence" then all the dependent +#' terms are tapered. If the character variable "all" then all terms are tapered. +#' It can also be the RHS of a formula giving the terms to be tapered. +#' @param response {Name of the edge attribute whose value is to be +#' modeled in the valued ERGM framework. Defaults to \code{NULL} for +#' simple presence or absence, modeled via a binary ERGM.} +#' @param constraints {A formula specifying one or more constraints +#' on the support of the distribution of the networks being modeled, +#' using syntax similar to the \code{formula} argument, on the +#' right-hand side. See \link[=ergm-constraints]{ERGM constraints}.} +#' @param reference {A one-sided formula specifying +#' the reference measure (\eqn{h(y)}) to be used. +#' See help for \link[=ergm-references]{ERGM reference measures} implemented in the +#' \strong{\link[=ergm-package]{ergm}} package for details.} +#' @param verbose A `logical`: if this is +#' \code{TRUE}, the program will print out additional +#' information about the progress of estimation. +#' @param ... Additional arguments to \code{\link{ergm}}. +#' @returns +#' An formula containing a specification of the Tapered ERGM model. In addition to all of the ergm items, +#' this object contains tapering.centers, tapering.coef and orig.formula. tapering.centers are the centers for the tapering term. +#' tapering.coef are the tapering coefficients = 1/ beta^2. orig.formula is the formula passed into ergm.tapered. +#' @importFrom stats var as.formula +#' @references \itemize{ +#' * Fellows, I. and M. S. Handcock (2017), +#' Removing Phase Transitions from Gibbs Measures. Volume 54 of +#' Proceedings of Machine Learning Research, Fort Lauderdale, +#' FL, USA, pp. 289–297. PMLR. +#' * Blackburn, B. and M. S. Handcock (2022), +#' Practical Network Modeling via Tapered Exponential-family Random Graph Models. +#' Journal of Computational and Graphical Statistics +#' \doi{10.1080/10618600.2022.2116444}. +#' +#' } +#' @examples +#' \dontrun{ +#' data(sampson) +#' fit <- ergm.tapered.formula(samplike ~ edges + triangles()) +#' summary(fit) +#' } +#' @export +ergm.tapered.formula <- function(formula, r=2, beta=NULL, tau=NULL, tapering.centers=NULL, target.stats=NULL, + family="taper", taper.terms="all", + response=NULL, constraints=~., reference=~Bernoulli, + verbose=FALSE, ...){ + + # Determine the dyadic independence terms + nw <- ergm.getnetwork(formula) + m<-ergm_model(formula, nw, response=response, ...) + + if(is.null(tapering.centers)) tapering.centers <- target.stats + + # set tapering terms + if(is.character(taper.terms) & length(taper.terms)==1){ + if(taper.terms=="dependent"){ + a <- sapply(m$terms, function(term){is.null(term$dependence) || term$dependence}) + taper.terms <- list_rhs.formula(formula) + tmp <- taper.terms + taper.terms <- NULL + for(i in seq_along(tmp)){if(a[i]){taper.terms <- c(taper.terms,tmp[[i]])}} + taper_formula <- append_rhs.formula(~.,taper.terms, env=NULL) + trimmed_formula=suppressWarnings(filter_rhs.formula(formula, function(term,taper.terms){all(term != taper.terms)}, taper.terms)) + reformula <- append_rhs.formula(trimmed_formula,taper_formula, env=NULL) + }else{if(taper.terms=="all"){ + taper.terms <- list_rhs.formula(formula) + taper_formula <- append_rhs.formula(~.,taper.terms, env=NULL) + trimmed_formula=suppressWarnings(filter_rhs.formula(formula, function(term,taper.terms){all(term != taper.terms)}, taper.terms)) + reformula <- formula + }else{ + if(!inherits(taper.terms,"formula")){ + stop('taper.terms must be "dependent", "all" or a formula of terms.') + } + taper.terms <- list_rhs.formula(taper.terms) + taper_formula <- append_rhs.formula(~.,taper.terms, env=NULL) + trimmed_formula=suppressWarnings(filter_rhs.formula(formula, function(term,taper.terms){all(term != taper.terms)}, taper.terms)) + reformula <- append_rhs.formula(trimmed_formula,taper_formula, env=NULL) + }} + }else{ + taper_formula <- taper.terms + taper.terms <- list_rhs.formula(taper.terms) + trimmed_formula=suppressWarnings(filter_rhs.formula(formula, function(term,taper.terms){all(term != taper.terms)}, taper.terms)) + reformula <- append_rhs.formula(trimmed_formula,taper_formula, env=NULL) + } + attr(taper.terms,"env") <- NULL + if(is.null(target.stats)){ + taper.stats <- summary(append_rhs.formula(nw ~.,taper.terms), response=response, ...) + }else{ + taper.stats <- target.stats + } + + # set tapering coefficient + tau <- switch(family, + "stereo"={ + if(is.null(beta) & is.null(tau)){ + 1 + }else{ + if(is.null(tau)){ + beta + }else{ + tau + } + }}, + {if(is.null(tau)){ + if(is.null(beta)){ + 1 / (r^2 * pmax(1,abs(taper.stats))) + }else{ + 1 / beta^2 + }}else{tau}} + ) + if(family=="stereo"){ + names(tau) <- "stereo" + }else{ + names(tau) <- names(taper.stats) + } + + taper_terms <- switch(paste0(family,ifelse(fixed,"_fixed","_notfixed")), + "stereo_fixed"=statnet.common::nonsimp_update.formula(taper_formula,.~Stereo(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE), + "stereo_notfixed"=statnet.common::nonsimp_update.formula(taper_formula,.~Stereo(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE), + "taper_fixed"=statnet.common::nonsimp_update.formula(taper_formula,.~Taper(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE), + "taper_notfixed"=statnet.common::nonsimp_update.formula(taper_formula,.~Kpenalty(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE), + statnet.common::nonsimp_update.formula(taper_formula,.~Taper(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE) + ) + + if(length(list_rhs.formula(formula))==length(taper.terms)){ + newformula <- switch(paste0(family,ifelse(fixed,"_fixed","_notfixed")), + "stereo_fixed"=statnet.common::nonsimp_update.formula(formula,.~Stereo(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE), + "stereo_notfixed"=statnet.common::nonsimp_update.formula(formula,.~Stereo(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE), + "taper_fixed"=statnet.common::nonsimp_update.formula(formula,.~Taper(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE), + "taper_notfixed"=statnet.common::nonsimp_update.formula(formula,.~Kpenalty(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE), + statnet.common::nonsimp_update.formula(formula,.~Taper(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE) + ) + + }else{ + newformula <- append_rhs.formula(trimmed_formula,taper_terms, env=NULL) + } + ostats <- summary(reformula, response=response, ...) + if(!is.null(tapering.centers)){ + tmp <- match(names(ostats), names(tapering.centers)) + if(any(is.na(tmp))){ + stop(paste("tapering.centers needs to have a named center for each statistic in the model:", + names(ostats))) + } + ostats <- tapering.centers[tmp] + } + npar <- length(ostats) + + env <- new.env(parent=environment(formula)) + env$.taper.center <- taper.stats + env$.taper.coef <- tau + environment(newformula) <- env + + return(newformula) +} diff --git a/R/simulate_ergm.tapered.R b/R/simulate_ergm.tapered.R new file mode 100755 index 0000000..d932f03 --- /dev/null +++ b/R/simulate_ergm.tapered.R @@ -0,0 +1,221 @@ +#' Simulates from a Tapered ERGM formula +#' @param formula An ergm formula to fit +#' @param coef Vector of parameter values for the model from which the +#' sample is to be drawn. +#' @param nsim Number of networks to be randomly drawn from the given +#' distribution on the set of all networks, returned by the Metropolis-Hastings +#' algorithm. +#' @param r The scaling factor to use for the heuristic of setting beta equal to r standard deviations of the observed statistics +#' @param beta The tapering parameters, expressed as in Fellows and Handcock (2017). If not NULL, these override the heuristics (r). +#' @param tau The tapering parameters, expressed as natural parameters. If not NULL, these override the beta and the heuristics (r). +#' @param tapering.centers The centers of the tapering terms. If null, these are taken to be the mean value parameters. +#' Equivalently, this vector is the mean-value parameter values for the +#' model. If this is given, the algorithm finds the natural +#' parameter values corresponding to these mean-value parameters. +#' If \code{NULL}, the mean-value parameters used are the observed +#' statistics of the network in the formula. +#' } +#' @param family The type of tapering used. This should either be the \code{stereo} or \code{taper}, the +#' tapering model of Fellows and Handcock (2017). +#' @param taper.terms Specification of the tapering used. If the character variable "dependence" then all the dependent +#' terms are tapered. If the character variable "all" then all terms are tapered. +#' It can also be the RHS of a formula giving the terms to be tapered. +#' @param response {Name of the edge attribute whose value is to be +#' modeled in the valued ERGM framework. Defaults to \code{NULL} for +#' simple presence or absence, modeled via a binary ERGM.} +#' @param reference {A one-sided formula specifying +#' the reference measure (\eqn{h(y)}) to be used. +#' See help for \link[=ergm-references]{ERGM reference measures} implemented in the +#' \strong{\link[=ergm-package]{ergm}} package.} +#' @param constraints {A formula specifying one or more constraints +#' on the support of the distribution of the networks being modeled, +#' using syntax similar to the \code{formula} argument, on the +#' right-hand side. See \link[=ergm-constraints]{ERGM constraints} +#' \code{\link{ergm}} for details.} +#' @param seed the random seed to use; see [simulate()]. +#' @param control An object of class control.ergm. Passed to the ergm function. +#' @param verbose A `logical`: if this is +#' \code{TRUE}, the program will print out additional +#' information about the progress of estimation. +#' @param ... Additional arguments to \code{\link{ergm}}. +#' @returns +#' An object of class c('tapered.ergm','ergm') containing the fit model. In addition to all of the ergm items, +#' this object contains tapering.centers, tapering.coef and orig.formula. tapering.centers are the centers for the tapering term. +#' tapering.coef are the tapering coefficients = 1/ beta^2. orig.formula is the formula passed into ergm.tapered. +#' @importFrom stats var as.formula +#' @references \itemize{ +#' * Fellows, I. and M. S. Handcock (2017), +#' Removing Phase Transitions from Gibbs Measures. Volume 54 of +#' Proceedings of Machine Learning Research, Fort Lauderdale, +#' FL, USA, pp. 289–297. PMLR. +#' * Blackburn, B. and M. S. Handcock (2022), +#' Practical Network Modeling via Tapered Exponential-family Random Graph Models. +#' Journal of Computational and Graphical Statistics +#' \doi{10.1080/10618600.2022.2116444}. +#' +#' } +#' @examples +#' \dontrun{ +#' data(sampson) +#' fit <- ergm.tapered(samplike ~ edges + triangles()) +#' summary(fit) +#' } +#' @export +simulate_ergm.tapered <- function(formula, coef, nsim=1, r=2, beta=NULL, tau=NULL, tapering.centers=NULL, + family="taper", taper.terms="all", + response=NULL, constraints=~., reference=~Bernoulli, seed=NULL, + control = control.simulate.formula(), verbose=FALSE, ...){ + + # Needed by ergm.getMCMCsample + match.llik.arg.pars <- c("MCMLE.metric","MCMLE.termination","MCMC.esteq.exclude.statistics",STATIC_TAPERING_CONTROLS) + for(arg in match.llik.arg.pars) + control$loglik[arg]<-list(control[[arg]]) + + # Determine the dyadic independence terms + nw <- ergm.getnetwork(formula) + m<-ergm_model(formula, nw, response=response, ...) + + # set tapering terms + if(is.character(taper.terms) & length(taper.terms)==1){ + if(taper.terms=="dependent"){ + a <- sapply(m$terms, function(term){is.null(term$dependence) || term$dependence}) + taper.terms <- list_rhs.formula(formula) + tmp <- taper.terms + taper.terms <- NULL + for(i in seq_along(tmp)){if(a[i]){taper.terms <- c(taper.terms,tmp[[i]])}} + taper_formula <- append_rhs.formula(~.,taper.terms, env=NULL) + trimmed_formula=suppressWarnings(filter_rhs.formula(formula, function(term,taper.terms){all(term != taper.terms)}, taper.terms)) + reformula <- append_rhs.formula(trimmed_formula,taper_formula, env=NULL) + }else{if(taper.terms=="all"){ + taper.terms <- list_rhs.formula(formula) + taper_formula <- append_rhs.formula(~.,taper.terms, env=NULL) + trimmed_formula=suppressWarnings(filter_rhs.formula(formula, function(term,taper.terms){all(term != taper.terms)}, taper.terms)) + reformula <- formula + }else{ + if(!inherits(taper.terms,"formula")){ + stop('taper.terms must be "dependent", "all" or a formula of terms.') + } + taper.terms <- list_rhs.formula(taper.terms) + taper_formula <- append_rhs.formula(~.,taper.terms, env=NULL) + trimmed_formula=suppressWarnings(filter_rhs.formula(formula, function(term,taper.terms){all(term != taper.terms)}, taper.terms)) + reformula <- append_rhs.formula(trimmed_formula,taper_formula, env=NULL) + }} + }else{ + taper_formula <- taper.terms + taper.terms <- list_rhs.formula(taper.terms) + trimmed_formula=suppressWarnings(filter_rhs.formula(formula, function(term,taper.terms){all(term != taper.terms)}, taper.terms)) + reformula <- append_rhs.formula(trimmed_formula,taper_formula, env=NULL) + } + attr(taper.terms,"env") <- NULL + if(is.null(tapering.centers)){ + taper.stats <- summary(append_rhs.formula(nw ~.,taper.terms), response=response, ...) + }else{ + taper.stats <- tapering.centers + } + + # set tapering coefficient + tau <- switch(family, + "stereo"={ + if(is.null(beta) & is.null(tau)){ + 1 + }else{ + if(is.null(tau)){ + beta + }else{ + tau + } + }}, + {if(is.null(tau)){ + if(is.null(beta)){ + 1 / (r^2 * pmax(1,abs(taper.stats))) + }else{ + 1 / beta^2 + }}else{tau}} + ) + if(family=="stereo"){ + names(tau) <- "stereo" + }else{ + names(tau) <- names(taper.stats) + } + + fixed <- TRUE + taper_terms <- switch(paste0(family,ifelse(fixed,"_fixed","_notfixed")), + "stereo_fixed"=statnet.common::nonsimp_update.formula(taper_formula,.~Stereo(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE), + "stereo_notfixed"=statnet.common::nonsimp_update.formula(taper_formula,.~Stereo(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE), + "taper_fixed"=statnet.common::nonsimp_update.formula(taper_formula,.~Taper(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE), + "taper_notfixed"=statnet.common::nonsimp_update.formula(taper_formula,.~Kpenalty(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE), + statnet.common::nonsimp_update.formula(taper_formula,.~Taper(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE) + ) + + if(length(list_rhs.formula(formula))==length(taper.terms)){ + newformula <- switch(paste0(family,ifelse(fixed,"_fixed","_notfixed")), + "stereo_fixed"=statnet.common::nonsimp_update.formula(formula,.~Stereo(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE), + "stereo_notfixed"=statnet.common::nonsimp_update.formula(formula,.~Stereo(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE), + "taper_fixed"=statnet.common::nonsimp_update.formula(formula,.~Taper(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE), + "taper_notfixed"=statnet.common::nonsimp_update.formula(formula,.~Kpenalty(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE), + statnet.common::nonsimp_update.formula(formula,.~Taper(~.,coef=.taper.coef,m=.taper.center), + environment(), from.new=TRUE) + ) + + }else{ + newformula <- append_rhs.formula(trimmed_formula,taper_terms, env=NULL) + } + ostats <- summary(reformula, response=response, ...) + if(!is.null(tapering.centers)){ + tmp <- match(names(ostats), names(tapering.centers)) + if(any(is.na(tmp))){ + stop(paste("tapering.centers needs to have a named center for each statistic in the model:", + names(ostats))) + } + ostats <- tapering.centers[tmp] + } + npar <- length(ostats) + + env <- new.env(parent=environment(formula)) + env$.taper.center <- taper.stats + env$.taper.coef <- tau + environment(newformula) <- env + + if(verbose){ + message(sprintf("The tapering formula is:\n %s", paste(deparse(newformula), sep="\n", collapse = "\n"))) + message("The (natural) tapering parameters are:") + for(i in seq_along(tau)){ + message(sprintf(" %s : %f",names(tau)[i],tau[i])) + } + message("\n") + } + + re.names <- names(summary(newformula)) + + # fit ergm + fit <- simulate(newformula, nsim=nsim, seed=seed, coef=coef, control=control, + response=response, constraints=constraints, reference=reference, verbose=verbose, ...) + + if(is.matrix(fit)){ + a <- grep("Taper_Penalty",colnames(fit),fixed=TRUE) + if(length(a) > 0){ + fit0 <- fit[,-a,drop=FALSE] + attr(fit0, "monitored") <- attr(fit0, "monitored")[-a] + attr(fit0, "formula") <- attr(fit0, "formula") + attr(fit0, "constraints") <- attr(fit0, "constraints") + attr(fit0, "reference") <- attr(fit0, "reference") + fit <- fit0 + } + } + attr(fit, "tapering.centers") <- taper.stats + attr(fit, "tapering.centers.o") <- ostats + attr(fit, "tapering.centers.t") <- taper.terms + attr(fit, "tapering.coef") <- tau + attr(fit, "r") <- r + + fit +} diff --git a/README.md b/README.md new file mode 100644 index 0000000..8e5b74c --- /dev/null +++ b/README.md @@ -0,0 +1,251 @@ +# `ergm.tapered`: Tapered Exponential-Family Models for Networks + +RDS network + +[![rstudio mirror downloads](https://cranlogs.r-pkg.org/badges/ergm.tapered?color=2ED968)](https://cranlogs.r-pkg.org/) +[![cran version](https://www.r-pkg.org/badges/version/ergm.tapered)](https://cran.r-project.org/package=ergm.tapered) +[![Coverage status](https://codecov.io/gh/statnet/ergm.tapered/branch/master/graph/badge.svg)](https://codecov.io/github/statnet/ergm.tapered?branch=master) +[![R build status](https://github.com/statnet/ergm.tapered/workflows/R-CMD-check/badge.svg)](https://github.com/statnet/ergm.tapered/actions) + +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 +behavior. + +Each ERGM has a corresponding Tapered ERGM. Indeed, the `ergm.tapered` package fits any `ergm` as it is based on `ergm` itself. + +# Installation + + + + + + + +To install the latest development version from github, you can also use: + +```{r} +# If devtools is not installed: +# install.packages("devtools") + +devtools::install_github("statnet/ergm.tapered") +``` + + +For now, this will install a variant of the current version of `ergm` (the `tapered` branch) that is needed. In a bit the CRAN version of `ergm` will have it. + +If you want to install that variant directly, you can use + +```{r} +devtools::install_github("statnet/ergm", branch="tapered") +``` + +# Implementation + +Load package and example data + +```{r} +library(ergm.tapered) +data(sampson) +``` + +Sampson (1969) recorded the social interactions among a group of monks while he was a resident as an experimenter at the cloister. + Of particular interest are the data on positive affect relations + ("liking," in which each monk was asked if he had positive relations + to each of the other monks. Each monk ranked only his top three + choices (or four, in the case of ties) on "liking". Here, we + consider a directed edge from monk A to monk B to exist if A + nominated B among these top choices at any one of three time points during the year. + For details see: + +``` +help(sampson) +``` + +We can make a quick visualization of the network + +```{r} +plot(sampson) +``` + +![](.github/figures/samplikeplot.png) + +A natural model is one that includes a term measuring the transitivity +of triples in the network, defined as a set of edges {(i,j), (j,k), +(i,k)}. + +``` r + fit <- ergm(samplike ~ edges + ttriple) +``` + + Starting maximum pseudolikelihood estimation (MPLE): + + ... + + Optimizing with step length 0.0054. + The log-likelihood improved by 27.9793. + Estimating equations are not within tolerance region. + Iteration 3 of at most 60: + Error in ergm.MCMLE(init, nw, model, initialfit = (initialfit <- NULL), : + Unconstrained MCMC sampling did not mix at all. Optimization cannot continue. + Calls: ergm -> ergm.MCMLE + In addition: Warning message: + In ergm_MCMC_sample(s, control, theta = mcmc.init, verbose = max(verbose - : + Unable to reach target effective size in iterations alotted. + Calls: ergm -> ergm.MCMLE -> ergm_MCMC_sample + + +This fit fails to converge computationally as the model is near +degenerate. We could try to get it to fit by working on the +computational algorithm. However, `ergm.tapered` considers a variant of +the ERGM that reflects our prior belief that the true generating process +is non-degenerate: + +``` r +fit <- ergm.tapered(samplike ~ edges + ttriple) +``` + + Starting maximum pseudolikelihood estimation (MPLE): + + ... + + Iteration 2 of at most 60: + Optimizing with step length 1.0000. + The log-likelihood improved by 0.0047. + Precision adequate twice. Stopping. + Finished MCMLE. + Evaluating log-likelihood at the estimate. Fitting the dyad-independent submodel... + Bridging between the dyad-independent submodel and the full model... + Setting up bridge sampling... + Using 16 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 . + Bridging finished. + This model was fit using MCMC. To examine model diagnostics and check + for degeneracy, use the mcmc.diagnostics() function. + +``` r +summary(fit) +``` + + Results: + + Estimate Std. Error MCMC % z value Pr(>|z|) + edges -1.83992 0.27737 0 -6.634 < 1e-04 *** + ttriple 0.21810 0.05816 0 3.750 0.000177 *** + --- + Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 + The estimated tapering scaling factor is 2. + + Null Deviance: 424.2 on 306 degrees of freedom + Residual Deviance: 358.0 on 304 degrees of freedom + + AIC: 362 BIC: 369.4 (Smaller is better. MC Std. Err. = 0) + +This tapered ERGM fits. The summary indicates that the coefficient on +the transitive triple term is positive (about 0.22) and statistically +above zero. This coefficient has the same interpretation as those in a standard ERGM. + +This model fixes the tapering parameter at 2 units. Let’s try to taper +less by increasing the tapering parameter to 3 (`r=3`): + +``` r +fit <- ergm.tapered(samplike ~ edges + ttriple, r=3) +``` + + Starting maximum pseudolikelihood estimation (MPLE): + + ... + + Iteration 2 of at most 60: + This model was fit using MCMC. To examine model diagnostics and check + for degeneracy, use the mcmc.diagnostics() function. + +``` r +summary(fit) +``` + + Results: + + Estimate Std. Error MCMC % z value Pr(>|z|) + edges -1.82890 0.27246 0 -6.712 <1e-04 *** + ttriple 0.21174 0.05255 0 4.029 <1e-04 *** + --- + Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 + The estimated tapering scaling factor is 3. + + Null Deviance: 424.2 on 306 degrees of freedom + Residual Deviance: 358.6 on 304 degrees of freedom + + AIC: 362.6 BIC: 370.1 (Smaller is better. MC Std. Err. = 0) + +It does not effect it much. + +The software allows the tapering to be estimated based on the shape of +the distributions of the model statistics. Let’s try that: + +``` r +fit <- ergm.tapered(samplike ~ edges + ttriple, fixed=FALSE) +``` + + Starting maximum pseudolikelihood estimation (MPLE): + + Iteration 4 of at most 60: + Optimizing with step length 1.0000. + The log-likelihood improved by 0.0003. + Precision adequate twice. Stopping. + Finished MCMLE. + +``` r +summary(fit) +``` + + Results: + + Estimate Std. Error MCMC % z value Pr(>|z|) + edges -1.83372 0.28426 0 -6.451 <1e-04 *** + ttriple 0.21236 0.09938 0 2.137 0.0326 * + --- + Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 + The estimated tapering scaling factor is 2.813. + + Null Deviance: 424.2 on 306 degrees of freedom + Residual Deviance: 359.0 on 304 degrees of freedom + + AIC: 365 BIC: 376.2 (Smaller is better. MC Std. Err. = 0) + +The estimated tapering parameter is about `2.8` (It is printed under the +coefficient table). This is between the default value and the second +guess. The coefficients of the ERGM terms are about the same as before. + +Enjoy trying `ergm.tapering`! + + + +See the following papers for more information and examples: + +#### Statistical Methodology + +* Fellows, Ian E. and Handcock, Mark S. (2017) [Removing Phase Transitions from Gibbs Measures](https://proceedings.mlr.press/v54/fellows17a/fellows17a.pdf), *Proceedings of the 20th International Conference on Artificial Intelligence and Statistics (AISTATS)s*, Volume 54. +* Blackburn, Bart and Handcock, Mark S. (2022) [Practical Network Modeling via Tapered Exponential-family Random Graph Models](https://doi.org/10.1080/10618600.2022.2116444), *Journal of Computational and Graphical Statistics*. + + + +## Public and Private repositories + +To facilitate open development of the package while giving the core developers an opportunity to publish on their developments before opening them up for general use, this project comprises two repositories: +* A public repository `statnet/ergm.tapered` +* A private repository `statnet/ergm.tapered-private` + +The intention is that all developments in `statnet/ergm.tapered-private` will eventually make their way into `statnet/ergm.tapered` and onto CRAN. + +Developers and Contributing Users to the Statnet Project should read https://statnet.github.io/private/ for information about the relationship between the public and the private repository and the workflows involved. + +## Latest Windows and MacOS binaries + +A set of binaries is built after every commit to the repository. We strongly encourage testing against them before filing a bug report, as they may contain fixes that have not yet been sent to CRAN. They can be downloaded through the following links: + +* [MacOS binary (a `.tgz` file in a `.zip` file)](https://nightly.link/statnet/ergm.tapered/workflows/R-CMD-check.yaml/master/macOS-rrelease-binaries.zip) +* [Windows binary (a `.zip` file in a `.zip` file)](https://nightly.link/statnet/ergm.tapered/workflows/R-CMD-check.yaml/master/Windows-rrelease-binaries.zip) + +You will need to extract the MacOS `.tgz` or the Windows `.zip` file from the outer `.zip` file before installing. These binaries are usually built under the latest version of R and their operating system and may not work under other versions. + +You may also want to install the corresponding latest binaries for packages on which `ergm.tapered` depends, in particular [`statnet.common`](https://github.com/statnet/statnet.common), [`network`](https://github.com/statnet/network), and [`ergm`](https://github.com/statnet/ergm). diff --git a/man/control.ergm.tapered.Rd b/man/control.ergm.tapered.Rd index a9e8db3..f0fb275 100644 --- a/man/control.ergm.tapered.Rd +++ b/man/control.ergm.tapered.Rd @@ -13,6 +13,7 @@ control.ergm.tapered( MCMLE.kurtosis.penalty = 2, MCMLE.termination = c("precision", "confidence", "Hummel", "Hotelling", "none"), MCMLE.MCMC.precision = 0.005, + estimate.tapered.bias = TRUE, loglik = control.logLik.ergm.tapered(), ... ) @@ -47,6 +48,9 @@ errors are below the precision bound, and the Hummel step length is 1 for two consecutive iterations. This is an experimental feature \code{\link[=ergm]{ergm()}}. The default value in \code{\link[=ergm.tapered]{ergm.tapered()}} is 0.15, higher than in \code{\link[=ergm]{ergm()}}.} +\item{estimate.tapered.bias}{logical. It \code{TRUE} the bias of the estimated estimates due to tapering is estimated. +If this is \code{FALSE} the MPLE of the untapered model is not computed, saving computational time.} + \item{loglik}{list List of additional control arguments for the loglik. See \code{\link{control.logLik.ergm.tapered}}.} \item{\dots}{Additional arguments, passed to other functions This argument diff --git a/man/ergm.tapered.Rd b/man/ergm.tapered.Rd index 6580a6c..afcb876 100644 --- a/man/ergm.tapered.Rd +++ b/man/ergm.tapered.Rd @@ -57,6 +57,16 @@ It can also be the RHS of a formula giving the terms to be tapered.} modeled in the valued ERGM framework. Defaults to \code{NULL} for simple presence or absence, modeled via a binary ERGM.}} +\item{constraints}{{A formula specifying one or more constraints +on the support of the distribution of the networks being modeled, +using syntax similar to the \code{formula} argument, on the +right-hand side. See \link[=ergm-constraints]{ERGM constraints}.}} + +\item{reference}{{A one-sided formula specifying +the reference measure (\eqn{h(y)}) to be used. +See help for \link[=ergm-references]{ERGM reference measures} implemented in the +\strong{\link[=ergm-package]{ergm}} package for details.}} + \item{estimate}{{If "MPLE," then the maximum pseudolikelihood estimator is returned. If "MLE" (the default), then an approximate maximum likelihood estimator is returned. For certain models, the MPLE and MLE are equivalent, @@ -66,7 +76,7 @@ likelihood calculation even when the MLE and MPLE are the same, see the the Monte-Carlo contrastive divergence estimate is returned. ) }} -\item{control}{An object of class control.ergm. Passed to the ergm function.} +\item{control}{An object of class control.erg.tapered. Passed to the ergm function.} \item{fixed}{A \code{logical}: if this is \code{TRUE}, the tapering is fixed at the passed values. If this is \code{FALSE}, the tapering is estimated using a kurtosis penalized likelihood. diff --git a/man/ergm.tapered.formula.Rd b/man/ergm.tapered.formula.Rd new file mode 100644 index 0000000..03e6ae4 --- /dev/null +++ b/man/ergm.tapered.formula.Rd @@ -0,0 +1,101 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ergm.tapered.formula.R +\name{ergm.tapered.formula} +\alias{ergm.tapered.formula} +\title{Create an Tapered ERGM formula specifying an Tapered ERGM model as a standard ERGM} +\usage{ +ergm.tapered.formula( + formula, + r = 2, + beta = NULL, + tau = NULL, + tapering.centers = NULL, + target.stats = NULL, + family = "taper", + taper.terms = "all", + response = NULL, + constraints = ~., + reference = ~Bernoulli, + verbose = FALSE, + ... +) +} +\arguments{ +\item{formula}{An ergm formula to fit} + +\item{r}{The scaling factor to use for the heuristic of setting beta equal to r standard deviations of the observed statistics} + +\item{beta}{The tapering parameters, expressed as in Fellows and Handcock (2017). If not NULL, these override the heuristics (r).} + +\item{tau}{The tapering parameters, expressed as natural parameters. If not NULL, these override the beta and the heuristics (r).} + +\item{tapering.centers}{The centers of the tapering terms. If null, these are taken to be the mean value parameters.} + +\item{target.stats}{{vector of "observed network statistics," +if these statistics are for some reason different than the +actual statistics of the network on the left-hand side of +\code{formula}. +Equivalently, this vector is the mean-value parameter values for the +model. If this is given, the algorithm finds the natural +parameter values corresponding to these mean-value parameters. +If \code{NULL}, the mean-value parameters used are the observed +statistics of the network in the formula. +}} + +\item{family}{The type of tapering used. This should either be the \code{stereo} or \code{taper}, the +tapering model of Fellows and Handcock (2017).} + +\item{taper.terms}{Specification of the tapering used. If the character variable "dependence" then all the dependent +terms are tapered. If the character variable "all" then all terms are tapered. +It can also be the RHS of a formula giving the terms to be tapered.} + +\item{response}{{Name of the edge attribute whose value is to be +modeled in the valued ERGM framework. Defaults to \code{NULL} for +simple presence or absence, modeled via a binary ERGM.}} + +\item{constraints}{{A formula specifying one or more constraints +on the support of the distribution of the networks being modeled, +using syntax similar to the \code{formula} argument, on the +right-hand side. See \link[=ergm-constraints]{ERGM constraints}.}} + +\item{reference}{{A one-sided formula specifying +the reference measure (\eqn{h(y)}) to be used. +See help for \link[=ergm-references]{ERGM reference measures} implemented in the +\strong{\link[=ergm-package]{ergm}} package for details.}} + +\item{verbose}{A \code{logical}: if this is +\code{TRUE}, the program will print out additional +information about the progress of estimation.} + +\item{...}{Additional arguments to \code{\link{ergm}}.} +} +\value{ +An formula containing a specification of the Tapered ERGM model. In addition to all of the ergm items, +this object contains tapering.centers, tapering.coef and orig.formula. tapering.centers are the centers for the tapering term. +tapering.coef are the tapering coefficients = 1/ beta^2. orig.formula is the formula passed into ergm.tapered. +} +\description{ +Create an Tapered ERGM formula specifying an Tapered ERGM model as a standard ERGM +} +\examples{ +\dontrun{ +data(sampson) +fit <- ergm.tapered.formula(samplike ~ edges + triangles()) +summary(fit) +} +} +\references{ +\itemize{ +\itemize{ +\item Fellows, I. and M. S. Handcock (2017), +Removing Phase Transitions from Gibbs Measures. Volume 54 of +Proceedings of Machine Learning Research, Fort Lauderdale, +FL, USA, pp. 289–297. PMLR. +\item Blackburn, B. and M. S. Handcock (2022), +Practical Network Modeling via Tapered Exponential-family Random Graph Models. +Journal of Computational and Graphical Statistics +\doi{10.1080/10618600.2022.2116444}. +} + +} +} diff --git a/man/simulate_ergm.tapered.Rd b/man/simulate_ergm.tapered.Rd new file mode 100644 index 0000000..152da20 --- /dev/null +++ b/man/simulate_ergm.tapered.Rd @@ -0,0 +1,103 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/simulate_ergm.tapered.R +\name{simulate_ergm.tapered} +\alias{simulate_ergm.tapered} +\title{Simulates from a Tapered ERGM formula} +\usage{ +simulate_ergm.tapered( + formula, + coef, + nsim = 1, + r = 2, + beta = NULL, + tau = NULL, + tapering.centers = NULL, + family = "taper", + taper.terms = "all", + response = NULL, + constraints = ~., + reference = ~Bernoulli, + seed = NULL, + control = control.simulate.formula(), + verbose = FALSE, + ... +) +} +\arguments{ +\item{formula}{An ergm formula to fit} + +\item{coef}{Vector of parameter values for the model from which the +sample is to be drawn.} + +\item{nsim}{Number of networks to be randomly drawn from the given +distribution on the set of all networks, returned by the Metropolis-Hastings +algorithm.} + +\item{r}{The scaling factor to use for the heuristic of setting beta equal to r standard deviations of the observed statistics} + +\item{beta}{The tapering parameters, expressed as in Fellows and Handcock (2017). If not NULL, these override the heuristics (r).} + +\item{tau}{The tapering parameters, expressed as natural parameters. If not NULL, these override the beta and the heuristics (r).} + +\item{family}{The type of tapering used. This should either be the \code{stereo} or \code{taper}, the +tapering model of Fellows and Handcock (2017).} + +\item{taper.terms}{Specification of the tapering used. If the character variable "dependence" then all the dependent +terms are tapered. If the character variable "all" then all terms are tapered. +It can also be the RHS of a formula giving the terms to be tapered.} + +\item{response}{{Name of the edge attribute whose value is to be +modeled in the valued ERGM framework. Defaults to \code{NULL} for +simple presence or absence, modeled via a binary ERGM.}} + +\item{constraints}{{A formula specifying one or more constraints +on the support of the distribution of the networks being modeled, +using syntax similar to the \code{formula} argument, on the +right-hand side. See \link[=ergm-constraints]{ERGM constraints} +\code{\link{ergm}} for details.}} + +\item{reference}{{A one-sided formula specifying +the reference measure (\eqn{h(y)}) to be used. +See help for \link[=ergm-references]{ERGM reference measures} implemented in the +\strong{\link[=ergm-package]{ergm}} package.}} + +\item{seed}{the random seed to use; see \code{\link[=simulate]{simulate()}}.} + +\item{control}{An object of class control.ergm. Passed to the ergm function.} + +\item{verbose}{A \code{logical}: if this is +\code{TRUE}, the program will print out additional +information about the progress of estimation.} + +\item{...}{Additional arguments to \code{\link{ergm}}.} +} +\value{ +An object of class c('tapered.ergm','ergm') containing the fit model. In addition to all of the ergm items, +this object contains tapering.centers, tapering.coef and orig.formula. tapering.centers are the centers for the tapering term. +tapering.coef are the tapering coefficients = 1/ beta^2. orig.formula is the formula passed into ergm.tapered. +} +\description{ +Simulates from a Tapered ERGM formula +} +\examples{ +\dontrun{ +data(sampson) +fit <- ergm.tapered(samplike ~ edges + triangles()) +summary(fit) +} +} +\references{ +\itemize{ +\itemize{ +\item Fellows, I. and M. S. Handcock (2017), +Removing Phase Transitions from Gibbs Measures. Volume 54 of +Proceedings of Machine Learning Research, Fort Lauderdale, +FL, USA, pp. 289–297. PMLR. +\item Blackburn, B. and M. S. Handcock (2022), +Practical Network Modeling via Tapered Exponential-family Random Graph Models. +Journal of Computational and Graphical Statistics +\doi{10.1080/10618600.2022.2116444}. +} + +} +}