Regression in Python: Statsmodels Workflow, Formulas, Diagnostics and Interpretation

Code: Expand Only the Software You Need

Standardized Effects in Code

Python checkpoint. how the script compares predictors measured on different scales. G2 has standardized beta 0.797983, G1 has 0.121852 and the remaining absolute beta values are below 0.044 in the saved model.

Statistical reading. The standardized display is a transformation of the fitted slopes and observed standard deviations, not a separate regression.

Decision rule. Use standardized betas to discuss relative modeled strength while retaining B for practical unit changes.

Verification step 7. Compare the chart values with the absolute_standardized_beta column rather than sorting on raw B.

Code and Formula Blocks

import pandas as pd
import statsmodels.api as sm

df = pd.read_csv("dataset.csv")
cols = ["G1","G2","studytime","failures","absences","age","Medu","Fedu"]
X = sm.add_constant(df[cols], has_constant="add")
fit = sm.OLS(df["G3"], X, missing="drop").fit()
pred = fit.predict(X)
resid = df.loc[pred.index, "G3"] - pred
print(fit.summary())

robust = fit.get_robustcov_results(cov_type="HC3")
print(robust.summary())

Verified Python Core Implementation

def prepare_data(
    data: pd.DataFrame,
) -> tuple[pd.DataFrame, str, str, list[str]]:
    outcome = choose_outcome(data)
    simple_predictor = choose_simple_predictor(data, outcome)
    multiple_predictors = choose_multiple_predictors(
        data,
        outcome,
        simple_predictor,
    )

    keep_columns = [outcome] + list(
        dict.fromkeys([simple_predictor] + multiple_predictors)
    )

    analysis_data = data[keep_columns].copy()

    for column in analysis_data.columns:
        analysis_data[column] = pd.to_numeric(
            analysis_data[column],
            errors="coerce",
        )

    analysis_data = (
        analysis_data
        .replace([np.inf, -np.inf], np.nan)
        .dropna()
        .reset_index(drop=True)
    )

    if len(analysis_data) <= len(multiple_predictors) + 10:
        raise ValueError(
            "Not enough complete observations for the requested regression models."
        )

    return (
        analysis_data,
        outcome,
        simple_predictor,
        multiple_predictors,
    )


def fit_ols(
    data: pd.DataFrame,
    outcome: str,
    predictors: list[str],
):
    y = data[outcome].astype(float)
    x = sm.add_constant(
        data[predictors].astype(float),
        has_constant="add",
    )
    return sm.OLS(y, x).fit()


def make_model_summary(
    result,
    model_name: str,
    outcome: str,
    predictors: list[str],
) -> pd.DataFrame:
    residuals = np.asarray(result.resid, dtype=float)

    return pd.DataFrame(
        [
            {
                "model": model_name,
                "outcome": outcome,
                "predictors": ", ".join(predictors),
                "observations": int(result.nobs),
                "number_of_predictors": len(predictors),
                "r_squared": float(result.rsquared),
                "adjusted_r_squared": float(result.rsquared_adj),
                "f_statistic": float(result.fvalue),
                "f_p_value": float(result.f_pvalue),
                "aic": float(result.aic),
                "bic": float(result.bic),
                "rmse": float(np.sqrt(np.mean(residuals ** 2))),
                "mae": float(np.mean(np.abs(residuals))),
                "residual_standard_error": float(
                    np.sqrt(np.sum(residuals ** 2) / result.df_resid)
                ),
                "durbin_watson": float(durbin_watson(residuals)),
                "decision_alpha_0_05": (
                    "Model is statistically significant"
                    if result.f_pvalue < ALPHA
                    else "Model is not statistically significant"
                ),
            }
        ]
    )


def make_coefficient_table(
    result,
    model_name: str,
    outcome_standard_deviation: float,
    predictor_standard_deviations: dict[str, float],
) -> pd.DataFrame:
    confidence_intervals = result.conf_int(alpha=ALPHA)

    rows: list[dict] = []

    for term in result.params.index:
        coefficient = float(result.params[term])
        standard_error = float(result.bse[term])
        t_value = float(result.tvalues[term])
        p_value = float(result.pvalues[term])
        lower = float(confidence_intervals.loc[term, 0])
        upper = float(confidence_intervals.loc[term, 1])

        display_term = "Intercept" if term == "const" else term

        if term == "const":
            standardized_beta = np.nan
        else:
            standardized_beta = (
                coefficient
                * predictor_standard_deviations[term]
                / outcome_standard_deviation
            )

        rows.append(
            {
                "model": model_name,
                "term": display_term,
                "coefficient_b": coefficient,
                "standard_error": standard_error,
                "t_value": t_value,
                "p_value": p_value,
                "ci95_lower": lower,
                "ci95_upper": upper,
                "standardized_beta": standardized_beta,
                "absolute_standardized_beta": (
                    abs(standardized_beta)
                    if np.isfinite(standardized_beta)
                    else np.nan
                ),
                "decision_alpha_0_05": (
                    "Statistically significant"
                    if p_value < ALPHA
                    else "Not statistically significant"
                ),
            }
        )

    return pd.DataFrame(rows)


