# How to interpret coefficients of logistic regression using natural cubic splines?

I'm trying to experiment with using natural cubic splines for logistic regression in R, and I'm having trouble interpreting and making use of the results.

For a variety of reasons, I need to be able to plot (and possibly further analyze) the results in Python. The problem is that while I can do this with standard logistic regression, I don't know how to do this for the case of logistic regression using (natural cubic) splines.

An example: doing logistic regression on the frequency of stellar bars in spiral galaxies, as a function of the total stellar mass of individual galaxies. (Here's a link to the data file, if you want to try fitting things yourself.)

If I do standard logistic regression with a quadratic logit function -- as I did with an earlier version of the dataset in Erwin 2018 [Github link to Python and R notebooks: https://github.com/perwin/s4g_barfractions], I get something like this in R (stripping away some of the output):

ff <- "barpresence_vs_logmstar_for_R_sp_w30_updateddist.txt"

logMstarFitSp_quad <- glm(bar ~ logmstar + I(logmstar^2),
family = binomial, logmstarBarSpTable)

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   -19.2096    26.0781  -0.737    0.461
logmstar        4.2945     5.3723   0.799    0.424
I(logmstar^2)  -0.2332     0.2760  -0.845    0.398
AIC: 828.4


I can easily plot the result in Python, since I know the "Coefficients" values correspond to coefficients in a quadratic equation: logit = beta_0 + beta_1*x + beta_2*x^2 where beta_0 = "(Intercept)", beta_1 = "logmstar", and beta_2 = "I(logmstar^2)", so the actual expression becomes logit = -19.2096 + 4.2945*x - 0.2332*(x**2) and then I can generate probability values with Python code like this: probability = 1.0 / (1.0 + np.exp(-logit))

Here's the result of asking R to plot the fit, via plot_model(logMstarFitSp_quad, type = "pred", terms = c("logmstar"))

and here's the result of doing the same in Python, which agrees nicely:

def logit_quad( x, params=[-19.2096, 4.2945, -0.2332] ):
return params[0] + params[1]*x + params[2]*(x**2)
xx = np.arange(8.5,11,0.01)
plt.plot(xx, 1.0 / (1.0 + np.exp(-logit_quad(xx)))


If on the other hand I try doing a simple natural-cubic-spline logistic regression in R (for simplicity, starting off with just one internal knot, so df=2) I get:

require("splines")
logMstarFitSp_spline2df <- glm(bar ~ ns(logmstar, 2),
family = binomial, logmstarBarSpTable)
summary(logMstarFitSp_spline2df)

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)        0.5049     0.2285   2.210   0.0271 *
ns(logmstar, 2)1  -0.2487     0.5061  -0.492   0.6231
ns(logmstar, 2)2  -0.7645     0.4276  -1.788   0.0738 .
AIC: 828.16


For what it's worth, I can figure out that the single internal knot is located at logmstar = 9.5575, via

str(attributes(summ$terms)$predvars)
language list(bar, ns(logmstar, knots = c(50% = 9.5575), Boundary.knots = c(8.759, 11.024), intercept = FALSE))


... and I have no idea how to interpret this, or what the logit expression would look like in Python. What are the "Coefficients" in this model?

In principle, I think they should be the coefficients that multiply the basis functions of the spline, so that

   logit = beta_0*h_1(x) + beta_1*h_2(x) + beta_1*h_3(x)


where h_1, h_2, and h_3 are the three basis functions (and beta_0 = "(Intercept)", etc.).

But what are the basis functions? (My attempts to read the documentation for the R function ns have left me baffled; and I have not found much help in various spline-requested questions on this site.)

• If you're willing to use rms::Glm instead of stats::glm, there is a much easier way to obtain the equations, by using Function, as described here Commented Aug 14 at 16:21
• Btw I provide context here regarding what the basis functions and coefficients from splines actually mean @PeterErwin Commented Aug 14 at 17:33
• Two more possibilities: (1) compute the basis function yourself (e.g. using splines::ns(...) outside of your logistic fit), and export that as a CSV file (maybe not what you want); (2) generate equivalent spline bases in Python, e.g. using the patsy library as in this answer (it may take a bit of fussing to get all the arguments/defaults to match up with R's); (3) use a truncated power basis (as in rms::rcs or rms::rcspline.eval), which is not as numerically nice but is analytically simpler Commented Aug 14 at 17:39
• If you want the formula, then you will definitely save yourself headaches by using the truncated power basis rather than the B-spline basis (both can be formulated in a 'natural' way, i.e. with derivatives >1 going to zero at the boundaries, but the truncated power is much simpler) Commented Aug 14 at 20:47

If you're willing to use rms::Glm instead of stats::glm, there is a much easier way to obtain the equations, by using Function, as described here.

If you do, there is another useful post here that explains the difference between restricted cubic splines fit with rms::rcs and splines::ns, and a post here with general guidance on understanding spline basis functions..

Conversely, if you don't care much for the mathematical expression and just want to be able to plot, just create a sequence of new $$x$$-values that span the range of the data, and then use predict:

x_seq  <- seq(8.5, 11, l = 1000)
y_pred <- predict(logMstarFitSp_spline2df,
newdata = data.frame(logmstar = x_seq), type = "response")


With a bit more work, you can also obtain the confidence band, by using predict(..., se.fit = TRUE), then using a normal approximation on the linear scale, then converting back to the original scale:

x_seq <- seq(8.5, 11, l = 1000)
pred <- predict(GLMlogMstarFitSp_spline2dfnewdata = data.frame(logmstar = x_seq), se.fit = TRUE)

# Compute a 95% CI, using a normal approximation on the linear scale
lwr  <- pred$$fit + qnorm(0.025) * pred$$se.fit
upr  <- pred$$fit + qnorm(0.975) * pred$$se.fit

# Convert back to the original scale
il <- binomial()$$linkinv newy <- cbind(il(pred$$fit), il(lwr), il(upr))
plot(bar ~ logmstar, data = logmstarBarSpTable, pch = 16, bty = "n",
xlab = "logmstar", ylab = bquote(P["bar"]))
lines(newy[, 1] ~ newx, lwd = 3, col = "red")
lines(newy[, 2] ~ newx, lwd = 1, col = "red", lty = "dashed")
lines(newy[, 3] ~ newx, lwd = 1, col = "red", lty = "dashed")

• (+1) It might not be clear from the linked pages that ns() and rcs() differ not only in default knot positions but also in the choice of spline basis functions. Frank Harrell explains the choice for rcs() in Section 2.4.5 of Regression Modeling Strategies. That choice will provide formulas that are perhaps easier to associate with the modeled coefficients.
– EdM
Commented Aug 14 at 17:49