Methodology

Baseline Characteristics / Descriptive Statistics

Purpose and scope

The descriptive statistics module summarizes metric and categorical variables for the complete dataset and, optionally, for every non-missing level of one group variable. It is intended for baseline-characteristics tables, exploratory summaries, data-quality checks, and transparent reporting before group-comparison analyses. The module focuses on descriptive summaries and additionally reports an inferential 95% confidence interval for the mean. It does not calculate group-comparison p-values or standardized mean differences.

Missing values and numeric parsing

Empty cells, non-finite numbers, and the case-insensitive tokens NA, N/A, NaN, null, none, and missing are treated as missing. Numeric parsing respects the selected point or comma decimal separator. In automatic mode, when commas occur without points, the final comma is treated as the decimal mark and preceding commas as grouping separators. If both comma and point occur, the final separator is treated as the decimal mark. A non-missing value that cannot be parsed as a complete number is reported separately as an invalid numeric value and is not silently coerced.

Grouping

With No group, one summary is produced for the complete cohort. When a group variable is selected, Sigma reports Overall plus one column for every observed non-missing group level. Rows with a missing group value remain part of Overall but are excluded from all group-specific summaries. The selected group variable cannot simultaneously be analyzed as a parameter.

Automatic variable type and manual override

Type detection is performed on all non-missing values of a parameter. A variable with at most two distinct values is classified as categorical, including numeric binary codes such as 0/1. Otherwise, a variable is classified as metric when at least 80% of its non-missing values are parseable numbers; all remaining variables are categorical. A high-cardinality warning is shown when at least 20 distinct levels comprise 90% or more of the non-missing observations.

Every variable can be overridden as Metric or Categorical. The option Show all statistics for all parameters calculates both numeric statistics and frequency distributions for every selected variable. For text-only variables, numeric quantities are explicitly reported as N/A.

Metric statistics

Let \(x_1,\ldots,x_n\) denote the valid numeric observations after excluding missing and invalid values. Sigma reports total observations, valid numeric n, missing and invalid counts and percentages, mean, sample variance, sample standard deviation (SD), standard error of the mean (SEM), a 95% confidence interval for the mean, median, quartiles, interquartile range (IQR), minimum, maximum, range, adjusted skewness, adjusted excess kurtosis, and user-selected percentiles.

\[ \bar{x}=\frac{1}{n}\sum_{i=1}^{n}x_i,\qquad s^2=\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\bar{x})^2,\qquad SEM=\frac{s}{\sqrt{n}}. \]

Variance, SD, SEM, and the mean confidence interval are undefined for \(n<2\). For \(n\ge2\), the two-sided 95% confidence interval is

\[ \bar{x}\pm t_{0.975,n-1}\cdot SEM. \]

This interval assumes independent observations. Its nominal finite-sample coverage is exact when the observations are sampled from a normal population; for non-normal data, its interpretation relies on the sampling distribution of the mean being adequately approximated by a normal distribution, which is generally more plausible in larger samples. It describes uncertainty in the population mean and should not be interpreted as an interval containing 95% of individual observations.

Percentiles and IQR

Percentiles use the same linear interpolation as R's quantile(..., type = 7). For sorted observations and probability \(p\), the zero-based position is \(h=(n-1)p\); Sigma linearly interpolates between the observations at \(\lfloor h\rfloor\) and \(\lceil h\rceil\). The first quartile, median, and third quartile correspond to \(p=0.25,0.5,0.75\), and \(IQR=Q_3-Q_1\).

Skewness and excess kurtosis

Define the central moments \(m_k=n^{-1}\sum_i(x_i-\bar{x})^k\). Sigma reports adjusted Fisher–Pearson skewness for \(n\ge3\),

\[ G_1=\frac{\sqrt{n(n-1)}}{n-2}\frac{m_3}{m_2^{3/2}}, \]

and adjusted excess kurtosis for \(n\ge4\),

\[ G_2=\frac{n-1}{(n-2)(n-3)}\left[(n+1)\left(\frac{m_4}{m_2^2}-3\right)+6\right]. \]

Both shape measures are undefined for a constant sample because \(m_2=0\).

Categorical statistics

Sigma reports total, valid, and missing observations, the number of distinct levels, all modal levels, and the frequency of every level. Level percentages use the number of non-missing observations as their denominator. The missing percentage uses all observations in the corresponding cohort as its denominator. Levels are displayed in first-observed order.

Output and limitations

The baseline table combines Mean ± SD and Median [Q1–Q3] rows for metric parameters with n (%) rows for categorical levels. Detailed tables expose all calculations, percentiles, frequencies, and type-detection diagnostics, rounded on screen according to the selected display precision. The detailed and grouped Excel sheets retain full-precision numeric values, whereas the publication-style Excel baseline table and all Word tables use the selected display rounding. When a group variable is selected, an additional publication layout uses a two-row header: every group is one horizontally merged header cell spanning all reported statistics beneath it. Metric and categorical parameters are presented in separate grouped tables for direct transfer into manuscripts. Word and Excel exports can include these grouped tables, the baseline table, numeric details and percentiles, and frequencies and data quality.

• Summaries describe the observed sample and do not establish causal or inferential group differences.
• Very small samples yield unstable distribution-shape estimates and wide confidence intervals.
• Numeric category codes may require a manual categorical override.
• High-cardinality variables such as identifiers can create very long frequency tables.

Normality Testing

Purpose and scope

The Sigma normality testing module assesses whether a univariate numeric sample is compatible with a Gaussian (normal) distribution. Normality testing is typically used to support methodological choices (e.g., whether a parametric model is plausible) and should be interpreted together with graphical diagnostics and the study context.

Important: “Failing to reject normality” (e.g., p ≥ 0.05) does not prove that the data are normal; it indicates that the observed deviations are not strong enough to reject normality at the chosen α. Conversely, with large sample sizes even small and clinically irrelevant deviations may yield p < 0.05.

Input data and preprocessing

The module accepts one or more numeric variables, which are analyzed independently. Numeric fields are parsed strictly according to the selected decimal-separator setting; missing, non-numeric, partially numeric, and non-finite values are excluded separately for each parameter. If a group variable is used, rows with a missing group value are also excluded rather than being treated as a separate group.

All tests are computed on the resulting sample x = (x_1, \ldots, x_n) of size n. Where required, the sample mean \(\bar{x}\) and sample standard deviation \(s = \sqrt{\sum_i(x_i-\bar{x})^2/(n-1)}\) are estimated from the data. In grouped analyses, Sigma reports the pooled eligible sample and each observed group level separately. If all observations are identical or a required variance or transformation is undefined, the respective test is not reported.

Implemented tests and statistical definitions

Shapiro–Wilk test

The Shapiro–Wilk test evaluates normality by measuring how closely the ordered sample resembles the expected order statistics of a normal distribution. Let \(x_{(1)} \le \cdots \le x_{(n)}\) denote the ordered sample and \(\bar{x}\) the sample mean. The test statistic is:

\[ W = \dfrac{\left(\sum_{i=1}^{n} a_i x_{(i)}\right)^2}{\sum_{i=1}^{n} (x_i - \bar{x})^2}, \]

where the coefficients \(a_i\) are derived from expected values of normal order statistics and their covariance structure. The p-value is obtained using Royston-type approximations based on transformations of \(1-W\), following the standard Shapiro–Wilk computational approach.

Sample size requirement: Sigma computes the Shapiro–Wilk test only for 3 ≤ n ≤ 5000.

Anderson–Darling test

The Anderson–Darling (AD) test is an EDF (empirical distribution function) goodness-of-fit test with relatively high sensitivity in the distribution tails. Let \(F\) denote the fitted normal CDF (using sample mean and standard deviation) and \(x_{(i)}\) the ordered observations. The base AD statistic is:

\[ A^2 = -n - \frac{1}{n}\sum_{i=1}^{n}\left((2i-1)\left[\ln F(x_{(i)}) + \ln\left(1-F(x_{(n+1-i)})\right)\right]\right). \]

Sigma reports the finite-sample adjusted statistic

\[ A^{2*} = A^2\left(1 + \frac{0.75}{n} + \frac{2.25}{n^2}\right), \]

and computes the p-value using the same Stephens-type piecewise approximation used by nortest::ad.test. For A²* ≥ 10, the approximation returns the documented extreme upper-tail value 3.7 × 10−24; it must not be extrapolated beyond its final interval.

Sample size requirements: the AD statistic can be computed for n ≥ 3. P-value reporting is based on the standard approximation and is provided for n ≥ 8.

Lilliefors-corrected Kolmogorov–Smirnov normality test

Sigma implements the one-sample Kolmogorov–Smirnov normality statistic against a fitted normal distribution, i.e. with mean and standard deviation estimated from the sample (composite null hypothesis). Let \(F\) be the fitted normal CDF and \(F_n\) the empirical CDF. The test statistic is:

\[ D = \sup_x \left|F_n(x) - F(x)\right|. \]

The p-value follows the approximation used by nortest::lillie.test. For n ≤ 100, define Dd = D and nd = n; for n > 100, define Dd = D(n/100)0.49 and nd = 100. The initial tail approximation is

\[ p_0 = \exp\left[-7.01256D_d^2(n_d+2.78019) +2.99587D_d\sqrt{n_d+2.78019}-0.122119 +\frac{0.974598}{\sqrt{n_d}}+\frac{1.67997}{n_d}\right]. \]

If p0 > 0.1, the upper range is refined using the published piecewise polynomial approximation in K = (√n − 0.01 + 0.85/√n)D. The reported value is therefore a Lilliefors p-value for a normal distribution whose mean and variance were estimated from the same sample—not the classical KS p-value for a fully specified distribution.

Sample size requirement: this test is reported for n ≥ 5.

D’Agostino–Pearson omnibus K² test

The D’Agostino–Pearson omnibus test combines skewness and kurtosis into a single normality statistic. Sigma computes biased sample skewness and Pearson kurtosis from central moments, applies the finite-sample transformations of D’Agostino and Pearson to obtain approximately standard-normal scores, and then forms:

\[ K^2 = Z_{\text{skew}}^2 + Z_{\text{kurt}}^2. \]

Under the null hypothesis of normality, \(K^2\) is evaluated against a \(\chi^2\) distribution with 2 degrees of freedom, so that p = exp(−K²/2). The implementation follows the definitions used by scipy.stats.normaltest, including the factor 0.5 in the logarithmic skewness transformation.

Sample size requirement: due to the moment-based approximations, the test is computed for n ≥ 8.

Hypotheses, significance level, and reported outputs

For each selected test, the null hypothesis is \(H_0\): the data are drawn from a normal distribution (with parameters estimated from the sample where applicable). The alternative \(H_1\) is that the distribution deviates from normality.

Sigma reports for each selected test:

• sample size n used by the test,
• the test statistic (e.g., W, D, A^{2*}, K^2),
• the corresponding p-value, if available for the implemented approximation regime.

Default interpretation uses α = 0.05 unless specified otherwise: p < α indicates evidence against normality, while p ≥ α indicates that the sample is compatible with normality at that α level.

Assumptions and limitations

• Independence: Observations are assumed independent.
• Sensitivity to n: With large n, normality tests may reject for small, practically irrelevant deviations; with small n, power can be limited.
• Discrete / tied data: Strong discreteness or many ties may affect EDF-based tests and moment-based approximations.
• Degenerate data: Constant or near-constant data can make one or more tests undefined.
• Normality is not model adequacy: Passing a normality test does not guarantee that a parametric model is appropriate; conversely, rejecting normality does not automatically invalidate parametric methods.

Warnings, edge cases, and safeguards

Sigma only reports a test if its minimum sample size requirement is met and the computation is numerically meaningful. Typical situations where results may be omitted or flagged include:

• n below the test-specific minimum,
• Shapiro–Wilk sample sizes above 5000,
• constant or near-constant data (variance approximately zero),
• invalid numeric input after preprocessing (e.g., all values missing/non-numeric).

Numerical implementation and reference equivalence

Computations are performed in JavaScript using double-precision floating point arithmetic (IEEE-754). Central moments are evaluated after shifting and scaling the observations, and the normal log-CDF is used for Anderson–Darling tail terms to avoid overflow, cancellation, and log(0) errors. Established p-value approximations are reproduced directly; very small floating-point differences can nevertheless occur because the underlying languages and elementary functions are not bitwise identical.

