Logistic Regression

import numpy as np
import plotnine as pn
import pandas as pd
import statsmodels.api as sm
from sklearn.metrics import confusion_matrix
# carnivora
df = pd.read_csv("https://raw.githubusercontent.com/roualdes/data/refs/heads/master/carnivora.csv")
df = df.dropna(subset = ["SuperFamily", "LY"])
df["dog"] = (df["SuperFamily"] == "Caniformia").astype(np.int64)
fit = sm.GLM.from_formula("dog ~ LY", 
                          data = df, 
                          family=sm.families.Binomial() ).fit()
fit.summary()
Generalized Linear Model Regression Results
Dep. Variable: dog No. Observations: 49
Model: GLM Df Residuals: 47
Model Family: Binomial Df Model: 1
Link Function: Logit Scale: 1.0000
Method: IRLS Log-Likelihood: -33.400
Date: Tue, 28 Apr 2026 Deviance: 66.799
Time: 13:11:54 Pearson chi2: 48.8
No. Iterations: 4 Pseudo R-squ. (CS): 0.02237
Covariance Type: nonrobust
coef std err z P>|z| [0.025 0.975]
Intercept -0.8570 0.842 -1.017 0.309 -2.508 0.794
LY 0.0045 0.004 1.025 0.306 -0.004 0.013
x = np.linspace(np.min(df["LY"]), np.max(df["LY"]), 101)
ndf = pd.DataFrame({"LY": x})
ndf["phat"] = fit.predict(ndf)
pn.ggplot() + \
    pn.geom_point(df, pn.aes("LY", "dog")) + \
    pn.geom_line(ndf, pn.aes("LY", "phat"))

df["p_dog"] = fit.predict() > 0.5
#np.mean((df["dog"] == 1) & (df["p_dog"] == 1)) # True Positive
#np.mean((df["dog"] == 1) & (df["p_dog"] == 0))  # False Negative
np.mean((df["dog"] == 0) & (df["p_dog"] == 1)) # False Positive
#np.mean((df["dog"] == 0) & (df["p_dog"] == 0)) # True Negative
0.14285714285714285
confusion_matrix(df["dog"], df["p_dog"], normalize = "all")
array([[0.36734694, 0.14285714],
       [0.34693878, 0.14285714]])
Predicted 0 Predicted 1
Actual 0 TN 😀 FP 🙁
Actual 1 FN 🙁 TP 😀
df = pd.read_csv("http://roualdes.sfo3.digitaloceanspaces.com/data/maize.csv")
df = df.dropna(subset = ["yield", "plantheight"])
df["lots"] = (df["yield"] > np.median(df["yield"])).astype(np.int64)
pn.ggplot() + \
    pn.geom_point(df, pn.aes("plantheight", "lots"))

fit = sm.GLM.from_formula("lots ~ plantheight", 
                          data = df, 
                          family=sm.families.Binomial() ).fit()
df["p_lots"] = fit.predict()
pn.ggplot() + \
    pn.geom_point(df, pn.aes("plantheight", "lots")) + \
    pn.geom_line(df, pn.aes("plantheight", "p_lots"))

def bootstrap(arr, T, R = 1_000):
    N = np.shape(arr)[0]
    Ts = np.zeros(R)
    rng = np.random.default_rng()
    for r in range(R):
        idx = rng.integers(N, size = N)
        if type(arr) is np.ndarray:
            Ts[r] = T(arr[idx])
        else:
            Ts[r] = T(arr.iloc[idx])
    return Ts
m = np.mean(df["plantheight"])
def logistic_ci(_data):
    fit = sm.GLM.from_formula("lots ~ plantheight", 
                          data = _data, 
                          family=sm.families.Binomial() ).fit()
    ndf = pd.DataFrame({"plantheight": [m, m + 1]})
    return np.diff(fit.predict(ndf))[0]
b = bootstrap(df, logistic_ci)
np.quantile(b, [0.025, 0.975])
array([0.00948544, 0.01030603])