import pandas as pd
import numpy as np
import scipy.stats as spicy
import plotnine as pn
import statsmodels.api as smLinear Regression
df = pd.read_csv("https://raw.githubusercontent.com/roualdes/data/refs/heads/master/penguins.csv")Unique Intercepts and Unique Slopes
fit = sm.OLS.from_formula("body_mass_g ~ sex * flipper_length_mm", data = df).fit()
fit.summary()| Dep. Variable: | body_mass_g | R-squared: | 0.806 |
| Model: | OLS | Adj. R-squared: | 0.804 |
| Method: | Least Squares | F-statistic: | 455.2 |
| Date: | Tue, 14 Apr 2026 | Prob (F-statistic): | 1.04e-116 |
| Time: | 12:38:33 | Log-Likelihood: | -2427.2 |
| No. Observations: | 333 | AIC: | 4862. |
| Df Residuals: | 329 | BIC: | 4878. |
| Df Model: | 3 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
| Intercept | -5443.9607 | 440.283 | -12.365 | 0.000 | -6310.086 | -4577.836 |
| sex[T.male] | 406.8015 | 587.303 | 0.693 | 0.489 | -748.541 | 1562.144 |
| flipper_length_mm | 47.1527 | 2.226 | 21.179 | 0.000 | 42.773 | 51.532 |
| sex[T.male]:flipper_length_mm | -0.2942 | 2.924 | -0.101 | 0.920 | -6.047 | 5.458 |
| Omnibus: | 0.262 | Durbin-Watson: | 1.711 |
| Prob(Omnibus): | 0.877 | Jarque-Bera (JB): | 0.375 |
| Skew: | 0.051 | Prob(JB): | 0.829 |
| Kurtosis: | 2.871 | Cond. No. | 8.24e+03 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 8.24e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
Two qualitative Explanatory Variables
fit = sm.OLS.from_formula("body_mass_g ~ sex + species * flipper_length_mm", data = df).fit()
fit.summary()| Dep. Variable: | body_mass_g | R-squared: | 0.869 |
| Model: | OLS | Adj. R-squared: | 0.867 |
| Method: | Least Squares | F-statistic: | 360.2 |
| Date: | Tue, 14 Apr 2026 | Prob (F-statistic): | 1.44e-140 |
| Time: | 12:38:33 | Log-Likelihood: | -2361.8 |
| No. Observations: | 333 | AIC: | 4738. |
| Df Residuals: | 326 | BIC: | 4764. |
| Df Model: | 6 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
| Intercept | -60.0183 | 734.674 | -0.082 | 0.935 | -1505.319 | 1385.283 |
| sex[T.male] | 521.5581 | 38.145 | 13.673 | 0.000 | 446.518 | 596.599 |
| species[T.Chinstrap] | 927.8374 | 1222.338 | 0.759 | 0.448 | -1476.828 | 3332.503 |
| species[T.Gentoo] | -1066.4650 | 1165.273 | -0.915 | 0.361 | -3358.869 | 1225.939 |
| flipper_length_mm | 18.4395 | 3.888 | 4.742 | 0.000 | 10.790 | 26.089 |
| species[T.Chinstrap]:flipper_length_mm | -5.1393 | 6.300 | -0.816 | 0.415 | -17.533 | 7.254 |
| species[T.Gentoo]:flipper_length_mm | 8.9573 | 5.637 | 1.589 | 0.113 | -2.133 | 20.048 |
| Omnibus: | 1.300 | Durbin-Watson: | 2.083 |
| Prob(Omnibus): | 0.522 | Jarque-Bera (JB): | 1.377 |
| Skew: | 0.110 | Prob(JB): | 0.502 |
| Kurtosis: | 2.774 | Cond. No. | 2.08e+04 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.08e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
ndf = pd.DataFrame({"sex": ["female", "male"],
"species": ["Adelie", "Adelie"],
"flipper_length_mm": [200, 200]})
ndf| sex | species | flipper_length_mm | |
|---|---|---|---|
| 0 | female | Adelie | 200 |
| 1 | male | Adelie | 200 |
fit.predict(ndf)0 3627.886070
1 4149.444186
dtype: float64
Non-linear predictors
fit = sm.OLS.from_formula("body_mass_g ~ bill_depth_mm * flipper_length_mm", data = df).fit()
fit.summary()| Dep. Variable: | body_mass_g | R-squared: | 0.787 |
| Model: | OLS | Adj. R-squared: | 0.785 |
| Method: | Least Squares | F-statistic: | 416.2 |
| Date: | Tue, 14 Apr 2026 | Prob (F-statistic): | 4.20e-113 |
| Time: | 12:38:33 | Log-Likelihood: | -2507.3 |
| No. Observations: | 342 | AIC: | 5023. |
| Df Residuals: | 338 | BIC: | 5038. |
| Df Model: | 3 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
| Intercept | -3.61e+04 | 4636.271 | -7.786 | 0.000 | -4.52e+04 | -2.7e+04 |
| bill_depth_mm | 1771.7958 | 273.003 | 6.490 | 0.000 | 1234.796 | 2308.795 |
| flipper_length_mm | 196.0737 | 22.603 | 8.675 | 0.000 | 151.613 | 240.534 |
| bill_depth_mm:flipper_length_mm | -8.5964 | 1.340 | -6.414 | 0.000 | -11.233 | -5.960 |
| Omnibus: | 8.083 | Durbin-Watson: | 2.205 |
| Prob(Omnibus): | 0.018 | Jarque-Bera (JB): | 8.056 |
| Skew: | 0.372 | Prob(JB): | 0.0178 |
| Kurtosis: | 3.107 | Cond. No. | 7.97e+05 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 7.97e+05. This might indicate that there are
strong multicollinearity or other numerical problems.
fit = sm.OLS.from_formula("body_mass_g ~ flipper_length_mm + I(flipper_length_mm ** 2)", data = df).fit()
fit.summary()| Dep. Variable: | body_mass_g | R-squared: | 0.775 |
| Model: | OLS | Adj. R-squared: | 0.774 |
| Method: | Least Squares | F-statistic: | 584.8 |
| Date: | Tue, 14 Apr 2026 | Prob (F-statistic): | 1.27e-110 |
| Time: | 12:38:33 | Log-Likelihood: | -2516.5 |
| No. Observations: | 342 | AIC: | 5039. |
| Df Residuals: | 339 | BIC: | 5050. |
| Df Model: | 2 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
| Intercept | 1.677e+04 | 4559.552 | 3.678 | 0.000 | 7802.304 | 2.57e+04 |
| flipper_length_mm | -173.2311 | 44.999 | -3.850 | 0.000 | -261.743 | -84.719 |
| I(flipper_length_mm ** 2) | 0.5482 | 0.111 | 4.956 | 0.000 | 0.331 | 0.766 |
| Omnibus: | 10.640 | Durbin-Watson: | 2.281 |
| Prob(Omnibus): | 0.005 | Jarque-Bera (JB): | 10.849 |
| Skew: | 0.432 | Prob(JB): | 0.00441 |
| Kurtosis: | 3.123 | Cond. No. | 9.06e+06 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 9.06e+06. This might indicate that there are
strong multicollinearity or other numerical problems.
ndf = pd.DataFrame({"flipper_length_mm": np.linspace(170, 230, 101)})
ndf["yhat"] = fit.predict(ndf)
pn.ggplot() + \
pn.geom_point(dfn, pn.aes("flipper_length_mm", "body_mass_g"), alpha = 0.1) + \
pn.geom_line(ndf, pn.aes("flipper_length_mm", "yhat"), size = 2)Polynomials of high order are dangerous
rng = np.random.default_rng()
x = np.arange(4)
y = 3 + 5 * x + rng.normal(size = 4) * 5
df = pd.DataFrame({"x": x, "y": y})
fit = sm.OLS.from_formula("y ~ x + I(x * x) + I(x * x * x)", data = df).fit()
df['yhat'] = fit.predict()
pn.ggplot() + \
pn.geom_point(df, pn.aes("x", "y")) + \
pn.geom_line(df, pn.aes("x", "yhat"))Nonlinearly transforming the response
df = pd.read_csv("https://raw.githubusercontent.com/roualdes/data/refs/heads/master/carnivora.csv")
df["brain_weight"] = df["SB"]
df["body_weight"] = df["SW"]pn.ggplot() + \
pn.geom_point(df, pn.aes("brain_weight", "body_weight")) + \
pn.labs("Brain weight", "Body Weight")df["log10_body_weight"] = np.log10(df["body_weight"])
df["log10_brain_weight"] = np.log10(df["brain_weight"])
fit = sm.OLS.from_formula("np.log10(body_weight) ~ np.log10(brain_weight)", data = df).fit()
fit.summary()| Dep. Variable: | np.log10(body_weight) | R-squared: | 0.919 |
| Model: | OLS | Adj. R-squared: | 0.918 |
| Method: | Least Squares | F-statistic: | 1244. |
| Date: | Tue, 14 Apr 2026 | Prob (F-statistic): | 8.70e-62 |
| Time: | 13:24:36 | Log-Likelihood: | 17.769 |
| No. Observations: | 112 | AIC: | -31.54 |
| Df Residuals: | 110 | BIC: | -26.10 |
| Df Model: | 1 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
| Intercept | -1.4394 | 0.062 | -23.076 | 0.000 | -1.563 | -1.316 |
| np.log10(brain_weight) | 1.4080 | 0.040 | 35.268 | 0.000 | 1.329 | 1.487 |
| Omnibus: | 27.817 | Durbin-Watson: | 1.650 |
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 59.485 |
| Skew: | 0.974 | Prob(JB): | 1.21e-13 |
| Kurtosis: | 5.992 | Cond. No. | 6.83 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
x = np.linspace(np.min(df["brain_weight"]) , np.max(df["brain_weight"]), 101)
ndf = pd.DataFrame({"brain_weight": x})
ndf["yhat"] = 10 ** fit.predict(ndf)pn.ggplot() + \
pn.geom_point(df, pn.aes("brain_weight", "body_weight")) + \
pn.geom_line(ndf, pn.aes("brain_weight", "yhat")) + \
pn.labs("Brain weight", "Body Weight")\(R^2\)
pn.ggplot(df, pn.aes("brain_weight", "body_weight")) + \
pn.geom_point() + \
pn.geom_smooth(method = "lm") + \
pn.labs("Brain weight", "Body Weight")spicy.stats.pearsonr(df["body_weight"], df["brain_weight"])[0]** 2 # correlation/var/folders/md/4bj70ky11855zs8hgz7gb_0w0000gn/T/ipykernel_35749/2097628686.py:1: DeprecationWarning: Please import `pearsonr` from the `scipy.stats` namespace; the `scipy.stats.stats` namespace is deprecated and will be removed in SciPy 2.0.0.
0.8540655598474169
fit = sm.OLS.from_formula("body_weight ~ brain_weight", data = df).fit()
fit.summary()| Dep. Variable: | body_weight | R-squared: | 0.854 |
| Model: | OLS | Adj. R-squared: | 0.853 |
| Method: | Least Squares | F-statistic: | 643.8 |
| Date: | Tue, 14 Apr 2026 | Prob (F-statistic): | 8.76e-48 |
| Time: | 13:47:05 | Log-Likelihood: | -492.71 |
| No. Observations: | 112 | AIC: | 989.4 |
| Df Residuals: | 110 | BIC: | 994.9 |
| Df Model: | 1 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
| Intercept | -15.9456 | 2.354 | -6.773 | 0.000 | -20.611 | -11.280 |
| brain_weight | 0.6385 | 0.025 | 25.372 | 0.000 | 0.589 | 0.688 |
| Omnibus: | 43.983 | Durbin-Watson: | 2.268 |
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 429.253 |
| Skew: | 0.894 | Prob(JB): | 6.15e-94 |
| Kurtosis: | 12.423 | Cond. No. | 117. |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.