The production calculation core has been compared directly with these independent reference implementations:

• stats::shapiro.test for Shapiro–Wilk,
• nortest::lillie.test for Lilliefors-corrected KS normality testing,
• nortest::ad.test for Anderson–Darling normality testing, and
• scipy.stats.normaltest for D’Agostino–Pearson K².

In Sigma, the result labelled “Lilliefors” is the fitted-normal variant of the KS statistic, not the classical KS test with completely fixed distribution parameters.

References

1. Royston, J. P. "W test for normality: algorithm AS 181." Applied Statistics 31.2 (1982): 176-180.
2. Royston, P. "Remark AS R94: A Remark on Algorithm AS 181: The W-test for Normality." Applied Statistics 44.4 (1995): 547-551.
3. Stephens, Michael A. "EDF statistics for goodness of fit and some comparisons." Journal of the American Statistical Association 69.347 (1974): 730-737.
4. Lilliefors, Hubert W. "On the Kolmogorov-Smirnov test for normality with mean and variance unknown." Journal of the American Statistical Association 62.318 (1967): 399-402.
5. D’Agostino, Ralph B., Albert Belanger, and Ralph B. D’Agostino Jr. "A suggestion for using powerful and informative tests of normality." The American Statistician 44.4 (1990): 316-321.
6. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.
7. Thode, Henry C. Testing for normality. CRC Press, 2002.

Chi-square (χ²) Test (2×2) and Fisher’s Exact Test

Purpose and scope

This module analyzes associations between a binary group variable (e.g., Control vs Treatment) and one or more binary outcome variables (parameters) using a 2×2 contingency table. For each outcome, Sigma reports both an approximate large-sample test (Pearson χ² with Yates’ continuity correction) and an exact test (Fisher’s exact test), together with common effect measures (odds ratio, risk ratio, risk difference).

Input data, variable coding, and preprocessing

The analysis requires:

• a group column with exactly two observed (non-missing) levels, and
• one or more parameter columns, each with exactly two observed (non-missing) levels.

Rows with missing values in the group column or the current parameter column are excluded for that parameter (complete-case evaluation per parameter). Outcome levels are determined from these complete cases within the two selected group levels. If the group column does not have exactly two levels, or a parameter is not binary after this filtering, the corresponding analysis is not performed.

The user selects which group level is labeled Control and which is labeled Treatment. For every selected parameter, the user explicitly chooses which observed outcome level is the event; the other level becomes the reference. Effect measures (OR/RR/RD) are reported for “Treatment vs Control” and “Event vs Reference”. A visible deterministic default is proposed, but it can be changed before analysis.

Note: The direction (sign) and magnitude of OR/RR/RD depend on which category is treated as “Event” versus “Reference”. Verify the displayed event selection before interpreting or exporting effect estimates.

Contingency table notation

For each binary parameter, Sigma constructs the 2×2 table:

\[ \begin{array}{c|cc|c} & \text{Event} & \text{Reference} & \text{Total} \\ \hline \text{Treatment} & a & b & a+b \\ \text{Control} & c & d & c+d \\ \hline \text{Total} & a+c & b+d & N \end{array} \]

with total sample size \(N = a+b+c+d\).

Pearson chi-square test with Yates’ continuity correction

Under the null hypothesis of independence (no association between group and outcome), the expected counts are:

\[ E_{ij}=\frac{(\text{row total})_i\cdot(\text{column total})_j}{N}. \]

Sigma computes the Pearson chi-square statistic for a 2×2 table using Yates’ continuity correction:

\[ \chi^2_Y = \sum_{cells}\frac{\left(\max\left(0,\;|O-E|-0.5\right)\right)^2}{E}, \]

where \(O\) and \(E\) denote observed and expected counts for each cell. Degrees of freedom for a 2×2 table are \(df=1\).

The two-sided p-value is obtained from the chi-square distribution with 1 degree of freedom. For \(df=1\), Sigma uses the equivalent normal relationship:

\[ p = 2\left(1-\Phi\left(\sqrt{\chi^2_Y}\right)\right), \]

where \(\Phi\) is the standard normal CDF. Numerically, Sigma evaluates the equivalent survival function directly as \(\operatorname{erfc}(\sqrt{\chi^2_Y/2})\), avoiding upper-tail cancellation.

Non-estimable margins

Pearson’s χ² statistic is not reported when a row total or column total is zero, because at least one expected count is then zero and the statistic is not estimable. Sigma marks this condition explicitly instead of returning χ² = 0 or p = 1. Fisher’s two-sided p-value remains available when its fixed-margin distribution is defined; effect estimates are omitted when a group or outcome margin is zero.

Fisher’s exact test (two-sided)

Fisher’s exact test evaluates the null hypothesis of independence conditional on the fixed row and column margins. Sigma computes the probability of the observed table under the hypergeometric model and reports a two-sided p-value by summing probabilities of all tables with the same margins whose probabilities are less than or equal to that of the observed table (standard two-sided definition for Fisher’s exact test). Probabilities are compared and accumulated on the log scale with a small relative floating-point tolerance.

In small samples or when expected counts are low, Fisher’s exact test is often preferred because it does not rely on large-sample chi-square approximations.

Effect measures and confidence intervals

Odds ratio (OR)

The odds ratio compares the odds of the event in Treatment versus Control:

\[ OR=\frac{a/b}{c/d}=\frac{ad}{bc}. \]

If any cell count is zero while all row and column margins remain positive, Sigma applies the Haldane–Anscombe correction by adding 0.5 to all four cells before computing \(OR\) and its confidence interval. If a margin is zero, the odds ratio is not estimable and is omitted.

A 95% confidence interval is computed using the Wald method on the log scale:

\[ SE(\log OR)=\sqrt{\frac{1}{a}+\frac{1}{b}+\frac{1}{c}+\frac{1}{d}},\quad CI_{95\%}=\exp\left(\log OR \pm 1.96\cdot SE(\log OR)\right), \]

where the corrected cell counts are used if the Haldane–Anscombe adjustment was applied.

Risk ratio (RR)

The risk ratio compares event probabilities:

\[ risk_T=\frac{a}{a+b},\quad risk_C=\frac{c}{c+d},\quad RR=\frac{risk_T}{risk_C}. \]

A 95% confidence interval is computed on the log scale (Wald method) when both event counts are nonzero:

\[ SE(\log RR)=\sqrt{\left(\frac{1}{a}-\frac{1}{a+b}\right)+\left(\frac{1}{c}-\frac{1}{c+d}\right)},\quad CI_{95\%}=\exp\left(\log RR \pm 1.96\cdot SE(\log RR)\right). \]

If \(a=0\) or \(c=0\), \(\log RR\) is undefined; in this case Sigma may omit the RR CI.

Risk difference (RD)

The risk difference is:

\[ RD=risk_T-risk_C. \]

A 95% Wald confidence interval is computed as:

\[ SE(RD)=\sqrt{\frac{risk_T(1-risk_T)}{a+b}+\frac{risk_C(1-risk_C)}{c+d}},\quad CI_{95\%}=RD \pm 1.96\cdot SE(RD). \]

Low expected counts flag

Sigma flags low expected counts if at least one expected cell count is < 5. This is a common rule-of-thumb indicating that large-sample chi-square approximations may be unreliable, and that Fisher’s exact test (or alternative methods) should be emphasized.

Assumptions and limitations

• Independence: observations are assumed independent; the test is not designed for paired or clustered data.
• 2×2 only: Sigma’s χ² module evaluates only binary group vs. binary outcome (no general r×c tables).
• Direction depends on coding: OR/RR/RD depend on which outcome category is treated as “Event”. Ensure coding matches the intended interpretation.
• Multiple testing: when many parameters are tested, Sigma does not automatically apply multiplicity corrections (e.g., Bonferroni or FDR); interpretation should consider the number of tests.

Numerical implementation and equivalence to R

Calculations are performed in JavaScript using double-precision floating point arithmetic (IEEE-754). χ² tail probabilities are evaluated directly through the complementary error function, and Fisher probabilities are evaluated with log-gamma and log-sum-exp calculations. Conceptually, the procedures correspond to common reference implementations in R:

• stats::chisq.test (Pearson χ²; note that Sigma applies Yates’ correction for 2×2 tables),
• stats::fisher.test (two-sided Fisher’s exact test).

Small numerical differences between software packages may occur due to rounding and implementation details, but the underlying statistical definitions are standard.

References

1. Pearson, K. (1900). On the criterion that a given system of deviations from the probable… Philosophical Magazine 50(302): 157–175.
2. Yates, F. (1934). Contingency tables involving small numbers and the chi-square test. Journal of the Royal Statistical Society (Suppl.) 1: 217–235.
3. Fisher, R. A. (1922). On the interpretation of χ² from contingency tables, and the calculation of P. Journal of the Royal Statistical Society 85(1): 87–94.
4. Woolf, B. (1955). On estimating the relation between blood group and disease. Annals of Human Genetics 19: 251–253.
5. Haldane, J. B. S. (1955). The estimation and significance of the logarithm of a ratio of frequencies. Annals of Human Genetics 20: 309–311.
6. Anscombe, F. J. (1956). On estimating binomial response relations. Biometrika 43(3/4): 461–464.
7. Agresti, A. Categorical Data Analysis. Wiley (current edition).

Unpaired (Independent Samples) t-Test

Purpose and scope

The unpaired (independent samples) t-test compares the means of a continuous variable between two independent groups (e.g., Control vs Treatment). Sigma reports both the Welch t-test (unequal variances; default for inference) and the Student t-test with pooled variance (equal variances assumed), along with confidence intervals and standardized effect sizes.

Input data and preprocessing

The analysis requires a group column with exactly two observed (non-missing) levels and one or more numeric parameter columns analyzed independently. For each parameter, observations with missing or non-numeric entries in either the group or the parameter are excluded for that parameter (complete-case evaluation per parameter). Numeric strings are accepted only when the entire normalized value is syntactically valid under the selected decimal-separator setting. Non-finite values and partially numeric strings are treated as invalid and excluded.

The user assigns group labels Control and Treatment. Throughout the module, the mean difference is defined as:

\[ \Delta = \bar{x}_T - \bar{x}_C, \]

where \(\bar{x}_T\) and \(\bar{x}_C\) are the sample means in Treatment and Control, respectively.

Descriptive statistics

For each group and parameter, Sigma reports sample size \(n\), mean \(\bar{x}\), sample variance \(s^2\) (denominator \(n-1\)), standard deviation \(s\), standard error \(SE=s/\sqrt{n}\), and distribution summaries (min/max, and quartiles). Quartiles are computed from sorted values using a standard linear interpolation approach.

Welch t-test (unequal variances)

The Welch t-test does not assume equal variances between groups. Let \(s_T^2\) and \(s_C^2\) denote the sample variances. The standard error of the mean difference is:

\[ SE_W(\Delta)=\sqrt{\frac{s_T^2}{n_T}+\frac{s_C^2}{n_C}}. \]

The Welch t-statistic is:

\[ t_W=\frac{\Delta}{SE_W(\Delta)}. \]

Degrees of freedom are approximated using the Welch–Satterthwaite equation:

\[ df_W= \frac{\left(\frac{s_T^2}{n_T}+\frac{s_C^2}{n_C}\right)^2}{ \frac{\left(\frac{s_T^2}{n_T}\right)^2}{n_T-1}+ \frac{\left(\frac{s_C^2}{n_C}\right)^2}{n_C-1} }. \]

The two-sided p-value is computed from the Student t distribution with \(df_W\) degrees of freedom:

\[ p = 2\left(1 - F_{t,df_W}\left(|t_W|\right)\right), \]

where \(F_{t,df}\) denotes the t-distribution CDF.

Student t-test with pooled variance (equal variances assumed)

Under the assumption of equal variances, Sigma additionally computes the pooled-variance t-test. The pooled variance is:

\[ s_p^2=\frac{(n_T-1)s_T^2+(n_C-1)s_C^2}{n_T+n_C-2}, \quad s_p=\sqrt{s_p^2}. \]

The standard error of the mean difference is:

