Extracting variance component parameters from fitted mgcv::gam object

In their seminal book "Semiparametric Regression", Ruppert, Wand and Carroll advocate for Monte Carlo simulation of the critical values of likelihood ratio statistic when using the likelihood ratio procedure for hypothesis testing in the framework of mixed models (see 4.8.2 pages 105-7 for details). As the authors put it, "The idea is to set the values of all fixed effects and variance component parameters [$\epsilon$-values and random effects] equal to their estimates under the null distribution and then to simulate the distribution of the likelihood ratio test under the null model at the parameters and with the covariates equal to their observed values."

Suppose we have a time course data comprised of 39 observations collected at 5 time points (0, 3, 7, 14, and 21 days) from 39 patients that are separated in two groups (A and B). Further assume we want to fit cubic smoothing splines to this data using mgcv::gam and perform a likelihood ratio test to determine if temporal trends are different between the two groups:

# load data
data <- c(4.61,4.61,4.36,4.42,4.47,4.39,4.52,4.39,4.62,5.06,4.72,5.09,4.48,4.81,4.57,4.65,4.85,4.82,4.71,4.77,4.66,4.33,4.51,4.47,4.35,4.53,4.51,4.54,4.25,4.28,4.35,4.34,4.52,4.29,4.37,4.42,4.57,4.44,4.33)
time <- c(rep(0,8),rep(3,8),rep(7,8),rep(14,7),rep(21,8))
groupB <- c(1,1,1,1,0,0,0,0,1,1,1,1,0,0,0,0,1,1,1,1,0,0,0,0,1,1,1,1,0,0,0,1,1,1,1,0,0,0,0)

# fit models
require(mgcv)
alternative.fit <- gam(data ~ groupB + s(time, k=5, bs="cr") + s(time, k=5, bs="cr", by=groupB), family=gaussian(), method="ML")
null.fit <- gam(data ~ groupB + s(time, k=5, bs="cr"), family=gaussian(), method="ML")

# extract log-likelihoods and calculate LR-statistic
LR.stat <- (-2) * (logLik(null.fit) - logLik(alternative.fit))


Figure below shows smooth curves for each group obtained from the alternative fit (code not shown for brevity):

One can assume that the LR statistic follows a chi-squared distribution under the null and compare the two models using 'anova' function:

anova(null.fit,alternative.fit, test = "Chisq")

Analysis of Deviance Table

Model 1: data ~ groupB + s(time, k = 5, bs = "cr")
Model 2: data ~ groupB + s(time, k = 5, bs = "cr") + s(time, k = 5, bs = "cr", by = groupB)
Resid. Df Resid. Dev    Df Deviance Pr(>Chi)
1    32.815    0.62492
2    29.494    0.45838 3.321  0.16654  0.01399 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


This gives us an approximate p-value of 0.01, suggesting a (marginally) significant interaction effect in agreement with our visual assessment. Now suppose we want to verify this result using Monte Carlo simulation as outlined above. The estimate of fixed effects (i.e. intercept and groupB) is readily available via 'summary':

summary(null.fit)

Family: gaussian

Formula:
data ~ groupB + s(time, k = 5, bs = "cr")

Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  4.46419    0.03140 142.156  < 2e-16 ***
groupB       0.14384    0.04388   3.278  0.00245 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
edf Ref.df     F  p-value
s(time) 3.599  3.892 10.44 3.26e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.543   Deviance explained = 59.8%
-ML = -19.717  Scale est. = 0.01871   n = 39


My question is how to extract variance component parameters associated with $\epsilon$-values and random effects, and how to use their estimates to simulate data under the null.

My initial thought was to extract residuals (i.e. $\epsilon$-values) and calculate their variance (let's call it var.epsilon), which will then allow us to generate random noise from N(0,var.epsilon). As for the other variance component parameter, I thought about using 'gam.vcomp' to obtain the variance of random effect associated with smoothing (let's call it var.smooth), and similarly generate random effects from N(0,var.smooth). Finally, to simulate data under the null we simply add all terms and repeat this process a large number of times:

# variance of epsilon values (unmodeled noise)
var.epsilon <- var(null.fit\$residuals)
# variance of smoothing random effect introduced via method="ML"
var.smooth <- (gam.vcomp(null.fit)[1,1])^2
# adding all components to simulate data under the null
simData <- (4.46419) + (0.14384)*(groupB) + rnorm(39,0,sqrt(var.smooth)) + rnorm(39,0,sqrt(var.epsilon))


Is this procedure correct? If we had additional random effects in the null model, for example random patient effect introduced by longitudinal sampling, would the above procedure still be valid as long as we include all variance components in our simulation?