{"spec_id":"diagnostic-regression-panel","library":"plotnine","language":"python","code":"\"\"\" anyplot.ai\ndiagnostic-regression-panel: Regression Diagnostic Panel (Four-Plot Display)\nLibrary: plotnine 0.15.4 | Python 3.13.13\nQuality: 85/100 | Created: 2026-05-13\n\"\"\"\n\nimport os\nimport site\nimport sys\n\n\n# 'python plotnine.py' adds this file's directory to sys.path[0], shadowing\n# the installed plotnine library. Insert site-packages at the front so the\n# library is found before this file. ruff recognises sys.path.insert() as a\n# valid pre-import path setup and does not raise E402 for subsequent imports.\nsys.path.insert(0, site.getsitepackages()[0])\n\nimport numpy as np\nimport pandas as pd\nfrom plotnine import (\n    aes,\n    element_blank,\n    element_line,\n    element_rect,\n    element_text,\n    facet_wrap,\n    geom_hline,\n    geom_line,\n    geom_point,\n    geom_text,\n    ggplot,\n    labs,\n    scale_color_manual,\n    theme,\n)\nfrom scipy.stats import probplot\nfrom sklearn.linear_model import LinearRegression\nfrom statsmodels.nonparametric.smoothers_lowess import lowess\n\n\n# Theme tokens\nTHEME = os.getenv(\"ANYPLOT_THEME\", \"light\")\nPAGE_BG = \"#FAF8F1\" if THEME == \"light\" else \"#1A1A17\"\nELEVATED_BG = \"#FFFDF6\" if THEME == \"light\" else \"#242420\"\nINK = \"#1A1A17\" if THEME == \"light\" else \"#F0EFE8\"\nINK_SOFT = \"#4A4A44\" if THEME == \"light\" else \"#B8B7B0\"\nBRAND = \"#009E73\"\nACCENT = \"#C475FD\"\n\n# Data — materials-testing regression: tensile strength vs temperature and load\n# Three high-leverage observations at extreme temperatures reveal influential points\nnp.random.seed(42)\nn_main = 97\nn_extreme = 3\n\ntemp_main = np.random.uniform(20, 80, n_main)\ntemp_extreme = np.array([180.0, 195.0, 210.0])\ntemperature = np.concatenate([temp_main, temp_extreme])\n\nload_main = np.random.normal(50, 10, n_main)\nload_extreme = np.random.normal(50, 10, n_extreme)\nload = np.concatenate([load_main, load_extreme])\n\nn = n_main + n_extreme\nX = np.column_stack([np.ones(n), temperature, load])\np_full = X.shape[1]  # includes intercept column\n\n# True relationship: strength declines with temperature, rises with load\n# Heteroscedastic errors: variance grows with temperature\nerrors = np.random.normal(0, 2 + 0.04 * temperature, n)\nstrength = 120 - 0.3 * temperature + 0.8 * load + errors\n\n# Add two response outliers at ordinary x-values\nstrength[20] += 25.0\nstrength[55] -= 28.0\n\n# Fit OLS (without intercept column — sklearn adds intercept internally)\nX_fit = np.column_stack([temperature, load])\nmodel = LinearRegression()\nmodel.fit(X_fit, strength)\nfitted_vals = model.predict(X_fit)\nresiduals = strength - fitted_vals\nk = X_fit.shape[1]  # number of predictors (excludes intercept)\np_hat = k + 1  # number of estimated parameters (including intercept)\n\n# Hat matrix for leverage (include intercept column)\nXtXinv = np.linalg.inv(X.T @ X)\nH = X @ XtXinv @ X.T\nleverage = np.diag(H)\n\n# Internally studentized (standardized) residuals\nsigma2 = np.sum(residuals**2) / (n - p_hat)\nsigma = np.sqrt(sigma2)\nstd_residuals = residuals / (sigma * np.sqrt(np.clip(1 - leverage, 1e-10, None)))\n\n# Cook's distance\ncooks_d = (std_residuals**2 * leverage) / (p_hat * (1 - leverage))\ntop3_idx = set(np.argsort(cooks_d)[-3:])\n\n# Q-Q data\n(qq_theoretical, qq_observed), (qq_slope, qq_intercept, _) = probplot(std_residuals, dist=\"norm\")\n\n# LOWESS smoothers for panels 1 and 3\nlw1 = lowess(residuals, fitted_vals, frac=0.55, return_sorted=True)\nsqrt_abs_std = np.sqrt(np.abs(std_residuals))\nlw3 = lowess(sqrt_abs_std, fitted_vals, frac=0.55, return_sorted=True)\n\n# Panel names\nPANELS = [\"Residuals vs Fitted\", \"Normal Q-Q\", \"Scale-Location\", \"Residuals vs Leverage\"]\n\n# Main scatter data (long form)\nobs_ids = list(range(n))\ninfluence_flags = [\"high\" if i in top3_idx else \"normal\" for i in obs_ids]\n\ndf_main = pd.DataFrame(\n    {\n        \"x_val\": np.concatenate([fitted_vals, qq_theoretical, fitted_vals, leverage]),\n        \"y_val\": np.concatenate([residuals, qq_observed, sqrt_abs_std, std_residuals]),\n        \"panel\": pd.Categorical(\n            PANELS[0:1] * n + PANELS[1:2] * n + PANELS[2:3] * n + PANELS[3:4] * n, categories=PANELS, ordered=True\n        ),\n        \"obs_id\": obs_ids * 4,\n        \"influence\": pd.Categorical(influence_flags * 4, categories=[\"normal\", \"high\"]),\n    }\n)\n\n# LOWESS overlay — panels 1 and 3 only\ndf_lowess = pd.DataFrame(\n    {\n        \"x_val\": np.concatenate([lw1[:, 0], lw3[:, 0]]),\n        \"y_val\": np.concatenate([lw1[:, 1], lw3[:, 1]]),\n        \"panel\": pd.Categorical(PANELS[0:1] * len(lw1) + PANELS[2:3] * len(lw3), categories=PANELS, ordered=True),\n    }\n)\n\n# Q-Q reference line — panel 2 only\nqq_x_ends = np.array([qq_theoretical.min(), qq_theoretical.max()])\nqq_y_ends = qq_slope * qq_x_ends + qq_intercept\ndf_qq_ref = pd.DataFrame(\n    {\"x_val\": qq_x_ends, \"y_val\": qq_y_ends, \"panel\": pd.Categorical(PANELS[1:2] * 2, categories=PANELS, ordered=True)}\n)\n\n# Zero reference lines — panels 1 and 4\ndf_hline = pd.DataFrame(\n    {\"yintercept\": [0.0, 0.0], \"panel\": pd.Categorical([PANELS[0], PANELS[3]], categories=PANELS, ordered=True)}\n)\n\n# Cook's distance contours — panel 4\n# With high-leverage extreme-temperature observations, h can reach ~0.3,\n# making D=0.5 and D=1.0 contours visible within |std_resid| ≤ 3\nh_max_plot = leverage.max() * 1.3\nh_grid = np.linspace(1e-4, min(h_max_plot, 0.995), 500)\nsr_clip = max(np.abs(std_residuals).max() * 1.1, 3.5)\ncook_segs = []\nfor gid, (level, sign) in enumerate([(0.5, 1), (0.5, -1), (1.0, 1), (1.0, -1)]):\n    sr = sign * np.sqrt(level * p_hat * (1 - h_grid) / h_grid)\n    mask = np.abs(sr) <= sr_clip\n    if mask.sum() > 1:\n        cook_segs.append(\n            pd.DataFrame(\n                {\n                    \"x_val\": h_grid[mask],\n                    \"y_val\": sr[mask],\n                    \"panel\": pd.Categorical([PANELS[3]] * mask.sum(), categories=PANELS, ordered=True),\n                    \"cook_group\": [gid] * mask.sum(),\n                }\n            )\n        )\n\nif cook_segs:\n    df_cook = pd.concat(cook_segs, ignore_index=True)\nelse:\n    df_cook = pd.DataFrame(\n        {\n            \"x_val\": pd.Series(dtype=float),\n            \"y_val\": pd.Series(dtype=float),\n            \"panel\": pd.Categorical([], categories=PANELS, ordered=True),\n            \"cook_group\": pd.Series(dtype=int),\n        }\n    )\n\n# Labels for top influential points (panels 1, 3, 4 — not Q-Q)\nlabel_panels = [PANELS[0], PANELS[2], PANELS[3]]\ndf_labels = df_main[(df_main[\"influence\"] == \"high\") & df_main[\"panel\"].isin(label_panels)].copy()\ndf_labels[\"label\"] = df_labels[\"obs_id\"].astype(str)\n\n# Theme\nanyplot_theme = theme(\n    figure_size=(16, 9),\n    plot_background=element_rect(fill=PAGE_BG, color=PAGE_BG),\n    panel_background=element_rect(fill=PAGE_BG),\n    panel_grid_major=element_line(color=INK, size=0.3, alpha=0.10),\n    panel_grid_minor=element_blank(),\n    panel_border=element_rect(color=INK_SOFT, fill=None),\n    axis_title=element_text(color=INK, size=16),\n    axis_text=element_text(color=INK_SOFT, size=13),\n    plot_title=element_text(color=INK, size=20),\n    strip_background=element_rect(fill=ELEVATED_BG, color=INK_SOFT),\n    strip_text=element_text(color=INK, size=15),\n    legend_position=\"none\",\n)\n\n# Plot\nplot = (\n    ggplot(df_main, aes(x=\"x_val\", y=\"y_val\"))\n    + geom_hline(data=df_hline, mapping=aes(yintercept=\"yintercept\"), color=INK_SOFT, linetype=\"dashed\", size=0.7)\n    + geom_line(data=df_lowess, mapping=aes(x=\"x_val\", y=\"y_val\"), color=ACCENT, size=1.0)\n    + geom_line(data=df_qq_ref, mapping=aes(x=\"x_val\", y=\"y_val\"), color=INK_SOFT, size=0.8, linetype=\"dashed\")\n    + geom_line(\n        data=df_cook,\n        mapping=aes(x=\"x_val\", y=\"y_val\", group=\"cook_group\"),\n        color=INK_SOFT,\n        size=0.6,\n        linetype=\"dotted\",\n        alpha=0.7,\n    )\n    + geom_point(aes(color=\"influence\"), size=2.5, alpha=0.7, stroke=0.3)\n    + scale_color_manual(values={\"normal\": BRAND, \"high\": ACCENT})\n    + geom_text(data=df_labels, mapping=aes(label=\"label\"), va=\"bottom\", size=9, color=INK_SOFT)\n    + facet_wrap(\"~panel\", scales=\"free\", ncol=2)\n    + labs(title=\"diagnostic-regression-panel · plotnine · anyplot.ai\", x=\"\", y=\"\")\n    + anyplot_theme\n)\n\n# Save\nplot.save(f\"plot-{THEME}.png\", dpi=300)\n"}