Added an estimate of the bias in the coefficients introduced by tapering

Added an estimate of the bias in the coefficients introduced by tapering.

The standard ERGM coefficients have a conditional log-odds interpretation. Tapered coefficients do not have this interpretation, but are close to the conditional log-odds interpretation. This "bias" term measures the average difference between the coefficient and the conditional log-odds for the tapered ERGM model. That is, positive bias indicates the coefficient is above the conditional log-odds.

DESCRIPTION

@@ -22,4 +22,4 @@ Imports:
 License: GPL-3
 Roxygen: list(markdown = TRUE)
 Encoding: UTF-8
-RoxygenNote: 7.1.1
+RoxygenNote: 7.1.2

R/ergm.tapered.R

@@ -217,6 +217,10 @@ 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 <- ergm(formula, 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, ...)
@@ -265,6 +269,10 @@ ergm.tapered <- function(formula, r=2, beta=NULL, tau=NULL, tapering.centers=NUL
   fit$tapering.coef <- tau
   fit$r <- r
   fit$orig.formula <- formula
+
+  fit$taudelta.mean <- apply((2*fit.MPLE$glm$model[,1]-1)*sweep(fit.MPLE$xmat.full,2,fit$tapering.coef,"*"),2,weighted.mean,weight=fit.MPLE$glm$prior.weights)
+  fit$taudelta.median <- apply((2*fit.MPLE$glm$model[,1]-1)*sweep(fit.MPLE$xmat.full,2,fit$tapering.coef,"*"),2,wtd.median,weight=fit.MPLE$glm$prior.weights)
+
   class(fit) <- c("ergm.tapered",family,class(fit))
 
 # fit$covar <- fit$est.cov

R/print.summary.ergm.tapered.R

@@ -36,7 +36,7 @@
 #' @export
 print.summary.ergm.tapered <- function (x, 
               digits = max(3, getOption("digits") - 3),
-              correlation=x$correlation, covariance=x$covariance,
+              correlation=x$correlation, covariance=x$covariance, extended=x$extended,
               signif.stars= getOption("show.signif.stars"),
               eps.Pvalue=0.0001, print.formula=FALSE, print.fitinfo=TRUE, print.coefmat=TRUE, print.message=TRUE, print.deviances=TRUE, print.drop=TRUE, print.offset=TRUE, print.call=TRUE,...){
 
@@ -56,9 +56,16 @@ print.summary.ergm.tapered <- function (x,
   }
 
   if(print.coefmat){
-    printCoefmat(coef(x), digits=digits, signif.stars=signif.stars,
+    if(extended == TRUE){
+     coef.x <- cbind(coef(x)[,1:2], tau=x$tapering.coef, bias=as.numeric(sprintf("%.4f", -x$taudelta.mean)), coef(x)[,4:5])
+    #digits <- c(3,3,3,3,3,3)
+    }else{
+     coef.x <- coef(x)
+    #digits <- c(3,3,3,3)
+    }
+    printCoefmat(coef.x, digits=digits, signif.stars=signif.stars,
                  P.values=TRUE, has.Pvalue=TRUE, na.print="NA",
-                 eps.Pvalue=eps.Pvalue, cs.ind=1:2, tst.ind=4L,...)
+                 eps.Pvalue=eps.Pvalue, cs.ind=c(1:2,4), tst.ind=5L,...)
   }
 
   writeLines(c(strwrap(paste0("The estimated tapering scaling factor is ", format(x$r, digits = digits),".")),''))

R/summary.ergm.tapered.R

@@ -17,6 +17,9 @@
 #'   the likelihood variation and the MCMC variation. If \code{FALSE}
 #'   only the likelihood varuation is used. The \eqn{p}-values are
 #'   based on this source of variation.
+#' @param extended logical; if \code{TRUE}, the natural parameters of the tapered terms
+#' are printed as an appended column of the summary table. Also printed are the bias
+#' in the parameter values due to the tapering
 #' @param \dots Arguments to \code{\link{logLik.ergm}}
 #' @return The function \code{\link{summary.ergm.tapered}} computes and
 #'   returns a list of summary statistics of the fitted
@@ -40,7 +43,7 @@
 #' @export
 summary.ergm.tapered <- function (object, ..., 
                           correlation=FALSE, covariance=FALSE,
-                          total.variation=TRUE)
+                          total.variation=TRUE, extended=FALSE)
 {
   if("digits" %in% names(list(...))) warning("summary.ergm.tapered() no longer takes a digits= argument.")
 
@@ -61,13 +64,16 @@ summary.ergm.tapered <- function (object, ...,
 
   ans <- list(formula=object$formula,
               correlation=correlation,
+              extended = extended,
               degeneracy.value = object$degeneracy.value,
               offset = object$offset[cnames],
               drop = NVL(object$drop[cnames], rep(0,length(object$offset[cnames]))),
               estimable = NVL(object$estimable[cnames], rep(TRUE,length(object$offset[cnames]))),
               covariance=covariance,
               pseudolikelihood=pseudolikelihood,
               independence=independence,
+              tapering.coef=object$tapering.coef,
+              taudelta.mean=object$taudelta.mean,
               estimate=object$estimate,
               control=object$control)
 
@@ -191,6 +197,7 @@ summary.ergm.tapered <- function (object, ...,
     a[is.nan(a) | a < 0] <- 0
     ans$r <- mean(sqrt(a))
   }
+
   class(ans) <- "summary.ergm.tapered"
   ans
 }