For a listing of all routines in the mgcv
package type:
library(help="mgcv")
mgcv(y,X,sp,S,off,C=NULL,w=rep(1,length(y)),H=NULL,
scale=1,gcv=TRUE,control=mgcv.control())
ncol(X)
must give the number of model parameters, while nrow(X)
should give the number of data.control$fixed==TRUE
then these are taken as being the
smoothing parameters. Otherwise any positive values are assumed to be initial estimates and negative values to
signal auto-initialization.off[i]
indicates the element of the parameter vector that
S[[i]][1,1]
relates to.p
is the model parameter vector and k=nrow(S[[i]])-1
, then the ith penalty is given by
t(p[off[i]:(off[i
y
).S
. mgcv
cannot mix fixed and estimated penalties.gcv
is TRUE then smoothing parameters are estimated by GCV,
otherwise UBRE is used.mgcv.control
.g
above - i.e. the leading diagonal of the Hessian.}
TRUE
if the second smoothing parameter guess improved the GCV/UBRE score. (Please report examples
where this is FALSE
)}
TRUE
if the algorithm terminated by failing to improve the GCV/UBRE score rather than by "converging".
Not necessarily a problem, but check the above derivative information quite carefully.}The method operates by alternating very efficient direct searches for
$\rho$
with Newton or steepest descent updates of the logs of the $\lambda_i$.
Because the GCV/UBRE scores are flat w.r.t. very large or very small $\lambda_i$,
it's important to get good starting parameters, and to be careful not to step into a flat region
of the smoothing parameter space. For this reason the algorithm rescales any Newton step that
would result in a $log(\lambda_i)$ change of more than 5. Newton steps are only used
if the Hessian of the GCV/UBRE is postive definite, otherwise steepest descent is used. Similarly steepest
descent is used if the Newton step has to be contracted too far (indicating that the quadratic model
underlying Newton is poor). All initial steepest descent steps are scaled so that their largest component is
1. However a step is calculated, it is never expanded if it is successful (to avoid flat portions of the objective),
but steps are successively halved if they do not decrease the GCV/UBRE score, until they do, or the direction is deemed to have
failed. M$conv
provides some convergence diagnostics.
The method is coded in C
and is intended to be portable. It should be
noted that seriously ill conditioned problems (i.e. with close to column rank
deficiency in the design matrix) may cause problems, especially if weights vary
wildly between observations.
Wood, S.N. (2000) Modelling and Smoothing Parameter Estimation with Multiple Quadratic Penalties. J.R.Statist.Soc.B 62(2):413-428
gam
,
magic
library(help="mgcv") # listing of all routines
set.seed(1);n<-400;sig2<-4
x0 <- runif(n, 0, 1);x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1);x3 <- runif(n, 0, 1)
f <- 2 * sin(pi * x0)
f <- f + exp(2 * x1) - 3.75887
f <- f+0.2*x2^11*(10*(1-x2))^6+10*(10*x2)^3*(1-x2)^10-1.396
e <- rnorm(n, 0, sqrt(sig2))
y <- f + e
# set up additive model
G<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),fit=FALSE)
# fit using mgcv
mgfit<-mgcv(G$y,G$X,G$sp,G$S,G$off,C=G$C)
Run the code above in your browser using DataLab