-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Add target.stats capability and control$estimate.tapered.bias
Add the ability to use target.stats and not double fit the MPLE when estimate="MPLE"
- Loading branch information
Showing
12 changed files
with
954 additions
and
19 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -95,3 +95,4 @@ dkms.conf | |
| .Rproj.user | ||
| R/ergm.tapered.R.withcomments | ||
| R/ergm.tapered.R.nosancall | ||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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) | ||
| } |
Oops, something went wrong.