update edf calculations

R/compute_edf.R

@@ -1,7 +1,8 @@
 #' Compute approximate EDFs of smooths
 #' @importFrom stats fitted
 #'@noRd
-compute_edf = function(mgcv_model, object, rho_names, sigma_raw_names){
+compute_edf = function(mgcv_model, object, rho_names, sigma_raw_names,
+                       conservative = FALSE){
 
   if(length(mgcv_model$smooth) > 0){
 
@@ -88,8 +89,11 @@ compute_edf = function(mgcv_model, object, rho_names, sigma_raw_names){
         mu[which(mu_variance == 0)]
     }
 
-    w <- as.numeric(mgcv_model$family$mu.eta(eta)^2 / mu_variance)
-    XWX <- t(X) %*% (w * X)
+    if(!conservative){
+      w <- as.numeric(mgcv_model$family$mu.eta(as.vector(eta))^2 / mu_variance)
+      XWX <- t(X) %*% (w * X)
+    } else XWX <- t(X) %*% X
+
     lambda <- mgcv_model$sp
     XWXS <- XWX
 
@@ -104,7 +108,18 @@ compute_edf = function(mgcv_model, object, rho_names, sigma_raw_names){
       names(edf) <- names(coef(mgcv_model))
     }
     mgcv_model$edf <- edf
+    mgcv_model$edf1 <- edf
+    mgcv_model$edf2 <- edf
   }
 
+  # Add frequentist version of parameter covariance estimators
+  # for calculation of smooth term p-values;
+  # rV <- mroot(mgcv_model$Vp)
+  # V <- tcrossprod(rV)
+  # XWX <- crossprod(mgcv_model$R)
+  # Ve_hat <- V %*% XWX
+  # mgcv_model$Ve <- Ve_hat %*% V
+
+  mgcv_model$Ve <- mgcv_model$Vp
   return(mgcv_model)
 }

R/mvgam.R

@@ -2168,7 +2168,7 @@ mvgam = function(formula,
   # also necessary for computing EDFs and approximate p-values of smooths
   if(!lfo){
     V <- cov(mcmc_chains(out_gam_mod, 'b'))
-    ss_gam$Ve <- ss_gam$Vp <- ss_gam$Vc <- V
+    ss_gam$Vp <- ss_gam$Vc <- V
 
     # Add the posterior median coefficients
     p <- mcmc_summary(out_gam_mod, 'b',
@@ -2182,7 +2182,7 @@ mvgam = function(formula,
     # Repeat for any trend-specific mgcv model
     if(!missing(trend_formula)){
       V <- cov(mcmc_chains(out_gam_mod, 'b_trend'))
-      trend_mgcv_model$Ve <- trend_mgcv_model$Vp <- trend_mgcv_model$Vc <- V
+      trend_mgcv_model$Vp <- trend_mgcv_model$Vc <- V
       p <- mcmc_summary(out_gam_mod, 'b_trend',
                         variational = algorithm %in% c('meanfield',
                                                        'fullrank',

R/residuals.mvgam.R

@@ -55,15 +55,15 @@ residuals.mvgam <- function(object,
   }
 
   if(summary){
-    Qupper <- apply(resid_matrix, 2, quantile, probs = max(probs), na.rm = TRUE)
-    Qlower <- apply(resid_matrix, 2, quantile, probs = min(probs), na.rm = TRUE)
+    Qupper <- apply(resid_matrix, 1, quantile, probs = max(probs), na.rm = TRUE)
+    Qlower <- apply(resid_matrix, 1, quantile, probs = min(probs), na.rm = TRUE)
 
     if(robust){
-      estimates <- apply(resid_matrix, 2, median, na.rm = TRUE)
-      errors <- apply(abs(resid_matrix - estimates), 2, median, na.rm = TRUE)
+      estimates <- apply(resid_matrix, 1, median, na.rm = TRUE)
+      errors <- apply(abs(resid_matrix - estimates), 1, median, na.rm = TRUE)
     } else {
-      estimates <- apply(resid_matrix, 2, mean, na.rm = TRUE)
-      errors <- apply(resid_matrix, 2, sd, na.rm = TRUE)
+      estimates <- apply(resid_matrix, 1, mean, na.rm = TRUE)
+      errors <- apply(resid_matrix, 1, sd, na.rm = TRUE)
     }
 
     out <- cbind(estimates, errors, Qlower, Qupper)

R/summary.mvgam.R

@@ -8,18 +8,20 @@
 #'Defaults to \code{TRUE} but use \code{FALSE} for a more concise summary
 #'@param smooth_test Logical. Compute estimated degrees of freedom and approximate
 #'p-values for smooth terms? Defaults to \code{TRUE}, but users may wish to set
-#'to \code{FALSE} for complex models with many smooth terms
+#'to \code{FALSE} for complex models with many smooth or random effect terms
 #'@param digits The number of significant digits for printing out the summary;
 #'  defaults to \code{2}.
 #'@param ... Ignored
 #'@author Nicholas J Clark
 #'@details `summary.mvgam` and `summary.mvgam_prefit` return brief summaries of the model's call, along with posterior intervals for
 #'some of the key parameters in the model. Note that some smooths have extra penalties on the null space,
 #'so summaries for the \code{rho} parameters may include more penalty terms than the number of smooths in
-#'the original model formula. Approximate p-values for smooth terms are also returned, with methods used for their
+#'the original model formula. Approximate p-values for smooth terms are also returned,
+#'with methods used for their
 #'calculation following those used for `mgcv` equivalents (see \code{\link[mgcv]{summary.gam}} for details).
 #'The Estimated Degrees of Freedom (edf) for smooth terms is computed using
-#'`edf.type = 1` as described in the documentation for \code{\link[mgcv]{jagam}}. Experiments suggest
+#'either `edf.type = 1` for models with no trend component, or `edf.type = 0` for models with
+#'trend components. These are described in the documentation for \code{\link[mgcv]{jagam}}. Experiments suggest
 #'these p-values tend to be more conservative than those that might be returned from an equivalent
 #'model fit with \code{\link[mgcv]{summary.gam}} using `method = 'REML'`
 #'
@@ -38,10 +40,16 @@ summary.mvgam = function(object,
 
 #### Smooth tests ####
   if(smooth_test){
+    if(inherits(object$trend_model, 'mvgam_trend')){
+      trend_model <- object$trend_model$label
+    } else {
+      trend_model <- object$trend_model
+    }
     object$mgcv_model <- compute_edf(object$mgcv_model,
                                      object,
                                      'rho',
-                                     'sigma_raw')
+                                     'sigma_raw',
+                                     conservative = trend_model == 'None')
 
     if(!is.null(object$trend_call)){
       object$trend_mgcv_model <- compute_edf(object$trend_mgcv_model,
@@ -287,7 +295,7 @@ if(any(!is.na(object$sp_names)) & smooth_test){
       if(!is.null(object$trend_call)){
         cat("\nApproximate significance of GAM observation smooths:\n")
       } else {
-        cat("\nApproximate significance of GAM observation smooths:\n")
+        cat("\nApproximate significance of GAM smooths:\n")
       }
       suppressWarnings(printCoefmat(gam_sig_table,
                                     digits = min(3, digits + 1),
@@ -299,7 +307,7 @@ if(any(!is.na(object$sp_names)) & smooth_test){
     if(!is.null(object$trend_call)){
       cat("\nApproximate significance of GAM observation smooths:\n")
     } else {
-      cat("\nApproximate significance of GAM observation smooths:\n")
+      cat("\nApproximate significance of GAM smooths:\n")
     }
     suppressWarnings(printCoefmat(gam_sig_table,
                                   digits = min(3, digits + 1),

README.Rmd

@@ -113,7 +113,7 @@ lynx_mvgam <- mvgam(data = lynx_train,
                newdata = lynx_test,
                formula = population ~ s(season, bs = 'cc', k = 12),
                knots = list(season = c(0.5, 19.5)),
-               family = 'poisson',
+               family = poisson(),
                trend_model = AR(p = 3),
                use_stan = TRUE)
 ```

README.md

@@ -327,29 +327,29 @@ summary(lynx_mvgam)
 #> 
 #> 
 #> GAM coefficient (beta) estimates:
-#>                2.5%    50%  97.5% Rhat n_eff
-#> (Intercept)   6.100  6.600  7.000 1.00   658
-#> s(season).1  -0.630  0.058  0.730 1.00   829
-#> s(season).2  -0.150  0.750  1.800 1.00   489
-#> s(season).3  -0.056  1.100  2.400 1.01   448
-#> s(season).4  -0.550  0.410  1.400 1.01   695
-#> s(season).5  -1.200 -0.059  0.960 1.00   474
-#> s(season).6  -1.100  0.058  1.100 1.00   560
-#> s(season).7  -0.740  0.350  1.400 1.00   865
-#> s(season).8  -1.000  0.160  1.700 1.00   473
-#> s(season).9  -1.100 -0.330  0.600 1.00   572
-#> s(season).10 -1.300 -0.640 -0.031 1.00   827
+#>                 2.5%      50%  97.5% Rhat n_eff
+#> (Intercept)   6.1000  6.60000  7.000 1.01   758
+#> s(season).1  -0.6400  0.02600  0.690 1.00   810
+#> s(season).2  -0.2400  0.81000  1.800 1.02   227
+#> s(season).3  -0.0084  1.20000  2.400 1.02   219
+#> s(season).4  -0.5000  0.44000  1.300 1.00   680
+#> s(season).5  -1.2000 -0.10000  0.950 1.02   444
+#> s(season).6  -1.1000  0.00088  1.100 1.01   564
+#> s(season).7  -0.7300  0.36000  1.500 1.00   673
+#> s(season).8  -0.9800  0.24000  1.800 1.02   337
+#> s(season).9  -1.2000 -0.29000  0.680 1.02   450
+#> s(season).10 -1.4000 -0.66000 -0.025 1.01   451
 #> 
-#> Approximate significance of GAM observation smooths:
+#> Approximate significance of GAM smooths:
 #>            edf Ref.df Chi.sq p-value
-#> s(season) 9.93     10  16564    0.29
+#> s(season) 9.97     10  19379    0.24
 #> 
 #> Latent trend AR parameter estimates:
 #>           2.5%   50% 97.5% Rhat n_eff
-#> ar1[1]    0.79  1.10 1.500    1   751
-#> ar2[1]   -0.84 -0.41 0.048    1  1406
-#> ar3[1]   -0.49 -0.14 0.250    1   680
-#> sigma[1]  0.40  0.51 0.650    1  1074
+#> ar1[1]    0.75  1.10  1.40 1.01   703
+#> ar2[1]   -0.85 -0.40  0.04 1.00  1839
+#> ar3[1]   -0.47 -0.13  0.31 1.01   393
+#> sigma[1]  0.41  0.51  0.64 1.00  1027
 #> 
 #> Stan MCMC diagnostics:
 #> n_eff / iter looks reasonable for all parameters
@@ -358,7 +358,7 @@ summary(lynx_mvgam)
 #> 0 of 2000 iterations saturated the maximum tree depth of 12 (0%)
 #> E-FMI indicated no pathological behavior
 #> 
-#> Samples were drawn using NUTS(diag_e) at Wed Mar 27 6:38:25 PM 2024.
+#> Samples were drawn using NUTS(diag_e) at Thu Mar 28 8:44:38 PM 2024.
 #> For each parameter, n_eff is a crude measure of effective sample size,
 #> and Rhat is the potential scale reduction factor on split MCMC chains
 #> (at convergence, Rhat = 1)
@@ -480,7 +480,7 @@ plot(lynx_mvgam, type = 'forecast', newdata = lynx_test)
 <img src="man/figures/README-unnamed-chunk-21-1.png" alt="Plotting forecast distributions using mvgam in R" width="60%" style="display: block; margin: auto;" />
 
     #> Out of sample CRPS:
-    #> [1] 2941.682
+    #> [1] 2942.292
 
 And the estimated latent trend component, again using the more flexible
 `plot_mvgam_...()` option to show first derivatives of the estimated
@@ -615,7 +615,7 @@ summary(mod, include_betas = FALSE)
 #>              2.5%  50% 97.5% Rhat n_eff
 #> (Intercept) -0.19 0.19  0.44    1   822
 #> 
-#> Approximate significance of GAM observation smooths:
+#> Approximate significance of GAM smooths:
 #>                            edf Ref.df Chi.sq p-value    
 #> s(season)                4.639      5  35.50  <2e-16 ***
 #> s(season):seriesseries_1 0.887      4   0.53    0.98    
@@ -641,7 +641,7 @@ summary(mod, include_betas = FALSE)
 #> 0 of 2000 iterations saturated the maximum tree depth of 12 (0%)
 #> E-FMI indicated no pathological behavior
 #> 
-#> Samples were drawn using NUTS(diag_e) at Wed Mar 27 6:40:30 PM 2024.
+#> Samples were drawn using NUTS(diag_e) at Thu Mar 28 8:46:19 PM 2024.
 #> For each parameter, n_eff is a crude measure of effective sample size,
 #> and Rhat is the potential scale reduction factor on split MCMC chains
 #> (at convergence, Rhat = 1)