0
$\begingroup$

I am attempting to run a relative risk model using the MCGV package in R with multiple predictors. The syntax for my model in R is the following:

model1=gam(datacount~s(trend,k=6*2)+s(T,k=5),family=quasipoisson)

where datacount is my variable of interest, trend is continuous time and T is the temperature. While I have five knots defined for the temperature, a plot of my results show only one knot within the spline in the y-axis label.

enter image description here

How can I force things to make sure that I actually am ending up with five knots in my temperature spline? Any info would be greatly appreciated.

$\endgroup$
1
  • $\begingroup$ I don't know why someone downvoted this? I suppose it is somewhat OT as it pertains to a specific implementation of a GAM and not GAMs in general. But if you think that it should have a vote to close. $\endgroup$ Commented Aug 31, 2017 at 15:37

1 Answer 1

2
$\begingroup$

The thin plate regression spline (TPRS) basis you are using (that's the default in s()) has a "knot" for each unique data point. This is not very efficient, so what Simon Wood does in the mgcv package is take the full TPRS, so all the $N$ basis functions and performs an eigen decomposition of the full space spanned by the basis. Then the first k eigenvectors of the decomposed basis are retained for use as a new set of basis functions that maintain most of the signal of the full TPRS basis but at much reduced cost.

Because of identifiability constraints on the splines k may not be exactly the k you specify; for the TPRS the final basis that is actually used in the model will be of dimension k-1 given all the defaults. Take that into account when fixing the basis dimension, k, at some known size.

To achieve what you want, you can instruct gam() to not do any form of smoothness selection. In that case, the model is really just a GLM using a basis expansion of the covariate. To fix the degree of smoothness at whatever dimension you specify, add fx = TRUE to the s(...) term, hence you model would be

gam(datacount ~ s(trend, k=6*2) + s(T, k=6, fx=TRUE), family=quasipoisson)

Remember, this is not using 5 knots but rather 5 basis functions. The plot is not showing that there is only a single knot. What this plot shows to me is that the model used 1 effective degree of freedom to represent the effect of T.

What you may confuse as a knot is the point where the confidence band tends to zero as the estimated effect tends also to zero. This is a result of the identifiability constraint applied to the basis expansion of T — the TPRS basis expansion of T contains two (by default, for this model) basis functions that are perfectly smooth

  1. a constant, flat line, basis function, and
  2. a linear basis function

These are in the null space of the penalty matrix for this basis. The first of these is not identifiable when your model has an intercept term; I could add a constant to the coefficient for the intercept and subtract that same constant from the coefficient for that basis function and I would not change the fit of the model at all. Hence there are infinite solutions to the model. What gam() does is apply a constraint that the spline sums to zero (that is why the spline is centred about 0 in the plot) over the range of T. This causes problems for the coverage properties of the confidence interval in such cases, which can be addressed by adding seWithMean = TRUE to the plot() call that generated the plot you show. However, this is not the result of using only 1 knot.

(Note that these identifiability constraints do nothing with the second basis function in the penalty null space, the linear function, and it remains in the basis used to fit the model; this is why the GAM can't smooth back to a null or constant effect (i.e. remove a term entirely from the model), because the smoothness penalty does not apply to functions in the penalty null space.)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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