def make_anova_table(result, model_name: str) -> pd.DataFrame:
    observed = np.asarray(result.model.endog, dtype=float)
    fitted = np.asarray(result.fittedvalues, dtype=float)
    residuals = np.asarray(result.resid, dtype=float)

    total_sum_squares = float(
        np.sum((observed - np.mean(observed)) ** 2)
    )
    regression_sum_squares = float(
        np.sum((fitted - np.mean(observed)) ** 2)
    )
    residual_sum_squares = float(
        np.sum(residuals ** 2)
    )

    regression_df = int(result.df_model)
    residual_df = int(result.df_resid)
    total_df = regression_df + residual_df

    regression_mean_square = (
        regression_sum_squares / regression_df
        if regression_df > 0
        else np.nan
    )
    residual_mean_square = (
        residual_sum_squares / residual_df
        if residual_df > 0
        else np.nan
    )

    return pd.DataFrame(
        [
            {
                "model": model_name,
                "source": "Regression",
                "sum_of_squares": regression_sum_squares,
                "df": regression_df,
                "mean_square": regression_mean_square,
                "f_value": float(result.fvalue),
                "p_value": float(result.f_pvalue),
            },
            {
                "model": model_name,
                "source": "Residual",
                "sum_of_squares": residual_sum_squares,
                "df": residual_df,
                "mean_square": residual_mean_square,
                "f_value": np.nan,
                "p_value": np.nan,
            },
            {
                "model": model_name,
                "source": "Total",
                "sum_of_squares": total_sum_squares,
                "df": total_df,
                "mean_square": np.nan,
                "f_value": np.nan,
                "p_value": np.nan,
            },
        ]
    )


def make_vif_table(
    data: pd.DataFrame,
    predictors: list[str],
) -> pd.DataFrame:
    design = sm.add_constant(
        data[predictors].astype(float),
        has_constant="add",
    )

    rows: list[dict] = []

    for index, term in enumerate(design.columns):
        if term == "const":
            continue

        vif = float(
            variance_inflation_factor(
                design.to_numpy(dtype=float),
                index,
            )
        )

        rows.append(
            {
                "predictor": term,
                "vif": vif,
                "tolerance": 1 / vif if vif != 0 else np.nan,
                "diagnostic": (
                    "Potential multicollinearity"
                    if vif >= 5
                    else "Acceptable"
                ),
            }
        )

    return pd.DataFrame(rows)


def make_diagnostic_tests(
    result,
    model_name: str,
) -> pd.DataFrame:
    residuals = np.asarray(result.resid, dtype=float)
    exogenous = np.asarray(result.model.exog, dtype=float)

    if len(residuals) >= 3:
        shapiro_statistic, shapiro_p_value = stats.shapiro(
            residuals[:5000]
        )
    else:
        shapiro_statistic, shapiro_p_value = np.nan, np.nan

    jarque_bera = stats.jarque_bera(residuals)

    try:
        bp_statistic, bp_p_value, bp_f_statistic, bp_f_p_value = (
            het_breuschpagan(
                residuals,
                exogenous,
            )
        )
    except Exception:
        bp_statistic = np.nan
        bp_p_value = np.nan
        bp_f_statistic = np.nan
        bp_f_p_value = np.nan

    return pd.DataFrame(
        [
            {
                "model": model_name,
                "diagnostic": "Shapiro-Wilk residual normality",
                "statistic": shapiro_statistic,
                "p_value": shapiro_p_value,
                "interpretation": (
                    "Potential non-normality"
                    if shapiro_p_value < ALPHA
                    else "No strong evidence against normality"
                ),
            },
            {
                "model": model_name,
                "diagnostic": "Jarque-Bera residual normality",
                "statistic": jarque_bera.statistic,
                "p_value": jarque_bera.pvalue,
                "interpretation": (
                    "Potential non-normality"
                    if jarque_bera.pvalue < ALPHA
                    else "No strong evidence against normality"
                ),
            },
            {
                "model": model_name,
                "diagnostic": "Breusch-Pagan heteroscedasticity",
                "statistic": bp_statistic,
                "p_value": bp_p_value,
                "interpretation": (
                    "Potential heteroscedasticity"
                    if bp_p_value < ALPHA
                    else "No strong evidence of heteroscedasticity"
                ),
            },
            {
                "model": model_name,
                "diagnostic": "Breusch-Pagan F test",
                "statistic": bp_f_statistic,
                "p_value": bp_f_p_value,
                "interpretation": (
                    "Potential heteroscedasticity"
                    if bp_f_p_value < ALPHA
                    else "No strong evidence of heteroscedasticity"
                ),
            },
            {
                "model": model_name,
                "diagnostic": "Durbin-Watson residual independence",
                "statistic": float(durbin_watson(residuals)),
                "p_value": np.nan,
                "interpretation": (
                    "Values near 2 indicate approximate independence"
                ),
            },
        ]
    )


