Analysis¶

Statistical analysis of collected experimental results.

ANOVA¶

industrialstats.analysis.anova ¶

ANOVA analysis for experimental designs.

ANOVAAnalysis ¶

ANOVAAnalysis(data: DataFrame, response_column: str)

Perform ANOVA analysis on experimental data.

Initialize an ANOVA analysis instance.

Parameters:

Name	Type	Description	Default
`data`	`DataFrame`	Experimental dataset.	required
`response_column`	`str`	Name of the response variable column.	required

Raises:

Type	Description
`ValueError`	If `response_column` is missing or no rows remain after filtering.

Source code in src/industrialstats/analysis/anova.py

def __init__(self, data: pd.DataFrame, response_column: str):
    """Initialize an ANOVA analysis instance.

    Parameters
    ----------
    data : pd.DataFrame
        Experimental dataset.
    response_column : str
        Name of the response variable column.

    Raises
    ------
    ValueError
        If ``response_column`` is missing or no rows remain after filtering.
    """
    if response_column not in data.columns:
        raise ValueError(f"Response column '{response_column}' not found in data")

    self.data = data.copy()
    self.response = response_column
    self.model: sm.regression.linear_model.RegressionResultsWrapper | None = None
    self.anova_table: pd.DataFrame | None = None

    # Remove any rows with missing response values
    self.data = self.data.dropna(subset=[response_column])

    if len(self.data) == 0:
        raise ValueError("No valid data rows after removing missing values")

fit_model ¶

fit_model(formula: str) -> RegressionResultsWrapper

Fit a linear model using an R-style formula.

Parameters:

Name	Type	Description	Default
`formula`	`str`	Formula string such as `"response ~ factor1 * factor2"`.	required

Returns:

Type	Description
`RegressionResultsWrapper`	Fitted model instance.

Raises:

Type	Description
`ValueError`	If the model fails to fit.

Source code in src/industrialstats/analysis/anova.py

def fit_model(
    self, formula: str
) -> sm.regression.linear_model.RegressionResultsWrapper:
    """Fit a linear model using an R-style formula.

    Parameters
    ----------
    formula : str
        Formula string such as ``"response ~ factor1 * factor2"``.

    Returns
    -------
    RegressionResultsWrapper
        Fitted model instance.

    Raises
    ------
    ValueError
        If the model fails to fit.
    """
    try:
        self.model = ols(formula, data=self.data).fit()
        return self.model
    except (ValueError, np.linalg.LinAlgError) as e:
        raise ValueError(
            f"Error fitting model with formula '{formula}': {e!s}"
        ) from e

anova_table_calculation ¶

anova_table_calculation(typ: int = 2) -> DataFrame

Compute the ANOVA table for the fitted model.

Parameters:

Name	Type	Description	Default
`typ`	`int`	ANOVA type (1, 2, or 3). Defaults to 2.	`2`

Returns:

Type	Description
`DataFrame`	Table with sums of squares, degrees of freedom and statistics.

Raises:

Type	Description
`ValueError`	If no model has been fitted or calculation fails.

Source code in src/industrialstats/analysis/anova.py

def anova_table_calculation(self, typ: int = 2) -> pd.DataFrame:
    """Compute the ANOVA table for the fitted model.

    Parameters
    ----------
    typ : int, optional
        ANOVA type (1, 2, or 3). Defaults to 2.

    Returns
    -------
    pd.DataFrame
        Table with sums of squares, degrees of freedom and statistics.

    Raises
    ------
    ValueError
        If no model has been fitted or calculation fails.
    """
    if self.model is None:
        raise ValueError("Model not fitted. Call fit_model() first.")

    try:
        self.anova_table = anova_lm(self.model, typ=typ)

        # Add additional columns for clarity
        self.anova_table["Mean_Square"] = (
            self.anova_table["sum_sq"] / self.anova_table["df"]
        )

        # Calculate effect sizes (eta-squared and partial eta-squared)
        total_ss = self.anova_table["sum_sq"].sum()
        self.anova_table["Eta_Squared"] = self.anova_table["sum_sq"] / total_ss

        # Partial eta-squared (for Type II and III)
        if typ in [2, 3]:
            residual_ss = (
                self.anova_table.loc["Residual", "sum_sq"]
                if "Residual" in self.anova_table.index
                else 0.0
            )
            # pandas-stubs types a .loc scalar as a broad union covering
            # dates and bytes; an ANOVA sum of squares is always numeric.
            error_ss = float(residual_ss)  # type: ignore[arg-type]
            self.anova_table["Partial_Eta_Squared"] = self.anova_table["sum_sq"] / (
                self.anova_table["sum_sq"] + error_ss
            )

        # Add significance stars
        def add_significance_stars(p_value):
            if pd.isna(p_value):
                return ""
            if p_value < 0.001:
                return "***"
            if p_value < 0.01:
                return "**"
            if p_value < 0.05:
                return "*"
            if p_value < 0.1:
                return "."
            return ""

        self.anova_table["Significance"] = self.anova_table["PR(>F)"].apply(
            add_significance_stars
        )

        return self.anova_table

    except (ValueError, np.linalg.LinAlgError) as e:
        raise ValueError(f"Error calculating ANOVA table: {e!s}") from e

multiple_comparisons ¶

multiple_comparisons(factor: str, method: str = 'tukey', alpha: float = 0.05) -> DataFrame

Run multiple comparison tests on a factor.

Parameters:

Name	Type	Description	Default
`factor`	`str`	Factor name for pairwise comparisons.	required
`method`	`str`	Comparison method (`"tukey"`, `"bonferroni"`, `"holm"`). Defaults to `"tukey"`.	`'tukey'`
`alpha`	`float`	Family-wise error rate. Defaults to 0.05.	`0.05`

Returns:

Type	Description
`DataFrame`	Pairwise comparison results.

Raises:

Type	Description
`ValueError`	If the factor is not present in the data or the method is unsupported.

Source code in src/industrialstats/analysis/anova.py

def multiple_comparisons(
    self, factor: str, method: str = "tukey", alpha: float = 0.05
) -> pd.DataFrame:
    """Run multiple comparison tests on a factor.

    Parameters
    ----------
    factor : str
        Factor name for pairwise comparisons.
    method : str, optional
        Comparison method (``"tukey"``, ``"bonferroni"``, ``"holm"``). Defaults to
        ``"tukey"``.
    alpha : float, optional
        Family-wise error rate. Defaults to 0.05.

    Returns
    -------
    pd.DataFrame
        Pairwise comparison results.

    Raises
    ------
    ValueError
        If the factor is not present in the data or the method is unsupported.
    """
    if factor not in self.data.columns:
        raise ValueError(f"Factor '{factor}' not found in data")

    if method.lower() == "tukey":
        from statsmodels.stats.multicomp import pairwise_tukeyhsd

        mc_result = pairwise_tukeyhsd(
            self.data[self.response], self.data[factor], alpha=alpha
        )

        # Convert to DataFrame
        mc_df = pd.DataFrame(
            {
                "Group1": mc_result.groupsunique[
                    mc_result._multicomp.pairindices[0]
                ],
                "Group2": mc_result.groupsunique[
                    mc_result._multicomp.pairindices[1]
                ],
                "Mean_Diff": mc_result.meandiffs,
                "P_adj": mc_result.pvalues,
                "Lower_CI": mc_result.confint[:, 0],
                "Upper_CI": mc_result.confint[:, 1],
                "Reject_H0": mc_result.reject,
            }
        )

        return mc_df

    if method.lower() == "bonferroni":
        # Bonferroni correction
        groups = self.data[factor].unique()
        n_comparisons = len(groups) * (len(groups) - 1) // 2
        bonferroni_alpha = alpha / n_comparisons

        results = []
        for i, group1 in enumerate(groups):
            for group2 in groups[i + 1 :]:
                data1 = self.data[self.data[factor] == group1][self.response]
                data2 = self.data[self.data[factor] == group2][self.response]

                # Perform t-test
                _t_stat, p_val = stats.ttest_ind(data1, data2)

                # Bonferroni adjusted p-value
                p_adj = min(p_val * n_comparisons, 1.0)

                # Confidence interval for difference
                diff = data1.mean() - data2.mean()
                pooled_se = np.sqrt(
                    (data1.var() / len(data1)) + (data2.var() / len(data2))
                )
                df = len(data1) + len(data2) - 2
                t_critical = stats.t.ppf(1 - bonferroni_alpha / 2, df)
                margin_error = t_critical * pooled_se

                results.append(
                    {
                        "Group1": group1,
                        "Group2": group2,
                        "Mean_Diff": diff,
                        "P_adj": p_adj,
                        "Lower_CI": diff - margin_error,
                        "Upper_CI": diff + margin_error,
                        "Reject_H0": p_adj < alpha,
                    }
                )

        return pd.DataFrame(results)

    raise NotImplementedError(f"Method '{method}' not implemented")

residual_analysis ¶

residual_analysis() -> dict[str, ndarray]

Perform comprehensive residual analysis.

Returns:

Type	Description
`Dict[str, ndarray]`	Dictionary containing various residual statistics.

Source code in src/industrialstats/analysis/anova.py

def residual_analysis(self) -> dict[str, np.ndarray]:
    """Perform comprehensive residual analysis.

    Returns
    -------
    Dict[str, np.ndarray]
        Dictionary containing various residual statistics.
    """
    if self.model is None:
        raise ValueError("Model not fitted. Call fit_model() first.")

    residuals = self.model.resid
    fitted_values = self.model.fittedvalues

    # Standardized residuals
    mse = self.model.mse_resid
    standardized_residuals = residuals / np.sqrt(mse)

    # Studentized residuals
    influence = self.model.get_influence()
    leverage = influence.hat_matrix_diag
    studentized_residuals = residuals / (np.sqrt(mse * (1 - leverage)))

    # Externally studentized residuals
    externally_studentized = influence.resid_studentized_external

    # Cook's distance
    cooks_d = influence.cooks_distance[0]

    return {
        "residuals": residuals,
        "fitted_values": fitted_values,
        "standardized_residuals": standardized_residuals,
        "studentized_residuals": studentized_residuals,
        "externally_studentized_residuals": externally_studentized,
        "leverage": leverage,
        "cooks_distance": cooks_d,
    }

assumptions_tests ¶

assumptions_tests() -> dict[str, dict[str, Any]]

Test ANOVA assumptions.

Returns:

Type	Description
`Dict[str, Dict[str, Any]]`	Test results for normality, homogeneity of variance, and independence.

Source code in src/industrialstats/analysis/anova.py

def assumptions_tests(self) -> dict[str, dict[str, Any]]:
    """Test ANOVA assumptions.

    Returns
    -------
    Dict[str, Dict[str, Any]]
        Test results for normality, homogeneity of variance, and independence.
    """
    if self.model is None:
        raise ValueError("Model not fitted. Call fit_model() first.")

    residuals = self.model.resid
    results = {}

    # 1. Normality test (Shapiro-Wilk)
    try:
        shapiro_stat, shapiro_p = stats.shapiro(residuals)
        results["normality"] = {
            "test": "Shapiro-Wilk",
            "statistic": shapiro_stat,
            "p_value": shapiro_p,
            "assumption_met": shapiro_p > 0.05,
            "interpretation": (
                "Residuals are normally distributed"
                if shapiro_p > 0.05
                else "Residuals may not be normally distributed"
            ),
        }
    except (ValueError, np.linalg.LinAlgError) as e:
        results["normality"] = {
            "test": "Shapiro-Wilk",
            "error": str(e),
            "assumption_met": None,
        }

    # 2. Homogeneity of variance (Levene's test)
    try:
        # Get factor columns (exclude response and metadata)
        factor_cols = [
            col
            for col in self.data.columns
            if col not in [self.response, "RunID", "RunOrder", "Replicate"]
        ]

        if factor_cols:
            # Use first factor for Levene's test
            main_factor = factor_cols[0]
            groups = [
                self.data[self.data[main_factor] == level][self.response].values
                for level in self.data[main_factor].unique()
            ]

            levene_stat, levene_p = stats.levene(*groups)
            results["homogeneity"] = {
                "test": "Levene",
                "statistic": levene_stat,
                "p_value": levene_p,
                "assumption_met": levene_p > 0.05,
                "interpretation": (
                    "Variances are homogeneous"
                    if levene_p > 0.05
                    else "Variances may be heterogeneous"
                ),
            }
        else:
            results["homogeneity"] = {
                "test": "Levene",
                "error": "No factors found for testing",
                "assumption_met": None,
            }
    except (ValueError, np.linalg.LinAlgError) as e:
        results["homogeneity"] = {
            "test": "Levene",
            "error": str(e),
            "assumption_met": None,
        }

    # 3. Independence test (Durbin-Watson)
    try:
        from statsmodels.stats.diagnostic import durbin_watson

        dw_stat = durbin_watson(residuals)

        # DW statistic interpretation
        if 1.5 <= dw_stat <= 2.5:
            independence_met = True
            interpretation = "No evidence of autocorrelation"
        else:
            independence_met = False
            interpretation = "Possible autocorrelation detected"

        results["independence"] = {
            "test": "Durbin-Watson",
            "statistic": dw_stat,
            "assumption_met": independence_met,
            "interpretation": interpretation,
            "note": "Values around 2 indicate no autocorrelation",
        }
    except (ValueError, ImportError) as e:
        results["independence"] = {
            "test": "Durbin-Watson",
            "error": str(e),
            "assumption_met": None,
        }

    return results

model_summary ¶

model_summary() -> dict[str, Any]

Get comprehensive model summary.

Returns:

Type	Description
`dict`	Model fit statistics and summary information.

Source code in src/industrialstats/analysis/anova.py

def model_summary(self) -> dict[str, Any]:
    """Get comprehensive model summary.

    Returns
    -------
    dict
        Model fit statistics and summary information.
    """
    if self.model is None:
        raise ValueError("Model not fitted. Call fit_model() first.")

    return {
        "r_squared": self.model.rsquared,
        "adj_r_squared": self.model.rsquared_adj,
        "f_statistic": self.model.fvalue,
        "f_pvalue": self.model.f_pvalue,
        "mse": self.model.mse_resid,
        "rmse": np.sqrt(self.model.mse_resid),
        "aic": self.model.aic,
        "bic": self.model.bic,
        "n_observations": self.model.nobs,
        "df_residuals": self.model.df_resid,
        "df_model": self.model.df_model,
    }

contrast_analysis ¶

contrast_analysis(contrasts: dict[str, list[float]], factor: str) -> DataFrame

Perform contrast analysis.

Parameters:

Name	Type	Description	Default
`contrasts`	`dict[str, list[float]]`	Mapping of contrast names to coefficient vectors.	required
`factor`	`str`	Factor name for the contrasts.	required

Returns:

Type	Description
`DataFrame`	Contrast analysis results.

