| Title: | Clinical Inference Engine for Physiological Signal Studies |
|---|---|
| Description: | Statistical inference for rehabilitation studies: mixed-effects / MMRM longitudinal models, single-case (N-of-1) designs, estimands with multiple imputation, and per-estimate uncertainty. Heavy modelling backends are optional (guarded) Suggests. |
| Authors: | Yusuke Matsui [aut, cre] |
| Maintainer: | Yusuke Matsui <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.3.1 |
| Built: | 2026-10-03 13:43:52 UTC |
| Source: | https://github.com/x-biosignal/PhysioClinStats |
End-to-end analysis of a declared defineEstimand estimand:
multiply-impute the dropout-missing longitudinal response (in wide form, so
the within-subject correlation is respected), fit an MMRM
(fitMMRM) to each completed data set, and pool the fixed
effects with Rubin's rules (poolEstimates()). Treatment-policy and
hypothetical strategies use the standard MAR imputation; other strategies
warn that they need a bespoke imputation model.
analyseEstimand( data, estimand, response, treatment, time, subject, covariates = NULL, m = 20, method = "norm", seed = NULL, covariance = "unstructured" )analyseEstimand( data, estimand, response, treatment, time, subject, covariates = NULL, m = 20, method = "norm", seed = NULL, covariance = "unstructured" )
data |
Long-format data (subject x time rows) with the response, treatment, time, subject, and any (subject-level) covariates. |
estimand |
An |
response, treatment, time, subject
|
Column names. |
covariates |
Optional subject-level covariate column names. |
m |
Number of imputations (default 20). |
method |
mice method (default |
seed |
Optional integer seed. |
covariance |
MMRM covariance structure (default |
An AnalysisResult (type "estimand_analysis") carrying
the estimand attributes and the pooled fixed-effect table.
ICH E9(R1); Rubin 1987; Mallinckrodt 2008 (MMRM).
defineEstimand(), multipleImputation(), poolEstimates()
For a numeric sample, returns the analytic Normal-Normal conjugate posterior
for the mean: a normal prior N(prior_mean, prior_sd^2) combined with a
normal likelihood of known scale sigma gives a normal posterior in closed
form (a flat prior, prior_sd = Inf, reproduces the sampling posterior
N(mean(y), sigma^2/n)). For a model formula the estimation is delegated
to rstanarm/brms when installed, and errors informatively otherwise – the
analytic numeric-vector path always runs.
bayesEstimate( y, data = NULL, prior_mean = 0, prior_sd = Inf, sigma = NULL, level = 0.95, seed = NULL )bayesEstimate( y, data = NULL, prior_mean = 0, prior_sd = Inf, sigma = NULL, level = 0.95, seed = NULL )
y |
A numeric vector (analytic path) or a model |
data |
Data frame for the |
prior_mean, prior_sd
|
Normal prior mean and SD (default 0 and |
sigma |
Known likelihood SD; |
level |
Credible level (default 0.95). |
seed |
Optional integer random seed recorded in the provenance. |
An AnalysisResult (from PhysioCore) with type = "bayes" whose
result holds posterior_mean, posterior_sd, ci_lower, ci_upper,
the level/method, plus a provenance log with the seed.
Gelman, A. et al. (2013). Bayesian Data Analysis, 3rd ed.
credibleInterval(), conformalInterval()
set.seed(1) bayesEstimate(rnorm(30, mean = 5), prior_mean = 0, prior_sd = 10)set.seed(1) bayesEstimate(rnorm(30, mean = 5), prior_mean = 0, prior_sd = 10)
Calls mediation::mediate() after checking that the mediator and outcome
models describe the same ordered observations. Estimates rely on sequential
ignorability, consistency, positivity, absence of a treatment-induced
mediator-outcome confounder, and correct mediator and outcome models. An
indirect association, even when statistically significant, is not evidence
that a biological mechanism has been established.
causalMediation( model_m, model_y, treat, mediator, control.value = 0, treat.value = 1, covariates = NULL, sims = 1000L, boot = FALSE, boot.ci.type = "perc", conf.level = 0.95, robustSE = FALSE, sensitivity = TRUE, seed = NULL, ... )causalMediation( model_m, model_y, treat, mediator, control.value = 0, treat.value = 1, covariates = NULL, sims = 1000L, boot = FALSE, boot.ci.type = "perc", conf.level = 0.95, robustSE = FALSE, sensitivity = TRUE, seed = NULL, ... )
model_m, model_y
|
Supported fitted mediator and outcome models, or
|
treat, mediator
|
Treatment and mediator variable names. |
control.value, treat.value
|
Treatment contrast supplied to the backend. |
covariates |
Optional named list fixing covariates for a conditional mediation estimand. |
sims |
Positive integer number of simulations or bootstrap replicates. |
boot |
Use the nonparametric bootstrap. |
boot.ci.type |
Bootstrap confidence-interval type. |
conf.level |
Confidence level in |
robustSE |
Request heteroskedasticity-consistent uncertainty for
supported |
sensitivity |
Run |
seed |
Optional integer seed. The caller's global RNG state is restored. |
... |
Further arguments passed to |
An AnalysisResult with a complete ACME/ADE/total/proportion table,
the backend object, settings, formulas, assumptions, and sensitivity
diagnostics.
The inference functions load heavy modelling packages only when present
(guarded via requireNamespace). This reports availability so callers
can fail early with a helpful message.
clinStatsBackends()clinStatsBackends()
A named logical vector for mmrm, emmeans, lme4,
and SingleCaseES.
clinStatsBackends()clinStatsBackends()
Distribution-free predictive intervals with finite-sample marginal coverage
at least 1 - alpha, requiring only exchangeability. "split" conformal
(Lei et al. 2018) scores an independent calibration set with the fitted
model and takes the conformal quantile of the absolute residuals;
"jackknife_plus" (Barber et al. 2021) refits the model leave-one-out over
the calibration data and combines the leave-one-out predictions and
residuals.
conformalInterval( model, calib, newx, alpha = 0.1, type = c("split", "jackknife_plus"), seed = NULL )conformalInterval( model, calib, newx, alpha = 0.1, type = c("split", "jackknife_plus"), seed = NULL )
model |
A fitted model supporting |
calib |
A data frame with the response and predictors: an independent
calibration set for |
newx |
A data frame of predictor values to predict for. |
alpha |
Miscoverage level (default 0.1 for 90% intervals). |
type |
|
seed |
Optional integer random seed recorded in the provenance. |
The coverage guarantee requires the calibration data to be exchangeable with,
and (for "split") independent of, the data the model was trained on. When
the model uses a transformed response (e.g. log(y) ~ x) the interval is on
that transformed scale.
An AnalysisResult (from PhysioCore) with type = "conformal" whose
result holds point, lower, upper (one per newx row), the
alpha/method/quantile, and a provenance log capturing the seed.
Lei, J. et al. (2018). JASA 113(523):1094-1111. Barber, R.F. et al. (2021). Ann Statist 49(1):486-507.
conformalPredict(), bayesEstimate()
set.seed(1) train <- data.frame(x = rnorm(60)); train$y <- 2 * train$x + rnorm(60) calib <- data.frame(x = rnorm(60)); calib$y <- 2 * calib$x + rnorm(60) fit <- lm(y ~ x, data = train) res <- conformalInterval(fit, calib, data.frame(x = c(-1, 0, 1))) PhysioExperiment::resultValue(res)$lowerset.seed(1) train <- data.frame(x = rnorm(60)); train$y <- 2 * train$x + rnorm(60) calib <- data.frame(x = rnorm(60)); calib$y <- 2 * calib$x + rnorm(60) fit <- lm(y ~ x, data = train) res <- conformalInterval(fit, calib, data.frame(x = c(-1, 0, 1))) PhysioExperiment::resultValue(res)$lower
Convenience wrapper around split-conformal conformalInterval() for one new
observation, returning the guaranteed-coverage predicted band.
conformalPredict(model, calib, newx, alpha = 0.1, seed = NULL)conformalPredict(model, calib, newx, alpha = 0.1, seed = NULL)
model |
A fitted model (see |
calib |
The calibration data frame (response + predictors). |
newx |
A one-row data frame of predictors for the new patient. |
alpha |
Miscoverage level (default 0.1). |
seed |
Optional integer random seed. |
An AnalysisResult with type = "conformal" for the single patient.
set.seed(1) d <- data.frame(x = rnorm(80)); d$y <- 1.5 * d$x + rnorm(80) fit <- lm(y ~ x, data = d[1:40, ]) conformalPredict(fit, d[41:80, ], data.frame(x = 0.5))set.seed(1) d <- data.frame(x = rnorm(80)); d$y <- 1.5 * d$x + rnorm(80) fit <- lm(y ~ x, data = d[1:40, ]) conformalPredict(fit, d[41:80, ], data.frame(x = 0.5))
Fits a Cox proportional-hazards model with survival::coxph, tests the
proportional-hazards assumption with survival::cox.zph, and warns when the
global test is significant. Returns hazard ratios with confidence intervals.
coxModel( data, formula = NULL, time = "time", event = "event", covariates = NULL, conf_level = 0.95, ph_alpha = 0.05 )coxModel( data, formula = NULL, time = "time", event = "event", covariates = NULL, conf_level = 0.95, ph_alpha = 0.05 )
data |
A data.frame. |
formula |
A model formula whose response is |
time, event, covariates
|
Alternative to |
conf_level |
Confidence level (default 0.95). |
ph_alpha |
Significance level for the global PH warning (default 0.05). |
An AnalysisResult (type = "cox_model") whose estimate is the
hazard-ratio vector, with result$coefficients (coef, HR, SE, z, p, CI),
result$ph_test (the cox.zph table), result$ph_violated and
result$fit.
Cox DR (1972). Regression models and life-tables. JRSS B, 34(2). Grambsch PM, Therneau TM (1994). Proportional hazards tests and diagnostics based on weighted residuals. Biometrika, 81(3).
survivalFit(), milestoneHazard()
if (requireNamespace("survival", quietly = TRUE)) { coxModel(survival::lung, time = "time", event = "status", covariates = c("age", "sex", "ph.ecog")) }if (requireNamespace("survival", quietly = TRUE)) { coxModel(survival::lung, time = "time", event = "status", covariates = c("age", "sex", "ph.ecog")) }
The highest-density interval (HDI, the narrowest interval carrying level
probability) or the equal-tailed quantile interval, from a vector of
posterior draws.
credibleInterval(x, level = 0.95, method = c("hdi", "quantile"))credibleInterval(x, level = 0.95, method = c("hdi", "quantile"))
x |
Numeric vector of posterior samples. |
level |
Credible level (default 0.95). |
method |
|
Numeric c(lower, upper).
set.seed(1) credibleInterval(rnorm(2000), method = "hdi")set.seed(1) credibleInterval(rnorm(2000), method = "hdi")
Constructs a validated estimand from its five ICH E9(R1) attributes. The
returned object is a fully-named list (class "estimand") that
round-trips through the estimand slot of an
AnalysisResult.
defineEstimand( treatment, population, endpoint, intercurrent = list(event = NA_character_, strategy = "treatment-policy"), summary_measure = "difference in means" )defineEstimand( treatment, population, endpoint, intercurrent = list(event = NA_character_, strategy = "treatment-policy"), summary_measure = "difference in means" )
treatment |
The treatment condition(s) being compared. |
population |
The target population. |
endpoint |
The endpoint / variable of interest. |
intercurrent |
A named list with |
summary_measure |
The population-level summary (e.g.
|
An "estimand" object (a named list of the five attributes plus
the strategy's analysis recipe).
ICH E9(R1) addendum on estimands and sensitivity analysis (2019).
analyseEstimand(), multipleImputation()
defineEstimand( treatment = "rehab protocol A vs B", population = "post-stroke", endpoint = "6-month gait speed", intercurrent = list(event = "treatment discontinuation", strategy = "treatment-policy"), summary_measure = "difference in means")defineEstimand( treatment = "rehab protocol A vs B", population = "post-stroke", endpoint = "6-month gait speed", intercurrent = list(event = "treatment discontinuation", strategy = "treatment-policy"), summary_measure = "difference in means")
Fits an outcome (or change from baseline) as a function of a continuous / ordinal therapy dose, with optional confounder adjustment, subject random intercepts and a non-linear dose form. Returns the fitted dose effect, a predicted dose-response curve and the coefficient table.
doseResponse( data, outcome, dose, covariates = NULL, subject = NULL, baseline = NULL, form = c("linear", "log", "spline", "emax"), df = 3, n_grid = 50 )doseResponse( data, outcome, dose, covariates = NULL, subject = NULL, baseline = NULL, form = c("linear", "log", "spline", "emax"), df = 3, n_grid = 50 )
data |
A |
outcome |
Name of the outcome column. |
dose |
Name of the therapy-dose column (hours / repetitions / sessions). |
covariates |
Optional character vector of confounder columns to adjust for (added linearly). |
subject |
Optional subject-id column; when given, a random intercept per
subject is fit via |
baseline |
Optional baseline column; when given it is added as an adjustment covariate (ANCOVA-style change analysis). |
form |
Dose form: |
df |
Spline degrees of freedom for |
n_grid |
Number of dose points in the predicted curve (default 50). |
A PhysioExperiment::AnalysisResult of type = "dose_response": estimate
is the dose effect (slope, spline coefficients, or E0/Emax/ED50);
result holds coefficients, the curve (dose, predicted), the fit
and the form.
Ruberg SJ (1995) Dose response studies. J Biopharm Stat 5:1-14.
fitMixedModel(), proportionalRecoveryRule()
set.seed(1) d <- data.frame(dose = rep(c(0, 5, 10, 20), each = 15)) d$change <- 0.4 * d$dose + rnorm(nrow(d), 0, 2) dr <- doseResponse(d, outcome = "change", dose = "dose", form = "linear") PhysioExperiment::resultValue(dr)$coefficientsset.seed(1) d <- data.frame(dose = rep(c(0, 5, 10, 20), each = 15)) d$change <- 0.4 * d$dose + rnorm(nrow(d), 0, 2) dr <- doseResponse(d, outcome = "change", dose = "dose", form = "linear") PhysioExperiment::resultValue(dr)$coefficients
Valid ICH E9(R1) intercurrent-event strategies
ESTIMAND_STRATEGIESESTIMAND_STRATEGIES
An object of class character of length 5.
Computes estimated marginal (least-squares) means from a fitted mixed or MMRM
model with emmeans, optionally with a follow-up contrast, returning the
marginal-mean table (and contrasts) in an AnalysisResult.
estimatedMarginalMeans(model, specs, contrasts = NULL, level = 0.95, ...)estimatedMarginalMeans(model, specs, contrasts = NULL, level = 0.95, ...)
model |
An |
specs |
A formula or character spec for the marginal means, e.g.
|
contrasts |
Optional contrast method passed to |
level |
Confidence level for the intervals (default 0.95). |
... |
Further arguments forwarded to |
An AnalysisResult (type = "emmeans") with result$emmeans (the
marginal-mean table) and, when requested, result$contrasts.
Lenth RV (2024). emmeans: Estimated Marginal Means. R package.
pairwiseContrasts(), fitMMRM()
if (requireNamespace("mmrm", quietly = TRUE) && requireNamespace("emmeans", quietly = TRUE)) { data(fev_data, package = "mmrm") fit <- fitMMRM(fev_data, "FEV1", "ARMCD", "AVISIT", "USUBJID", covariates = c("RACE", "SEX")) estimatedMarginalMeans(fit, ~ ARMCD | AVISIT) }if (requireNamespace("mmrm", quietly = TRUE) && requireNamespace("emmeans", quietly = TRUE)) { data(fev_data, package = "mmrm") fit <- fitMMRM(fev_data, "FEV1", "ARMCD", "AVISIT", "USUBJID", covariates = c("RACE", "SEX")) estimatedMarginalMeans(fit, ~ ARMCD | AVISIT) }
Fits a linear mixed model with lme4::lmer and returns tidy fixed-effect and
random-effect tables in an AnalysisResult. Denominator degrees of freedom
(and the associated t tests) use lmerTest when it is installed - either the
Satterthwaite or the Kenward-Roger approximation.
fitMixedModel( data, formula, random = NULL, method = c("REML", "ML"), df = c("satterthwaite", "kenward-roger") )fitMixedModel( data, formula, random = NULL, method = c("REML", "ML"), df = c("satterthwaite", "kenward-roger") )
data |
A data.frame in long format. |
formula |
A model formula. It may carry the random-effects terms
directly (e.g. |
random |
Optional one-sided random-effects formula (e.g. |
method |
|
df |
Denominator-df method for the fixed-effect tests: |
An AnalysisResult (type = "mixed_model") whose estimate is the
fixed-effect coefficient vector, with result$fixed (the tidy fixed-effect
table), result$random (variance components), result$sigma, result$fit
(the fitted model, for estimatedMarginalMeans()) and result$formula.
Bates D, Maechler M, Bolker B, Walker S (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1). Kuznetsova A, Brockhoff PB, Christensen RHB (2017). lmerTest package. Journal of Statistical Software, 82(13).
fitMMRM(), estimatedMarginalMeans()
if (requireNamespace("lme4", quietly = TRUE)) { data(sleepstudy, package = "lme4") fitMixedModel(sleepstudy, Reaction ~ Days + (Days | Subject)) }if (requireNamespace("lme4", quietly = TRUE)) { data(sleepstudy, package = "lme4") fitMixedModel(sleepstudy, Reaction ~ Days + (Days | Subject)) }
Fits the Mallinckrodt MMRM for a longitudinal randomised trial: fixed effects
for treatment, categorical time and their interaction (plus any covariates),
with a within-subject covariance over the repeated time points. Uses the
mmrm package when available - with Kenward-Roger or Satterthwaite adjusted
degrees of freedom - and falls back to nlme::gls otherwise. The fallback
reproduces the mmrm fixed effects for the unstructured covariance (its
primary use) and approximates the homogeneous cs/ar1/toeplitz
structures; it reports between-within (containment) degrees of freedom and
emits a message, since Kenward-Roger / Satterthwaite df need mmrm.
fitMMRM( data, response, treatment, time, subject, covariates = NULL, covariance = c("unstructured", "ar1", "compound-symmetry", "toeplitz"), df = c("kenward-roger", "satterthwaite") )fitMMRM( data, response, treatment, time, subject, covariates = NULL, covariance = c("unstructured", "ar1", "compound-symmetry", "toeplitz"), df = c("kenward-roger", "satterthwaite") )
data |
A long-format data.frame with one row per subject-time. |
response, treatment, time, subject
|
Column names of the outcome, the treatment factor, the categorical time factor and the subject id. |
covariates |
Optional character vector of extra fixed-effect covariate columns (e.g. baseline value, stratifiers). |
covariance |
Within-subject covariance: |
df |
Denominator-df method: |
An AnalysisResult (type = "mmrm") whose estimate is the
fixed-effect coefficient vector, with result$coefficients (estimate, SE,
df, t, p), result$backend ("mmrm" or "gls"), result$fit and
result$formula.
Mallinckrodt CH et al. (2008). Recommendations for the primary analysis of continuous endpoints in longitudinal clinical trials. Drug Information Journal, 42. Sabanes Bove D et al. (2023). mmrm: Mixed Models for Repeated Measures. R package.
fitMixedModel(), estimatedMarginalMeans()
if (requireNamespace("mmrm", quietly = TRUE)) { data(fev_data, package = "mmrm") fitMMRM(fev_data, "FEV1", "ARMCD", "AVISIT", "USUBJID", covariates = c("RACE", "SEX")) }if (requireNamespace("mmrm", quietly = TRUE)) { data(fev_data, package = "mmrm") fitMMRM(fev_data, "FEV1", "ARMCD", "AVISIT", "USUBJID", covariates = c("RACE", "SEX")) }
Fits a finite mixture of polynomial growth trajectories to panel data (each subject belongs to one latent class), selecting the number of classes by BIC. Returns per-subject class assignments, posterior class probabilities, and the normalised entropy of the classification.
latentClassGrowth( data, subject, time, outcome, n_classes = 1:3, degree = 1, seed = NULL )latentClassGrowth( data, subject, time, outcome, n_classes = 1:3, degree = 1, seed = NULL )
data |
A long-format data frame. |
subject, time, outcome
|
Column names for the grouping factor, the time variable, and the response. |
n_classes |
Integer vector of class counts to compare (default
|
degree |
Polynomial degree of the per-class growth curve (default 1, linear). |
seed |
Optional RNG seed (flexmix uses random starts). |
An AnalysisResult (type "latent_class_growth") whose
result holds the selected n_classes, the bic table, the
per-subject assignment, the posterior probabilities, the
entropy, and the fitted flexmix model.
Nagin 2005 (group-based trajectory modelling); Gruen & Leisch 2008 (flexmix).
recoveryTrajectoryLME(), proportionalRecoveryRule()
set.seed(1) df <- do.call(rbind, lapply(1:40, function(s) { fast <- s <= 20; t <- 0:6 y <- (if (fast) 5 * t else 0.5 * t) + rnorm(7, 0, 1.5) data.frame(subject = s, time = t, y = y) })) latentClassGrowth(df, "subject", "time", "y", n_classes = 1:3)set.seed(1) df <- do.call(rbind, lapply(1:40, function(s) { fast <- s <= 20; t <- 0:6 y <- (if (fast) 5 * t else 0.5 * t) + rnorm(7, 0, 1.5) data.frame(subject = s, time = t, y = y) })) latentClassGrowth(df, "subject", "time", "y", n_classes = 1:3)
Fits a Cox model for the time-to-milestone endpoint against a grouping variable and returns the group hazard ratio - the relative rate of attaining the milestone - with its confidence interval.
milestoneHazard(milestone, group, conf_level = 0.95)milestoneHazard(milestone, group, conf_level = 0.95)
milestone |
A time-to-milestone data.frame from |
group |
Column name of the grouping variable. |
conf_level |
Confidence level (default 0.95). |
An AnalysisResult (type = "milestone_hazard") as coxModel().
Wraps mice::mice for missing-at-random multiple imputation (the
default), or delegates to rbmi for reference-based (jump-to-reference)
imputation. Imputations are reproducible under a fixed seed.
multipleImputation( data, method = "pmm", m = 5, predictors = NULL, seed = NULL, reference_based = FALSE, ... )multipleImputation( data, method = "pmm", m = 5, predictors = NULL, seed = NULL, reference_based = FALSE, ... )
data |
A data frame with missing values ( |
method |
Imputation method passed to mice (default |
m |
Number of imputations (default 5). |
predictors |
Optional: a character vector of variables to use as
predictors for every imputed column, or a full mice
|
seed |
Optional integer seed for reproducibility. |
reference_based |
Logical; if |
... |
Further arguments passed to |
A mice "mids" object (MAR path) with m and
method recorded.
van Buuren & Groothuis-Oudshoorn 2011 (mice); Carpenter et al. 2013 (reference-based MI); rbmi.
poolEstimates(), analyseEstimand()
if (requireNamespace("mice", quietly = TRUE)) { imp <- multipleImputation(mice::nhanes, m = 5, seed = 1) }if (requireNamespace("mice", quietly = TRUE)) { imp <- multipleImputation(mice::nhanes, m = 5, seed = 1) }
Convenience wrapper computing all pairwise differences of the estimated
marginal means for specs.
pairwiseContrasts(model, specs, adjust = "tukey", level = 0.95, ...)pairwiseContrasts(model, specs, adjust = "tukey", level = 0.95, ...)
model |
An |
specs |
A formula or character spec for the marginal means, e.g.
|
adjust |
Multiplicity adjustment passed to |
level |
Confidence level for the intervals (default 0.95). |
... |
Further arguments forwarded to |
An AnalysisResult (type = "emmeans_contrasts") with
result$contrasts.
Flattens an epoched or longitudinal PhysioExperiment into a long
data.frame(subject, time, channel, value) suitable for fitMixedModel(),
fitMMRM() and the marginal-means helpers.
peToLong(pe, assay = NULL, value_name = "value", subjects = NULL)peToLong(pe, assay = NULL, value_name = "value", subjects = NULL)
pe |
A |
assay |
Assay name (default the first). |
value_name |
Name for the value column (default |
subjects |
Optional character/factor labels for the third dimension
(subjects/trials); defaults to the |
A long data.frame with columns subject, time, channel and the
value column.
Combines the coefficient estimates from a list of models fitted to multiply
imputed data sets using Rubin's rules with the Barnard-Rubin degrees-of-
freedom adjustment. The pooled estimate, total variance, df, confidence
interval, and fraction of missing information reproduce mice::pool.
poolEstimates(fits, conf_level = 0.95)poolEstimates(fits, conf_level = 0.95)
fits |
A list of fitted models (each with |
conf_level |
Confidence level for the pooled interval (default 0.95). |
An AnalysisResult (type "pooled_estimates") whose
result$estimates is a data frame of pooled estimate,
std.error, df, statistic, p.value,
conf.low, conf.high, and fmi per term.
Rubin 1987; Barnard & Rubin 1999. mice::pool.
multipleImputation(), analyseEstimand()
set.seed(1) fits <- lapply(1:5, function(i) lm(mpg ~ hp + wt, data = mtcars[sample(nrow(mtcars), replace = TRUE), ])) poolEstimates(fits)set.seed(1) fits <- lapply(1:5, function(i) lm(mpg ~ hp + wt, data = mtcars[sample(nrow(mtcars), replace = TRUE), ])) poolEstimates(fits)
Tests whether recovery follows the proportional recovery rule
(default
), classifies fitters vs non-fitters, handles ceiling effects,
and runs two diagnostics against the known statistical artefact: a
mixture-vs-single-line comparison and a regression-to-the-mean (label-shuffle)
control.
proportionalRecoveryRule( initial, followup, max_score, expected_prop = 0.7, fitter_tol = 0.2, ceiling_frac = 0.9, n_shuffle = 999, seed = NULL )proportionalRecoveryRule( initial, followup, max_score, expected_prop = 0.7, fitter_tol = 0.2, ceiling_frac = 0.9, n_shuffle = 999, seed = NULL )
initial |
Numeric baseline scores. |
followup |
Numeric follow-up scores (same length). |
max_score |
Scalar maximum attainable score (the recovery ceiling). |
expected_prop |
Expected recovery proportion (default 0.7). |
fitter_tol |
Half-width around |
ceiling_frac |
Subjects with |
n_shuffle |
Label-shuffles for the regression-to-the-mean control (default 999). |
seed |
Optional RNG seed. |
A genuine proportional rule implies the follow-up regresses on the initial
score with slope ; a pure regression-to-the-mean / coupling
artefact instead reproduces the -vs-potential slope when the
follow-up scores are shuffled across subjects. When the observed slope is not
distinguishable from that shuffled null, artefact_suspected is set.
An AnalysisResult (type "proportional_recovery") whose
result holds the fitted slope (+ CI and test vs
expected_prop), per-subject classification, the
mixture_preferred flag, and artefact_suspected.
Prabhakaran 2008; Winters 2015; Hawe, Scott & Dukelow 2019 (proportional-recovery artefact).
recoveryTrajectoryLME(), latentClassGrowth()
set.seed(1) init <- runif(60, 5, 55) fu <- init + 0.7 * (66 - init) + rnorm(60, 0, 3) # genuine 0.7 rule proportionalRecoveryRule(init, fu, max_score = 66)set.seed(1) init <- runif(60, 5, 55) fu <- init + 0.7 * (66 - init) + rnorm(60, 0, 3) # genuine 0.7 rule proportionalRecoveryRule(init, fu, max_score = 66)
Fits a nonlinear mixed-effects recovery curve to longitudinal panel data with subject-level random effects (partial pooling), and returns the population fixed effects together with per-subject predicted asymptote, rate, and time-to-90%-recovery.
recoveryTrajectoryLME( data, subject, time, outcome, model = c("exponential", "asymptotic", "logistic") )recoveryTrajectoryLME( data, subject, time, outcome, model = c("exponential", "asymptotic", "logistic") )
data |
A long-format data frame. |
subject, time, outcome
|
Column names (character) for the grouping factor, the time variable, and the response. |
model |
|
Models: "exponential"/"asymptotic" use the self-starting
asymptotic curve (SSasymp);
"logistic" uses SSlogis.
An AnalysisResult (type "recovery_trajectory") whose
result holds the fixed_effects, random_effects, a
per-subject data frame (asymptote, rate, time_to_90),
and the fitted nlme object.
Pinheiro & Bates 2000 (nlme); Lindstrom & Bates 1990.
proportionalRecoveryRule(), latentClassGrowth()
set.seed(1) df <- do.call(rbind, lapply(1:12, function(s) { A <- 60 + rnorm(1, 0, 5); rate <- 0.3 * exp(rnorm(1, 0, 0.2)) t <- 0:8; data.frame(subject = s, time = t, y = A * (1 - exp(-rate * t)) + rnorm(9, 0, 2)) })) fit <- recoveryTrajectoryLME(df, "subject", "time", "y")set.seed(1) df <- do.call(rbind, lapply(1:12, function(s) { A <- 60 + rnorm(1, 0, 5); rate <- 0.3 * exp(rnorm(1, 0, 0.2)) t <- 0:8; data.frame(subject = s, time = t, y = A * (1 - exp(-rate * t)) + rnorm(9, 0, 2)) })) fit <- recoveryTrajectoryLME(df, "subject", "time", "y")
Require an optional backend or stop with guidance
requireBackend(pkg)requireBackend(pkg)
pkg |
Backend package name. |
Invisibly TRUE if available; otherwise an error.
Orchestrates the SCED analysis of a reversal design: it separates the phases
(by default the four phases A1, B1, A2, B2), and for each intervention
contrast (B1 vs A1 and B2 vs A2) computes the non-overlap effect
sizes, the 2-SD-band decision and per-phase celeration lines.
scedABAB( x, phases = c("A1", "B1", "A2", "B2"), contrasts = list(c("A1", "B1"), c("A2", "B2")), improvement = c("increase", "decrease"), value = NULL, phase = NULL, time = NULL, assay = NULL, channel = 1L )scedABAB( x, phases = c("A1", "B1", "A2", "B2"), contrasts = list(c("A1", "B1"), c("A2", "B2")), improvement = c("increase", "decrease"), value = NULL, phase = NULL, time = NULL, assay = NULL, channel = 1L )
x |
A |
phases |
Ordered character vector of the phase labels (default
|
contrasts |
A list of length-2 |
improvement |
|
value, phase, time, assay, channel
|
Column/assay selectors passed to the
extractor when |
An object of class sced_abab (a list) with phase_data,
contrasts (per contrast: NAP/Tau/Tau-U/PND/PEM results, the 2-SD decision
and both phase celeration lines) and a summary data.frame, with print
and plot methods.
scedNAP(), scedTwoSDBand(), scedCelerationLine()
df <- data.frame( value = c(10, 11, 9, 18, 19, 20, 10, 12, 11, 21, 22, 20), phase = rep(c("A1", "B1", "A2", "B2"), each = 3)) res <- scedABAB(df) resdf <- data.frame( value = c(10, 11, 9, 18, 19, 20, 10, 12, 11, 21, 22, 20), phase = rep(c("A1", "B1", "A2", "B2"), each = 3)) res <- scedABAB(df) res
Fits a trend line to one phase, either by the White & Haring split-middle method (the standard SCED hand method) or by ordinary least squares.
scedCelerationLine( value, time = NULL, method = c("split_middle", "ols"), ratio_period = 7 )scedCelerationLine( value, time = NULL, method = c("split_middle", "ols"), ratio_period = 7 )
value |
Numeric phase observations, or pass |
time |
Optional numeric session index (default |
method |
|
ratio_period |
Sessions spanned by the reported celeration ratio (default 7, a weekly ratio for daily data). |
An AnalysisResult (type = "sced_celeration") whose estimate is
the slope (celeration per session), with result$intercept,
result$fitted, result$bounce (the ratio of the largest positive to the
largest negative residual magnitude) and result$celeration_ratio (the
multiplicative change over ratio_period sessions).
White OR, Haring NG (1980). Exceptional Teaching, 2nd ed. Columbus, OH: Merrill.
scedCelerationLine(c(4, 6, 5, 8, 7, 10, 9))scedCelerationLine(c(4, 6, 5, 8, 7, 10, 9))
The probability that a randomly chosen intervention observation exceeds a randomly chosen baseline observation, with ties counted as one half (Parker & Vannest 2009) - equivalent to the area under the ROC curve. The confidence interval is the score interval of Newcombe (2006) and the standard error the unbiased estimator, both as in SingleCaseES.
scedNAP( A_data, B_data = NULL, improvement = c("increase", "decrease"), confidence = 0.95 )scedNAP( A_data, B_data = NULL, improvement = c("increase", "decrease"), confidence = 0.95 )
A_data, B_data
|
Numeric baseline (A) and intervention (B) observations,
or pass a data.frame to the phase orchestrator |
improvement |
|
confidence |
Confidence level for the interval (default 0.95). |
An AnalysisResult (type = "sced_nap") whose estimate is the NAP,
with an analytic uncertainty interval and result$se, result$p_value.
Parker RI, Vannest KJ (2009). An improved effect size for single-case research: non-overlap of all pairs. Behavior Therapy, 40(4), 357-367. Newcombe RG (2006). Confidence intervals for an effect size measure based on the Mann-Whitney statistic. Statistics in Medicine, 25.
scedNAP(c(20, 20, 26, 25), c(28, 25, 30, 29))scedNAP(c(20, 20, 26, 25), c(28, 25, 30, 29))
The proportion of intervention-phase observations beyond the baseline median, counting ties as one half (Ma 2006).
scedPEM(A_data, B_data = NULL, improvement = c("increase", "decrease"))scedPEM(A_data, B_data = NULL, improvement = c("increase", "decrease"))
A_data, B_data
|
Numeric baseline (A) and intervention (B) observations,
or pass a data.frame to the phase orchestrator |
improvement |
|
An AnalysisResult (type = "sced_pem") whose estimate is the PEM.
Ma H-H (2006). An alternative method for quantitative synthesis of single-subject research: percentage of data points exceeding the median. Behavior Modification, 30(5), 598-617.
scedPEM(c(20, 20, 26, 25), c(28, 25, 30, 29))scedPEM(c(20, 20, 26, 25), c(28, 25, 30, 29))
The proportion of intervention-phase observations that exceed the most extreme baseline observation (Scruggs, Mastropieri & Casto 1987).
scedPND(A_data, B_data = NULL, improvement = c("increase", "decrease"))scedPND(A_data, B_data = NULL, improvement = c("increase", "decrease"))
A_data, B_data
|
Numeric baseline (A) and intervention (B) observations,
or pass a data.frame to the phase orchestrator |
improvement |
|
An AnalysisResult (type = "sced_pnd") whose estimate is the PND
in [0, 1].
Scruggs TE, Mastropieri MA, Casto G (1987). The quantitative synthesis of single-subject research. Remedial and Special Education, 8(2).
scedPEM(), scedNAP(), scedTauU()
scedPND(c(20, 20, 26, 25), c(28, 25, 30, 29))scedPND(c(20, 20, 26, 25), c(28, 25, 30, 29))
The rank-correlation non-overlap effect size,
(Parker et al. 2011), with the SE and interval derived from scedNAP().
scedTau( A_data, B_data = NULL, improvement = c("increase", "decrease"), confidence = 0.95 )scedTau( A_data, B_data = NULL, improvement = c("increase", "decrease"), confidence = 0.95 )
A_data, B_data
|
Numeric baseline (A) and intervention (B) observations,
or pass a data.frame to the phase orchestrator |
improvement |
|
confidence |
Confidence level for the interval (default 0.95). |
An AnalysisResult (type = "sced_tau") whose estimate is Tau in
[-1, 1].
scedTau(c(20, 20, 26, 25), c(28, 25, 30, 29))scedTau(c(20, 20, 26, 25), c(28, 25, 30, 29))
Parker, Vannest, Davis & Sauber's (2011) Tau-U combines the A-vs-B non-overlap with a correction for baseline (phase A) trend. Two reference conventions differ on the denominator and this function exposes both:
scedTauU( A_data, B_data = NULL, improvement = c("increase", "decrease"), method = c("parker", "scan"), confidence = 0.95 )scedTauU( A_data, B_data = NULL, improvement = c("increase", "decrease"), method = c("parker", "scan"), confidence = 0.95 )
A_data, B_data
|
Numeric baseline (A) and intervention (B) observations,
or pass a data.frame to the phase orchestrator |
improvement |
|
method |
Denominator convention, |
confidence |
Confidence level for the interval (default 0.95). |
method = "parker" (default)matches SingleCaseES::Tau_U:
, where is the A-vs-B
Kendall S and the within-baseline Kendall S. Can exceed
1 when a baseline trend runs counter to improvement.
method = "scan"matches scan::tau_u()'s "A vs. B - Trend A" row:
the same S numerator over a tie-corrected pair count , giving a
Kendall tau-b in [-1, 1].
The Z statistic and p-value come from the tie-corrected variance of the
Kendall S numerator and are identical under both methods (they depend on S,
not the denominator), matching scan's Z and p.
An AnalysisResult (type = "sced_tau_u") whose estimate is Tau-U,
with result$S, result$se_S, result$z, result$p_value and the
component S values.
Parker RI, Vannest KJ, Davis JL, Sauber SB (2011). Combining nonoverlap and trend for single-case research: Tau-U. Behavior Therapy, 42(2), 284-299.
scedTauU(c(20, 20, 26, 25, 22, 23), c(28, 25, 24, 27, 30, 30, 29, 28))scedTauU(c(20, 20, 26, 25, 22, 23), c(28, 25, 24, 27, 30, 30, 29, 28))
Draws the Nelson (1984) 2-SD band from the baseline (A) phase - the baseline
mean plus and minus two baseline standard deviations - and flags the
intervention (B) phase for a systematic shift: the rule fires when a run of
consecutive or more successive B points falls on the same side beyond the
band.
scedTwoSDBand( A_data, B_data = NULL, improvement = c("increase", "decrease"), k = 2, consecutive = 2L )scedTwoSDBand( A_data, B_data = NULL, improvement = c("increase", "decrease"), k = 2, consecutive = 2L )
A_data, B_data
|
Numeric baseline and intervention observations. |
improvement |
|
k |
Band half-width in baseline SDs (default 2). |
consecutive |
Number of successive out-of-band points that trigger the flag (default 2, Nelson's rule). |
An AnalysisResult (type = "sced_2sd") whose estimate is the
logical decision, with result$mean, result$sd, result$upper,
result$lower, result$outside (per-B-point side: -1/0/+1) and
result$first_run_at (index of the first triggering point, or NA).
Nelson LS (1984). The Shewhart control chart - tests for special causes. Journal of Quality Technology, 16(4), 237-239. Applied to SCED in Gast & Ledford (2014).
scedCelerationLine(), scedABAB()
scedTwoSDBand(c(10, 12, 11, 9, 10), c(15, 16, 17, 16, 18))scedTwoSDBand(c(10, 12, 11, 9, 10), c(15, 16, 17, 16, 18))
Fits a Kaplan-Meier estimator with survival::survfit, optionally stratified
by a grouping variable, and returns the median and quartile attainment times
with their confidence intervals.
survivalFit( data, time = "time", event = "event", group = NULL, conf_level = 0.95 )survivalFit( data, time = "time", event = "event", group = NULL, conf_level = 0.95 )
data |
A data.frame with the survival columns. |
time, event
|
Column names of the follow-up time and the event indicator (1 = event, 0 = censored). |
group |
Optional grouping-variable column name for stratified curves. |
conf_level |
Confidence level (default 0.95). |
An object of class physio_km (a list) with fit (the survfit),
quantiles (median / Q1 / Q3 attainment times with CIs, per group) and
summary table, plus print and plot methods.
Kaplan EL, Meier P (1958). Nonparametric estimation from incomplete observations. JASA, 53(282), 457-481.
if (requireNamespace("survival", quietly = TRUE)) { survivalFit(survival::lung, "time", "status", group = "sex") }if (requireNamespace("survival", quietly = TRUE)) { survivalFit(survival::lung, "time", "status", group = "sex") }
Applies eligibility at the declared time zero and estimates either an intention-to-treat contrast or a per-protocol clone-censor-weight contrast. The result remains conditional on consistency, exchangeability, positivity, correct treatment/censoring/outcome models, and a correctly aligned time zero. A Cox hazard ratio is non-collapsible and is not a marginal risk ratio.
targetTrialEmulate( data, protocol, id, time, treatment, outcome, event = NULL, baseline_covariates = NULL, time_varying_covariates = NULL, estimand = c("per_protocol", "intention_to_treat"), weight_backend = c("internal", "WeightIt", "ipw"), stabilized = TRUE, numerator_covariates = NULL, weight_truncation = c(0.01, 0.99), max_weight = Inf, conf.level = 0.95 )targetTrialEmulate( data, protocol, id, time, treatment, outcome, event = NULL, baseline_covariates = NULL, time_varying_covariates = NULL, estimand = c("per_protocol", "intention_to_treat"), weight_backend = c("internal", "WeightIt", "ipw"), stabilized = TRUE, numerator_covariates = NULL, weight_truncation = c(0.01, 0.99), max_weight = Inf, conf.level = 0.95 )
data |
Long person-period data. |
protocol |
A |
id, time, treatment, outcome
|
Column names. For survival outcomes,
|
event |
Optional binary event-status column. |
baseline_covariates, time_varying_covariates
|
Declared predictors for treatment/adherence models. |
estimand |
|
weight_backend |
Propensity backend. Optional backends are validation paths and never silent fallbacks. |
stabilized |
Use time-based numerator probabilities. |
numerator_covariates |
Additional numerator-model covariates. |
weight_truncation |
Lower and upper quantiles applied after raw cumulative weights are retained. |
max_weight |
Additional positive upper cap applied after quantile truncation. |
conf.level |
Confidence level. |
An AnalysisResult containing effects, models, clone rows, both raw
and analysis weights, diagnostics, censoring counts, and provenance.
Defines the seven protocol components that must be fixed before emulating a
target trial. time_zero is a separate alignment rule and is never inferred
from the observed data. Static strategies are scalar sustained treatment
values. Dynamic strategy functions receive the current row as data and
the participant's rows through the current time as history; they must not
access later rows.
targetTrialProtocol( eligibility, treatment_strategies, assignment, time_zero, follow_up, outcome, causal_contrast, analysis_plan, protocol_id = NULL, version = "1.0.0" )targetTrialProtocol( eligibility, treatment_strategies, assignment, time_zero, follow_up, outcome, causal_contrast, analysis_plan, protocol_id = NULL, version = "1.0.0" )
eligibility |
Eligibility rule. Use a function of baseline data or a one-sided formula for an executable rule. |
treatment_strategies |
Named list of at least two static scalar or dynamic function strategies. |
assignment, outcome, causal_contrast, analysis_plan
|
Non-empty, serializable protocol declarations. |
time_zero |
Explicit scalar time value or alignment rule. |
follow_up |
Non-empty follow-up declaration. A finite positive scalar
is used as the fixed horizon by |
protocol_id |
Optional stable identifier. When omitted, one is derived from the protocol declarations. |
version |
Semantic protocol version. |
An object of class target_trial_protocol.
protocol <- targetTrialProtocol( eligibility = function(data) data$eligible, treatment_strategies = list(never = 0, always = 1), assignment = "cloning at eligibility", time_zero = 0, follow_up = 12, outcome = "binary recovery by week 12", causal_contrast = "always versus never", analysis_plan = "clone-censor-weight" ) protocolprotocol <- targetTrialProtocol( eligibility = function(data) data$eligible, treatment_strategies = list(never = 0, always = 1), assignment = "cloning at eligibility", time_zero = 0, follow_up = 12, outcome = "binary recovery by week 12", causal_contrast = "always versus never", analysis_plan = "clone-censor-weight" ) protocol
From a long-format longitudinal series, finds for each subject the first time
its outcome crosses a clinical threshold, producing a right-censored
(time, event) pair per subject: event = 1 at the first crossing, or
event = 0 (right-censored) at the last observation - or at censor_at - for
subjects who never attain the milestone.
timeToMilestone( data, id, time, value, threshold, direction = c("increase", "decrease"), censor_at = NULL )timeToMilestone( data, id, time, value, threshold, direction = c("increase", "decrease"), censor_at = NULL )
data |
A long-format data.frame with one row per subject-visit. |
id, time, value
|
Column names of the subject id, the visit time and the outcome value. |
threshold |
The clinical milestone value. |
direction |
|
censor_at |
Optional common administrative censoring time for non-attainers; defaults to each subject's last observation time. |
A data.frame with columns id, time, event (1 attained /
0 censored) and attained, carrying a Surv object as the "surv"
attribute.
survivalFit(), milestoneHazard()
d <- data.frame( id = rep(c("a", "b"), each = 3), time = rep(1:3, 2), value = c(10, 30, 55, 12, 20, 28)) timeToMilestone(d, "id", "time", "value", threshold = 50)d <- data.frame( id = rep(c("a", "b"), each = 3), time = rep(1:3, 2), value = c(10, 30, 55, 12, 20, 28)) timeToMilestone(d, "id", "time", "value", threshold = 50)