\[ SE_P(\Delta)=s_p\sqrt{\frac{1}{n_T}+\frac{1}{n_C}}, \]

and the pooled t-statistic is:

\[ t_P=\frac{\Delta}{SE_P(\Delta)}, \quad df_P=n_T+n_C-2. \]

The two-sided p-value is computed as \(p = 2(1 - F_{t,df_P}(|t_P|))\).

Interpretation: Welch’s test is generally preferred when variances differ; Sigma reports both results for transparency. In the summary output, the Welch two-sided p-value is used as the primary inference quantity.

Confidence intervals for the mean difference

Sigma reports a two-sided 95% confidence interval for \(\Delta=\bar{x}_T-\bar{x}_C\) using the relevant standard error and critical t-value:

\[ CI_{95\%} = \Delta \pm t_{1-\alpha/2,df}\cdot SE(\Delta), \quad \alpha=0.05, \]

where \(df=df_W\) and \(SE=SE_W\) for Welch, and \(df=df_P\) and \(SE=SE_P\) for the pooled-variance test.

Effect sizes

Sigma reports standardized mean differences based on the pooled standard deviation \(s_p\). Cohen’s \(d\) is computed as:

\[ d=\frac{\Delta}{s_p}. \]

Hedges’ \(g\) applies the exact gamma-function small-sample correction:

\[ J(df_P)=\frac{\Gamma(df_P/2)}{\sqrt{df_P/2}\,\Gamma((df_P-1)/2)}, \quad g = J\cdot d. \]

Assumptions and limitations

• Independence: observations are assumed independent between groups.
• Scale and distribution: the outcome variable should be continuous (or approximately continuous). Normality within each group supports exact t-based inference; in moderate-to-large samples the tests are relatively robust to mild deviations from normality.
• Variance structure: Welch’s test accommodates unequal variances; the pooled-variance test assumes equal variances and can be misleading when this assumption is violated.
• Multiple testing: when multiple parameters are tested, Sigma does not automatically apply multiplicity corrections (e.g., Bonferroni/FDR).

Warnings, edge cases, and safeguards

• Tests are only reported when both groups contain at least two finite observations and the corresponding standard error is finite and greater than zero.
• Constant or numerically non-estimable data are explicitly reported as non-estimable; no p-value or confidence interval is fabricated.
• The sign of \(\Delta\), \(d\), and \(g\) depends on the Control/Treatment labeling.

Numerical implementation and equivalence to R

Computations are performed in JavaScript using double-precision floating point arithmetic (IEEE-754). Student t tail probabilities are evaluated through the regularized incomplete beta function. This provides direct, stable two-sided tail probabilities without subtracting a CDF from one. Critical values are obtained by root-finding against the same distribution implementation. Small numerical differences across software may occur due to floating-point rounding.

Conceptually, the procedures correspond to R’s stats::t.test for independent samples: Welch’s test (var.equal = FALSE) and the pooled-variance test (var.equal = TRUE).

References

1. Student (Gosset). (1908). The probable error of a mean. Biometrika.
2. Welch, B. L. (1947). The generalization of “Student’s” problem when several different population variances are involved. Biometrika.
3. Satterthwaite, F. E. (1946). An approximate distribution of estimates of variance components. Biometrics Bulletin.
4. Cohen, J. (1988). Statistical Power Analysis for the Behavioral Sciences (2nd ed.).
5. Hedges, L. V., & Olkin, I. (1985). Statistical Methods for Meta-Analysis.

Paired t-Test

Purpose and scope

The Sigma paired t-test module performs a classical paired-samples (dependent-samples) t-test, mathematically equivalent to common implementations in R (t.test(..., paired = TRUE)) and standard statistical software. The test compares two measurements obtained from the same subjects (Time 1 and Time 2) and bases inference on the within-subject differences.

Input data, pairing modes, and preprocessing

The analysis is defined by one or more pairs of numeric columns representing measurements at Time 1 and Time 2. Sigma supports two pairing modes:

• Manual pairing: the user explicitly selects the Time 1 column and Time 2 column for each pair (and assigns a label for reporting).
• Automatic pairing: Sigma can automatically detect pairs by matching column names using a user-specified suffix pattern (e.g., _pre / _post) or prefix pattern (e.g., pre_ / post_) for Time 1 and Time 2. Detected pairs can be reviewed and converted into the manual list for reproducible reporting.

For each pair, rows are included if and only if both Time 1 and Time 2 values are present and can be parsed as finite numeric values and their difference is finite. Rows with a missing, non-numeric, or non-finite value in either column are excluded for that pair (pairwise complete-case evaluation). The CSV reader supports quoted fields, escaped quotes, comma, semicolon, and tab delimiters, and selectable decimal conventions. Numeric conversion is strict: a field such as 12abc is not accepted as the number 12.

Numeric column acceptance (70% rule)

To reduce accidental selection of text columns, Sigma uses a numeric eligibility heuristic. A column is offered for analysis if at least 70% of all non-missing entries can be parsed as finite numbers. Columns failing this rule are not analyzed. This classification rule is separate from the row-level complete-case evaluation: within an eligible column, every individual invalid value is still excluded from its pair.

The 70% threshold is a software heuristic, not a statistical assumption. Users should verify that the selected columns represent the intended continuous measurements and review the reported inclusion counts.

Definition of within-subject differences

Let \(x_{1,i}\) denote the measurement at Time 1 and \(x_{2,i}\) the measurement at Time 2 for subject \(i\), for all subjects with valid paired observations \(i=1,\ldots,n\). Sigma defines the within-subject difference as:

\[ d_i = x_{2,i} - x_{1,i}. \]

All subsequent statistics are computed on the vector of differences \(d=(d_1,\ldots,d_n)\).

Descriptive statistics

For Time 1, Time 2, and the differences, Sigma reports:

• sample size \(n\),
• mean \(\bar{x}\) (or \(\bar{d}\)),
• standard deviation \(s=\sqrt{s^2}\),
• standard error \(SE=s/\sqrt{n}\).

Means and variances are accumulated with a numerically stable online algorithm rather than from the potentially unstable identity \(\sum x_i^2-n\bar{x}^2\).

Paired t-test

Let \(\bar{d}\) be the mean of the differences and \(s_d\) their sample standard deviation:

\[ \bar{d}=\frac{1}{n}\sum_{i=1}^{n} d_i,\qquad s_d^2=\frac{1}{n-1}\sum_{i=1}^{n}(d_i-\bar{d})^2,\qquad SE_d=\frac{s_d}{\sqrt{n}}. \]

The paired t-statistic is:

\[ t=\frac{\bar{d}}{SE_d}=\frac{\bar{d}}{s_d/\sqrt{n}}, \qquad df=n-1. \]

The two-sided null hypothesis is \(H_0:\mu_d=0\), where \(\mu_d\) denotes the population mean of differences. Sigma computes the two-sided p-value from the Student t distribution. For numerical stability in the tails, it evaluates the equivalent regularized incomplete-beta expression directly:

\[ p=I_{\frac{df}{df+t^2}}\left(\frac{df}{2},\frac{1}{2}\right). \]

This avoids subtracting a CDF rounded close to one and therefore preserves very small non-zero p-values.

Confidence interval for the mean difference

Sigma reports a two-sided 95% confidence interval for \(\mu_d\) using \(\alpha=0.05\) and the 97.5% t-quantile \(t_{crit}=t_{1-\alpha/2,df}\):

\[ CI_{95\%}=\bar{d}\pm t_{crit}\cdot SE_d. \]

Effect sizes

Sigma reports standardized effect sizes based on the variability of within-subject differences. Cohen’s \(d_z\) for paired samples is:

\[ d_z=\frac{\bar{d}}{s_d}. \]

The corresponding small-sample corrected effect size is labelled Hedges’ \(g_z\). Sigma uses the exact gamma-function correction for \(df>1\):

\[ J(df)=\frac{\Gamma(df/2)}{\sqrt{df/2}\,\Gamma((df-1)/2)}, \qquad g_z=J(df)\,d_z. \]

For two valid pairs (\(df=1\)), the paired t-test remains defined but this exact correction is not; consequently, Sigma leaves \(g_z\) blank and reports a warning. Effect-size definitions for paired data are not interchangeable: \(d_z\) uses the standard deviation of the difference scores and should be named explicitly when reported.

Assumptions and limitations

• Independence at the subject level: pairs are assumed to come from independent subjects; the test does not model clustering beyond the pairing.
• Approximate normality of differences: exact t-based inference assumes the differences \(d_i\) are approximately normally distributed in the population. With increasing sample size, the test is generally robust to moderate deviations.
• Measurement scale: outcomes should be continuous (or approximately continuous) and measured consistently at both time points.
• Correct pairing and influential differences: the two cells in each row must belong to the same subject. Strong outliers in the difference scores can materially affect the mean, standard deviation, test statistic, and confidence interval.
• Complete-case interpretation: inference describes the observations with two usable measurements. Systematic missingness can bias the result and is not corrected by the paired t-test.

Warnings, edge cases, and safeguards

Sigma reports paired t-test results only when at least two valid pairs are available and the differences are not constant:

• \(n \ge 2\) valid paired observations are required.
• If the variance of differences is zero (i.e., \(s_d=0\)), no t-statistic is reported.
• Numerically almost-constant differences are rejected because reliable t-based inference cannot be formed.
• Pairs are rejected if one or both columns fail the numeric acceptance rule (70% rule).
• Non-finite differences, including subtraction overflow, are excluded and counted separately.

When several column pairs are analyzed, Sigma reports a multiple-testing warning. It does not automatically adjust p-values; any multiplicity strategy must be selected according to the study design. The module performs a two-sided classical paired t-test and does not substitute a non-parametric test when assumptions are doubtful.

Numerical implementation and equivalence to R

Computations are performed in JavaScript using double-precision floating point arithmetic (IEEE-754). Student-t probabilities are obtained from the regularized incomplete beta function; the 97.5% quantile is found by bracketed numerical inversion of the t CDF. The almost-constant-data check follows the numerical safeguard used by R. The procedure corresponds to R’s stats::t.test with paired = TRUE, alternative = "two.sided", mu = 0, and conf.level = 0.95. Minor residual differences can occur from floating-point evaluation.

References

1. Student (Gosset). (1908). The probable error of a mean. Biometrika.
2. Cohen, J. (1988). Statistical Power Analysis for the Behavioral Sciences (2nd ed.).
3. Hedges, L. V., & Olkin, I. (1985). Statistical Methods for Meta-Analysis.
4. Lakens, D. (2013). Calculating and reporting effect sizes to facilitate cumulative science: a practical primer for t-tests and ANOVAs. Frontiers in Psychology, 4, 863.
5. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.

One-Way ANOVA — Welch, Classical, Tukey, and Games–Howell

Purpose and scope

The One-Way ANOVA module compares the mean of a continuous outcome across two or more mutually exclusive, independent groups defined by one categorical factor. Several continuous outcomes may be selected in one run, but each outcome is analyzed separately. Welch ANOVA is the recommended default; classical equal-variance ANOVA is available as the primary or alternative omnibus method. Sigma also reports group descriptives, classical effect sizes, the median-centered Brown–Forsythe variance diagnostic, and multiplicity-adjusted pairwise comparisons.

The module is intended for one-way between-subject designs. It does not fit repeated-measures, factorial, mixed, multilevel, ANCOVA, or regression-adjusted models, and it does not create graphics.

Input structure, numeric detection, and exclusions

Input is a CSV file in long format, with one independent observation per row. The user selects one group column and one or more continuous outcome columns. Every distinct non-missing group value defines a level, retained in order of first appearance. An outcome is eligible for selection when at least 70% of its non-missing cells can be parsed as finite numbers under the selected decimal convention.

For a given outcome, rows with a missing group, missing outcome, or non-numeric/non-finite outcome are excluded and counted separately. No imputation is performed. A group level with no usable outcome is dropped and reported. At least two populated groups and at least two usable observations per group are required for the standard inferential workflow. Because outcomes are processed independently, their complete-case cohorts may differ.

Group descriptives

For each group Sigma reports the usable sample size \(n_j\), mean, sample standard deviation, standard error, median, and interquartile range. The mean confidence interval is

