Is there a way to plot the tuning parameter selection in mgcv? Take a simple analysis using generalized additive models in mgcv: 
 x <- seq(0, 1, length=500)
 y <- sin(2*pi*x)^2 + rnorm(500, sd=0.2)
 m <- gam(y ~ s(x))

I know how to see the final tuning parameter selected (type m$sp). Is there a way to plot the GCV vs. smoothing parameter to graphically see how it was selected? I know you can do stuff like this in glmnet when you choose the elastic net penalty size. Anything like this for mgcv? Thanks
 A: The way to do this is to fit the GAM for known values of the smoothness parameter $\lambda$, and set the smoothness parameter for each fit. You can then extract the GCV score from each model and plot them against the sequence of $\lambda$ values.
To fix the smoothness parameter you use the sp argument. For your example data this means setting sp to be a length one vector. If you are doing this for multiple parameters you'll need to fit models over a grid of 2 or more smoothness parameters.
Here is an example using your dummy data set:
## packages
library("mgcv")
library("ggplot2")

## data
set.seed(42)
x <- seq(0, 1, length=500)
y <- sin(2*pi*x)^2 + rnorm(500, sd=0.2)
df <- data.frame(x = x, y = y)

## Generate GCV traces
lambda <- 1e-6
gcv <- numeric(length = 100)

## loop over values of lambda, fixing sp
for (i in seq_along(gcv)) {
    m <- gam(y ~ s(x, k = 20), data = df, method = "GCV.Cp", sp = lambda)
    gcv[i] <- m$gcv.ubre
    lambda <- lambda * 1.1
}

## gather results
res <- data.frame(lambda = 1e-6 * 1.1^{0:99}, gcv = gcv)

## plot
ggplot(res, aes(x = lambda, y = gcv)) +
    geom_line() +
    scale_x_log10() +
    xlab(expression(lambda)) +
    ylab("Score")

This produces

You can do this for the other fitting methods available in mgcv also; the score is always in $gcv.ubre, despite the name of that component.
You'll notice that it take a lot longer to fit those models than it does to find the optimum value of $\lambda$ using GCV if you fitted without setting sp. mgcv has very efficient algorithms for selecting the optimal lambda and the steps it takes along the GCV surface (curev) will not produce a nice curve and typically aren't worth saving in the model object as a result. Hence the need to do the fitting iteratively, as shown above, to yield a nice curve of values over $\lambda$.
The above methods requires a little tuning to get right; you need to set lambda small enough initially and not increase it took quickly to home in on the minimum of the curve.
