Skip to contents

kumarsurv

Companion R package to:

Ospina, R., Pimentel, J. S. and Cribari-Neto, F. Copula-Based Kumaraswamy Regression for Anomaly Detection in Doubly Bounded Time Series (manuscript submitted for publication).

A marginal regression model for time series in the unit interval, with a Kumaraswamy marginal distribution parametrized directly by a quantile of prescribed level and serial dependence induced by a Gaussian copula with ARMA structure. Anomaly detection combines bootstrap-calibrated goodness-of-fit tests, one-step prediction bands and a cumulative sum (CUSUM) chart calibrated through the in-control average run length. Applied to weekly acute diarrheal disease incidence in Platina, SP (Brazil).

Site: https://raydonal.github.io/kumarsurv/

Installation 

# install.packages("remotes")
remotes::install_github("Raydonal/kumarsurv")

Usage 

library(kumarsurv)
data(platina)

X <- cbind(intercept = 1, trend = platina$trend, sin = platina$sin, cos = platina$cos)
Z <- cbind(intercept = 1, trend = platina$trend, sin = platina$sin, cos = platina$cos)
y <- platina$incidence
y[y == 0] <- 1e-7  # boundary replacement, see the paper's Section 6.1

fit <- fit_copula(y, X, Z, p = 1, q = 1, family = "kuma", tau = 0.5)
r <- qresiduals(fit)$r  # predictive quantile residuals

For tau other than the median, use fit_copula_multi_tau() instead of a bare fit_copula() call: the likelihood surface can have more than one interior stationary point at a given tau, and fit_copula_multi_tau() runs the monotonicity/invariance diagnostic the paper’s Section 6.3 and Table 3 rely on to pick the right one. See vignette("reproducing-the-paper").

Reproducing the paper

Every table, figure and numerical value in the paper and its supplementary material is generated by a script, all 53 of them browsable directly at inst/scripts/ on GitHub (also available at system.file("scripts", package = "kumarsurv") once the package is installed). inst/scripts/run_all.R runs them in dependency order:

scripts_dir <- system.file("scripts", package = "kumarsurv")
# fast scripts + figures:
source(file.path(scripts_dir, "run_all.R"))
# + the M = 300 / B = 300 Monte Carlo / bootstrap runs (slow, ~2h45):
Sys.setenv(RUN_PESADOS = "1")
source(file.path(scripts_dir, "run_all.R"))

inst/scripts/copula_ml.R and R/copula_ml.R are deliberately two separate files, not one kept in sync with the other. inst/scripts/copula_ml.R is a frozen, byte-for-byte archival copy of the original script that generated every number in the paper: it is never edited, so the scripts in inst/scripts/ (which source("copula_ml.R") it directly, not the package’s namespace) always reproduce exactly what was published, regardless of any later change to the package. R/copula_ml.R is the same model-fitting code documented and exported for general use; it has since received two small robustness fixes (parameter-name fallback when X/Z lack column names; suppressing, not changing, spurious warnings from out-of-domain BFGS excursions) that do not exist in the frozen script because they do not need to: both were verified, before being committed, to leave every fitted value (par, se, loglik, aic) bit-for-bit unchanged. Concretely, running inst/scripts/01_dados_reais_SEM.R from a clean checkout of this repository, untouched by anything in R/, reproduces Table 4 of the paper exactly: AIC -5064.54 (Kumaraswamy) and -5086.40 (beta) under ARMA(1,1), to the last printed digit. config_paths.R resolves the scripts’ three legacy path conventions regardless of the working directory.

Result → script

Full mapping, every table and figure in both the paper and its supplementary material (this is the complete table that the paper’s own supplementary material points here for, not a subset):

Paper

