fix standard errors predict + newdata

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: gslnls
 Type: Package
 Title: GSL Multi-Start Nonlinear Least-Squares Fitting
-Version: 1.3.2
-Date: 2024-04-30
+Version: 1.3.3
+Date: 2024-05-15
 Authors@R: person("Joris", "Chau", email = "[email protected]", role = c("aut", "cre"))
 Description: An R interface to nonlinear least-squares optimization with the GNU Scientific Library (GSL), see M. Galassi et al. (2009, ISBN:0954612078). The available trust region methods include the Levenberg-Marquardt algorithm with and without geodesic acceleration, the Steihaug-Toint conjugate gradient algorithm for large systems and several variants of Powell's dogleg algorithm. Multi-start optimization based on quasi-random samples is implemented using a modified version of the algorithm in Hickernell and Yuan (1997, OR Transactions). Bindings are provided to tune a number of parameters affecting the low-level aspects of the trust region algorithms. The interface mimics R's nls() function and returns model objects inheriting from the same class.
 BugReports: https://github.com/JorisChau/gslnls/issues

NEWS.md

@@ -1,3 +1,7 @@
+# gslnls 1.3.3
+
+* Fix standard errors `predict()` when using `newdata`
+
 # gslnls 1.3.2
 
 * Reverted to static Makevars.win (supplied by T. Kalibera)

R/nls_methods.R

@@ -272,17 +272,16 @@ print.gsl_nls <- function(x, digits = max(3L, getOption("digits") - 3L), ...) {
 #' @export
 predict.gsl_nls <- function(object, newdata, scale = NULL, interval = c("none", "confidence", "prediction"), level = 0.95, ...) {
   interval <- match.arg(interval, c("none", "confidence", "prediction"))
+
   if (missing(newdata)) {
     fit <- as.vector(fitted(object))
     if(interval != "none") {
       Fdot <- object$m$gradient()
-      Rmat <- object$m$Rmat()
     }
   } else {
     fit <- object$m$predict(newdata)
     if(interval != "none") {
       Fdot <- object$m$gradient1(newdata)
-      Rmat <- qr.R(qr(Fdot))
       if(!is.null(object$weights))
         warning("unweighted Jacobian matrix used to calculate standard errors, evaluate predictions without 'newdata' argument to use weighted Jacobian.")
     }
@@ -291,7 +290,7 @@ predict.gsl_nls <- function(object, newdata, scale = NULL, interval = c("none",
     if(is.null(scale))
       scale <- sigma(object)
     a <- c((1 - level) / 2, (1 + level) / 2)
-    ses <- scale * sqrt(1 * (interval == "prediction") + rowSums(Fdot %*% chol2inv(Rmat) * Fdot))
+    ses <- scale * sqrt(1 * (interval == "prediction") + rowSums(Fdot %*% chol2inv(object$m$Rmat()) * Fdot))
     ci <- fit + ses %o% qt(a, df.residual(object))
     cimat <- cbind(fit = fit, lwr = ci[, 1], upr = ci[, 2])
     return(cimat)