checks sandwich estimator

Author: jchiquet

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: PLNmodels
 Title: Poisson Lognormal Models
-Version: 1.2.1
+Version: 1.2.2
 Authors@R: c(
     person("Julien", "Chiquet", , "julien.chiquet@inrae.fr", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-3629-3429")),

NEWS.md

@@ -1,5 +1,6 @@
 # Current version
 
+* fix sandwich variance estimation (PR #140)
 * fix use of native pipe to ensure compatibility with R 3.6 (merge PR #125, fix #124)
 
 # PLNmodels 1.2.0 (2024-03-05)

R/PLNfit-S3methods.R

@@ -156,7 +156,7 @@ sigma.PLNfit <- function(object, ...) {
 #' @description Extracts univariate standard errors for the estimated coefficient of B. Standard errors are computed from the (approximate) Fisher information matrix.
 #'
 #' @param object an R6 object with class PLNfit
-#' @param type string describing the type of variance approximation: "variational", "jackknife", "sandwich" (only for fixed covariance). Default is "variational".
+#' @param type string describing the type of variance approximation: "variational", "jackknife", "sandwich". Default is "sandwich".
 #' @param parameter string describing the target parameter: either B (regression coefficients) or Omega (inverse residual covariance)
 #'
 #' @seealso [vcov.PLNfit()] for the complete variance covariance estimation of the coefficient
@@ -169,13 +169,13 @@ sigma.PLNfit <- function(object, ...) {
 #'               control = PLN_param(config_post = list(variational_var = TRUE)))
 #' standard_error(myPLN)
 #' @export
-standard_error <- function(object, type = c("variational", "jackknife", "sandwich"), parameter = c("B", "Omega")) {
+standard_error <- function(object, type = c("sandwich", "variational", "jackknife"), parameter = c("B", "Omega")) {
   UseMethod("standard_error", object)
 }
 
 #' @describeIn standard_error Component-wise standard errors of B in [`PLNfit`]
 #' @export
-standard_error.PLNfit <- function(object, type = c("variational", "sandwich", "jackknife", "bootstrap"), parameter = c("B", "Omega")) {
+standard_error.PLNfit <- function(object, type = c("sandwich", "variational", "jackknife", "bootstrap"), parameter = c("B", "Omega")) {
   type <- match.arg(type)
   par  <- match.arg(parameter)
   if (type == "variational" & is.null(attr(object$model_par$B, "variance_variational")))
@@ -191,7 +191,7 @@ standard_error.PLNfit <- function(object, type = c("variational", "sandwich", "j
 
 #' @describeIn standard_error Component-wise standard errors of B in [`PLNfit_fixedcov`]
 #' @export
-standard_error.PLNfit_fixedcov <- function(object, type = c("variational", "sandwich", "jackknife", "bootstrap", "sandwich"), parameter = c("B", "Omega")) {
+standard_error.PLNfit_fixedcov <- function(object, type = c("sandwich", "variational", "jackknife", "bootstrap"), parameter = c("B", "Omega")) {
   type <- match.arg(type)
   par  <- match.arg(parameter)
   if (par == "Omega")

man/standard_error.Rd

@@ -4,7 +4,7 @@ data(trichoptera)
 trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
 
 ## Error messages ---------------------
-test_that("Check that fisher and standard_error return objects with proper dimensions and sign",  {
+test_that("Variational Fisher: Check that fisher and standard_error return objects with proper dimensions and sign",  {
 
   myPLN_cov <- PLN(Abundance ~ Wind + offset(log(Offset)), data = trichoptera, control = PLN_param(config_post = list(variational_var = TRUE)))
   expect_is(myPLN_cov, "PLNfit")
@@ -31,20 +31,19 @@ test_that("Check that fisher and standard_error return objects with proper dimen
 myPLN <- PLN(Abundance ~ 1 + offset(log(Offset)), data = trichoptera, control = PLN_param(config_post = list(variational_var = TRUE)))
 
 ## Consistency -----------------------
-test_that("Check internal consistency of Fisher matrix for PLN models with no covariates",  {
+test_that("Variational Fisher: Check internal consistency of Fisher matrix for PLN models with no covariates",  {
   tol <- 1e-8
 
   n <- nrow(myPLN$fitted)
   ## Consistency of the standard error matrix
-  sem <- (sqrt(n) * standard_error(myPLN)) %>% as.numeric()
+  sem <- (sqrt(n) * standard_error(myPLN, "variational")) %>% as.numeric()
   manual.sem <- 1/colMeans(myPLN$fitted) %>% sqrt()
 
   ## Internal consistency
   expect_equal(sem, manual.sem, tolerance = tol)
 
 })
 
-
 test_that("Check temporal consistency of Fisher matrix for PLN models with no covariates",  {
   tol <- 1e-2
 
@@ -73,6 +72,31 @@ test_that("Check temporal consistency of Fisher matrix for PLN models with no co
 
 })
 
+
+## Error messages ---------------------
+test_that("Sandwich: check that fisher and standard_error return objects with proper dimensions and sign",  {
+
+  myPLN_cov <- PLN(Abundance ~ Wind + offset(log(Offset)), data = trichoptera, control = PLN_param(config_post = list(sandwich_var = TRUE)))
+  expect_is(myPLN_cov, "PLNfit")
+  p <- myPLN_cov$p
+  d <- myPLN_cov$d
+
+
+  sem <- standard_error(myPLN_cov)
+  ## Dimensions
+  expect_equal(dim(sem), c(d, p))
+
+  ## Names
+  expect_equal(rownames(sem), rownames(coef(myPLN_cov)))
+  expect_equal(colnames(sem), colnames(coef(myPLN_cov)))
+
+  ## Standard errors are all positive
+  for (i in 1:(p*d)) {
+    expect_gte(sem[i], 0)
+  }
+
+})
+
 params <- PLNmodels:::create_parameters(n = 30, p = 3, d = 1, depths = 1e3)
 B <- params$B
 X <- as.matrix(params$X)