Monte Carlo test for comparing curvature of binomial response surfaces from effective degrees of freedom of GCV-fitted splines

Question

I would like to compare the curvature of two response surfaces, each of the form:

binomial ~ continuous variables1-5

I think it would be appropriate to use the effective degrees of freedom of a thin-plate spline fit to the data with generalized cross-validation as an estimate of the curvature. Is this correct? One of the surfaces is complex (effective df = 40), so parametric approaches are not practical here.

To test for a significant difference in curvature, I calculated the observed difference in effective df for thin-plate splines fit to the data using generalized cross-validation (using the Tps function in the Fields package in R). I then calculated this difference for 1000 samples with the binomial response variables randomly reassigned for each dataset to generate an empirical null distribution and then counted results >= the observed difference to get a one-tailed P-value.

My problem is this: the null distribution of effective degrees of freedom is incredibly large, spanning -10 to 120 df for the randomized datasets, even though the observed difference in df was only 20 df (40 and 20 df, respectively)!

Is this something I should be worried about? Why is GCV fitting such complex curves to randomized response surfaces? This problem seems to be inherent to resampling smoothing splines, GAMs, thin-plate splines, REML methods, and all variations of bootstrapping my various datasets: random binomial response surfaces often result in dramatically over-fitted splines.

Is there a better nonparametric way to compare curvature of two surfaces?

Thanks!

Answer 1

I'm working on a similar problem but I'm also not sure how to proceed, so maybe we can help each other. Along the way I've found Simon Wood's book "Generalized Additive Models: An Introduction with R" very helpful - have you looked at that?

I'm using a GAM with cubic smoothing splines to model two time series, using the mgcv package in R (letting GCV identify the "best" smoothing parameter). The data are population measurements of seabirds on two islands and I wanted to compare the two populations in terms of stability in time. In other words I want to compare wiggliness, or curvature of the smooths across the time range of my data.

My thought was to do a bootstrapping test like yours, except I was going to compare the integral of the squared second derivative of the spline function, say J(f), where f is the fitted spline. J(f) approximates the total curvature of the spline, so I thought it was a good metric for wiggliness. It also appears in the least squares objective function which is minimized when fitting the splines:

SSE + lambda*J(f)

where lambda is the smoothing parameter. J(f) has a nice representation as a quadratic in the vector of coefficients for the fitted spline (at least using the formulation of basis functions used in Wood, 2006). Have you considered using this, or some other approximation of total curvature, rather than EDF?

I imagine that when you resample from both datasets, you sometimes obtain points that are very close together in the space of your covariates, but very far apart in terms of your response variable. This would force your smooths to curve sharply in some spots to minimize the squared error, at the expense of an increase in EDF.

This is my worry with my own dataset. I haven't tried the bootstrap yet but I will let you know how it goes if you're interested.