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.

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.