def make_prediction_table(
    data: pd.DataFrame,
    outcome: str,
    simple_predictor: str,
    multiple_predictors: list[str],
    simple_result,
    multiple_result,
) -> pd.DataFrame:
    multiple_influence = multiple_result.get_influence()
    prediction_frame = multiple_result.get_prediction().summary_frame(
        alpha=ALPHA
    )

    output = data[
        [outcome] + list(
            dict.fromkeys([simple_predictor] + multiple_predictors)
        )
    ].copy()

    output.insert(
        0,
        "case_number",
        np.arange(1, len(output) + 1),
    )

    output["simple_predicted_value"] = (
        simple_result.fittedvalues.to_numpy()
    )
    output["simple_residual"] = (
        simple_result.resid.to_numpy()
    )

    output["multiple_predicted_value"] = (
        multiple_result.fittedvalues.to_numpy()
    )
    output["multiple_residual"] = (
        multiple_result.resid.to_numpy()
    )
    output["multiple_standardized_residual"] = (
        multiple_influence.resid_studentized_internal
    )
    output["multiple_external_studentized_residual"] = (
        multiple_influence.resid_studentized_external
    )
    output["multiple_leverage"] = (
        multiple_influence.hat_matrix_diag
    )
    output["multiple_cooks_distance"] = (
        multiple_influence.cooks_distance[0]
    )
    output["multiple_dffits"] = (
        multiple_influence.dffits[0]
    )

    output["mean_ci_lower_95"] = prediction_frame["mean_ci_lower"].to_numpy()
    output["mean_ci_upper_95"] = prediction_frame["mean_ci_upper"].to_numpy()
    output["prediction_interval_lower_95"] = (
        prediction_frame["obs_ci_lower"].to_numpy()
    )
    output["prediction_interval_upper_95"] = (
        prediction_frame["obs_ci_upper"].to_numpy()
    )

    return output


def add_title_and_subtitle(
    figure,
    title: str,
    subtitle: str,
) -> None:
    figure.suptitle(
        title,
        x=0.06,
        y=0.975,
        ha="left",
        fontsize=18,
        fontweight="bold",
    )
    figure.text(
        0.06,
        0.925,
        subtitle,
        ha="left",
        fontsize=10.5,
    )
    figure.subplots_adjust(
        top=0.82,
        left=0.11,
        right=0.95,
        bottom=0.15,
    )


def save_figure(
    figure,
    filename: str,
) -> Path:
    path = PNG_DIR / filename
    figure.savefig(
        path,
        dpi=300,
        bbox_inches="tight",
    )
    plt.close(figure)
    return path


