Handle SE calculation for GAM fits

Author: baddstats

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: spatstat.model
-Version: 3.2-1.006
-Date: 2023-03-15
+Version: 3.2-1.007
+Date: 2023-03-16
 Title: Parametric Statistical Modelling and Inference for the 'spatstat' Family
 Authors@R: c(person("Adrian", "Baddeley", 
                     role = c("aut", "cre", "cph"),

NEWS

@@ -1,5 +1,5 @@
 
-        CHANGES IN spatstat.model VERSION 3.2-1.006
+        CHANGES IN spatstat.model VERSION 3.2-1.007
 
 OVERVIEW
 
@@ -9,6 +9,20 @@ OVERVIEW
 
     o Minor changes to documentation.
 
+    o Bug fixes.
+
+BUG FIXES
+
+    o predict.ppm
+    Standard error calculation ('se=TRUE') crashed if the fitted model
+    was a generalised additive model (fitted with 'use.gam=TRUE').
+    Fixed.
+
+    o effectfun
+    Standard error calculation ('se.fit=TRUE') crashed if the fitted model
+    was a generalised additive model (fitted with 'use.gam=TRUE').
+    Fixed.
+    
         CHANGES IN spatstat.model VERSION 3.2-1
 
 OVERVIEW

R/predict.ppm.R

@@ -490,30 +490,43 @@ predict.ppm <- local({
       if(type == "intensity") 
         z <- PoisSaddle(z, fitin(model))
 
-      ##
       if(needSE) {
-        ## extract variance-covariance matrix of parameters
-        vc <- vcov(model, new.coef=new.coef)
-        ## compute model matrix
-        fmla <- rhs.of.formula(formula(glmfit))
+        ## compute standard error, needed for some purpose
+        if(inherits(glmfit, "gam")) {
+          ## compute SE using predict.gam
+          if(!is.null(new.coef))
+            warning("new.coef ignored in standard error calculation")
+          SE <- predict(glmfit, newdata=newdata, type="response",
+                        se.fit=TRUE)[[2]]
+        } else {
+          ## Use vcov.ppm
+          ## extract variance-covariance matrix of parameters
+          vc <- vcov(model, new.coef=new.coef)
+          ## compute model matrix
+          fmla <- rhs.of.formula(formula(glmfit))
 #        mf <- model.frame(fmla, newdata, ..., na.action=na.pass)
 #        mm <- model.matrix(fmla, mf, ..., na.action=na.pass)
-        mf <- model.frame(fmla, newdata, na.action=na.pass)
-        mm <- model.matrix(fmla, mf, na.action=na.pass)
-        if(nrow(mm) != nrow(newdata))
-          stop("Internal error: row mismatch in SE calculation")
-        ## compute relative variance = diagonal of quadratic form
-        if(ncol(mm) != ncol(vc))
-          stop("Internal error: column mismatch in SE calculation")
-        vv <- quadform(mm, vc)
-        ## standard error
-        SE <- lambda * sqrt(vv)
+          mf <- model.frame(fmla, newdata, na.action=na.pass)
+          mm <- model.matrix(fmla, mf, na.action=na.pass)
+          if(nrow(mm) != nrow(newdata))
+            stop("Internal error: row mismatch in SE calculation")
+          ## compute relative variance = diagonal of quadratic form
+          if(ncol(mm) != ncol(vc))
+            stop("Internal error: column mismatch in SE calculation")
+          vv <- quadform(mm, vc)
+          ## standard error
+          SE <- lambda * sqrt(vv)
+        }
+
+        ## Save desired quantities depending on SE
         if(se) 
           zse <- SE
+
         if(interval == "confidence") {
           z <- lambda + outer(SE, ci.q, "*")
           colnames(z) <- ci.names
-        } 
+        }
+        
       } 
 
       ## ############################################################  

inst/doc/packagesizes.txt

@@ -4,4 +4,4 @@ date version nhelpfiles nobjects ndatasets Rlines srclines
 "2023-01-26" "3.1-2" 267 706 0 38848 1155
 "2023-02-14" "3.2-0" 269 709 0 38930 1155
 "2023-02-15" "3.2-1" 269 709 0 38930 1155
-"2023-03-15" "3.2-1.006" 269 709 0 38931 1155
+"2023-03-16" "3.2-1.007" 269 709 0 38944 1155