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

import plotnine as pn
import statsmodels.api as sm
# same dataset from hypothesis tests, take 01
# penguins
df = pd.read_csv("https://raw.githubusercontent.com/roualdes/data/refs/heads/master/penguins.csv")
dfn = df.dropna(subset = ["sex", "body_mass_g"])
sdx = dfn["sex"] == "male"
m = dfn.loc[sdx, "body_mass_g"]
f = dfn.loc[~sdx, "body_mass_g"]
spicy.ttest_ind(m, f) # Two sample t-test
TtestResult(statistic=8.541720337994516, pvalue=4.897246751596224e-16, df=331.0)
pn.ggplot(dfn) + \
    pn.geom_jitter(pn.aes("sex", "body_mass_g", color = "sex"), alpha = 0.25) + \
    pn.stat_summary(pn.aes("sex", "body_mass_g"), fun_data = "mean_cl_boot") + \
    pn.labs(y = "Body mass (g)") + \
    pn.theme_minimal()

# ordinary least squares
fit = sm.OLS.from_formula("body_mass_g ~ sex", data = df).fit()
fit.summary()
OLS Regression Results
Dep. Variable: body_mass_g R-squared: 0.181
Model: OLS Adj. R-squared: 0.178
Method: Least Squares F-statistic: 72.96
Date: Mon, 06 Apr 2026 Prob (F-statistic): 4.90e-16
Time: 13:39:29 Log-Likelihood: -2667.0
No. Observations: 333 AIC: 5338.
Df Residuals: 331 BIC: 5346.
Df Model: 1
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 3862.2727 56.829 67.963 0.000 3750.481 3974.064
sex[T.male] 683.4118 80.009 8.542 0.000 526.022 840.801
Omnibus: 128.992 Durbin-Watson: 0.410
Prob(Omnibus): 0.000 Jarque-Bera (JB): 26.595
Skew: 0.402 Prob(JB): 1.68e-06
Kurtosis: 1.873 Cond. No. 2.63


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
3862.2727 + 683.4118
4545.6845
dfn.groupby("sex").aggregate(mn = ("body_mass_g", np.mean))
mn
sex
female 3862.272727
male 4545.684524
# Analysis of Variance output
sm.stats.anova_lm(fit)
df sum_sq mean_sq F PR(>F)
sex 1.0 3.887890e+07 3.887890e+07 72.960986 4.897247e-16
Residual 331.0 1.763808e+08 5.328724e+05 NaN NaN
8.541720337994516 ** 2
72.96098633250914
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)

fit = sm.OLS.from_formula("body_mass_g ~ flipper_length_mm", data = df).fit()
fit.summary()
OLS Regression Results
Dep. Variable: body_mass_g R-squared: 0.759
Model: OLS Adj. R-squared: 0.758
Method: Least Squares F-statistic: 1071.
Date: Mon, 06 Apr 2026 Prob (F-statistic): 4.37e-107
Time: 13:39:29 Log-Likelihood: -2528.4
No. Observations: 342 AIC: 5061.
Df Residuals: 340 BIC: 5069.
Df Model: 1
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -5780.8314 305.815 -18.903 0.000 -6382.358 -5179.305
flipper_length_mm 49.6856 1.518 32.722 0.000 46.699 52.672
Omnibus: 5.634 Durbin-Watson: 2.190
Prob(Omnibus): 0.060 Jarque-Bera (JB): 5.585
Skew: 0.313 Prob(JB): 0.0613
Kurtosis: 3.019 Cond. No. 2.89e+03


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.89e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
-5780.8314 + 49.6856 * 200
4156.288600000001
fit = sm.OLS.from_formula("body_mass_g ~ island", data = df).fit()
fit.summary()
OLS Regression Results
Dep. Variable: body_mass_g R-squared: 0.394
Model: OLS Adj. R-squared: 0.390
Method: Least Squares F-statistic: 110.0
Date: Mon, 06 Apr 2026 Prob (F-statistic): 1.52e-37
Time: 13:40:02 Log-Likelihood: -2686.2
No. Observations: 342 AIC: 5378.
Df Residuals: 339 BIC: 5390.
Df Model: 2
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 4716.0180 48.468 97.301 0.000 4620.682 4811.354
island[T.Dream] -1003.1147 74.249 -13.510 0.000 -1149.162 -857.067
island[T.Torgersen] -1009.6454 100.207 -10.076 0.000 -1206.752 -812.539
Omnibus: 9.164 Durbin-Watson: 1.683
Prob(Omnibus): 0.010 Jarque-Bera (JB): 9.126
Skew: -0.370 Prob(JB): 0.0104
Kurtosis: 3.303 Cond. No. 3.52


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.