Running examples in ‘GeDS-Ex.R’ failed
The error most likely occurred in:

> ### Name: GGeDS
> ### Title: Generalized Geometrically Designed Spline regression estimation
> ### Aliases: GGeDS
> 
> ### ** Examples
> 
> ######################################################################
> # Generate a data sample for the response variable Y and the covariate X
> # assuming Poisson distributed error and log link function
> # See section 4.1 in Dimitrova et al. (2023)
> set.seed(123)
> N <- 500
> f_1 <- function(x) (10*x/(1+100*x^2))*4+4
> X <- sort(runif(N, min = -2, max = 2))
> # Specify a model for the mean of Y to include only a component
> # non-linear in X, defined by the function f_1
> means <- exp(f_1(X))
> 
> #############
> ## POISSON ##
> #############
> # Generate Poisson distributed Y according to the mean model
> Y <- rpois(N, means)
> 
> # Fit a Poisson GeDS regression using GGeDS
> (Gmod <- GGeDS(Y ~ f(X), beta = 0.2, phi = 0.99, q = 2, family = poisson(),
+                Xextr = c(-2,2)))

Call:
GGeDS(formula = Y ~ f(X), family = poisson(), beta = 0.2, phi = 0.99, 
    q = 2, Xextr = c(-2, 2))

Number of internal knots of the second order spline: 13
Deviances:
Order 2  Order 3  Order 4  
  490.6    483.1    511.7  

> # Plot the quadratic and cubic GeDS fits
> plot(X, log(Y), xlab = "x", ylab = expression(f[1](x)))
> lines(X, sapply(X, f_1), lwd = 2)
> lines(Gmod, n = 3, col = "red")
> lines(Gmod, n = 4, col = "blue", lty = 2)
> legend("topleft",
+        legend = expression(f[1](x), "Quadratic", "Cubic"),
+        col = c("black", "red", "blue"),
+        lty = c(1, 1, 2),
+        lwd = c(2, 1, 1),
+        bty = "n")
> 
> # Generate GeDS prediction at X=0, first on the response scale and then on
> # the predictor scale
> predict(Gmod, n = 3, newdata = data.frame(X = 0))
         [,1]
[1,] 434.4222
> predict(Gmod, n = 3, newdata = data.frame(X = 0), type = "link")
         [,1]
[1,] 6.074017
> 
> # Apply some of the other available methods, e.g.
> # knots, coefficients and deviance extractions for the
> # quadratic GeDS fit
> knots(Gmod)
 [1] -2.000000000 -2.000000000 -2.000000000 -1.524015922 -1.085910542
 [6] -0.699338583 -0.278677095 -0.077787009 -0.007524582  0.040184661
[11]  0.078467796  0.132250944  0.189800300  0.356034439  0.605776373
[16]  2.000000000  2.000000000  2.000000000
> coef(Gmod)
      N1       N2       N3       N4       N5       N6       N7       N8 
3.831712 3.738951 3.682653 3.541883 3.433803 2.273420 2.043467 4.805209 
      N9      N10      N11      N12      N13      N14      N15 
5.849636 6.074017 5.728570 5.284927 4.700625 4.185276 4.261842 
> deviance(Gmod)
[1] 483.0868
> 
> # the same but for the cubic GeDS fit
> knots(Gmod, n = 4)
 [1] -2.00000000 -2.00000000 -2.00000000 -2.00000000 -1.33555637 -0.87062032
 [7] -0.50533049 -0.19853797 -0.04412116  0.01403654  0.06004872  0.10722022
[13]  0.15979247  0.29226465  0.47547609  2.00000000  2.00000000  2.00000000
[19]  2.00000000
> coef(Gmod, n = 4)
      N1       N2       N3       N4       N5       N6       N7       N8 
3.818281 3.767209 3.674072 3.588658 3.384239 2.817902 1.356245 4.366175 
      N9      N10      N11      N12      N13      N14      N15 
6.028346 6.029461 5.634932 5.071845 4.089328 4.391411 4.214580 
> deviance(Gmod, n = 4)
[1] 511.7263
> 
> ###########
> ## GAMMA ##
> ###########
> # Generate Gamma distributed Y according to the mean model
> Y <- rgamma(N, shape = means, rate = 0.1)
> # Fit a Gamma GeDS regression using GGeDS
> Gmod <- GGeDS(Y ~ f(X), beta = 0.1, phi = 0.99, q = 2, family = Gamma(log),
+               Xextr = c(-2,2))
> plot(Gmod, f = function(x) exp(f_1(x))/0.1)
> 
> ##############
> ## BINOMIAL ##
> ##############
> # Generate Binomial distributed Y according to the mean model
> eta <- f_1(X) - 4
> means <- exp(eta)/(1+exp(eta))
> Y <- rbinom(N, size = 50, prob = means) / 50
> # Fit a Binomial GeDS regression using GGeDS
> Gmod <- GGeDS(Y ~ f(X), beta = 0.1, phi = 0.99, family = "quasibinomial",
+               Xextr = c(-2,2))
> plot(Gmod, f = function(x) exp(f_1(x) - 4)/(1 + exp(f_1(x) - 4)))
> 
> 
> ##########################################
> # A real data example
> # See Dimitrova et al. (2023), Section 4.2
> 
> data("coalMining")
> (Gmod2 <- GGeDS(formula = accidents ~ f(years), beta = 0.1, phi = 0.98,
+                  family = poisson(), data = coalMining))

Call:
GGeDS(formula = accidents ~ f(years), family = poisson(), data = coalMining, 
    beta = 0.1, phi = 0.98)

Number of internal knots of the second order spline: 5
Deviances:
Order 2  Order 3  Order 4  
  117.2    118.0    119.7  

