sginv() calls to compute covariance matrices now use tol=0 to avoid u…

…nderestimating standard errors of parameters corresponding to highly correlated statistics.
statnet · Dec 9, 2023 · f6b2f2d · f6b2f2d
1 parent e645c8a
commit f6b2f2d
Show file tree

Hide file tree

Showing 8 changed files with 10 additions and 10 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: ergm
-Version: 4.6-7275
-Date: 2023-11-29
+Version: 4.6-7276
+Date: 2023-12-09
 Title: Fit, Simulate and Diagnose Exponential-Family Models for Networks
 Authors@R: c(
              person(c("Mark", "S."), "Handcock", role=c("aut"), email="handcock@stat.ucla.edu"),

diff --git a/R/approx.hotelling.diff.test.R b/R/approx.hotelling.diff.test.R
@@ -166,7 +166,7 @@ approx.hotelling.diff.test<-function(x,y=NULL, mu0=0, assume.indep=FALSE, var.eq
 
   if(p==0) stop("data are essentially constant")
 
-  ivcov.d <-sginv(vcov.d[!novar,!novar,drop=FALSE])
+  ivcov.d <-sginv(vcov.d[!novar,!novar,drop=FALSE], tol=0)
 
   method <- paste0("Hotelling's ",
                   NVL2(y, "Two", "One"),

diff --git a/R/ergm.MCMLE.R b/R/ergm.MCMLE.R
@@ -321,7 +321,7 @@ ergm.MCMLE <- function(init, s, s.obs,
       # if some statistic has a variance of exactly 0.
       novar <- diag(Vm) == 0
       Vm[!novar,!novar] <- snearPD(Vm[!novar,!novar,drop=FALSE], posd.tol=0, base.matrix=TRUE)$mat
-      iVm <- sginv(Vm)
+      iVm <- sginv(Vm, tol=0)
       diag(Vm)[novar] <- sqrt(.Machine$double.xmax) # Virtually any nonzero difference in estimating functions will map to a very large number.
       d2 <- xTAx(estdiff, iVm)
       if(d2<2) last.adequate <- TRUE

diff --git a/R/ergm.estimate.R b/R/ergm.estimate.R
@@ -225,7 +225,7 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
     # or more statistics.
     if(inherits(Lout$par,"try-error")){
       Lout$par <- try(eta0 
-                      + sginv(-Lout$hessian) %*%
+                      + sginv(-Lout$hessian, tol=0) %*%
                       xobs,
                       silent=TRUE)
     }
@@ -339,7 +339,7 @@ ergm.estimate<-function(init, model, statsmatrices, statsmatrices.obs=NULL,
     }
 
     covar <- matrix(NA, ncol=length(theta), nrow=length(theta))
-    covar[!model$etamap$offsettheta,!model$etamap$offsettheta ] <- sginv(-Lout$hessian)
+    covar[!model$etamap$offsettheta,!model$etamap$offsettheta ] <- sginv(-Lout$hessian, tol=0)
     dimnames(covar) <- list(names(theta),names(theta))
     He <- matrix(NA, ncol=length(theta), nrow=length(theta))
     He[!model$etamap$offsettheta,!model$etamap$offsettheta ] <- Lout$hessian

diff --git a/R/ergm.logitreg.R b/R/ergm.logitreg.R
@@ -137,7 +137,7 @@ ergm.logitreg <- function(x, y, wt = rep(1, length(y)),
   names(fit$coefficients) <- dn[if(!is.null(m)) !m$etamap$offsettheta else TRUE]
   fit$deviance <- -2*fit$value
   fit$iter <- fit$iterations
-  asycov <- sginv(-fit$hessian)
+  asycov <- sginv(-fit$hessian, tol=0)
   fit$cov.unscaled <- asycov
   if(!fit$converged)
       message("Trust region algorithm did not converge.")

diff --git a/R/ergm.mple.R b/R/ergm.mple.R
@@ -160,7 +160,7 @@ ergm.mple<-function(s, s.obs, init=NULL,
   covar[!is.na(theta)&!m$etamap$offsettheta,
         !is.na(theta)&!m$etamap$offsettheta] <- real.cov
   hess[!is.na(theta)&!m$etamap$offsettheta,
-        !is.na(theta)&!m$etamap$offsettheta] <- if(length(real.cov)) -sginv(real.cov) else matrix(0,0,0)
+        !is.na(theta)&!m$etamap$offsettheta] <- if(length(real.cov)) -sginv(real.cov, tol=0) else matrix(0,0,0)
 #
   iteration <-  mplefit$iter 
 

diff --git a/R/ergm.pen.glm.R b/R/ergm.pen.glm.R
@@ -101,7 +101,7 @@ ergm.pen.glm <- function(formula,
    iter <- iter + 1
    XW2 <- sweep(x, 1, (weights*(pi * (1 - pi)))^0.5, "*") #### X' (W ^ 1/2)
    Fisher <- crossprod(XW2)  #### X' W  X
-   covs <- sginv(Fisher)  ### (X' W  X) ^ -1
+   covs <- sginv(Fisher, tol=0)  ### (X' W  X) ^ -1
 #  H <- crossprod(XW2, covs) %*% XW2
 #  H <- XW2 %*% covs %*% t(XW2)
    diagH <- pi

diff --git a/R/vcov.ergm.R b/R/vcov.ergm.R
@@ -39,7 +39,7 @@ vcov.ergm <- function(object, sources=c("all","model","estimation"), ...){
     if(is.null(object$hessian) && is.null(object$covar)){
       object$covar <- matrix(NA, p, p)
     }
-    v.mod <- NVL(object$covar, sginv(-object$hessian))
+    v.mod <- NVL(object$covar, sginv(-object$hessian, tol=0))
     v.mod[is.na(diag(v.mod))|diag(v.mod)<0|is.infinite(coef(object)),] <- NA
     v.mod[,is.na(diag(v.mod))|diag(v.mod)<0|is.infinite(coef(object))] <- NA
     v.mod[object$offset,] <- 0