1
$\begingroup$

I'm learning about modeling B-Splines using PyMC3. The design matrix of splines (apparently) can become quite complicated, so it's easier to delegate this construction to an API, Patsy. In the context of B-Splines, I don't understand what the formula-like string is supposed to do.

The tutorial I'm following (see section 4.74 in https://github.com/pymc-devs/resources/blob/master/Rethinking_2/Chp_04.ipynb) uses Patsy's dmatrix in order to construct its design matrix. I don't quite "get it."

I've created some synthetic data, namely, columns, year and value. In the below synthetic data, I've used one knot per observation and degree of 0 (no observation is influenced by more than 1 spline) to hopefully keep things simple.

Synthetic data

fake = pd.DataFrame({'year':[1,2,3],'value':[1,1,1]})
k = [1,2,3]

Without intercept:

dmatrix(
    "bs(year, knots=knots, degree=0, include_intercept=False) ",
    {"year": fake.year.values, "knots": k},
).view()
>>>
array([[1., 1., 0., 0.],
       [1., 0., 1., 0.],
       [1., 0., 0., 1.]])

Without intercept -1:

dmatrix(
    "bs(year, knots=knots, degree=0, include_intercept=False) -1",
    {"year": fake.year.values, "knots": k},
).view()
>>>
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

With intercept:

dmatrix(
    "bs(year, knots=knots, degree=0, include_intercept=True) ",
    {"year": fake.year.values, "knots": k},
).view()
>>>
array([[1., 0., 1., 0., 0.],
       [1., 0., 0., 1., 0.],
       [1., 0., 0., 0., 1.]])

With intercept -1:

dmatrix(
    "bs(year, knots=knots, degree=0, include_intercept=True) -1",
    {"year": fake.year.values, "knots": k},
).view()
>>>
array([[0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])

API references: https://patsy.readthedocs.io/en/latest/spline-regression.html https://patsy.readthedocs.io/en/latest/API-reference.html

From what I can gather, the default behavior is the first column containing all 1's. By including -1 at the end of the "formula-like" string object, this column of ones if removed. Likewise, by default, the second column is all zeroes, and include_intercept=False removes this column of zeroes.

I'd like to know, for the context of spline regression: (A) what is the columns of zeroes (intercept?) for and why is it all zeroes? Likewise, (B) what is the column of all ones for? (I would have guessed that ones represented the intercept, but I guess this is wrong.)

$\endgroup$
1
$\begingroup$

The column that is giving you all zeros doesn't always give all zeros. I never use it, but I see that if I change year to include 0, the bs-intercept column has a one in it:

patsy.dmatrix(
    "bs(year, knots=knots, degree=0, include_intercept=True) ",
    {"year": [0,1,2,3], "knots": k}).view()

array([[1., 1., 0., 0., 0.],
       [1., 0., 1., 0., 0.],
       [1., 0., 0., 1., 0.],
       [1., 0., 0., 0., 1.]])

The first column is the more traditional "intercept" term in regression modeling, where if you didn't have the spline part making things confusing, you would be predicting $y \sim \beta_0 + \beta_1x_1 + \ldots \beta_p x_p$.

It helps me to think of bs as a fancy alternative to including nonlinear transformations of $x$---instead of $x^2$ and $x^3$ to try to represent a nonlinear response, I can use bs(x) to do a better job with less typing. For me a picture is worth a thousand words, so here is one that I find useful:

xx = np.linspace(-5,5,100)
df = patsy.dmatrix(
        "bs(xx, knots=[0,2], degree=3, include_intercept=False)",
        return_type='dataframe')
df.index = xx

df.plot()
plt.legend(loc=(1.01, .01));

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.