def create_charts(
    data: pd.DataFrame,
    outcome: str,
    simple_predictor: str,
    multiple_predictors: list[str],
    simple_result,
    multiple_result,
    simple_summary: pd.DataFrame,
    multiple_summary: pd.DataFrame,
    multiple_coefficients: pd.DataFrame,
    predictions: pd.DataFrame,
) -> list[Path]:
    charts: list[Path] = []

    figure, axis = plt.subplots(figsize=(9.5, 6.8))
    axis.hist(
        data[outcome],
        bins=16,
        edgecolor="white",
        alpha=0.9,
        color=CHART_COLORS[0],
    )
    axis.axvline(
        data[outcome].mean(),
        linewidth=2,
        color=CHART_COLORS[2],
        label=f"Mean = {data[outcome].mean():.2f}",
    )
    axis.set_xlabel(outcome)
    axis.set_ylabel("Frequency")
    axis.legend(frameon=False)
    axis.grid(axis="y", alpha=0.3)
    add_title_and_subtitle(
        figure,
        "Outcome Distribution",
        "Distribution of the dependent variable used in the regression analyses.",
    )
    charts.append(
        save_figure(
            figure,
            "01_regression_in_python_outcome_distribution.png",
        )
    )

    figure, axis = plt.subplots(figsize=(9.5, 6.8))
    axis.scatter(
        data[simple_predictor],
        data[outcome],
        alpha=0.65,
        color=CHART_COLORS[1],
        edgecolor="white",
        linewidth=0.5,
        label="Observed",
    )
    order = np.argsort(data[simple_predictor].to_numpy())
    axis.plot(
        data[simple_predictor].to_numpy()[order],
        simple_result.fittedvalues.to_numpy()[order],
        linewidth=2.5,
        color=CHART_COLORS[3],
        label="Regression line",
    )
    axis.set_xlabel(simple_predictor)
    axis.set_ylabel(outcome)
    axis.legend(frameon=False)
    axis.grid(alpha=0.3)
    add_title_and_subtitle(
        figure,
        "Simple Linear Regression",
        f"Observed {outcome} values and fitted regression line for {simple_predictor}.",
    )
    charts.append(
        save_figure(
            figure,
            "02_simple_regression_scatter_with_line.png",
        )
    )

    figure, axis = plt.subplots(figsize=(8.5, 7))
    axis.scatter(
        data[outcome],
        multiple_result.fittedvalues,
        alpha=0.65,
        color=CHART_COLORS[2],
        edgecolor="white",
        linewidth=0.5,
    )
    minimum = min(
        data[outcome].min(),
        multiple_result.fittedvalues.min(),
    )
    maximum = max(
        data[outcome].max(),
        multiple_result.fittedvalues.max(),
    )
    axis.plot(
        [minimum, maximum],
        [minimum, maximum],
        linestyle="--",
        linewidth=2,
        color=CHART_COLORS[4],
    )
    axis.set_xlabel(f"Observed {outcome}")
    axis.set_ylabel(f"Predicted {outcome}")
    axis.grid(alpha=0.3)
    add_title_and_subtitle(
        figure,
        "Observed versus Predicted Values",
        "Points closer to the diagonal line indicate stronger prediction accuracy.",
    )
    charts.append(
        save_figure(
            figure,
            "03_multiple_regression_observed_vs_predicted.png",
        )
    )

    figure, axis = plt.subplots(figsize=(9.5, 6.8))
    axis.scatter(
        multiple_result.fittedvalues,
        multiple_result.resid,
        alpha=0.65,
        color=CHART_COLORS[3],
        edgecolor="white",
        linewidth=0.5,
    )
    axis.axhline(
        0,
        linestyle="--",
        linewidth=2,
        color=CHART_COLORS[4],
    )
    axis.set_xlabel("Fitted values")
    axis.set_ylabel("Residuals")
    axis.grid(alpha=0.3)
    add_title_and_subtitle(
        figure,
        "Residuals versus Fitted Values",
        "A random band around zero supports linearity and constant variance.",
    )
    charts.append(
        save_figure(
            figure,
            "04_multiple_regression_residuals_vs_fitted.png",
        )
    )

    figure, axis = plt.subplots(figsize=(9.5, 6.8))
    axis.hist(
        multiple_result.resid,
        bins=18,
        edgecolor="white",
        color=CHART_COLORS[4],
        alpha=0.9,
    )
    axis.axvline(
        0,
        linestyle="--",
        linewidth=2,
        color=CHART_COLORS[2],
    )
    axis.set_xlabel("Residual")
    axis.set_ylabel("Frequency")
    axis.grid(axis="y", alpha=0.3)
    add_title_and_subtitle(
        figure,
        "Residual Distribution",
        "A centered and approximately symmetric distribution supports residual normality.",
    )
    charts.append(
        save_figure(
            figure,
            "05_multiple_regression_residual_histogram.png",
        )
    )

    figure = plt.figure(figsize=(8.5, 7))
    axis = figure.add_subplot(111)
    stats.probplot(
        multiple_result.resid,
        dist="norm",
        plot=axis,
    )
    axis.grid(alpha=0.3)
    add_title_and_subtitle(
        figure,
        "Q-Q Plot of Residuals",
        "Points close to the reference line support approximate normality.",
    )
    charts.append(
        save_figure(
            figure,
            "06_multiple_regression_residual_qq_plot.png",
        )
    )

    coefficient_plot = (
        multiple_coefficients[
            multiple_coefficients["term"] != "Intercept"
        ]
        .copy()
        .sort_values("standardized_beta")
    )

    figure, axis = plt.subplots(
        figsize=(
            10.5,
            max(5.5, 0.55 * len(coefficient_plot) + 2.5),
        )
    )
    positions = np.arange(len(coefficient_plot))
    axis.barh(
        positions,
        coefficient_plot["standardized_beta"],
        color=CHART_COLORS[5],
    )
    axis.set_yticks(
        positions,
        coefficient_plot["term"],
    )
    axis.axvline(
        0,
        linestyle="--",
        linewidth=1.5,
        color=CHART_COLORS[2],
    )
    axis.set_xlabel("Standardized beta coefficient")
    axis.grid(axis="x", alpha=0.3)
    add_title_and_subtitle(
        figure,
        "Standardized Beta Coefficients",
        "Absolute standardized beta values compare relative predictor importance.",
    )
    charts.append(
        save_figure(
            figure,
            "07_multiple_regression_standardized_betas.png",
        )
    )

    figure, axis = plt.subplots(figsize=(9.5, 6.8))
    model_labels = ["Simple regression", "Multiple regression"]
    r_squared_values = [
        simple_summary.loc[0, "r_squared"],
        multiple_summary.loc[0, "r_squared"],
    ]
    adjusted_values = [
        simple_summary.loc[0, "adjusted_r_squared"],
        multiple_summary.loc[0, "adjusted_r_squared"],
    ]
    x_positions = np.arange(2)
    width = 0.36
    axis.bar(
        x_positions - width / 2,
        r_squared_values,
        width,
        label="R-squared",
        color=CHART_COLORS[0],
    )
    axis.bar(
        x_positions + width / 2,
        adjusted_values,
        width,
        label="Adjusted R-squared",
        color=CHART_COLORS[1],
    )
    axis.set_xticks(
        x_positions,
        model_labels,
    )
    axis.set_ylim(
        0,

Report assembly and reproducible export

        min(1.05, max(r_squared_values) + 0.15),
    )
    axis.set_ylabel("Model fit")
    axis.legend(frameon=False)
    axis.grid(axis="y", alpha=0.3)
    add_title_and_subtitle(
        figure,
        "Simple and Multiple Regression Fit",
        "Adjusted R-squared accounts for the number of predictors.",
    )
    charts.append(
        save_figure(
            figure,
            "08_simple_vs_multiple_regression_fit.png",
        )
    )

    figure, axis = plt.subplots(figsize=(9.5, 6.8))
    axis.scatter(
        predictions["multiple_leverage"],
        predictions["multiple_cooks_distance"],
        alpha=0.65,
        color=CHART_COLORS[6],
        edgecolor="white",
        linewidth=0.5,
    )
    leverage_threshold = (
        2 * (len(multiple_predictors) + 1) / len(predictions)
    )
    cook_threshold = 4 / len(predictions)
    axis.axvline(
        leverage_threshold,
        linestyle="--",
        linewidth=1.8,
        label="Leverage 2p/n",
        color=CHART_COLORS[3],
    )
    axis.axhline(
        cook_threshold,
        linestyle=":",
        linewidth=1.8,
        label="Cook 4/n",
        color=CHART_COLORS[2],
    )
    axis.set_xlabel("Leverage")
    axis.set_ylabel("Cook's distance")
    axis.legend(frameon=False)
    axis.grid(alpha=0.3)
    add_title_and_subtitle(
        figure,
        "Leverage and Cook's Distance",
        "Cases beyond the reference lines should be investigated for influence.",
    )
    charts.append(
        save_figure(
            figure,
            "09_multiple_regression_influence_diagnostics.png",
        )
    )

    figure, axis = plt.subplots(figsize=(11, 9))
    correlation_matrix = data[
        [outcome] + multiple_predictors
    ].corr()

    image = axis.imshow(
        correlation_matrix,
        vmin=-1,
        vmax=1,
        cmap="coolwarm",
    )
    axis.set_xticks(
        np.arange(len(correlation_matrix.columns)),
        correlation_matrix.columns,
        rotation=45,
        ha="right",
    )
    axis.set_yticks(
        np.arange(len(correlation_matrix.index)),
        correlation_matrix.index,
    )
    figure.colorbar(
        image,
        ax=axis,
        fraction=0.046,
        pad=0.04,
        label="Correlation",
    )
    add_title_and_subtitle(
        figure,
        "Regression Variable Correlation Matrix",
        "Strong predictor correlations may affect coefficient stability.",
    )
    charts.append(
        save_figure(
            figure,
            "10_regression_variable_correlation_heatmap.png",
        )
    )

    return charts


