From f6b2f2d582a4d7e220c81d3d1190e14e4990c952 Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Sat, 9 Dec 2023 12:17:30 +1100 Subject: [PATCH] sginv() calls to compute covariance matrices now use tol=0 to avoid underestimating standard errors of parameters corresponding to highly correlated statistics. --- DESCRIPTION | 4 ++-- R/approx.hotelling.diff.test.R | 2 +- R/ergm.MCMLE.R | 2 +- R/ergm.estimate.R | 4 ++-- R/ergm.logitreg.R | 2 +- R/ergm.mple.R | 2 +- R/ergm.pen.glm.R | 2 +- R/vcov.ergm.R | 2 +- 8 files changed, 10 insertions(+), 10 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 07e12b44f..01dce0e0c 100644 --- 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 index 0c215846e..072fa3e3c 100644 --- 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 index 8e6019175..e680a1537 100644 --- 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 index 53ab3ddaf..6a5b2110f 100644 --- 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 index cc2ca4a15..a8df4cf54 100644 --- 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 index 206719e8a..3a02ef3cf 100644 --- 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 index 94dcdc0f8..0a339de1d 100644 --- 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 index aa760f6b5..78738c727 100644 --- 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