\[ \bar x_j \pm t_{0.975,\,n_j-1}\frac{s_j}{\sqrt{n_j}}. \]

These are separate pointwise 95% intervals for the group means. They are not simultaneous intervals and should not be substituted for the adjusted pairwise intervals described below.

Classical one-way ANOVA

Let \(k\) be the number of groups, \(N=\sum_j n_j\), \(\bar x_j\) the mean in group \(j\), and \(\bar x\) the overall mean. The classical decomposition is

\[ SS_B=\sum_{j=1}^{k}n_j(\bar x_j-\bar x)^2, \qquad SS_W=\sum_{j=1}^{k}\sum_{i=1}^{n_j}(x_{ij}-\bar x_j)^2, \qquad SS_T=SS_B+SS_W. \]

With \(MS_B=SS_B/(k-1)\) and \(MS_W=SS_W/(N-k)\), the omnibus statistic is

\[ F=\frac{MS_B}{MS_W}, \qquad df_1=k-1, \qquad df_2=N-k. \]

Under the null hypothesis all population group means are equal. Classical F inference assumes independent observations, normally distributed within-group errors, and a common population variance. Moderate non-normality is often tolerable in sufficiently large, reasonably balanced samples, but strong heteroskedasticity combined with unequal group sizes can distort the test.

Welch one-way ANOVA

Welch ANOVA does not pool the group variances. Define \(w_j=n_j/s_j^2\), \(W=\sum_j w_j\), the weighted mean \(\bar x_w=\sum_jw_j\bar x_j/W\), and

\[ B=\sum_{j=1}^{k}\frac{(1-w_j/W)^2}{n_j-1}. \]

Sigma calculates

\[ F_W= \frac{\displaystyle \frac{1}{k-1}\sum_{j=1}^{k}w_j(\bar x_j-\bar x_w)^2} {\displaystyle 1+\frac{2(k-2)}{k^2-1}B}, \qquad df_1=k-1, \qquad df_2=\frac{k^2-1}{3B}. \]

Welch ANOVA permits unequal variances and unequal group sizes and is therefore the default primary method. It still assumes independent observations and uses an approximate F distribution. Every group must have a finite, strictly positive sample variance. Selecting Welch by default is not a license to ignore data quality, extreme outliers, dependence, or inappropriate outcome scales.

Brown–Forsythe variance diagnostic

For each observation Sigma forms the absolute deviation from its group median, \(z_{ij}=|x_{ij}-\tilde x_j|\), and applies the classical one-way ANOVA calculation to these deviations. The resulting Brown–Forsythe F test evaluates equality of group variances using a median-centered form of Levene's test. Median centering is less sensitive to non-normality than mean-centered Levene testing.

The diagnostic is reported separately and never switches the primary inference method automatically. Choosing classical or Welch inference by first testing variances and then selecting a method can alter the overall type-I error behavior. Method choice should instead reflect the design, plausible variance structure, balance, and prespecified analysis plan. A non-significant Brown–Forsythe result is not proof that variances are equal, particularly in small samples.

Classical effect-size estimates

Sigma reports eta squared and omega squared from the classical sums-of-squares decomposition:

\[ \eta^2=\frac{SS_B}{SS_T}, \qquad \omega^2=\frac{SS_B-(k-1)MS_W}{SS_T+MS_W}. \]

Eta squared is the observed proportion of total variation assigned to between-group differences and is upward biased as a population effect-size estimator. Omega squared applies a finite-sample correction and can be negative; Sigma does not silently truncate it to zero. Both quantities are tied to the classical model decomposition. They are descriptive effect-size estimates, not robust Welch effect sizes, and their interpretation depends on the sampled groups, outcome scale, and study design.

Tukey HSD and the Tukey–Kramer extension

Tukey post-hoc inference uses the pooled \(MS_W\) and residual degrees of freedom from classical ANOVA. For groups \(a\) and \(b\), Sigma reports the contrast in the explicit direction \(\Delta_{ba}=\bar x_b-\bar x_a\), its ordinary difference standard error

\[ SE(\Delta_{ba})=\sqrt{MS_W\left(\frac{1}{n_a}+\frac{1}{n_b}\right)}, \]

and the studentized-range statistic and simultaneous interval

\[ q=\frac{|\Delta_{ba}|}{SE(\Delta_{ba})/\sqrt{2}}, \qquad \Delta_{ba}\pm q_{0.95;\,k,\,N-k}\frac{SE(\Delta_{ba})}{\sqrt{2}}. \]

With unequal sample sizes this is the Tukey–Kramer extension, although the interface uses the familiar label “Tukey HSD.” Adjusted p-values and simultaneous 95% intervals control the family-wise error rate for all pairwise comparisons within that outcome under the classical equal-variance model.

Games–Howell pairwise comparisons

Games–Howell is the variance-heterogeneous companion to Welch ANOVA. For each pair define

\[ SE_{ab}=\sqrt{\frac{s_a^2}{n_a}+\frac{s_b^2}{n_b}}, \qquad \nu_{ab}= \frac{\left(s_a^2/n_a+s_b^2/n_b\right)^2} {\frac{(s_a^2/n_a)^2}{n_a-1}+\frac{(s_b^2/n_b)^2}{n_b-1}}. \]

The reported statistic and simultaneous interval are

\[ q_{ab}=\frac{\sqrt{2}|\Delta_{ba}|}{SE_{ab}}, \qquad \Delta_{ba}\pm q_{0.95;\,k,\,\nu_{ab}}\frac{SE_{ab}}{\sqrt{2}}. \]

Each pair has its own Welch–Satterthwaite degrees of freedom. The studentized-range distribution uses the total number of analyzed groups \(k\), not merely the two groups in the pair. Games–Howell does not assume a common variance, but very small groups and unstable variance estimates can still produce wide intervals and low effective degrees of freedom.

Post-hoc selection, omnibus tests, and multiplicity

In automatic mode Sigma pairs Games–Howell with Welch ANOVA and Tukey with classical ANOVA. The user may request either method, both, or no post-hoc analysis. Pairwise results are calculated even when the corresponding omnibus p-value is not below 0.05; this permits prespecified comparisons and avoids imposing an automatic omnibus gate. Researchers should nevertheless distinguish planned contrasts from exploratory post-hoc testing and interpret the complete family of results.

Family-wise adjustment applies separately within each outcome and selected post-hoc method. If several outcomes are analyzed, Sigma does not adjust omnibus or post-hoc results across outcomes and displays a warning. Any broader multiplicity strategy must be defined from the scientific testing hierarchy.

Numerical implementation and safeguards

Calculations use IEEE-754 double precision. A common translation is removed internally before means and sums of squares are accumulated, preserving the statistics while protecting analyses whose observations share a very large numerical offset. F probabilities use the regularized incomplete beta function. Studentized-range probabilities are evaluated by numerical integration over the normal-range and chi-square mixing distributions; quantiles are obtained by bracketed inversion. For two groups Sigma uses the exact identity between the studentized range and the two-sided Student t distribution.

Classical ANOVA is reported as infinite only when between-group variation is positive and pooled within-group variance is exactly zero; Tukey intervals are then not estimable. Welch and Games–Howell require a positive variance in every group. Constant outcomes, insufficient group sizes, missing residual degrees of freedom, and non-finite intermediate values produce controlled diagnostic messages instead of fabricated estimates.

Assumptions and limitations

• Independence: observations must be independent within and between groups. Clustered, paired, or repeated observations require another model.
• Outcome scale: the mean and variance must be scientifically meaningful for the selected outcome.
• Distribution: omnibus and pairwise inference rely on normal-theory approximations; small samples, strong skewness, and influential outliers require particular caution.
• Variance structure: classical ANOVA and Tukey assume equal population variances; Welch and Games–Howell relax this assumption but do not solve dependence or severe model misspecification.
• One factor only: interactions, confounder adjustment, factorial effects, random effects, and repeated measures are outside the module.
• No automatic scientific decision: Brown–Forsythe, p-values, and warning flags support judgment but do not select a scientifically appropriate model automatically.
• No graphics: distributional and residual graphics should be created and assessed separately, for example in Phi Statistics Drawer.

References

1. Fisher RA. Statistical Methods for Research Workers. Oliver and Boyd; 1925.
2. Welch BL. On the comparison of several mean values: an alternative approach. Biometrika. 1951;38:330–336.
3. Brown MB, Forsythe AB. Robust tests for the equality of variances. Journal of the American Statistical Association. 1974;69:364–367.
4. Tukey JW. Comparing individual means in the analysis of variance. Biometrics. 1949;5:99–114.
5. Kramer CY. Extension of multiple range tests to group means with unequal numbers of replications. Biometrics. 1956;12:307–310.
6. Games PA, Howell JF. Pairwise multiple comparison procedures with unequal n's and/or variances: a Monte Carlo study. Journal of Educational Statistics. 1976;1:113–125.
7. Lakens D. Calculating and reporting effect sizes to facilitate cumulative science: a practical primer for t-tests and ANOVAs. Frontiers in Psychology. 2013;4:863.

Univariable Logistic Regression

Purpose and scope

The Sigma univariable logistic regression module fits classical binary logistic regression models with an intercept and one predictor at a time. For each selected predictor, Sigma fits a separate univariable model using maximum likelihood estimation (MLE) via an iteratively reweighted least squares (IRLS) algorithm (Newton–Raphson for the logistic log-likelihood). The approach is mathematically equivalent to standard implementations such as R’s glm(..., family = binomial) for single-predictor models.

Data structure, outcome definition, and predictors

The tool expects a binary outcome variable and one or more candidate predictor variables. The user selects:

• a binary outcome column,
• the outcome level treated as the event (coded as 1),
• the outcome level treated as the non-event (coded as 0),
• one or more numeric predictors to be analyzed in separate univariable models.

Internally, the outcome is recoded to \(y_i \in \{0,1\}\), where \(y_i=1\) denotes the user-defined event and \(y_i=0\) denotes the non-event. Each predictor \(x_i\) is treated as continuous (no automatic categorization or spline/nonlinear terms).

The analysis can be started only when the outcome contains exactly two non-missing levels. This prevents accidental omission of a third outcome category.

Missing data handling and numeric checks

For each predictor, Sigma performs a complete-case analysis on the pair \((y_i, x_i)\):

• rows with missing outcome values are excluded,
• rows with missing or non-numeric predictor values are excluded,
• all remaining rows with valid \((y_i,x_i)\) are included for that predictor’s model.

As a data-quality safeguard, a column is offered as numeric only if at least 70% of up to the first 30 non-empty entries are valid finite numbers. Parsing honors the selected decimal separator and requires the complete cell to be a valid numeric representation; malformed strings such as 12abc are not accepted as partial numbers.

For each predictor, Sigma reports:

• \(n_{total}\): number of rows with valid outcome and predictor,
• \(n_{events}\): number of rows with \(y=1\),
• \(n_{non\text{-}events}\): number of rows with \(y=0\),
• event rate (%): \(100 \cdot n_{events}/n_{total}\).

Logistic regression model

For each predictor, Sigma fits the binary logistic regression model:

\[ Y_i \sim \text{Bernoulli}(p_i), \qquad \text{logit}(p_i)=\log\!\left(\frac{p_i}{1-p_i}\right)=\beta_0+\beta_1 x_i, \]

with fitted probability \(p_i = \frac{1}{1+\exp\!\left(-(\beta_0+\beta_1 x_i)\right)}\). Here, \(\beta_1\) is the log-odds change per unit increase in \(x\).

The log-likelihood for a sample of size \(n\) is:

\[ \ell(\beta_0,\beta_1)=\sum_{i=1}^{n}\left[y_i\log(p_i)+(1-y_i)\log(1-p_i)\right]. \]

The MLEs \((\hat\beta_0,\hat\beta_1)\) maximize \(\ell\) and solve the corresponding score equations.

Estimation via IRLS (Newton–Raphson)

To make the numerical calculation invariant to the measurement unit, Sigma centers and scales the predictor internally. Coefficients and their covariance matrix are transformed back to the original unit before display. The reported odds ratio therefore always refers to a one-unit increase in the original predictor.

Sigma initializes the intercept from the observed event odds and the slope at zero. It then uses IRLS with a maximum of 100 iterations and a relative parameter-update tolerance of 10⁻¹⁰, together with a scaled-score convergence check. At iteration \(t\), define:

