| gam.models {mgcv} | R Documentation |
This page is intended to provide some more information on how to specify GAMs. A GAM is a GLM in which the linear predictor depends, in part, on a sum of smooth functions of predictors and (possibly) linear functionals of smooth functions of (possibly dummy) predictors.
Specifically let y_i denote an independent random variable with mean mu_i and an exponential family distribution, or failing that a known mean variance relationship suitable for use of quasi-likelihood methods. Then the the linear predictor of a GAM has a structure something like
g(mu_i)=X_i b + f_1(x_1i,x_2i) + f_2(x_3i) + L_i f_3(x_4) + ...
where g is a known smooth monotonic `link' function, X_i b is the parametric part of the linear predictor, the x_j are predictor variables, the f_j are smooth functions and L_i is some linear functional of f_3. There may of course be multiple linear functional terms, or none.
The key idea here is that the dependence of the response on the predictors can be represented as a parametric sub-model plus the sum of some (functionals of) smooth functions of one or more of the predictor variables. Thus the model is quite flexible relative to strictly parametric linear or generalized linear models, but still has much more structure than the completely general model that says that the response is just some smooth function of all the covariates.
Note one important point. In order for the model to be identifiable
the smooth functions usually have to be constrained to have zero mean (usually
taken over the set of covariate values). The constraint is needed if the term involving the
smooth includes a constant function in its span. gam always applies such constraints
unless there is a by variable present, in which case an assessment is made of whether
the constraint is needed or not (see below).
The following sections discuss specifying model structures for gam.
Specification of the distribution and link function is done using the family
argument to gam and works in the same way as for glm.
This page therefore concentrates on the model formula for gam.
Consider the example model.
g(mu_i) = b_0 + b_1 x_1i + b_2 x_2i + f1(x_3i) + f2(x_4i,x_5i)
where the response variables y_i has expectation mu_i and g is a link function.
The gam formula for this would be
y ~ x1 + x2 + s(x3) + s(x4,x5).
This would use the default basis for the smooths (a thin plate
regression spline basis for each), with automatic selection of the
effective degrees of freedom for both smooths. The dimension of the
smoothing basis is given a default value as well (the dimension of the
basis sets an upper limit on the maximum possible degrees of
freedom for the basis - the limit is typically one less than basis
dimension). Full details of how to control smooths are given in
s and te, and further discussion of basis
dimension choice can be found in choose.k.
For the moment suppose that we would like to change
the basis of the first smooth to a cubic regression spline basis with
a dimension of 20, while fixing the second term at 25 degrees of
freedom. The appropriate formula would be:
y ~ x1 + x2 + s(x3,bs="cr",k=20) + s(x4,x5,k=26,fx=TRUE).
The above assumes that x_4 and x_5 are naturally on
similar scales (e.g. they might be co-ordinates), so that isotropic smoothing
is appropriate. If this assumption is false then tensor product smoothing might be
better (see te).
y ~ x1 + x2 + s(x3) + te(x4,x5)
would generate a tensor product smooth of x_4 and x_5.
By default this smooth would have basis dimension 25 and use cubic regression spline marginals.
Varying the defaults is easy. For example
y ~ x1 + x2 + s(x3) + te(x4,x5,bs=c("cr","ps"),k=c(6,7))
specifies that the tensor product should use a rank 6 cubic regression spline marginal
and a rank 7 P-spline marginal to create a smooth with basis dimension 42.
Another common modelling task is to decide if a model like:
E(y)=f(x,z)
is really necessary or whether:
E(y)=f1(x)+f2(z)
would do just as well. One possibility is to examine the results
of fitting:
y ~ s(x) + s(z) + s(x,z).
gam automatically generates side conditions to make this model
identifiable. You can also estimate `overlapping' models like:
y ~ s(x,z) + s(z,v).
by variables are the means for constructing `varying-coefficient models' (geographic regression models) and
for letting smooths `interact' with factors or parametric terms. They are also the key to specifying general linear
functionals of smooths.
The s and te terms used to specify smooths accept an argument by,
which is a numeric or factor variable of the same dimension as the covariates of the smooth.
If a by variable is numeric, then its ith element multiples the ith
row of the model matrix corresponding to the smooth term concerned. If a by variable is a factor then it generates
an indicator vector for each level of the factor. The model matrix for the smooth term is then replicated for each factor level,
and each copy has its rows multiplied by the corresponding rows of its
indicator variable. The smoothness penalties are also duplicated for each
factor level. In short a different smooth is generated
for each factor level (the id argument to s and te can be used to force all
such smooths to have the same smoothing parameter).
As an example, consider the model
E(y_i) = b_0 + f(x_i)z_i
where f is a smooth function, and z_i is a numeric variable.
The appropriate formula is:
y ~ s(x,by=z)
- the by argument ensures that the smooth function gets multiplied by
covariate z. Note that when using factor variables, then the presence of centering constraints usually
means that the by variable should be included as a parametric term, as well.
The example code below also illustrates the use of factor by variables.
by variables may be supplied as numeric matrices as part of specifying general linear functional terms.
If a by variable is present and numeric (rather than a factor) then the corresponding smooth is only subjected
to an identifiability constraint if (i) the by variable is a constant vector, or, (ii) for a matrix
by variable, L, if L%*%rep(1,ncol(L)) is constant or (iii) if a user defined smooth constructor
supplies an identifiability constraint explicitly.
It is sometimes desirable to insist that different smooth terms have the same degree of smoothness.
This can be done by using the id argument to s or te terms. Smooths
which share an id will have the same smoothing parameter. Really this only makes sense if the
smooths use the same basis functions, and the default behaviour is to force this to happen: all smooths
sharing an id have the same basis functions as the first smooth occurring with that id. Note
that if you want exactly the same function for each smooth, then this is best achieved by making use of the
summation convention covered under `linear functional terms'.
As an example suppose that E(y_i)=mu_i and
g(mu_i) = f1(x_1i) + f2(x_2i,x_3i) + f3(x_4i)
but that f1 and f3 should have the same smoothing parameters (and x_2
and x_3 are on different scales). Then
the gam formula
y ~ s(x1,id=1) + te(x_2,x3) + s(x4,id=1)
would achieve the desired result. id can be numbers or character strings. Giving an id to a
term with a factor by variable causes the smooths at each level of the factor to have the same smoothing
parameter.
Smooth term ids are not supported by gamm.
General linear functional terms have a long history in the spline literature including in the penalized
GLM context (see e.g. Wahba (sC,M
q8sUh G@
1@Nh`|TN6X+.
ujc"Wb 1Ailۮ_w ~[M"{A)i 7֔:@u!k z.=w5ES 4:KCpK8\+h4g||J[H`E* Faq QL+7п碵p
7NHT위..w.!4L
j88It0~ǵ* JZ%7b_<<8*(zlU{
N(mͬETF Y]VpHIvs[ʼnwH e^&T camP :h+Xx3RZ`{/o>ḑȌݍ6HB*ey]twԉnϗL37ƥl=XևQ?FKls3$8]ǴG0nÚG9/\i @ ƛ6eq 1~?w'L7SC)w/s"p]>k~O+"o&-Û縑ptwsWOo<5T 9ӻf3\E=Xqoa#A=*.jBN,QWeBY/$=
[#fYm{.K&
ɂ5,DZ4U3;.3 㴅33+4ZfD-z_½ѳGD_NbNDl5(Jh?3=FX2O,gPt8mF1aI4ܾi5E]7ۮ'