> (Gmod3 <- GGeDS(formula = accidents ~ f(years), beta = 0.1, phi = 0.985,
+                  family = poisson(), data = coalMining))

Call:
GGeDS(formula = accidents ~ f(years), family = poisson(), data = coalMining, 
    beta = 0.1, phi = 0.985)

Number of internal knots of the second order spline: 14
Deviances:
Order 2  Order 3  Order 4  
   99.4    102.8    105.8  

> plot(coalMining$years, coalMining$accidents, type = "h", xlab = "Years",
+      ylab = "Accidents")
> lines(Gmod2, tr = exp, n = 4, col = "red")
> lines(Gmod3, tr = exp, n = 4, col = "blue", lty = 2)
> legend("topright", c("phi = 0.98","phi = 0.985"), col = c("red", "blue"),
+        lty=c(1, 2))
> 
> 
> ## Not run: 
> ##D ##########################################
> ##D # The same regression in the example of GeDS
> ##D # but assuming Gamma and Poisson responses
> ##D # See Dimitrova et al. (2023), Section 4.2
> ##D 
> ##D data('BaFe2As2')
> ##D (Gmod4 <- GGeDS(intensity ~ f(angle), data = BaFe2As2, beta = 0.6, phi = 0.995, q = 3,
> ##D                 family = Gamma(log), stoptype = "RD"))
> ##D plot(Gmod4)
> ##D 
> ##D (Gmod5 <- GGeDS(intensity ~ f(angle), data = BaFe2As2, beta = 0.1, phi = 0.995, q = 3,
> ##D                 family = poisson(), stoptype = "SR"))
> ##D plot(Gmod5)
> ## End(Not run)
> 
> ##########################################
> # Life tables
> # See Dimitrova et al. (2023), Section 4.2
> 
> data(EWmortality)
> attach(EWmortality)
> (M1 <- GGeDS(formula = Deaths ~ f(Age) + offset(log(Exposure)),
+               family = quasipoisson(), phi = 0.99, beta = 0.1, q = 3,
+               stoptype = "LR"))

Call:
GGeDS(formula = Deaths ~ f(Age) + offset(log(Exposure)), family = quasipoisson(), 
    beta = 0.1, phi = 0.99, q = 3, stoptype = "LR")

Number of internal knots of the second order spline: 24
Deviances:
Order 2  Order 3  Order 4  
  644.1    805.7    802.3  

> 
> Exposure_init <- Exposure + 0.5 * Deaths
> Rate <- Deaths / Exposure_init
> (M2 <- GGeDS(formula = Rate ~ f(Age), weights = Exposure_init,
+               family = quasibinomial(), phi = 0.99, beta = 0.1,
+               q = 3, stoptype = "LR"))

Call:
GGeDS(formula = Rate ~ f(Age), family = quasibinomial(), weights = Exposure_init, 
    beta = 0.1, phi = 0.99, q = 3, stoptype = "LR")

Number of internal knots of the second order spline: 20
Deviances:
Order 2  Order 3  Order 4  
  739.8    937.7    724.3  

> 
> 
> op <- par(mfrow=c(2,2))
> plot(Age, Deaths/Exposure, ylab = expression(mu[x]), xlab = "Age")
> lines(M1, n = 3, tr = exp, lwd = 1, col = "red")
> plot(Age, Rate, ylab = expression(q[x]), xlab = "Age")
> lines(M2, n = 3, tr = quasibinomial()$linkinv, lwd = 1, col = "red")
> plot(Age, log(Deaths/Exposure), ylab = expression(log(mu[x])), xlab = "Age")
> lines(M1, n = 3, lwd = 1, col = "red")
> plot(Age, quasibinomial()$linkfun(Rate), ylab = expression(logit(q[x])), xlab = "Age")
> lines(M2, n = 3, lwd = 1, col = "red")
> par(op)
> 
> #########################################
> # bivariate example
> set.seed(123)
> doublesin <- function(x) {
+ # Adjusting the output to ensure it's positive
+ exp(sin(2*x[,1]) + sin(2*x[,2]))
+ }
> X <- round(runif(400, min = 0, max = 3), 2)
> Y <- round(runif(400, min = 0, max = 3), 2)
> # Calculate lambda for Poisson distribution
> lambda <- doublesin(cbind(X,Y))
> # Generate Z from Poisson distribution
> Z <- rpois(400, lambda)
> data <- data.frame(X, Y, Z)
> 
> # Fit a Poisson GeDS regression using GGeDS
> BivGeDS <- GGeDS(Z ~ f(X,Y), beta = 0.2, phi = 0.99, family = "poisson")
> 
> # Poisson mean deviance w.r.t data
> deviance(BivGeDS, n = 2) # or sum(poisson()$dev.resids(Z, BivGeDS$Linear.Fit$Predicted, wt = 1))
[1] 395.044
> deviance(BivGeDS, n = 3)
[1] 374.905
> deviance(BivGeDS, n = 4)
[1] 373.0677
> 
> # Poisson mean deviance w.r.t true function#' 
> f_XY <- apply(cbind(X, Y), 1, function(row) doublesin(matrix(row, ncol = 2)))
> mean(poisson()$dev.resids(f_XY, BivGeDS$Linear.Fit$Predicted, wt = 1))
[1] 0.1113454
> mean(poisson()$dev.resids(f_XY, BivGeDS$Quadratic.Fit$Predicted, wt = 1))
[1] 0.05282198
> mean(poisson()$dev.resids(f_XY, BivGeDS$Cubic.Fit$Predicted, wt = 1))
[1] 0.05211246
> 
> # Surface plot of the generating function (doublesin)
> plot(BivGeDS, f = doublesin)
Error in others$xlab %||% "X" : could not find function "%||%"
Calls: plot -> plot -> .local
Execution halted