7 min read · May 25, 2018
--
If you, like me, are learning how to handle generalised additive models, it may have crossed your mind to perform some automated variable selection, doing things like varying the degrees of freedom of particular splines and such in your model.
I would caution you that unless have you have been forced to do stepwise variable selection or are very careful about interpreting it, it is generally considered a bad idea, and you may want to consider alternatives. Here is a good place to start your search if so.
We will be using the gam package, which is a good place to start learning GAMs. If you are looking to use GAMs in a more mature capacity in R, you should investigate the mgcv package instead, which is favored in the community.
Disclaimer: do not use the mgcv package and gam at the same time. Their namespaces are not compatible and overlap. This will lead to tears and horrible bugs.
The function to do this, of course, is called step.Gam
. The documentation is only somewhat helpful — access it via help(step.Gam)
. It includes a rather minimalistic example to help you start off:
Gam.object <- gam(y~x+z, data=gam.data)
step.object <- step.Gam(Gam.object, scope=list("x"=~1+x+s(x,4)+s(x,6)+s(x,12),"z"=~1+z+s(z,4)))
Here, Gam.object
is the starting model, which will be used as the starting point in exploring our variable space and building other models. AIC is used to pick better models internally and is reported in the function output at every step.
One of the most important and confusing parts is the scope=list
component of the step.Gam
function call — this list is crucial as it defines which variables to explore while building up the function. You need to provide a named formula, named exactly like the predictor variable you want to stepwise. Thus in this function call, we have a list component called "x"
which corresponds to the formula ~1 + x + s(x, 4) + s(x, 6) + s(x, 12)
and another component called "z"
which corresponds to the formula ~1 + z + s(z, 4)
. These components correspond exactly to the predictors in the Gam.object
formula, which means this scope will work with this starting formula. step.Gam
will use these variables separated by plus signs as the building blocks of its models — so if you want to add extra variables for the stepwise selection to consider, simply add it after a plus sign — these can be, for example, natural splines, ns(x, knots=2)
.
Let’s do a slightly more advanced example to show you further how this works. The code can be found here. We are using a dataset called prostate
which is a few variables and about a hundred observations, trying to find a good model to predict the Cscore
response variable. We want to try fitting a smoothing spline model and select how many degrees of freedom we require for a good fit, judged by AIC. We will do backwards and forwards selection to illustrate — by default, step.Gam
scans in both directions.
First, let us define the scope list:
scope_list = list(
"lcavol" = ~1 + lcavol + s(lcavol, df=2) + s(lcavol, df=3) + s(lcavol, df =4) + s(lcavol, df=5),
"lweight" = ~1 + lweight + s(lweight, df=2) + s(lweight, df=3) + s(lweight, df=4) + s(lweight, df=5),
"age" = ~1 + age + s(age, df=2) + s(age, df=3) + s(age, df=4) + s(age, df=5),
"lbph" = ~1 + lbph + s(lbph, df=2) + s(lbph, df=3) + s(lbph, df=4) + s(lbph, df=5),
"lcp" = ~1 + lcp + s(lcp, df=2) + s(lcp, df=3) + s(lcp, df=4) + s(lcp, df=5),
"lpsa" = ~1 + lpsa + s(lpsa, df=2) + s(lpsa, df=3) + s(lpsa, df=4) + s(lpsa, df=5)
)
One can see that we’re only considering splines up to five degrees of freedom.
First, let’s pick a very random model to demonstrate how the function goes both ways by default.
start_model_rand = gam(Cscore~lpsa+s(lcavol, df=3), data=prostate)
And run the selection by typing:
step.Gam(start_model_rand, scope_list)
Which should produce the following output:
Start: Cscore ~ lpsa + s(lcavol, df = 3); AIC= 969
Step:1 Cscore ~ s(lcavol, df = 3) + s(lpsa, df = 2) ; AIC= 924
Step:2 Cscore ~ s(lcavol, df = 3) + s(lpsa, df = 3) ; AIC= 906
Step:3 Cscore ~ s(lcavol, df = 3) + lcp + s(lpsa, df = 3) ; AIC= 899
Step:4 Cscore ~ s(lcavol, df = 3) + s(lcp, df = 2) + s(lpsa, df = 3) ; AIC= 895
Step:5 Cscore ~ s(lcavol, df = 3) + s(lcp, df = 3) + s(lpsa, df = 3) ; AIC= 889
Step:6 Cscore ~ s(lcavol, df = 3) + s(lcp, df = 3) + s(lpsa, df = 4) ; AIC= 885
Step:7 Cscore ~ s(lcavol, df = 3) + age + s(lcp, df = 3) + s(lpsa, df = 4) ; AIC= 882
Step:8 Cscore ~ s(lcavol, df = 3) + age + s(lcp, df = 4) + s(lpsa, df = 4) ; AIC= 879
Step:9 Cscore ~ s(lcavol, df = 3) + age + s(lcp, df = 4) + s(lpsa, df = 5) ; AIC= 878
Step:10 Cscore ~ s(lcavol, df = 2) + age + s(lcp, df = 4) + s(lpsa, df = 5) ; AIC= 877
Step:11 Cscore ~ lcavol + age + s(lcp, df = 4) + s(lpsa, df = 5) ; AIC= 876
Call:
gam(formula = Cscore ~ lcavol + age + s(lcp, df = 4) + s(lpsa,
df = 5), data = prostate, trace = FALSE)Degrees of Freedom: 96 total; 85 Residual
Residual Deviance: 36186
Thus it returns the “best” model, best meaning it has the lowest AIC in the space it has searched — in this case being Cscore ~ lcavol + age + s(lcp, df = 4) + s(lpsa, df = 5)
.
Now let’s attempt forward stepwise selection. For this, we can use a somewhat minimalistic starting model that includes each variable (lpsa + lcavol
etc), using the dot formula operator to fill these in for us.
start_model_minimal = gam(Cscore~., data=prostate)
step.Gam(start_model_minimal, scope_list, direction = "forward")
Which produces:
Start: Cscore ~ .; AIC= 974
Step:1 Cscore ~ svi + lcavol + lweight + age + lbph + lcp + s(lpsa, df = 2) ; AIC= 917
Step:2 Cscore ~ svi + lcavol + lweight + age + lbph + lcp + s(lpsa, df = 3) ; AIC= 898
Step:3 Cscore ~ svi + lcavol + lweight + age + lbph + s(lcp, df = 2) + s(lpsa, df = 3) ; AIC= 894
Step:4 Cscore ~ svi + lcavol + lweight + age + lbph + s(lcp, df = 3) + s(lpsa, df = 3) ; AIC= 888
Step:5 Cscore ~ svi + lcavol + lweight + age + lbph + s(lcp, df = 3) + s(lpsa, df = 4) ; AIC= 884
Step:6 Cscore ~ svi + lcavol + lweight + age + lbph + s(lcp, df = 4) + s(lpsa, df = 4) ; AIC= 882
Step:7 Cscore ~ svi + lcavol + lweight + age + lbph + s(lcp, df = 4) + s(lpsa, df = 5) ; AIC= 880
Call:
gam(formula = Cscore ~ svi + lcavol + lweight + age + lbph +
s(lcp, df = 4) + s(lpsa, df = 5), data = prostate, trace = FALSE)Degrees of Freedom: 96 total; 82 Residual
Residual Deviance: 35649
Quite a different model this time — again a good reminder to be careful with interpreting the output of stepwise selection.
We can try backward selection too. In this case let’s start from a beefy case so the function will have something to cut off. We can also specify trace=2
to the function to see some more information about the selection process happening.
start_model = gam(Cscore~s(lcavol, df=5)+s(lweight, df=5)+s(age, df=5)+s(lbph, df=5)+s(lcp, df=5)+s(lpsa, df=5), data = prostate)step.Gam.mod(start_model, scope = scope_list, direction = "backward", trace=2)
Particularly, this verbose trace option will print every model the selection procedure has considered — very useful to see which places one can still explore manually. As the output is quite long, I provide only the last two steps here:
Step:18 Cscore ~ lcavol + age + s(lcp, df = 5) + s(lpsa, df = 5) ; AIC= 876
Trial: Cscore ~ 1 + 1 + age + 1 + s(lcp, df = 5) + s(lpsa, df = 5) ; AIC= 880
Trial: Cscore ~ lcavol + 1 + 1 + 1 + s(lcp, df = 5) + s(lpsa, df = 5) ; AIC= 880
Trial: Cscore ~ lcavol + 1 + age + 1 + s(lcp, df = 4) + s(lpsa, df = 5) ; AIC= 876
Trial: Cscore ~ lcavol + 1 + age + 1 + s(lcp, df = 5) + s(lpsa, df = 4) ; AIC= 878
Step:19 Cscore ~ lcavol + age + s(lcp, df = 4) + s(lpsa, df = 5) ; AIC= 876
Trial: Cscore ~ 1 + 1 + age + 1 + s(lcp, df = 4) + s(lpsa, df = 5) ; AIC= 879
Trial: Cscore ~ lcavol + 1 + 1 + 1 + s(lcp, df = 4) + s(lpsa, df = 5) ; AIC= 880
Trial: Cscore ~ lcavol + 1 + age + 1 + s(lcp, df = 3) + s(lpsa, df = 5) ; AIC= 878
Trial: Cscore ~ lcavol + 1 + age + 1 + s(lcp, df = 4) + s(lpsa, df = 4) ; AIC= 877
Call:
gam(formula = Cscore ~ lcavol + age + s(lcp, df = 4) + s(lpsa,
df = 5), data = prostate, trace = FALSE)Degrees of Freedom: 96 total; 85 Residual
Residual Deviance: 36186
I hope this step-by-step guide helped you understand how to work with this seemingly rarely used function. One important thing to note is that step.Gam
is written entirely in R — which means don’ t be afraid to put breakpoints into it and use your debugger to inspect variables at any given point in the code to understand what it wants and expects from you.