## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  fig.align = "center"
)

## ----library------------------------------------------------------------------
library(betaregscale)

## ----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

## ----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))

## ----scoreprob----------------------------------------------------------------
round(brs_predict_scoreprob(fit, newdata = nd), 3)
head(brs_predict_scoreprob(fit, newdata = nd, scores = 0:2, format = "long"))

## ----marginal-----------------------------------------------------------------
brs_marginal(fit, nd, type = "exceed", k = 4)
brs_marginal(fit, nd, type = "expected_score",
             contrast = rbind(c(2, 1), c(3, 1)))

## ----marginal-sim-------------------------------------------------------------
set.seed(1)
brs_marginal(fit, nd[1:2, ], type = "expected_score", contrast = c(2, 1),
             se_method = "simulation", n_sim = 200)

## ----ame----------------------------------------------------------------------
set.seed(2)
ame <- brs_marginaleffects(fit, n_sim = 100, keep_draws = TRUE)
ame

## ----ame-plot, fig.height = 3-------------------------------------------------
autoplot.brs_marginaleffects(ame, type = "forest")

## ----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)))

## ----plot, fig.height = 6-----------------------------------------------------
plot(fit, gg = TRUE)

## ----autoplot, fig.height = 4-------------------------------------------------
autoplot.brs(fit, type = "calibration")
autoplot.brs(fit, type = "score_dist")

## ----cens, fig.height = 5-----------------------------------------------------
brs_cens(fit, gg = TRUE)

## ----bootstrap----------------------------------------------------------------
set.seed(3)
boot <- brs_bootstrap(fit, R = 49, keep_draws = TRUE)
boot

## ----bootstrap-plot, fig.height = 4-------------------------------------------
autoplot.brs_bootstrap(boot, type = "ci_forest")

## ----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))

## ----tables-------------------------------------------------------------------
brs_est(fit)
fit_time <- brs(y ~ time, data = nrs, ncuts = 10)
brs_table(full = fit, time_only = fit_time, sort_by = "AIC")