Result Script
Table 1 (tail calibration under misspecification) 18_calib_cauda.R
Table 2 (estimates, Phase I) 01_dados_reais_SEM.R (estimates); analise_covariaveis.R (standard errors, tests, Wald interval, alternative trends)
Table 3 (multi-tau structure) analise_multiquantil.R; multitau_controle.R (in-control coverage); calib_SEM.R and teste_multi_tau.R (note)
Table 4 (Kumaraswamy vs. beta fit) 01_dados_reais_SEM.R; 13_bootstrap_gof.R (W2W^2, A2A^2 and pp-values)
Table 5 (calibration of control limits) calib_controle_SEM.R
Figures (series, control limits, CUSUM, detection) fig3_horizontal.R; 06_figuras_SEM.R; fig4_cusum_SEM.R; 07_fig_deteccao_SEM.R
Simulation: comparison between marginals/parametrizations 17_crps_fix_B300.R (CRPS); 15_sim_comp_B300.R (pinball loss); 16_simcomp_recovery_B300.R (biases)
Simulation: serial dependence and ARL calibration 09_efeito_correlacao.R; 10_banda_simulada.R; 11_arl.R
Application: ARMA selection, Ljung-Box, AR profile tab_arma_selecao.R; analise_inflacao_diagnostico.R; perfil_ar1_exato.R
Application: out-of-sample, sensitivity to ε\varepsilon 08_vantagem_SEM.R; sensib_zeros_controle.R
Application: Farrington comparison and the chart 02_farrington_SEM.R
Hurdle model (preliminary, Concluding remarks) analise_harmonicos_inflacao.R

Supplementary material

Result Script
Parameter recovery, design of Varin (2014) 14_mc_validacao_B300.R
Parameter recovery, application design, and tied values simulacao_desenho_aplicado_300.R; simaplic_empates.R
Parameter recovery under each marginal distribution 16_simcomp_recovery_B300.R
Effect of ignoring correlation (table) + prediction-band figure 09_efeito_correlacao.R; 10_banda_simulada.R
ARL figure (chart-limit calibration curves) 11_arl.R
ARMA order selection (full AIC table) tab_arma_selecao.R
Harmonic periodicity + AR likelihood profile tables perfil_ar1_exato.R; analise_harmonicos_inflacao.R (nested submodel)
Sensitivity to the replacement of zeros, bootstrap stability sensib_zeros_controle.R; 13_bootstrap_gof.R (stability)
Out-of-sample predictive performance (full table) 08_vantagem_SEM.R
Calibration of control limits, with/without trend calib_controle_SEM.R
Density and quantile-curve figures fig_densidades.R; fig_quantis_corrigida.R

Verification scripts (inst/scripts/verificacao_R47/, verificacao_R50/, verificacao_R51/)

What it checks Script
Exact likelihood (Kalman vs. Durbin-Levinson) R47_valida_lib.R; R47_vero_exata.R
Numerical Hessian step size R47_teste_hessiana.R
Definition of the outbreak weeks R47_calib_robustez.R
The two stationary points at tau = 0.90 R50_tau090_pontos.R
Alternative trends, restricted maxima R50_tendencias_semestral.R
Constant-shape model: likelihood-ratio test of the shape coefficients, sensitivity of the control limits and of the CUSUM (Section S6.7) R51_forma_constante_sensib.R

Two known, environment-dependent non-reproducibilities, documented in the paper’s supplementary material: the Kumaraswamy bootstrap pp-value (13_bootstrap_gof.R; Section S6.4, “Sensitivity to the replacement of zeros and stability of the bootstrap”; five runs and their Monte Carlo standard error are reported), and fits where the likelihood is nearly flat (tau = 0.90, near the ARMA stationarity boundary; Section S6.6, “Calibration of the control limits”, paragraph “The level tau = 0.90”), where the last digit of floating-point precision can select a different stationary point.

Citation

The paper is submitted to Biostatistics (Oxford University Press) and not yet published; cite it as a manuscript in preparation until it is accepted.

@unpublished{ospina2026kumarsurv,
  title  = {Copula-Based Kumaraswamy Regression for Anomaly Detection in Doubly Bounded Time Series},
  author = {Ospina, Raydonal and Pimentel, Jonatha S. and Cribari-Neto, Francisco},
  year   = {2026},
  note   = {Manuscript submitted for publication}
}

License

MIT © Raydonal Ospina, Jonatha S. Pimentel and Francisco Cribari-Neto