diff --git a/R/compute_edf.R b/R/compute_edf.R index 4f3b005e..ea46d0c1 100644 --- a/R/compute_edf.R +++ b/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) } diff --git a/R/mvgam.R b/R/mvgam.R index 35412457..5b76c6e8 100644 --- a/R/mvgam.R +++ b/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', diff --git a/R/residuals.mvgam.R b/R/residuals.mvgam.R index b6c86012..933251ff 100644 --- a/R/residuals.mvgam.R +++ b/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) diff --git a/R/summary.mvgam.R b/R/summary.mvgam.R index fffff26e..54eb2b4b 100644 --- a/R/summary.mvgam.R +++ b/R/summary.mvgam.R @@ -8,7 +8,7 @@ #'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 @@ -16,10 +16,12 @@ #'@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), diff --git a/README.Rmd b/README.Rmd index 0076178f..9f386cf7 100644 --- a/README.Rmd +++ b/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) ``` diff --git a/README.md b/README.md index 045d4a69..a9273108 100644 --- a/README.md +++ b/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) Plotting forecast distributions using mvgam in R #> 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) diff --git a/docs/reference/Rplot002.png b/docs/reference/Rplot002.png index 6d219198..c4d06b58 100644 Binary files a/docs/reference/Rplot002.png and b/docs/reference/Rplot002.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-13-1.png b/docs/reference/figures/README-unnamed-chunk-13-1.png index 11dd9bdb..09f75cdb 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-13-1.png and b/docs/reference/figures/README-unnamed-chunk-13-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-14-1.png b/docs/reference/figures/README-unnamed-chunk-14-1.png index fb9b521f..36bc2232 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-14-1.png and b/docs/reference/figures/README-unnamed-chunk-14-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-15-1.png b/docs/reference/figures/README-unnamed-chunk-15-1.png index b9314215..0aef7e7b 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-15-1.png and b/docs/reference/figures/README-unnamed-chunk-15-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-16-1.png b/docs/reference/figures/README-unnamed-chunk-16-1.png index 2bb0b7c2..78dc3d14 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-16-1.png and b/docs/reference/figures/README-unnamed-chunk-16-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-17-1.png b/docs/reference/figures/README-unnamed-chunk-17-1.png index 7581f43b..6afa071d 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-17-1.png and b/docs/reference/figures/README-unnamed-chunk-17-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-18-1.png b/docs/reference/figures/README-unnamed-chunk-18-1.png index 747a5f34..c799fe95 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-18-1.png and b/docs/reference/figures/README-unnamed-chunk-18-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-19-1.png b/docs/reference/figures/README-unnamed-chunk-19-1.png index 2f6d21eb..a66e6d6d 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-19-1.png and b/docs/reference/figures/README-unnamed-chunk-19-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-20-1.png b/docs/reference/figures/README-unnamed-chunk-20-1.png index 45c8febf..5b25c223 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-20-1.png and b/docs/reference/figures/README-unnamed-chunk-20-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-21-1.png b/docs/reference/figures/README-unnamed-chunk-21-1.png index f5d5ea3b..40dd2644 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-21-1.png and b/docs/reference/figures/README-unnamed-chunk-21-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-22-1.png b/docs/reference/figures/README-unnamed-chunk-22-1.png index e96f84e0..59250c51 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-22-1.png and b/docs/reference/figures/README-unnamed-chunk-22-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-23-1.png b/docs/reference/figures/README-unnamed-chunk-23-1.png index b5260a34..24340f18 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-23-1.png and b/docs/reference/figures/README-unnamed-chunk-23-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-24-1.png b/docs/reference/figures/README-unnamed-chunk-24-1.png index e59e02c0..9ba85217 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-24-1.png and b/docs/reference/figures/README-unnamed-chunk-24-1.png differ diff --git a/docs/reference/index.html b/docs/reference/index.html index 0ad03cd1..866c1d3e 100644 --- a/docs/reference/index.html +++ b/docs/reference/index.html @@ -10,7 +10,7 @@ mvgam - 1.0.91 + 1.1.0