Source code in src/industrialstats/analysis/anova.py

def contrast_analysis(
    self, contrasts: dict[str, list[float]], factor: str
) -> pd.DataFrame:
    """Perform contrast analysis.

    Parameters
    ----------
    contrasts : dict[str, list[float]]
        Mapping of contrast names to coefficient vectors.
    factor : str
        Factor name for the contrasts.

    Returns
    -------
    pd.DataFrame
        Contrast analysis results.
    """
    if self.model is None:
        raise ValueError("Model not fitted. Call fit_model() first.")

    if factor not in self.data.columns:
        raise ValueError(f"Factor '{factor}' not found in data")

    # Get factor levels and means
    factor_levels = sorted(self.data[factor].unique())
    group_means = [
        self.data[self.data[factor] == level][self.response].mean()
        for level in factor_levels
    ]
    group_ns = [
        len(self.data[self.data[factor] == level]) for level in factor_levels
    ]

    mse = self.model.mse_resid
    results = []

    for contrast_name, coefficients in contrasts.items():
        if len(coefficients) != len(factor_levels):
            raise ValueError(
                f"Contrast '{contrast_name}' must have {len(factor_levels)} coefficients"
            )

        # Calculate contrast value
        contrast_value = sum(
            c * m for c, m in zip(coefficients, group_means, strict=True)
        )

        # Calculate standard error
        se_squared = mse * sum(
            c**2 / n for c, n in zip(coefficients, group_ns, strict=True)
        )
        se = np.sqrt(se_squared)

        # Calculate t-statistic and p-value
        t_stat = contrast_value / se
        df = self.model.df_resid
        p_value = 2 * (1 - stats.t.cdf(abs(t_stat), df))

        # Confidence interval
        t_critical = stats.t.ppf(0.975, df)
        ci_lower = contrast_value - t_critical * se
        ci_upper = contrast_value + t_critical * se

        results.append(
            {
                "Contrast": contrast_name,
                "Value": contrast_value,
                "SE": se,
                "t_statistic": t_stat,
                "p_value": p_value,
                "CI_Lower": ci_lower,
                "CI_Upper": ci_upper,
                "Significant": p_value < 0.05,
            }
        )

    return pd.DataFrame(results)

power_analysis_post_hoc ¶

power_analysis_post_hoc(alpha: float = 0.05) -> dict[str, float]

Calculate observed power for each effect in the model.

Parameters:

Name	Type	Description	Default
`alpha`	`float`	Significance level. Defaults to 0.05.	`0.05`

Returns:

Type	Description
`dict[str, float]`	Observed power for each effect.

Source code in src/industrialstats/analysis/anova.py

def power_analysis_post_hoc(self, alpha: float = 0.05) -> dict[str, float]:
    """Calculate observed power for each effect in the model.

    Parameters
    ----------
    alpha : float, optional
        Significance level. Defaults to 0.05.

    Returns
    -------
    dict[str, float]
        Observed power for each effect.
    """
    if self.anova_table is None:
        raise ValueError(
            "ANOVA table not calculated. Call anova_table_calculation() first."
        )

    from scipy.stats import ncf

    powers = {}

    for effect in self.anova_table.index:
        if effect != "Residual" and "F" in self.anova_table.columns:
            f_stat = self.anova_table.loc[effect, "F"]
            df1 = self.anova_table.loc[effect, "df"]
            df2 = (
                self.anova_table.loc["Residual", "df"]
                if "Residual" in self.anova_table.index
                else 1
            )

            if not pd.isna(f_stat) and f_stat > 0:
                # Non-centrality parameter
                lambda_nc = f_stat * df1

                # Critical F-value
                f_critical = stats.f.ppf(1 - alpha, df1, df2)

                # Observed power
                power = 1 - ncf.cdf(f_critical, df1, df2, lambda_nc)
                powers[effect] = power

    return powers

mixed_effects_model ¶

mixed_effects_model(fixed_effects: list[str], random_effects: list[str], nested_effects: list[str] | None = None) -> dict[str, Any]

Fit a mixed-effects model with optional nesting.

Parameters:

Name	Type	Description	Default
`fixed_effects`	`list[str]`	Factors treated as fixed effects in the model.	required
`random_effects`	`list[str]`	Random-effect factors. The first entry specifies the grouping variable used for the random intercept.	required
`nested_effects`	`list[str]`	Random effects nested within the main grouping factor. Each entry is treated as a variance component.	`None`

Returns:

Type	Description
`dict`	Contains the following keys: `aic` Model Akaike Information Criterion. `params` Estimated fixed-effect parameters. `random_effects_var` Estimated variances for random effects. `lrt` Likelihood-ratio test statistics for each random effect.

Raises:

Type	Description
`ValueError`	If a specified factor is not present in the data.

Source code in src/industrialstats/analysis/anova.py

def mixed_effects_model(
    self,
    fixed_effects: list[str],
    random_effects: list[str],
    nested_effects: list[str] | None = None,
) -> dict[str, Any]:
    """Fit a mixed-effects model with optional nesting.

    Parameters
    ----------
    fixed_effects : list[str]
        Factors treated as fixed effects in the model.
    random_effects : list[str]
        Random-effect factors. The first entry specifies the grouping
        variable used for the random intercept.
    nested_effects : list[str], optional
        Random effects nested within the main grouping factor. Each entry is
        treated as a variance component.

    Returns
    -------
    dict
        Contains the following keys:

        ``aic``
            Model Akaike Information Criterion.
        ``params``
            Estimated fixed-effect parameters.
        ``random_effects_var``
            Estimated variances for random effects.
        ``lrt``
            Likelihood-ratio test statistics for each random effect.

    Raises
    ------
    ValueError
        If a specified factor is not present in the data.
    """

    for eff in fixed_effects + random_effects + (nested_effects or []):
        if eff not in self.data.columns:
            raise ValueError(f"Factor '{eff}' not found in data")

    import statsmodels.api as sm

    formula = f"{self.response} ~ " + " + ".join(fixed_effects)
    group_col = random_effects[0]

    vc_formula: dict[str, str] = {}
    if nested_effects:
        for effect in nested_effects:
            vc_formula[effect] = f"0 + C({effect})"

    model = sm.MixedLM.from_formula(
        formula,
        self.data,
        groups=group_col,
        re_formula="1",
        vc_formula=vc_formula or None,
    )
    result = model.fit()

    random_vars: dict[str, float] = {
        random_effects[0]: float(result.cov_re.iloc[0, 0])
    }
    if nested_effects:
        for idx, effect in enumerate(nested_effects):
            random_vars[effect] = float(result.vcomp[idx])

    lrt_results: dict[str, dict[str, float]] = {}

    ols_model = sm.OLS.from_formula(formula, self.data).fit()
    lr_stat = 2 * (result.llf - ols_model.llf)
    df = result.df_modelwc - (ols_model.df_model + 1)
    pvalue = stats.chi2.sf(lr_stat, df)
    lrt_results[random_effects[0]] = {
        "lr_stat": float(lr_stat),
        "p_value": float(pvalue),
        "df": int(df),
    }

    if nested_effects:
        for effect in nested_effects:
            reduced_vc = {k: v for k, v in vc_formula.items() if k != effect}
            reduced_model = sm.MixedLM.from_formula(
                formula,
                self.data,
                groups=group_col,
                re_formula="1",
                vc_formula=reduced_vc or None,
            )
            reduced_result = reduced_model.fit()
            lr_stat = 2 * (result.llf - reduced_result.llf)
            df = result.df_modelwc - reduced_result.df_modelwc
            pvalue = stats.chi2.sf(lr_stat, df)
            lrt_results[effect] = {
                "lr_stat": float(lr_stat),
                "p_value": float(pvalue),
                "df": int(df),
            }

    return {
        "aic": float(result.aic),
        "params": result.params.to_dict(),
        "random_effects_var": random_vars,
        "lrt": lrt_results,
    }

unbalanced_anova ¶

unbalanced_anova() -> dict[str, Any]

Perform Type II ANOVA for unbalanced designs.

Returns:

Type	Description
`dict`	Dictionary containing the ANOVA table.

Raises:

Type	Description
`ValueError`	If no model has been fitted.

Source code in src/industrialstats/analysis/anova.py

def unbalanced_anova(self) -> dict[str, Any]:
    """Perform Type II ANOVA for unbalanced designs.

    Returns
    -------
    dict
        Dictionary containing the ANOVA table.

    Raises
    ------
    ValueError
        If no model has been fitted.
    """

    if self.model is None:
        raise ValueError("Model not fitted. Call fit_model() first.")

    table = anova_lm(self.model, typ=2)
    return {"anova_table": table}

nested_anova ¶

nested_anova(nesting_structure: dict[str, str]) -> dict[str, Any]

Perform nested ANOVA for hierarchical designs.

Parameters:

Name	Type	Description	Default
`nesting_structure`	`dict[str, str]`	Mapping of nested factor to its parent factor.	required

Returns:

Type	Description
`dict`	Dictionary with the ANOVA table.

Source code in src/industrialstats/analysis/anova.py

def nested_anova(self, nesting_structure: dict[str, str]) -> dict[str, Any]:
    """Perform nested ANOVA for hierarchical designs.

    Parameters
    ----------
    nesting_structure : dict[str, str]
        Mapping of nested factor to its parent factor.

    Returns
    -------
    dict
        Dictionary with the ANOVA table.
    """

    terms = []
    for child, parent in nesting_structure.items():
        if child not in self.data.columns or parent not in self.data.columns:
            raise ValueError("Factors specified in nesting_structure missing")
        terms.append(f"C({parent})/C({child})")

    formula = f"{self.response} ~ " + " + ".join(terms)
    self.fit_model(formula)
    table = self.anova_table_calculation(typ=2)
    return {"anova_table": table}

repeated_measures_anova ¶

repeated_measures_anova(subject_column: str, within_factors: list[str]) -> dict[str, Any]

Analyze repeated measures designs.

Parameters:

Name	Type	Description	Default
`subject_column`	`str`	Identifier for each experimental unit.	required
`within_factors`	`list[str]`	Factors measured repeatedly.	required

Returns:

Type	Description
`dict`	Dictionary with the ANOVA table.

Raises:

Type	Description
`ValueError`	If specified columns are not in the data.

Source code in src/industrialstats/analysis/anova.py

def repeated_measures_anova(
    self, subject_column: str, within_factors: list[str]
) -> dict[str, Any]:
    """Analyze repeated measures designs.

    Parameters
    ----------
    subject_column : str
        Identifier for each experimental unit.
    within_factors : list[str]
        Factors measured repeatedly.

    Returns
    -------
    dict
        Dictionary with the ANOVA table.

    Raises
    ------
    ValueError
        If specified columns are not in the data.
    """

    if subject_column not in self.data.columns:
        raise ValueError(f"Subject column '{subject_column}' not found")
    for fac in within_factors:
        if fac not in self.data.columns:
            raise ValueError(f"Factor '{fac}' not found in data")

    from statsmodels.stats.anova import AnovaRM

    rm = AnovaRM(
        self.data,
        depvar=self.response,
        subject=subject_column,
        within=within_factors,
    )
    res = rm.fit()
    return {"anova_table": res.anova_table}

Split-plot inference¶

Split-plot experiments have two randomization stages and therefore two error strata. Whole-plot-only treatment terms must be tested against whole-plot error; terms containing subplot factors must be tested against subplot error.

Split-plot randomization, error strata, and inference structure

Source: ../diagrams/split_plot_inference.dot

The classical stratum-specific F tests and the random-intercept mixed model answer complementary questions. A generic OLS residual denominator is not valid for whole-plot treatment effects, and mixed-model Wald tests are not substituted for the classical balanced split-plot F tests.

industrialstats.analysis.split_plot ¶

Inference helpers for balanced complete split-plot experiments.

SplitPlotErrorStrata `dataclass` ¶

SplitPlotErrorStrata(whole_plot_treatments: int, subplot_treatments: int, replicates: int, whole_plots: int, runs: int, whole_plot_error_df: int, subplot_error_df: int)

Balanced split-plot experimental-unit and error-stratum summary.

SplitPlotAnalysis ¶

SplitPlotAnalysis(data: DataFrame, response_column: str, whole_plot_factors: list[str], subplot_factors: list[str], *, whole_plot_column: str = 'WholePlot', replicate_column: str = 'Replicate')

Validate and analyse a balanced complete split-plot experiment.

Source code in src/industrialstats/analysis/split_plot.py

def __init__(
    self,
    data: pd.DataFrame,
    response_column: str,
    whole_plot_factors: list[str],
    subplot_factors: list[str],
    *,
    whole_plot_column: str = "WholePlot",
    replicate_column: str = "Replicate",
) -> None:
    if not whole_plot_factors:
        raise ValueError("At least one whole-plot factor is required")
    if not subplot_factors:
        raise ValueError("At least one subplot factor is required")
    overlap = set(whole_plot_factors) & set(subplot_factors)
    if overlap:
        raise ValueError(
            "Factors cannot be both whole-plot and subplot factors: "
            + ", ".join(sorted(overlap))
        )

    required = {
        response_column,
        whole_plot_column,
        replicate_column,
        *whole_plot_factors,
        *subplot_factors,
    }
    missing = sorted(required - set(data.columns))
    if missing:
        raise ValueError("Missing required column(s): " + ", ".join(missing))

    self.data = data.dropna(subset=[response_column]).copy()
    if self.data.empty:
        raise ValueError("No valid data rows after removing missing responses")

    self.response = response_column
    self.whole_plot_factors = list(whole_plot_factors)
    self.subplot_factors = list(subplot_factors)
    self.whole_plot_column = whole_plot_column
    self.replicate_column = replicate_column

error_strata ¶

error_strata() -> SplitPlotErrorStrata

Return balanced whole-plot and subplot error-stratum degrees of freedom.

Source code in src/industrialstats/analysis/split_plot.py