\[ \eta_i=\beta_0^{(t)}+\beta_1^{(t)}x_i,\qquad p_i=\frac{1}{1+\exp(-\eta_i)},\qquad w_i=p_i(1-p_i). \]

Probabilities and the Bernoulli log-likelihood are evaluated with sign-aware logistic and log1p-based formulas. The linear predictor is not artificially clamped.

The score components for \(\beta=(\beta_0,\beta_1)\) are:

\[ U_0=\sum_{i}(y_i-p_i),\qquad U_1=\sum_{i}(y_i-p_i)x_i. \]

The observed Fisher information matrix for the two-parameter model is:

\[ I(\beta)= \begin{pmatrix} I_{00} & I_{01}\\ I_{01} & I_{11} \end{pmatrix}, \quad I_{00}=\sum_i w_i,\quad I_{01}=\sum_i w_i x_i,\quad I_{11}=\sum_i w_i x_i^2. \]

The Newton update solves \(I(\beta)\Delta\beta = U(\beta)\) by Cholesky factorization. A likelihood-monitored step-halving procedure shortens a Newton step if the full step would reduce the log-likelihood. Iteration stops only when both the relative update and scaled score meet their tolerances. If convergence is not achieved, Sigma suppresses the coefficient, odds ratio, confidence interval, and p-value rather than reporting an interim iterate.

Standard errors, Wald test, odds ratio, and confidence interval

After convergence, Sigma recomputes the Fisher information matrix at the final maximum-likelihood estimate and obtains the covariance matrix by Cholesky-based inversion. For the slope \(\hat\beta_1\):

\[ \mathrm{Var}(\hat\beta_1)=\left[I(\hat\beta)^{-1}\right]_{11}, \qquad SE(\hat\beta_1)=\sqrt{\mathrm{Var}(\hat\beta_1)}. \]

A Wald z-statistic is computed as \(z=\hat\beta_1/SE(\hat\beta_1)\), and the two-sided p-value is:

\[ p = 2\left(1-\Phi(|z|)\right), \]

where \(\Phi\) is the standard normal CDF. The two-sided tail probability is evaluated directly through a high-accuracy complementary-error-function calculation, avoiding cancellation from subtracting a CDF value close to one. Sigma uses the critical value \(z_{0.975}=1.9599639845\) for 95% Wald intervals. The odds ratio and its 95% confidence interval are:

\[ OR=\exp(\hat\beta_1),\qquad CI_{95\%}(OR)=\left[\exp(\hat\beta_1-z_{0.975}SE),\ \exp(\hat\beta_1+z_{0.975}SE)\right]. \]

Warnings and quality safeguards

To improve interpretability in pathological data constellations, Sigma attaches warning codes to each fitted univariable model and displays these in the results table and detailed view.

• Warning 1 — Non-estimable coefficient: Issued for a single-level outcome, a constant predictor, or complete/quasi-complete separation. Separation is checked for numeric predictors generally, not only binary predictors, by comparing the observed predictor ranges of events and non-events. In these situations the usual finite MLE is absent or the slope is not identifiable. Sigma suppresses the coefficient, standard error, odds ratio, confidence interval, and p-value.
• Warning 2 — Very few events or non-events: Issued when \(n_{events} < 5\) or \(n_{non\text{-}events} < 5\) in the complete-case subset for the predictor. Estimates are computed, but inference may be unstable and confidence intervals wide.
• Warning 3 — Numerical instability: Indicates strong ill-conditioning, extreme fitted coefficients, or lack of strict IRLS convergence. Results from a converged but flagged model require caution; if the model does not converge, inferential quantities are suppressed.

Models are only fitted when both outcome categories are present (at least one event and one non-event) and the predictor shows non-zero variability in the complete-case subset. Otherwise, the model is not fitted and an appropriate warning is issued.

Assumptions and limitations

• Independence: observations are assumed independent at the subject level.
• Linearity on the logit scale: the model assumes a linear relationship between the predictor and the log-odds of the event; nonlinearity is not modeled (no splines/polynomials).
• Large-sample inference: p-values and confidence intervals rely on Wald (normal) approximations; with small samples or rare events, Wald inference can be unstable.
• No penalized/bias-reduced methods: the current implementation does not include Firth correction or likelihood-ratio-based confidence intervals.

Numerical implementation and equivalence to R

Calculations use IEEE-754 double-precision arithmetic. Centering and scaling are numerical devices only and are reversed analytically for all reported quantities. Positive-definite systems are solved by Cholesky factorization with relative, scale-aware pivot checks. Two-sided normal probabilities are computed from the complementary error function.

Conceptually, Sigma corresponds to R’s stats::glm with family = binomial (logit link). Small numerical differences across software may occur due to rounding, convergence criteria, and implementation details, especially in separation scenarios where Sigma intentionally suppresses non-identifiable estimates.

References

1. Cox, D. R. (1958). The regression analysis of binary sequences. Journal of the Royal Statistical Society. Series B.
2. McCullagh, P., & Nelder, J. A. (1989). Generalized Linear Models (2nd ed.). Chapman & Hall/CRC.
3. Agresti, A. Categorical Data Analysis. Wiley (current edition).
4. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.

Multivariable Logistic Regression

Overview

The Sigma multivariable logistic regression module implements a classical binary logistic regression model with one intercept and multiple predictors, mathematically equivalent to glm(..., family = binomial) in R and standard implementations in SPSS or Stata for multivariable models. Parameters are estimated by maximum likelihood using an iteratively reweighted least squares (IRLS) / Newton–Raphson algorithm in IEEE-754 double precision.

In contrast to the univariable module, all selected predictors enter the same model simultaneously. Predictors with zero variance or exact linear dependence in the complete-case sample are detected, reported, and excluded from the fitted model.

Data Structure, Outcome Definition and Predictors

The tool requires a binary outcome and at least one numeric predictor. The user specifies:

• a binary outcome column (e.g. 0 = non-event, 1 = event),
• which outcome value is treated as the “event” category,
• which outcome value is treated as the “non-event” category,
• one or more numeric predictor columns to be included together in a single multivariable model.

Internally, the outcome is recoded to yi ∈ {0, 1}, with yi = 1 representing the user-defined event and yi = 0 the non-event. Predictors are treated as continuous without automatic categorization, transformation, or interaction terms.

The analysis can be started only when the outcome contains exactly two non-missing levels. This prevents a third outcome category from being silently omitted.

Missing Data, Numeric Checks and Listwise Deletion

Sigma applies listwise deletion for the multivariable model: a row is included only if the outcome and all selected predictors are present and numeric. Complete cases are formed before variance and rank checks, so excluded rows cannot influence whether a predictor is treated as constant or collinear.

Numeric columns are identified conservatively: a column is considered numeric only if at least 70% of up to the first 30 non-empty values are valid finite numbers (using the selected decimal separator setting). The complete cell must be numeric; partial parses such as 12abc are rejected. Non-numeric columns are not offered as predictors.

After listwise deletion, let n be the number of complete cases used in the model and k the number of predictors actually fitted (excluding the intercept, constants, and exactly aliased predictors). The module reports n, number of events nevents, non-events nnon-events, and the event rate:

Event rate (%) = 100 · nevents / n.

The events-per-variable ratio (EPV) is computed as:

EPV = nevents / k.

Handling of Constant and Exactly Collinear Predictors

Predictors with zero variance (all values identical) in the analysis sample cannot be estimated separately from the intercept. Such predictors are flagged, excluded from the regression fit, and listed as “dropped predictors” in the model summary. In the results table, these predictors are shown with counts and a warning code indicating that no meaningful coefficient / odds ratio is available. Sigma also applies a rank-revealing, two-pass modified Gram–Schmidt check to the standardized design matrix. If a predictor is an exact linear combination of predictors already retained (relative residual threshold 10⁻¹⁰), its individual coefficient is aliased and is left blank; the remaining identifiable model is still fitted. Predictor order determines which member of an exactly dependent set is retained.

Multivariable Logistic Regression Model

For subject i, let Yi be the binary outcome and Xi· = (xi1, …, xik) the predictor values. The model is:

Yi ~ Bernoulli(pi)
logit(pi) = log( pi / (1 − pi) ) = β0 + Σj=1k βj · xij

where βj is the log-odds ratio per one-unit increase in predictor Xj, adjusted for all other predictors.

In matrix form, with design matrix X of dimension n × (k+1) (first column = intercept), parameter vector β, linear predictor η = Xβ, and probabilities pi = 1/(1+exp(−ηi)), the log-likelihood is:

ℓ(β) = Σi=1n { yi · log(pi) + (1 − yi) · log(1 − pi) }.

Estimation via IRLS (Newton–Raphson)

Each retained predictor is centered and scaled internally to improve numerical conditioning. After fitting, coefficients and the full covariance matrix are transformed back to the original predictor units. Sigma uses IRLS with a maximum of 100 iterations, a relative parameter-update tolerance of 10⁻¹⁰, and a scaled-score convergence check. At iteration t:

• η(t) = X β(t)
• pi(t) = 1/(1+exp(−ηi(t)))
• wi(t) = pi(t)(1−pi(t))

The score vector and observed Fisher information are computed as:

U(β(t)) = Xᵀ (y − p(t))
I(β(t)) = Xᵀ W(t) X

and the Newton step solves I Δβ = U by Cholesky factorization with scale-aware pivot checks. Parameters are updated via β(t+1) = β(t) + Δβ. A likelihood-monitored step-halving procedure shortens the update if a full Newton step would reduce the log-likelihood.

Probabilities and log-likelihood contributions use sign-aware logistic and log1p-based formulas; the linear predictor is not clamped. A model is declared converged only when both update and score criteria are met. If convergence fails, Sigma leaves all inferential coefficient fields blank.

Standard Errors, Wald Tests and Confidence Intervals

After convergence, the information matrix is recomputed at the final estimate and inverted by Cholesky factorization to obtain the covariance matrix V = I(β̂)−1. Standard errors are:

SE(β̂j) = √Vjj.

Wald statistics and two-sided p-values are computed as:

zj = β̂j / SE(β̂j)
pj = 2 · (1 − Φ(|zj|))

The two-sided normal tail is evaluated directly through a high-accuracy complementary-error-function calculation. This avoids cancellation from subtracting a CDF value close to one.

Odds ratios and 95% Wald confidence intervals (for predictors; intercept excluded) are:

ORj = exp(β̂j)
CIlower(ORj) = exp(β̂j − z0.975 · SE(β̂j))
CIupper(ORj) = exp(β̂j + z0.975 · SE(β̂j))

with z0.975 ≈ 1.96. The intercept is used for model fitting and diagnostics but its OR is not displayed.

Model-Level Diagnostics and Warnings

The model summary reports outcome labels, n, events/non-events, event rate, predictors fitted/selected, EPV, convergence status, constant predictors, and exactly aliased predictors.

Model-level warnings include:

• Warning A: very few events or non-events (e.g. <5) → estimates may be unstable.
• Warning B: low EPV (EPV < 5 high risk; 5 ≤ EPV < 10 caution).
• Warning C: lack of convergence within the iteration limit → inferential coefficient fields are suppressed.
• Warning D: indications of separation or strong multicollinearity (e.g. near-singular information matrix, failure to solve/invert, extreme coefficient/OR/CI behavior) → OR/CI may be inflated or undefined.

Predictor-Level Warning Codes

• Warning 1: no meaningful estimate (single-level outcome, constant or exactly aliased predictor, or complete/quasi-complete separation) → β/SE/OR/CI/p suppressed.
• Warning 2: very few events or non-events (<5) → unstable estimates, wide CIs likely.
• Warning 3: strong or quasi-separation / numerical instability → estimates may be inflated; interpret with caution.

Assumptions and Limitations

• Independent observations at the subject level.
• Correct logit link specification; no automatic modeling of non-linearity (no splines, fractional polynomials) and no interactions.
• Wald inference relies on large-sample approximations; for small samples/rare events results may be unstable.
• No penalized/bias-reduced methods (e.g. Firth), no exact logistic regression, and no likelihood-ratio-based intervals.

Equivalence to R and numerical notes

