Alexander März
2013-Apr-20 15:00 UTC
[R] Calculate confidence intervals in mgcv for unconditional on the, smoothing parameters
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]]