`as.matrix.projection()` for GAMs: Duplicated linear term

[fweber144]

On branch develop (commit 96d8e11), this reprex:

options(mc.cores = parallel::detectCores(logical = FALSE))
data("kidiq", package = "rstanarm")
kidiq <- head(kidiq, 45)
fit_gauss <- rstanarm::stan_gamm4(
  kid_score ~ s(mom_iq),
  data = kidiq,
  iter = 500,
  seed = 734572
)

library(projpred)
prj <- project(fit_gauss, solution_terms = c("mom_iq", "s(mom_iq)"),
               ndraws = 25)
colnames(as.matrix(prj))

gives

 [1] "b_Intercept"   "b_mom_iq"      "b_s(mom_iq).1" "b_s(mom_iq).2"
 [5] "b_s(mom_iq).3" "b_s(mom_iq).4" "b_s(mom_iq).5" "b_s(mom_iq).6"
 [9] "b_s(mom_iq).7" "b_s(mom_iq).8" "b_s(mom_iq).9" "sigma"

(so the output matrix includes a linear term for the smoothed predictor) whereas

colnames(as.matrix(fit_gauss))

(correctly) gives

 [1] "(Intercept)"           "s(mom_iq).1"          
 [3] "s(mom_iq).2"           "s(mom_iq).3"          
 [5] "s(mom_iq).4"           "s(mom_iq).5"          
 [7] "s(mom_iq).6"           "s(mom_iq).7"          
 [9] "s(mom_iq).8"           "s(mom_iq).9"          
[11] "sigma"                 "smooth_sd[s(mom_iq)1]"
[13] "smooth_sd[s(mom_iq)2]"

Perhaps the problem is not really with as.matrix.projection(), but rather with project() which would need to project onto the submodel with the smoothed term only (and not an additional linear term), even if solution_terms = c("mom_iq", "s(mom_iq)")?

There's another issue with as.matrix.projection() for GAMs visible here (the smooth_sd parameters) which I'll open a new issue for.

[fweber144]

There's another issue with as.matrix.projection() for GAMs visible here (the smooth_sd parameters) which I'll open a new issue for.

This is issue #151 now.

[AlejandroCatalina]

This has to do with the way we split the original fit object's formula to do the forward search. We defined some conventions that seemed reasonable but they may give raise to stuff like this. In this case, for instance, it was decided that we also want to project the main effect of the smoothing term, which serves as a way of asking "do I need the smoothing term or is the linear term enough in terms of predictive accuracy?".

[fweber144]

In this case, for instance, it was decided that we also want to project the main effect of the smoothing term, which serves as a way of asking "do I need the smoothing term or is the linear term enough in terms of predictive accuracy?".

That definitely makes sense, but I think we need to differentiate between the (forward) search and the projection onto a pre-specified submodel. My concern here is the latter.

In any case, a submodel shouldn't contain both, the nonsmoothed term (e.g. mom_iq) as well as the smoothed term (e.g. s(mom_iq)), right? For example, if you fit a model with both of these terms in rstanarm, you'll get a lot of convergence issues. So including both terms leads to a degenerate model. Therefore, I think projpred should catch that case (when projecting onto a pre-specified submodel) and either

• throw an error or
• exclude the nonsmoothed term and throw a warning.

EDIT: An earlier version referred to mgcv where fitting a model with a nonsmoothed as well as a smoothed term would give 0 reference degrees of freedom for the smoothed term and an F-statistic as well as a p-value which are NaN. However, that doesn't seem to be necessarily the case, see e.g.:

library(mgcv)
set.seed(2)
dat <- gamSim(1, n = 40)
b <- gam(y ~ x2 + s(x2), data = dat)
summary(b)

When debugging the projection of the first draw in the reprex above, however, we do get 0 reference degrees of freedom for the smoothed term and an F-statistic as well as a p-value which are NaN.

[AlejandroCatalina]

As far as I remember if the smooth is included first the linear term can’t be added, but if the linear term is added first the smooth can be included later. So a model with both could exist and from the modeling perspective it could make sense in the right model, but I agree they have to be properly constrained and so on. Probably @ Get Outlook for iOS<https://aka.ms/o0ukef>

…

________________________________ From: Frank Weber ***@***.***> Sent: Monday, July 5, 2021 6:46:39 PM To: stan-dev/projpred ***@***.***> Cc: Alejandro Catalina ***@***.***>; Comment ***@***.***> Subject: Re: [stan-dev/projpred] `as.matrix.projection()` for GAMs: Duplicated linear term (#150) In this case, for instance, it was decided that we also want to project the main effect of the smoothing term, which serves as a way of asking "do I need the smoothing term or is the linear term enough in terms of predictive accuracy?". That definitely makes sense, but I think we need to differentiate between the (forward) search and the projection onto a pre-specified submodel. My concern here is the latter. In any case, a submodel shouldn't contain both, the nonsmoothed term (e.g. mom_iq) as well as the smoothed term (e.g. s(mom_iq)), right? For example, if you fit a model with both of these terms in rstanarm, you'll get a lot of convergence issues. And if you fit such a model in mgcv, you'll get 0 reference degrees of freedom (Ref.df in summary.gam()'s output) for the smoothed term and an F-statistic as well as a p-value which are NaN. So including both terms leads to a degenerate model. Therefore, I think projpred should catch that case (when projecting onto a pre-specified submodel) and either * throw an error or * exclude the nonsmoothed term and throw a warning. — You are receiving this because you commented. Reply to this email directly, view it on GitHub<#150 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABZ5FH3COG6A7T6SVMV5T4TTWHHV7ANCNFSM46Q6WWEQ>.

[fweber144]

The models I have tested with rstanarm and mgcv had the linear term included first, so that order didn't help.