For well-behaved datasets (adequate events, no extreme separation, no near-singular design), results are expected to agree closely with R’s glm(..., family = binomial) up to rounding differences. Differences may occur in edge cases because Sigma explicitly flags instability and leaves quantities undefined when no converged finite MLE or invertible final information matrix is obtained. Because Sigma evaluates Fisher information at the final MLE, tiny standard-error differences from implementations that retain weights from their last IRLS working step can occur.

Propensity Score Matching (PSM)

Purpose and scope

The Sigma Propensity Score Matching (PSM) module implements a streamlined propensity score matching workflow for observational comparisons of a binary exposure or treatment T (treated vs. control). Its goal is to reduce confounding by constructing a matched sample with improved balance in measured baseline covariates. The module is intended primarily as a design-stage tool: it estimates propensity scores, performs matching under the implemented fixed settings, and reports overlap and covariate balance diagnostics. Outcome analyses are conceptually performed after matching using methods appropriate for paired matched data.

Propensity score methods rely on the assumptions of (i) conditional exchangeability / no unmeasured confounding given the selected covariates, and (ii) overlap / positivity, meaning that subjects with similar covariate patterns have a non-zero probability of belonging to either group. These assumptions cannot be verified from observed data alone; the reported diagnostics evaluate only measured covariates and empirical overlap.

Notation and estimand

Let T ∈ {0,1} denote treatment assignment, with T=1 for treated and T=0 for controls, and let X denote the vector of measured baseline covariates. The propensity score is defined as:

\[ e(X) = \Pr(T = 1 \mid X). \]

Matching on e(X) aims to create treated and control groups that are comparable with respect to the observed covariates X. In Sigma, the implemented matching procedure is ATT-oriented: treated subjects are the index units, and control subjects are selected as nearest neighbors when suitable matches exist. Consequently, the matched sample is most naturally interpreted as targeting the average treatment effect on the treated (ATT) within the region of empirical overlap.

Inputs and variable types

• Treatment indicator: the user selects one column and specifies which level represents T=1 (“treated”) and which represents T=0 (“control”). Exactly two non-missing observed levels are required; otherwise preparation and matching are blocked.
• Covariates: the user selects baseline covariates to be included in the propensity score model. Numeric covariates enter the model as continuous predictors. Binary covariates enter as a single 0/1-coded predictor. Categorical covariates are dummy-coded using reference coding: for a factor with m observed levels, Sigma creates m−1 indicator variables, with the first deterministically sorted level used as the reference category. Any two-level covariate is automatically recognized as binary; its deterministically sorted first level is coded 0 and its second level 1. The final mapping is shown in the preprocessing audit and exports, and every automatically assigned covariate type can be overridden by the user.

Missing data handling (as implemented)

In the current Sigma PSM tool, missing data are handled by complete-case analysis only. Rows are excluded before propensity score estimation and matching if any of the following is true:

• the treatment value is missing,
• any selected covariate is missing, or
• a covariate specified as numeric cannot be parsed as a finite numeric value.

Accordingly, the matched analysis pertains to the subset of observations with complete information on treatment and all selected covariates. Complete-case analysis is transparent and reproducible, but it may reduce sample size and may induce bias if missingness is not completely at random or conditionally ignorable given the retained data.

Propensity score model

Sigma estimates propensity scores using a logistic regression model with intercept:

\[ T_i \sim \text{Bernoulli}(p_i), \qquad \text{logit}(p_i)=\log\!\left(\frac{p_i}{1-p_i}\right)=\beta_0 + X_i^\top\beta. \]

The fitted propensity score is:

\[ \hat e(X_i)=\hat p_i = \frac{1}{1+\exp\left(-(\hat\beta_0 + X_i^\top\hat\beta)\right)}. \]

Estimation (maximum likelihood via IRLS)

Parameters are estimated by maximum likelihood using an iteratively reweighted least squares (IRLS) / Newton–Raphson procedure in double precision. At iteration t, with current parameter vector \(\beta^{(t)}\),

\[ \eta_i^{(t)} = \beta_0^{(t)} + X_i^\top \beta^{(t)}, \quad p_i^{(t)} = \frac{1}{1+\exp(-\eta_i^{(t)})}, \quad w_i^{(t)} = p_i^{(t)}(1-p_i^{(t)}). \]

Using the classical working-response formulation,

\[ z_i^{(t)} = \eta_i^{(t)} + \frac{T_i - p_i^{(t)}}{w_i^{(t)}}, \]

the updated parameter vector is obtained by solving

\[ (X^\top W^{(t)} X)\,\beta^{(t+1)} = X^\top W^{(t)} z^{(t)}. \]

Before estimation, non-constant predictor columns are internally centered and scaled; coefficients are transformed back to the original units after fitting. This makes fitted probabilities invariant to a change of measurement units while leaving all reported balance values on the original scale. Weighted least-squares steps are solved with a column-pivoted, rank-aware QR decomposition rather than normal equations. Constant and aliased columns are identified and reported. Iteration uses the same outcome-based binomial initialization principle as R glm(..., family = binomial) and stops by relative deviance change or at 100 iterations. A small lower numerical bound is applied only to working weights to prevent division by zero. No ridge penalty is added and fitted propensity scores are not artificially clamped. Non-convergence blocks matching; extreme fitted probabilities and rank deficiency are reported as diagnostics because they may indicate separation or an unstable model specification.

Distance measure for matching

In the current implementation, matching is performed on the propensity score scale itself. For subject i, the matching score is

\[ s_i = \hat p_i. \]

The distance between a treated subject i and a control subject j is the absolute difference in fitted propensity scores:

\[ d_{ij} = \left| \hat p_i - \hat p_j \right|. \]

Matching algorithm (nearest neighbor; 1:1; without replacement)

Sigma implements greedy nearest-neighbor matching with a fixed 1:1 ratio and without replacement. Each treated subject is considered as a candidate for matching to at most one control, and each control subject can be used at most once.

Treated subjects are processed in descending order of the chosen matching score (largest propensity score first), with deterministic tie-breaking by row order. For each treated subject, Sigma searches the set of currently available controls and selects the control with the smallest absolute distance \(d_{ij}\), again using deterministic tie-breaking. If no eligible control satisfies the caliper restriction, the treated subject remains unmatched.

Because matching is greedy and without replacement, the final matched sample depends on the sequential allocation of controls. This is methodologically standard for nearest-neighbor matching and is consistent with an ATT-oriented design, but it is not equivalent to optimal global matching.

Caliper matching

Sigma supports caliper-restricted matching. A candidate treated–control pair is accepted only if

\[ d_{ij} \le c, \]

where c is the caliper threshold. In the current implementation, the user may either specify an absolute caliper directly on the propensity score scale or use a standardized caliper. For the standardized option, Sigma interprets the entered value \(\gamma\) as a multiple of the empirical standard deviation of the matching distance scale in the full analysis sample. Since the implemented distance scale is the propensity score itself, this becomes

\[ c = \gamma \cdot \mathrm{SD}(\hat p). \]

Thus, with the default standardized setting and caliper value 0.2, Sigma uses \(c = 0.2 \times \mathrm{SD}(\hat p)\) in the complete-case analysis sample.

A smaller caliper improves match closeness but may leave more treated subjects unmatched. Therefore, the practical estimand remains ATT-like but increasingly restricted to the subset of treated subjects lying within empirical overlap.

Overlap, common support, and unmatched units

Sigma reports overlap diagnostics on the propensity score scale. These include the pre- and post-matching propensity score ranges in treated and control subjects, their overlap interval, selected quantiles, and the proportion of units lying outside the range of the opposite group. These diagnostics are descriptive and warning-oriented; they do not alter the matching process except through the user-specified caliper. A zero-width intersection is explicitly distinguished from non-overlap: it is reported as support limited to a single point.

Unmatched treated units indicate limited overlap and imply that the matched analysis pertains to a subset of treated subjects for whom suitable controls were available under the implemented matching constraints.

Balance diagnostics (pre- and post-matching)

Sigma reports covariate balance before and after matching. The primary diagnostic is the standardized mean difference (SMD), computed on the design matrix actually used for propensity score estimation, i.e., numeric predictors, binary predictors, and dummy-coded categorical predictors. For categorical covariates, Sigma additionally reports the reference-level contrast reconstructed from the dummy coding so that all observed levels are represented in the balance output.

SMD definition used in Sigma

Sigma uses an ATT-style standardization. For each balance feature Z, let \(\bar Z_T\) and \(\bar Z_C\) denote the treated and control means (or proportions for binary indicators), and let \(s_{T,\mathrm{pre}}\) denote the standard deviation of Z in the treated group before matching. Then Sigma reports

\[ \mathrm{SMD} = \frac{\bar Z_T - \bar Z_C}{s_{T,\mathrm{pre}}}. \]

The denominator \(s_{T,\mathrm{pre}}\) is held constant for the pre- and post-matching SMD of the same feature. This yields an ATT-oriented balance metric and makes post-matching changes directly interpretable relative to the original variability in the treated group.

Pre- and post-matching samples

Pre-matching balance is computed in the complete-case analysis sample used for propensity score estimation. Post-matching balance is computed in the matched sample. Because the current implementation uses 1:1 matching without replacement, each matched treated and each matched control contributes weight 1 in the post-matching balance calculations.

Additional descriptive balance metric

Sigma also reports z-diff as an additional descriptive measure. With group-specific weighted sample variances \(s_T^2\) and \(s_C^2\), and Kish effective sample sizes \(n_{\mathrm{eff},T}\) and \(n_{\mathrm{eff},C}\), it is calculated as

\[ z_{\mathrm{diff}} = \frac{\bar Z_T-\bar Z_C} {\sqrt{s_T^2/n_{\mathrm{eff},T}+s_C^2/n_{\mathrm{eff},C}}}, \qquad n_{\mathrm{eff}}=\frac{(\sum_i w_i)^2}{\sum_i w_i^2}. \]

In the current 1:1 design without replacement all retained weights are 1, so the effective sample sizes equal the matched group counts. The statistic is not used as a hypothesis test: after matching, paired observations are dependent, whereas this descriptive standard error does not model within-pair correlation. SMD remains the primary balance diagnostic.

Interpreting balance

There is no universal cutoff, but in many applied settings absolute SMD values below approximately 0.1 are considered compatible with good covariate balance. Balance should be assessed across all included covariates rather than relying on a single summary value alone.

Outputs and exports

Sigma reports:

• the number of rows retained after complete-case preprocessing, exclusion reasons, and binary mappings,
• propensity score model diagnostics (convergence, iterations, deviance, estimable rank, and score range),
• matching diagnostics (matched and unmatched treated/control counts, distance scale, caliper, effective caliper),
• overlap / common-support diagnostics on the propensity score scale,
• covariate balance tables with pre- and post-matching SMD and z-diff values.

The module can additionally export the matched cohort and the match pairs, including metadata for downstream analyses. For the current 1:1 without-replacement design, Sigma also supports paired wide-format export.

After matching: outcome analysis (principles)

Because the current implementation produces 1:1 matched pairs without replacement, downstream outcome analyses should generally account for the paired design. Depending on the outcome type, this may include paired analyses, conditional models, or regression methods that explicitly respect the matching structure. Residual covariate adjustment after matching may be considered when small residual imbalances remain.

Sigma treats matching primarily as a design-stage procedure and emphasizes transparent reporting of overlap and balance. Correct outcome analysis remains the responsibility of the analyst and depends on the endpoint and estimand of interest.

Assumptions and limitations

• No unmeasured confounding: validity depends on whether all important confounders were measured and included in the propensity score model.
• Overlap: limited overlap can leave treated subjects unmatched and changes the target population to treated subjects within the overlap region.
• Model specification: propensity scores are estimated by logistic regression. In propensity score work, the key criterion is achieved covariate balance rather than predictive discrimination alone.
• Complete-case analysis: excluding incomplete observations may reduce sample size and may bias results if missingness is informative.
• Greedy matching: nearest-neighbor matching without replacement is deterministic and practical, but not globally optimal.
• Inference after matching: valid outcome analyses should account for the matched-pair structure rather than treating the matched sample as an ordinary independent sample.

References

