2
$\begingroup$

Consider the following very simple example:

set.seed( 1 )
SimData <- data.frame( X = runif( 1000, 0, 1 ) )
SimData$Y <- rnbinom( nrow( SimData ), mu = 100*sin( SimData$X*2*pi )+100, size = 10 )
plot( gam( Y ~ s( X ), data = SimData, family = nb( link = log ) ) )

Here is the result:

enter image description here

The effective degrees of freedom is very-very close to 9, giving rise to the suspicion that the default basis dimension is not large enough. Let's increase it:

plot( gam( Y ~ s( X, k = 20 ), data = SimData, family = nb( link = log ) ) )

enter image description here

Hm. The overall picture is essentially the same, yet, the EDF suggests that the basis dimension is still not large enough!

Even if we increase k to 30, the EDF still gets larger (24.8), so, essentially EDF seems to simply follow the k limit, which is pretty bizarre... (especially for such a simple functional form, and especially that is was already well captured by the default model).

EDIT (07 Sep, 2017): According to an answer to the original question, the application of adaptive splines (bs="ad") can be the solution to this problem. However...

Let's take another simple example:

set.seed( 1 )
SimData <- data.frame( X = runif( 1000, 0, 1 ) )
SimData$Y <- rnbinom( nrow( SimData ), mu = 100*sin( SimData$X*2*pi*3 )+1000, size = 10 )
plot( gam( Y ~ s( X ), data = SimData, family = nb( link = log ), method = "REML" ) )

enter image description here

Seems perfect! Let's now "spoil" it:

SimData$Y[ SimData$X<=0.05|SimData$X>=0.95 ] <- 0

This gives rise to the original problem: EDF is 8.92 with default k, 18.91 if k=20, 46.74 if k=50 etc. As an illustration, here is the k=50 case:

enter image description here

So indeed, we have the original problem. Let's try therefore bs="ad":

plot( gam( Y ~ s( X, bs = "ad" ), data = SimData, family = nb( link = log ), method = "REML" ) )

enter image description here

So unfortunately the problem remained, even with bs="ad"! (It's really the same situation: increasing k to 50 gives an EDF of 44.72, k=100 gives 78.03. That's why I decided to edit this question instead of starting a new one: this seems to be somehow the same story...)

$\endgroup$
  • 1
    $\begingroup$ As a further comment to my answer, it occurs to me that on the log scale, your data is reminiscent of the head acceleration example data (available as data set mcycle in MASS IIRC) that is often used to motivate adaptive splines. Not quite the same (there the response is continuous, not integer), but the true model on the log scale in your example is similar in shape to the true model on the raw data scale in the mcycle data set. $\endgroup$ – Reinstate Monica - G. Simpson Sep 7 '17 at 16:36
2
$\begingroup$

Fitting on the log scale is causing the problem; on that scale, as shown in your plot, the degree of smoothness is not really the same across the function; at low values of x there is little curvature in the true function, but there is a lot of curvature around the minimum at ~ x = 7.5. As, by default, the spline wiggliness is controlled by a single smoothness parameter, this has the effect of assuming the same degree of smoothness or curvature everywhere. As you increase k, the model is able to adapt to the data around the rapid rise and fall, but in doing so it also adapts, because it has too given the smoothness parameter/penalty, to sampling noise elsewhere in the data.

With these data I was able to achieve an acceptable fit by changing the link function to the sqrt link via family = nb(link = "sqrt"). On this scale the rise and fall is not as sharp or constrained to a narrow range of x as it is on the log scale.

m <-  gam(Y ~ s(X, k = 20), data = SimData, family = nb(link = "sqrt"), method = "REML")

The fitted spline is:

enter image description here

There is still some evidence of adaptation to sampling noise in the region where x is low, but this is quite trivial and by any reasonable metric would be indistinguishable from the true function I would guess.

The basis dimension test implemented in gam.check() also seems reasonably happy with this — note that I increased k to 20 but this appears sufficient for this model, whereas k = 10 wasn't:

> gam.check(m)

Method: REML   Optimizer: outer newton
full convergence after 4 iterations.
Gradient range [2.520842e-06,0.0003287577]
(score 4439.901 & scale 1).
Hessian positive definite, eigenvalue range [5.53589,345.0859].
Model rank =  20 / 20 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

       k'  edf k-index p-value
s(X) 19.0 14.9    1.01    0.76

The diagnostic plots are also pretty reasonable:

enter image description here

If you want to fit on the log scale, an adaptive spline can be used, which also appears to give good results for these data. An adaptive spline uses multiple smoothness penalties (parameters) for each spline and in effect is like fitting several separate splines to parts of data, each with their own smoothness parameter. Such splines require a lot more data to achieve good results however.

mad <-  gam(Y ~ s(X, k = 20, bs = "ad"), data = SimData, family = nb, method = "REML")

The fitted spline is now much smoother for low values of x, but there is perhaps a little overfitting to sampling noise at large x. On the scale of the response this is going to be trivial however.

enter image description here

and gam.check() is also happy that k = 20 is sufficient for this basis/model/data:

> gam.check(mad)

Method: REML   Optimizer: outer newton
full convergence after 11 iterations.
Gradient range [-0.0002868044,0.001512152]
(score 4461.377 & scale 1).
Hessian positive definite, eigenvalue range [0.0005701379,344.5224].
Model rank =  20 / 20 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

       k'  edf k-index p-value
s(X) 19.0 12.9    0.97    0.24

and again, the diagnostics don't look that bad either

enter image description here

although there is a hint of non-constant variance in residuals in the plot to the top right.

I'll also add that the test in gam.check() for the dimension of the basis is not infallible. It is a guide at best. If the model looks good in other respects, I would place less emphasis on the dimension test results. However, for the default k = 10 model you fitted, the model isn't really sufficiently wiggly to capture the low counts around the dip and it is quite badly biased there.

enter image description here

I'd argue that whilst just increasing k wasn't helping here, the initial model was also not adequate for these data.

$\endgroup$
  • $\begingroup$ Fantastic Gavin, thank you very much!! I don't really want to change the link function, as my original problem is medical, and at that field they really like coefficients that have a well-known interpretation (through exponentiation, in this case)... Adaptive splines, however, seem extremely promising! My only problem is that while it solves this problem, I could come up with another one, where the problem (the same one) remains - even with adaptive splines! I edited the question, I'd be really interested to hear your opinion. Thank you again! $\endgroup$ – Tamas Ferenci Sep 7 '17 at 20:14
  • $\begingroup$ @TamasFerenci I would argue that a GAM (or at least a simple one such as that which you fitted to the example data) is a poor choice for the data generating function in the case of your second example; the mean shift at the margins of the covariate are hard thresholds and as such are incompatible with the assumption of the GAM that truth is smooth or at least can be adequately approximated by a smooth function. $\endgroup$ – Reinstate Monica - G. Simpson Sep 7 '17 at 22:52
  • 1
    $\begingroup$ OK, got it! You're right. I accepted your answer, thanks again! $\endgroup$ – Tamas Ferenci Sep 8 '17 at 4:49

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.