def create_pdf_report(
    report_lines: list[str],
    tables: dict[str, pd.DataFrame],
    chart_paths: list[Path],
) -> Path:
    pdf_path = PDF_DIR / "Regression-in-Python-Report.pdf"

    with PdfPages(pdf_path) as pdf:
        figure = plt.figure(figsize=(11, 8.5))
        figure.text(
            0.06,
            0.94,
            TEST_NAME,
            fontsize=22,
            fontweight="bold",
        )

        y_position = 0.88
        for line in report_lines[:28]:
            figure.text(
                0.06,
                y_position,
                line,
                fontsize=10,
            )
            y_position -= 0.034

        pdf.savefig(
            figure,
            bbox_inches="tight",
        )
        plt.close(figure)

        for title, table_data in tables.items():
            figure, axis = plt.subplots(figsize=(11, 8.5))
            axis.axis("off")

            display = table_data.copy().head(14)

            for column in display.columns:
                if pd.api.types.is_numeric_dtype(display[column]):
                    display[column] = display[column].map(
                        lambda value: (
                            ""
                            if pd.isna(value)
                            else f"{value:.5f}"
                        )
                    )

            axis.set_title(
                title,
                fontsize=16,
                fontweight="bold",
                loc="left",
                pad=20,
            )

            table = axis.table(
                cellText=display.values,
                colLabels=display.columns,
                cellLoc="left",
                colLoc="left",
                loc="center",
            )
            table.auto_set_font_size(False)
            table.set_fontsize(6.8)
            table.scale(1, 1.2)

            pdf.savefig(
                figure,
                bbox_inches="tight",
            )
            plt.close(figure)

        for chart_path in chart_paths:
            image = plt.imread(chart_path)

            figure, axis = plt.subplots(figsize=(11, 8.5))
            axis.imshow(image)
            axis.axis("off")
            axis.set_title(
                chart_path.name,
                fontsize=13,
                fontweight="bold",
                loc="left",
            )

            pdf.savefig(
                figure,
                bbox_inches="tight",
            )
            plt.close(figure)

    return pdf_path