def error_strata(self) -> SplitPlotErrorStrata:
    """Return balanced whole-plot and subplot error-stratum degrees of freedom."""
    self._validate_whole_plot_units()

    whole_treatment_counts = (
        self.data.groupby(self.whole_plot_factors, observed=True)[
            self.whole_plot_column
        ]
        .nunique()
        .to_numpy()
    )
    if len({int(value) for value in whole_treatment_counts}) != 1:
        raise ValueError(
            "Every whole-plot treatment combination must have the same number "
            "of independent whole-plot replicates"
        )

    replicates = int(whole_treatment_counts[0])
    whole_plot_treatments = int(
        self.data[self.whole_plot_factors].drop_duplicates().shape[0]
    )
    subplot_levels = [
        int(self.data[factor].nunique(dropna=False))
        for factor in self.subplot_factors
    ]
    subplot_treatments = prod(subplot_levels)
    expected_subplot_combinations = int(
        self.data[self.subplot_factors].drop_duplicates().shape[0]
    )
    if expected_subplot_combinations != subplot_treatments:
        raise ValueError(
            "Subplot factors do not form a complete factorial treatment set"
        )

    for whole_plot_id, frame in self.data.groupby(
        self.whole_plot_column, observed=True, sort=False
    ):
        subplot_combinations = frame[self.subplot_factors].drop_duplicates()
        if (
            len(frame) != subplot_treatments
            or len(subplot_combinations) != subplot_treatments
        ):
            raise ValueError(
                f"Whole plot {whole_plot_id!r} must contain exactly one complete "
                "subplot factorial"
            )

    whole_plots = whole_plot_treatments * replicates
    runs = whole_plots * subplot_treatments
    if len(self.data) != runs:
        raise ValueError(
            "Observed run count is inconsistent with a balanced split-plot"
        )

    return SplitPlotErrorStrata(
        whole_plot_treatments=whole_plot_treatments,
        subplot_treatments=subplot_treatments,
        replicates=replicates,
        whole_plots=whole_plots,
        runs=runs,
        whole_plot_error_df=whole_plot_treatments * (replicates - 1),
        subplot_error_df=(
            whole_plot_treatments * (replicates - 1) * (subplot_treatments - 1)
        ),
    )

anova_table ¶

anova_table() -> DataFrame

Return the classical balanced split-plot ANOVA table.

Whole-plot-only treatment terms are tested against WholePlot Error. Terms containing at least one subplot factor are tested against Subplot Error. Sums of squares use orthogonal Helmert-contrast projections, so the decomposition is invariant to row order and to arbitrary treatment labels.

Returns:

Type	Description
`DataFrame`	ANOVA table with source, randomization stratum, degrees of freedom, sums of squares, mean squares, F statistics, p-values, and the denominator error term used for each fixed effect.

Source code in src/industrialstats/analysis/split_plot.py

def anova_table(self) -> pd.DataFrame:
    """Return the classical balanced split-plot ANOVA table.

    Whole-plot-only treatment terms are tested against ``WholePlot Error``.
    Terms containing at least one subplot factor are tested against
    ``Subplot Error``. Sums of squares use orthogonal Helmert-contrast
    projections, so the decomposition is invariant to row order and to
    arbitrary treatment labels.

    Returns
    -------
    pandas.DataFrame
        ANOVA table with source, randomization stratum, degrees of freedom,
        sums of squares, mean squares, F statistics, p-values, and the
        denominator error term used for each fixed effect.
    """
    strata = self.error_strata()
    if strata.whole_plot_error_df <= 0:
        raise ValueError(
            "Whole-plot treatment combinations require replication to estimate "
            "whole-plot error"
        )
    if strata.subplot_error_df <= 0:
        raise ValueError(
            "The subplot stratum requires replication to estimate subplot error"
        )

    all_factors = [*self.whole_plot_factors, *self.subplot_factors]
    level_map = {
        factor: list(pd.unique(self.data[factor])) for factor in all_factors
    }
    all_terms = [
        term
        for order in range(1, len(all_factors) + 1)
        for term in combinations(all_factors, order)
    ]
    whole_plot_factor_set = set(self.whole_plot_factors)
    whole_plot_terms = [
        term for term in all_terms if set(term) <= whole_plot_factor_set
    ]
    subplot_terms = [term for term in all_terms if term not in whole_plot_terms]

    grouped = self.data.groupby(
        self.whole_plot_column,
        observed=True,
        sort=False,
    )
    whole_plot_frame = grouped[self.whole_plot_factors].first().reset_index()
    whole_plot_frame[self.response] = grouped[self.response].mean().to_numpy()

    whole_plot_response = whole_plot_frame[self.response].to_numpy(dtype=float)
    whole_plot_centered = whole_plot_response - whole_plot_response.mean()
    whole_plot_total_ss = (
        float(whole_plot_centered @ whole_plot_centered) * strata.subplot_treatments
    )

    rows: list[dict[str, Any]] = []
    whole_plot_model_ss = 0.0
    for term in whole_plot_terms:
        matrix = self._term_matrix(whole_plot_frame, term, level_map)
        sum_sq = (
            self._projection_sum_of_squares(whole_plot_centered, matrix)
            * strata.subplot_treatments
        )
        whole_plot_model_ss += sum_sq
        rows.append(
            {
                "Source": self._term_name(term),
                "Stratum": "whole_plot",
                "df": int(matrix.shape[1]),
                "sum_sq": sum_sq,
            }
        )

    whole_plot_error_ss = whole_plot_total_ss - whole_plot_model_ss
    if whole_plot_error_ss < -1e-10:
        raise ValueError("Whole-plot decomposition produced negative error SS")
    whole_plot_error_ss = max(0.0, whole_plot_error_ss)
    whole_plot_error_ms = whole_plot_error_ss / strata.whole_plot_error_df

    response = self.data[self.response].to_numpy(dtype=float)
    whole_plot_means = (
        grouped[self.response].transform("mean").to_numpy(dtype=float)
    )
    within_response = response - whole_plot_means
    subplot_total_ss = float(within_response @ within_response)

    subplot_model_ss = 0.0
    subplot_model_df = 0
    for term in subplot_terms:
        matrix = self._term_matrix(self.data, term, level_map)
        sum_sq = self._projection_sum_of_squares(within_response, matrix)
        term_df = int(matrix.shape[1])
        subplot_model_ss += sum_sq
        subplot_model_df += term_df
        rows.append(
            {
                "Source": self._term_name(term),
                "Stratum": "subplot",
                "df": term_df,
                "sum_sq": sum_sq,
            }
        )

    derived_subplot_error_df = (
        len(self.data) - strata.whole_plots - subplot_model_df
    )
    if derived_subplot_error_df != strata.subplot_error_df:
        raise ValueError(
            "Subplot error degrees of freedom disagree with the validated "
            "balanced-design identity"
        )

    subplot_error_ss = subplot_total_ss - subplot_model_ss
    if subplot_error_ss < -1e-10:
        raise ValueError("Subplot decomposition produced negative error SS")
    subplot_error_ss = max(0.0, subplot_error_ss)
    subplot_error_ms = subplot_error_ss / strata.subplot_error_df

    rows.extend(
        [
            {
                "Source": "WholePlot Error",
                "Stratum": "whole_plot_error",
                "df": strata.whole_plot_error_df,
                "sum_sq": whole_plot_error_ss,
            },
            {
                "Source": "Subplot Error",
                "Stratum": "subplot_error",
                "df": strata.subplot_error_df,
                "sum_sq": subplot_error_ss,
            },
        ]
    )

    table = pd.DataFrame(rows)
    table["mean_sq"] = table["sum_sq"] / table["df"]
    table["F"] = np.nan
    table["PR(>F)"] = np.nan
    table["Denominator"] = pd.NA

    for index, row in table.iterrows():
        stratum = row["Stratum"]
        if stratum == "whole_plot":
            denominator_ms = whole_plot_error_ms
            denominator_df = strata.whole_plot_error_df
            denominator_name = "WholePlot Error"
        elif stratum == "subplot":
            denominator_ms = subplot_error_ms
            denominator_df = strata.subplot_error_df
            denominator_name = "Subplot Error"
        else:
            continue

        numerator_ms = float(row["mean_sq"])
        if denominator_ms == 0.0:
            f_statistic = np.inf if numerator_ms > 0.0 else np.nan
            p_value = 0.0 if np.isinf(f_statistic) else np.nan
        else:
            f_statistic = numerator_ms / denominator_ms
            p_value = float(stats.f.sf(f_statistic, int(row["df"]), denominator_df))

        table.loc[index, "F"] = f_statistic
        table.loc[index, "PR(>F)"] = p_value
        table.loc[index, "Denominator"] = denominator_name

    return table

expected_mean_squares ¶

expected_mean_squares(*, whole_plot_variance: float, residual_variance: float) -> dict[str, float]

Return error-stratum EMS values for the random-intercept split-plot model.

Source code in src/industrialstats/analysis/split_plot.py

def expected_mean_squares(
    self,
    *,
    whole_plot_variance: float,
    residual_variance: float,
) -> dict[str, float]:
    """Return error-stratum EMS values for the random-intercept split-plot model."""
    if whole_plot_variance < 0 or residual_variance < 0:
        raise ValueError("Variance components must be non-negative")
    strata = self.error_strata()
    return {
        "whole_plot_error": residual_variance
        + strata.subplot_treatments * whole_plot_variance,
        "subplot_error": residual_variance,
    }

fit_mixed_model ¶

fit_mixed_model(*, reml: bool = True, method: str = 'lbfgs') -> dict[str, Any]

Fit the full fixed-treatment split-plot model with a random whole-plot intercept.

Source code in src/industrialstats/analysis/split_plot.py

def fit_mixed_model(
    self,
    *,
    reml: bool = True,
    method: str = "lbfgs",
) -> dict[str, Any]:
    """Fit the full fixed-treatment split-plot model with a random whole-plot intercept."""
    strata = self.error_strata()
    all_factors = [*self.whole_plot_factors, *self.subplot_factors]
    fixed_terms = " * ".join(self._categorical_term(name) for name in all_factors)
    formula = f"{self._quote(self.response)} ~ {fixed_terms}"

    model = sm.MixedLM.from_formula(
        formula,
        data=self.data,
        groups=self.whole_plot_column,
        re_formula="1",
    )
    result = model.fit(reml=reml, method=method)

    whole_plot_variance = float(result.cov_re.iloc[0, 0])
    residual_variance = float(result.scale)
    return {
        "formula": formula,
        "converged": bool(result.converged),
        "reml": reml,
        "fixed_effects": result.fe_params.to_dict(),
        "whole_plot_variance": whole_plot_variance,
        "residual_variance": residual_variance,
        "error_strata": asdict(strata),
        "expected_mean_squares": self.expected_mean_squares(
            whole_plot_variance=whole_plot_variance,
            residual_variance=residual_variance,
        ),
        "log_likelihood": float(result.llf),
    }

Effects¶

industrialstats.analysis.effects ¶

Effect analysis with canonical two-level factorial contrast semantics.

EffectsAnalysis ¶

EffectsAnalysis(design_matrix: DataFrame, response_data: list[float])

Bases: EffectsAnalysis

Calculate factorial effects with one canonical two-level convention.

Complete balanced two-level factorials use the same -1/+1 orthogonal contrast engine as :meth:industrialstats.designs.factorial.FactorialDesign.calculate_effects. Multi-level and incomplete layouts retain the established analysis paths.

Initialize effect analysis and exclude known design metadata columns.

Source code in src/industrialstats/analysis/effects.py

def __init__(self, design_matrix: pd.DataFrame, response_data: list[float]):
    """Initialize effect analysis and exclude known design metadata columns."""
    stored_orders = design_matrix.attrs.get("factor_level_orders")
    super().__init__(design_matrix, response_data)
    self.factor_names = [
        column
        for column in design_matrix.columns
        if column not in _METADATA_COLUMNS
    ]
    if not self.factor_names:
        raise ValueError("No factor columns found in design matrix")

    self._factor_level_orders: dict[str, list[Any]] | None = None
    if isinstance(stored_orders, dict) and all(
        name in stored_orders for name in self.factor_names
    ):
        self._factor_level_orders = {
            name: list(stored_orders[name]) for name in self.factor_names
        }

calculate_main_effects ¶

calculate_main_effects() -> dict[str, float]

Calculate main effects, using canonical contrasts for complete 2^k data.

Source code in src/industrialstats/analysis/effects.py

def calculate_main_effects(self) -> dict[str, float]:
    """Calculate main effects, using canonical contrasts for complete ``2^k`` data."""
    if not self._is_complete_balanced_two_level_factorial():
        return super().calculate_main_effects()

    return calculate_two_level_factorial_effects(
        self.design_matrix,
        self.response_data,
        self.factor_names,
        max_order=1,
        level_orders=self._factor_level_orders,
    )

calculate_interaction_effects ¶

calculate_interaction_effects(max_order: int = 2) -> dict[str, float]

Calculate interaction effects under the canonical factorial convention.

Source code in src/industrialstats/analysis/effects.py

def calculate_interaction_effects(self, max_order: int = 2) -> dict[str, float]:
    """Calculate interaction effects under the canonical factorial convention."""
    if max_order < 2:
        return {}
    if not self._is_complete_balanced_two_level_factorial():
        return super().calculate_interaction_effects(max_order=max_order)

    effects = calculate_two_level_factorial_effects(
        self.design_matrix,
        self.response_data,
        self.factor_names,
        max_order=max_order,
        level_orders=self._factor_level_orders,
    )
    return {name: effect for name, effect in effects.items() if "*" in name}

Model fitting¶

industrialstats.analysis.model_fitting ¶

Advanced model fitting and selection for experimental data.

ModelFitting ¶

ModelFitting(data: DataFrame, response_column: str)

Advanced model fitting with automatic term selection and validation.

Initialize model fitting.

Parameters:

Name	Type	Description	Default
`data`	`DataFrame`	Experimental data.	required
`response_column`	`str`	Name of the response variable.	required

Source code in src/industrialstats/analysis/model_fitting.py

def __init__(self, data: pd.DataFrame, response_column: str):
    """Initialize model fitting.

    Parameters
    ----------
    data : pandas.DataFrame
        Experimental data.
    response_column : str
        Name of the response variable.
    """
    if response_column not in data.columns:
        raise ValueError(f"Response column '{response_column}' not found")

    self.data = data.copy()
    self.response = response_column
    self.factor_columns = [
        col
        for col in data.columns
        if col
        not in [response_column, "RunID", "RunOrder", "Replicate", "DesignPoint"]
    ]

    # Remove missing values
    self.data = self.data.dropna(subset=[response_column])

    if len(self.data) == 0:
        raise ValueError("No valid data after removing missing values")

    self.fitted_models: dict[str, Any] = {}
    self.model_comparison: pd.DataFrame | None = None

stepwise_selection ¶

stepwise_selection(entry_threshold: float = 0.05, removal_threshold: float = 0.1, max_terms: int | None = None) -> dict[str, Any]

Perform stepwise model selection.

Parameters:

Name	Type	Description	Default
`entry_threshold`	`float`	P-value threshold for entering terms, by default `0.05`.	`0.05`
`removal_threshold`	`float`	P-value threshold for removing terms, by default `0.10`.	`0.1`
`max_terms`	`int`	Maximum number of terms in the model.	`None`

Returns:

Type	Description
`dict`	Stepwise selection results.

Source code in src/industrialstats/analysis/model_fitting.py

