I want to evaluate the significance of including two smooth terms in a GAM by comparing that model to a baseline model that does not include those predictors. A simulation of this case can be reproduced using the code below:
library(mgcv)
dat <- gamSim(1,n=1000,dist="normal",scale=.1)
model <- gam(y ~ s(x0, bs = "cr", k = 20) + s(x1, bs = "cr", k = 20) + te(x2, x3, bs = "cr"), data = dat)
baseline <- gam(y ~ te(x2, x3, bs = "cr"), data = dat)
anova(baseline, model)
However, I would like to do this evaluation using a cross-validation sheme. Fitting the model on part of the data and evaluating the log-likelihood on a left out fold. For this I would calculate the log-likelihood on the left out fold manually. To calculate the p-value associated with this difference in log-likelihood I need the degrees of freedom. How do I determine the correct degrees of freedom for this comparison?
set.seed(123)
nfolds <- 5
dat$folds <- sample(1:nfolds, size = nrow(dat), replace = TRUE)
train_data <- dat[dat$folds != 1, ]
test_data <- dat[dat$folds == 1, ]
test_y <- test_data$y
model <- gam(y ~ s(x0, bs = "cr", k = 20) + s(x1, bs = "cr", k = 20) + te(x2, x3, bs = "cr"), data = train_data)
baseline <- gam(y ~ te(x2, x3, bs = "cr"), data = train_data)
pred_model <- predict(model, test_data)
pred_baseline <- predict(baseline, test_data)
stdev_model <- sigma(model)
stdev_baseline <- sigma(baseline)
loglik_model <- sum(dnorm(test_y, pred_model, stdev_model, log=TRUE))
loglik_baseline <- sum(dnorm(test_y, pred_baseline, stdev_baseline, log=TRUE))
delta <- loglik_model - loglik_baseline
teststat <- -2 * delta
pval <- pchisq(teststat, df = , lower.tail = FALSE)