1. Rosenbaum PR, Rubin DB. The central role of the propensity score in observational studies for causal effects. Biometrika. 1983.
2. Rosenbaum PR, Rubin DB. Constructing a control group using multivariate matched sampling methods that incorporate the propensity score. The American Statistician. 1985.
3. Stuart EA. Matching methods for causal inference: A review and a look forward. Statistical Science. 2010.
4. Austin PC. An introduction to propensity score methods for reducing the effects of confounding in observational studies. Multivariate Behavioral Research. 2011.
5. Austin PC. Balance diagnostics for comparing the distribution of baseline covariates between treatment groups in propensity-score matched samples. Statistics in Medicine. 2009.
6. Ho DE, Imai K, King G, Stuart EA. MatchIt: Nonparametric preprocessing for parametric causal inference. Journal of Statistical Software. 2011.

Survival Analysis — Kaplan–Meier, Log-Rank Test, and Cox Regression

Purpose and scope

Sigma analyzes classical right-censored time-to-event endpoints using Kaplan–Meier estimation, the unweighted log-rank test, and Cox proportional hazards regression. The module is table-focused and does not create survival graphics. It currently accepts CSV input and supports one observation per subject, a binary event indicator, optional groups, and time-fixed metric or categorical Cox covariates.

Input, event coding, and exclusions

Follow-up time must be numeric, finite, and non-negative. The status variable must contain exactly two non-missing levels; the user explicitly assigns one level to the event and the other to right censoring. Rows with missing or invalid time/status values are excluded from all survival calculations. A missing group value remains in the Overall Kaplan–Meier estimate but is excluded from group-specific estimates and the log-rank test.

Before risk sets are formed, times that differ only by floating-point round-off are consolidated using the same default tolerance rule as R survival::aeqSurv. Consecutive unique times are treated as a near-tie when their absolute difference, or that difference divided by the mean absolute unique time, is no greater than the square root of machine precision. The left-most time represents the resulting block. This prevents mathematically tied times created by preceding calculations from being split by insignificant binary floating-point differences.

Cox models use a common complete-case cohort across all selected covariates so that univariable and multivariable estimates are based on the same observations. Metric covariates require a complete numeric value. Categorical covariates require a non-missing level and are represented by indicator variables relative to the selected reference level. A variable with at most two distinct non-missing levels is automatically treated as categorical; otherwise it is metric when at least 80% of non-missing values are numeric. Every covariate type can be overridden manually. When Cox stratification or cluster-robust inference is requested, the corresponding strata and cluster variables must also be complete on this common analysis cohort.

Kaplan–Meier estimator

At an ordered event time \(t_j\), let \(n_j\) be the number at risk immediately before \(t_j\) and \(d_j\) the number of events at that time. Events and censorings tied at \(t_j\) are all included in \(n_j\); the event block is applied before those observations leave the risk set. The product-limit estimator is

\[ \hat S(t)=\prod_{t_j\le t}\left(1-\frac{d_j}{n_j}\right). \]

Greenwood's cumulative variance component and the standard error of \(\hat S(t)\) are

\[ \hat G(t)=\sum_{t_j\le t}\frac{d_j}{n_j(n_j-d_j)},\qquad SE\{\hat S(t)\}=\hat S(t)\sqrt{\hat G(t)}. \]

Pointwise 95% confidence intervals use the complementary log-log transformation. With \(\eta=\log[-\log\{\hat S(t)\}]\), the transformed standard error is \(\sqrt{\hat G(t)}/|\log\{\hat S(t)\}|\); the interval is back-transformed to the survival scale. Boundary estimates of 0 or 1 have no finite complementary log-log transform; their confidence limits are therefore reported as not available. These are pointwise intervals, not simultaneous bands.

Median survival is the time at which \(\hat S(t)\) reaches 0.5. If the curve equals exactly 0.5 over an interval, Sigma follows the standard survival convention and reports the midpoint of that interval. Its confidence limits are obtained by inverting the pointwise confidence limits. If a curve or one of its confidence limits never reaches 0.5, the corresponding median is reported as not available. Median follow-up is estimated by reverse Kaplan–Meier: censorings are treated as events and endpoint events as censored observations.

Fixed-time survival tables

At each user-selected time, Sigma reports the right-continuous Kaplan–Meier estimate, its standard error and 95% confidence interval, the number still at risk immediately before that time, and the cumulative numbers of endpoint events and censorings through that time. Values beyond the last observed follow-up are not extrapolated.

Log-rank test

When at least two non-missing groups are available, Sigma performs the unweighted omnibus log-rank test. At each pooled event time, the expected events in group \(g\) are \(e_{gj}=d_j n_{gj}/n_j\). Observed-minus-expected scores are accumulated over time, and their hypergeometric covariance contributions are

\[ V_{gh,j}=\frac{d_j(n_j-d_j)}{n_j-1} \left[\mathbf{1}(g=h)p_{gj}-p_{gj}p_{hj}\right], \qquad p_{gj}=\frac{n_{gj}}{n_j}. \]

The quadratic form has an asymptotic chi-square distribution with \(k-1\) degrees of freedom for \(k\) groups. Sigma additionally reports group-specific observed and expected event counts. The test assumes independent subjects and non-informative censoring. It is most sensitive when group hazards are approximately proportional and can have limited power for crossing survival curves.

Restricted mean survival time

For a user-selected restriction time \(\tau>0\), restricted mean survival time is the area under the Kaplan–Meier curve up to that common horizon,

\[ \widehat{RMST}(\tau)=\int_0^\tau \hat S(t)\,dt, \qquad \widehat{RMTL}(\tau)=\tau-\widehat{RMST}(\tau). \]

Sigma evaluates the integral exactly over the constant Kaplan–Meier steps. The restriction time may not exceed the maximum observed follow-up in any cohort included in the table. With \(A_j(\tau)=\int_{t_j}^{\tau}\hat S(u)\,du\), the Greenwood plug-in variance is

\[ \widehat{Var}\{\widehat{RMST}(\tau)\} =\sum_{t_j\le\tau} A_j(\tau)^2\frac{d_j}{n_j(n_j-d_j)}. \]

Pointwise 95% intervals use the normal approximation and are not truncated to the logical interval \([0,\tau]\). For independent groups, Sigma reports each non-reference group minus the selected reference: the difference standard error is the square root of the sum of the two RMST variances. It also reports the corresponding two-sided normal-approximation p-value. Sigma additionally reports the RMST ratio; its interval uses the delta method on the log-ratio scale. A common restriction time is essential because RMST values at different horizons are not directly comparable. These comparison formulas assume disjoint, independent groups; they are not paired or cluster-robust RMST comparisons.

Cox proportional hazards model

For covariate vector \(X_i\), Sigma fits

\[ h(t\mid X_i)=h_0(t)\exp(\beta^\top X_i). \]

Coefficients are estimated by maximizing the Cox partial likelihood with Newton–Raphson iteration. Covariates are centered and scaled internally to improve numerical stability; reported coefficients, standard errors, hazard ratios, and confidence intervals are transformed back to the original covariate units. A continuous hazard ratio corresponds to a one-unit increase. Categorical coefficients compare one level with the selected reference level.

Efron's approximation is the default for tied event times. If \(D_j\) is a tied block of \(d_j\) events and \(R_j\) is the risk set, the block contribution to the partial-likelihood denominator successively removes fractions \(l/d_j\), \(l=0,\ldots,d_j-1\), of the tied-event risk score. Alternatively, the user can select Breslow handling, which uses the unchanged risk-set denominator for all \(d_j\) events in the block. Efron is generally preferable when ties are present; the selectable Breslow method supports established analysis specifications and reproducibility requirements.

Sigma reports \(\hat\beta\), standard error, Wald \(z\), two-sided Wald p-value, \(HR=\exp(\hat\beta)\), and the 95% Wald interval \(\exp\{\hat\beta\pm1.959964\,SE(\hat\beta)\}\). It fits one univariable model per selected covariate and one multivariable model containing all selected covariates. With model-based variance, multi-level categorical terms receive an overall likelihood-ratio test obtained by comparing the full model with the otherwise identical model excluding that term. With robust variance, their overall test is instead a joint robust Wald test of all indicator coefficients for that term.

Stratified baseline hazards and robust variance

With an optional strata variable, Sigma fits \(h_s(t\mid X_i)=h_{0s}(t)\exp(\beta^\top X_i)\). Each stratum has its own unspecified baseline hazard, while regression coefficients are shared. Risk sets and tied-event blocks are formed separately within strata and their log partial likelihoods, scores, and information matrices are summed. The strata variable therefore does not receive a hazard ratio; effects are identified only by comparisons within the same stratum.

Model-based covariance is the inverse observed information. When robust sandwich variance is selected, let \(U_i\) denote the subject-level score residual and let \(U_c=\sum_{i\in c}U_i\) for cluster \(c\). Sigma uses

\[ \widehat{Var}_{robust}(\hat\beta) =I(\hat\beta)^{-1}\left(\sum_c U_cU_c^\top\right)I(\hat\beta)^{-1}. \]

If no cluster variable is selected, every row is treated as its own independent cluster. Coefficient standard errors, confidence intervals, coefficient Wald tests, the model Wald test, and overall factor Wald tests then use this sandwich covariance. The partial likelihood, likelihood-ratio test, score test, and AIC are unchanged and remain model-based quantities. Robust inference relies on a sufficiently large number of independent clusters; Sigma warns when fewer than 20 clusters are available.

Model-level statistics and concordance

The multivariable summary contains the null and fitted log partial likelihood, the likelihood-ratio, Wald, and score chi-square tests, partial-likelihood AIC \(-2\ell+2p\), convergence information, the tie and variance methods, strata and cluster counts, and the number of events per fitted parameter. Events per parameter are descriptive only and are not used as an automatic validity threshold.

Harrell's concordance statistic is the proportion of comparable subject pairs in which the subject with the earlier observed event has the higher Cox linear predictor; tied predictions receive one half credit. An event tied in time with a censored observation remains comparable, whereas pairs of events tied in time are not ordered. In a stratified model, only pairs within the same stratum are comparable. Sigma reports the infinitesimal-jackknife standard error obtained from subject-level pair influences; with clustered robust inference these influences are summed by cluster before the standard error is calculated. Concordance assesses ranking/discrimination and does not measure calibration or establish causal validity.

Proportional hazards diagnostics

For the multivariable Cox model, Sigma tests the proportional hazards assumption for every model term and globally. Event time is transformed using the pooled Kaplan–Meier estimate and centered over events; Cox risk sets continue to respect any selected strata. The diagnostic is an efficient score test for adding covariate-by-transformed-time terms to the fitted Cox model, which is the formal Grambsch–Therneau/Schoenfeld-residual approach. A small p-value is evidence that the corresponding log hazard ratio changes over time; it is not proof of a particular time-varying shape. The test should be interpreted alongside sample size, event count, subject-matter knowledge, and—in a separate graphics application—residual or survival plots.

Warnings and current limitations

• Right censoring is assumed independent/non-informative conditional on modeled covariates.
• Each subject contributes one row and at most one endpoint event; dependence must be represented through the optional cluster-robust variance.
• The Cox model assumes proportional hazards and a linear effect on the log-hazard scale for metric covariates.
• Robust or cluster-robust variance changes uncertainty estimates but does not correct non-proportional hazards, non-linearity, confounding, or coefficient bias.
• Extreme coefficients, singular information, sparse categories, few events, and non-convergence are reported as warnings or controlled model failures.
• Competing event types must not be collapsed silently; they belong in Sigma's Competing Risk Analysis tool.
• The module does not include delayed entry, start-stop data, time-varying covariates, interactions, splines, adjusted survival profiles, or frailty.
• No survival graphics are produced in Sigma Statistics; graphical reporting is handled by Phi Statistics Drawer.

References

1. Kaplan EL, Meier P. Nonparametric estimation from incomplete observations. JASA. 1958.
2. Mantel N. Evaluation of survival data and two new rank order statistics arising in its consideration. Cancer Chemotherapy Reports. 1966.
3. Cox DR. Regression models and life-tables. JRSS Series B. 1972.
4. Efron B. The efficiency of Cox's likelihood function for censored data. JASA. 1977.
5. Lin DY, Wei LJ. The robust inference for the Cox proportional hazards model. JASA. 1989.
6. Grambsch PM, Therneau TM. Proportional hazards tests and diagnostics based on weighted residuals. Biometrika. 1994.
7. Schemper M, Smith TL. A note on quantifying follow-up in studies of failure time. Controlled Clinical Trials. 1996.
8. Royston P, Parmar MKB. Restricted mean survival time: an alternative to the hazard ratio for the design and analysis of randomized trials with a time-to-event outcome. BMC Medical Research Methodology. 2013.

