Dear R-Help members,

I am using Simon Wood`s mgcv package version1.7-22and R version 3.0.0 
(2013-04-03) for fitting a GAM-Model to the LIDAR Data contained in the 
"SemiPar" package. Here is the code for fitting the model and for 
plotting the result:

data("lidar")
attach(lidar)

###
# mgcv fitting
###
gam_fit <- gam(logratio ~ s(range, k = 40, bs = "cr"), gamma = 1.4, 
method = "REML")
plot(gam_fit, res = TRUE, shade = TRUE, pch = "°")

Since the confidence interval in the above plot is pointwise, I would 
like to calculate the critical value needed in order to construct a 
simultaneous confidence band. Here isthe code(it`s a bit lengthy and has 
to be tidied up but I wanted to to some manual calculations on the 
objects afterwards)

###
# Simulate 100000 coefficient vectors from the posterior for ß and 
calculate the critical value
###
cov_beta_diff <- vcov(gam_fit)
basis_matrix <- model.matrix(gam_fit)
reml_sigma_squared <- gam.vcomp(gam_fit)
sigma_2_epsilon <- reml_sigma_squared[2,1]^2
smooth_matrix <- 
basis_matrix%*%vcov(gam_fit,dispersion=1)%*%t(basis_matrix)
stdv_f_x <- sqrt(diag(sigma_2_epsilon * smooth_matrix))

set.seed(123)
R_sim <- 100000
mean_val <- rep(0, nrow(cov_beta_diff))
coef_diff <- mvrnorm(n = R_sim, mu = mean_val , Sigma = cov_beta_diff)
Xp <- predict(gam_fit, type = "lpmatrix")
sim_bet <- Xp %*% t(coef_diff)
abs_ratio <- abs(sim_bet / stdv_f_x)

ratio_max <- apply(abs_ratio, 2, max)
cv <- quantile(ratio_max, probs = 0.95)

So everything is fine until this point. My problem is:everything has 
been presented conditional on the smoothing parameter estimated in 
gam_fit. How can I "calculate" a critical value for thesimultaneous 
confidence band that is not conditioned on the estimated smoothing 
parameter? Simon Wood seems to offer a solution in his book "Generalized 
Additive Models: An introduction with R". But I couldn`t figure where to 
look for it and out how to implement it using R.

Any help is appreciated

Best
Alex


	[[alternative HTML version deleted]]

