rmcelreath / rethinking

Statistical Rethinking course and book package
2.1k stars 596 forks source link

Smooth functions of a variable - modeling syntax #400

Open wesleyburr opened 1 year ago

wesleyburr commented 1 year ago

Chapter 4.5.2 of rethinking, 2nd edition - the use of splines. The model example in 4.5.2 uses *B %% w to represent the smooth function of some variable, pre-computed as a B-spline basis matrix B, with parameter vector w. All well and good. The example implemented on page 119 runs fine, using start = to initialize the w** vector as all 0s to start.

What is not clear from this, or at least was not sufficiently that it clicked for me to flag it for my students while teaching this to them, is that without initializing w, the model throws an unconformable objects error out of quap (specifically, from the fit call on line 531 of map-quap.R).

This probably deserves an Overthinking box to mention to the interested reader that without initializing w via start, the model parser has no idea what dimension w should be, and will not set it up properly. It might also be possible to simply force w to be the dimension of the column dimension of the matrix B, by detecting the %*% operator and then putting in special logic; or at least, this might need to be flagged as a more interpretable error message for the readers?

Context: I assigned my students an assignment consisting of a Gaussian-link GLM with a smooth function of one parameter, a categorical variable, and an additive linear term. Several students were not getting things to work, and it took me a few minutes to figure out that it was because they had not included the start initialization in their quap call.

rmcelreath commented 1 year ago

Yeah good point. It's an awkward problem either way, given that typed variables are not so natural to so many students. I could give your idea about %*% a try. Should we easy to add a check for uncormforbale objects while the linear model is parsed.

wesleyburr commented 1 year ago

Update on this: I asked the students to try to use ulam() this week, and ran into a similar issue. It seems like the way you parse the formulas into Stan code is a bit different for ulam() [seems = it is, as you clearly state in the docs, because of the more powerful support for GLMMs and other features].

In an ulam() call, setting the start values is insufficient to properly set up the vector size for the parameters for a smooth function. The workaround I found is to manually drop into Stan-like syntax and pre-pend the prior statements for those particular (vector) priors, e.g., vector[9]:w ~ rnorm(0, 1) instead of the more textbook-described w ~ rnorm(0, 1), where the parser takes care of it.

I have some students who seem to like to push the boundaries ....