Linear Regression

import pandas as pd
import numpy as np
import scipy.stats as spicy

import plotnine as pn
import statsmodels.api as sm
df = pd.read_csv("https://raw.githubusercontent.com/roualdes/data/refs/heads/master/penguins.csv")

Unique Intercepts and a shared Slope

fit = sm.OLS.from_formula("body_mass_g ~ sex + flipper_length_mm", data = df).fit()
fit.summary()
OLS Regression Results
Dep. Variable: body_mass_g R-squared: 0.806
Model: OLS Adj. R-squared: 0.805
Method: Least Squares F-statistic: 684.8
Date: Tue, 14 Apr 2026 Prob (F-statistic): 3.53e-118
Time: 12:38:33 Log-Likelihood: -2427.2
No. Observations: 333 AIC: 4860.
Df Residuals: 330 BIC: 4872.
Df Model: 2
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -5410.3002 285.798 -18.931 0.000 -5972.515 -4848.085
sex[T.male] 347.8503 40.342 8.623 0.000 268.491 427.209
flipper_length_mm 46.9822 1.441 32.598 0.000 44.147 49.817
Omnibus: 0.262 Durbin-Watson: 1.710
Prob(Omnibus): 0.877 Jarque-Bera (JB): 0.376
Skew: 0.051 Prob(JB): 0.829
Kurtosis: 2.870 Cond. No. 2.95e+03


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.95e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

For this model the coefficients would be named as

\[\beta_0, \beta_1, \beta_2\]

where \(\beta_2\) refers to the slope on flipper length.

dfn = df.dropna(subset = ["sex", "body_mass_g"])

pn.ggplot(dfn) + \
    pn.geom_point(pn.aes("flipper_length_mm", "body_mass_g"))  + \
    pn.stat_smooth(pn.aes("flipper_length_mm", "body_mass_g"), method = "lm", se = False)

Unique Intercepts and Unique Slopes

fit = sm.OLS.from_formula("body_mass_g ~ sex * flipper_length_mm", data = df).fit()
fit.summary()
OLS Regression Results
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()
OLS Regression Results
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()
OLS Regression Results
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()
OLS Regression Results
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()
OLS Regression Results
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()
OLS Regression Results
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.