def main() -> None:

Advanced Regression in Python Topics

Pandas Import and Row Integrity

Python checkpoint. how data enter the pipeline without silent coercion. The saved dataset overview, variable summary and analysis-data tables document the rows and model columns used by the script.

Statistical reading. A reproducible fit depends on numeric conversion, missing-value handling and column selection occurring before the design matrix is handed to statsmodels.

Decision rule. Fail loudly when a required column is absent or contains unexpected text instead of allowing a reduced model to run.

Verification step 2. Compare the 649-row count across the overview, coefficient model metadata and predictions table.

Statsmodels OLS Design Matrix

Python checkpoint. why the intercept must be added deliberately. The coefficient output contains an intercept for both simple and multiple equations, with -0.501155 in the latter. The model includes eight slopes.

Statistical reading. An omitted constant changes residual centering, R-squared interpretation and every coefficient; adding it twice can also corrupt the intended design.

Decision rule. Construct X in an explicit column order and call add_constant once before fitting.

Verification step 3. Print exogenous names from the fitted object and reconcile them with the saved coefficient table.

Variance Inflation Factors

Python checkpoint. whether the Python design shows damaging predictor redundancy. The largest VIF values are 4.164832 for G1 and 4.120795 for G2; all named predictor rows are labeled acceptable in the saved table.

Statistical reading. These values reveal substantial overlap between the two grade predictors but remain below the workflow’s closer-review threshold of five.

Decision rule. Describe the overlap and inspect stability rather than claiming that an acceptable label proves coefficients are immune to collinearity.

Verification step 8. Recompute each VIF on the same eight-column matrix and exclude the intercept from substantive interpretation.

Model Comparison without Overclaiming

Python checkpoint. how AIC, BIC and fit metrics frame the two equations. The simple model has AIC 2162.270813 and BIC 2171.221679; the multiple model has AIC 2146.324659 and BIC 2186.603553.

Statistical reading. AIC favors the multiple fit while BIC penalizes its additional parameters more strongly, demonstrating why one ranking criterion should be specified in advance.

Decision rule. Use out-of-sample validation for predictive claims instead of choosing a model solely from in-sample criteria.

Verification step 13. Keep the same rows and outcome when comparing criteria, as the saved pipeline does.

Verified Chart Stories and Same-Topic Evidence

Every figure below comes from the dedicated Regression in Python folder.

Discussion of What the Python Findings Mean

