Fit smoothing spline ANOVA models in non-Gaussian regression. The
symbolic model specification via formula follows the same
rules as in lm and glm.
gssanova0(formula, family, type=NULL, data=list(), weights, subset,
offset, na.action=na.omit, partial=NULL, method=NULL,
varht=1, nu=NULL, prec=1e-7, maxiter=30)gssanova1(formula, family, type=NULL, data=list(), weights, subset,
offset, na.action=na.omit, partial=NULL, method=NULL,
varht=1, alpha=1.4, nu=NULL, id.basis=NULL, nbasis=NULL,
seed=NULL, random=NULL, skip.iter=FALSE)
gssanova0 returns a list object of class
c("gssanova0","ssanova0","gssanova").
gssanova1 returns a list object of class
c("gssanova","ssanova").
The method summary.gssanova0 or
summary.gssanova can be used to obtain summaries of
the fits. The method predict.ssanova0 or
predict.ssanova can be used to evaluate the fits at
arbitrary points along with standard errors, on the link scale. The
methods residuals.gssanova and
fitted.gssanova extract the respective traits from the
fits.
Symbolic description of the model to be fit.
Description of the error distribution. Supported
are exponential families "binomial", "poisson",
"Gamma", "inverse.gaussian", and
"nbinomial". Also supported are accelerated life model
families "weibull", "lognorm", and
"loglogis". Further more, proportional odds logistic
regression "polr" for ordinal response is also
supported.
List specifying the type of spline for each variable.
See mkterm for details.
Optional data frame containing the variables in the model.
Optional vector of weights to be used in the fitting process.
Optional vector specifying a subset of observations to be used in the fitting process.
Optional offset term with known parameter 1.
Function which indicates what should happen when the data contain NAs.
Optional symbolic description of parametric terms in partial spline models.
Score used to drive the performance-oriented
iteration. Supported are method="v" for GCV,
method="m" for GML, and method="u" for Mallows' CL.
Dispersion parameter needed for method="u".
Ignored when method="v" or method="m" are
specified.
Inverse scale parameter in accelerated life model families. Ignored for exponential families.
Precision requirement for the iterations.
Maximum number of iterations allowed for performance-oriented iteration, and for inner-loop multiple smoothing parameter selection when applicable.
Tuning parameter modifying GCV or Mallows' CL.
Index designating selected "knots".
Number of "knots" to be selected. Ignored when
id.basis is supplied.
Seed for reproducible random selection of "knots".
Ignored when id.basis is supplied.
Input for parametric random effects in nonparametric
mixed-effect models. See mkran for details.
Flag indicating whether to use initial values of
theta and skip theta iteration. See ssanova for
notes on skipping theta iteration.
For family="binomial", the response can be specified either
as two columns of counts or as a column of sample proportions plus a
column of total counts entered through the argument weights,
as in glm.
For family="nbinomial", the response may be specified as two
columns with the second being the known sizes, or simply as a single
column with the common unknown size to be estimated through the
maximum likelihood.
For family="weibull", "lognorm", or "loglogis",
the response consists of three columns, with the first giving the
follow-up time, the second the censoring status, and the third the
left-truncation time. For data with no truncation, the third column
can be omitted.
For family="polr", the response should be an ordered factor.
The model specification via formula is intuitive. For
example, y~x1*x2 yields a model of the form
$$
y = C + f_{1}(x1) + f_{2}(x2) + f_{12}(x1,x2) + e
$$
with the terms denoted by "1", "x1", "x2", and
"x1:x2".
The model terms are sums of unpenalized and penalized terms. Attached to every penalized term there is a smoothing parameter, and the model complexity is largely determined by the number of smoothing parameters.
Only one link is implemented for each family. It is the
logit link for "binomial", and the log link for
"poisson", "Gamma", and "inverse.gaussian".
For "nbinomial", the working parameter is the logit of the
probability \(p\); see NegBinomial. For
"weibull", "lognorm", and "loglogis", it is the
location parameter for the log lifetime.
The models are fitted by penalized likelihood method through the
performance-oriented iteration as described in the reference. For
family="binomial", "poisson", "nbinomial",
"weibull", "lognorm", and "loglogis", the score
driving the performance-oriented iteration defaults to
method="u" with varht=1. For family="Gamma"
and "inverse.gaussian", the default is method="v".
gssanova0 uses the algorithm of ssanova0 for
the iterated penalized least squares problems, whereas
gssanova1 uses the algorithm of ssanova.
In gssanova1, a subset of the observations are selected as
"knots." Unless specified via id.basis or nbasis, the
number of "knots" \(q\) is determined by \(max(30,10n^{2/9})\),
which is appropriate for the default cubic splines for numerical
vectors.
Gu, C. (1992), Cross-validating non Gaussian data. Journal of Computational and Graphical Statistics, 1, 169-179.
Gu, C. (2013), Smoothing Spline ANOVA Models (2nd Ed). New York: Springer-Verlag.
GU, C. (2014), Smoothing Spline ANOVA Models: R Package gss. Journal of Statistical Software, 58(5), 1-25. URL http://www.jstatsoft.org/v58/i05/.
## Fit a cubic smoothing spline logistic regression model
test <- function(x)
{.3*(1e6*(x^11*(1-x)^6)+1e4*(x^3*(1-x)^10))-2}
x <- (0:100)/100
p <- 1-1/(1+exp(test(x)))
y <- rbinom(x,3,p)
logit.fit <- gssanova0(cbind(y,3-y)~x,family="binomial")
## The same fit
logit.fit1 <- gssanova0(y/3~x,"binomial",weights=rep(3,101))
## Obtain estimates and standard errors on a grid
est <- predict(logit.fit,data.frame(x=x),se=TRUE)
## Plot the fit and the Bayesian confidence intervals
plot(x,y/3,ylab="p")
lines(x,p,col=1)
lines(x,1-1/(1+exp(est$fit)),col=2)
lines(x,1-1/(1+exp(est$fit+1.96*est$se)),col=3)
lines(x,1-1/(1+exp(est$fit-1.96*est$se)),col=3)
## Clean up
if (FALSE) rm(test,x,p,y,logit.fit,logit.fit1,est)
dev.off()
Run the code above in your browser using DataLab