Competing Risk Analysis — CIFs, Gray Test, and Fine & Gray Regression

Purpose and scope

The Sigma competing risks module provides (i) non-parametric estimation of cumulative incidence functions (CIFs) for each cause with pointwise confidence intervals, (ii) an unadjusted two-group comparison of CIFs using Gray’s test (ρ = 0) implemented to match cmprsk::cuminc, and (iii) regression modeling of covariate effects on the CIF via the Fine & Gray proportional subdistribution hazards model (effects reported as subdistribution hazard ratios, SHR).

In competing risks settings, multiple mutually exclusive event types can occur; occurrence of one type prevents occurrence of the others. CIFs quantify absolute risk for a specific cause over time. Gray’s test compares CIF curves between groups. Fine & Gray regression models covariate effects on the CIF through the subdistribution hazard.

Data structure and preprocessing

The input must contain at least:

• A time variable timei (numeric, non-negative), representing follow-up time for individual i.
• A status variable statusi with status = 0 indicating censoring and positive integers status ∈ {1,2,…} indicating mutually exclusive causes of failure.
• Optionally, a grouping variable. Group-specific CIFs can be estimated for every observed level; Gray’s test and Fine & Gray regression in the current version require exactly two non-missing levels.
• For Fine & Gray regression, the user selects the reference group, the cause or causes of interest, and whether the censoring distribution is estimated separately by group or pooled.

Records with missing time or status are excluded. Non-numeric or negative times and non-integer or negative status codes stop the analysis with an input error instead of being silently reinterpreted. The censoring code is fixed to 0; positive integers are distinct causes. Rows with a missing group remain in the overall CIF and are excluded from group-specific CIFs, Gray’s test, and Fine & Gray regression. The current regression model contains one binary group indicator and no additional covariates.

Event-time grid and basic risk set quantities

Let t1 < t2 < … < tm be the ordered distinct times at which at least one event (any cause > 0) occurs. For each event time tj, define:

• nj: number at risk just before tj,
• drj: number of events from cause r at tj,
• dj = Σr drj: total events (all causes) at tj.

Ties (multiple events at the same time) are handled by aggregating counts at each distinct time. Individuals with time = tj are considered at risk just prior to tj and are removed from the risk set after processing the tied block.

Overall survival for “any event”

Sigma computes the Kaplan–Meier survival function for the composite endpoint “any event” (any status > 0), overall and (if grouping is present) within each group. With S(t0)=1, the product-limit estimator on the event-time grid is:

\[ \hat S(t_j) = \hat S(t_{j-1})\left(1-\frac{d_j}{n_j}\right), \qquad \hat S(t_j-) = \hat S(t_{j-1}). \]

Cumulative incidence functions (CIFs)

For cause r, the CIF is the absolute risk of failing from cause r by time t:

\[ I_r(t)=\Pr(T\le t,\;R=r). \]

Sigma uses the standard non-parametric Aalen–Johansen–type estimator on the event-time grid:

\[ \hat I_r(t_\ell)=\sum_{j:\,t_j\le t_\ell} \hat S(t_j-)\,\frac{d_{rj}}{n_j}. \]

The estimators satisfy \(\sum_{r=1}^k \hat I_r(t)=1-\hat S(t)\), i.e. the sum of cause-specific cumulative incidence equals one minus survival for “any event”.

Variance and pointwise confidence intervals for CIFs

Sigma reports pointwise standard errors and confidence intervals for each \(\hat I_r(t)\). Standard errors use the asymptotic Aalen variance returned by cmprsk::cuminc. Sigma follows the finite-sample and tied-event recursions in cmprsk's cinc routine, including separate contributions from the cause of interest and competing events. The recurrence is evaluated without display rounding; only negligible negative values caused by floating-point arithmetic are truncated to zero. Standard errors are SE = √Var.

Note: Sigma reports pointwise intervals at each time point. These are not simultaneous confidence bands.

Pointwise confidence intervals are constructed using a log(-log) transformation applied directly to the CIF estimate. For a fixed cause \(r\) and time \(t\), let \(\hat I=\hat I_r(t)\) and \(SE=\sqrt{\widehat{\mathrm{Var}}(\hat I)}\). If \(0<\hat I<1\), define:

\[ \eta=\ln\!\bigl(-\ln(\hat I)\bigr). \]

Using the delta method, the standard error on the transformed scale is approximated by:

\[ SE_\eta \approx \left|\frac{1}{\hat I\,\ln(\hat I)}\right|\,SE \;=\;\frac{SE}{-\hat I\,\ln(\hat I)} \quad (\text{since } \ln(\hat I)<0). \]

With confidence level \(1-\alpha\) (default 95%), let \(z=\Phi^{-1}(1-\alpha/2)\). The CI is computed by forming \(\eta_L=\eta-zSE_\eta\) and \(\eta_U=\eta+zSE_\eta\) and back-transforming. An equivalent closed form used operationally is:

\[ \hat I_L = \hat I^{\exp\!\left(-z\,\frac{SE}{\hat I\ln(\hat I)}\right)}, \qquad \hat I_U = \hat I^{\exp\!\left(+z\,\frac{SE}{\hat I\ln(\hat I)}\right)}. \]

If \(\hat I\) is at the boundary (0 or 1) so the transform is undefined and the estimated variance is positive, Sigma falls back to a normal approximation on the original scale:

\[ \hat I \pm z\,SE, \]

with truncation to [0,1]. If the boundary estimate has zero variance, the corresponding interval is degenerate at the estimate.

Gray’s test for group differences in CIFs (two groups; ρ = 0)

If a two-level grouping variable is provided, Sigma computes group-specific CIFs and performs Gray’s test for each cause to compare CIF curves between the two groups (unadjusted comparison).

Null hypothesis

For a given cause \(j\), Gray’s test assesses:

\[ H_0:\; I_{j,1}(t)=I_{j,2}(t)\;\;\text{for all }t, \]

Implementation aligned to cmprsk::cuminc

Sigma implements Gray’s test in the same computational form used by the R package cmprsk (cuminc), i.e. the original algorithm as implemented in the Fortran routines commonly referred to as crstm/crst. The test is evaluated on tied event-time blocks and uses the recursive construction of the score and its variance that accounts for (i) the cause of interest and (ii) competing-event contributions, including tie corrections when multiple events occur at the same time.

The weight parameter is fixed to the unweighted case (ρ = 0), corresponding to a constant weight function. The resulting test statistic for each cause is reported as a one-degree-of-freedom chi-square:

The upper-tail p-value is evaluated directly as erfc(sqrt(chi-square / 2)). Direct tail evaluation avoids the cancellation and approximation error that occurs when subtracting an approximate normal CDF from one.

\[ \chi^2 = \frac{U^2}{\widehat{\mathrm{Var}}(U)} \;\;\sim\;\; \chi^2(1)\quad \text{(asymptotically under }H_0\text{)}. \]

Scope note: Gray’s test in Sigma is restricted to exactly two groups and the unweighted case (ρ = 0). Like most large-sample tests, the \(\chi^2\) approximation can be unreliable with very small numbers of events for the cause of interest.

Fine & Gray regression (subdistribution hazard model)

In addition to unadjusted CIF estimation and Gray’s test, Sigma implements the Fine & Gray proportional subdistribution hazards model for a chosen cause of interest j. The model targets covariate effects on the CIF via the subdistribution hazard.

Subdistribution hazard and model

\[ \lambda_j^*(t\mid Z) = \lambda_{j0}^*(t)\,\exp(\beta^\top Z), \]

where Z is the covariate vector and \(\exp(\beta)\) are subdistribution hazard ratios (SHR).

IPCW for right censoring (with optional censoring strata)

Fine & Gray estimation handles right censoring using inverse probability of censoring weights (IPCW) based on the Kaplan–Meier estimator of the censoring survival function \(\hat G_C(t-)\), constructed by treating censoring as the event and all observed failures as censored in the censoring model. The left-continuous value immediately before each event time is used, including when failures and censorings are tied. The user can estimate the censoring distribution separately within the two modeled groups (cengroup = group) or use one pooled estimate, analogous to cmprsk::crr.

To match the censoring curves constructed by cmprsk through survival::survfit, numerically near-equal censoring-model times are coalesced using aeqSurv's √machine epsilon rule before the censoring Kaplan–Meier estimate is formed. For group-specific censoring distributions, the interpolated survival is set to zero beyond the final fitted follow-up time of the corresponding censoring group, following crr's convention. The observed event-time grid itself is not rounded by this compatibility rule.

\[ Y_i^*(t)= \mathbf{1}(X_i\ge t) +\mathbf{1}(X_i<t,\,\Delta_i\notin\{0,j\}) \frac{\hat G_C(t-)}{\hat G_C(X_i-)}, \]

where j is the target cause, 0 denotes censoring, and Yi*(t) is the weighted subdistribution risk-set contribution. Subjects who experience a competing event remain in the risk set after their event with the displayed IPCW weight; subjects censored before t and subjects who already experienced the target event do not.

Undefined censoring weights cause a controlled model failure and are not replaced by an arbitrary positive constant. The validity of IPCW requires independent censoring within the selected pooled or group-specific censoring model.

Estimation and inference

Sigma minimizes the Fine–Gray negative pseudo-log-likelihood using Newton updates, a scaled-gradient convergence criterion of 10−6, at most ten updates, and step-halving when a full Newton step does not provide sufficient improvement. These conventions follow cmprsk::crr. The robust sandwich covariance includes both subject-level estimating-equation contributions and the influence of estimating the censoring distribution. After convergence, Sigma reports \(\hat\beta\), robust standard errors, Wald z statistics and p-values, SHR, and 95% Wald confidence intervals.

Before fitting, Sigma checks target-cause events in both groups. If either group has no event of the selected cause, the group coefficient is classified as non-estimable because of complete separation. Failed convergence, singular information, invalid censoring weights, and invalid robust variance suppress coefficients and inferential output rather than returning a numerically extreme estimate.

Outputs and exports

• Cause-specific CIF tables (overall and, if group provided, by group) on the observed event-time grid,
• Point estimates, standard errors, and pointwise confidence intervals for CIFs,
• Gray’s test statistic and p-value per cause (two groups; ρ = 0, aligned to cmprsk::cuminc),
• Fine & Gray regression results for the selected cause or all causes, including reference group, events by group, convergence status, censoring method, SHR, robust SE, Wald test, and CI,
• Optional exports of CIF tables and model summaries for documentation and downstream analyses.

Assumptions and limitations

• Competing events are assumed mutually exclusive and correctly coded (one cause per subject or censoring).
• Right censoring is assumed independent/non-informative; Fine & Gray uses IPCW based on this assumption.
• Pointwise CIF confidence intervals are provided (not simultaneous bands).
• Gray’s test is implemented only for exactly two groups and only for the unweighted case (ρ = 0).
• Fine & Gray regression assumes proportional subdistribution hazards; diagnostics/time-varying effects are not included.
• Fine & Gray regression currently supports one binary group indicator; multivariable and time-varying covariates are not included.
• With few target-cause events, Wald inference and the proportional subdistribution-hazards model remain large-sample procedures even when numerical convergence is achieved.

References

1. Gray RJ. A class of K-sample tests for comparing the cumulative incidence of a competing risk. Annals of Statistics. 1988.
2. Fine JP, Gray RJ. A proportional hazards model for the subdistribution of a competing risk. JASA. 1999.
3. Choudhury JB. Non-parametric confidence interval estimation for competing risks analysis. Statistics in Medicine. 2002.
4. Scrucca L, Santucci A, Aversa F. Competing risk analysis using R: an easy guide for clinicians. Bone Marrow Transplantation. 2007.
5. cmprsk R package documentation (cuminc, crr) as methodological reference for unadjusted CIF/Gray output and Fine & Gray regression reporting.