def stepwise_selection(
    self,
    entry_threshold: float = 0.05,
    removal_threshold: float = 0.10,
    max_terms: int | None = None,
) -> dict[str, Any]:
    """Perform stepwise model selection.

    Parameters
    ----------
    entry_threshold : float, optional
        P-value threshold for entering terms, by default ``0.05``.
    removal_threshold : float, optional
        P-value threshold for removing terms, by default ``0.10``.
    max_terms : int, optional
        Maximum number of terms in the model.

    Returns
    -------
    dict
        Stepwise selection results.
    """
    # Generate candidate terms
    candidate_terms = self._generate_candidate_terms()

    def _extract_p_value(term: str, pvals: dict[str, float]) -> float | None:
        """Map simplified term names to statsmodels parameter keys."""
        if term in pvals:
            return pvals[term]
        pattern = (
            ":".join([f"C({t})" for t in term.split("*")])
            if "*" in term
            else f"C({term})"
        )
        for key, val in pvals.items():
            if key.startswith(pattern):
                return val
        return None

    if max_terms is None:
        max_terms = min(len(candidate_terms), len(self.data) // 3)

    # Start with intercept-only model
    current_terms = ["Intercept"]
    selection_history = []

    while True:
        improved = False

        # Forward step: try adding terms
        best_addition = None
        best_p_value = float("inf")

        for term in candidate_terms:
            if term not in current_terms and len(current_terms) < max_terms:
                trial_terms = [*current_terms, term]
                try:
                    model_result = self._fit_terms(trial_terms)

                    p_value = _extract_p_value(term, model_result["p_values"])
                    if (
                        p_value is not None
                        and p_value < entry_threshold
                        and p_value < best_p_value
                    ):
                        best_addition = term
                        best_p_value = p_value
                except (ValueError, np.linalg.LinAlgError) as e:
                    logger.debug("Failed to fit trial terms %s: %s", trial_terms, e)
                    continue

        # Add best term if found
        if best_addition is not None:
            current_terms.append(best_addition)
            selection_history.append(
                {
                    "action": "add",
                    "term": best_addition,
                    "p_value": best_p_value,
                    "current_terms": current_terms.copy(),
                }
            )
            improved = True

        # Backward step: try removing terms
        worst_removal = None
        worst_p_value = 0

        for term in current_terms[1:]:  # Skip intercept
            trial_terms = [t for t in current_terms if t != term]
            try:
                model_result = self._fit_terms(current_terms)

                p_value = _extract_p_value(term, model_result["p_values"])
                if (
                    p_value is not None
                    and p_value > removal_threshold
                    and p_value > worst_p_value
                ):
                    worst_removal = term
                    worst_p_value = p_value
            except (ValueError, np.linalg.LinAlgError) as e:
                logger.debug("Failed to evaluate term %s: %s", term, e)
                continue

        # Remove worst term if found
        if worst_removal is not None:
            current_terms.remove(worst_removal)
            selection_history.append(
                {
                    "action": "remove",
                    "term": worst_removal,
                    "p_value": worst_p_value,
                    "current_terms": current_terms.copy(),
                }
            )
            improved = True

        # Stop if no improvement
        if not improved:
            break

    # Fit final model
    final_model = self._fit_terms(current_terms)

    return {
        "selected_terms": current_terms,
        "selection_history": selection_history,
        "final_model": final_model,
        "entry_threshold": entry_threshold,
        "removal_threshold": removal_threshold,
    }

hierarchical_fitting ¶

hierarchical_fitting(max_order: int = 3, significance_level: float = 0.05) -> dict[str, Any]

Fit hierarchical models while respecting effect hierarchy.

Terms are added in increasing order of interaction degree. A term is only considered if all of its lower-order components are already present in the model, enforcing the principle described by Montgomery [1]_. Each candidate term is fit and retained when its p-value is below significance_level.

Parameters:

Name	Type	Description	Default
`max_order`	`int`	Maximum interaction order. For example, `2` fits up to two-factor interactions. Default is `3`.	`3`
`significance_level`	`float`	Significance level for term inclusion. Default is `0.05`.	`0.05`

Returns:

Type	Description
`dict`	Dictionary containing selected terms and fitted model statistics.

all_subsets_selection ¶

all_subsets_selection(criterion: str = 'AIC') -> dict[str, Any]

Perform all possible subsets selection.

Parameters:

Name	Type	Description	Default
`criterion`	`str`	Selection criterion: `'AIC'`, `'BIC'`, `'R2'`, `'R2_adj'`.	`"AIC"`

Returns:

Type	Description
`Dict[str, Any]`	All subsets results.

Source code in src/industrialstats/analysis/model_fitting.py

def all_subsets_selection(self, criterion: str = "AIC") -> dict[str, Any]:
    """Perform all possible subsets selection.

    Parameters
    ----------
    criterion : str, default="AIC"
        Selection criterion: ``'AIC'``, ``'BIC'``, ``'R2'``, ``'R2_adj'``.

    Returns
    -------
    Dict[str, Any]
        All subsets results.
    """
    candidate_terms = self._generate_candidate_terms()
    max_terms = min(len(candidate_terms), len(self.data) // 4)  # Conservative limit

    if len(candidate_terms) > 20:
        warnings.warn(
            "Large number of candidate terms. Consider using stepwise selection.",
            stacklevel=2,
        )

    best_models_by_size = {}
    all_models = []

    # Try all subset sizes
    for subset_size in range(1, max_terms + 1):
        best_criterion = (
            float("inf") if criterion in ["AIC", "BIC"] else float("-inf")
        )
        best_model = None
        best_terms = None

        # Try all combinations of this size
        for term_combination in combinations(candidate_terms, subset_size):
            terms = ["Intercept", *list(term_combination)]

            try:
                model_result = self._fit_terms(terms)
                criterion_value = model_result["model_metrics"][criterion]

                all_models.append(
                    {
                        "terms": terms,
                        "n_terms": len(terms),
                        "criterion_value": criterion_value,
                        "r_squared": model_result["model_metrics"]["R2"],
                        "model_result": model_result,
                    }
                )

                # Check if best for this size
                if criterion in ["AIC", "BIC"]:
                    is_better = criterion_value < best_criterion
                else:
                    is_better = criterion_value > best_criterion

                if is_better:
                    best_criterion = criterion_value
                    best_model = model_result
                    best_terms = terms

            except (ValueError, np.linalg.LinAlgError) as e:
                logger.debug("Failed to fit subset %s: %s", terms, e)
                continue

        if best_model is not None:
            best_models_by_size[subset_size] = {
                "terms": best_terms,
                "model": best_model,
                "criterion_value": best_criterion,
            }

    # Find overall best model
    overall_best = None
    overall_best_criterion = (
        float("inf") if criterion in ["AIC", "BIC"] else float("-inf")
    )

    for _size, model_info in best_models_by_size.items():
        criterion_value = model_info["criterion_value"]

        if criterion in ["AIC", "BIC"]:
            is_better = criterion_value < overall_best_criterion
        else:
            is_better = criterion_value > overall_best_criterion

        if is_better:
            overall_best = model_info
            overall_best_criterion = criterion_value

    return {
        "best_models_by_size": best_models_by_size,
        "overall_best": overall_best,
        "all_models": sorted(all_models, key=lambda x: x["criterion_value"]),
        "criterion": criterion,
    }

cross_validation ¶

cross_validation(model_terms: list[str], k_folds: int = 5, random_state: int | None = None) -> dict[str, Any]

Perform k-fold cross-validation.

Parameters:

Name	Type	Description	Default
`model_terms`	`list of str`	Model terms to validate.	required
`k_folds`	`int`	Number of folds, by default 5.	`5`
`random_state`	`int`	Seed for reproducible splitting.	`None`

Returns:

Type	Description
`dict`	Cross-validation results.

Source code in src/industrialstats/analysis/model_fitting.py

def cross_validation(
    self,
    model_terms: list[str],
    k_folds: int = 5,
    random_state: int | None = None,
) -> dict[str, Any]:
    """Perform k-fold cross-validation.

    Parameters
    ----------
    model_terms : list of str
        Model terms to validate.
    k_folds : int, optional
        Number of folds, by default 5.
    random_state : int, optional
        Seed for reproducible splitting.

    Returns
    -------
    dict
        Cross-validation results.
    """
    from sklearn.model_selection import KFold

    kf = KFold(n_splits=k_folds, shuffle=True, random_state=random_state)

    cv_results = {"fold_results": [], "predictions": [], "actuals": []}

    for fold, (train_idx, test_idx) in enumerate(kf.split(self.data)):
        train_data = self.data.iloc[train_idx]
        test_data = self.data.iloc[test_idx]

        # Fit model on training data
        try:
            # Create temporary ModelFitting object for training data
            train_fitter = ModelFitting(train_data, self.response)
            train_model = train_fitter._fit_terms(model_terms)

            # Predict on test data
            test_predictions = self._predict_with_model(
                train_model, test_data, model_terms
            )
            test_actuals = test_data[self.response].values

            # Calculate fold metrics
            fold_rmse = np.sqrt(np.mean((test_predictions - test_actuals) ** 2))
            fold_mae = np.mean(np.abs(test_predictions - test_actuals))
            fold_r2 = 1 - np.sum((test_actuals - test_predictions) ** 2) / np.sum(
                (test_actuals - np.mean(test_actuals)) ** 2
            )

            cv_results["fold_results"].append(
                {
                    "fold": fold + 1,
                    "rmse": fold_rmse,
                    "mae": fold_mae,
                    "r2": fold_r2,
                    "n_train": len(train_data),
                    "n_test": len(test_data),
                }
            )

            cv_results["predictions"].extend(test_predictions)
            cv_results["actuals"].extend(test_actuals)

        except (ValueError, np.linalg.LinAlgError) as e:
            logger.debug("Cross-validation fold %s failed: %s", fold + 1, e)
            cv_results["fold_results"].append({"fold": fold + 1, "error": str(e)})

    # Calculate overall CV metrics
    if cv_results["predictions"]:
        all_predictions = np.array(cv_results["predictions"])
        all_actuals = np.array(cv_results["actuals"])

        cv_results["overall_rmse"] = np.sqrt(
            np.mean((all_predictions - all_actuals) ** 2)
        )
        cv_results["overall_mae"] = np.mean(np.abs(all_predictions - all_actuals))
        cv_results["overall_r2"] = 1 - np.sum(
            (all_actuals - all_predictions) ** 2
        ) / np.sum((all_actuals - np.mean(all_actuals)) ** 2)

        # Calculate mean and std of fold metrics
        valid_folds = [f for f in cv_results["fold_results"] if "error" not in f]
        if valid_folds:
            cv_results["mean_rmse"] = np.mean([f["rmse"] for f in valid_folds])
            cv_results["std_rmse"] = np.std([f["rmse"] for f in valid_folds])
            cv_results["mean_r2"] = np.mean([f["r2"] for f in valid_folds])
            cv_results["std_r2"] = np.std([f["r2"] for f in valid_folds])

    return cv_results

bootstrap_validation ¶

bootstrap_validation(model_terms: list[str], n_bootstrap: int = 100, random_state: int | None = None) -> dict[str, Any]

Perform bootstrap validation.

Parameters:

Name	Type	Description	Default
`model_terms`	`list of str`	Model terms to validate.	required
`n_bootstrap`	`int`	Number of bootstrap samples, by default 100.	`100`
`random_state`	`int`	Seed for reproducible resampling.	`None`

Returns:

Type	Description
`dict`	Bootstrap validation results.

Source code in src/industrialstats/analysis/model_fitting.py

def bootstrap_validation(
    self,
    model_terms: list[str],
    n_bootstrap: int = 100,
    random_state: int | None = None,
) -> dict[str, Any]:
    """Perform bootstrap validation.

    Parameters
    ----------
    model_terms : list of str
        Model terms to validate.
    n_bootstrap : int, optional
        Number of bootstrap samples, by default 100.
    random_state : int, optional
        Seed for reproducible resampling.

    Returns
    -------
    dict
        Bootstrap validation results.
    """
    rng = np.random.default_rng(random_state)
    n_samples = len(self.data)
    bootstrap_results = {
        "coefficients": {term: [] for term in model_terms},
        "r_squared": [],
        "rmse": [],
        "predictions": [],
    }

    # Original model for comparison
    original_model = self._fit_terms(model_terms)

    for _bootstrap_idx in range(n_bootstrap):
        # Create bootstrap sample
        bootstrap_indices = rng.integers(0, n_samples, size=n_samples)
        bootstrap_data = self.data.iloc[bootstrap_indices].reset_index(drop=True)

        try:
            # Fit model on bootstrap sample
            bootstrap_fitter = ModelFitting(bootstrap_data, self.response)
            bootstrap_model = bootstrap_fitter._fit_terms(model_terms)

            # Store coefficients
            for term in model_terms:
                if term in bootstrap_model["coefficients"]:
                    bootstrap_results["coefficients"][term].append(
                        bootstrap_model["coefficients"][term]
                    )
                else:
                    bootstrap_results["coefficients"][term].append(np.nan)

            # Store model metrics
            bootstrap_results["r_squared"].append(
                bootstrap_model["model_metrics"]["R2"]
            )
            bootstrap_results["rmse"].append(
                bootstrap_model["model_metrics"]["RMSE"]
            )

            # Predict on original data
            predictions = self._predict_with_model(
                bootstrap_model, self.data, model_terms
            )
            bootstrap_results["predictions"].append(predictions)

        except (ValueError, np.linalg.LinAlgError) as e:
            logger.debug("Bootstrap iteration failed: %s", e)
            for term in model_terms:
                bootstrap_results["coefficients"][term].append(np.nan)
            bootstrap_results["r_squared"].append(np.nan)
            bootstrap_results["rmse"].append(np.nan)
            bootstrap_results["predictions"].append(np.full(n_samples, np.nan))

    # Calculate bootstrap statistics
    bootstrap_stats = {}

    # Coefficient statistics
    for term in model_terms:
        coeff_values = [
            c for c in bootstrap_results["coefficients"][term] if not np.isnan(c)
        ]
        if coeff_values:
            bootstrap_stats[f"{term}_mean"] = np.mean(coeff_values)
            bootstrap_stats[f"{term}_std"] = np.std(coeff_values)
            bootstrap_stats[f"{term}_ci_lower"] = np.percentile(coeff_values, 2.5)
            bootstrap_stats[f"{term}_ci_upper"] = np.percentile(coeff_values, 97.5)

            # Bias calculation
            original_coeff = original_model["coefficients"].get(term, 0)
            bootstrap_stats[f"{term}_bias"] = np.mean(coeff_values) - original_coeff

    # Model performance statistics
    valid_r2 = [r2 for r2 in bootstrap_results["r_squared"] if not np.isnan(r2)]
    valid_rmse = [rmse for rmse in bootstrap_results["rmse"] if not np.isnan(rmse)]

    if valid_r2:
        bootstrap_stats["r2_mean"] = np.mean(valid_r2)
        bootstrap_stats["r2_std"] = np.std(valid_r2)
        bootstrap_stats["r2_ci_lower"] = np.percentile(valid_r2, 2.5)
        bootstrap_stats["r2_ci_upper"] = np.percentile(valid_r2, 97.5)

    if valid_rmse:
        bootstrap_stats["rmse_mean"] = np.mean(valid_rmse)
        bootstrap_stats["rmse_std"] = np.std(valid_rmse)
        bootstrap_stats["rmse_ci_lower"] = np.percentile(valid_rmse, 2.5)
        bootstrap_stats["rmse_ci_upper"] = np.percentile(valid_rmse, 97.5)

    return {
        "bootstrap_results": bootstrap_results,
        "bootstrap_stats": bootstrap_stats,
        "original_model": original_model,
        "n_bootstrap": n_bootstrap,
        "success_rate": len(valid_r2) / n_bootstrap,
    }

regularized_fitting ¶

regularized_fitting(method: str = 'lasso', alphas: Sequence[float] | None = None, cv: int = 5, l1_ratio: float = 0.5, plot_path: bool = False, random_state: int | None = None) -> dict[str, Any]

Fit linear models with regularization and cross-validation.

Parameters:

Name	Type	Description	Default
`method`	`(lasso, ridge, elasticnet)`	Regularization technique to use. Default is `"lasso"`.	`"lasso"`
`alphas`	`sequence of float`	Grid of regularization strengths to evaluate. If `None`, scikit-learn chooses an appropriate set.	`None`
`cv`	`int`	Number of cross-validation folds. Default is `5`.	`5`
`l1_ratio`	`float`	Elastic net mixing parameter, with `0` = ridge and `1` = lasso. Used only when `method="elasticnet"`. Default is `0.5`.	`0.5`
`plot_path`	`bool`	If `True`, plot coefficient paths across regularization strengths.	`False`
`random_state`	`int`	Seed for reproducible cross-validation splits.	`None`

Returns:

Type	Description
`dict`	Results containing fitted model, selected features, and path data.

References

.. [1] Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B. .. [2] Zou, H., & Hastie, T. (2005). Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society: Series B.

Source code in src/industrialstats/analysis/model_fitting.py

def regularized_fitting(
    self,
    method: str = "lasso",
    alphas: Sequence[float] | None = None,
    cv: int = 5,
    l1_ratio: float = 0.5,
    plot_path: bool = False,
    random_state: int | None = None,
) -> dict[str, Any]:
    """Fit linear models with regularization and cross-validation.

    Parameters
    ----------
    method : {"lasso", "ridge", "elasticnet"}, optional
        Regularization technique to use. Default is ``"lasso"``.
    alphas : sequence of float, optional
        Grid of regularization strengths to evaluate. If ``None``,
        scikit-learn chooses an appropriate set.
    cv : int, optional
        Number of cross-validation folds. Default is ``5``.
    l1_ratio : float, optional
        Elastic net mixing parameter, with ``0`` = ridge and ``1`` = lasso.
        Used only when ``method="elasticnet"``. Default is ``0.5``.
    plot_path : bool, optional
        If ``True``, plot coefficient paths across regularization strengths.
    random_state : int, optional
        Seed for reproducible cross-validation splits.

    Returns
    -------
    dict
        Results containing fitted model, selected features, and path data.

    References
    ----------
    .. [1] Tibshirani, R. (1996). Regression shrinkage and selection via the
           lasso. Journal of the Royal Statistical Society: Series B.
    .. [2] Zou, H., & Hastie, T. (2005). Regularization and variable
           selection via the elastic net. Journal of the Royal Statistical
           Society: Series B.
    """

    from sklearn.linear_model import (
        ElasticNetCV,
        LassoCV,
        Ridge,
        RidgeCV,
        enet_path,
        lasso_path,
    )

    X = pd.get_dummies(self.data[self.factor_columns], drop_first=True)
    y = self.data[self.response].values

    # scikit-learn deprecated passing ``alphas=None`` explicitly: from 1.9
    # the default becomes an alpha count rather than None. Omitting the
    # argument entirely selects the library default on every supported
    # version and keeps the automatic grid behaviour.
    alpha_kwargs = {} if alphas is None else {"alphas": alphas}

    if method.lower() == "lasso":
        model = LassoCV(cv=cv, random_state=random_state, **alpha_kwargs).fit(X, y)
        path_alphas, coefs, _ = lasso_path(X, y, **alpha_kwargs)
    elif method.lower() == "ridge":
        if alphas is None:
            alphas = np.logspace(-6, 6, 100)
        model = RidgeCV(alphas=alphas, cv=cv).fit(X, y)
        path_alphas = np.array(alphas)
        coefs = []
        for a in path_alphas:
            coefs.append(Ridge(alpha=a).fit(X, y).coef_)
        coefs = np.array(coefs).T
    elif method.lower() == "elasticnet":
        model = ElasticNetCV(
            l1_ratio=l1_ratio,
            cv=cv,
            random_state=random_state,
            **alpha_kwargs,
        ).fit(X, y)
        path_alphas, coefs, _ = enet_path(X, y, l1_ratio=l1_ratio, **alpha_kwargs)
    else:
        raise ValueError("method must be 'lasso', 'ridge', or 'elasticnet'")

    coefficients = dict(zip(X.columns, model.coef_, strict=True))
    selected_features = [
        feat for feat, coef in coefficients.items() if not np.isclose(coef, 0.0)
    ]

    if plot_path:
        for idx, feat in enumerate(X.columns):
            plt.plot(path_alphas, coefs[idx], label=feat)
        plt.xscale("log")
        plt.gca().invert_xaxis()
        plt.xlabel("alpha")
        plt.ylabel("coefficient")
        plt.title(f"{method.title()} coefficient paths")
        plt.legend(loc="best")
        plt.tight_layout()

    return {
        "model": model,
        "best_alpha": getattr(model, "alpha_", None),
        "coefficients": coefficients,
        "selected_features": selected_features,
        "path_alphas": path_alphas,
        "path_coefficients": coefs,
        "method": method.lower(),
    }

model_comparison ¶

model_comparison(model_list: list[list[str]]) -> DataFrame

Compare multiple models using various criteria.

Parameters:

Name	Type	Description	Default
`model_list`	`list of list of str`	List of model term lists to compare.	required

Returns:

Type	Description
`DataFrame`	Model comparison table.

Source code in src/industrialstats/analysis/model_fitting.py

def model_comparison(self, model_list: list[list[str]]) -> pd.DataFrame:
    """Compare multiple models using various criteria.

    Parameters
    ----------
    model_list : list of list of str
        List of model term lists to compare.

    Returns
    -------
    pandas.DataFrame
        Model comparison table.
    """
    comparison_results = []

    for i, model_terms in enumerate(model_list):
        try:
            model_result = self._fit_terms(model_terms)

            comparison_results.append(
                {
                    "Model": f"Model_{i + 1}",
                    "Terms": " + ".join(model_terms),
                    "N_Terms": len(model_terms),
                    "R2": model_result["model_metrics"]["R2"],
                    "R2_Adj": model_result["model_metrics"]["R2_adj"],
                    "AIC": model_result["model_metrics"]["AIC"],
                    "BIC": model_result["model_metrics"]["BIC"],
                    "RMSE": model_result["model_metrics"]["RMSE"],
                    "F_Statistic": model_result["model_metrics"].get(
                        "F_statistic", np.nan
                    ),
                    "F_P_Value": model_result["model_metrics"].get(
                        "F_p_value", np.nan
                    ),
                }
            )

            # Store fitted model
            self.fitted_models[f"Model_{i + 1}"] = model_result

        except (ValueError, np.linalg.LinAlgError) as e:
            logger.debug("Model comparison failed for model %s: %s", i + 1, e)
            comparison_results.append(
                {
                    "Model": f"Model_{i + 1}",
                    "Terms": " + ".join(model_terms),
                    "N_Terms": len(model_terms),
                    "Error": str(e),
                }
            )

    self.model_comparison = pd.DataFrame(comparison_results)
    return self.model_comparison

residual_diagnostics ¶

residual_diagnostics(model_terms: list[str]) -> dict[str, Any]

Perform comprehensive residual diagnostics.

Parameters:

Name	Type	Description	Default
`model_terms`	`List[str]`	Model terms to diagnose.	required

Returns:

Type	Description
`Dict[str, Any]`	Diagnostic results.

Source code in src/industrialstats/analysis/model_fitting.py

def residual_diagnostics(self, model_terms: list[str]) -> dict[str, Any]:
    """Perform comprehensive residual diagnostics.

    Parameters
    ----------
    model_terms : List[str]
        Model terms to diagnose.

    Returns
    -------
    Dict[str, Any]
        Diagnostic results.
    """
    model_result = self._fit_terms(model_terms)

    residuals = model_result["residuals"]

    diagnostics = {}

    # Basic residual statistics
    diagnostics["residual_stats"] = {
        "mean": np.mean(residuals),
        "std": np.std(residuals),
        "min": np.min(residuals),
        "max": np.max(residuals),
        "range": np.max(residuals) - np.min(residuals),
    }

    # Normality tests
    try:
        shapiro_stat, shapiro_p = stats.shapiro(residuals)
        diagnostics["normality_test"] = {
            "shapiro_wilk_statistic": shapiro_stat,
            "shapiro_wilk_p_value": shapiro_p,
            "normal_assumption": shapiro_p > 0.05,
        }
    except (ValueError, np.linalg.LinAlgError) as e:
        logger.debug("Normality test failed: %s", e)
        diagnostics["normality_test"] = {
            "error": f"Unable to perform normality test: {e}"
        }

    # Homoscedasticity tests
    try:
        # Breusch-Pagan test
        from statsmodels.stats.diagnostic import het_breuschpagan

        bp_stat, bp_p, _bp_f_stat, _bp_f_p = het_breuschpagan(
            residuals, model_result["model_object"].model.exog
        )

        diagnostics["homoscedasticity_test"] = {
            "breusch_pagan_statistic": bp_stat,
            "breusch_pagan_p_value": bp_p,
            "homoscedastic_assumption": bp_p > 0.05,
        }
    except (ValueError, np.linalg.LinAlgError, ImportError) as e:
        logger.debug("Homoscedasticity test failed: %s", e)
        diagnostics["homoscedasticity_test"] = {
            "error": f"Unable to perform homoscedasticity test: {e}"
        }

    # Independence test
    try:
        from statsmodels.stats.diagnostic import durbin_watson

        dw_stat = durbin_watson(residuals)

        diagnostics["independence_test"] = {
            "durbin_watson_statistic": dw_stat,
            "independent_assumption": 1.5 <= dw_stat <= 2.5,
        }
    except (ValueError, ImportError) as e:
        logger.debug("Independence test failed: %s", e)
        diagnostics["independence_test"] = {
            "error": f"Unable to perform independence test: {e}"
        }

    # Outlier detection
    standardized_residuals = residuals / np.std(residuals)
    outliers = np.abs(standardized_residuals) > 2.5

    diagnostics["outlier_analysis"] = {
        "n_outliers": np.sum(outliers),
        "outlier_indices": np.where(outliers)[0].tolist(),
        "max_standardized_residual": np.max(np.abs(standardized_residuals)),
    }

    # Leverage and influence
    try:
        influence = model_result["model_object"].get_influence()
        leverage = influence.hat_matrix_diag
        cooks_d = influence.cooks_distance[0]

        high_leverage = leverage > 2 * len(model_terms) / len(self.data)
        high_influence = cooks_d > 4 / len(self.data)

        diagnostics["leverage_influence"] = {
            "max_leverage": np.max(leverage),
            "n_high_leverage": np.sum(high_leverage),
            "high_leverage_indices": np.where(high_leverage)[0].tolist(),
            "max_cooks_d": np.max(cooks_d),
            "n_high_influence": np.sum(high_influence),
            "high_influence_indices": np.where(high_influence)[0].tolist(),
        }
    except (ValueError, np.linalg.LinAlgError) as e:
        logger.debug("Leverage and influence calculation failed: %s", e)
        diagnostics["leverage_influence"] = {
            "error": f"Unable to calculate leverage and influence: {e}"
        }

    return diagnostics

lack_of_fit_test ¶

lack_of_fit_test(model_terms: list[str]) -> dict[str, Any]

Perform lack-of-fit test for models with replicates.

Parameters:

Name	Type	Description	Default
`model_terms`	`List[str]`	Model terms to test.	required

Returns:

Type	Description
`Dict[str, Any]`	Lack-of-fit test results.

Source code in src/industrialstats/analysis/model_fitting.py

def lack_of_fit_test(self, model_terms: list[str]) -> dict[str, Any]:
    """Perform lack-of-fit test for models with replicates.

    Parameters
    ----------
    model_terms : List[str]
        Model terms to test.

    Returns
    -------
    Dict[str, Any]
        Lack-of-fit test results.
    """
    # Check if we have replicates
    factor_combinations = self.data[self.factor_columns].drop_duplicates()

    if len(factor_combinations) == len(self.data):
        return {"error": "No replicates found for lack-of-fit test"}

    # Fit the model
    model_result = self._fit_terms(model_terms)

    # Calculate pure error and lack-of-fit
    pure_error_ss = 0
    pure_error_df = 0

    for _, combination in factor_combinations.iterrows():
        # Find all replicates for this combination
        mask = True
        for factor in self.factor_columns:
            mask &= self.data[factor] == combination[factor]

        replicates = self.data[mask]

        if len(replicates) > 1:
            # Calculate pure error for this combination
            replicate_responses = replicates[self.response].values
            replicate_mean = np.mean(replicate_responses)

            pure_error_ss += np.sum((replicate_responses - replicate_mean) ** 2)
            pure_error_df += len(replicates) - 1

    if pure_error_df == 0:
        return {"error": "Insufficient replicates for lack-of-fit test"}

    # Calculate lack-of-fit
    total_error_ss = np.sum(model_result["residuals"] ** 2)
    total_error_df = len(self.data) - len(model_terms)

    lof_ss = total_error_ss - pure_error_ss
    lof_df = total_error_df - pure_error_df

    if lof_df <= 0:
        return {"error": "Model is saturated - cannot test lack-of-fit"}

    # Calculate F-statistic
    lof_ms = lof_ss / lof_df
    pure_error_ms = pure_error_ss / pure_error_df

    f_statistic = lof_ms / pure_error_ms
    p_value = 1 - stats.f.cdf(f_statistic, lof_df, pure_error_df)

    return {
        "lack_of_fit_ss": lof_ss,
        "lack_of_fit_df": lof_df,
        "lack_of_fit_ms": lof_ms,
        "pure_error_ss": pure_error_ss,
        "pure_error_df": pure_error_df,
        "pure_error_ms": pure_error_ms,
        "f_statistic": f_statistic,
        "p_value": p_value,
        "adequate_fit": p_value > 0.05,
        "n_unique_combinations": len(factor_combinations),
        "total_observations": len(self.data),
    }

Diagnostics¶

Model diagnostics are a decision process, not a single goodness-of-fit number. The package combines formal assumption checks, influence diagnostics, explicit outlier thresholds, residual plots, and an adequacy summary before producing remediation guidance.

industrialstats model diagnostics and remediation flow

Source: ../diagrams/model_diagnostics.dot

Formal assumption tests should be interpreted together with residual plots and influence measures. A model can pass a normality or variance test and still contain observations with enough leverage or Cook's distance to destabilize inference; conversely, a flagged point should be investigated rather than deleted automatically.

industrialstats.analysis.diagnostics ¶

Comprehensive regression diagnostics.

This module implements residual diagnostics described by Cook & Weisberg (1982) and extends classical assumption checks with actionable guidance. The :class:ModelDiagnostics class operates on dictionaries returned by the high-level fitting utilities in :mod:industrialstats.analysis.model_fitting and requires the original design data to contextualize the diagnostics.

Examples:

>>> import pandas as pd
>>> import numpy as np
>>> import statsmodels.api as sm
>>> from industrialstats.analysis.diagnostics import ModelDiagnostics
>>> rng = np.random.default_rng(42)
>>> x1 = rng.normal(size=120)
>>> x2 = rng.normal(size=120)
>>> y = 1.5 + 2.0 * x1 - 1.2 * x2 + rng.normal(size=120)
>>> data = pd.DataFrame({"y": y, "x1": x1, "x2": x2})
>>> model = sm.OLS(data["y"], sm.add_constant(data[["x1", "x2"]])).fit()
>>> model_result = {
...     "model_object": model,
...     "residuals": model.resid,
...     "fitted_values": model.fittedvalues,
...     "model_metrics": {"R2": model.rsquared},
... }
>>> diagnostics = ModelDiagnostics(model_result, data)
>>> summary = diagnostics.assumption_tests()
>>> round(summary["normality"]["shapiro"]["p_value"], 3) >= 0.05
True

ModelDiagnostics ¶

ModelDiagnostics(model_result: dict[str, Any], data: DataFrame)

Diagnostic analytics for linear models.

The diagnostics follow the influence framework of Cook & Weisberg [1]_ and are compatible with dictionaries returned by :class:~industrialstats.analysis.model_fitting.ModelFitting.

Parameters:

Name	Type	Description	Default
`model_result`	`Dict[str, Any]`	Dictionary containing the fitted `statsmodels` object under the key `"model_object"` alongside residual vectors and fitted values. The minimal expected keys are `{"model_object", "residuals", "fitted_values"}`.	required
`data`	`DataFrame`	Original dataset used during model fitting. The frame is copied to avoid inadvertent mutation during diagnostics.	required

Raises:

Type	Description
`TypeError`	If `model_result` is not a dictionary or lacks the required `statsmodels` interfaces.
`KeyError`	When mandatory keys are absent from `model_result`.
`ValueError`	If residual and fitted vector lengths are inconsistent with `data`.

References

.. [1] Cook, R. D., & Weisberg, S. (1982). Residuals and Influence in Regression. Chapman & Hall/CRC.

Source code in src/industrialstats/analysis/diagnostics.py

def __init__(self, model_result: dict[str, Any], data: pd.DataFrame) -> None:
    if not isinstance(model_result, dict):
        raise TypeError("model_result must be a dictionary of model outputs")

    required_keys = {"model_object", "residuals", "fitted_values"}
    missing_keys = required_keys.difference(model_result.keys())
    if missing_keys:
        raise KeyError(
            "model_result is missing required keys: "
            + ", ".join(sorted(missing_keys))
        )

    model_object = model_result["model_object"]
    if not hasattr(model_object, "get_influence"):
        raise TypeError(
            "model_result['model_object'] must expose statsmodels influence diagnostics"
        )

    self.model_result = model_result
    self.model = model_object
    self.data = data.copy(deep=True)

    self.residuals = np.asarray(model_result["residuals"], dtype=float)
    self.fitted_values = np.asarray(model_result["fitted_values"], dtype=float)

    if self.residuals.ndim != 1 or self.fitted_values.ndim != 1:
        raise ValueError(
            "residuals and fitted_values must be one-dimensional arrays"
        )

    n_obs = len(self.data)
    if (
        len(self.residuals) != len(self.fitted_values)
        or len(self.residuals) != n_obs
    ):
        raise ValueError(
            "Length of residuals, fitted values, and data rows must match"
        )

    self._assumption_cache: dict[str, dict[str, Any]] | None = None
    self._influence_cache: dict[str, np.ndarray] | None = None
    self._outlier_cache: dict[str, list[int]] | None = None

assumption_tests ¶

assumption_tests() -> dict[str, dict[str, Any]]

Evaluate classical regression assumptions.

The procedure combines the Shapiro-Wilk and Anderson-Darling tests for normality, Levene and Bartlett tests for homoscedasticity across fitted quantile groups, and the Durbin-Watson statistic for independence.

Returns:

Type	Description
`Dict[str, Dict[str, Any]]`	Nested mapping summarising each assumption. For example, `result["normality"]["passes"]` indicates whether both normality tests are satisfied at the 5% level. The Anderson-Darling entry reports whichever evidence the installed SciPy can supply: `"p_value"` on SciPy 1.17 and newer, or `"critical_value_5pct"` on older releases. Both keys are always present and the unavailable one is `None`.

Examples:

>>> tests = diagnostics.assumption_tests()
>>> sorted(tests.keys())
['homoscedasticity', 'independence', 'normality']
>>> tests["independence"]["passes"]
True

Source code in src/industrialstats/analysis/diagnostics.py

def assumption_tests(self) -> dict[str, dict[str, Any]]:
    """Evaluate classical regression assumptions.

    The procedure combines the Shapiro-Wilk and Anderson-Darling tests for
    normality, Levene and Bartlett tests for homoscedasticity across fitted
    quantile groups, and the Durbin-Watson statistic for independence.

    Returns
    -------
    Dict[str, Dict[str, Any]]
        Nested mapping summarising each assumption. For example,
        ``result["normality"]["passes"]`` indicates whether both normality
        tests are satisfied at the 5% level.

        The Anderson-Darling entry reports whichever evidence the installed
        SciPy can supply: ``"p_value"`` on SciPy 1.17 and newer, or
        ``"critical_value_5pct"`` on older releases. Both keys are always
        present and the unavailable one is ``None``.

    Examples
    --------
    >>> tests = diagnostics.assumption_tests()
    >>> sorted(tests.keys())
    ['homoscedasticity', 'independence', 'normality']
    >>> tests["independence"]["passes"]
    True
    """

    if self._assumption_cache is not None:
        return self._assumption_cache

    if len(self.residuals) < 8:
        raise ValueError(
            "At least 8 observations are required for assumption tests"
        )

    # Normality diagnostics
    shapiro_stat, shapiro_p = stats.shapiro(self.residuals)
    anderson_stat, anderson_p, anderson_crit, anderson_pass = (
        self._anderson_normality(self.residuals)
    )
    normality_pass = (shapiro_p > 0.05) and anderson_pass

    # Homoscedasticity diagnostics via fitted quantile groups
    df = pd.DataFrame({"fitted": self.fitted_values, "resid": self.residuals})
    try:
        df["group"] = pd.qcut(
            df["fitted"],
            q=min(4, max(2, df["fitted"].nunique())),
            duplicates="drop",
        )
    except ValueError as exc:  # occurs when data are constant
        raise ValueError("Cannot form groups for homoscedasticity checks") from exc

    grouped = [
        grp["resid"].to_numpy() for _, grp in df.groupby("group", observed=True)
    ]
    if len(grouped) < 2:
        raise ValueError("Need at least two groups to run homoscedasticity tests")

    levene_stat, levene_p = stats.levene(*grouped, center="median")
    bartlett_stat, bartlett_p = stats.bartlett(*grouped)
    homoscedastic_pass = (levene_p > 0.05) and (bartlett_p > 0.05)

    # Independence diagnostics via Durbin-Watson
    dw_stat = float(durbin_watson(self.residuals))
    independence_pass = 1.5 <= dw_stat <= 2.5

    self._assumption_cache = {
        "normality": {
            "passes": bool(normality_pass),
            "shapiro": {
                "statistic": float(shapiro_stat),
                "p_value": float(shapiro_p),
            },
            "anderson": {
                "statistic": anderson_stat,
                "p_value": anderson_p,
                "critical_value_5pct": anderson_crit,
            },
        },
        "homoscedasticity": {
            "passes": bool(homoscedastic_pass),
            "levene": {
                "statistic": float(levene_stat),
                "p_value": float(levene_p),
            },
            "bartlett": {
                "statistic": float(bartlett_stat),
                "p_value": float(bartlett_p),
            },
        },
        "independence": {
            "passes": bool(independence_pass),
            "durbin_watson": float(dw_stat),
        },
    }
    return self._assumption_cache

influence_analysis ¶

influence_analysis() -> dict[str, ndarray]

Compute influence diagnostics under the Cook & Weisberg framework.

Returns:

Type	Description
`Dict[str, ndarray]`	Arrays of studentized residuals, Cook's distances, DFFITS, leverage, and DFBETAS for each observation.

Examples:

>>> influence = diagnostics.influence_analysis()
>>> {k: v.shape for k, v in influence.items()}["leverage"]
(120,)

Source code in src/industrialstats/analysis/diagnostics.py

def influence_analysis(self) -> dict[str, np.ndarray]:
    """Compute influence diagnostics under the Cook & Weisberg framework.

    Returns
    -------
    Dict[str, numpy.ndarray]
        Arrays of studentized residuals, Cook's distances, DFFITS, leverage,
        and DFBETAS for each observation.

    Examples
    --------
    >>> influence = diagnostics.influence_analysis()
    >>> {k: v.shape for k, v in influence.items()}["leverage"]
    (120,)
    """

    if self._influence_cache is not None:
        return self._influence_cache

    influence = self.model.get_influence()
    studentized = influence.resid_studentized_external
    cooks_d = influence.cooks_distance[0]
    dffits = influence.dffits[0]
    leverage = influence.hat_matrix_diag
    dfbetas = influence.dfbetas

    self._influence_cache = {
        "studentized_residuals": np.asarray(studentized, dtype=float),
        "cooks_distance": np.asarray(cooks_d, dtype=float),
        "dffits": np.asarray(dffits, dtype=float),
        "leverage": np.asarray(leverage, dtype=float),
        "dfbetas": np.asarray(dfbetas, dtype=float),
    }
    return self._influence_cache

outlier_detection ¶

outlier_detection() -> dict[str, list[int]]

Identify influential observations with multiple criteria.

Returns:

Type	Description
`Dict[str, List[int]]`	Observation indices flagged by studentized residual, Cook's distance and DFBETAS thresholds. Indices are returned in ascending order.

Examples:

>>> diagnostics.outlier_detection()["cooks_distance"]
[]

Source code in src/industrialstats/analysis/diagnostics.py

def outlier_detection(self) -> dict[str, list[int]]:
    """Identify influential observations with multiple criteria.

    Returns
    -------
    Dict[str, List[int]]
        Observation indices flagged by studentized residual, Cook's distance
        and DFBETAS thresholds. Indices are returned in ascending order.

    Examples
    --------
    >>> diagnostics.outlier_detection()["cooks_distance"]
    []
    """

    if self._outlier_cache is not None:
        return self._outlier_cache

    influence = self.influence_analysis()
    n_obs = len(self.residuals)
    p_params = int(getattr(self.model, "df_model", len(self.model.params) - 1)) + 1

    studentized = influence["studentized_residuals"]
    cooks_d = influence["cooks_distance"]
    dfbetas = influence["dfbetas"]

    student_threshold = 3.0
    cook_threshold = 4.0 / max(n_obs, 1)
    dfbetas_threshold = 2.0 / math.sqrt(max(n_obs, 1))

    student_idx = np.where(np.abs(studentized) > student_threshold)[0]
    cooks_idx = np.where(cooks_d > cook_threshold)[0]
    dfbetas_idx = np.unique(np.where(np.abs(dfbetas) > dfbetas_threshold)[0])

    leverage = influence["leverage"]
    leverage_threshold = 2.0 * p_params / max(n_obs, 1)
    leverage_idx = np.where(leverage > leverage_threshold)[0]

    self._outlier_cache = {
        "studentized_residuals": student_idx.astype(int).tolist(),
        "cooks_distance": cooks_idx.astype(int).tolist(),
        "dfbetas": dfbetas_idx.astype(int).tolist(),
        "leverage": leverage_idx.astype(int).tolist(),
    }
    return self._outlier_cache

model_adequacy ¶

model_adequacy() -> dict[str, Any]

Summarise overall adequacy of the fitted model.

The summary merges assumption test outcomes, influence diagnostics, and model fit statistics. Diagnostic plots for residual behaviour and Cook's distance are included to support expert review.

Returns:

Type	Description
`Dict[str, Any]`	Dictionary with keys `assumptions`, `outliers`, `influence`, `model_metrics`, `overall_pass`, and `plots`.

Examples:

>>> adequacy = diagnostics.model_adequacy()
>>> sorted(adequacy["plots"].keys())
['cook_distance', 'qq_plot', 'residuals_vs_fitted']

Source code in src/industrialstats/analysis/diagnostics.py

def model_adequacy(self) -> dict[str, Any]:
    """Summarise overall adequacy of the fitted model.

    The summary merges assumption test outcomes, influence diagnostics, and
    model fit statistics. Diagnostic plots for residual behaviour and Cook's
    distance are included to support expert review.

    Returns
    -------
    Dict[str, Any]
        Dictionary with keys ``assumptions``, ``outliers``, ``influence``,
        ``model_metrics``, ``overall_pass``, and ``plots``.

    Examples
    --------
    >>> adequacy = diagnostics.model_adequacy()
    >>> sorted(adequacy["plots"].keys())
    ['cook_distance', 'qq_plot', 'residuals_vs_fitted']
    """

    assumptions = self.assumption_tests()
    outliers = self.outlier_detection()
    influence = self.influence_analysis()
    metrics = self.model_result.get("model_metrics", {})

    assumption_pass = all(result["passes"] for result in assumptions.values())
    outlier_count = sum(len(indices) for indices in outliers.values())
    influence_summary = {
        key: float(np.nanmax(np.abs(values))) for key, values in influence.items()
    }

    try:
        plots = self._generate_diagnostic_plots(influence)
    except RuntimeError:
        plots = {}

    overall_pass = (
        assumption_pass
        and outlier_count == 0
        and influence_summary["cooks_distance"] < (4.0 / len(self.residuals))
    )

    return {
        "assumptions": assumptions,
        "outliers": outliers,
        "influence": influence_summary,
        "model_metrics": metrics,
        "overall_pass": bool(overall_pass),
        "plots": plots,
    }

recommendation_system ¶

recommendation_system() -> list[str]

Produce actionable recommendations based on diagnostics.

Recommendations interpret assumption violations and influential point detections to guide remedial strategies such as variance-stabilising transformations, robust regression, or data review.

Returns:

Type	Description
`List[str]`	Human-readable recommendations ordered by severity.

Examples:

>>> diagnostics.recommendation_system()
['No major issues detected. Consider validating on a holdout set.']

Source code in src/industrialstats/analysis/diagnostics.py

def recommendation_system(self) -> list[str]:
    """Produce actionable recommendations based on diagnostics.

    Recommendations interpret assumption violations and influential point
    detections to guide remedial strategies such as variance-stabilising
    transformations, robust regression, or data review.

    Returns
    -------
    List[str]
        Human-readable recommendations ordered by severity.

    Examples
    --------
    >>> diagnostics.recommendation_system()  # doctest: +SKIP
    ['No major issues detected. Consider validating on a holdout set.']
    """

    suggestions: list[str] = []
    assumptions = self.assumption_tests()
    outliers = self.outlier_detection()
    influence = self.influence_analysis()

    if not assumptions["normality"]["passes"]:
        suggestions.append(
            "Residuals deviate from normality; consider Box-Cox transformations or non-parametric approaches."
        )
    if not assumptions["homoscedasticity"]["passes"]:
        suggestions.append(
            "Variance heterogeneity detected; weighted least squares or modelling variance as a function of predictors is recommended."
        )
    if not assumptions["independence"]["passes"]:
        suggestions.append(
            "Residual autocorrelation present; incorporate lag terms or mixed-effects structures to address dependence."
        )

    if any(outliers.values()):
        suggestions.append(
            "Investigate high-influence observations flagged by studentized residuals, Cook's distance, or DFBETAS before finalising conclusions."
        )

    cooks_peak = float(np.nanmax(np.abs(influence["cooks_distance"])))
    if cooks_peak > (4.0 / len(self.residuals)):
        suggestions.append(
            "Cook's distance exceeds the 4/n heuristic; reassess the modelling assumptions for the flagged runs."
        )

    if not suggestions:
        suggestions.append(
            "No major issues detected. Consider validating on a holdout set to confirm predictive adequacy."
        )

    return suggestions

Power analysis¶

industrialstats.analysis.power_analysis ¶

Power analysis and sample size determination for experimental designs.

PowerAnalysisResult `dataclass` ¶

PowerAnalysisResult(effect_size: float, alpha: float, power: float, sample_size: int, test_type: str, additional_info: dict[str, Any])

Container for power analysis results.

PowerAnalysis ¶

PowerAnalysis()

Comprehensive power analysis for experimental designs.

Supports power calculations for t-tests, ANOVA, factorial designs, and regression models.

Initialize power analysis.

Source code in src/industrialstats/analysis/power_analysis.py

def __init__(self):
    """Initialize power analysis."""
    self.results_history: list[PowerAnalysisResult] = []

t_test_power ¶

t_test_power(effect_size: float | None = None, alpha: float = 0.05, power: float | None = None, sample_size: int | None = None, test_type: str = 'two_sample') -> PowerAnalysisResult

Power analysis for t-tests.

Parameters:

Name	Type	Description	Default
`effect_size`	`float`	Cohen's d effect size.	`None`
`alpha`	`float`	Type I error rate. Defaults to 0.05.	`0.05`
`power`	`float`	Statistical power (`1 - eta`).	`None`
`sample_size`	`int`	Sample size per group.	`None`
`test_type`	`str`	Type of t-test (`"one_sample"`, `"two_sample"`, or `"paired"`). Defaults to `"two_sample"`.	`'two_sample'`

Returns:

Type	Description
`PowerAnalysisResult`	Power analysis results.

Source code in src/industrialstats/analysis/power_analysis.py

def t_test_power(
    self,
    effect_size: float | None = None,
    alpha: float = 0.05,
    power: float | None = None,
    sample_size: int | None = None,
    test_type: str = "two_sample",
) -> PowerAnalysisResult:
    """Power analysis for t-tests.

    Parameters
    ----------
    effect_size : float, optional
        Cohen's d effect size.
    alpha : float, optional
        Type I error rate. Defaults to 0.05.
    power : float, optional
        Statistical power (``1 - \beta``).
    sample_size : int, optional
        Sample size per group.
    test_type : str, optional
        Type of t-test (``"one_sample"``, ``"two_sample"``, or ``"paired"``).
        Defaults to ``"two_sample"``.

    Returns
    -------
    PowerAnalysisResult
        Power analysis results.
    """
    # Validate inputs
    non_none_params = sum(x is not None for x in [effect_size, power, sample_size])
    if non_none_params != 2:
        raise ValueError(
            "Exactly two of effect_size, power, sample_size must be specified"
        )

    if test_type not in ["one_sample", "two_sample", "paired"]:
        raise ValueError(
            "test_type must be 'one_sample', 'two_sample', or 'paired'"
        )

    # Calculate missing parameter
    if effect_size is None:
        effect_size = self._solve_for_effect_size_t_test(
            alpha, power, sample_size, test_type
        )
    elif power is None:
        power = self._calculate_power_t_test(
            effect_size, alpha, sample_size, test_type
        )
    elif sample_size is None:
        sample_size = self._solve_for_sample_size_t_test(
            effect_size, alpha, power, test_type
        )

    # Additional calculations
    additional_info = {
        "critical_value": stats.t.ppf(1 - alpha / 2, sample_size - 1),
        "degrees_of_freedom": (
            sample_size - 1 if test_type != "two_sample" else 2 * sample_size - 2
        ),
        "minimum_detectable_difference": effect_size,
        "test_description": self._get_test_description(test_type),
    }

    result = PowerAnalysisResult(
        effect_size=effect_size,
        alpha=alpha,
        power=power,
        sample_size=sample_size,
        test_type=f"t_test_{test_type}",
        additional_info=additional_info,
    )

    self.results_history.append(result)
    return result

anova_power ¶

anova_power(effect_size: float | None = None, alpha: float = 0.05, power: float | None = None, sample_size: int | None = None, n_groups: int = 3) -> PowerAnalysisResult

Power analysis for one-way ANOVA.

Parameters:

Name	Type	Description	Default
`effect_size`	`float`	Cohen's `f` effect size.	`None`
`alpha`	`float`	Type I error rate. Defaults to 0.05.	`0.05`
`power`	`float`	Desired power.	`None`
`sample_size`	`int`	Sample size per group.	`None`
`n_groups`	`int`	Number of groups. Defaults to 3.	`3`

Returns:

Type	Description
`PowerAnalysisResult`	Power analysis results.

Source code in src/industrialstats/analysis/power_analysis.py

def anova_power(
    self,
    effect_size: float | None = None,
    alpha: float = 0.05,
    power: float | None = None,
    sample_size: int | None = None,
    n_groups: int = 3,
) -> PowerAnalysisResult:
    """Power analysis for one-way ANOVA.

    Parameters
    ----------
    effect_size : float, optional
        Cohen's ``f`` effect size.
    alpha : float, optional
        Type I error rate. Defaults to 0.05.
    power : float, optional
        Desired power.
    sample_size : int, optional
        Sample size per group.
    n_groups : int, optional
        Number of groups. Defaults to 3.

    Returns
    -------
    PowerAnalysisResult
        Power analysis results.
    """
    # Validate inputs
    non_none_params = sum(x is not None for x in [effect_size, power, sample_size])
    if non_none_params != 2:
        raise ValueError(
            "Exactly two of effect_size, power, sample_size must be specified"
        )

    if n_groups < 2:
        raise ValueError("n_groups must be at least 2")

    # Calculate missing parameter
    if effect_size is None:
        effect_size = self._solve_for_effect_size_anova(
            alpha, power, sample_size, n_groups
        )
    elif power is None:
        power = self._calculate_power_anova(
            effect_size, alpha, sample_size, n_groups
        )
    elif sample_size is None:
        sample_size = self._solve_for_sample_size_anova(
            effect_size, alpha, power, n_groups
        )

    # Calculate additional metrics
    df_between = n_groups - 1
    df_within = n_groups * (sample_size - 1)
    total_n = n_groups * sample_size

    additional_info = {
        "n_groups": n_groups,
        "total_sample_size": total_n,
        "df_between": df_between,
        "df_within": df_within,
        "critical_f": stats.f.ppf(1 - alpha, df_between, df_within),
        "eta_squared": effect_size**2 / (1 + effect_size**2),
    }

    result = PowerAnalysisResult(
        effect_size=effect_size,
        alpha=alpha,
        power=power,
        sample_size=sample_size,
        test_type="one_way_anova",
        additional_info=additional_info,
    )

    self.results_history.append(result)
    return result

factorial_power ¶

factorial_power(effect_size: float | None = None, alpha: float = 0.05, power: float | None = None, replicates: int | None = None, factor_levels: list[int] | None = None, effect: tuple[int, ...] | None = None) -> PowerAnalysisResult

Power analysis for factorial designs.

Parameters:

Name	Type	Description	Default
`effect_size`	`float`	Cohen's `f` effect size for main effects.	`None`
`alpha`	`float`	Type I error rate. Defaults to 0.05.	`0.05`
`power`	`float`	Statistical power.	`None`
`replicates`	`int`	Number of replicates.	`None`
`factor_levels`	`list of int`	Number of levels for each factor. Defaults to `[2, 2]`.	`None`
`effect`	`tuple of int`	Indices of factors forming the effect of interest. `(0,)` specifies the main effect for the first factor, `(0, 1)` the two-way interaction between the first and second factors.	`None`

Returns:

Type	Description
`PowerAnalysisResult`	Power analysis results.

Source code in src/industrialstats/analysis/power_analysis.py

def factorial_power(
    self,
    effect_size: float | None = None,
    alpha: float = 0.05,
    power: float | None = None,
    replicates: int | None = None,
    factor_levels: list[int] | None = None,
    effect: tuple[int, ...] | None = None,
) -> PowerAnalysisResult:
    """Power analysis for factorial designs.

    Parameters
    ----------
    effect_size : float, optional
        Cohen's ``f`` effect size for main effects.
    alpha : float, optional
        Type I error rate. Defaults to 0.05.
    power : float, optional
        Statistical power.
    replicates : int, optional
        Number of replicates.
    factor_levels : list of int, optional
        Number of levels for each factor. Defaults to ``[2, 2]``.
    effect : tuple of int, optional
        Indices of factors forming the effect of interest. ``(0,)``
        specifies the main effect for the first factor, ``(0, 1)`` the
        two-way interaction between the first and second factors.

    Returns
    -------
    PowerAnalysisResult
        Power analysis results.
    """
    # Validate inputs
    non_none_params = sum(x is not None for x in [effect_size, power, replicates])
    if non_none_params != 2:
        raise ValueError(
            "Exactly two of effect_size, power, replicates must be specified"
        )

    factor_levels = list(factor_levels) if factor_levels is not None else [2, 2]

    if len(factor_levels) < 1:
        raise ValueError("At least one factor required")

    # Calculate design parameters
    n_treatment_combinations = int(np.prod(factor_levels))

    if effect is None:
        effect = (0,)

    if any(i >= len(factor_levels) or i < 0 for i in effect):
        raise ValueError("effect indices must correspond to factor_levels")

    # Calculate missing parameter
    if effect_size is None:
        effect_size = self._solve_for_effect_size_factorial(
            alpha, power, replicates, factor_levels, effect
        )
    elif power is None:
        power = self._calculate_power_factorial(
            effect_size, alpha, replicates, factor_levels, effect
        )
    elif replicates is None:
        replicates = self._solve_for_replicates_factorial(
            effect_size, alpha, power, factor_levels, effect
        )

    # Degrees of freedom for the specified effect
    df_effect = int(np.prod([factor_levels[i] - 1 for i in effect]))

    # Calculate error degrees of freedom
    df_error = n_treatment_combinations * (replicates - 1)
    total_n = n_treatment_combinations * replicates

    additional_info = {
        "factor_levels": factor_levels,
        "n_factors": len(factor_levels),
        "n_treatment_combinations": n_treatment_combinations,
        "total_sample_size": total_n,
        "df_effect": df_effect,
        "df_error": df_error,
        "design_type": f"{len(factor_levels)}-factor factorial",
        "design_notation": "x".join(map(str, factor_levels)),
        "effect": effect,
    }

    result = PowerAnalysisResult(
        effect_size=effect_size,
        alpha=alpha,
        power=power,
        sample_size=replicates,
        test_type="factorial_design",
        additional_info=additional_info,
    )

    self.results_history.append(result)
    return result

regression_power ¶

regression_power(effect_size: float | None = None, alpha: float = 0.05, power: float | None = None, sample_size: int | None = None, n_predictors: int = 1) -> PowerAnalysisResult

Power analysis for multiple regression.

Parameters:

Name	Type	Description	Default
`effect_size`	`float`	Cohen's :math:`f^2` effect size.	`None`
`alpha`	`float`	Type I error rate. Defaults to 0.05.	`0.05`
`power`	`float`	Statistical power.	`None`
`sample_size`	`int`	Total sample size.	`None`
`n_predictors`	`int`	Number of predictors in the model. Defaults to 1.	`1`

Returns:

Type	Description
`PowerAnalysisResult`	Power analysis results.

Source code in src/industrialstats/analysis/power_analysis.py

def regression_power(
    self,
    effect_size: float | None = None,
    alpha: float = 0.05,
    power: float | None = None,
    sample_size: int | None = None,
    n_predictors: int = 1,
) -> PowerAnalysisResult:
    """Power analysis for multiple regression.

    Parameters
    ----------
    effect_size : float, optional
        Cohen's :math:`f^2` effect size.
    alpha : float, optional
        Type I error rate. Defaults to 0.05.
    power : float, optional
        Statistical power.
    sample_size : int, optional
        Total sample size.
    n_predictors : int, optional
        Number of predictors in the model. Defaults to 1.

    Returns
    -------
    PowerAnalysisResult
        Power analysis results.
    """
    # Validate inputs
    non_none_params = sum(x is not None for x in [effect_size, power, sample_size])
    if non_none_params != 2:
        raise ValueError(
            "Exactly two of effect_size, power, sample_size must be specified"
        )

    if n_predictors < 1:
        raise ValueError("n_predictors must be at least 1")

    # Calculate missing parameter
    if effect_size is None:
        effect_size = self._solve_for_effect_size_regression(
            alpha, power, sample_size, n_predictors
        )
    elif power is None:
        power = self._calculate_power_regression(
            effect_size, alpha, sample_size, n_predictors
        )
    elif sample_size is None:
        sample_size = self._solve_for_sample_size_regression(
            effect_size, alpha, power, n_predictors
        )

    # Calculate additional metrics
    df_model = n_predictors
    df_error = sample_size - n_predictors - 1

    if df_error <= 0:
        raise ValueError("Sample size too small for the number of predictors")

    r_squared = effect_size / (1 + effect_size)

    additional_info = {
        "n_predictors": n_predictors,
        "df_model": df_model,
        "df_error": df_error,
        "r_squared": r_squared,
        "adjusted_r_squared": 1 - (1 - r_squared) * (sample_size - 1) / df_error,
        "critical_f": stats.f.ppf(1 - alpha, df_model, df_error),
    }

    result = PowerAnalysisResult(
        effect_size=effect_size,
        alpha=alpha,
        power=power,
        sample_size=sample_size,
        test_type="multiple_regression",
        additional_info=additional_info,
    )

    self.results_history.append(result)
    return result

power_curve ¶

power_curve(test_type: str, fixed_params: dict[str, Any], varying_param: str, param_range: list[float]) -> dict[str, Any]

Generate power curve by varying one parameter.

Parameters:

Name	Type	Description	Default
`test_type`	`str`	Type of test (`"t_test"`, `"anova"`, `"factorial"`, `"regression"`).	required
`fixed_params`	`dict`	Fixed parameters for the analysis.	required
`varying_param`	`str`	Parameter to vary (`"effect_size"`, `"sample_size"`, `"alpha"`, `"power"`).	required
`param_range`	`list of float`	Range of values for the varying parameter.	required

Returns:

Type	Description
`Dict[str, Any]`	Power curve data and the generated plot.

Source code in src/industrialstats/analysis/power_analysis.py

def power_curve(
    self,
    test_type: str,
    fixed_params: dict[str, Any],
    varying_param: str,
    param_range: list[float],
) -> dict[str, Any]:
    """Generate power curve by varying one parameter.

    Parameters
    ----------
    test_type : str
        Type of test (``"t_test"``, ``"anova"``, ``"factorial"``, ``"regression"``).
    fixed_params : dict
        Fixed parameters for the analysis.
    varying_param : str
        Parameter to vary (``"effect_size"``, ``"sample_size"``, ``"alpha"``, ``"power"``).
    param_range : list of float
        Range of values for the varying parameter.

    Returns
    -------
    Dict[str, Any]
        Power curve data and the generated plot.
    """
    results = []

    for param_value in param_range:
        # Set up parameters
        params = fixed_params.copy()
        params[varying_param] = param_value

        try:
            # Call appropriate power analysis method
            if test_type == "t_test":
                result = self.t_test_power(**params)
            elif test_type == "anova":
                result = self.anova_power(**params)
            elif test_type == "factorial":
                result = self.factorial_power(**params)
            elif test_type == "regression":
                result = self.regression_power(**params)
            else:
                raise ValueError(f"Unknown test_type: {test_type}")

            results.append(
                {
                    varying_param: param_value,
                    "effect_size": result.effect_size,
                    "alpha": result.alpha,
                    "power": result.power,
                    "sample_size": result.sample_size,
                }
            )

        except ValueError as e:
            logger.debug(
                "Skipping invalid parameter combination (%s=%s): %s",
                varying_param,
                param_value,
                e,
            )
            continue

    if not results:
        raise ValueError("No valid parameter combinations found")

    # Create plot
    fig, ax = plt.subplots(figsize=(10, 6))

    x_values = [r[varying_param] for r in results]

    if varying_param == "power":
        y_values = [r["effect_size"] for r in results]
        ax.set_ylabel("Effect Size")
    elif (
        varying_param == "effect_size"
        or varying_param == "sample_size"
        or varying_param == "alpha"
    ):
        y_values = [r["power"] for r in results]
        ax.set_ylabel("Statistical Power")
    else:
        y_values = [r["power"] for r in results]
        ax.set_ylabel("Statistical Power")

    ax.plot(x_values, y_values, "b-", linewidth=2, marker="o", markersize=4)
    ax.set_xlabel(varying_param.replace("_", " ").title())
    ax.set_title(f"Power Curve: {test_type.replace('_', ' ').title()}")
    ax.grid(True, alpha=0.3)

    # Add reference lines
    if varying_param != "power":
        ax.axhline(
            y=0.8, color="red", linestyle="--", alpha=0.7, label="Power = 0.8"
        )
        ax.axhline(
            y=0.9, color="orange", linestyle="--", alpha=0.7, label="Power = 0.9"
        )
        ax.legend()

    plt.tight_layout()

    return {
        "results": results,
        "figure": fig,
        "varying_param": varying_param,
        "test_type": test_type,
    }

factorial_power_curve ¶

factorial_power_curve(effect_sizes: list[float], alpha: float = 0.05, replicates: int = 1, factor_levels: list[int] | None = None, effect: tuple[int, ...] | None = None) -> dict[str, Any]

Generate power curve for factorial designs over effect sizes.

Parameters:

Name	Type	Description	Default
`effect_sizes`	`list of float`	Effect sizes to evaluate.	required
`alpha`	`float`	Type I error rate. Defaults to 0.05.	`0.05`
`replicates`	`int`	Number of replicates per treatment combination. Defaults to `1`.	`1`
`factor_levels`	`list of int`	Number of levels for each factor. Defaults to `[2, 2]`.	`None`
`effect`	`tuple of int`	Indices of factors forming the effect of interest. Defaults to `(0,)`.	`None`

Returns:

Type	Description
`Dict[str, Any]`	Power curve data and the generated figure.

Source code in src/industrialstats/analysis/power_analysis.py

def factorial_power_curve(
    self,
    effect_sizes: list[float],
    alpha: float = 0.05,
    replicates: int = 1,
    factor_levels: list[int] | None = None,
    effect: tuple[int, ...] | None = None,
) -> dict[str, Any]:
    """Generate power curve for factorial designs over effect sizes.

    Parameters
    ----------
    effect_sizes : list of float
        Effect sizes to evaluate.
    alpha : float, optional
        Type I error rate. Defaults to 0.05.
    replicates : int, optional
        Number of replicates per treatment combination. Defaults to ``1``.
    factor_levels : list of int, optional
        Number of levels for each factor. Defaults to ``[2, 2]``.
    effect : tuple of int, optional
        Indices of factors forming the effect of interest. Defaults to
        ``(0,)``.

    Returns
    -------
    Dict[str, Any]
        Power curve data and the generated figure.
    """
    factor_levels = list(factor_levels) if factor_levels is not None else [2, 2]

    powers = []
    for es in effect_sizes:
        result = self.factorial_power(
            effect_size=es,
            alpha=alpha,
            replicates=replicates,
            factor_levels=factor_levels,
            effect=effect,
        )
        powers.append(result.power)

    fig, ax = plt.subplots(figsize=(10, 6))
    ax.plot(effect_sizes, powers, marker="o")
    ax.set_xlabel("Effect Size")
    ax.set_ylabel("Statistical Power")
    ax.set_title("Factorial Design Power Curve")
    ax.grid(True, alpha=0.3)

    ax.axhline(0.8, color="red", linestyle="--", alpha=0.7, label="Power = 0.8")
    ax.axhline(0.9, color="orange", linestyle="--", alpha=0.7, label="Power = 0.9")
    ax.legend()

    return {
        "effect_sizes": effect_sizes,
        "powers": powers,
        "figure": fig,
        "factor_levels": factor_levels,
        "effect": effect if effect is not None else (0,),
    }

sample_size_table ¶

sample_size_table(test_type: str, effect_sizes: list[float], powers: list[float] | None = None, alpha: float = 0.05, **kwargs) -> DataFrame

Generate sample size table for different effect sizes and powers.

Parameters:

Name	Type	Description	Default
`test_type`	`str`	Type of test.	required
`effect_sizes`	`list of float`	Effect sizes to include.	required
`powers`	`list of float`	Power levels to include. Defaults to `[0.8, 0.9, 0.95]`.	`None`
`alpha`	`float`	Type I error rate. Defaults to 0.05.	`0.05`
`**kwargs`		Additional parameters for specific tests.	`{}`

Returns:

Type	Description
`DataFrame`	Sample size table.

Source code in src/industrialstats/analysis/power_analysis.py

def sample_size_table(
    self,
    test_type: str,
    effect_sizes: list[float],
    powers: list[float] | None = None,
    alpha: float = 0.05,
    **kwargs,
) -> pd.DataFrame:
    """Generate sample size table for different effect sizes and powers.

    Parameters
    ----------
    test_type : str
        Type of test.
    effect_sizes : list of float
        Effect sizes to include.
    powers : list of float, optional
        Power levels to include. Defaults to ``[0.8, 0.9, 0.95]``.
    alpha : float, optional
        Type I error rate. Defaults to 0.05.
    **kwargs
        Additional parameters for specific tests.

    Returns
    -------
    pd.DataFrame
        Sample size table.
    """
    powers = list(powers) if powers is not None else [0.8, 0.9, 0.95]

    table_data = []

    for effect_size in effect_sizes:
        row = {"Effect_Size": effect_size}

        for power in powers:
            try:
                # Calculate sample size
                if test_type == "t_test":
                    result = self.t_test_power(
                        effect_size=effect_size, alpha=alpha, power=power, **kwargs
                    )
                elif test_type == "anova":
                    result = self.anova_power(
                        effect_size=effect_size, alpha=alpha, power=power, **kwargs
                    )
                elif test_type == "factorial":
                    result = self.factorial_power(
                        effect_size=effect_size, alpha=alpha, power=power, **kwargs
                    )
                elif test_type == "regression":
                    result = self.regression_power(
                        effect_size=effect_size, alpha=alpha, power=power, **kwargs
                    )
                else:
                    raise ValueError(f"Unknown test_type: {test_type}")

                row[f"Power_{power}"] = result.sample_size

            except ValueError as e:
                logger.debug("Power calculation failed for power=%s: %s", power, e)
                row[f"Power_{power}"] = np.nan

        table_data.append(row)

    return pd.DataFrame(table_data)

minimum_detectable_effect ¶

minimum_detectable_effect(test_type: str, alpha: float = 0.05, power: float = 0.8, sample_size: int = 20, **kwargs) -> PowerAnalysisResult

Calculate minimum detectable effect for a given design.

Parameters:

Name	Type	Description	Default
`test_type`	`str`	Type of test.	required
`alpha`	`float`	Type I error rate. Defaults to 0.05.	`0.05`
`power`	`float`	Statistical power. Defaults to 0.8.	`0.8`
`sample_size`	`int`	Sample size. Defaults to 20.	`20`
`**kwargs`		Additional parameters for specific tests.	`{}`

Returns:

Type	Description
`PowerAnalysisResult`	Analysis with minimum detectable effect.

Source code in src/industrialstats/analysis/power_analysis.py

def minimum_detectable_effect(
    self,
    test_type: str,
    alpha: float = 0.05,
    power: float = 0.8,
    sample_size: int = 20,
    **kwargs,
) -> PowerAnalysisResult:
    """Calculate minimum detectable effect for a given design.

    Parameters
    ----------
    test_type : str
        Type of test.
    alpha : float, optional
        Type I error rate. Defaults to 0.05.
    power : float, optional
        Statistical power. Defaults to 0.8.
    sample_size : int, optional
        Sample size. Defaults to 20.
    **kwargs
        Additional parameters for specific tests.

    Returns
    -------
    PowerAnalysisResult
        Analysis with minimum detectable effect.
    """
    if test_type == "t_test":
        return self.t_test_power(
            alpha=alpha, power=power, sample_size=sample_size, **kwargs
        )
    if test_type == "anova":
        return self.anova_power(
            alpha=alpha, power=power, sample_size=sample_size, **kwargs
        )
    if test_type == "factorial":
        return self.factorial_power(
            alpha=alpha, power=power, replicates=sample_size, **kwargs
        )
    if test_type == "regression":
        return self.regression_power(
            alpha=alpha, power=power, sample_size=sample_size, **kwargs
        )
    raise ValueError(f"Unknown test_type: {test_type}")

summary_report ¶

summary_report() -> str

Generate summary report of all power analyses performed.

Source code in src/industrialstats/analysis/power_analysis.py

def summary_report(self) -> str:
    """Generate summary report of all power analyses performed."""
    if not self.results_history:
        return "No power analyses performed yet."

    report = "POWER ANALYSIS SUMMARY REPORT\n"
    report += "=" * 50 + "\n\n"

    for i, result in enumerate(self.results_history, 1):
        report += f"Analysis {i}: {result.test_type.replace('_', ' ').title()}\n"
        report += "-" * 30 + "\n"
        report += f"Effect Size: {result.effect_size:.4f}\n"
        report += f"Alpha: {result.alpha}\n"
        report += f"Power: {result.power:.4f}\n"
        report += f"Sample Size: {result.sample_size}\n"

        # Add test-specific information
        if "n_groups" in result.additional_info:
            report += f"Number of Groups: {result.additional_info['n_groups']}\n"
        if "n_factors" in result.additional_info:
            report += f"Number of Factors: {result.additional_info['n_factors']}\n"
        if "n_predictors" in result.additional_info:
            report += (
                f"Number of Predictors: {result.additional_info['n_predictors']}\n"
            )

        report += "\n"

    return report

Analysis¶

ANOVA¶

industrialstats.analysis.anova ¶

ANOVAAnalysis ¶

fit_model ¶

anova_table_calculation ¶

multiple_comparisons ¶

residual_analysis ¶

assumptions_tests ¶

model_summary ¶

contrast_analysis ¶

power_analysis_post_hoc ¶

mixed_effects_model ¶

unbalanced_anova ¶

nested_anova ¶

repeated_measures_anova ¶

Split-plot inference¶

industrialstats.analysis.split_plot ¶

SplitPlotErrorStrata dataclass ¶

SplitPlotAnalysis ¶

error_strata ¶

anova_table ¶

expected_mean_squares ¶

fit_mixed_model ¶

Effects¶

industrialstats.analysis.effects ¶

EffectsAnalysis ¶

calculate_main_effects ¶

calculate_interaction_effects ¶

Model fitting¶

industrialstats.analysis.model_fitting ¶

ModelFitting ¶

stepwise_selection ¶

hierarchical_fitting ¶

all_subsets_selection ¶

cross_validation ¶

bootstrap_validation ¶

regularized_fitting ¶

model_comparison ¶

residual_diagnostics ¶

lack_of_fit_test ¶

Diagnostics¶

industrialstats.analysis.diagnostics ¶

ModelDiagnostics ¶

assumption_tests ¶

influence_analysis ¶

outlier_detection ¶

model_adequacy ¶

recommendation_system ¶

Power analysis¶

industrialstats.analysis.power_analysis ¶

PowerAnalysisResult dataclass ¶

PowerAnalysis ¶

t_test_power ¶

anova_power ¶

factorial_power ¶

regression_power ¶

power_curve ¶

factorial_power_curve ¶

sample_size_table ¶

minimum_detectable_effect ¶

summary_report ¶

SplitPlotErrorStrata `dataclass` ¶

PowerAnalysisResult `dataclass` ¶