gam object produced by gam() and produces various useful
summaries from it. (See sink to divert output to a file.)## S3 method for class 'gam':
summary(object, dispersion=NULL, freq=FALSE, p.type = 0, ...)## S3 method for class 'summary.gam':
print(x,digits = max(3, getOption("digits") - 3),
signif.stars = getOption("show.signif.stars"),...)
gam object as produced by gam().summary.gam object produced by summary.gam().NULL to use estimate or
default (e.g. 1 for Poisson).TRUE then
the frequentist covariance matrix of the parameters is used instead.summary.gam produces a list of summary information for a fitted gam object.p.coeff's divided by their standard errors.r.sq does not include any offset in the one parameter model.dev.expl can be substantially lower than r.sq when an offset is present.freq=TRUE), divided
by scale parameter.freq=TRUE).P-values for terms penalized via `paraPen' are unlikely to be correct.
print.summary.gam tries to print various bits of summary information useful for term selection in a pretty way.
Unless p.type=5, p-values for smooth terms are usually based on a
test statistic motivated by an extension of Nychka's (1988) analysis of the frequentist properties
of Bayesian confidence intervals for smooths.
These have better frequentist performance (in terms of power and distribution under the null)
than the alternative strictly frequentist approximation. When the Bayesian intervals have good
across the function properties then the p-values have close to the correct null distribution
and reasonable power (but there are no optimality results for the power). Full details are in Wood (2013),
although what is computed is actually a slight variant in which the components of the test statistic are weighted by the iterative fitting weights.
Let $\bf f$ denote the vector of values of a smooth term evaluated at the original covariate values and let ${\bf V}_f$ denote the corresponding Bayesian covariance matrix. Let ${\bf V}_f^{r-}$ denote the rank $r$ pseudoinverse of ${\bf V}_f$, where $r$ is the EDF for the term. The statistic used is then $$T = {\bf f}^T {\bf V}_f^{r-}{\bf f}$$ (this can be calculated efficiently without forming the pseudoinverse explicitly). $T$ is compared to an approximation to an appropriate weighted sum of chi-squared random variables.
The non-integer rank truncated inverse is constructed to give an
approximation varying smoothly between the bounding integer rank approximations, while yielding test statistics with the correct mean and variance under the null. Alternatively (p.type==1) $r$ is obtained by
biased rounding of the EDF: values less than .05 above the preceding integer are rounded down, while other values are rounded up. Another option (p.type==-1) uses a statistic of formal rank given by the number of coefficients for the smooth, but with its terms weighted by the eigenvalues of the covariance matrix, so that penalized terms are down-weighted, but the null distribution requires simulation. Other options for p.type are 2 (naive rounding), 3 (round up), 4 (numerical rank determination): these are poor options for theoretically known reasons, and will generate a warning.
Note that for terms with no unpenalized terms the Nychka (1988) requirement for smoothing bias to be substantially less than variance breaks down (see e.g. appendix of Marra and Wood, 2012), and this results in incorrect null distribution for p-values computed using the above approach. In this case it is necessary to use an alternative approach designed for random effects variance components, and this is done.
In this zero-dimensional null space/random effects case, the p-values are again conditional on the smoothing parameters/variance component estimates, and may therefore be somewhat too low when these are subject to large uncertainty. The idea is to condition on the smoothing parameter estimates, and then to use the likelihood ratio test statistic conditional on those estimates. The distribution of this test statistic under the null is computable as a weighted sum of chi-squared random variables.
In simulations the p-values have best behaviour under ML smoothness selection, with REML coming second. In general the p-values behave well, but conditioning on the smoothing parameters means that they may be somewhat too low when smoothing parameters are highly uncertain. High uncertainty happens in particular when smoothing parameters are poorly identified, which can occur with nested smooths or highly correlated covariates (high concurvity).
If p.type=5 then the frequentist approximation for p-values of smooth terms described in section
4.8.5 of Wood (2006) is used. The approximation is not great. If ${\bf p}_i$
is the parameter vector for the ith smooth term, and this term has estimated
covariance matrix ${\bf V}_i$ then the
statistic is ${\bf p}_i^\prime {\bf V}_i^{k-} {\bf
p}_i$, where ${\bf V}^{k-}_i$ is the rank k
pseudo-inverse of ${\bf V_i}$, and k is estimated rank of
${\bf V_i}$. p-values are obtained as follows. In the case of
known dispersion parameter, they are obtained by comparing the chi.sq statistic to the
chi-squared distribution with k degrees of freedom, where k is the estimated
rank of ${\bf V_i}$. If the dispersion parameter is unknown (in
which case it will have been estimated) the statistic is compared
to an F distribution with k upper d.f. and lower d.f. given by the residual degrees of freedom for the model.
Typically the p-values will be somewhat too low.
By default the p-values for parametric model terms are also based on Wald tests using the Bayesian
covariance matrix for the coefficients. This is appropriate when there are "re" terms present, and is
otherwise rather similar to the results using the frequentist covariance matrix (freq=TRUE), since
the parametric terms themselves are usually unpenalized. Default P-values for parameteric terms that are
penalized using the paraPen argument will not be good. However if such terms represent conventional
random effects with full rank penalties, then setting freq=TRUE is appropriate.
Marra, G and S.N. Wood (2012) Coverage Properties of Confidence Intervals for Generalized Additive Model Components. Scandinavian Journal of Statistics, 39(1), 53-74.
Nychka (1988) Bayesian Confidence Intervals for Smoothing Splines. Journal of the American Statistical Association 83:1134-1143.
Wood S.N. (2006) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press.
gam, predict.gam,
gam.check, anova.gam, gam.vcomp, sp.vcovlibrary(mgcv)
set.seed(0)
dat <- gamSim(1,n=200,scale=2) ## simulate data
b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat)
plot(b,pages=1)
summary(b)
## now check the p-values by using a pure regression spline.....
b.d <- round(summary(b)$edf)+1 ## get edf per smooth
b.d <- pmax(b.d,3) # can't have basis dimension less than 3!
bc<-gam(y~s(x0,k=b.d[1],fx=TRUE)+s(x1,k=b.d[2],fx=TRUE)+
s(x2,k=b.d[3],fx=TRUE)+s(x3,k=b.d[4],fx=TRUE),data=dat)
plot(bc,pages=1)
summary(bc)
## Example where some p-values are less reliable...
dat <- gamSim(6,n=200,scale=2)
b <- gam(y~s(x0,m=1)+s(x1)+s(x2)+s(x3)+s(fac,bs="re"),data=dat)
## Here s(x0,m=1) can be penalized to zero, so p-value approximation
## cruder than usual...
summary(b)
## Now force summary to report approximate p-value for "re"
## terms as well. In this case the p-value is OK, since the
## random effect structure is so simple.
summary(b,all.p=TRUE)
## p-value check - increase k to make this useful!
k<-20;n <- 200;p <- rep(NA,k)
for (i in 1:k)
{ b<-gam(y~te(x,z),data=data.frame(y=rnorm(n),x=runif(n),z=runif(n)),
method="ML")
p[i]<-summary(b)$s.p[1]
}
plot(((1:k)-0.5)/k,sort(p))
abline(0,1,col=2)
ks.test(p,"punif") ## how close to uniform are the p-values?Run the code above in your browser using DataLab