fix vcov for models with offsets

NEWS.md

@@ -3,6 +3,7 @@ develop
 
 1. add Schoenfeld residual output for Cox models
 2. check for negative curvature before computing CIs
+3. fix `vcov` when model has an offset
 
 Cyclops v3.5.0
 ==============

R/ModelFit.R

@@ -1124,6 +1124,9 @@ vcov.cyclopsFit <- function(object, control, overrideNoRegularization = FALSE, .
     .checkInterface(object$cyclopsData, testOnly = TRUE)
     .setControl(object$cyclopsData$cyclopsInterfacePtr, control)
     fisherInformation <- .cyclopsGetFisherInformation(object$cyclopsData$cyclopsInterfacePtr, NULL)
+    if (.cyclopsGetHasOffset(object$cyclopsData)) {
+        fisherInformation <- fisherInformation[2:nrow(fisherInformation), 2:ncol(fisherInformation)]
+    }
     vcov <- solve(fisherInformation)
     if (!is.null(object$coefficientNames)) {
         rownames(vcov) <- object$coefficientNames

tests/testthat/test-smallCLR.R

@@ -1,5 +1,6 @@
 library("testthat")
 library("survival")
+library("gnm")
 
 context("test-smallCLR.R")
 
@@ -42,5 +43,21 @@ test_that("Small conditional logistic regression", {
 # 	expect_equal(predict(cyclopsFit), predict(gold, type = "risk"), tolerance = tolerance)
 })
 
+test_that("Small conditional poisson regression with an offset", {
+
+    gold.cp <- gnm(event ~ exgr + agegr + offset(loginterval),
+                   family = poisson, eliminate = indiv,
+                   data = Cyclops::oxford)
+    vcov(gold.cp)
+
+    dataPtr <- createCyclopsData(event ~ exgr + agegr + strata(indiv) + offset(loginterval),
+                                 data = Cyclops::oxford,
+                                 modelType = "cpr")
+    cyclopsFit <- fitCyclopsModel(dataPtr,
+                                  prior = createPrior("none"))
+
+    expect_equal(coef(cyclopsFit), vcov(gold.cp))
+    expect_equal(vcov(cyclopsFit), vcov(gold.cp))
+})