The most important substantive pattern is continuity between assessment periods. G2 alone captures roughly 84 percent of fitted G3 variation, and G1 adds a smaller conditional contribution once G2 is present. This is plausible because the variables measure related academic performance at successive times. It also means that high fit partly reflects temporal proximity and overlapping measurement, not necessarily broad explanatory coverage of motivation, teaching quality or home conditions.

The multiple equation improves fitted RMSE by only about 0.029 grade points relative to the simple model, despite adding seven variables. That comparison favors parsimony when transparent prediction from G2 is sufficient, but it does not make the extra terms useless. The multiple fit answers a different question about conditional associations, reveals suppression in absences, and provides a richer diagnostic surface. Choice between the two should follow intended use rather than a mechanical preference for either simplicity or maximum R-squared.

G1 and G2 VIF values near four deserve discussion because their slopes partition shared prior-grade information. A coefficient can move noticeably if one of the pair is removed even while predictions remain stable. For scientific interpretation, stability plots across resamples or alternate prior-grade specifications would be more revealing than a single VIF cutoff. For operational prediction, the question is whether both variables improve validated error enough to justify collecting and maintaining them.

Residual tests point to a mismatch between the convenient Gaussian OLS story and the observed tail. Zero final grades create departures that are substantively meaningful; they may reflect dropout, noncompletion or a grading rule rather than random noise. A sensitivity analysis could compare robust covariance, robust regression, a bounded-outcome strategy or a two-part mechanism. Each alternative changes the estimand, so it should be motivated by how zeros arise rather than selected merely to make a normality test nonsignificant.

The positive adjusted absence coefficient is an example of why automated prose is dangerous. Its magnitude is small, its interval barely excludes zero, and its sign differs from the bivariate direction shown in the coefficient-focused outputs. A responsible discussion notes the conditioning set and shared associations, checks coding and influential observations, and avoids saying that more absence improves grades. The regression describes a conditional pattern; it does not justify an intervention recommendation.

Software reproducibility adds epistemic value only when the saved artifacts remain connected. Ten CSV tables expose exact numbers, ten PNG files provide graphical views, and the PDF offers a stable review copy. If a chart is regenerated after code changes but the report is not, apparent agreement can become stale. Hashes, timestamps, a run manifest or one scripted output directory would make it easier to prove that every derivative came from the same execution.

The verified run is best understood as a strong descriptive baseline. It establishes the dominant role of recent prior grade, quantifies smaller conditional terms, and reveals residual limitations that should shape future modeling. It does not establish causal pathways, transportability to another school system or fairness across subgroups. Those questions need new design information, subgroup diagnostics and validation data rather than more decimals from the same fit.

Python Diagnostics, Influence Review and Validation Decisions

Residuals versus fitted values should be read before formal tests because the plot indicates what kind of misspecification may be present. Curvature suggests a mean-form problem; a funnel suggests changing variance; isolated vertical extremes suggest unusual outcomes. In these grade data, discrete bands and a cluster of severe negative errors matter. A single omnibus p value cannot distinguish those mechanisms, so the visual and the case export serve different but complementary roles.

The normal Q Q plot is expected to bend in the lower tail when observed zeros lie far below ordinary predictions. Shapiro-Wilk and Jarque-Bera confirm departure, but the next decision is not automatically to transform G3. A transformation may damage the interpretation of grade points and cannot make a bounded outcome unbounded. Bootstrap intervals, robust covariance or a model designed for the data-generating boundary can be considered while retaining an explicit description of what each remedy changes.

Breusch-Pagan evidence supports an HC3 sensitivity analysis for coefficient uncertainty. In statsmodels, get_robustcov_results(cov_type=’HC3′) returns a result wrapper with unchanged OLS point estimates and leverage-adjusted covariance. The robust intervals should be reported alongside, not silently substituted into, the classical table. If conclusions differ, that divergence is substantive evidence that variance misspecification affects inference.

Leverage is determined by the rows of X, so it should be inspected together with predictor values rather than outcome alone. Case 1 combines G1 equal to zero with later grade information and has leverage 0.1102142948. A sensitivity fit excluding that row can show changes in coefficients, predictions and diagnostics, but the primary analysis should retain the case unless source verification or the research protocol supplies a reason to exclude it.

Cook distance and DFFITS answer related but nonidentical questions. Cook distance summarizes joint movement of the fitted coefficient vector after deletion; DFFITS measures the deletion effect on that case’s own fitted value. Ranking both, then adding externally studentized residual and leverage, produces an interpretable review table. Thresholds are triage devices whose practical importance depends on observed coefficient change and the validity of the source record.

