Extracting knot points from glm when using bs in R as a variable

Question

I fit a model using the following:

mymodel <- glm(LS ~ bs(LA, df = 8) + bs(IN, df = 7),
               family = binomial, data = mydata, na.action = na.omit)

No problem, I have the model fit now I am trying to extract the knot points used. I followed a post on extracting knot points using attr and str. That was for a model that was just a spline. I think the knots are somewhere in the structure in terms

I called str(mymodel$terms) there are ..-attr(*, "variables"). I am having trouble going further with attr but I am relatively certain that this is basically what I need to do. Any guidance to get the knots is appreciated.

Answer 1

You can use

eval(attr(mymodel$terms, "predvars"))

which evaluates the language object contained in the predvars attribute of the terms component of the fitted model.

Here is an example with a silly fitted model

mod <- glm(rnorm(length(women$height)) ~ bs(women$height, df = 5))

from which we can evaluate the required part of the terms component of mod

> eval(attr(mod$terms, "predvars"))
[[1]]
 [1] -1.20088330 -0.46267556 -0.04791518 -1.42748340  2.32896914  0.07858849
 [7]  2.16635328 -0.78670562 -1.68737883  0.71389437 -0.64123154 -0.04891306
[13] -0.07260125  0.71263717 -2.63426761

[[2]]
                 1           2           3            4           5
 [1,] 0.000000e+00 0.000000000 0.000000000 0.000000e+00 0.000000000
 [2,] 4.534439e-01 0.059857872 0.001639942 0.000000e+00 0.000000000
 [3,] 5.969388e-01 0.203352770 0.013119534 0.000000e+00 0.000000000
 [4,] 5.338010e-01 0.376366618 0.044278426 0.000000e+00 0.000000000
 [5,] 3.673469e-01 0.524781341 0.104956268 0.000000e+00 0.000000000
 [6,] 2.001640e-01 0.595025510 0.204719388 9.110787e-05 0.000000000
 [7,] 9.110787e-02 0.566326531 0.336734694 5.830904e-03 0.000000000
 [8,] 3.125000e-02 0.468750000 0.468750000 3.125000e-02 0.000000000
 [9,] 5.830904e-03 0.336734694 0.566326531 9.110787e-02 0.000000000
[10,] 9.110787e-05 0.204719388 0.595025510 2.001640e-01 0.000000000
[11,] 0.000000e+00 0.104956268 0.524781341 3.673469e-01 0.002915452
[12,] 0.000000e+00 0.044278426 0.376366618 5.338010e-01 0.045553936
[13,] 0.000000e+00 0.013119534 0.203352770 5.969388e-01 0.186588921
[14,] 0.000000e+00 0.001639942 0.059857872 4.534439e-01 0.485058309
[15,] 0.000000e+00 0.000000000 0.000000000 0.000000e+00 1.000000000
attr(,"degree")
[1] 3
attr(,"knots")
33.33333% 66.66667% 
 62.66667  67.33333 
attr(,"Boundary.knots")
[1] 58 72
attr(,"intercept")
[1] FALSE
attr(,"class")
[1] "bs"     "basis"  "matrix"

In the resulting list, the first and second components are the response and predictor data respectively. To this list a number of attributes are attached to the second component, the bs data. You need to extract those.

ll <- eval(attr(mod$terms, "predvars"))
attr(ll[[2]], "knots")
attr(ll[[2]], "Boundary.knots")

> attr(ll[[2]], "knots")
33.33333% 66.66667% 
 62.66667  67.33333 
> attr(ll[[2]], "Boundary.knots")
[1] 58 72