--- title: "Analyst tools: predictions, residuals and validation" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Analyst tools: predictions, residuals and validation} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5, fig.align = "center" ) ``` This vignette shows the post-estimation tools for `brs()` fits: predictions on the mean and score scales, contrasts with standard errors, marginal effects, residuals, plots, the parametric bootstrap, cross-validation and summary tables. The model itself is described in `vignette("brs-intro")`; the mixed-model versions of these tools are in `vignette("brs-mm")`. Small numbers of replicates and folds keep the vignette fast, so the results are illustrative. ```{r library} library(betaregscale) ``` ## Data and fit Simulated NRS-11 scores of 240 patients at three post-operative times, with age and a treatment indicator: ```{r data} set.seed(2023) nrs <- data.frame( time = factor(rep(c("6h", "12h", "24h"), each = 80), levels = c("6h", "12h", "24h")), age = round(rnorm(240, 50, 12)), trt = rbinom(240, 1, 0.5) ) nrs$age_c <- (nrs$age - 50) / 10 eta <- -0.6 + c(0, 0.5, -0.3)[nrs$time] + 0.2 * nrs$age_c - 0.5 * nrs$trt shp <- brs_repar(mu = plogis(eta), phi = 0.2, repar = 2) nrs$y <- round(10 * rbeta(240, shp$shape1, shp$shape2)) fit <- brs(y ~ time + age_c + trt, data = nrs, ncuts = 10) fit ``` ## Predictions `predict()` returns, by `type`: | `type` | value | |:--|:--| | `"response"` (default) | the mean $E(Y) = \mu$ on $(0, 1)$ | | `"link"` | the linear predictor of the mean | | `"precision"` | the second parameter on its own scale (the dispersion $\sigma$ under `repar = 2`) | | `"variance"` | $\mathrm{Var}(Y)$ | | `"quantile"` | quantiles of $Y$ at the probabilities `at` | | `"score"` | the latent score at $E(Y)$ ($K\,E(Y)$ under `"mid"`) | | `"expected_score"` | the expected recorded score $\sum_s s\,P(S = s)$ | ```{r predict} nd <- data.frame(time = c("6h", "12h", "24h"), age_c = 0, trt = 0) cbind(nd, mu = predict(fit, nd), score = predict(fit, nd, type = "score"), expected_score = predict(fit, nd, type = "expected_score")) predict(fit, nd, type = "quantile", at = c(0.25, 0.5, 0.75)) ``` `brs_predict_scoreprob()` gives the probability of each score, as a matrix or in long format: ```{r scoreprob} round(brs_predict_scoreprob(fit, newdata = nd), 3) head(brs_predict_scoreprob(fit, newdata = nd, scores = 0:2, format = "long")) ``` ## Score-scale quantities with standard errors `brs_marginal()` returns score probabilities (`type = "prob"`), expected scores (`"expected_score"`) or exceedance probabilities $P(S \ge k)$ (`"exceed"`) with delta-method standard errors and Wald intervals, computed on the logit scale for probabilities and expected scores so they stay in range (contrasts use Wald intervals on the natural scale). `contrast` gives differences between pairs of rows of `newdata`; `c(i, j)` is row `i` minus row `j`: ```{r marginal} brs_marginal(fit, nd, type = "exceed", k = 4) brs_marginal(fit, nd, type = "expected_score", contrast = rbind(c(2, 1), c(3, 1))) ``` With `se_method = "simulation"`, the standard errors and percentile intervals come from `n_sim` draws of the parameters (set a seed first): ```{r marginal-sim} set.seed(1) brs_marginal(fit, nd[1:2, ], type = "expected_score", contrast = c(2, 1), se_method = "simulation", n_sim = 200) ``` ## Average marginal effects `brs_marginaleffects()` gives the average marginal effect of each numeric covariate on the mean (`type = "response"`) or on the linear predictor (`"link"`): a central difference for continuous covariates and the discrete change for 0/1 covariates, averaged over the rows. Intervals use `n_sim` draws from the asymptotic normal law of the estimates. Factors are not included; code a factor as 0/1 dummies to get its contrast, or use `brs_marginal()` contrasts as above. ```{r ame} set.seed(2) ame <- brs_marginaleffects(fit, n_sim = 100, keep_draws = TRUE) ame ``` Effects on $E(Y)$ are on $(0, 1)$; multiply by $K$ for the score scale under `"mid"`. ```{r ame-plot, fig.height = 3} autoplot.brs_marginaleffects(ame, type = "forest") ``` ## Residuals `residuals()` offers the types below. The randomized quantile residuals (`"rqr"`, Dunn and Smyth, 1996) are standard normal under the model whatever the censoring and are the recommended type for diagnostics. | `type` | definition | |:--|:--| | `"response"` (default) | `yt` minus the fitted mean | | `"pearson"` | response residual over the fitted SD | | `"deviance"` | beta deviance residual at `yt` | | `"rqr"` | $\Phi^{-1}(u_i)$, $u_i$ uniform on $(F(l_i), F(u_i))$, with 0 and 1 at the borders | | `"weighted"`, `"sweighted"` | residuals on the $\mathrm{logit}(Y)$ scale (Espinheira et al., 2008) | | `"cell_deviance"` | $\mathrm{sign}_i\sqrt{-2\log P_i}$ from the probability of the observed cell | All types except `"rqr"` and `"cell_deviance"` use one point of the cell, `yt`: $s/K$ under `"mid"` (so the border scores sit at $10^{-5}$ and $1 - 10^{-5}$) and the cell center under `"right"`/`"left"`. They are unreliable at the border scores. The cell deviance residuals are signed by the observed minus the expected score; for score data their squares add up to $-2\log L$: ```{r residuals} set.seed(1) r_q <- residuals(fit, type = "rqr") c(mean = mean(r_q), sd = sd(r_q)) r_d <- residuals(fit, type = "cell_deviance") c(sum_sq = sum(r_d^2), minus2logLik = -2 * as.numeric(logLik(fit))) ``` ## Plots `plot()` draws up to six diagnostic panels of quantile residuals (four by default) in base graphics or, with `gg = TRUE`, ggplot2: ```{r plot, fig.height = 6} plot(fit, gg = TRUE) ``` `autoplot.brs()` adds a calibration plot, the predicted score distribution, fitted CDFs and residuals by censoring type: ```{r autoplot, fig.height = 4} autoplot.brs(fit, type = "calibration") autoplot.brs(fit, type = "score_dist") ``` `brs_cens()` summarizes the censoring types of a fit: ```{r cens, fig.height = 5} brs_cens(fit, gg = TRUE) ``` ## Parametric bootstrap `brs_bootstrap()` draws new responses at the fitted shapes, re-coarsens them on the fit's cells and refits the original formula, `R` times. It returns percentile, basic, normal or BCa intervals next to the Wald intervals; failed refits are counted and printed. Use `R` of at least 199 in practice; 49 is used here for speed. ```{r bootstrap} set.seed(3) boot <- brs_bootstrap(fit, R = 49, keep_draws = TRUE) boot ``` ```{r bootstrap-plot, fig.height = 4} autoplot.brs_bootstrap(boot, type = "ci_forest") ``` ## Cross-validation `brs_cv()` runs repeated k-fold cross-validation. For each held-out fold it reports `log_score`, the mean log probability of the held-out observations (higher is better), and the RMSE and MAE of `yt` around the predicted mean: ```{r cv} set.seed(4) cv_full <- brs_cv(y ~ time + age_c + trt, data = nrs, k = 5, ncuts = 10) set.seed(4) cv_time <- brs_cv(y ~ time, data = nrs, k = 5, ncuts = 10) cv_full c(full = mean(cv_full$log_score), time_only = mean(cv_time$log_score)) ``` Resetting the seed gives both models the same folds. ## Tables `brs_est()` gives the coefficient table, `brs_gof()` the fit criteria and `brs_table()` a comparison of several fits with their censoring counts: ```{r tables} brs_est(fit) fit_time <- brs(y ~ time, data = nrs, ncuts = 10) brs_table(full = fit, time_only = fit_time, sort_by = "AIC") ``` ## References - Dunn, P. K., and Smyth, G. K. (1996). Randomized quantile residuals. *Journal of Computational and Graphical Statistics*, 5(3), 236–244. - Espinheira, P. L., Ferrari, S. L. P., and Cribari-Neto, F. (2008). On beta regression residuals. *Journal of Applied Statistics*, 35(4), 407–419. - Lopes, J. E. (2023). *Modelos de regressão beta para dados de escala*. Master's dissertation, Universidade Federal do Paraná, Curitiba.