The correlation heatmap and VIF table are design diagnostics, not residual diagnostics. They explain why G1 and G2 compete for unique slope information, while residual plots evaluate how the final equation misses G3. Keeping those categories separate prevents a common mistake in which multicollinearity is blamed for heteroscedasticity or a non-normal tail. Different problems require different evidence and different responses.

A defensible validation extension would use nested resampling if any predictor selection, transformation or tuning is introduced. Each training fold would learn preprocessing, fit the candidate equation and produce predictions for its untouched fold. The resulting error distribution could be compared with 1.247020 fitted RMSE. Without this separation, a lower training error after adding complexity may reflect adaptation to these rows rather than improved future performance.

Python Limitations, Common Errors and Scope Boundaries

A pandas dataframe can conceal type problems until a modeling call. Numeric-looking strings, locale commas, sentinel values and mixed missing markers may coerce unpredictably. The verified script should validate ranges and dtypes explicitly, then preserve a rejected-row log if cleaning removes data. Because the current report says all 649 rows were used, an unnoticed coercion that drops even one record is a material divergence from the published result.

Notebook execution order is another risk. A cell may display results from an old X while later cells operate on a modified dataframe. Restarting the kernel and running the complete analysis top to bottom is the minimum check. A standalone script is easier to audit because object creation, fitting and export occur in one declared sequence, but it still needs guarded paths and should fail rather than reuse output files from a previous topic.

Statsmodels summary text is convenient for people but brittle for programmatic extraction. Parsing aligned console columns can misread scientific notation or lose term names. The package correctly stores structured CSVs; future updates should continue to access attributes directly and format only at the publication boundary. A PDF belongs at the end of the pipeline, never as the numerical input to a new article or calculator.

The model treats G3 as continuous and unbounded even though grades occupy a restricted scale and include exact zeros. OLS can predict outside the permitted range at extreme covariate combinations. Clipping those predictions after fitting changes operational output but does not repair the statistical model or interval calculations. If boundary behavior matters, compare a scientifically appropriate alternative and state the different target it estimates.

Classical p values assume the specified covariance and a prespecified analysis. Multiple exploratory specifications, subgroup searches or repeated outcome checks would require transparent multiplicity handling. The saved folder supports the two declared formulas only; it contains no record that dozens of alternatives were tried. The article should not imply that the reported probabilities account for searches that are absent from the evidence.

Causal interpretation is outside the fitted design. Prior grades, absences and studytime may share causes with final grade, and conditioning on intermediate variables can change the meaning of coefficients. Temporal order for G1 and G2 improves their predictive relevance but does not create random assignment. Language such as associated with, conditional on and predicts in this sample is more accurate than caused or produced.

Transportability remains unknown. The 649 records may differ from future students in grading rules, curriculum, missingness or support systems. Before deployment, coefficients and error distributions should be checked on a genuinely new cohort, with calibration examined across relevant subgroups. Model monitoring should focus on prediction error and residual patterns, not merely whether the old R-squared can be reproduced.

Statsmodels Summary Output Audit

The saved statsmodels summary output is a compact view of one fitted result object, not a substitute for checking the design matrix. In this run, the multiple model contains an intercept and eight declared predictor columns for 649 complete observations. The reconstruction confirms that the coefficient ordering in the summary table follows the explicit column list rather than an alphabetical or automatically discovered order.

The headline fit row reports multiple R-squared 0.85077719824595 and adjusted R-squared 0.848911913224024. Their small separation reflects the penalty for eight slopes at this sample size; it does not turn the adjusted value into a validation score. The F statistic 456.111097363367 and probability 1.45489804091868e-258 test the joint zero-slope null under the fitted OLS assumptions.

The coefficient portion of the summary output must be read in original units. G2 has B = 0.884807 with the other columns held constant, while G1 has B = 0.143397 under the same conditioning. Standardized magnitudes are calculated in a separate table from sample standard deviations, so the summary’s raw slopes and the beta ranking answer different numerical questions.

The diagnostic header and companion exports are reconciled rather than read independently. RMSE 1.24702023001007 is computed from case residuals with an n denominator, whereas the regression standard error uses residual degrees of freedom. Jarque-Bera 9923.74188901137 and Breusch-Pagan 35.7228073187284 are warning statistics saved by dedicated tests, not labels inferred from the appearance of a chart.

A reproducible summary requires the same row filter, response vector, constant column and predictor order. Changing any one of those inputs can change coefficient estimates while leaving valid-looking Python output. The article therefore treats the dataset overview, formula record, coefficient CSV and diagnostic tables as a connected audit trail for this exact model.