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]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
