Restricted Cubic Spline

n_knots = 3
model = smf.glm(
    formula = f'mortality_30d ~ bs(PAR, df={n_knots}, include_intercept=False)',
    data=df,
    family=sm.families.Binomial(),
    #exposure=df['duration_30d']
)
res = model.fit()

model_linear = smf.glm(formula='mortality_30d ~ PAR', data = df, family=sm.families.Binomial())
res_linear = model_linear.fit()

LR = 2 * (res.llf - res_linear.llf)
df_diff = res.df_model - res_linear.df_model
p_nonlinear = stats.chi2.sf(LR, df_diff)

x_pred = np.linspace(df["PAR"].min(), df["PAR"].max(), 300)
df_pred = pd.DataFrame(data = {'PAR': x_pred})
pred = res.get_prediction(df_pred)
pred_mean = pred.predicted_mean
ci = pred.conf_int()

fig, ax = plt.subplots(figsize=(16,8))
ax.plot(x_pred, pred_mean, label="Fitted", color='red')
ax.fill_between(x_pred, ci[:,0], ci[:,1], color='red', alpha=0.2)
ax.axhline(y=0.5, color='black', linestyle='--', linewidth=0.8)
ax.axvline(x=df["PAR"].median(), color='gray', linestyle='--', linewidth=0.7)
# Add histogram
ax2 = ax.twinx()
sns.histplot(df["PAR"], ax=ax2, bins=40, color='blue', alpha=0.3, stat="density")
ax2.set_yticks([])
# Set labels
ax.set_xlabel("PAR")
ax.set_ylabel(f"Mortality 30 days probability")
ax.set_title(f"Logistic Regression 30 days")
text = f"P-nonlinear {format_p(p_nonlinear)}"
ax.text(0.7, 0.88, text, transform=ax.transAxes, fontsize=10, verticalalignment='top', horizontalalignment='right')
plt.show()