| Title: | Shape Constrained Additive Models |
|---|---|
| Description: | Generalized additive models under shape constraints on the component functions of the linear predictor. Models can include multiple shape-constrained (univariate and bivariate) and unconstrained terms. Routines of the package 'mgcv' are used to set up the model matrix, print, and plot the results. Multiple smoothing parameter estimation by the Generalized Cross Validation or similar. See Pya and Wood (2015) <doi:10.1007/s11222-013-9448-7> for an overview. A broad selection of shape-constrained smoothers, linear functionals of smooths with shape constraints, and Gaussian models with AR1 residuals. |
| Authors: | Natalya Pya <[email protected]> |
| Maintainer: | Natalya Pya <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 1.2-17 |
| Built: | 2024-12-17 06:31:44 UTC |
| Source: | CRAN |
scam provides functions for generalized additive modelling under shape constraints on the component functions of the linear predictor of the GAM. Models can contain multiple univariate and bivariate shape constrained terms, unconstrained terms and parametric terms. A wide variety of shape constrained smooths covered in shape.constrained.smooth.terms are provided.
The model set-up is similar to that of gam() of the package mgcv, so unconstrained smooths of one or more variables of the mgcv can be included in SCAMs. User-defined smooths can be added as well. SCAM is estimated by penalized log likelihood maximization and provides automatic smoothness selection by minimizing generalized cross validation or similar. A Bayesian approach is used to obtain a covariance matrix of the model coefficients and credible intervals for each smooth. Linear functionals of smooth functions with shape constraints, parametric model terms, simple linear random effects terms, bivariate interaction smooths with increasing/decreasing constraints (smooth ANOVA), and identity link Gaussian models with AR1 residuals are supported.
scam provides generalized additive modelling under shape constraints functions scam, summary.scam, plot.scam, scam.check, predict.scam, anova.scam, and vis.scam. These are based on the functions of the unconstrained GAM of the package mgcv and are similar in use.
The use of scam() is much like the use of gam(), except that within a scam model formula, shape constrained smooths of one or two predictors can be specified using s terms with a type of shape constraints used specified as a letter character string of the argument bs, e.g. s(x, bs="mpi") for smooth subject to increasing constraint. See shape.constrained.smooth.terms for a complete overview of what is available. scam model estimation is performed by penalized likelihood maximization, with smoothness selection by GCV, UBRE/AIC criteria. See scam, linear.functional.terms for a short discussion of model specification and some examples. See scam arguments optimizer and optim.method, and scam.control for detailed control of scam model fitting. For checking and visualization, see scam.check, plot.scam, qq.scam and vis.scam. For extracting fitting results, see summary.scam and anova.scam.
A Bayesian approach to smooth modelling is used to obtain covariance matrix of the model coefficients and credible intervals for each smooth. Vp element of the fitted object of class scam returns the Bayesian covariance matrix, Ve returns the frequentist estimated covariance matrix for the parameter estimators. The
frequentist estimated covariance matrix for the reparametrized parameter estimators (obtained using the delta method) is returned in Ve.t, which is particularly useful for testing individual smooth
terms for equality to the zero function (not so useful for CI's as smooths are usually biased).
Vp.t returns the Bayesian covariance matrix for the reparametrized parameters. Frequentist approximations can be used for hypothesis testing based on model comparison; see anova.scam and summary.scam for info on hypothesis testing.
For a complete list of functions type library(help=scam).
Natalya Pya <[email protected]> based partly on mgcv by Simon Wood
Maintainer: Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
Wood S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapman and Hall/CRC Press
Wood, S.N. (2008) Fast stable direct fitting and smoothness selection for generalized additive models. Journal of the Royal Statistical Society (B) 70(3):495-518
Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1):3-36
The package was part supported by EPSRC grants EP/I000917/1, EP/K005251/1 and the Science Committee of the Ministry of Science and Education of the Republic of Kazakhstan grant #2532/GF3.
## see examples for scam## see examples for scam
Performs hypothesis tests relating to one or more fitted scam objects.
The function is a clone of anova.gam of the mgcv package.
The documentation below is similar to that of object anova.gam.
## S3 method for class 'scam' anova(object, ..., dispersion = NULL, test = NULL, freq = FALSE,p.type=0) ## S3 method for class 'anova.scam' print(x, digits = max(3, getOption("digits") - 3),...)## S3 method for class 'scam' anova(object, ..., dispersion = NULL, test = NULL, freq = FALSE,p.type=0) ## S3 method for class 'anova.scam' print(x, digits = max(3, getOption("digits") - 3),...)
object, ...
|
fitted model objects of class |
x |
an |
dispersion |
a value for the dispersion parameter: not normally used. |
test |
what sort of test to perform for a multi-model call. One of
|
freq |
whether to use frequentist or Bayesian approximations for parametric term
p-values. See |
p.type |
selects exact test statistic to use for single smooth term p-values. See
|
digits |
number of digits to use when printing output. |
see anova.gam for details.
In the multi-model case anova.scam produces output identical to
anova.glm, which it in fact uses.
In the single model case an object of class anova.scam is produced,
which is in fact an object returned from summary.scam.
print.anova.scam simply produces tabulated output.
If models 'a' and 'b' differ only in terms with no un-penalized components then p values from anova(a,b) are unreliable, and usually much too low.
Default P-values will usually be wrong for parametric terms penalized using ‘paraPen’: use freq=TRUE to obtain better p-values when the penalties are full rank and represent conventional random effects.
For a single model, interpretation is similar to drop1, not anova.lm.
Simon N. Wood [email protected]
Scheipl, F., Greven, S. and Kuchenhoff, H. (2008) Size and power of tests for a zero random effect variance or polynomial regression in additive and linear mixed models. Comp. Statist. Data Anal. 52, 3283-3299
Wood, S.N. (2013a) On p-values for smooth components of an extended generalized additive model. Biometrika 100:221-228
Wood, S.N. (2013b) A simple test for random effects in regression models. Biometrika 100:1005-1010
scam, predict.scam,
scam.check, summary.scam, anova.gam
library(scam) set.seed(0) fac <- rep(1:4,20) x1 <- runif(80)*5 x2 <- runif(80,-1,2) x3 <- runif(80, 0, 1) y <- fac+log(x1)/5 y <- y + exp(-1.3*x2) +rnorm(80)*0.1 fac <- factor(fac) b <- scam(y ~ fac+s(x1,bs="mpi")+s(x2,bs="mpd")+s(x3)) b1 <- scam(y ~ fac+s(x1,bs="mpi")+s(x2,bs="mpd")) anova(b,b1,test="F") ## b2 <- scam(y ~ fac +s(x1)+s(x2)+te(x1,x2))library(scam) set.seed(0) fac <- rep(1:4,20) x1 <- runif(80)*5 x2 <- runif(80,-1,2) x3 <- runif(80, 0, 1) y <- fac+log(x1)/5 y <- y + exp(-1.3*x2) +rnorm(80)*0.1 fac <- factor(fac) b <- scam(y ~ fac+s(x1,bs="mpi")+s(x2,bs="mpd")+s(x3)) b1 <- scam(y ~ fac+s(x1,bs="mpi")+s(x2,bs="mpd")) anova(b,b1,test="F") ## b2 <- scam(y ~ fac +s(x1)+s(x2)+te(x1,x2))
Function to efficiently estimate smoothing parameters of SCAM by GCV/UBRE score optimization. The procedure is outer to the model fitting by the Newton-Raphson method. The function uses the BFGS method where the Hessian matrix is updated iteratively at each step. Backtracking is included to satisfy the sufficient decrease condition.
The function is not normally called directly, but rather service routines for scam.
bfgs_gcv.ubre(fn=gcv.ubre_grad, rho, ini.fd=TRUE, G, env, n.pen=length(rho), typx=rep(1,n.pen), typf=1, control)bfgs_gcv.ubre(fn=gcv.ubre_grad, rho, ini.fd=TRUE, G, env, n.pen=length(rho), typx=rep(1,n.pen), typf=1, control)
fn |
GCV/UBRE Function which returs the GCV/UBRE value and its derivative wrt log smoothing parameter. |
rho |
log of the initial values of the smoothing parameters. |
ini.fd |
If TRUE, a finite difference to the Hessian is used to find the initial inverse Hessian, otherwise the initial inverse Hessian is a diagonal matrix ‘100*I’. |
G |
A list of items needed to fit a SCAM. |
env |
Get the enviroment for the model coefficients, their derivatives and the smoothing parameter. |
n.pen |
Smoothing parameter dimension. |
typx |
A vector whose component is a positive scalar specifying the typical magnitude of sp. |
typf |
A positive scalar estimating the magnitude of the gcv near the minimum. |
control |
Control option list as returned by |
A list is returned with the following items:
gcv.ubre |
The optimal value of GCV/UBRE. |
rho |
The best value of the log smoothing parameter. |
dgcv.ubre |
The gradient of the GCV/UBRE. |
iterations |
The number of iterations taken until convergence. |
conv.bfgs |
Convergence information indicating why the BFGS terminated (given below). |
termcode |
An integer code indicating why the optimization process terminated. 1: relative gradient is close to zero, current iterate probably is a solution. 2: scaled distance between last two steps less than ‘steptol’, current iterate probably is a local minimizer, but it's possible that the algorithm is making very slow progress, or ‘steptol’ is too large. 3: last global step failed to locate a point lower than
estimate. Either estimate is an approximate
local minimum of the function or 4: iteration limit exceeded. 5: five consecutive steps of length |
object |
A list of elements returned by the fitting procedure |
dgcv.ubre.check |
If |
check.grad |
If |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society: Series B. 73(1): 1–34
This function calculates the finite-difference approximation of the GCV/UBRE gradient for the fitted model and compares it with the analytical gradient.
check.analytical(object, data, del=1e-6,control)check.analytical(object, data, del=1e-6,control)
object |
A fitted |
data |
An original data frame or list containing the model response variable and covariates. |
del |
A positive scalar (default is 1e-6) giving an increment for finite difference approximation. |
control |
Control option list as returned by |
A list is returned with the following items:
dgcv.ubre.fd |
The finite-difference approximation of the gradient. |
check.grad |
The relative difference in percentage between the analytical and finite differenced derivatives. |
Natalya Pya <[email protected]>
Function to get derivatives of the smooth model terms (currently only of the univariate smooths). Analytical derivatives for SCOP-splines (shape constrained P-splines), finite difference approximation is used for all others
derivative.scam(object,smooth.number=1,deriv=1)derivative.scam(object,smooth.number=1,deriv=1)
object |
fitted scam object |
smooth.number |
ordered number of the smooth model term (1,2,...), ordered as in the formula, which derivative is needed to be calculated. |
deriv |
either 1 if the 1st derivative is required, or 2 if the 2nd |
d |
values of the derivative of the smooth term. |
se.d |
standard errors of the derivative. |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
set.seed(2) n <- 200 x1 <- runif(n)*4-1; f1 <- exp(4*x1)/(1+exp(4*x1)) # monotone increasing smooth x2 <- sort(runif(n)*3-1) # decreasing smooth f2 <- exp(-1.3*x2) f <- f1+ f2 y <- f+ rnorm(n)*0.2 ## fit model, results, and plot... b <- scam(y~ s(x1,k=20,bs="mpi")+s(x2,k=15,bs="mpd")) d1 <- derivative.scam(b,smooth.number=1,deriv=1) par(mfrow=c(1,2)) xx <- sort(x1,index=TRUE) plot(xx$x,d1$d[xx$ix],type="l",xlab=expression(x[1]), ylab=expression(df[1]/dx[1])) d2 <- derivative.scam(b,smooth.number=2,deriv=1) xx <- sort(x2,index=TRUE) plot(xx$x,d2$d[xx$ix],type="l",xlab=expression(x[2]), ylab=expression(df[2]/dx[2]))set.seed(2) n <- 200 x1 <- runif(n)*4-1; f1 <- exp(4*x1)/(1+exp(4*x1)) # monotone increasing smooth x2 <- sort(runif(n)*3-1) # decreasing smooth f2 <- exp(-1.3*x2) f <- f1+ f2 y <- f+ rnorm(n)*0.2 ## fit model, results, and plot... b <- scam(y~ s(x1,k=20,bs="mpi")+s(x2,k=15,bs="mpd")) d1 <- derivative.scam(b,smooth.number=1,deriv=1) par(mfrow=c(1,2)) xx <- sort(x1,index=TRUE) plot(xx$x,d1$d[xx$ix],type="l",xlab=expression(x[1]), ylab=expression(df[1]/dx[1])) d2 <- derivative.scam(b,smooth.number=2,deriv=1) xx <- sort(x2,index=TRUE) plot(xx$x,d2$d[xx$ix],type="l",xlab=expression(x[2]), ylab=expression(df[2]/dx[2]))
Description of scam formula (see gam of the mgcv package for Details), and how to extract it from a fitted scam object.
The function is a clone of formula.gam of the mgcv package.
## S3 method for class 'scam' formula(x,...)## S3 method for class 'scam' formula(x,...)
x |
fitted model objects of class |
... |
un-used in this case |
see formula.gam for details.
Returns the model formula, x$formula. Provided so that anova methods
print an appropriate description of the model.
For the estimation of the SCAM smoothing parameters the GCV/UBRE score is optimized outer to the Newton-Raphson procedure of the model fitting. This function returns the value of the GCV/UBRE score and calculates its first derivative with respect to the log smoothing parameter using the method of Wood (2009).
The function is not normally called directly, but rather service routines for bfgs_gcv.ubre.
gcv.ubre_grad(rho, G, env, control)gcv.ubre_grad(rho, G, env, control)
rho |
log of the initial values of the smoothing parameters. |
G |
a list of items needed to fit a SCAM. |
env |
Get the enviroment for the model coefficients, their derivatives and the smoothing parameter. |
control |
A list of fit control parameters as returned by |
A list is returned with the following items:
dgcv.ubre |
The value of GCV/UBRE gradient. |
gcv.ubre |
The GCV/UBRE score value. |
scale.est |
The value of the scale estimate. |
object |
The elements of the fitting procedure |
dgcv.ubre.check |
If |
check.grad |
If |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
Wood S.N. (2006) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press.
Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society: Series B. 73(1): 1–34
Since scam uses the model setup of gam of the mgcv package, in the same way as in gam scam allows the response variable to depend on linear
functionals of smooth terms in the s with additional shape constraints.
See linear.functional.terms(mgcv).
## Not run: ########################################### ## similar to a "signal" regression ## example from mgcv() ... ########################################### library(scam) ## decreasing smooth... set.seed(4) rf <- function(x=seq(-1,3,length=100)) { ## taken from mgcv package ## generates random functions... m <- ceiling(runif(1)*5) ## number of components f <- x*0; mu <- runif(m,min(x),max(x)); sig <- (runif(m)+.5)*(max(x)-min(x))/10 for (i in 1:m) f <- f+ dnorm(x,mu[i],sig[i]) f } ## simulate 200 functions and store in rows of L... L <- matrix(NA,200,100) for (i in 1:200) L[i,] <- rf() ## simulate the functional predictors x <- seq(-1,3,length=100) ## evaluation points f2 <- function(x) { ## the coefficient function -4*exp(4*x)/(1+exp(4*x)) } f <- f2(x) plot(x,f ,type="l") y <- L%*%f + rnorm(200)*20 ## simulated response data X <- matrix(x,200,100,byrow=TRUE) b <- scam(y~s(X,by=L,k=20,bs="mpdBy")) par(mfrow=c(1,2)) plot(b,shade=TRUE);lines(x,f,col=2); ## compare with gam() of mgcv package... g <- gam(y~s(X,by=L,k=20)) plot(g,shade=TRUE);lines(x,f,col=2) ## increasing smooth.... L <- matrix(NA,200,100) for (i in 1:200) L[i,] <- rf() ## simulate the functional predictors x <- seq(-1,3,length=100) ## evaluation points f2 <- function(x) { ## the coefficient function 4*exp(4*x)/(1+exp(4*x)) } f <- f2(x) plot(x,f ,type="l") y <- L%*%f + rnorm(200)*20 ## simulated response data X <- matrix(x,200,100,byrow=TRUE) b <- scam(y~s(X,by=L,k=20,bs="mpiBy")) par(mfrow=c(1,2)) plot(b,shade=TRUE);lines(x,f,col=2); ## compare with unconstrained fit... g <- scam(y~s(X,by=L,k=20)) plot(g,shade=TRUE);lines(x,f,col=2) ## convex smooth... ## simulate 200 functions and store in rows of L... set.seed(4) L <- matrix(NA,200,100) for (i in 1:200) L[i,] <- rf(x=sort(2*runif(100)-1)) ## simulate the functional predictors x <- sort(runif(100,-1,1)) ## evaluation points f2 <- function(x){4*x^2 } ## the coefficient function f <- f2(x) plot(x,f ,type="l") y <- L%*%f + rnorm(200)*30 ## simulated response data X <- matrix(x,200,100,byrow=TRUE) b <- scam(y~s(X,by=L,k=20,bs="cxBy")) par(mfrow=c(1,2)) plot(b,shade=TRUE);lines(x,f,col=2); g <- scam(y~s(X,by=L,k=20)) plot(g,shade=TRUE);lines(x,f,col=2) ## End(Not run)## Not run: ########################################### ## similar to a "signal" regression ## example from mgcv() ... ########################################### library(scam) ## decreasing smooth... set.seed(4) rf <- function(x=seq(-1,3,length=100)) { ## taken from mgcv package ## generates random functions... m <- ceiling(runif(1)*5) ## number of components f <- x*0; mu <- runif(m,min(x),max(x)); sig <- (runif(m)+.5)*(max(x)-min(x))/10 for (i in 1:m) f <- f+ dnorm(x,mu[i],sig[i]) f } ## simulate 200 functions and store in rows of L... L <- matrix(NA,200,100) for (i in 1:200) L[i,] <- rf() ## simulate the functional predictors x <- seq(-1,3,length=100) ## evaluation points f2 <- function(x) { ## the coefficient function -4*exp(4*x)/(1+exp(4*x)) } f <- f2(x) plot(x,f ,type="l") y <- L%*%f + rnorm(200)*20 ## simulated response data X <- matrix(x,200,100,byrow=TRUE) b <- scam(y~s(X,by=L,k=20,bs="mpdBy")) par(mfrow=c(1,2)) plot(b,shade=TRUE);lines(x,f,col=2); ## compare with gam() of mgcv package... g <- gam(y~s(X,by=L,k=20)) plot(g,shade=TRUE);lines(x,f,col=2) ## increasing smooth.... L <- matrix(NA,200,100) for (i in 1:200) L[i,] <- rf() ## simulate the functional predictors x <- seq(-1,3,length=100) ## evaluation points f2 <- function(x) { ## the coefficient function 4*exp(4*x)/(1+exp(4*x)) } f <- f2(x) plot(x,f ,type="l") y <- L%*%f + rnorm(200)*20 ## simulated response data X <- matrix(x,200,100,byrow=TRUE) b <- scam(y~s(X,by=L,k=20,bs="mpiBy")) par(mfrow=c(1,2)) plot(b,shade=TRUE);lines(x,f,col=2); ## compare with unconstrained fit... g <- scam(y~s(X,by=L,k=20)) plot(g,shade=TRUE);lines(x,f,col=2) ## convex smooth... ## simulate 200 functions and store in rows of L... set.seed(4) L <- matrix(NA,200,100) for (i in 1:200) L[i,] <- rf(x=sort(2*runif(100)-1)) ## simulate the functional predictors x <- sort(runif(100,-1,1)) ## evaluation points f2 <- function(x){4*x^2 } ## the coefficient function f <- f2(x) plot(x,f ,type="l") y <- L%*%f + rnorm(200)*30 ## simulated response data X <- matrix(x,200,100,byrow=TRUE) b <- scam(y~s(X,by=L,k=20,bs="cxBy")) par(mfrow=c(1,2)) plot(b,shade=TRUE);lines(x,f,col=2); g <- scam(y~s(X,by=L,k=20)) plot(g,shade=TRUE);lines(x,f,col=2) ## End(Not run)
Function to extract the log-likelihood for a fitted scam
model (fitted by penalized likelihood maximization).
Used by AIC.
The function is a clone of logLik.gam of the mgcv package.
The documentation below is similar to that of object logLik.gam.
## S3 method for class 'scam' logLik(object,...)## S3 method for class 'scam' logLik(object,...)
object |
fitted model objects of class |
... |
unused in this case |
see logLik.gam for details.
Standard logLik object: see logLik.
Hastie and Tibshirani, 1990, Generalized Additive Models.
Wood, S.N. (2008) Fast stable direct fitting and smoothness selection for generalized additive models. J.R.Statist. Soc. B 70(3):495-518
This function returns the marginal model matrices and the list of penalty matrices for the tensor product bivariate smooth with the single concavity or convexity restriction along the second covariate. The marginal smooth functions of both covariates are constructed using the B-spline basis functions.
marginal.matrices.tescv.ps(x, z, xk, zk, m, q1, q2)marginal.matrices.tescv.ps(x, z, xk, zk, m, q1, q2)
x |
A numeric vector of the values of the first covariate at which to evaluate the B-spline marginal functions.
The values in |
z |
A numeric vector of the values of the second covariate at which to evaluate the B-spline marginal functions.
The values in |
xk |
A numeric vector of knot positions for the first covariate, |
zk |
A numeric vector of knot positions for the second covariate, |
m |
A pair of two numbers where |
q1 |
A number denoting the basis dimension of the first marginal smooth. |
q2 |
A number denoting the basis dimension of the second marginal smooth. |
The function is not called directly, but is rather used internally by the constructor
smooth.construct.tescv.smooth.spec and smooth.construct.tescx.smooth.spec .
X1 |
Marginal model matrix for the first unconstrained marginal smooth. |
X2 |
Marginal model matrix for the second monotonic marginal smooth. |
S |
A list of penalty matrices for this tensor product smooth. |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
smooth.construct.tescv.smooth.spec,
smooth.construct.tescx.smooth.spec,
marginal.matrices.tesmi1.ps,
smooth.construct.tesmd1.smooth.spec,
smooth.construct.tesmd2.smooth.spec
This function returns the marginal model matrices and the list of penalty matrices for the tensor product bivariate smooth with the single monotone increasing or decreasing restriction along the first covariate. The marginal smooth functions of both covariates are constructed using the B-spline basis functions.
marginal.matrices.tesmi1.ps(x, z, xk, zk, m, q1, q2)marginal.matrices.tesmi1.ps(x, z, xk, zk, m, q1, q2)
x |
A numeric vector of the values of the first covariate at which to evaluate the B-spline marginal functions.
The values in |
z |
A numeric vector of the values of the second covariate at which to evaluate the B-spline marginal functions.
The values in |
xk |
A numeric vector of knot positions for the first covariate, |
zk |
A numeric vector of knot positions for the second covariate, |
m |
A pair of two numbers where |
q1 |
A number denoting the basis dimension of the first marginal smooth. |
q2 |
A number denoting the basis dimension of the second marginal smooth. |
The function is not called directly, but is rather used internally by the constructor
smooth.construct.tesmi1.smooth.spec and smooth.construct.tesmd1.smooth.spec .
X1 |
Marginal model matrix for the first monotonic marginal smooth. |
X2 |
Marginal model matrix for the second unconstrained marginal smooth. |
S |
A list of penalty matrices for this tensor product smooth. |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
Wood S.N. (2006) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press.
smooth.construct.tesmi1.smooth.spec,
smooth.construct.tesmi2.smooth.spec,
marginal.matrices.tesmi2.ps,
smooth.construct.tesmd1.smooth.spec,
smooth.construct.tesmd2.smooth.spec
This function returns the marginal model matrices and the list of penalty matrices for the tensor product bivariate smooth with the single monotone increasing or decreasing restriction along the second covariate. The marginal smooth functions of both covariates are constructed using the B-spline basis functions.
marginal.matrices.tesmi2.ps(x, z, xk, zk, m, q1, q2)marginal.matrices.tesmi2.ps(x, z, xk, zk, m, q1, q2)
x |
A numeric vector of the values of the first covariate at which to evaluate the B-spline marginal functions.
The values in |
z |
A numeric vector of the values of the second covariate at which to evaluate the B-spline marginal functions.
The values in |
xk |
A numeric vector of knot positions for the first covariate, |
zk |
A numeric vector of knot positions for the second covariate, |
m |
A pair of two numbers where |
q1 |
A number denoting the basis dimension of the first marginal smooth. |
q2 |
A number denoting the basis dimension of the second marginal smooth. |
The function is not called directly, but is rather used internally by the constructor
smooth.construct.tesmi2.smooth.spec and
smooth.construct.tesmd2.smooth.spec .
X1 |
Marginal model matrix for the first unconstrained marginal smooth. |
X2 |
Marginal model matrix for the second monotonic marginal smooth. |
S |
A list of penalty matrices for this tensor product smooth. |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
Wood S.N. (2006) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press.
smooth.construct.tesmi1.smooth.spec,
smooth.construct.tesmi2.smooth.spec,
marginal.matrices.tesmi1.ps,
smooth.construct.tesmd1.smooth.spec,
smooth.construct.tesmd2.smooth.spec
The function is a clone of the plot.gam of the mgcv package with the differences
in the construction of the Bayesian confidence intervals of the shape constrained smooth terms. The function
takes a fitted scam object produced by scam() and plots the
component smooth functions that make it up, on the scale of the linear
predictor. Optionally produces term plots for parametric model components
as well.
Note: The fitted shape constrained smooth functions are centred when plotted, which is done in order to be in line with plots of unconstrained smooths (as in gam()). Although 'zeroed intercept' constraints are applied to deal with identifiability of the scop-splines.
## S3 method for class 'scam' plot(x,residuals=FALSE,rug=TRUE,se=TRUE,pages=0,select=NULL,scale=-1, n=100,n2=40,pers=FALSE,theta=30,phi=30,jit=FALSE,xlab=NULL, ylab=NULL,main=NULL,ylim=NULL,xlim=NULL,too.far=0.1, all.terms=FALSE,shade=FALSE,shade.col="gray80", shift=0,trans=I,seWithMean=FALSE,unconditional = FALSE, by.resids = FALSE,scheme=0,...)## S3 method for class 'scam' plot(x,residuals=FALSE,rug=TRUE,se=TRUE,pages=0,select=NULL,scale=-1, n=100,n2=40,pers=FALSE,theta=30,phi=30,jit=FALSE,xlab=NULL, ylab=NULL,main=NULL,ylim=NULL,xlim=NULL,too.far=0.1, all.terms=FALSE,shade=FALSE,shade.col="gray80", shift=0,trans=I,seWithMean=FALSE,unconditional = FALSE, by.resids = FALSE,scheme=0,...)
The list of the arguments is the same as in plot.gam of the mgcv package.
x |
a fitted |
residuals |
If |
rug |
when TRUE (default) then the covariate to which the plot applies is displayed as a rug plot at the foot of each plot of a 1-d smooth, and the locations of the covariates are plotted as points on the contour plot representing a 2-d smooth. |
se |
when TRUE (default) upper and lower lines are added to the
1-d plots at 2 standard errors
above and below the estimate of the smooth being plotted while for
2-d plots, surfaces at +1 and -1 standard errors are contoured
and overlayed on the contour plot for the estimate. If a
positive number is supplied then this number is multiplied by
the standard errors when calculating standard error curves or
surfaces. See also |
pages |
(default 0) the number of pages over which to spread the output. For example,
if |
select |
Allows the plot for a single model term to be selected for printing. e.g. if you just want the plot for the second smooth term set |
scale |
set to -1 (default) to have the same y-axis scale for each plot, and to 0 for a
different y axis for each plot. Ignored if |
n |
number of points used for each 1-d plot - for a nice smooth plot this needs to be several times the estimated degrees of freedom for the smooth. Default value 100. |
n2 |
Square root of number of points used to grid estimates of 2-d functions for contouring. |
pers |
Set to |
theta |
One of the perspective plot angles. |
phi |
The other perspective plot angle. |
jit |
Set to TRUE if you want rug plots for 1-d terms to be jittered. |
xlab |
If supplied then this will be used as the x label for all plots. |
ylab |
If supplied then this will be used as the y label for all plots. |
main |
Used as title (or z axis label) for plots if supplied. |
ylim |
If supplied then this pair of numbers are used as the y limits for each plot. |
xlim |
If supplied then this pair of numbers are used as the x limits for each plot. |
too.far |
If greater than 0 then this is used to determine when a location is too
far from data to be plotted when plotting 2-D smooths. This is useful since smooths tend to go wild away from data.
The data are scaled into the unit square before deciding what to exclude, and |
all.terms |
if set to |
shade |
Set to |
shade.col |
define the color used for shading confidence bands. |
shift |
constant to add to each smooth (on the scale of the linear
predictor) before plotting. Can be useful for some diagnostics, or with |
trans |
function to apply to each smooth (after any shift), before
plotting. |
seWithMean |
if |
unconditional |
if |
by.resids |
Should partial residuals be plotted for terms with |
scheme |
Integer (0,1 or 2) or integer vector selecting a plotting scheme for each plot.
|
... |
other graphics parameters to pass on to plotting commands. |
The function generates plots.
Natalya Pya <[email protected]> based on the plot.gam of the mgcv by Simon Wood
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
Wood S.N. (2006) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press.
## simulating data... require(scam) n <- 200 set.seed(1) x0 <- rep(1:4,50) x1 <- runif(n)*6-3 f1 <- 3*exp(-x1^2) # unconstrained smooth term x2 <- runif(n)*4-1; f2 <- exp(4*x2)/(1+exp(4*x2)) # monotone increasing smooth x3 <- runif(n)*5; f3 <- -log(x3)/5 # monotone decreasing smooth f <- f1+f2+f3 y <- 2*x0 + f + rnorm(n)*.3 x0 <- factor(x0) ## fit the model and plot ... b <- scam(y~x0+s(x1,k=15,bs="cr")+s(x2,k=30,bs="mpi")+s(x3,k=30,bs="mpd")) plot(b,pages=1,residuals=TRUE,all.terms=TRUE,shade=TRUE,shade.col=3) ## same with EFS and BFGS methods for smoothing parameter and models coefficients estimations... b <- scam(y~x0+s(x1,k=15,bs="cr")+s(x2,k=30,bs="mpi")+s(x3,k=30,bs="mpd"),optimizer=c("efs","bfgs")) plot(b,pages=1,residuals=TRUE,all.terms=TRUE,shade=TRUE,shade.col=3) ## Not run: ## example with 2-d plots... ## simulating data... set.seed(2) n <- 30 x0 <- rep(1:9,100) x1 <- sort(runif(n)*4-1) x2 <- sort(runif(n)) x3 <- runif(n*n, 0, 1) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) { f[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i]))+2*sin(pi*x2[j])} f1 <- as.vector(t(f)) f2 <- x3*0 e <- rnorm(length(f1))*.1 y <- 2*x0 + f1 + f2 + e x0 <- factor(x0) x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x0=x0,x1=x11,x2=x22,x3=x3,y=y) ## fit model and plot ... b <- scam(y~x0+s(x1,x2,k=c(10,10),bs=c("tesmd1","ps"),m=2)+s(x3),data=dat,optimizer="efs") op <- par(mfrow=c(2,2)) plot(b,all.terms=TRUE) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data",pch=19,cex=.3) par(op) ## and use of schemes... op <- par(mfrow=c(2,2)) plot(b,all.terms=TRUE,scheme=1) par(op) op <- par(mfrow=c(2,2)) plot(b,all.terms=TRUE,scheme=c(2,1)) par(op) ## End(Not run)## simulating data... require(scam) n <- 200 set.seed(1) x0 <- rep(1:4,50) x1 <- runif(n)*6-3 f1 <- 3*exp(-x1^2) # unconstrained smooth term x2 <- runif(n)*4-1; f2 <- exp(4*x2)/(1+exp(4*x2)) # monotone increasing smooth x3 <- runif(n)*5; f3 <- -log(x3)/5 # monotone decreasing smooth f <- f1+f2+f3 y <- 2*x0 + f + rnorm(n)*.3 x0 <- factor(x0) ## fit the model and plot ... b <- scam(y~x0+s(x1,k=15,bs="cr")+s(x2,k=30,bs="mpi")+s(x3,k=30,bs="mpd")) plot(b,pages=1,residuals=TRUE,all.terms=TRUE,shade=TRUE,shade.col=3) ## same with EFS and BFGS methods for smoothing parameter and models coefficients estimations... b <- scam(y~x0+s(x1,k=15,bs="cr")+s(x2,k=30,bs="mpi")+s(x3,k=30,bs="mpd"),optimizer=c("efs","bfgs")) plot(b,pages=1,residuals=TRUE,all.terms=TRUE,shade=TRUE,shade.col=3) ## Not run: ## example with 2-d plots... ## simulating data... set.seed(2) n <- 30 x0 <- rep(1:9,100) x1 <- sort(runif(n)*4-1) x2 <- sort(runif(n)) x3 <- runif(n*n, 0, 1) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) { f[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i]))+2*sin(pi*x2[j])} f1 <- as.vector(t(f)) f2 <- x3*0 e <- rnorm(length(f1))*.1 y <- 2*x0 + f1 + f2 + e x0 <- factor(x0) x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x0=x0,x1=x11,x2=x22,x3=x3,y=y) ## fit model and plot ... b <- scam(y~x0+s(x1,x2,k=c(10,10),bs=c("tesmd1","ps"),m=2)+s(x3),data=dat,optimizer="efs") op <- par(mfrow=c(2,2)) plot(b,all.terms=TRUE) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data",pch=19,cex=.3) par(op) ## and use of schemes... op <- par(mfrow=c(2,2)) plot(b,all.terms=TRUE,scheme=1) par(op) op <- par(mfrow=c(2,2)) plot(b,all.terms=TRUE,scheme=c(2,1)) par(op) ## End(Not run)
The various built in smooth classes for use with scam have associate
Predict.matrix method functions to enable prediction from the fitted model.
## S3 method for class 'mpi.smooth' Predict.matrix(object, data) ## S3 method for class 'lmpi.smooth' Predict.matrix(object, data) ## S3 method for class 'miso.smooth' Predict.matrix(object, data) ## S3 method for class 'mifo.smooth' Predict.matrix(object, data) ## S3 method for class 'mpd.smooth' Predict.matrix(object, data) ## S3 method for class 'cv.smooth' Predict.matrix(object, data) ## S3 method for class 'cx.smooth' Predict.matrix(object, data) ## S3 method for class 'micx.smooth' Predict.matrix(object, data) ## S3 method for class 'micv.smooth' Predict.matrix(object, data) ## S3 method for class 'mdcx.smooth' Predict.matrix(object, data) ## S3 method for class 'mdcv.smooth' Predict.matrix(object, data) ## S3 method for class 'po.smooth' Predict.matrix(object, data) ## S3 method for class 'mpdBy.smooth' Predict.matrix(object, data) ## S3 method for class 'mpiBy.smooth' Predict.matrix(object, data) ## S3 method for class 'cxBy.smooth' Predict.matrix(object, data) ## S3 method for class 'cvBy.smooth' Predict.matrix(object, data) ## S3 method for class 'mdcxBy.smooth' Predict.matrix(object, data) ## S3 method for class 'mdcvBy.smooth' Predict.matrix(object, data) ## S3 method for class 'micxBy.smooth' Predict.matrix(object, data) ## S3 method for class 'micvBy.smooth' Predict.matrix(object, data) ## S3 method for class 'tedmd.smooth' Predict.matrix(object, data) ## S3 method for class 'tedmi.smooth' Predict.matrix(object, data) ## S3 method for class 'tesmd1.smooth' Predict.matrix(object, data) ## S3 method for class 'tesmd2.smooth' Predict.matrix(object, data) ## S3 method for class 'tesmi1.smooth' Predict.matrix(object, data) ## S3 method for class 'tesmi2.smooth' Predict.matrix(object, data) ## S3 method for class 'temicx.smooth' Predict.matrix(object, data) ## S3 method for class 'temicv.smooth' Predict.matrix(object, data) ## S3 method for class 'tedecx.smooth' Predict.matrix(object, data) ## S3 method for class 'tedecv.smooth' Predict.matrix(object, data) ## S3 method for class 'tescx.smooth' Predict.matrix(object, data) ## S3 method for class 'tescv.smooth' Predict.matrix(object, data) ## S3 method for class 'tecvcv.smooth' Predict.matrix(object, data) ## S3 method for class 'tecxcv.smooth' Predict.matrix(object, data) ## S3 method for class 'tecxcx.smooth' Predict.matrix(object, data) ## S3 method for class 'tismi.smooth' Predict.matrix(object, data) ## S3 method for class 'tismd.smooth' Predict.matrix(object, data)## S3 method for class 'mpi.smooth' Predict.matrix(object, data) ## S3 method for class 'lmpi.smooth' Predict.matrix(object, data) ## S3 method for class 'miso.smooth' Predict.matrix(object, data) ## S3 method for class 'mifo.smooth' Predict.matrix(object, data) ## S3 method for class 'mpd.smooth' Predict.matrix(object, data) ## S3 method for class 'cv.smooth' Predict.matrix(object, data) ## S3 method for class 'cx.smooth' Predict.matrix(object, data) ## S3 method for class 'micx.smooth' Predict.matrix(object, data) ## S3 method for class 'micv.smooth' Predict.matrix(object, data) ## S3 method for class 'mdcx.smooth' Predict.matrix(object, data) ## S3 method for class 'mdcv.smooth' Predict.matrix(object, data) ## S3 method for class 'po.smooth' Predict.matrix(object, data) ## S3 method for class 'mpdBy.smooth' Predict.matrix(object, data) ## S3 method for class 'mpiBy.smooth' Predict.matrix(object, data) ## S3 method for class 'cxBy.smooth' Predict.matrix(object, data) ## S3 method for class 'cvBy.smooth' Predict.matrix(object, data) ## S3 method for class 'mdcxBy.smooth' Predict.matrix(object, data) ## S3 method for class 'mdcvBy.smooth' Predict.matrix(object, data) ## S3 method for class 'micxBy.smooth' Predict.matrix(object, data) ## S3 method for class 'micvBy.smooth' Predict.matrix(object, data) ## S3 method for class 'tedmd.smooth' Predict.matrix(object, data) ## S3 method for class 'tedmi.smooth' Predict.matrix(object, data) ## S3 method for class 'tesmd1.smooth' Predict.matrix(object, data) ## S3 method for class 'tesmd2.smooth' Predict.matrix(object, data) ## S3 method for class 'tesmi1.smooth' Predict.matrix(object, data) ## S3 method for class 'tesmi2.smooth' Predict.matrix(object, data) ## S3 method for class 'temicx.smooth' Predict.matrix(object, data) ## S3 method for class 'temicv.smooth' Predict.matrix(object, data) ## S3 method for class 'tedecx.smooth' Predict.matrix(object, data) ## S3 method for class 'tedecv.smooth' Predict.matrix(object, data) ## S3 method for class 'tescx.smooth' Predict.matrix(object, data) ## S3 method for class 'tescv.smooth' Predict.matrix(object, data) ## S3 method for class 'tecvcv.smooth' Predict.matrix(object, data) ## S3 method for class 'tecxcv.smooth' Predict.matrix(object, data) ## S3 method for class 'tecxcx.smooth' Predict.matrix(object, data) ## S3 method for class 'tismi.smooth' Predict.matrix(object, data) ## S3 method for class 'tismd.smooth' Predict.matrix(object, data)
object |
A smooth object, usually generated by a |
data |
A data frame containing the values of the named covariates at which the smooth term is to be evaluated. |
A matrix mapping the coefficients for the smooth term to its values at the supplied data values.
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
Wood S.N. (2006) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press.
This function is a clone of the mgcv library code predict.gam with some modifications
to adopt shape preserving smooth terms.
It takes a fitted scam object produced by scam()
and produces predictions given a new set of values for the model covariates
or the original values used for the model fit. Predictions can be accompanied
by standard errors, based on the posterior distribution of the model
coefficients.
It now alows prediction outside the range of knots, and use linear extrapolation in this case.
## S3 method for class 'scam' predict(object,newdata,type="link",se.fit=FALSE,terms=NULL,exclude=NULL, block.size=NULL,newdata.guaranteed=FALSE,na.action=na.pass,...)## S3 method for class 'scam' predict(object,newdata,type="link",se.fit=FALSE,terms=NULL,exclude=NULL, block.size=NULL,newdata.guaranteed=FALSE,na.action=na.pass,...)
object |
a fitted |
newdata |
A data frame or list containing the values of the model covariates at which predictions
are required. If this is not provided then predictions corresponding to the
original data are returned. If |
type |
When this has the value |
se.fit |
when this is TRUE (not default) standard error estimates are returned for each prediction. |
terms |
if |
exclude |
if |
block.size |
maximum number of predictions to process per call to underlying code: larger is quicker, but more memory intensive. Set to < 1 to use total number of predictions as this. |
newdata.guaranteed |
Set to |
na.action |
what to do about |
... |
other arguments. |
See predict.gam for details.
If type=="lpmatrix" then a matrix is returned which will
give a vector of linear predictor values (minus any offest) at the supplied covariate
values, when applied to the model coefficient vector.
Otherwise, if se.fit is TRUE then a 2 item list is returned with items (both arrays) fit
and se.fit containing predictions and associated standard error estimates, otherwise an
array of predictions is returned. The dimensions of the returned arrays depends on whether
type is "terms" or not: if it is then the array is 2 dimensional with each
term in the linear predictor separate, otherwise the array is 1 dimensional and contains the
linear predictor/predicted values (or corresponding s.e.s). The linear predictor returned termwise will
not include the offset or the intercept.
newdata can be a data frame, list or model.frame: if it's a model frame
then all variables must be supplied.
Natalya Pya <[email protected]> based partly on mgcv by Simon Wood
Chambers and Hastie (1993) Statistical Models in S. Chapman & Hall.
Wood S.N. (2006) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press.
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
## Not run: library(scam) set.seed(2) n <- 200 x1 <- runif(n)*6-3 f1 <- 3*exp(-x1^2) # unconstrained term x2 <- runif(n)*4-1; f2 <- exp(4*x2)/(1+exp(4*x2)) # monotone increasing smooth f <- f1+f2 y <- f+rnorm(n)*0.2 dat <- data.frame(x1=x1,x2=x2,y=y) b <- scam(y~s(x1,k=15,bs="cr")+s(x2,k=30,bs="mpi"),data=dat) newd <- data.frame(x1=seq(-3,3,length.out=20),x2=seq(-1,3,length.out=20)) pred <- predict(b,newd) pred predict(b,newd,type="terms",se=TRUE) ## prediction on the original data... pr <- predict(b,type="terms") x<-sort(x1,index=TRUE) old.par <- par(mfrow=c(2,2)) plot(x$x,(pr[,1])[x$ix],type="l",col=3,xlab="x1") z<-sort(x2,index=TRUE) plot(z$x,(pr[,2])[z$ix],type="l",col=3,xlab="x2") plot(b,select=1,scale=0,se=FALSE) plot(b,select=2,scale=0,se=FALSE) par(old.par) ## linear extrapolation with predict.scam()... set.seed(3) n <- 100 x <- sort(runif(n)*3-1) f <- exp(-1.3*x) y <- rpois(n,exp(f)) dat <- data.frame(x=x,y=y) b <- scam(y~s(x,k=15,bs="mpd"),family=poisson(link="log"),data=dat) newd <- data.frame(x=c(2.3,2.7,3.2)) fe <- predict(b,newd,type="link",se=TRUE) ylim<- c(min(y,exp(fe$fit)),max(y,exp(fe$fit))) plot(c(x,newd[[1]]),c(y,NA,NA,NA),ylim=ylim,ylab="y",xlab="x") lines(c(x,newd[[1]]),c(b$fitted,exp(fe$fit)),col=3) ## prediction on the original data... pr <- predict(b) plot(x,y) lines(x,exp(pr),col=3) ## Gaussian model .... ## simulating data... set.seed(2) n <- 200 x <- sort(runif(n)*4-1) f <- exp(4*x)/(1+exp(4*x)) # monotone increasing smooth y <- f+rnorm(n)*0.1 dat <- data.frame(x=x,y=y) b <- scam(y~ s(x,k=25,bs="mpi"),data=dat) newd <- data.frame(x=c(3.2,3.3,3.6)) fe <- predict(b,newd) plot(c(x,newd[[1]]),c(y,NA,NA,NA),ylab="y",xlab="x") lines(c(x,newd[[1]]),c(b$fitted,fe),col=3) ### passing observed data + new data... newd <- data.frame(x=c(x,3.2,3.3,3.6)) fe <- predict(b,newd,se=TRUE) plot(newd[[1]],c(y,NA,NA,NA),ylab="y",xlab="x") lines(newd[[1]],fe$fit,col=2) lines(newd[[1]],fe$fit+2*fe$se.fit,col=3) lines(newd[[1]],fe$fit-2*fe$se.fit,col=4) ## prediction with CI... newd <- data.frame(x=seq(-1.2,3.5,length.out=100)) fe <- predict(b,newd,se=TRUE) ylim<- c(min(y,fe$se.fit),max(y,fe$se.fit)) plot(newd[[1]],fe$fit,type="l",ylim=ylim,ylab="y",xlab="x") lines(newd[[1]],fe$fit+2*fe$se.fit,lty=2) lines(newd[[1]],fe$fit-2*fe$se.fit,lty=2) ## prediction on the original data... pr <- predict(b) plot(x,y) lines(x,pr,col=3) ## bivariate example... set.seed(2) n <- 30 x1 <- sort(runif(n)); x2 <- sort(runif(n)*4-1) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- 2*sin(pi*x1[i]) +exp(4*x2[j])/(1+exp(4*x2[j])) f <- as.vector(t(f)); y <- f+rnorm(length(f))*0.1 x11 <- matrix(0,n,n); x11[,1:n] <- x1; x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) b <- scam(y~s(x1,x2,k=c(10,10),bs="tesmi2"),data=dat,optimizer="efs") par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE); plot(b,pers=TRUE,theta = 80, phi = 40) n.out <- 20 xp <- seq(0,1.4,length.out=n.out) zp <- seq(-1,3.4,length.out=n.out) xp1 <- matrix(0,n.out,n.out); xp1[,1:n.out] <- xp xp1 <- as.vector(t(xp1)); xp2 <- rep(zp,n.out) newd <- data.frame(x1=xp1,x2=xp2) fe <- predict(b,newd) fc <- t(matrix(fe,n.out,n.out)) persp(xp,zp,fc,expand= 0.85,ticktype = "simple",xlab="x1", ylab="x2",zlab="f^",main="", theta = 80, phi = 40) ## obtaining a 'prediction matrix'... newd <- data.frame(x1=c(-2,-1),x2=c(0,1)) Xp <- predict(b,newdata=newd,type="lpmatrix") fv <- Xp%*% b$beta.t fv ## End(Not run)## Not run: library(scam) set.seed(2) n <- 200 x1 <- runif(n)*6-3 f1 <- 3*exp(-x1^2) # unconstrained term x2 <- runif(n)*4-1; f2 <- exp(4*x2)/(1+exp(4*x2)) # monotone increasing smooth f <- f1+f2 y <- f+rnorm(n)*0.2 dat <- data.frame(x1=x1,x2=x2,y=y) b <- scam(y~s(x1,k=15,bs="cr")+s(x2,k=30,bs="mpi"),data=dat) newd <- data.frame(x1=seq(-3,3,length.out=20),x2=seq(-1,3,length.out=20)) pred <- predict(b,newd) pred predict(b,newd,type="terms",se=TRUE) ## prediction on the original data... pr <- predict(b,type="terms") x<-sort(x1,index=TRUE) old.par <- par(mfrow=c(2,2)) plot(x$x,(pr[,1])[x$ix],type="l",col=3,xlab="x1") z<-sort(x2,index=TRUE) plot(z$x,(pr[,2])[z$ix],type="l",col=3,xlab="x2") plot(b,select=1,scale=0,se=FALSE) plot(b,select=2,scale=0,se=FALSE) par(old.par) ## linear extrapolation with predict.scam()... set.seed(3) n <- 100 x <- sort(runif(n)*3-1) f <- exp(-1.3*x) y <- rpois(n,exp(f)) dat <- data.frame(x=x,y=y) b <- scam(y~s(x,k=15,bs="mpd"),family=poisson(link="log"),data=dat) newd <- data.frame(x=c(2.3,2.7,3.2)) fe <- predict(b,newd,type="link",se=TRUE) ylim<- c(min(y,exp(fe$fit)),max(y,exp(fe$fit))) plot(c(x,newd[[1]]),c(y,NA,NA,NA),ylim=ylim,ylab="y",xlab="x") lines(c(x,newd[[1]]),c(b$fitted,exp(fe$fit)),col=3) ## prediction on the original data... pr <- predict(b) plot(x,y) lines(x,exp(pr),col=3) ## Gaussian model .... ## simulating data... set.seed(2) n <- 200 x <- sort(runif(n)*4-1) f <- exp(4*x)/(1+exp(4*x)) # monotone increasing smooth y <- f+rnorm(n)*0.1 dat <- data.frame(x=x,y=y) b <- scam(y~ s(x,k=25,bs="mpi"),data=dat) newd <- data.frame(x=c(3.2,3.3,3.6)) fe <- predict(b,newd) plot(c(x,newd[[1]]),c(y,NA,NA,NA),ylab="y",xlab="x") lines(c(x,newd[[1]]),c(b$fitted,fe),col=3) ### passing observed data + new data... newd <- data.frame(x=c(x,3.2,3.3,3.6)) fe <- predict(b,newd,se=TRUE) plot(newd[[1]],c(y,NA,NA,NA),ylab="y",xlab="x") lines(newd[[1]],fe$fit,col=2) lines(newd[[1]],fe$fit+2*fe$se.fit,col=3) lines(newd[[1]],fe$fit-2*fe$se.fit,col=4) ## prediction with CI... newd <- data.frame(x=seq(-1.2,3.5,length.out=100)) fe <- predict(b,newd,se=TRUE) ylim<- c(min(y,fe$se.fit),max(y,fe$se.fit)) plot(newd[[1]],fe$fit,type="l",ylim=ylim,ylab="y",xlab="x") lines(newd[[1]],fe$fit+2*fe$se.fit,lty=2) lines(newd[[1]],fe$fit-2*fe$se.fit,lty=2) ## prediction on the original data... pr <- predict(b) plot(x,y) lines(x,pr,col=3) ## bivariate example... set.seed(2) n <- 30 x1 <- sort(runif(n)); x2 <- sort(runif(n)*4-1) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- 2*sin(pi*x1[i]) +exp(4*x2[j])/(1+exp(4*x2[j])) f <- as.vector(t(f)); y <- f+rnorm(length(f))*0.1 x11 <- matrix(0,n,n); x11[,1:n] <- x1; x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) b <- scam(y~s(x1,x2,k=c(10,10),bs="tesmi2"),data=dat,optimizer="efs") par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE); plot(b,pers=TRUE,theta = 80, phi = 40) n.out <- 20 xp <- seq(0,1.4,length.out=n.out) zp <- seq(-1,3.4,length.out=n.out) xp1 <- matrix(0,n.out,n.out); xp1[,1:n.out] <- xp xp1 <- as.vector(t(xp1)); xp2 <- rep(zp,n.out) newd <- data.frame(x1=xp1,x2=xp2) fe <- predict(b,newd) fc <- t(matrix(fe,n.out,n.out)) persp(xp,zp,fc,expand= 0.85,ticktype = "simple",xlab="x1", ylab="x2",zlab="f^",main="", theta = 80, phi = 40) ## obtaining a 'prediction matrix'... newd <- data.frame(x1=c(-2,-1),x2=c(0,1)) Xp <- predict(b,newdata=newd,type="lpmatrix") fv <- Xp%*% b$beta.t fv ## End(Not run)
The default print method for a scam object. The code is a clone of print.gam
of the mgcv package with a slight simplification since only two methods of smoothing parameter
selection (by GCV or UBRE) was implemented for scam.
## S3 method for class 'scam' print(x,...)## S3 method for class 'scam' print(x,...)
x |
fitted model objects of class |
... |
other arguments. |
As for mgcv(gam) prints out the family, model formula, effective degrees of freedom for each smooth term,
and optimized value of the smoothness selection criterion used.
Natalya Pya <[email protected]>
Wood S.N. (2006) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press.
Takes a fitted scam object produced by scam() and produces
QQ plots of its residuals (conditional on the fitted model
coefficients and scale parameter). This is an adapted short version of qq.gam() of mgcv package of Simon N Wood.
qq.scam(object, rep=0, level=.9,s.rep=10, type=c("deviance","pearson","response"), pch=".", rl.col=3, rep.col="gray80", ...)qq.scam(object, rep=0, level=.9,s.rep=10, type=c("deviance","pearson","response"), pch=".", rl.col=3, rep.col="gray80", ...)
object |
a fitted |
rep |
How many replicate datasets to generate to simulate quantiles
of the residual distribution. |
level |
If simulation is used for the quantiles, then reference intervals can be provided for the QQ-plot, this specifies the level. 0 or less for no intervals, 1 or more to simply plot the QQ plot for each replicate generated. |
s.rep |
how many times to randomize uniform quantiles to data under direct computation. |
type |
what sort of residuals should be plotted? See
|
pch |
plot character to use. 19 is good. |
rl.col |
color for the reference line on the plot. |
rep.col |
color for reference bands or replicate reference plots. |
... |
extra graphics parameters to pass to plotting functions. |
QQ-plots of the the model residuals can be produced in one of two ways. The cheapest method generates reference quantiles by
associating a quantile of the uniform distribution with each datum, and feeding these uniform quantiles into the quantile function associated with each datum. The resulting quantiles are then used in place of each datum to generate approximate quantiles of residuals.
The residual quantiles are averaged over s.rep randomizations of the uniform quantiles to data.
The second method is to use direct simulatation. For each replicate, data are simulated from the fitted model, and the corresponding residuals computed. This is repeated rep times.
Quantiles are readily obtained from the empirical distribution of residuals so obtained. From this method reference bands are also computable.
Even if rep is set to zero, the routine will attempt to simulate quantiles if no quantile function is available for the family. If no random deviate generating function family is available (e.g. for the quasi families), then a normal QQ-plot is produced. The routine conditions on the fitted model coefficents and the scale parameter estimate.
The plots are very similar to those proposed in Ben and Yohai (2004), but are substantially cheaper to produce (the interpretation of residuals for binary data in Ben and Yohai is not recommended).
Simon N. Wood [email protected]
Natalya Pya [email protected] adapted for usage with scam
N.H. Augustin, E-A Sauleaub, S.N. Wood (2012) On quantile quantile plots for generalized linear models Computational Statistics & Data Analysis. 56(8), 2404-2409.
M.G. Ben and V.J. Yohai (2004) JCGS 13(1), 36-47.
This function is a clone of the mgcv library code residuals.gam.
It returns residuals for a fitted scam model
object. Pearson, deviance, working and response residuals are
available.
## S3 method for class 'scam' residuals(object, type = c("deviance", "pearson","scaled.pearson", "working", "response"),...)## S3 method for class 'scam' residuals(object, type = c("deviance", "pearson","scaled.pearson", "working", "response"),...)
object |
a |
type |
the type of residuals wanted. |
... |
other arguments. |
See residuals.gam for details.
An array of residuals.
Natalya Pya <[email protected]>
This function fits a SCAM to data. Various shape constrained smooths (SCOP-splines), including univariate smooths subject to monotonicity, convexity, or monotonicity plus convexity, bivariate smooths with double or single monotonicity are available as model terms. See shape.constrained.smooth.terms for a complete overview of what is available. Smoothness selection is estimated as part of the fitting. Confidence/credible intervals are available for each smooth term.
The shaped constrained smooths have been added to the gam() in package mgcv setup using the smooth.construct function. The routine calls a gam() function for the model set up, but there are separate functions for the model fitting, scam.fit, and smoothing parameter selection, bfgs_gcv.ubre. Any smooth available in the mgcv can be taken as a model term for SCAM. User-defined smooths can be included as well.
scam(formula, family = gaussian(), data = list(), gamma = 1, sp = NULL, weights = NULL, offset = NULL,optimizer=c("bfgs","newton"), optim.method=c("Nelder-Mead","fd"),scale = 0, knots=NULL, not.exp=FALSE, start= NULL, etastart=NULL,mustart= NULL, control=list(),AR1.rho=0, AR.start=NULL,drop.unused.levels=TRUE)scam(formula, family = gaussian(), data = list(), gamma = 1, sp = NULL, weights = NULL, offset = NULL,optimizer=c("bfgs","newton"), optim.method=c("Nelder-Mead","fd"),scale = 0, knots=NULL, not.exp=FALSE, start= NULL, etastart=NULL,mustart= NULL, control=list(),AR1.rho=0, AR.start=NULL,drop.unused.levels=TRUE)
formula |
A SCAM formula.
This is exactly like the formula for a GAM (see |
family |
A family object specifying the distribution and link to use in
fitting etc. See |
data |
A data frame or list containing the model response variable and
covariates required by the formula. By default the variables are taken
from
|
gamma |
A constant multiplier to inflate the model degrees of freedom in the GCV or UBRE/AIC score. |
sp |
A vector of smoothing parameters can be provided here. Smoothing parameters must be supplied in the order that
the smooth terms appear in the model formula. The default |
weights |
Prior weights on the data. |
offset |
Used to supply a model offset for use in fitting. Note that this offset will always be completely ignored when predicting, unlike an offset
included in |
optimizer |
An array specifying the numerical optimization methods
to optimize the smoothing parameter estimation criterion (specified in the first
element of |
optim.method |
In case of |
scale |
If this is positive then it is taken as the known scale parameter of the exponential family distribution.
Negative value indicates that the scale paraemter is unknown. 0 indicates that the scale parameter is 1 for Poisson and binomial
and unknown otherwise. This conforms to the behaviour of |
knots |
An optional list containing user specified knot values to be used for basis construction. Different terms can use different numbers of knots. |
not.exp |
if |
start |
Initial values for the model coefficients. |
etastart |
Initial values for the linear predictor. |
mustart |
Initial values for the expected values. |
control |
A list of fit control parameters to replace defaults returned by |
AR1.rho |
The AR1 correlation parameter. An AR1 error model can be used for the residuals of Gaussian-identity link models. Standardized residuals (approximately uncorrelated under correct model) returned in
|
AR.start |
logical variable of same length as data, |
drop.unused.levels |
as with |
A shape constrained additive model (SCAM) is a generalized linear model (GLM) in which the linear predictor is given by strictly parametric components plus a sum of smooth functions of the covariates where some of the functions are assumed to be shape constrained. For example,
where the independent response variables follow Poisson distribution with log link function,
, , and are smooth functions of the corresponding covariates, and
is subject to monotone increasing constraint.
Available shape constrained smooths are decsribed in shape.constrained.smooth.terms.
Residual auto-correlation with a simple AR1 correlation structure can be dealt with, for Gaussian models with identity link. Currently, the AR1 correlation parameter should be supplied (rather than estimated) in AR1.rho. AR.start input argument (logical) allows to set independent sections of AR1 correlation. Standardized residuals (approximately uncorrelated under correct model) are returned in std.rsd if AR1.rho is non zero. Use acf(model$std.rsd) for computing and plotting estimates of the autocorrelation function to check correlation.
The function returns an object of class "scam" with the following elements (this agrees with gamObject):
aic |
AIC of the fitted model: the degrees of freedom used to calculate this are the effective degrees of freedom of the model, and the likelihood is evaluated at the maximum of the penalized likelihood, not at the MLE. |
assign |
Array whose elements indicate which model term (listed in
|
bfgs.info |
If |
call |
the matched call. |
coefficients |
the coefficients of the fitted model. Parametric coefficients are first, followed by coefficients for each spline term in turn. |
coefficients.t |
the parametrized coefficients of the fitted model (exponentiated for the monotonic smooths). |
conv |
indicates whether or not the iterative fitting method converged. |
CPU.time |
indicates the real and CPU time (in seconds) taken by the fitting process in case of unknown smoothing parameters |
data |
the original supplied data argument. Only included if the |
deviance |
model deviance (not penalized deviance). |
df.null |
null degrees of freedom. |
df.residual |
effective residual degrees of freedom of the model. |
edf |
estimated degrees of freedom for each model parameter. Penalization means that many of these are less than 1. |
edf1 |
alternative estimate of edf. |
efs.info |
If |
family |
family object specifying distribution and link used. |
fitted.values |
fitted model predictions of expected value for each datum. |
formula |
the model formula. |
gcv.ubre |
the minimized GCV or UBRE score. |
dgcv.ubre |
the gradient of the GCV or UBRE score. |
iter |
number of iterations of the Newton-Raphson method taken to get convergence. |
linear.predictors |
fitted model prediction of link function of expected value for each datum. |
method |
|
min.edf |
Minimum possible degrees of freedom for whole model. |
model |
model frame containing all variables needed in original model fit. |
nlm.info |
If |
not.exp |
if |
nsdf |
number of parametric, non-smooth, model terms including the intercept. |
null.deviance |
deviance for single parameter model. |
offset |
model offset. |
optim.info |
If |
prior.weights |
prior weights on observations. |
pterms |
|
R |
Factor R from QR decomposition of weighted model matrix, unpivoted to be in same column order as model matrix. |
residuals |
the working residuals for the fitted model. |
scale.estimated |
|
sig2 |
estimated or supplied variance/scale parameter. |
smooth |
list of smooth objects, containing the basis information for each term in the
model formula in the order in which they appear. These smooth objects are returned by
the |
sp |
estimated smoothing parameters for the model. These are the underlying smoothing parameters, subject to optimization. |
std.rsd |
Standardized residuals (approximately uncorrelated under correct model) if |
termcode |
an integer indicating why the optimization process of the smoothness selection
terminated (see |
terms |
|
trA |
trace of the influence matrix, total number of the estimated degrees of freedom ( |
var.summary |
A named list of summary information on the predictor variables. See |
Ve |
frequentist estimated covariance matrix for the parameter estimators. |
Vp |
estimated covariance matrix for the parameters. This is a Bayesian posterior covariance matrix that results from adopting a particular Bayesian model of the smoothing process. |
Ve.t |
frequentist estimated covariance matrix for the reparametrized parameter estimators obtained using the delta method. Particularly useful for testing whether terms are zero. Not so useful for CI's as smooths are usually biased. |
Vp.t |
estimated covariance matrix for the reparametrized parameters obtained using the delta method. Paricularly useful for creating credible/confidence intervals. |
weights |
final weights used in the Newton-Raphson iteration. |
cmX |
column means of the model matrix (with elements corresponding to smooths set to zero). |
contrasts |
contrasts associated with a factor. |
xlevels |
levels of a factor variable used in the model. |
y |
response data. |
Natalya Pya <[email protected]> based partly on mgcv by Simon Wood
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1):3-36
Wood S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapman and Hall/CRC Press
Wood, S.N. (2006) On confidence intervals for generalized additive models based on penalized regression splines. Australian and New Zealand Journal of Statistics. 48(4): 445-464.
Wood, S.N. and M. Fasiolo (2017) A generalized Fellner-Schall method for smoothing parameter optimization with application to Tweedie location, scale and shape models. Biometrics 73 (4), 1071-1081
scam-package, shape.constrained.smooth.terms,
gam, s,
plot.scam, summary.scam, scam.check, predict.scam
##============================= ## Gaussian model, two smooth terms: unconstrained and increasing... ## simulating data... require(scam) set.seed(4) n <- 200 x1 <- runif(n)*6-3 f1 <- 3*exp(-x1^2) # unconstrained term x2 <- runif(n)*4-1; f2 <- exp(4*x2)/(1+exp(4*x2)) # monotone increasing smooth y <- f1+f2 +rnorm(n)*.5 dat <- data.frame(x1=x1,x2=x2,y=y) ## fit model, get results, and plot... b <- scam(y~s(x1,bs="cr")+s(x2,bs="mpi"),data=dat) b summary(b) plot(b,pages=1,shade=TRUE) ##=============================== ## Gaussian model, two smooth terms: increasing and mixed (decreasing and convex)... ## simulating data... set.seed(4) n <- 200 x1 <- runif(n)*4-1; f1 <- exp(4*x1)/(1+exp(4*x1)) # increasing smooth x2 <- runif(n)*3-1; f2 <- exp(-3*x2)/15 # decreasing and convex smooth y <- f1+f2 + rnorm(n)*.4 dat <- data.frame(x1=x1,x2=x2,y=y) ## fit model, results, and plot... b <- scam(y~ s(x1,bs="mpi")+s(x2, bs="mdcx"),data=dat) summary(b) plot(b,pages=1,scale=0,shade=TRUE) ##================================= ## Not run: ## using the extended Fellner-Schall method for smoothing parameter selection... b0 <- scam(y~ s(x1,bs="mpi")+s(x2,bs="mdcx"),data=dat,optimizer="efs") summary(b0) ## using the extended Fellner-Schall method for smoothing parameter selection, ## and BFGS for model coefficient estimation... b0 <- scam(y~ s(x1,bs="mpi")+s(x2,bs="mdcx"),data=dat,optimizer=c("efs","bfgs")) summary(b0) ## using optim() for smoothing parameter selection... b1 <- scam(y~ s(x1,bs="mpi")+s(x2,bs="mdcx"),data=dat,optimizer="optim") summary(b1) b2 <- scam(y~ s(x1,bs="mpi")+s(x2,bs="mdcx"),data=dat,optimizer="optim", optim.method=c("BFGS","fd")) summary(b2) ## using nlm()... b3 <- scam(y~ s(x1,bs="mpi")+s(x2,bs="mdcx"),data=dat,optimizer="nlm") summary(b3) ## End(Not run) ##=================================== ## Poisson model .... ## simulating data... set.seed(2) n <- 200 x1 <- runif(n)*6-3 f1 <- 3*exp(-x1^2) # unconstrained term x2 <- runif(n)*4-1; f2 <- exp(4*x2)/(1+exp(4*x2)) # monotone increasing smooth f <- f1+f2 y <- rpois(n,exp(f)) dat <- data.frame(x1=x1,x2=x2,y=y) ## fit model, get results, and plot... b <- scam(y~s(x1,bs="cr")+s(x2,bs="mpi"), family=poisson(link="log"),data=dat,optimizer=c("efs","bfgs")) summary(b) plot(b,pages=1,shade=TRUE) scam.check(b) ## Gamma model... ## simulating data... set.seed(6) n <- 300 x1 <- runif(n)*6-3 f1 <- 1.5*sin(1.5*x1) # unconstrained term x2 <- runif(n)*4-1; f2 <- 1.5/(1+exp(-10*(x2+.75)))+1.5/(1+exp(-5*(x2-.75))) # increasing smooth x3 <- runif(n)*6-3; f3 <- 3*exp(-x3^2) # unconstrained term f <- f1+f2+f3 y <- rgamma(n,shape=1,scale=exp(f)) dat <- data.frame(x1=x1,x2=x2,x3=x3,y=y) ## fit model, get results, and plot... b <- scam(y~s(x1,bs="ps")+s(x2,k=15,bs="mpi")+s(x3,bs="ps"), family=Gamma(link="log"),data=dat,optimizer=c("efs","bfgs")) b summary(b) par(mfrow=c(2,2)) plot(b,shade=TRUE) ## Not run: ## bivariate example... ## simulating data... set.seed(2) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(runif(n)) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) { f1[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i]))+2*sin(pi*x2[j])} f <- as.vector(t(f1)) y <- f+rnorm(length(f))*.2 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model and plot ... b <- scam(y~s(x1,x2,k=c(10,10),bs=c("tesmd1","ps")),data=dat,optimizer="efs") summary(b) old.par <- par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data",pch=".",cex=3) par(old.par) ## example with random effect smoother... set.seed(2) n <- 200 x1 <- runif(n)*6-3 f1 <- 3*exp(-x1^2) # unconstrained term x2 <- runif(n)*4-1; f2 <- exp(4*x2)/(1+exp(4*x2)) # increasing smooth f <- f1+f2 a <- factor(sample(1:10,200,replace=TRUE)) Xa <- model.matrix(~a-1) # random main effects y <- f + Xa%*%rnorm(length(levels(a)))*.5 + rnorm(n)*.4 dat <- data.frame(x1=x1,x2=x2,y=y,a=a) ## fit model and plot... b <- scam(y~s(x1,bs="cr")+s(x2,bs="mpi")+s(a,bs="re"), data=dat) summary(b) scam.check(b) plot(b,pages=1,shade=TRUE) ## example with AR1 errors... set.seed(8) n <- 500 x1 <- runif(n)*6-3 f1 <- 3*exp(-x1^2) # unconstrained term x2 <- runif(n)*4-1; f2 <- exp(4*x2)/(1+exp(4*x2)) # increasing smooth f <- f1+f2 e <- rnorm(n,0,sd=2) for (i in 2:n) e[i] <- .6*e[i-1] + e[i] y <- f + e dat <- data.frame(x1=x1,x2=x2,y=y) b <- scam(y~s(x1,bs="cr")+s(x2,k=25,bs="mpi"), data=dat, AR1.rho=.6, optimizer="efs") b ## Raw residuals still show correlation... acf(residuals(b)) ## But standardized are now fine... x11() acf(b$std.rsd) ## End(Not run)##============================= ## Gaussian model, two smooth terms: unconstrained and increasing... ## simulating data... require(scam) set.seed(4) n <- 200 x1 <- runif(n)*6-3 f1 <- 3*exp(-x1^2) # unconstrained term x2 <- runif(n)*4-1; f2 <- exp(4*x2)/(1+exp(4*x2)) # monotone increasing smooth y <- f1+f2 +rnorm(n)*.5 dat <- data.frame(x1=x1,x2=x2,y=y) ## fit model, get results, and plot... b <- scam(y~s(x1,bs="cr")+s(x2,bs="mpi"),data=dat) b summary(b) plot(b,pages=1,shade=TRUE) ##=============================== ## Gaussian model, two smooth terms: increasing and mixed (decreasing and convex)... ## simulating data... set.seed(4) n <- 200 x1 <- runif(n)*4-1; f1 <- exp(4*x1)/(1+exp(4*x1)) # increasing smooth x2 <- runif(n)*3-1; f2 <- exp(-3*x2)/15 # decreasing and convex smooth y <- f1+f2 + rnorm(n)*.4 dat <- data.frame(x1=x1,x2=x2,y=y) ## fit model, results, and plot... b <- scam(y~ s(x1,bs="mpi")+s(x2, bs="mdcx"),data=dat) summary(b) plot(b,pages=1,scale=0,shade=TRUE) ##================================= ## Not run: ## using the extended Fellner-Schall method for smoothing parameter selection... b0 <- scam(y~ s(x1,bs="mpi")+s(x2,bs="mdcx"),data=dat,optimizer="efs") summary(b0) ## using the extended Fellner-Schall method for smoothing parameter selection, ## and BFGS for model coefficient estimation... b0 <- scam(y~ s(x1,bs="mpi")+s(x2,bs="mdcx"),data=dat,optimizer=c("efs","bfgs")) summary(b0) ## using optim() for smoothing parameter selection... b1 <- scam(y~ s(x1,bs="mpi")+s(x2,bs="mdcx"),data=dat,optimizer="optim") summary(b1) b2 <- scam(y~ s(x1,bs="mpi")+s(x2,bs="mdcx"),data=dat,optimizer="optim", optim.method=c("BFGS","fd")) summary(b2) ## using nlm()... b3 <- scam(y~ s(x1,bs="mpi")+s(x2,bs="mdcx"),data=dat,optimizer="nlm") summary(b3) ## End(Not run) ##=================================== ## Poisson model .... ## simulating data... set.seed(2) n <- 200 x1 <- runif(n)*6-3 f1 <- 3*exp(-x1^2) # unconstrained term x2 <- runif(n)*4-1; f2 <- exp(4*x2)/(1+exp(4*x2)) # monotone increasing smooth f <- f1+f2 y <- rpois(n,exp(f)) dat <- data.frame(x1=x1,x2=x2,y=y) ## fit model, get results, and plot... b <- scam(y~s(x1,bs="cr")+s(x2,bs="mpi"), family=poisson(link="log"),data=dat,optimizer=c("efs","bfgs")) summary(b) plot(b,pages=1,shade=TRUE) scam.check(b) ## Gamma model... ## simulating data... set.seed(6) n <- 300 x1 <- runif(n)*6-3 f1 <- 1.5*sin(1.5*x1) # unconstrained term x2 <- runif(n)*4-1; f2 <- 1.5/(1+exp(-10*(x2+.75)))+1.5/(1+exp(-5*(x2-.75))) # increasing smooth x3 <- runif(n)*6-3; f3 <- 3*exp(-x3^2) # unconstrained term f <- f1+f2+f3 y <- rgamma(n,shape=1,scale=exp(f)) dat <- data.frame(x1=x1,x2=x2,x3=x3,y=y) ## fit model, get results, and plot... b <- scam(y~s(x1,bs="ps")+s(x2,k=15,bs="mpi")+s(x3,bs="ps"), family=Gamma(link="log"),data=dat,optimizer=c("efs","bfgs")) b summary(b) par(mfrow=c(2,2)) plot(b,shade=TRUE) ## Not run: ## bivariate example... ## simulating data... set.seed(2) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(runif(n)) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) { f1[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i]))+2*sin(pi*x2[j])} f <- as.vector(t(f1)) y <- f+rnorm(length(f))*.2 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model and plot ... b <- scam(y~s(x1,x2,k=c(10,10),bs=c("tesmd1","ps")),data=dat,optimizer="efs") summary(b) old.par <- par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data",pch=".",cex=3) par(old.par) ## example with random effect smoother... set.seed(2) n <- 200 x1 <- runif(n)*6-3 f1 <- 3*exp(-x1^2) # unconstrained term x2 <- runif(n)*4-1; f2 <- exp(4*x2)/(1+exp(4*x2)) # increasing smooth f <- f1+f2 a <- factor(sample(1:10,200,replace=TRUE)) Xa <- model.matrix(~a-1) # random main effects y <- f + Xa%*%rnorm(length(levels(a)))*.5 + rnorm(n)*.4 dat <- data.frame(x1=x1,x2=x2,y=y,a=a) ## fit model and plot... b <- scam(y~s(x1,bs="cr")+s(x2,bs="mpi")+s(a,bs="re"), data=dat) summary(b) scam.check(b) plot(b,pages=1,shade=TRUE) ## example with AR1 errors... set.seed(8) n <- 500 x1 <- runif(n)*6-3 f1 <- 3*exp(-x1^2) # unconstrained term x2 <- runif(n)*4-1; f2 <- exp(4*x2)/(1+exp(4*x2)) # increasing smooth f <- f1+f2 e <- rnorm(n,0,sd=2) for (i in 2:n) e[i] <- .6*e[i-1] + e[i] y <- f + e dat <- data.frame(x1=x1,x2=x2,y=y) b <- scam(y~s(x1,bs="cr")+s(x2,k=25,bs="mpi"), data=dat, AR1.rho=.6, optimizer="efs") b ## Raw residuals still show correlation... acf(residuals(b)) ## But standardized are now fine... x11() acf(b$std.rsd) ## End(Not run)
Takes a fitted scam object produced by scam() and produces some diagnostic information
about the fitting procedure and results. This function is almost the same as gam.check of the mgcv library. The default is to produce four residual plots and some information about the
convergence of the smoothness selection optimization.
scam.check(b,type=c("deviance","pearson","response"),old.style=FALSE, pch=".", rep=0, level=.9, rl.col=3, rep.col="gray80",...)scam.check(b,type=c("deviance","pearson","response"),old.style=FALSE, pch=".", rep=0, level=.9, rl.col=3, rep.col="gray80",...)
b |
a fitted |
old.style |
produces qq-norm plots as it was in scam versions < 1.2-15 when set to |
type |
type of residuals, see |
rep, level, rep.col
|
arguments passed to |
rl.col |
color for the reference line on the quantile-quantile plot. |
pch |
plot character to use for the quantile-quantile plot. |
... |
extra graphics parameters to pass to plotting functions. |
As for mgcv(gam) plots 4 standard diagnostic plots, and some other
convergence diagnostics. The
printed information relates to the optimization process used to select smoothing
parameters.
Natalya Pya [email protected] based partly on mgcv by Simon N Wood
Wood S.N. (2006) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press.
## Not run: library(scam) set.seed(2) n <- 200 x1 <- runif(n)*4-1; f1 <- exp(4*x1)/(1+exp(4*x1)) # monotone increasing smooth x2 <- runif(n)*3-1; f2 <- exp(-3*x2)/15 # monotone decreasing and convex smooth f <- f1+f2 y <- f+ rnorm(n)*0.2 dat <- data.frame(x1=x1,x2=x2,y=y) b <- scam(y~ s(x1,k=25,bs="mpi",m=2)+s(x2,k=25,bs="mdcx",m=2), family=gaussian(link="identity"),data=dat) plot(b,pages=1) scam.check(b) ## End(Not run)## Not run: library(scam) set.seed(2) n <- 200 x1 <- runif(n)*4-1; f1 <- exp(4*x1)/(1+exp(4*x1)) # monotone increasing smooth x2 <- runif(n)*3-1; f2 <- exp(-3*x2)/15 # monotone decreasing and convex smooth f <- f1+f2 y <- f+ rnorm(n)*0.2 dat <- data.frame(x1=x1,x2=x2,y=y) b <- scam(y~ s(x1,k=25,bs="mpi",m=2)+s(x2,k=25,bs="mdcx",m=2), family=gaussian(link="identity"),data=dat) plot(b,pages=1) scam.check(b) ## End(Not run)
This is an internal function of package scam which allows
control of the numerical options for fitting a SCAM.
scam.control(maxit = 200, maxHalf=30, devtol.fit=1e-7, steptol.fit=1e-7, keepData=FALSE,efs.lspmax=15,efs.tol=.1, nlm=list(),optim=list(), bfgs=list(), trace =FALSE, print.warn=FALSE,b.notexp=1, threshold.notexp=20)scam.control(maxit = 200, maxHalf=30, devtol.fit=1e-7, steptol.fit=1e-7, keepData=FALSE,efs.lspmax=15,efs.tol=.1, nlm=list(),optim=list(), bfgs=list(), trace =FALSE, print.warn=FALSE,b.notexp=1, threshold.notexp=20)
maxit |
Maximum number of IRLS iterations to perform used in |
maxHalf |
If a step of the BFGS optimization method leads
to a worse penalized deviance, then the step length of the model coefficients is halved. This is
the number of halvings to try before giving up used in |
devtol.fit |
A positive scalar giving the convergence control for the model fitting algorithm in |
steptol.fit |
A positive scalar giving the tolerance at which the scaled distance between
two successive iterates is considered close enough to zero to terminate the model fitting algorithm in |
keepData |
Should a copy of the original |
efs.lspmax |
maximum log smoothing parameters to allow under extended Fellner Schall smoothing parameter optimization. |
efs.tol |
change in GCV to count as negligible when testing for EFS convergence. If the step is small and the last 3 steps led to a GCV change smaller than this, then stop. |
nlm |
list of control parameters to pass to |
optim |
list of control parameters to pass to |
bfgs |
list of control parameters to pass to default BFGS optimizer used for outer estimation of log smoothing parameters. |
trace |
turns on or off some de-bugging information. |
print.warn |
when set to |
b.notexp |
parameter |
threshold.notexp |
parameter |
Outer iteration is used to estimate smoothing parameters of SCAM by GCV/UBRE score optimization. The default procedure is the built-in BFGS method which is controlled by the list bfgs with the following elements: steptol.bfgs (default 1e-7) is the relative convergence tolerance;
gradtol.bfgs (default 6.0554*1e-6) is a tolerance at which the gradient is considered to be close enougth to 0 to terminate the BFGS algorithm;
maxNstep is a positive scalar which gives the maximum allowable step length (default 5);
maxHalf gives the maximum number of step halving in "backtracking" to permit before giving up(default 30);
check.analytical is logical whether the analytical gradient of GCV/UBRE should be checked numerically (default FALSE);
del is an increment for finite differences when checking analytical gradients (default 1e-4).
If outer iteration using nlm is used for fitting, then the control list nlm stores control arguments
for calls to routine nlm. As in gam.control the list has the following named elements: ndigit is the number
of significant digits in the GCV/UBRE score; gradtol is the tolerance used to judge convergence of the gradient of the GCV/UBRE score to zero
(default 1e-6); stepmax is the maximum allowable log smoothing parameter step (default 2);
steptol is the minimum allowable step length (default 1e-4);
iterlim is the maximum number of optimization steps allowed (default 200);
check.analyticals indicates whether the built in exact derivative calculations should be checked numerically (default FALSE).
Any of these which are not supplied and named in the list are set to their default values.
Outer iteration using optim is controlled using list optim, which currently has one element: factr
which takes default value 1e7.
Natalya Pya Arnqvist [email protected] based partly on gam.control by Simon Wood
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1):3-36
This routine estimates SCAM coefficients given log smoothing parameters using the Newton-Raphson method.
The estimation of the smoothing parameters by the GCV/UBRE score optimization is outer to the model fitting. Routine
gcv.ubre_grad evaluates the first derivatives of the smoothness selection scores with respect to the
log smoothing parameters. Routine bfgs_gcv.ubre estimates the smoothing parameters using the BFGS method.
The function is not normally called directly, but rather service routines for scam.
scam.fit(G,sp, etastart=NULL, mustart=NULL, env=env, null.coef=rep(0,ncol(G$X)), control=scam.control())scam.fit(G,sp, etastart=NULL, mustart=NULL, env=env, null.coef=rep(0,ncol(G$X)), control=scam.control())
G |
A list of items needed to fit a SCAM. |
sp |
The vector of smoothing parameters. |
etastart |
Initial values for the linear predictor. |
mustart |
Initial values for the expected values. |
env |
Get the enviroment for the model coefficients, their derivatives and the smoothing parameter. |
null.coef |
coefficients for a null model, needed for an ability to check for immediate divergence. |
control |
A list of fit control parameters returned by |
The routine applies step halving to any step that increases the penalized deviance substantially.
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
Wood, S.N. (2008) Fast stable direct fitting and smoothness selection for generalized additive models. Journal of the Royal Statistical Society (B) 70(3):495-518
Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1):3-36
As in mgcv(gam), shape preserving smooth terms are specified in a scam formula using s
terms. All the shape constrained smooth terms (SCOP-splines) are constructed using the B-splines basis proposed by Eilers and Marx (1996) with a discrete penalty on the basis coefficients.
The univariate single penalty built-in shape constrained smooth classes are summarized as follows.
Monotone increasing SCOP-splines: bs="mpi". To achieve monotone increasing smooths this reparameterizes the coefficients so that they form an increasing sequence.
For details see smooth.construct.mpi.smooth.spec.
Monotone decreasing SCOP-splines: bs="mpd". To achieve monotone decreasing smooths this reparameterizes the coefficients
so that they form a decreasing sequence. A first order difference penalty applied to the basis coefficients starting with the
second is used for the monotone increasing and decreasing cases.
Convex SCOP-splines: bs="cx". This reparameterizes the coefficients so that the second order differences of the basis coefficients are greater than zero.
For details see smooth.construct.cx.smooth.spec.
Concave SCOP-splines: bs="cv". This reparameterizes the coefficients so that the second order differences of the basis coefficients are less than zero.
For details see smooth.construct.cv.smooth.spec.
Increasing and convex SCOP-splines: bs="micx". This reparameterizes the coefficients
so that the first and the second order differences of the basis coefficients are greater than zero.
For details see
smooth.construct.micx.smooth.spec.
Increasing and concave SCOP-splines: bs="micv". This reparameterizes the coefficients
so that the first order differences of the basis coefficients are greater than zero while the second order difference
are less than zero.
Decreasing and convex SCOP-splines: bs="mdcx". This reparameterizes the coefficients
so that the first order differences of the basis coefficients are less than zero while the second order difference
are greater.
For details see smooth.construct.mdcx.smooth.spec.
Decreasing and concave SCOP-splines: bs="mdcv". This reparameterizes the coefficients
so that the first and the second order differences of the basis coefficients are less than zero.
Increasing with an additional 'finish-at-zero' constraint SCOP-splines: bs="mifo". This sets the last (m+1) spline coefficients to zero. According to the B-spline basis functions properties, the value of the spline, f(x), is determined by m+2 non-zero basis functions, and only m+1 B-splines are non-zero at knots. Only m+2 B-splines are non-zero on any (k_i, k_{i+1}), and the sum of these m+2 basis functions is 1.
For details see smooth.construct.mifo.smooth.spec.
Increasing with an additional 'start-at-zero' constraint SCOP-spline: bs="miso". This sets the first (m+1) spline coefficients to zero. According to the B-spline basis functions properties, the value of the spline, f(x), is determined by m+2 non-zero basis functions, and only m+1 B-splines are non-zero at knots. Only m+2 B-splines are non-zero on any (k_i, k_{i+1}), and the sum of these m+2 basis functions is 1.
For details see smooth.construct.miso.smooth.spec.
SCOP-spline with positivity constraint: bs="po". This reparameterizes the coefficients so that there are positive. For details see smooth.construct.po.smooth.spec.
Decreasing/increasing SCOP-splines used with numeric 'by' variable: bs="mpdBy", bs="mpiBy". These work similar to mpd.smooth.spec, mpi.smooth.spec, but without applying an identifiability constraint ('zero intercept' constraint). Use when the smooth term has a numeric by variable that takes more than one value.
For details see smooth.construct.mpd.smooth.spec, smooth.construct.mpi.smooth.spec.
Convex/concave SCOP-splines used with numeric 'by' variable: bs="cxdBy", bs="cvBy". These work similar to cx.smooth.spec, cv.smooth.spec, but without applying an identifiability constraint ('zero intercept' constraint). Use when the smooth term has a numeric by variable that takes more than one value.
For details see smooth.construct.cx.smooth.spec,
Decreasing/increasing and convex/concave SCOP-splines used with numeric 'by' variable: bs="mdcxBy", bs="mdcvBy", bs="micxBy", bs="micvBy".
These work similar to mdcx.smooth.spec, mdcv.smooth.spec, micx.smooth.spec, micv.smooth.spec, but without applying an identifiability constraint ('zero intercept' constraint). Use when the smooth term has a numeric by variable that takes more than one value.
For details see smooth.construct.mdcx.smooth.spec,
smooth.construct.mdcv.smooth.spec,
Locally shape-constrained P-splines (LSCOP-splines): bs="lmpi". Locally increasing splines that are monotone increasing up to a specified change point and become unconstrained beyond that point. For details see smooth.construct.lmpi.smooth.spec.
For all types of the mixed constrained smoothing a first order difference penalty applied to the basis coefficients starting with the third one is used. Centring ('sum-to-zero') constraint has been applied to univariate SCOP-splines subject to monotonicity (convexity) constraints after implementing the 'zero intercept' identifiability constraint. This is achieved by dropping the first (constant) column of the spline model matrix and subtracting the corresponding column means from the elements of the remaining columns afterwards. 'Sum-to-zero' constraint orthogonalized the smooth to the model intercept term, thus avoiding confounding with the intercept. The standard errors of the estimated intercept become lower with the centring constraint.
Using the concept of the tensor product spline bases bivariate smooths under monotonicity constraint where monotonicity may be assumed on only one of the covariates (single monotonicity) or both of them (double monotonicity) are added as the smooth terms of the SCAM. Bivariate B-spline is constructed by expressing the coefficients of one of the marginal univariate B-spline bases as the B-spline of the other covariate. Double or single monotonicity is achieved by the corresponding re-parametrization of the bivariate basis coefficients to satisfy the sufficient conditions formulated in terms of the first order differences of the coefficients. The following explains the built in bivariate shape constrained smooth classes.
Double monotone increasing SCOP-splines: bs="tedmi".
See smooth.construct.tedmi.smooth.spec for details.
Double monotone decreasing SCOP-splines: bs="tedmd".
Single monotone increasing SCOP-splines along the first covariate direction: bs="tesmi1".
Single monotone increasing SCOP-splines along the second covariate direction: bs="tesmi2".
Single monotone decreasing SCOP-splines along the first covariate direction: bs="tesmd1".
Single monotone decreasing SCOP-splines along the second covariate direction: bs="tesmd2".
SCOP-splines with double concavity constraint: bs="tecvcv".
See smooth.construct.tecvcv.smooth.spec for details.
SCOP-splines with double convexity constraint: bs="tecxcx".
See smooth.construct.tecxcx.smooth.spec for details.
SCOP-splines with convexity wrt the first covariate and
concavity wrt the second covariate: bs="tecxcv". See smooth.construct.tecxcv.smooth.spec for details.
Decreasing along the first covariate and
concave along the second covariate SCOP-splines: bs="tedecv".
See smooth.construct.tedecv.smooth.spec for details.
Decreasing along the first covariate and convex along the second covariate SCOP-splines: bs="tedecx".
See smooth.construct.tedecx.smooth.spec for details.
Increasing along the first covariate and concave along the second covariate SCOP-splines: bs="temicv".
See smooth.construct.temicv.smooth.spec for details.
Increasing along the first covariate and convex along the second covariate SCOP-splines: bs="temicx".
See smooth.construct.temicx.smooth.spec for details.
Convex along the second covariate SCOP-splines: bs="tescx".
See smooth.construct.tescx.smooth.spec for details.
Concave along the second covariate SCOP-splines: bs="tescv".
See smooth.construct.tescv.smooth.spec for details.
Tensor product interaction with increasing constraint along the first
covariate and unconstrained along the second covariate: bs="tismi".
See smooth.construct.tismi.smooth.spec for details.
Tensor product interaction with decreasing constraint along the first
covariate and unconstrained along the second covariate: bs="tismd".
See smooth.construct.tismd.smooth.spec for details.
Double penalties for the shape constrained tensor product smooths are obtained from the penalties of the marginal smooths. For the bivariate SCOP-splines with monotonicity (convexity) constraints along one covariate, the 'sum-to-zero' constraints are applied after dropping the first columns of the model matrix of the constrained marginal smooth. The basis for the unconstrained marginal must be non-negative over the region where the marginal monotonicity (convexity) is to hold. For the bivariate interaction smooths "tismi" and "tismd" the following identifiability steps are implemented: i) dropped the first column of the "mpi" ("mpd") marginals, ii) applied 'sum-to-zero' constraints to the marginals and to the unconstrained B-spline basis, iii) tensor product constructed. The 'sum-to-zero' constraint is applied to the final tensor product model matrix afters removing its first column when constructing bivariate SCOP-splines with double monotonicity (convexity). These result in faster convergence of the optimization routines and more stable intercept estimates.
Also linear functionals of smooths with shape constraints (increasing/decreasing and convex/concave) are
supported. See linear.functional.terms.
Natalya Pya [email protected]
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
Eilers, P.H.C. and B.D. Marx (1996) Flexible Smoothing with B-splines and Penalties. Statistical Science, 11(2):89-121
Wood S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapman and Hall/CRC Press
Wood, S.N. (2006) Low rank scale invariant tensor product smooths for generalized additive mixed models. Biometrics 62(4):1025-1036
s, smooth.construct.mpi.smooth.spec,
smooth.construct.mpd.smooth.spec,
smooth.construct.cx.smooth.spec,
smooth.construct.cv.smooth.spec,
smooth.construct.micx.smooth.spec,
smooth.construct.micv.smooth.spec,
smooth.construct.mdcx.smooth.spec,
smooth.construct.mdcv.smooth.spec,
smooth.construct.tedmi.smooth.spec,
smooth.construct.tedmd.smooth.spec,
smooth.construct.tesmi1.smooth.spec,
smooth.construct.tesmi2.smooth.spec,
smooth.construct.tesmd1.smooth.spec,
smooth.construct.tesmd2.smooth.spec,
smooth.construct.tismi.smooth.spec,
smooth.construct.tismd.smooth.spec
## see examples for scam## see examples for scam
This is a special method function
for creating smooths subject to concavity constraint which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed using concave P-splines. This smooth is specified via model terms such as
s(x,k,bs="cv",m=2),
where k denotes the basis dimension and m+1 is the order of the B-spline basis.
cvBy.smooth.spec works similar to cv.smooth.spec but without applying an identifiability constraint ('zero intercept' constraint). cvBy.smooth.spec should be used when the smooth term has a numeric by variable that takes more than one value. In such cases, the smooth terms are fully identifiable without a 'zero intercept' constraint, so they are left unconstrained. This smooth is specified as
s(x,by=z,bs="cvBy"). See an example below.
However a factor by variable requires identifiability constraints, so s(x,by=fac,bs="cv") is used in this case.
## S3 method for class 'cv.smooth.spec' smooth.construct(object, data, knots) ## S3 method for class 'cvBy.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'cv.smooth.spec' smooth.construct(object, data, knots) ## S3 method for class 'cvBy.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the data required by this term,
with names given by |
knots |
An optional list containing the knots supplied for basis setup.
If it is |
An object of class "cv.smooth", "cvBy.smooth".
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
smooth.construct.cx.smooth.spec,
smooth.construct.mpi.smooth.spec,
smooth.construct.mdcv.smooth.spec,
smooth.construct.mdcx.smooth.spec,
smooth.construct.micx.smooth.spec,
smooth.construct.mpd.smooth.spec
## Not run: ## Concave P-splines example ## simulating data... require(scam) set.seed(1) n <- 100 x <- sort(2*runif(n)-1) f <- -4*x^2 y <- f + rnorm(n)*0.45 dat <- data.frame(x=x,y=y) b <- scam(y~s(x,k=15,bs="cv"),family=gaussian,data=dat,not.exp=FALSE) ## fit unconstrained model... b1 <- scam(y~s(x,k=15,bs="cr"),family=gaussian, data=dat,not.exp=FALSE) ## plot results ... plot(x,y,xlab="x",ylab="y",cex=.5) lines(x,f) ## the true function lines(x,b$fitted,col=2) ## constrained fit lines(x,b1$fitted,col=3) ## unconstrained fit ## Poisson version... y <- rpois(n,15*exp(f)) dat <- data.frame(x=x,y=y) ## fit model ... b <- scam(y~s(x,k=15,bs="cv"),family=poisson(link="log"),data=dat,not.exp=FALSE) ## fit unconstrained model... b1<-scam(y~s(x,k=15,bs="cr"),family=poisson(link="log"), data=dat,not.exp=FALSE) ## plot results ... plot(x,y,xlab="x",ylab="y",cex=.5) lines(x,15*exp(f)) ## the true function lines(x,b$fitted,col=2) ## constrained fit lines(x,b1$fitted,col=3) ## unconstrained fit ## plotting on log scale... plot(x,log(15*exp(f)),type="l",cex=.5) ## the true function lines(x,log(b$fitted),col=2) ## constrained fit lines(x,log(b1$fitted),col=3) ## unconstrained fit ## 'by' factor example... set.seed(9) n <- 400 x <- sort(runif(n,-.5,.5)) f1 <- -.7*x+cos(x)-3 f2 <- -20*x^2 par(mfrow=c(1,2)) plot(x,f1,type="l");plot(x,f2,type="l") e <- rnorm(n, 0, 1.5) fac <- as.factor(sample(1:2,n,replace=TRUE)) fac.1 <- as.numeric(fac==1) fac.2 <- as.numeric(fac==2) y <- f1*fac.1 + f2*fac.2 + e dat <- data.frame(y=y,x=x,fac=fac,f1=f1,f2=f2) b2 <- scam(y ~ fac+s(x,by=fac,bs="cv"),data=dat,optimizer="efs") plot(b2,pages=1,scale=0,shade=TRUE) summary(b2) x11() vis.scam(b2,theta=50,color="terrain") ## numeric 'by' variable example... set.seed(6) n <- 100 x <- sort(2*runif(n)-1) z <- runif(n,-2,3) f <- -4*x^2 y <- f*z + rnorm(n)*0.6 dat <- data.frame(x=x,z=z,y=y) b <- scam(y~s(x,k=15,by=z,bs="cvBy"),data=dat) summary(b) par(mfrow=c(1,2)) plot(b,shade=TRUE) ## unconstrained fit... b1 <- scam(y~s(x,k=15,by=z),data=dat) plot(b1,shade=TRUE) summary(b1) ## End(Not run)## Not run: ## Concave P-splines example ## simulating data... require(scam) set.seed(1) n <- 100 x <- sort(2*runif(n)-1) f <- -4*x^2 y <- f + rnorm(n)*0.45 dat <- data.frame(x=x,y=y) b <- scam(y~s(x,k=15,bs="cv"),family=gaussian,data=dat,not.exp=FALSE) ## fit unconstrained model... b1 <- scam(y~s(x,k=15,bs="cr"),family=gaussian, data=dat,not.exp=FALSE) ## plot results ... plot(x,y,xlab="x",ylab="y",cex=.5) lines(x,f) ## the true function lines(x,b$fitted,col=2) ## constrained fit lines(x,b1$fitted,col=3) ## unconstrained fit ## Poisson version... y <- rpois(n,15*exp(f)) dat <- data.frame(x=x,y=y) ## fit model ... b <- scam(y~s(x,k=15,bs="cv"),family=poisson(link="log"),data=dat,not.exp=FALSE) ## fit unconstrained model... b1<-scam(y~s(x,k=15,bs="cr"),family=poisson(link="log"), data=dat,not.exp=FALSE) ## plot results ... plot(x,y,xlab="x",ylab="y",cex=.5) lines(x,15*exp(f)) ## the true function lines(x,b$fitted,col=2) ## constrained fit lines(x,b1$fitted,col=3) ## unconstrained fit ## plotting on log scale... plot(x,log(15*exp(f)),type="l",cex=.5) ## the true function lines(x,log(b$fitted),col=2) ## constrained fit lines(x,log(b1$fitted),col=3) ## unconstrained fit ## 'by' factor example... set.seed(9) n <- 400 x <- sort(runif(n,-.5,.5)) f1 <- -.7*x+cos(x)-3 f2 <- -20*x^2 par(mfrow=c(1,2)) plot(x,f1,type="l");plot(x,f2,type="l") e <- rnorm(n, 0, 1.5) fac <- as.factor(sample(1:2,n,replace=TRUE)) fac.1 <- as.numeric(fac==1) fac.2 <- as.numeric(fac==2) y <- f1*fac.1 + f2*fac.2 + e dat <- data.frame(y=y,x=x,fac=fac,f1=f1,f2=f2) b2 <- scam(y ~ fac+s(x,by=fac,bs="cv"),data=dat,optimizer="efs") plot(b2,pages=1,scale=0,shade=TRUE) summary(b2) x11() vis.scam(b2,theta=50,color="terrain") ## numeric 'by' variable example... set.seed(6) n <- 100 x <- sort(2*runif(n)-1) z <- runif(n,-2,3) f <- -4*x^2 y <- f*z + rnorm(n)*0.6 dat <- data.frame(x=x,z=z,y=y) b <- scam(y~s(x,k=15,by=z,bs="cvBy"),data=dat) summary(b) par(mfrow=c(1,2)) plot(b,shade=TRUE) ## unconstrained fit... b1 <- scam(y~s(x,k=15,by=z),data=dat) plot(b1,shade=TRUE) summary(b1) ## End(Not run)
This is a special method function
for creating smooths subject to convexity constraint which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed using convex P-splines. This smooth is specified via model terms such as
s(x,k,bs="cx",m=2),
where k denotes the basis dimension and m+1 is the order of the B-spline basis.
cxBy.smooth.spec works similar to cx.smooth.spec but without applying an identifiability constraint ('zero intercept' constraint). cxBy.smooth.spec should be used when the smooth term has a numeric by variable that takes more than one value. In such cases, the smooth terms are fully identifiable without a 'zero intercept' constraint, so they are left unconstrained. This smooth is specified as
s(x,by=z,bs="cxBy"). See an example below.
However a factor by variable requires identifiability constraints, so s(x,by=fac,bs="cx") is used in this case.
## S3 method for class 'cx.smooth.spec' smooth.construct(object, data, knots) ## S3 method for class 'cxBy.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'cx.smooth.spec' smooth.construct(object, data, knots) ## S3 method for class 'cxBy.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the data required by this term,
with names given by |
knots |
An optional list containing the knots supplied for basis setup.
If it is |
An object of class "cx.smooth", "cxBy.smooth".
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
smooth.construct.cv.smooth.spec,
smooth.construct.mpi.smooth.spec,
smooth.construct.mdcv.smooth.spec,
smooth.construct.mdcx.smooth.spec,
smooth.construct.micv.smooth.spec,
smooth.construct.mpd.smooth.spec
## Not run: ## Convex SCOP-splines example... ## simulating data... require(scam) set.seed(16) n <- 100 x <- sort(2*runif(n)-1) f <- 4*x^2 y <- f + rnorm(n)*0.4 dat <- data.frame(x=x,y=y) b <- scam(y~s(x,k=15,bs="cx"),family=gaussian,data=dat) ## unconstrained fit... b1 <- scam(y~s(x,k=15),family=gaussian, data=dat) ## plot results ... plot(x,y,xlab="x",ylab="y") lines(x,f) ## the true function lines(x,b$fitted,col=2) ## constrained fit lines(x,b1$fitted,col=3) ## unconstrained fit ## Poisson version... set.seed(18) y <- rpois(n,exp(f)) dat <- data.frame(x=x,y=y) ## fit shape constrained model ... b <- scam(y~s(x,k=15,bs="cx"),family=poisson(link="log"),data=dat,optimizer="efs") ## unconstrained fit... b1 <- scam(y~s(x,k=15),family=poisson(link="log"), data=dat,optimizer="efs") ## plot results ... plot(x,y,xlab="x",ylab="y") lines(x,exp(f)) ## the true function lines(x,b$fitted,col=2) ## constrained fit lines(x,b1$fitted,col=3) ## unconstrained fit ## 'by' factor example... set.seed(9) n <- 400 x <- sort(runif(n,-.5,.5)) f1 <- .7*x-cos(x)+3 f2 <- 20*x^2 par(mfrow=c(1,2)) plot(x,f1,type="l");plot(x,f2,type="l") e <- rnorm(n, 0, 1.5) fac <- as.factor(sample(1:2,n,replace=TRUE)) fac.1 <- as.numeric(fac==1) fac.2 <- as.numeric(fac==2) y <- f1*fac.1 + f2*fac.2 + e dat <- data.frame(y=y,x=x,fac=fac,f1=f1,f2=f2) b2 <- scam(y ~ fac+s(x,by=fac,bs="cx"),data=dat,optimizer="efs") plot(b2,pages=1,scale=0) summary(b2) x11() vis.scam(b2,theta=50,color="terrain") ## numeric 'by' variable example... set.seed(6) n <- 100 x <- sort(2*runif(n)-1) z <- runif(n,-2,3) f <- 4*x^2 y <- f*z + rnorm(n)*.6 dat <- data.frame(x=x,z=z,y=y) b <- scam(y~s(x,k=15,by=z,bs="cxBy"),data=dat) summary(b) par(mfrow=c(1,2)) plot(b,shade=TRUE) ## unconstrained fit... b1 <- scam(y~s(x,k=15,by=z),data=dat) plot(b1,shade=TRUE) summary(b1) ## End(Not run)## Not run: ## Convex SCOP-splines example... ## simulating data... require(scam) set.seed(16) n <- 100 x <- sort(2*runif(n)-1) f <- 4*x^2 y <- f + rnorm(n)*0.4 dat <- data.frame(x=x,y=y) b <- scam(y~s(x,k=15,bs="cx"),family=gaussian,data=dat) ## unconstrained fit... b1 <- scam(y~s(x,k=15),family=gaussian, data=dat) ## plot results ... plot(x,y,xlab="x",ylab="y") lines(x,f) ## the true function lines(x,b$fitted,col=2) ## constrained fit lines(x,b1$fitted,col=3) ## unconstrained fit ## Poisson version... set.seed(18) y <- rpois(n,exp(f)) dat <- data.frame(x=x,y=y) ## fit shape constrained model ... b <- scam(y~s(x,k=15,bs="cx"),family=poisson(link="log"),data=dat,optimizer="efs") ## unconstrained fit... b1 <- scam(y~s(x,k=15),family=poisson(link="log"), data=dat,optimizer="efs") ## plot results ... plot(x,y,xlab="x",ylab="y") lines(x,exp(f)) ## the true function lines(x,b$fitted,col=2) ## constrained fit lines(x,b1$fitted,col=3) ## unconstrained fit ## 'by' factor example... set.seed(9) n <- 400 x <- sort(runif(n,-.5,.5)) f1 <- .7*x-cos(x)+3 f2 <- 20*x^2 par(mfrow=c(1,2)) plot(x,f1,type="l");plot(x,f2,type="l") e <- rnorm(n, 0, 1.5) fac <- as.factor(sample(1:2,n,replace=TRUE)) fac.1 <- as.numeric(fac==1) fac.2 <- as.numeric(fac==2) y <- f1*fac.1 + f2*fac.2 + e dat <- data.frame(y=y,x=x,fac=fac,f1=f1,f2=f2) b2 <- scam(y ~ fac+s(x,by=fac,bs="cx"),data=dat,optimizer="efs") plot(b2,pages=1,scale=0) summary(b2) x11() vis.scam(b2,theta=50,color="terrain") ## numeric 'by' variable example... set.seed(6) n <- 100 x <- sort(2*runif(n)-1) z <- runif(n,-2,3) f <- 4*x^2 y <- f*z + rnorm(n)*.6 dat <- data.frame(x=x,z=z,y=y) b <- scam(y~s(x,k=15,by=z,bs="cxBy"),data=dat) summary(b) par(mfrow=c(1,2)) plot(b,shade=TRUE) ## unconstrained fit... b1 <- scam(y~s(x,k=15,by=z),data=dat) plot(b1,shade=TRUE) summary(b1) ## End(Not run)
scam can use locally shape-constrained smooths that are monotone increasing up to a specified change point and become unconstrained beyond that point. This smooth is specified via model terms like
s(x,bs="lmpi",xt=list(xc=xc)), where xc sets a change point. The construction of the smooths uses B-spline bases with mildly non-linear re-parametrization of the basis coefficients over the shape-constrained interval. The 'wiggliness' of the smooths is controlled by discrete penalties applied directly to the basis coefficients, in particular, by penalizing the squared differences between adjacent spline coefficients.
The method is build by the mgcv constructor function for smooth terms, smooth.construct.
## S3 method for class 'lmpi.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'lmpi.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the data required by this term,
with names given by |
knots |
An optional list containing the knots supplied for basis setup.
If it is |
A smooth term of the form s(x,k,bs="lmpi",xt=list(xc=xc),m=2), where xc specifies a change point. k denotes the basis dimension and m+1 is the order of the B-spline basis. The default is a cubic spline m=2. The smooth uses a first order difference penalty on the coefficients starting from the second coefficient over the constrained part, and a second order difference penalty (a second order P-spline penalty) over the unconstrained part. The first order difference penalty over the constrained part shares with a second order P-spline penalty the feature of 'smoothing towards a straight line'.
The default basis dimension is k=10. The basis dimensions for constrained and unconstrained parts are set
proportionally to the share of each part from the total range of the covariate value, but with a minimum of 5 basis functions on each side.
The smooth currently does not support user-supplied knots. The knots are placed evenly throughout the constrained and unconstrained ranges of the covariate values (inner knots) and outside the covariate values range (outer knots). The change point xc is set as an inner knot m times. Having multiple knots at the change point guarantees the correct imposing of the shape constraints up to the change point, avoiding 'leakage' of the transformed basis functions of the constrained part into the unconstrained part.
Sum-to-zero ('centering') constraint is applied to the LSCOP-splines after imposing the shape-constrained model matrix transformation. Linear extrapolation is (only) used for prediction that requires extrapolation, i.e. prediction outside the range of the interior knots.
An object of class "lmpi.smooth".
Natalya Pya <[email protected]>
with contributions from Jens Lichter
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
smooth.construct.mpd.smooth.spec,
smooth.construct.cv.smooth.spec,
smooth.construct.cx.smooth.spec,
smooth.construct.mdcv.smooth.spec,
smooth.construct.mdcx.smooth.spec,
smooth.construct.micv.smooth.spec,
smooth.construct.micx.smooth.spec
## Local monotone increasing LSCOP-spline example... ## simulating data... require(scam) set.seed(4) n <- 200 x <- sort(runif(n)*6) f <- 4*tanh(2*x-5)-5*exp(-(x-4)^2) ## increasing up untill xt=2.978 xc <- 2.978 y <- f+rnorm(n)*.7 old.par <- par(mfrow=c(1,2)) plot(x,y,cex=.5) lines(x,f,lty=2); abline(v=xc,lty=2,col="red") ## fit model ... b <- scam(y~s(x,bs="lmpi",xt=list(xc=xc),k=15)) summary(b) lines(x,fitted(b),col=2,lwd=2) plot(b,shade=TRUE) par(old.par) ## unconstrained model to compare with... g <- scam(y~s(x)) plot(x,y,cex=.5) lines(x,f,lty=2); abline(v=xc,lty=2,col=2) lines(x,fitted(g),col=4,lwd=2) lines(x,fitted(b),col=2,lwd=1)## Local monotone increasing LSCOP-spline example... ## simulating data... require(scam) set.seed(4) n <- 200 x <- sort(runif(n)*6) f <- 4*tanh(2*x-5)-5*exp(-(x-4)^2) ## increasing up untill xt=2.978 xc <- 2.978 y <- f+rnorm(n)*.7 old.par <- par(mfrow=c(1,2)) plot(x,y,cex=.5) lines(x,f,lty=2); abline(v=xc,lty=2,col="red") ## fit model ... b <- scam(y~s(x,bs="lmpi",xt=list(xc=xc),k=15)) summary(b) lines(x,fitted(b),col=2,lwd=2) plot(b,shade=TRUE) par(old.par) ## unconstrained model to compare with... g <- scam(y~s(x)) plot(x,y,cex=.5) lines(x,f,lty=2); abline(v=xc,lty=2,col=2) lines(x,fitted(g),col=4,lwd=2) lines(x,fitted(b),col=2,lwd=1)
This is a special method function
for creating smooths subject to both monotone decreasing and concavity constraints which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed using mixed constrained P-splines. This smooth is specified via model terms such as
s(x,k,bs="mdcv",m=2),
where k denotes the basis dimension and m+1 is the order of the B-spline basis.
mdcvBy.smooth.spec works similar to mdcv.smooth.spec but without applying an identifiability constraint ('zero intercept' constraint). mdcvBy.smooth.spec should be used when the smooth term has a numeric by variable that takes more than one value. In such cases, the smooth terms are fully identifiable without a 'zero intercept' constraint, so they are left unconstrained. This smooth is specified as
s(x,by=z,bs="mdcvBy"). See an example below.
However a factor by variable requires identifiability constraints, so s(x,by=fac,bs="mdcv") is used in this case.
## S3 method for class 'mdcv.smooth.spec' smooth.construct(object, data, knots) ## S3 method for class 'mdcvBy.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'mdcv.smooth.spec' smooth.construct(object, data, knots) ## S3 method for class 'mdcvBy.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the data required by this term,
with names given by |
knots |
An optional list containing the knots supplied for basis setup.
If it is |
An object of class "mdcv.smooth", "mdcvBy.smooth".
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
smooth.construct.mpi.smooth.spec,
smooth.construct.mpd.smooth.spec,
smooth.construct.cx.smooth.spec,
smooth.construct.cv.smooth.spec,
smooth.construct.mdcx.smooth.spec,
smooth.construct.micx.smooth.spec,
smooth.construct.micv.smooth.spec
## Not run: ## Monotone decreasing and concave SCOP-splines example ## simulating data... require(scam) set.seed(2) n <- 100 x <- sort(runif(n)) f <- -x^4 y <- f+rnorm(n)*.2 dat <- data.frame(x=x,y=y) ## fit model ... b <- scam(y~s(x,bs="mdcv"),family=gaussian(),data=dat) ## fit unconstrained model ... b1 <- scam(y~s(x,bs="ps"),family=gaussian(),data=dat) ## plot results ... plot(x,y,xlab="x",ylab="y",cex=.5) lines(x,f) ## the true function lines(x,b$fitted.values,col=2) ## mixed constrained fit lines(x,b1$fitted.values,col=3) ## unconstrained fit ## numeric 'by' variable example... set.seed(6) n <- 100 x <- sort(runif(n)) z <- runif(n,-2,3) f <- -x^4 y <- f*z + rnorm(n)*0.4 dat <- data.frame(x=x,z=z,y=y) b <- scam(y~s(x,k=15,by=z,bs="mdcvBy"),data=dat) summary(b) par(mfrow=c(1,2)) plot(b,shade=TRUE) ## unconstrained fit... b1 <- scam(y~s(x,k=15,by=z),data=dat) plot(b1,shade=TRUE) summary(b1) ## End(Not run)## Not run: ## Monotone decreasing and concave SCOP-splines example ## simulating data... require(scam) set.seed(2) n <- 100 x <- sort(runif(n)) f <- -x^4 y <- f+rnorm(n)*.2 dat <- data.frame(x=x,y=y) ## fit model ... b <- scam(y~s(x,bs="mdcv"),family=gaussian(),data=dat) ## fit unconstrained model ... b1 <- scam(y~s(x,bs="ps"),family=gaussian(),data=dat) ## plot results ... plot(x,y,xlab="x",ylab="y",cex=.5) lines(x,f) ## the true function lines(x,b$fitted.values,col=2) ## mixed constrained fit lines(x,b1$fitted.values,col=3) ## unconstrained fit ## numeric 'by' variable example... set.seed(6) n <- 100 x <- sort(runif(n)) z <- runif(n,-2,3) f <- -x^4 y <- f*z + rnorm(n)*0.4 dat <- data.frame(x=x,z=z,y=y) b <- scam(y~s(x,k=15,by=z,bs="mdcvBy"),data=dat) summary(b) par(mfrow=c(1,2)) plot(b,shade=TRUE) ## unconstrained fit... b1 <- scam(y~s(x,k=15,by=z),data=dat) plot(b1,shade=TRUE) summary(b1) ## End(Not run)
This is a special method function
for creating smooths subject to both monotone decreasing and convexity constraints which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed using mixed constrained P-splines. This smooth is specified via model terms such as
s(x,k,bs="mdcx",m=2),
where k denotes the basis dimension and m+1 is the order of the B-spline basis.
mdcxBy.smooth.spec works similar to mdcx.smooth.spec but without applying an identifiability constraint ('zero intercept' constraint). mdcxBy.smooth.spec should be used when the smooth term has a numeric by variable that takes more than one value. In such cases, the smooth terms are fully identifiable without a 'zero intercept' constraint, so they are left unconstrained. This smooth is specified as
s(x,by=z,bs="mdcxBy"). See an example below.
However a factor by variable requires identifiability constraints, so s(x,by=fac,bs="mdcx") is used in this case.
## S3 method for class 'mdcx.smooth.spec' smooth.construct(object, data, knots) ## S3 method for class 'mdcxBy.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'mdcx.smooth.spec' smooth.construct(object, data, knots) ## S3 method for class 'mdcxBy.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the data required by this term,
with names given by |
knots |
An optional list containing the knots supplied for basis setup.
If it is |
An object of class "mdcx.smooth", "mdcxBy.smooth".
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
smooth.construct.mpi.smooth.spec,
smooth.construct.mpd.smooth.spec,
smooth.construct.cx.smooth.spec,
smooth.construct.cv.smooth.spec,
smooth.construct.mdcv.smooth.spec,
smooth.construct.micx.smooth.spec,
smooth.construct.micv.smooth.spec
## Not run: ## Monotone decreasing and convex SCOP-splines example ## simulating data... require(scam) set.seed(2) n <- 100 x <- sort(runif(n)*3-1) f <- (x-3)^6/1000 # monotone decreasing and convex smooth y <- f+rnorm(n)*.4 dat <- data.frame(x=x,y=y) ## fit model ... b <- scam(y~s(x,k=15,bs="mdcx"),family=gaussian(link="identity"),data=dat) ## fit unconstrained model ... b1 <- scam(y~s(x,k=15,bs="ps"),family=gaussian(link="identity"),data=dat) ## plot results ... plot(x,y,xlab="x",ylab="y") lines(x,f) ## the true function lines(x,b$fitted.values,col=2) ## mixed constrained fit lines(x,b1$fitted.values,col=3) ## unconstrained fit ## numeric 'by' variable example... set.seed(6) n <- 100 x <- sort(runif(n)*3-1) z <- runif(n,-2,3) f <- (x-3)^6/1000 y <- f*z + rnorm(n)*.4 dat <- data.frame(x=x,z=z,y=y) b <- scam(y~s(x,k=15,by=z,bs="mdcxBy"),data=dat) summary(b) par(mfrow=c(1,2)) plot(b,shade=TRUE) ## unconstrained fit... b1 <- scam(y~s(x,k=15,by=z),data=dat) plot(b1,shade=TRUE) summary(b1) ## End(Not run)## Not run: ## Monotone decreasing and convex SCOP-splines example ## simulating data... require(scam) set.seed(2) n <- 100 x <- sort(runif(n)*3-1) f <- (x-3)^6/1000 # monotone decreasing and convex smooth y <- f+rnorm(n)*.4 dat <- data.frame(x=x,y=y) ## fit model ... b <- scam(y~s(x,k=15,bs="mdcx"),family=gaussian(link="identity"),data=dat) ## fit unconstrained model ... b1 <- scam(y~s(x,k=15,bs="ps"),family=gaussian(link="identity"),data=dat) ## plot results ... plot(x,y,xlab="x",ylab="y") lines(x,f) ## the true function lines(x,b$fitted.values,col=2) ## mixed constrained fit lines(x,b1$fitted.values,col=3) ## unconstrained fit ## numeric 'by' variable example... set.seed(6) n <- 100 x <- sort(runif(n)*3-1) z <- runif(n,-2,3) f <- (x-3)^6/1000 y <- f*z + rnorm(n)*.4 dat <- data.frame(x=x,z=z,y=y) b <- scam(y~s(x,k=15,by=z,bs="mdcxBy"),data=dat) summary(b) par(mfrow=c(1,2)) plot(b,shade=TRUE) ## unconstrained fit... b1 <- scam(y~s(x,k=15,by=z),data=dat) plot(b1,shade=TRUE) summary(b1) ## End(Not run)
This is a special method function
for creating smooths subject to both monotone increasing and concavity constraints which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed using mixed constrained P-splines. This smooth is specified via model terms such as
s(x,k,bs="micv",m=2),
where k denotes the basis dimension and m+1 is the order of the B-spline basis.
micvBy.smooth.spec works similar to micv.smooth.spec but without applying an identifiability constraint ('zero intercept' constraint). micvBy.smooth.spec should be used when the smooth term has a numeric by variable that takes more than one value. In such cases, the smooth terms are fully identifiable without a 'zero intercept' constraint, so they are left unconstrained. This smooth is specified as
s(x,by=z,bs="micvBy"). See an example below.
However a factor by variable requires identifiability constraints, so s(x,by=fac,bs="micv") is used in this case.
## S3 method for class 'micv.smooth.spec' smooth.construct(object, data, knots) ## S3 method for class 'micvBy.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'micv.smooth.spec' smooth.construct(object, data, knots) ## S3 method for class 'micvBy.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the data required by this term,
with names given by |
knots |
An optional list containing the knots supplied for basis setup.
If it is |
An object of class "micv.smooth", "micvBy.smooth".
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
smooth.construct.mpi.smooth.spec,
smooth.construct.cx.smooth.spec,
smooth.construct.cv.smooth.spec,
smooth.construct.mdcv.smooth.spec,
smooth.construct.mdcx.smooth.spec,
smooth.construct.micx.smooth.spec,
smooth.construct.mpd.smooth.spec
## Not run: ## Monotone increasing and concave SCOP-splines example ## simulating data... set.seed(3) n <- 100 x <- sort(runif(n)*99+1) f <- log(x)/2 y <- f+rnorm(n)*.3 dat <- data.frame(x=x,y=y) ## fit model ... b <- scam(y~s(x,k=15,bs="micv"), data=dat) summary(b) # fit unconstrained model ... b1 <- scam(y~s(x,k=15,bs="ps"),data=dat) ## plot results ... plot(x,y,xlab="x",ylab="y",cex=.5) lines(x,f) ## the true function lines(x,b$fitted.values,col=2) ## mixed constrained fit lines(x,b1$fitted.values,col=3) ## unconstrained fit ## numeric 'by' variable example... set.seed(3) n <- 100 x <- sort(runif(n)*99+1) f <- log(x)/2 z <- runif(n,-2,3) y <- f*z + rnorm(n)*0.3 dat <- data.frame(x=x,z=z,y=y) b <- scam(y~s(x,k=15,by=z,bs="micvBy")-1,data=dat) summary(b) par(mfrow=c(1,2)) plot(b,shade=TRUE) ## unconstrained fit... b1 <- scam(y~s(x,k=15,by=z)-1,data=dat) plot(b1,shade=TRUE) summary(b1) ## End(Not run)## Not run: ## Monotone increasing and concave SCOP-splines example ## simulating data... set.seed(3) n <- 100 x <- sort(runif(n)*99+1) f <- log(x)/2 y <- f+rnorm(n)*.3 dat <- data.frame(x=x,y=y) ## fit model ... b <- scam(y~s(x,k=15,bs="micv"), data=dat) summary(b) # fit unconstrained model ... b1 <- scam(y~s(x,k=15,bs="ps"),data=dat) ## plot results ... plot(x,y,xlab="x",ylab="y",cex=.5) lines(x,f) ## the true function lines(x,b$fitted.values,col=2) ## mixed constrained fit lines(x,b1$fitted.values,col=3) ## unconstrained fit ## numeric 'by' variable example... set.seed(3) n <- 100 x <- sort(runif(n)*99+1) f <- log(x)/2 z <- runif(n,-2,3) y <- f*z + rnorm(n)*0.3 dat <- data.frame(x=x,z=z,y=y) b <- scam(y~s(x,k=15,by=z,bs="micvBy")-1,data=dat) summary(b) par(mfrow=c(1,2)) plot(b,shade=TRUE) ## unconstrained fit... b1 <- scam(y~s(x,k=15,by=z)-1,data=dat) plot(b1,shade=TRUE) summary(b1) ## End(Not run)
This is a special method function
for creating smooths subject to both monotone increasing and convexity constraints which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed using the mixed constrained P-splines. This smooth is specified via model terms such as
s(x,k,bs="micx",m=2),
where k denotes the basis dimension and m+1 is the order of the B-spline basis.
micxBy.smooth.spec works similar to micx.smooth.spec but without applying an identifiability constraint ('zero intercept' constraint). micxBy.smooth.spec should be used when the smooth term has a numeric by variable that takes more than one value. In such cases, the smooth terms are fully identifiable without a 'zero intercept' constraint, so they are left unconstrained. This smooth is specified as
s(x,by=z,bs="micvBy"). See an example below.
However a factor by variable requires identifiability constraints, so s(x,by=fac,bs="micx") is used in this case.
## S3 method for class 'micx.smooth.spec' smooth.construct(object, data, knots) ## S3 method for class 'micxBy.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'micx.smooth.spec' smooth.construct(object, data, knots) ## S3 method for class 'micxBy.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the data required by this term,
with names given by |
knots |
An optional list containing the knots supplied for basis setup.
If it is |
An object of class "micx.smooth", "micxBy.smooth".
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
smooth.construct.mpi.smooth.spec,
smooth.construct.cx.smooth.spec,
smooth.construct.cv.smooth.spec,
smooth.construct.mdcv.smooth.spec,
smooth.construct.mdcx.smooth.spec,
smooth.construct.micv.smooth.spec,
smooth.construct.mpd.smooth.spec
## Not run: ## Monotone increasing and convex SCOP-splines example ## simulating data... set.seed(1) n <- 100 x <- runif(n)*2 f <- 5*x^2/8 y <- rpois(n,exp(f)) dat <- data.frame(x=x,y=y) ## fit model ... b <- scam(y~s(x,bs="micx"),family=poisson,data=dat) ## fit unconstrained model ... b1 <- scam(y~s(x,bs="cr"),family=poisson,data=dat) ## plot results ... plot(x,y,xlab="x",ylab="y",cex=.5) x1 <- sort(x,index=TRUE) lines(x1$x,exp(f)[x1$ix]) ## the true function lines(x1$x,b$fitted.values[x1$ix],col=2) ## mixed constrained fit lines(x1$x,b1$fitted.values[x1$ix],col=3) ## unconstrained fit ## numeric 'by' variable example... set.seed(10) n <- 100 x <- runif(n)*2 f <- x^2 z <- runif(n,-2,3) y <- f*z + rnorm(n)*0.4 dat <- data.frame(x=x,z=z,y=y) b <- scam(y~s(x,by=z,bs="micxBy"),data=dat) summary(b) par(mfrow=c(1,2)) plot(b,shade=TRUE) ## unconstrained fit... b1 <- scam(y~s(x,by=z),data=dat) plot(b1,shade=TRUE) summary(b1) ## End(Not run)## Not run: ## Monotone increasing and convex SCOP-splines example ## simulating data... set.seed(1) n <- 100 x <- runif(n)*2 f <- 5*x^2/8 y <- rpois(n,exp(f)) dat <- data.frame(x=x,y=y) ## fit model ... b <- scam(y~s(x,bs="micx"),family=poisson,data=dat) ## fit unconstrained model ... b1 <- scam(y~s(x,bs="cr"),family=poisson,data=dat) ## plot results ... plot(x,y,xlab="x",ylab="y",cex=.5) x1 <- sort(x,index=TRUE) lines(x1$x,exp(f)[x1$ix]) ## the true function lines(x1$x,b$fitted.values[x1$ix],col=2) ## mixed constrained fit lines(x1$x,b1$fitted.values[x1$ix],col=3) ## unconstrained fit ## numeric 'by' variable example... set.seed(10) n <- 100 x <- runif(n)*2 f <- x^2 z <- runif(n,-2,3) y <- f*z + rnorm(n)*0.4 dat <- data.frame(x=x,z=z,y=y) b <- scam(y~s(x,by=z,bs="micxBy"),data=dat) summary(b) par(mfrow=c(1,2)) plot(b,shade=TRUE) ## unconstrained fit... b1 <- scam(y~s(x,by=z),data=dat) plot(b1,shade=TRUE) summary(b1) ## End(Not run)
This is a special method function
for creating smooths subject to a monotone increasing constraint plus the smooths should pass through zero at the right-end point of the covariate range. This is similar to the pc argument to s in mgcv(gam) when pc=max(x), where x is a covariate.
The smooth is built by the mgcv constructor function for smooth terms, smooth.construct. 'Zero intercept' identifiability constraints used for univariate SCOP-splines are substituted with a 'finish at zero' constraint here. This smooth is specified via model terms such as
s(x,k,bs="mifo",m=2),
where k denotes the basis dimension and m+1 is the order of the B-spline basis.
## S3 method for class 'mifo.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'mifo.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the data required by this term,
with names given by |
knots |
An optional list containing the knots supplied for basis setup.
If it is |
The constructor is not called directly, but as with gam(mgcv) is used internally.
A 'finish at zero' constraint is achieved by setting the last (m+1) spline coefficients to zero. According to the B-spline basis functions properties, the value of the spline, f(x), is determined by m+2 non-zero basis functions, and only m+1 B-splines are non-zero at knots. Only m+2 B-splines are non-zero on any [k_i, k_{i+1}), and the sum of these m+2 basis functions is 1.
If the knots of the spline are not supplied, then they are placed evenly throughout the covariate values with an exception of the m inner knots preceeding the last inner knot that are joined with that last knot. This is done in order to avoid an otherwise plateau fit at the right-end region. If the knots are supplied, then the number of supplied knots should be k+m+2, and the range of the middle k-m knots must include all the covariate values.
Note: when a plateau region is expected at the righ-end covariate region, the smooth might result in some decrease when approaching to zero.
An object of class "mifo.smooth".
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
smooth.construct.mpi.smooth.spec,
smooth.construct.miso.smooth.spec,
smooth.construct.mpd.smooth.spec,
smooth.construct.mdcv.smooth.spec,
smooth.construct.mdcx.smooth.spec, smooth.construct.micv.smooth.spec,
smooth.construct.micx.smooth.spec
## Monotone increasing SCOP-spline examples with a finish at zero constraint... set.seed(53) n <- 100;x <- runif(n);z <- runif(n) pc <- max(x) y <- exp(3*x)/10-exp(3*pc)/10 + z*(1-z)*5 + rnorm(100)*.4 m1 <- scam(y~s(x,bs='mifo')+s(z)) #,knots=knots) plot(m1,pages=1,scale=0) summary(m1) newd<- data.frame(x=pc,z=0) predict(m1,newd, type='terms')## Monotone increasing SCOP-spline examples with a finish at zero constraint... set.seed(53) n <- 100;x <- runif(n);z <- runif(n) pc <- max(x) y <- exp(3*x)/10-exp(3*pc)/10 + z*(1-z)*5 + rnorm(100)*.4 m1 <- scam(y~s(x,bs='mifo')+s(z)) #,knots=knots) plot(m1,pages=1,scale=0) summary(m1) newd<- data.frame(x=pc,z=0) predict(m1,newd, type='terms')
This is a special method function
for creating smooths subject to a monotone increasing constraint plus the smooths should pass through zero at the left-end point of the covariate range. This is similar to the pc argument to s in mgcv(gam) when pc=min(x), where x is a covariate.
The smooth is built by the mgcv constructor function for smooth terms, smooth.construct. 'Zero intercept' identifiability constraints used for univariate SCOP-splines are substituted with a 'start at zero' constraint here. This smooth is specified via model terms such as
s(x,k,bs="miso",m=2),
where k denotes the basis dimension and m+1 is the order of the B-spline basis.
## S3 method for class 'miso.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'miso.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the data required by this term,
with names given by |
knots |
An optional list containing the knots supplied for basis setup.
If it is |
The constructor is not called directly, but as with gam(mgcv) is used internally.
A 'start at zero' constraint is achieved by setting the first (m+1) spline coefficients to zero. According to the B-spline basis functions properties, the value of the spline, f(x), is determined by m+2 non-zero basis functions, and only m+1 B-splines are non-zero at knots. Only m+2 B-splines are non-zero on any [k_i, k_{i+1}), and the sum of these m+2 basis functions is 1.
If the knots of the spline are not supplied, then they are placed evenly throughout the covariate values with an exception of the m inner knots following the first inner knot that are joined with that first knot. This is done in order to avoid an otherwise plateau fit at the left-end region. If the knots are supplied, then the number of supplied knots should be k+m+2, and the range of the middle k-m knots must include all the covariate values.
An object of class "miso.smooth".
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
smooth.construct.mpi.smooth.spec,
smooth.construct.mifo.smooth.spec,
smooth.construct.mpd.smooth.spec,
smooth.construct.mdcv.smooth.spec,
smooth.construct.mdcx.smooth.spec, smooth.construct.micv.smooth.spec,
smooth.construct.micx.smooth.spec
## Monotone increasing SCOP-spline examples with a start at zero constraint... ## passing through 0 at -1... require(scam) set.seed(7) n <- 100; x <- c(-1,runif(n-1)*4-1); ## starting at -1 for a function to be zero at a start z <- runif(n) y <- exp(4*x)/(1+exp(4*x)) -0.01798621+ z*(1-z)*5 + rnorm(100)*.4 m1 <- scam(y~s(x,bs='miso')+s(z)) plot(m1,pages=1) newd<- data.frame(x=-1,z=0) predict(m1,newd, type='terms') ## Not run: ## passing through 0 at 0... set.seed(53) n <- 100; x <- c(0,runif(n-1)); ## starting at 0 for a function to be zero at a start z <- runif(n) y <- exp(3*x)/10-.1 + z*(1-z)*5 + rnorm(100)*.4 m2 <- scam(y~s(x,bs='miso')+s(z)) plot(m2,pages=1) newd<- data.frame(x=0,z=0) predict(m2,newd, type='terms') ## End(Not run)## Monotone increasing SCOP-spline examples with a start at zero constraint... ## passing through 0 at -1... require(scam) set.seed(7) n <- 100; x <- c(-1,runif(n-1)*4-1); ## starting at -1 for a function to be zero at a start z <- runif(n) y <- exp(4*x)/(1+exp(4*x)) -0.01798621+ z*(1-z)*5 + rnorm(100)*.4 m1 <- scam(y~s(x,bs='miso')+s(z)) plot(m1,pages=1) newd<- data.frame(x=-1,z=0) predict(m1,newd, type='terms') ## Not run: ## passing through 0 at 0... set.seed(53) n <- 100; x <- c(0,runif(n-1)); ## starting at 0 for a function to be zero at a start z <- runif(n) y <- exp(3*x)/10-.1 + z*(1-z)*5 + rnorm(100)*.4 m2 <- scam(y~s(x,bs='miso')+s(z)) plot(m2,pages=1) newd<- data.frame(x=0,z=0) predict(m2,newd, type='terms') ## End(Not run)
This is a special method function
for creating smooths subject to monotone decreasing constraints which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed using monotonic P-splines. This smooth is specified via model terms such as
s(x,k,bs="mpd",m=2),
where k denotes the basis dimension and m+1 is the order of the B-spline basis.
mpdBy.smooth.spec works similar to mpd.smooth.spec but without applying an identifiability constraint ('zero intercept' constraint). mpdBy.smooth.spec should be used when the smooth term has a numeric by variable that takes more than one value. In such cases, the smooth terms are fully identifiable without a 'zero intercept' constraint, so they are left unconstrained. This smooth is specified as
s(x,by=z,bs="mpdBy"). See an example below.
However a factor by variable requires identifiability constraints, so s(x,by=fac,bs="mpd") is used in this case.
## S3 method for class 'mpd.smooth.spec' smooth.construct(object, data, knots) ## S3 method for class 'mpdBy.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'mpd.smooth.spec' smooth.construct(object, data, knots) ## S3 method for class 'mpdBy.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the data required by this term,
with names given by |
knots |
An optional list containing the knots supplied for basis setup.
If it is |
An object of class "mpd.smooth", "mpdBy.smooth".
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
smooth.construct.mpi.smooth.spec,
smooth.construct.cx.smooth.spec,
smooth.construct.cv.smooth.spec,
smooth.construct.mdcv.smooth.spec,
smooth.construct.mdcx.smooth.spec,
smooth.construct.micv.smooth.spec,
smooth.construct.micx.smooth.spec
## Not run: ## Monotone decreasing SCOP-splines example... ## simulating data... require(scam) set.seed(3) n <- 100 x <- runif(n)*3-1 f <- exp(-1.3*x) y <- rpois(n,exp(f)) dat <- data.frame(x=x,y=y) ## fit model ... b <- scam(y~s(x,k=15,bs="mpd"),family=poisson(link="log"), data=dat) ## unconstrained model fit for comparison... b1 <- scam(y~s(x,k=15,bs="ps"),family=poisson(link="log"), data=dat) ## plot results ... plot(x,y,xlab="x",ylab="y",cex=.5) x1 <- sort(x,index=TRUE) lines(x1$x,exp(f)[x1$ix]) ## the true function lines(x1$x,b$fitted.values[x1$ix],col=2) ## decreasing fit lines(x1$x,b1$fitted.values[x1$ix],col=3) ## unconstrained fit ## 'by' factor example... set.seed(3) n <- 400 x <- runif(n, 0, 1) ## all three smooths are decreasing... f1 <- -log(x *5) f2 <- -exp(2 * x) + 4 f3 <- -5* sin(x) e <- rnorm(n, 0, 2) fac <- as.factor(sample(1:3,n,replace=TRUE)) fac.1 <- as.numeric(fac==1) fac.2 <- as.numeric(fac==2) fac.3 <- as.numeric(fac==3) y <- f1*fac.1 + f2*fac.2 + f3*fac.3 + e dat <- data.frame(y=y,x=x,fac=fac,f1=f1,f2=f2,f3=f3) b2 <- scam(y ~ fac+s(x,by=fac,bs="mpd"),data=dat) plot(b2,pages=1,scale=0,shade=TRUE) summary(b2) vis.scam(b2,theta=120,color="terrain") ## comparing with unconstrained fit... b3 <- scam(y ~ fac+s(x,by=fac),data=dat) x11() plot(b3,pages=1,scale=0,shade=TRUE) summary(b3) ## Note that since in scam() as in mgcv::gam() when using factor 'by' variables, 'centering' ## constraints are applied to the smooths, which usually means that the 'by' ## factor variable should be included as a parametric term, as well. ## numeric 'by' variable example... set.seed(3) n <- 100 x <- sort(runif(n,-1,2)) z <- runif(n,-2,3) f <- exp(-1.3*x) y <- f*z + rnorm(n)*0.4 dat <- data.frame(x=x,y=y,z=z) b <- scam(y~s(x,k=15,by=z,bs="mpdBy"),data=dat,optimizer="efs") plot(b,shade=TRUE) summary(b) ## unconstrained fit... b1 <- scam(y~s(x,k=15,by=z),data=dat) plot(b1,shade=TRUE) summary(b1) ## End(Not run)## Not run: ## Monotone decreasing SCOP-splines example... ## simulating data... require(scam) set.seed(3) n <- 100 x <- runif(n)*3-1 f <- exp(-1.3*x) y <- rpois(n,exp(f)) dat <- data.frame(x=x,y=y) ## fit model ... b <- scam(y~s(x,k=15,bs="mpd"),family=poisson(link="log"), data=dat) ## unconstrained model fit for comparison... b1 <- scam(y~s(x,k=15,bs="ps"),family=poisson(link="log"), data=dat) ## plot results ... plot(x,y,xlab="x",ylab="y",cex=.5) x1 <- sort(x,index=TRUE) lines(x1$x,exp(f)[x1$ix]) ## the true function lines(x1$x,b$fitted.values[x1$ix],col=2) ## decreasing fit lines(x1$x,b1$fitted.values[x1$ix],col=3) ## unconstrained fit ## 'by' factor example... set.seed(3) n <- 400 x <- runif(n, 0, 1) ## all three smooths are decreasing... f1 <- -log(x *5) f2 <- -exp(2 * x) + 4 f3 <- -5* sin(x) e <- rnorm(n, 0, 2) fac <- as.factor(sample(1:3,n,replace=TRUE)) fac.1 <- as.numeric(fac==1) fac.2 <- as.numeric(fac==2) fac.3 <- as.numeric(fac==3) y <- f1*fac.1 + f2*fac.2 + f3*fac.3 + e dat <- data.frame(y=y,x=x,fac=fac,f1=f1,f2=f2,f3=f3) b2 <- scam(y ~ fac+s(x,by=fac,bs="mpd"),data=dat) plot(b2,pages=1,scale=0,shade=TRUE) summary(b2) vis.scam(b2,theta=120,color="terrain") ## comparing with unconstrained fit... b3 <- scam(y ~ fac+s(x,by=fac),data=dat) x11() plot(b3,pages=1,scale=0,shade=TRUE) summary(b3) ## Note that since in scam() as in mgcv::gam() when using factor 'by' variables, 'centering' ## constraints are applied to the smooths, which usually means that the 'by' ## factor variable should be included as a parametric term, as well. ## numeric 'by' variable example... set.seed(3) n <- 100 x <- sort(runif(n,-1,2)) z <- runif(n,-2,3) f <- exp(-1.3*x) y <- f*z + rnorm(n)*0.4 dat <- data.frame(x=x,y=y,z=z) b <- scam(y~s(x,k=15,by=z,bs="mpdBy"),data=dat,optimizer="efs") plot(b,shade=TRUE) summary(b) ## unconstrained fit... b1 <- scam(y~s(x,k=15,by=z),data=dat) plot(b1,shade=TRUE) summary(b1) ## End(Not run)
This is a special method function
for creating smooths subject to a monotone increasing constraint which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed using monotonic P-splines. This smooth is specified via model terms such as
s(x,k,bs="mpi",m=2),
where k denotes the basis dimension and m+1 is the order of the B-spline basis.
mpiBy.smooth.spec works similar to mpi.smooth.spec but without applying an identifiability constraint ('zero intercept' constraint). mpiBy.smooth.spec should be used when the smooth term has a numeric by variable that takes more than one value. In such cases, the smooth terms are fully identifiable without a 'zero intercept' constraint, so they are left unconstrained. This smooth is specified as
s(x,by=z,bs="mpiBy"). See an example below.
However a factor by variable requires identifiability constraints, so s(x,by=fac,bs="mpi") is used in this case.
## S3 method for class 'mpi.smooth.spec' smooth.construct(object, data, knots) ## S3 method for class 'mpiBy.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'mpi.smooth.spec' smooth.construct(object, data, knots) ## S3 method for class 'mpiBy.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the data required by this term,
with names given by |
knots |
An optional list containing the knots supplied for basis setup.
If it is |
The constructor is not called directly, but as with gam(mgcv) is used internally.
If the knots of the spline are not supplied, then they are placed evenly throughout the covariate values. If the knots are supplied, then the number of supplied knots should be k+m+2, and the range of the middle k-m knots must include all the covariate values.
An object of class "mpi.smooth", "mpiBy.smooth".
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
smooth.construct.mpd.smooth.spec,
smooth.construct.cv.smooth.spec,
smooth.construct.cx.smooth.spec,
smooth.construct.mdcv.smooth.spec,
smooth.construct.mdcx.smooth.spec,
smooth.construct.micv.smooth.spec,
smooth.construct.micx.smooth.spec
## Monotone increasing P-splines example ## simulating data... require(scam) set.seed(12) n <- 100 x <- runif(n)*4-1 f <- 4*exp(4*x)/(1+exp(4*x)) y <- rpois(n,exp(f)) dat <- data.frame(x=x,y=y) ## fit model ... b <- scam(y~s(x,k=15,bs="mpi"),family=poisson(link="log"), data=dat) ## fit unconstrained model... b1 <- scam(y~s(x,k=15,bs="ps"),family=poisson(link="log"), data=dat) ## plot results ... plot(x,y,xlab="x",ylab="y") x1 <- sort(x,index=TRUE) lines(x1$x,exp(f)[x1$ix]) ## the true function lines(x1$x,b$fitted.values[x1$ix],col=2) ## monotone fit lines(x1$x,b1$fitted.values[x1$ix],col=3) ## unconstrained fit ## example with supplied knots... knots <- list(x=c (-1.5, -1.2, -.99, -.97, -.7, -.5, -.3, 0, 0.7, 0.9,1.1, 1.22,1.5,2.2,2.77,2.93,2.99, 3.2,3.6)) b2 <- scam(y~s(x,k=15,bs="mpi"),knots=knots, family=poisson(link="log"), data=dat) summary(b2) plot(b2,shade=TRUE) ## Not run: ## example with two terms... set.seed(0) n <- 200 x1 <- runif(n)*6-3 f1 <- 3*exp(-x1^2) # unconstrained term x2 <- runif(n)*4-1; f2 <- exp(4*x2)/(1+exp(4*x2)) # monotone increasing smooth f <- f1+f2 y <- f+rnorm(n)*.7 dat <- data.frame(x1=x1,x2=x2,y=y) knots <- list(x1=c(-4,-3.5,-2.99,-2.7,-2.5,-1.9,-1.1,-.9,-.3,0.3,.8,1.2,1.9,2.3, 2.7,2.99,3.5,4.1,4.5), x2=c(-1.5,-1.2,-1.1, -.89,-.69,-.5,-.3,0,0.7, 0.9,1.1,1.22,1.5,2.2,2.77,2.99,3.1, 3.2,3.6)) b3 <- scam(y~s(x1,k=15)+s(x2,bs="mpi", k=15), knots=knots,data=dat) summary(b3) plot(b3,pages=1,shade=TRUE) ## setting knots for f(x2) only... knots <- list(x2=c(-1.5,-1.2,-1.1, -.89,-.69,-.5,-.3, 0,0.7,0.9,1.1,1.22,1.5,2.2,2.77,2.99,3.1, 3.2,3.6)) b4 <- scam(y~s(x1,k=15,bs="bs")+s(x2,bs="mpi",k=15), knots=knots,data=dat) summary(b4) plot(b4,pages=1,shade=TRUE) ## 'by' factor example... set.seed(10) n <- 400 x <- runif(n, 0, 1) ## all three smooths are increasing... f1 <- log(x *5) f2 <- exp(2*x) - 4 f3 <- 5* sin(x) e <- rnorm(n, 0, 2) fac <- as.factor(sample(1:3,n,replace=TRUE)) fac.1 <- as.numeric(fac==1) fac.2 <- as.numeric(fac==2) fac.3 <- as.numeric(fac==3) y <- f1*fac.1 + f2*fac.2 + f3*fac.3 + e dat <- data.frame(y=y,x=x,fac=fac,f1=f1,f2=f2,f3=f3) b5 <- scam(y ~ fac+s(x,by=fac,bs="mpi"),data=dat) plot(b5,pages=1,scale=0,shade=TRUE) summary(b5) vis.scam(b5,theta=50,color="terrain") ## comparing with unconstrained fit... b6 <- scam(y ~ fac+s(x,by=fac),data=dat) x11() plot(b6,pages=1,scale=0,shade=TRUE) summary(b6) vis.scam(b6,theta=50,color="terrain") ## Note that since in scam() as in mgcv::gam() when using factor 'by' variables, 'centering' ## constraints are applied to the smooths, which usually means that the 'by' ## factor variable should be included as a parametric term, as well. ## numeric 'by' variable example... set.seed(3) n <- 200 x <- sort(runif(n,-1,2)) z <- runif(n,-2,3) f <- exp(1.3*x)-5 y <- f*z + rnorm(n)*2 dat <- data.frame(x=x,y=y,z=z) b <- scam(y~s(x,by=z,bs="mpiBy"),data=dat) plot(b,shade=TRUE) summary(b) ## unconstrained fit... b1 <- scam(y~s(x,k=15,by=z),data=dat) plot(b1,shade=TRUE) summary(b1) ## End(Not run)## Monotone increasing P-splines example ## simulating data... require(scam) set.seed(12) n <- 100 x <- runif(n)*4-1 f <- 4*exp(4*x)/(1+exp(4*x)) y <- rpois(n,exp(f)) dat <- data.frame(x=x,y=y) ## fit model ... b <- scam(y~s(x,k=15,bs="mpi"),family=poisson(link="log"), data=dat) ## fit unconstrained model... b1 <- scam(y~s(x,k=15,bs="ps"),family=poisson(link="log"), data=dat) ## plot results ... plot(x,y,xlab="x",ylab="y") x1 <- sort(x,index=TRUE) lines(x1$x,exp(f)[x1$ix]) ## the true function lines(x1$x,b$fitted.values[x1$ix],col=2) ## monotone fit lines(x1$x,b1$fitted.values[x1$ix],col=3) ## unconstrained fit ## example with supplied knots... knots <- list(x=c (-1.5, -1.2, -.99, -.97, -.7, -.5, -.3, 0, 0.7, 0.9,1.1, 1.22,1.5,2.2,2.77,2.93,2.99, 3.2,3.6)) b2 <- scam(y~s(x,k=15,bs="mpi"),knots=knots, family=poisson(link="log"), data=dat) summary(b2) plot(b2,shade=TRUE) ## Not run: ## example with two terms... set.seed(0) n <- 200 x1 <- runif(n)*6-3 f1 <- 3*exp(-x1^2) # unconstrained term x2 <- runif(n)*4-1; f2 <- exp(4*x2)/(1+exp(4*x2)) # monotone increasing smooth f <- f1+f2 y <- f+rnorm(n)*.7 dat <- data.frame(x1=x1,x2=x2,y=y) knots <- list(x1=c(-4,-3.5,-2.99,-2.7,-2.5,-1.9,-1.1,-.9,-.3,0.3,.8,1.2,1.9,2.3, 2.7,2.99,3.5,4.1,4.5), x2=c(-1.5,-1.2,-1.1, -.89,-.69,-.5,-.3,0,0.7, 0.9,1.1,1.22,1.5,2.2,2.77,2.99,3.1, 3.2,3.6)) b3 <- scam(y~s(x1,k=15)+s(x2,bs="mpi", k=15), knots=knots,data=dat) summary(b3) plot(b3,pages=1,shade=TRUE) ## setting knots for f(x2) only... knots <- list(x2=c(-1.5,-1.2,-1.1, -.89,-.69,-.5,-.3, 0,0.7,0.9,1.1,1.22,1.5,2.2,2.77,2.99,3.1, 3.2,3.6)) b4 <- scam(y~s(x1,k=15,bs="bs")+s(x2,bs="mpi",k=15), knots=knots,data=dat) summary(b4) plot(b4,pages=1,shade=TRUE) ## 'by' factor example... set.seed(10) n <- 400 x <- runif(n, 0, 1) ## all three smooths are increasing... f1 <- log(x *5) f2 <- exp(2*x) - 4 f3 <- 5* sin(x) e <- rnorm(n, 0, 2) fac <- as.factor(sample(1:3,n,replace=TRUE)) fac.1 <- as.numeric(fac==1) fac.2 <- as.numeric(fac==2) fac.3 <- as.numeric(fac==3) y <- f1*fac.1 + f2*fac.2 + f3*fac.3 + e dat <- data.frame(y=y,x=x,fac=fac,f1=f1,f2=f2,f3=f3) b5 <- scam(y ~ fac+s(x,by=fac,bs="mpi"),data=dat) plot(b5,pages=1,scale=0,shade=TRUE) summary(b5) vis.scam(b5,theta=50,color="terrain") ## comparing with unconstrained fit... b6 <- scam(y ~ fac+s(x,by=fac),data=dat) x11() plot(b6,pages=1,scale=0,shade=TRUE) summary(b6) vis.scam(b6,theta=50,color="terrain") ## Note that since in scam() as in mgcv::gam() when using factor 'by' variables, 'centering' ## constraints are applied to the smooths, which usually means that the 'by' ## factor variable should be included as a parametric term, as well. ## numeric 'by' variable example... set.seed(3) n <- 200 x <- sort(runif(n,-1,2)) z <- runif(n,-2,3) f <- exp(1.3*x)-5 y <- f*z + rnorm(n)*2 dat <- data.frame(x=x,y=y,z=z) b <- scam(y~s(x,by=z,bs="mpiBy"),data=dat) plot(b,shade=TRUE) summary(b) ## unconstrained fit... b1 <- scam(y~s(x,k=15,by=z),data=dat) plot(b1,shade=TRUE) summary(b1) ## End(Not run)
This is a special method function
for creating univariate smooths subject to a positivity constraint which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed using monotonic P-splines. This smooth is specified via model terms suach as
s(x,k,bs="po",m=2),
where k denotes the basis dimension and m+1 is the order of the B-spline basis.
Note: Models that include this smooth should not have an intercept. See examples below.
## S3 method for class 'po.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'po.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the data required by this term,
with names given by |
knots |
An optional list containing the knots supplied for basis setup.
If it is |
An object of class "po.smooth".
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
smooth.construct.mpd.smooth.spec,
smooth.construct.cv.smooth.spec,
smooth.construct.cx.smooth.spec,
smooth.construct.mdcv.smooth.spec,
smooth.construct.mdcx.smooth.spec, smooth.construct.micv.smooth.spec,
smooth.construct.micx.smooth.spec
## SCOP-splines example with positivity constraint... ## simulating data... ## Not run: require(scam) set.seed(3) n <- 100 x <- seq(-3,3,length.out=100) f <- dnorm(x) y <- f + rnorm(n)*0.1 b <- scam(y~s(x,bs="po")-1) b1 <- scam(y~s(x)) ## unconstrained model plot(x,y) lines(x,f) lines(x,fitted(b),col=2) lines(x,fitted(b1),col=3) ## two-term example... set.seed(3) n <- 200 x1 <- seq(-3,3,length.out=n) f1 <- 3*exp(-x1^2) ## positively constrained smooth x2 <- runif(n)*3-1; f2 <- exp(4*x2)/(1+exp(4*x2)) ## increasing smooth f <- f1+f2 y <- f+rnorm(n)*0.3 dat <- data.frame(x1=x1,x2=x2,y=y) ## fit model, results, and plot... b2 <- scam(y~s(x1,bs="po")+s(x2,bs="mpi")-1,data=dat) summary(b2) plot(b2,pages=1) b3 <- scam(y~s(x1,bs="ps")+s(x2,bs="ps"),data=dat) ## unconstrained model summary(b3) plot(b3,pages=1) ## End(Not run)## SCOP-splines example with positivity constraint... ## simulating data... ## Not run: require(scam) set.seed(3) n <- 100 x <- seq(-3,3,length.out=100) f <- dnorm(x) y <- f + rnorm(n)*0.1 b <- scam(y~s(x,bs="po")-1) b1 <- scam(y~s(x)) ## unconstrained model plot(x,y) lines(x,f) lines(x,fitted(b),col=2) lines(x,fitted(b1),col=3) ## two-term example... set.seed(3) n <- 200 x1 <- seq(-3,3,length.out=n) f1 <- 3*exp(-x1^2) ## positively constrained smooth x2 <- runif(n)*3-1; f2 <- exp(4*x2)/(1+exp(4*x2)) ## increasing smooth f <- f1+f2 y <- f+rnorm(n)*0.3 dat <- data.frame(x1=x1,x2=x2,y=y) ## fit model, results, and plot... b2 <- scam(y~s(x1,bs="po")+s(x2,bs="mpi")-1,data=dat) summary(b2) plot(b2,pages=1) b3 <- scam(y~s(x1,bs="ps")+s(x2,bs="ps"),data=dat) ## unconstrained model summary(b3) plot(b3,pages=1) ## End(Not run)
This is a special method function
for creating tensor product bivariate smooths subject to double concavity constraint, i.e. concavity constraint wrt both the first and the second covariates. This is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed from a pair of single penalty marginal smooths which are represented using the B-spline basis functions.
This tensor product is specified by model terms such as s(x1,x2,k=c(q1,q2),bs="tecvcv",m=c(2,2)),
where q1 and q2 denote the basis dimensions for the marginal smooths.
## S3 method for class 'tecvcv.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'tecvcv.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the values of the elements of |
knots |
An optional list containing the knots corresponding to |
An object of class "tecvcv.smooth". In addition to the usual
elements of a smooth class documented under smooth.construct of the mgcv library,
this object contains:
p.ident |
A vector of 0's and 1's for model parameter identification: 1's indicate parameters which will be exponentiated, 0's - otherwise. |
Zc |
A matrix of identifiability constraints. |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
smooth.construct.tedmd.smooth.spec
smooth.construct.temicx.smooth.spec
smooth.construct.tedecx.smooth.spec
smooth.construct.tecxcx.smooth.spec
smooth.construct.tecxcv.smooth.spec
## Not run: ## tensor product `tecvcv' example ## simulating data... set.seed(3) n <- 30 x1 <- sort(2*runif(n)-1) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- -4*(x1[i]^2+x2[j]^2) f <- as.vector(t(f1)) y <- f+rnorm(length(f))*.05 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tecvcv"), data=dat) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b, theta=30,phi=40) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- -4*(x1[i]^2+x2[j]^2) persp(x1,x2,f1,theta = 30, phi = 40) ## End(Not run)## Not run: ## tensor product `tecvcv' example ## simulating data... set.seed(3) n <- 30 x1 <- sort(2*runif(n)-1) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- -4*(x1[i]^2+x2[j]^2) f <- as.vector(t(f1)) y <- f+rnorm(length(f))*.05 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tecvcv"), data=dat) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b, theta=30,phi=40) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- -4*(x1[i]^2+x2[j]^2) persp(x1,x2,f1,theta = 30, phi = 40) ## End(Not run)
This is a special method function
for creating tensor product bivariate smooths subject to mixed constraints, convexity constraint wrt the first covariate and concavity wrt the second one. This is built by the mgcv constructor function for smooth terms, smooth.construct.
It is constructed from a pair of single penalty marginal smooths which are represented using the B-spline basis functions.
This tensor product is specified by model terms such as s(x1,x2,k=c(q1,q2),bs="tecxcv",m=c(2,2)),
where q1 and q2 denote the basis dimensions for the marginal smooths.
## S3 method for class 'tecxcv.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'tecxcv.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the values of the elements of |
knots |
An optional list containing the knots corresponding to |
An object of class "tecxcv.smooth". In addition to the usual
elements of a smooth class documented under smooth.construct of the mgcv library,
this object contains:
p.ident |
A vector of 0's and 1's for model parameter identification: 1's indicate parameters which will be exponentiated, 0's - otherwise. |
Zc |
A matrix of identifiability constraints. |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
smooth.construct.tedmd.smooth.spec
smooth.construct.tedecv.smooth.spec
smooth.construct.tedecx.smooth.spec
smooth.construct.tecvcv.smooth.spec
smooth.construct.tecxcx.smooth.spec
## Not run: ## tensor product `tecxcv' example ## simulating data... set.seed(5) n <- 30 x1 <- sort(2*runif(n)-1) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- 2*x1[i]^2 - 4*x2[j]^2 f <- as.vector(t(f1)) y <- f+rnorm(length(f))*.05 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tecxcv"), data=dat) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b,theta=30,phi=40) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- 2*x1[i]^2 - 4*x2[j]^2 persp(x1,x2,f1,theta = 30, phi = 40) ## End(Not run)## Not run: ## tensor product `tecxcv' example ## simulating data... set.seed(5) n <- 30 x1 <- sort(2*runif(n)-1) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- 2*x1[i]^2 - 4*x2[j]^2 f <- as.vector(t(f1)) y <- f+rnorm(length(f))*.05 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tecxcv"), data=dat) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b,theta=30,phi=40) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- 2*x1[i]^2 - 4*x2[j]^2 persp(x1,x2,f1,theta = 30, phi = 40) ## End(Not run)
This is a special method function
for creating tensor product bivariate smooths subject to double convexity constraint, convexity constraint wrt both the first and the second covariates. This is built by the mgcv constructor function for smooth terms, smooth.construct.
It is constructed from a pair of single penalty marginal smooths which are represented using the B-spline basis functions.
This tensor product is specified by model terms such as s(x1,x2,k=c(q1,q2),bs="tecxcx",m=c(2,2)),
where q1 and q2 denote the basis dimensions for the marginal smooths.
## S3 method for class 'tecxcx.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'tecxcx.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the values of the elements of |
knots |
An optional list containing the knots corresponding to |
An object of class "tecxcx.smooth". In addition to the usual
elements of a smooth class documented under smooth.construct of the mgcv library,
this object contains:
p.ident |
A vector of 0's and 1's for model parameter identification: 1's indicate parameters which will be exponentiated, 0's - otherwise. |
Zc |
A matrix of identifiability constraints. |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
smooth.construct.tedmd.smooth.spec
smooth.construct.tedecv.smooth.spec
smooth.construct.tedecx.smooth.spec
smooth.construct.tecvcv.smooth.spec
smooth.construct.tecxcv.smooth.spec
## Not run: ## tensor product `tecxcx' example ## simulating data... set.seed(2) n <- 30 x1 <- sort(2*runif(n)-1) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) { f1[i,j] <- 2*(x1[i]^2 + x2[j]^2)} f <- as.vector(t(f1)) y <- f+rnorm(length(f))*.05 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tecxcx"), data=dat) summary(b) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b,theta=20,phi=20) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- 2*(x1[i]^2 + x2[j]^2) persp(x1,x2,f1,theta = 30, phi = 40) ## End(Not run)## Not run: ## tensor product `tecxcx' example ## simulating data... set.seed(2) n <- 30 x1 <- sort(2*runif(n)-1) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) { f1[i,j] <- 2*(x1[i]^2 + x2[j]^2)} f <- as.vector(t(f1)) y <- f+rnorm(length(f))*.05 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tecxcx"), data=dat) summary(b) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b,theta=20,phi=20) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- 2*(x1[i]^2 + x2[j]^2) persp(x1,x2,f1,theta = 30, phi = 40) ## End(Not run)
This is a special method function
for creating tensor product bivariate smooths subject to mixed constraints, monotone decreasing constraint wrt the first covariate and concavity wrt the second one, which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed from a pair of single penalty marginal smooths which are represented using the B-spline basis functions.
This tensor product is specified by model terms such as s(x1,x2,k=c(q1,q2),bs="tedecv",m=c(2,2)),
where q1 and q2 denote the basis dimensions for the marginal smooths.
## S3 method for class 'tedecv.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'tedecv.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the values of the elements of |
knots |
An optional list containing the knots corresponding to |
An object of class "tedecv.smooth". In addition to the usual
elements of a smooth class documented under smooth.construct of the mgcv library,
this object contains:
p.ident |
A vector of 0's and 1's for model parameter identification: 1's indicate parameters which will be exponentiated, 0's - otherwise. |
Zc |
A matrix of identifiability constraints. |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
smooth.construct.tedmd.smooth.spec
smooth.construct.temicx.smooth.spec
smooth.construct.tedecx.smooth.spec
## Not run: ## tensor product `tedecv' example ## simulating data... set.seed(2) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) { f1[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i]))- 4*x2[j]^2} f <- as.vector(t(f1)) y <- f+rnorm(length(f))*0.1 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tedecv",m=2), data=dat) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b, theta=30) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i]))- 4*x2[j]^2 persp(x1,x2,f1,theta = 30, phi = 40) ## End(Not run)## Not run: ## tensor product `tedecv' example ## simulating data... set.seed(2) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) { f1[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i]))- 4*x2[j]^2} f <- as.vector(t(f1)) y <- f+rnorm(length(f))*0.1 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tedecv",m=2), data=dat) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b, theta=30) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i]))- 4*x2[j]^2 persp(x1,x2,f1,theta = 30, phi = 40) ## End(Not run)
This is a special method function
for creating tensor product bivariate smooths subject to mixed constraints, monotone decreasing constraint wrt the first covariate and convexity wrt the second one, which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed from a pair of single penalty marginal smooths which are represented using the B-spline basis functions.
This tensor product is specified by model terms such as s(x1,x2,k=c(q1,q2),bs="tedecx",m=c(2,2)),
where q1 and q2 denote the basis dimensions for the marginal smooths.
## S3 method for class 'tedecx.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'tedecx.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the values of the elements of |
knots |
An optional list containing the knots corresponding to |
An object of class "tedecx.smooth". In addition to the usual
elements of a smooth class documented under smooth.construct of the mgcv library,
this object contains:
p.ident |
A vector of 0's and 1's for model parameter identification: 1's indicate parameters which will be exponentiated, 0's - otherwise. |
Zc |
A matrix of identifiability constraints. |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
smooth.construct.tedmd.smooth.spec
smooth.construct.tedecv.smooth.spec
## Not run: ## tensor product `tedecx' example ## simulating data... set.seed(2) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i])) + 2*x2[j]^2 f <- as.vector(t(f1)) y <- f+rnorm(length(f))*0.05 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tedecx",m=2), not.exp=TRUE, data=dat) ## b1 <- scam(y~s(x1,bs="mpd",m=2)+s(x2,bs="cx",m=2), data=dat) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b,theta=20,phi=20) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i])) + 2*x2[j]^2 persp(x1,x2,f1,theta = 30, phi = 40) ## End(Not run)## Not run: ## tensor product `tedecx' example ## simulating data... set.seed(2) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i])) + 2*x2[j]^2 f <- as.vector(t(f1)) y <- f+rnorm(length(f))*0.05 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tedecx",m=2), not.exp=TRUE, data=dat) ## b1 <- scam(y~s(x1,bs="mpd",m=2)+s(x2,bs="cx",m=2), data=dat) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b,theta=20,phi=20) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i])) + 2*x2[j]^2 persp(x1,x2,f1,theta = 30, phi = 40) ## End(Not run)
This is a special method function
for creating tensor product bivariate smooths subject to double monotone decreasing constraints which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed from a pair of single penalty marginal smooths which are represented using the B-spline basis functions.
This tensor product is specified by model terms such as s(x1,x2,k=c(q1,q2),bs="tedmd",m=c(2,2)),
where q1 and q2 denote the basis dimensions for the marginal smooths.
From scam version 1.2-15, the sum-to-zero contraint is now applied to all bivariate SCOP-splines after imposing the scop-constraints (including scop identifiability constraint). This simply shifts the smooth vertically, leaving the shape of the smooths unchanged.
## S3 method for class 'tedmd.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'tedmd.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the values of the elements of |
knots |
An optional list containing the knots corresponding to |
An object of class "tedmd.smooth". In addition to the usual
elements of a smooth class documented under smooth.construct of the mgcv library,
this object contains:
p.ident |
A vector of 0's and 1's for model parameter identification: 1's indicate parameters which will be exponentiated, 0's - otherwise. |
Zc |
A matrix of identifiability constraints. |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
smooth.construct.tedmi.smooth.spec
## Not run: ## tensor product `tedmd' example ## simulating data... require(scam) set.seed(2) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(runif(n)) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) { f1[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i]))-2*exp(x2[j]-0.5)} f <- as.vector(t(f1)) y <- f+rnorm(length(f))*0.1 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tedmd"), data=dat) summary(b) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 80, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") ## End(Not run)## Not run: ## tensor product `tedmd' example ## simulating data... require(scam) set.seed(2) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(runif(n)) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) { f1[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i]))-2*exp(x2[j]-0.5)} f <- as.vector(t(f1)) y <- f+rnorm(length(f))*0.1 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tedmd"), data=dat) summary(b) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 80, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") ## End(Not run)
This is a special method function
for creating tensor product bivariate smooths subject to double monotone increasing constraints which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed from a pair of single penalty marginal smooths which are represented using the B-spline basis functions.
This tensor product is specified by model terms such as s(x1,x2,k=c(q1,q2),bs="tedmi",m=c(2,2)),
where q1 and q2 denote the basis dimensions for the marginal smooths.
## S3 method for class 'tedmi.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'tedmi.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the values of the elements of |
knots |
An optional list containing the knots corresponding to |
An object of class "tedmi.smooth". In addition to the usual
elements of a smooth class documented under smooth.construct of the mgcv library,
this object contains:
p.ident |
A vector of 0's and 1's for model parameter identification: 1's indicate parameters which will be exponentiated, 0's - otherwise. |
Zc |
A matrix of identifiability constraints. |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
smooth.construct.tedmd.smooth.spec
## Not run: ## tensor product `tedmi' example ## simulating data... set.seed(1) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(runif(n)) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) { f1[i,j] <- exp(4*x1[i])/(1+exp(4*x1[i]))+2*exp(x2[j]-0.5)} f <- as.vector(t(f1)) y <- f+rnorm(length(f))*0.1 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tedmi"), data=dat,optimizer="efs") ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") ## End(Not run)## Not run: ## tensor product `tedmi' example ## simulating data... set.seed(1) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(runif(n)) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) { f1[i,j] <- exp(4*x1[i])/(1+exp(4*x1[i]))+2*exp(x2[j]-0.5)} f <- as.vector(t(f1)) y <- f+rnorm(length(f))*0.1 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tedmi"), data=dat,optimizer="efs") ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") ## End(Not run)
This is a special method function
for creating tensor product bivariate smooths subject to mixed constraints, monotone increasing constraint wrt the first covariate and concavity wrt the second one, which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed from a pair of single penalty marginal smooths which are represented using the B-spline basis functions.
This tensor product is specified by model terms such as s(x1,x2,k=c(q1,q2),bs="temicv",m=c(2,2)),
where q1 and q2 denote the basis dimensions for the marginal smooths.
## S3 method for class 'temicv.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'temicv.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the values of the elements of |
knots |
An optional list containing the knots corresponding to |
An object of class "temicv.smooth". In addition to the usual
elements of a smooth class documented under smooth.construct of the mgcv library,
this object contains:
p.ident |
A vector of 0's and 1's for model parameter identification: 1's indicate parameters which will be exponentiated, 0's - otherwise. |
Zc |
A matrix of identifiability constraints. |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
smooth.construct.tedmd.smooth.spec
smooth.construct.temicx.smooth.spec
## Not run: ## tensor product `temicv' example ## simulating data... set.seed(4) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- exp(4*x1[i])/(1+exp(4*x1[i])) - 4*x2[j]^2 f <- as.vector(t(f1)) y <- f+rnorm(length(f))*0.1 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="temicv"), data=dat) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b, theta=30, phi = 40) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- exp(4*x1[i])/(1+exp(4*x1[i])) - 4*x2[j]^2 persp(x1,x2,f1,theta = 30, phi = 40) ## End(Not run)## Not run: ## tensor product `temicv' example ## simulating data... set.seed(4) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- exp(4*x1[i])/(1+exp(4*x1[i])) - 4*x2[j]^2 f <- as.vector(t(f1)) y <- f+rnorm(length(f))*0.1 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="temicv"), data=dat) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b, theta=30, phi = 40) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- exp(4*x1[i])/(1+exp(4*x1[i])) - 4*x2[j]^2 persp(x1,x2,f1,theta = 30, phi = 40) ## End(Not run)
This is a special method function
for creating tensor product bivariate smooths subject to mixed constraints, monotone increasing constraint wrt the first covariate and convexity wrt the second one, which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed from a pair of single penalty marginal smooths which are represented using the B-spline basis functions.
This tensor product is specified by model terms such as s(x1,x2,k=c(q1,q2),bs="temicx",m=c(2,2)),
where q1 and q2 denote the basis dimensions for the marginal smooths.
## S3 method for class 'temicx.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'temicx.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the values of the elements of |
knots |
An optional list containing the knots corresponding to |
An object of class "temicx.smooth". In addition to the usual
elements of a smooth class documented under smooth.construct of the mgcv library,
this object contains:
p.ident |
A vector of 0's and 1's for model parameter identification: 1's indicate parameters which will be exponentiated, 0's - otherwise. |
Zc |
A matrix of identifiability constraints. |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
smooth.construct.tedmd.smooth.spec
## Not run: ## tensor product `temicx' example ## simulating data... set.seed(1) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- exp(4*x1[i])/(1+exp(4*x1[i])) + 2*x2[j]^2 f <- as.vector(t(f1)) y <- f+rnorm(length(f))*0.1 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="temicx",m=2), data=dat) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b,theta = 30, phi = 40) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- exp(4*x1[i])/(1+exp(4*x1[i])) + 2*x2[j]^2 persp(x1,x2,f1,theta = 30, phi = 40) ## End(Not run)## Not run: ## tensor product `temicx' example ## simulating data... set.seed(1) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- exp(4*x1[i])/(1+exp(4*x1[i])) + 2*x2[j]^2 f <- as.vector(t(f1)) y <- f+rnorm(length(f))*0.1 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="temicx",m=2), data=dat) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b,theta = 30, phi = 40) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- exp(4*x1[i])/(1+exp(4*x1[i])) + 2*x2[j]^2 persp(x1,x2,f1,theta = 30, phi = 40) ## End(Not run)
This is a special method function
for creating tensor product bivariate smooths concave in the second covariate which is built by the mgcv constructor function for smooth terms, smooth.construct.
It is constructed from a pair of single penalty
marginal smooths. This tensor product is specified by model terms such as s(x1,x2,k=c(q1,q2),bs="tescv",m=c(2,2)),
where the basis for the first marginal smooth is specified in the second element of bs.
## S3 method for class 'tescv.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'tescv.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the values of the elements of |
knots |
An optional list containing the knots corresponding to |
An object of class "tescv.smooth". In addition to the usual
elements of a smooth class documented under smooth.construct of the mgcv library, this object contains:
p.ident |
A vector of 0's and 1's for model parameter identification: 1's indicate parameters which will be exponentiated, 0's - otherwise. |
Zc |
A matrix of identifiability constraints. |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
smooth.construct.temicv.smooth.spec
smooth.construct.temicx.smooth.spec
smooth.construct.tedecv.smooth.spec
smooth.construct.tedecx.smooth.spec
smooth.construct.tescx.smooth.spec
## Not run: ## tensor product `tescv' example ## simulating data... set.seed(5) n <- 30 x1 <- sort(runif(n)) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- sin(2*x1[i]) - 4*x2[j]^2 f1 <- as.vector(t(f1)) f <- (f1-min(f1))/(max(f1)-min(f1)) y <- f+rnorm(length(f))*0.1 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tescv",m=2), family=gaussian(), data=dat) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE, theta = 50, phi = 20) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b, theta = 50, phi = 20) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- sin(2*x1[i]) - 4*x2[j]^2 persp(x1,x2,f1,theta = 50, phi = 20) ## End(Not run)## Not run: ## tensor product `tescv' example ## simulating data... set.seed(5) n <- 30 x1 <- sort(runif(n)) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- sin(2*x1[i]) - 4*x2[j]^2 f1 <- as.vector(t(f1)) f <- (f1-min(f1))/(max(f1)-min(f1)) y <- f+rnorm(length(f))*0.1 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tescv",m=2), family=gaussian(), data=dat) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE, theta = 50, phi = 20) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b, theta = 50, phi = 20) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- sin(2*x1[i]) - 4*x2[j]^2 persp(x1,x2,f1,theta = 50, phi = 20) ## End(Not run)
This is a special method function
for creating tensor product bivariate smooths convex in the second covariate which is built by the mgcv constructor function for smooth terms, smooth.construct.
It is constructed from a pair of single penalty
marginal smooths. This tensor product is specified by model terms such as s(x1,x2,k=c(q1,q2),bs="tescx",m=c(2,2)),
where the basis for the first marginal smooth is specified in the second element of bs.
## S3 method for class 'tescx.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'tescx.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the values of the elements of |
knots |
An optional list containing the knots corresponding to |
An object of class "tescx.smooth". In addition to the usual
elements of a smooth class documented under smooth.construct of the mgcv library, this object contains:
p.ident |
A vector of 0's and 1's for model parameter identification: 1's indicate parameters which will be exponentiated, 0's - otherwise. |
Zc |
A matrix of identifiability constraints. |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
smooth.construct.temicv.smooth.spec
smooth.construct.temicx.smooth.spec
smooth.construct.tedecv.smooth.spec
smooth.construct.tedecx.smooth.spec
smooth.construct.tescv.smooth.spec
## Not run: ## tensor product `tescx' example ## simulating data... set.seed(2) n <- 30 x1 <- sort(runif(n)) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- sin(x1[i]) + 2*x2[j]^2 f1 <- as.vector(t(f1)) f <- (f1-min(f1))/(max(f1)-min(f1)) y <- f+rnorm(length(f))*0.1 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tescx",m=2), family=gaussian(), data=dat) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE, theta = 50, phi = 20) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b, theta = 50, phi = 20) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- sin(x1[i]) + 2*x2[j]^2 persp(x1,x2,f1,theta = 50, phi = 20) ## End(Not run)## Not run: ## tensor product `tescx' example ## simulating data... set.seed(2) n <- 30 x1 <- sort(runif(n)) x2 <- sort(2*runif(n)-1) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- sin(x1[i]) + 2*x2[j]^2 f1 <- as.vector(t(f1)) f <- (f1-min(f1))/(max(f1)-min(f1)) y <- f+rnorm(length(f))*0.1 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,k=c(10,10),bs="tescx",m=2), family=gaussian(), data=dat) ## plot results ... par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE, theta = 50, phi = 20) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") x11() vis.scam(b, theta = 50, phi = 20) ## plotting the truth... x11() x1 <- seq(min(x1),max(x1),length.out=30) x2 <- seq(min(x2),max(x2),length.out=30) f1 <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f1[i,j] <- sin(x1[i]) + 2*x2[j]^2 persp(x1,x2,f1,theta = 50, phi = 20) ## End(Not run)
This is a special method function
for creating tensor product bivariate smooths monotone decreasing in the first covariate which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed from a pair of single penalty marginal smooths.
This tensor product is specified by model terms such as s(x1,x2,k=c(q1,q2),bs="tesmd1",m=2).
The default basis for the second marginal smooth is P-spline. Cyclic cubic regression spline ("cc") is implemented in addition to the P-spline. See an example below on how to call for it.
## S3 method for class 'tesmd1.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'tesmd1.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the values of the elements of |
knots |
An optional list containing the knots corresponding to |
An object of class "tesmd1.smooth". In addition to the usual
elements of a smooth class documented under smooth.construct of the mgcv library,
this object contains:
p.ident |
A vector of 0's and 1's for model parameter identification: 1's indicate parameters which will be exponentiated, 0's - otherwise. |
Zc |
A matrix of identifiability constraints. |
margin.bs |
A two letter character string indicating the (penalized) smoothing basis to use for the second unconstrained marginal smooth. (eg |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
smooth.construct.tesmd2.smooth.spec
## Not run: ## tensor product `tesmd1' example ## simulating data... require(scam) set.seed(2) n <- 30 x1 <- sort(runif(n)*4-1); x2 <- sort(runif(n)) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i]))+2*sin(pi*x2[j]) f <- as.vector(t(f)) y <- f+rnorm(length(f))*.2 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,bs="tesmd1",k=10),data=dat) ## plot results ... old.par <- par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") par(old.par) vis.scam(b,theta=40,phi=20) ## example with cyclic cubic regression spline along the second covariate... set.seed(2) n <- 30 x1 <- sort(runif(n)*4-1); x2 <- sort(runif(n)) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i]))+sin(2*pi*x2[j]) f <- as.vector(t(f)) y <- f+rnorm(length(f))*.2 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b1 <- scam(y~s(x1,x2,bs="tesmd1",xt=list("cc"),k=10), data=dat) ## plot results ... old.par <- par(mfrow=c(2,2)) plot(b1,se=TRUE) plot(b1,pers=TRUE,theta = 30, phi = 40) plot(y,b1$fitted.values,xlab="Simulated data",ylab="Fitted data") par(old.par) vis.scam(b1,theta=40,phi=20) ## End(Not run)## Not run: ## tensor product `tesmd1' example ## simulating data... require(scam) set.seed(2) n <- 30 x1 <- sort(runif(n)*4-1); x2 <- sort(runif(n)) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i]))+2*sin(pi*x2[j]) f <- as.vector(t(f)) y <- f+rnorm(length(f))*.2 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,bs="tesmd1",k=10),data=dat) ## plot results ... old.par <- par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") par(old.par) vis.scam(b,theta=40,phi=20) ## example with cyclic cubic regression spline along the second covariate... set.seed(2) n <- 30 x1 <- sort(runif(n)*4-1); x2 <- sort(runif(n)) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i]))+sin(2*pi*x2[j]) f <- as.vector(t(f)) y <- f+rnorm(length(f))*.2 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b1 <- scam(y~s(x1,x2,bs="tesmd1",xt=list("cc"),k=10), data=dat) ## plot results ... old.par <- par(mfrow=c(2,2)) plot(b1,se=TRUE) plot(b1,pers=TRUE,theta = 30, phi = 40) plot(y,b1$fitted.values,xlab="Simulated data",ylab="Fitted data") par(old.par) vis.scam(b1,theta=40,phi=20) ## End(Not run)
This is a special method function
for creating tensor product bivariate smooths monotone decreasing in the second covariate which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed from a pair of single penalty
marginal smooths. This tensor product is specified by model terms such as s(x1,x2,k=c(q1,q2),bs="tesmd2",m=c(2,2)). The default basis for the first marginal smooth is P-spline. Cyclic cubic regression spline ("cc") is implemented in addition to the P-spline. See an example below on how to call for it.
## S3 method for class 'tesmd2.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'tesmd2.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the values of the elements of |
knots |
An optional list containing the knots corresponding to |
An object of class "tesmd2.smooth". In addition to the usual
elements of a smooth class documented under smooth.construct of the mgcv library,
this object contains:
p.ident |
A vector of 0's and 1's for model parameter identification: 1's indicate parameters which will be exponentiated, 0's - otherwise. |
Zc |
A matrix of identifiability constraints. |
margin.bs |
A two letter character string indicating the (penalized) smoothing basis to use for the first unconstrained marginal smooth. (eg |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
smooth.construct.tesmd1.smooth.spec
## Not run: ## tensor product `tesmd2' example ## simulating data... require(scam) set.seed(2) n <- 30 x1 <- sort(runif(n)); x2 <- sort(runif(n)*4-1) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- 2*sin(pi*x1[i])-exp(4*x2[j])/(1+exp(4*x2[j])) f <- as.vector(t(f)) y <- f+rnorm(length(f))*.2 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,bs="tesmd2",k=10),data=dat) ## plot results ... old.par <- par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,scheme=1,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") par(old.par) vis.scam(b,theta = 40, phi = 20) ## example with cyclic cubic regression spline along the 1st covariate... set.seed(4) n <- 30 x1 <- sort(runif(n)); x2 <- sort(runif(n)*4-1) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- sin(2*pi*x1[i])-exp(4*x2[j])/(1+exp(4*x2[j])) f <- as.vector(t(f)) y <- f+rnorm(length(f))*.2 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b1 <- scam(y~s(x1,x2,bs="tesmd2",xt=list("cc"),k=10), data=dat) ## plot results ... old.par <-par(mfrow=c(2,2)) plot(b1,se=TRUE) plot(b1,scheme=1,theta = 30, phi = 40) plot(y,b1$fitted.values,xlab="Simulated data",ylab="Fitted data") par(old.par) vis.scam(b1,theta=40,phi=20) ## End(Not run)## Not run: ## tensor product `tesmd2' example ## simulating data... require(scam) set.seed(2) n <- 30 x1 <- sort(runif(n)); x2 <- sort(runif(n)*4-1) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- 2*sin(pi*x1[i])-exp(4*x2[j])/(1+exp(4*x2[j])) f <- as.vector(t(f)) y <- f+rnorm(length(f))*.2 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,bs="tesmd2",k=10),data=dat) ## plot results ... old.par <- par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,scheme=1,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") par(old.par) vis.scam(b,theta = 40, phi = 20) ## example with cyclic cubic regression spline along the 1st covariate... set.seed(4) n <- 30 x1 <- sort(runif(n)); x2 <- sort(runif(n)*4-1) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- sin(2*pi*x1[i])-exp(4*x2[j])/(1+exp(4*x2[j])) f <- as.vector(t(f)) y <- f+rnorm(length(f))*.2 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b1 <- scam(y~s(x1,x2,bs="tesmd2",xt=list("cc"),k=10), data=dat) ## plot results ... old.par <-par(mfrow=c(2,2)) plot(b1,se=TRUE) plot(b1,scheme=1,theta = 30, phi = 40) plot(y,b1$fitted.values,xlab="Simulated data",ylab="Fitted data") par(old.par) vis.scam(b1,theta=40,phi=20) ## End(Not run)
This is a special method function
for creating tensor product bivariate smooths monotone increasing in the first covariate which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed from a pair of single penalty
marginal smooths. This tensor product is specified by model terms such as s(x1,x2,k=c(q1,q2),bs="tesmi1",m=2).
The basis for the second marginal smooth can be specified as a two letter character string of the argument xt (eg xt="cc" to specify cyclic cubic regression spline). See example below. The default basis for the second marginal smooth is P-spline. Cyclic cubic regression spline ("cc") is implemented in addition to the P-spline. See an example below on how to call for it.
## S3 method for class 'tesmi1.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'tesmi1.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the values of the elements of |
knots |
An optional list containing the knots corresponding to |
An object of class "tesmi1.smooth". In addition to the usual
elements of a smooth class documented under smooth.construct of the mgcv library,
this object contains:
p.ident |
A vector of 0's and 1's for model parameter identification: 1's indicate parameters which will be exponentiated, 0's - otherwise. |
Zc |
A matrix of identifiability constraints. |
margin.bs |
A two letter character string indicating the (penalized) smoothing basis to use for the second unconstrained marginal smooth. (eg |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
smooth.construct.tesmi2.smooth.spec
## Not run: ## tensor product `tesmi1' example... ## simulating data... require(scam) set.seed(2) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(runif(n)) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- exp(4*x1[i])/(1+exp(4*x1[i]))+2*sin(pi*x2[j]) f <- as.vector(t(f)) y <- f+rnorm(length(f))*.3 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,bs="tesmi1",k=c(10,10)), data=dat) ## plot results ... old.par<- par(mfrow=c(2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") par(old.par) vis.scam(b,theta=40,phi=20) ## example with cyclic cubic regression spline along the second covariate... set.seed(2) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(runif(n)) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- exp(4*x1[i])/(1+exp(4*x1[i]))+sin(2*pi*x2[j]) f <- as.vector(t(f)) y <- f+rnorm(length(f))*.2 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b1 <- scam(y~s(x1,x2,bs="tesmi1",xt=list("cc"),k=10), data=dat) ## plot results ... old.par<- par(mfrow=c(2,2)) plot(b1,se=TRUE) plot(b1,pers=TRUE,theta = 30, phi = 40) plot(y,b1$fitted.values,xlab="Simulated data",ylab="Fitted data") par(old.par) vis.scam(b1,theta=40,phi=20) ## End(Not run)## Not run: ## tensor product `tesmi1' example... ## simulating data... require(scam) set.seed(2) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(runif(n)) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- exp(4*x1[i])/(1+exp(4*x1[i]))+2*sin(pi*x2[j]) f <- as.vector(t(f)) y <- f+rnorm(length(f))*.3 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,bs="tesmi1",k=c(10,10)), data=dat) ## plot results ... old.par<- par(mfrow=c(2,2)) plot(b,se=TRUE) plot(b,pers=TRUE,theta = 30, phi = 40) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") par(old.par) vis.scam(b,theta=40,phi=20) ## example with cyclic cubic regression spline along the second covariate... set.seed(2) n <- 30 x1 <- sort(runif(n)*4-1) x2 <- sort(runif(n)) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- exp(4*x1[i])/(1+exp(4*x1[i]))+sin(2*pi*x2[j]) f <- as.vector(t(f)) y <- f+rnorm(length(f))*.2 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b1 <- scam(y~s(x1,x2,bs="tesmi1",xt=list("cc"),k=10), data=dat) ## plot results ... old.par<- par(mfrow=c(2,2)) plot(b1,se=TRUE) plot(b1,pers=TRUE,theta = 30, phi = 40) plot(y,b1$fitted.values,xlab="Simulated data",ylab="Fitted data") par(old.par) vis.scam(b1,theta=40,phi=20) ## End(Not run)
This is a special method function
for creating tensor product bivariate smooths monotone increasing in the second covariate which is built by
the mgcv constructor function for smooth terms, smooth.construct.
It is constructed from a pair of single penalty
marginal smooths. This tensor product is specified by model terms such as s(x1,x2,k=c(q1,q2),bs="tesmi2",m=c(2,2)). The default basis for the first marginal smooth is P-spline. Cyclic cubic regression spline ("cc") is implemented in addition to the P-spline. See an example below on how to call for it.
## S3 method for class 'tesmi2.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'tesmi2.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the values of the elements of |
knots |
An optional list containing the knots corresponding to |
An object of class "tesmi2.smooth". In addition to the usual
elements of a smooth class documented under smooth.construct of the mgcv library,
this object contains:
p.ident |
A vector of 0's and 1's for model parameter identification: 1's indicate parameters which will be exponentiated, 0's - otherwise. |
Zc |
A matrix of identifiability constraints. |
margin.bs |
A two letter character string indicating the (penalized) smoothing basis to use for the first unconstrained marginal smooth. (eg |
Natalya Pya <[email protected]>
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
smooth.construct.tesmi1.smooth.spec
## Not run: ## tensor product `tesmi2' example ## simulating data... require(scam) set.seed(2) n <- 30 x1 <- sort(runif(n)); x2 <- sort(runif(n)*4-1) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- 2*sin(pi*x1[i]) +exp(4*x2[j])/(1+exp(4*x2[j])) f <- as.vector(t(f)) y <- f+rnorm(length(f))*.3 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,bs="tesmi2",k=c(10,10)),data=dat) ## plot results ... old.par <- par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE, theta = 50, phi = 20) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") par(old.par) vis.scam(b,theta=50,phi=20) ## example with cyclic cubic regression spline along the 1st covariate... set.seed(2) n <- 30 x1 <- sort(runif(n)); x2 <- sort(runif(n)*4-1) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- sin(2*pi*x1[i])+ exp(4*x2[j])/(1+exp(4*x2[j])) f <- as.vector(t(f)) y <- f+rnorm(length(f))*.3 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b1 <- scam(y~s(x1,x2,bs="tesmi2",xt=list("cc"),k=10), data=dat) ## plot results ... old.par <- par(mfrow=c(2,2)) plot(b1,se=TRUE) plot(b1,pers=TRUE,theta = 50, phi = 20) plot(y,b1$fitted.values,xlab="Simulated data",ylab="Fitted data") par(old.par) vis.scam(b1,theta=40,phi=20) ## End(Not run)## Not run: ## tensor product `tesmi2' example ## simulating data... require(scam) set.seed(2) n <- 30 x1 <- sort(runif(n)); x2 <- sort(runif(n)*4-1) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- 2*sin(pi*x1[i]) +exp(4*x2[j])/(1+exp(4*x2[j])) f <- as.vector(t(f)) y <- f+rnorm(length(f))*.3 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b <- scam(y~s(x1,x2,bs="tesmi2",k=c(10,10)),data=dat) ## plot results ... old.par <- par(mfrow=c(2,2),mar=c(4,4,2,2)) plot(b,se=TRUE) plot(b,pers=TRUE, theta = 50, phi = 20) plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data") par(old.par) vis.scam(b,theta=50,phi=20) ## example with cyclic cubic regression spline along the 1st covariate... set.seed(2) n <- 30 x1 <- sort(runif(n)); x2 <- sort(runif(n)*4-1) f <- matrix(0,n,n) for (i in 1:n) for (j in 1:n) f[i,j] <- sin(2*pi*x1[i])+ exp(4*x2[j])/(1+exp(4*x2[j])) f <- as.vector(t(f)) y <- f+rnorm(length(f))*.3 x11 <- matrix(0,n,n) x11[,1:n] <- x1 x11 <- as.vector(t(x11)) x22 <- rep(x2,n) dat <- list(x1=x11,x2=x22,y=y) ## fit model ... b1 <- scam(y~s(x1,x2,bs="tesmi2",xt=list("cc"),k=10), data=dat) ## plot results ... old.par <- par(mfrow=c(2,2)) plot(b1,se=TRUE) plot(b1,pers=TRUE,theta = 50, phi = 20) plot(y,b1$fitted.values,xlab="Simulated data",ylab="Fitted data") par(old.par) vis.scam(b1,theta=40,phi=20) ## End(Not run)
This is a special method function
for creating tensor product bivariate interaction smooths with decreasing constraint in the first covariate,
appropriate when the main effects (and any lower interactions) are also present.
It is built by the mgcv constructor function for smooth terms, smooth.construct, and constructed from a pair of single penalty marginal smooths. This tensor product is specified by model terms such as s(x1,x2,k=c(q1,q2),bs="tismd").
See example below.
## S3 method for class 'tismd.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'tismd.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the values of the elements of |
knots |
An optional list containing the knots corresponding to |
In some cases, it is helpful to consider models with a main-effects + interactions structure, for example,
where and are smooth ‘main effects’ and is a smooth ‘interaction’ subject to decreasing constraint wrt ( can be subject to decreasing constraint).
Constructing such functional ANOVA decomposition recognises the fact that the tensor produc basis construction is exactly the same as the construction used for any interaction in a linear model. tismd produce tensor product interactions with decreasing constraint along the first covariate from which the main effects have been excluded, under the assumption that they will be included separately. For example, the ~ s(x) + s(z) + s(x,z,bs="tismd") would produce the above main effects + interaction structure. Specifically, the marginal smooths of a tensor product, tismd, are subject to identifiability constraints before constructing the tensor product basis. This results in the interaction smooths that do not include the corresponding main effects. tismd apply SCOP identifiability constraints to the first marginal and sum-to-zero constraints to the second unconstrained marginal. See Wood (2017, section 5.6.3) for ANOVA decompositions of unconstrained smooths.
An object of class "tismd.smooth". In addition to the usual
elements of a smooth class documented under smooth.construct of the mgcv library,
this object contains:
p.ident |
A vector of 0's and 1's for model parameter identification: 1's indicate parameters which will be exponentiated, 0's - otherwise. |
Zc |
A matrix of identifiability constraints. |
margin.bs |
A two letter character string indicating the (penalized) smoothing basis to use for the second unconstrained marginal smooth. (eg |
Natalya Pya [email protected]
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Wood S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapman and Hall/CRC Press.
smooth.construct.tesmd1.smooth.spec
smooth.construct.tismi.smooth.spec
## Not run: ## tensor product `tismd' example... ## simulating data... require(scam) test <- function(x,z){ -exp(4*x)/(1+exp(4*x))-2*sin(pi*z)-(x+1)^0.6*z } set.seed(7) n <- 600 x <- runif(n)*4-1 z <- runif(n) xs <- seq(-1,3,length=30); zs <- seq(0,1,length=30) pr <- data.frame(x=rep(xs,30),z=rep(zs,rep(30,30))) truth <- matrix(test(pr$x,pr$z),30,30) f <- test(x,z) y <- f + rnorm(n)*0.3 bi <- scam(y~ ti(x)+ti(z)+ s(x,z,bs="tismd")) summary(bi) old.par <- par(mfrow=c(2,2)) persp(xs,zs,truth);title("truth") vis.scam(bi);title("tismd") ## fitting with "tesmd1" instead... bc <- scam(y~s(x,z,bs="tesmd1")) vis.scam(bc);title("tesmd1") par(old.par) plot(bi,pages=1,scheme=2) plot(bi,select=3,scheme=1,zlim=c(-3,3)) ## End(Not run)## Not run: ## tensor product `tismd' example... ## simulating data... require(scam) test <- function(x,z){ -exp(4*x)/(1+exp(4*x))-2*sin(pi*z)-(x+1)^0.6*z } set.seed(7) n <- 600 x <- runif(n)*4-1 z <- runif(n) xs <- seq(-1,3,length=30); zs <- seq(0,1,length=30) pr <- data.frame(x=rep(xs,30),z=rep(zs,rep(30,30))) truth <- matrix(test(pr$x,pr$z),30,30) f <- test(x,z) y <- f + rnorm(n)*0.3 bi <- scam(y~ ti(x)+ti(z)+ s(x,z,bs="tismd")) summary(bi) old.par <- par(mfrow=c(2,2)) persp(xs,zs,truth);title("truth") vis.scam(bi);title("tismd") ## fitting with "tesmd1" instead... bc <- scam(y~s(x,z,bs="tesmd1")) vis.scam(bc);title("tesmd1") par(old.par) plot(bi,pages=1,scheme=2) plot(bi,select=3,scheme=1,zlim=c(-3,3)) ## End(Not run)
This is a special method function
for creating tensor product bivariate interaction smooths with increasing constraint in the first covariate,
appropriate when the main effects (and any lower interactions) are also present.
It is built by the mgcv constructor function for smooth terms, smooth.construct, and constructed from a pair of single penalty marginal smooths. This tensor product is specified by model terms such as s(x1,x2,k=c(q1,q2),bs="tismi").
See example below.
## S3 method for class 'tismi.smooth.spec' smooth.construct(object, data, knots)## S3 method for class 'tismi.smooth.spec' smooth.construct(object, data, knots)
object |
A smooth specification object, generated by an |
data |
A data frame or list containing the values of the elements of |
knots |
An optional list containing the knots corresponding to |
In some cases, it is helpful to consider models with a main-effects + interactions structure, for example,
where and are smooth ‘main effects’ and is a smooth ‘interaction’ subject to increasing constraint wrt ( can be subject to increasing constraint).
Constructing such functional ANOVA decomposition recognises the fact that the tensor produc basis construction is exactly the same as the construction used for any interaction in a linear model. tismi produce tensor product interactions with increasing constraint along the first covariate from which the main effects have been excluded, under the assumption that they will be included separately. For example, the ~ s(x) + s(z) + s(x,z,bs="tismi") would produce the above main effects + interaction structure. Specifically, the marginal smooths of a tensor product, tismi, are subject to identifiability constraints before constructing the tensor product basis. This results in the interaction smooths that do not include the corresponding main effects. tismi apply SCOP identifiability constraints to the first marginal and sum-to-zero constraints to the second unconstrained marginal. See Wood (2017, section 5.6.3) for ANOVA decompositions of unconstrained smooths.
An object of class "tismi.smooth". In addition to the usual
elements of a smooth class documented under smooth.construct of the mgcv library,
this object contains:
p.ident |
A vector of 0's and 1's for model parameter identification: 1's indicate parameters which will be exponentiated, 0's - otherwise. |
Zc |
A matrix of identifiability constraints. |
margin.bs |
A two letter character string indicating the (penalized) smoothing basis to use for the second unconstrained marginal smooth. (eg |
Natalya Pya [email protected]
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Wood S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapman and Hall/CRC Press
smooth.construct.tismd.smooth.spec
smooth.construct.tesmi1.smooth.spec
## Not run: ## tensor product `tismi' example... require(scam) test <- function(x,z){ exp(4*x)/(1+exp(4*x))+2*sin(pi*z)+(x+1)^0.6*z } set.seed(7) n <- 600 x <- runif(n)*4-1 z <- runif(n) xs <- seq(-1,3,length=30); zs <- seq(0,1,length=30) pr <- data.frame(x=rep(xs,30),z=rep(zs,rep(30,30))) truth <- matrix(test(pr$x,pr$z),30,30) f <- test(x,z) y <- f + rnorm(n)*0.3 bi <- scam(y~ ti(x)+ti(z)+ s(x,z,bs="tismi")) summary(bi) old.par <- par(mfrow=c(2,2)) persp(xs,zs,truth);title("truth") vis.scam(bi);title("tismi") ## fitting with "tesmi1"... bc <- scam(y~s(x,z,bs="tesmi1")) vis.scam(bc);title("tesmi1") par(old.par) plot(bi,pages=1,scheme=2) plot(bi,select=3,scheme=1,zlim=c(-5,5)) ## End(Not run)## Not run: ## tensor product `tismi' example... require(scam) test <- function(x,z){ exp(4*x)/(1+exp(4*x))+2*sin(pi*z)+(x+1)^0.6*z } set.seed(7) n <- 600 x <- runif(n)*4-1 z <- runif(n) xs <- seq(-1,3,length=30); zs <- seq(0,1,length=30) pr <- data.frame(x=rep(xs,30),z=rep(zs,rep(30,30))) truth <- matrix(test(pr$x,pr$z),30,30) f <- test(x,z) y <- f + rnorm(n)*0.3 bi <- scam(y~ ti(x)+ti(z)+ s(x,z,bs="tismi")) summary(bi) old.par <- par(mfrow=c(2,2)) persp(xs,zs,truth);title("truth") vis.scam(bi);title("tismi") ## fitting with "tesmi1"... bc <- scam(y~s(x,z,bs="tesmi1")) vis.scam(bc);title("tesmi1") par(old.par) plot(bi,pages=1,scheme=2) plot(bi,select=3,scheme=1,zlim=c(-5,5)) ## End(Not run)
Takes a fitted scam object produced by scam() and produces various useful
summaries from it. The same code as in summary.gam of the mgcv package is used here with slight modifications
to accept the exponentiated parameters of the shape constrained smooth terms and the corresponding covariance matrix.
## S3 method for class 'scam' summary(object, dispersion=NULL, freq=FALSE, ...) ## S3 method for class 'summary.scam' print(x,digits = max(3, getOption("digits") - 3), signif.stars = getOption("show.signif.stars"),...)## S3 method for class 'scam' summary(object, dispersion=NULL, freq=FALSE, ...) ## S3 method for class 'summary.scam' print(x,digits = max(3, getOption("digits") - 3), signif.stars = getOption("show.signif.stars"),...)
object |
a fitted |
x |
a |
dispersion |
A known dispersion parameter. |
freq |
By default p-values for individual terms are calculated using the Bayesian estimated covariance matrix of the parameter estimators. If this is set to TRUE then the frequentist covariance matrix of the parameters is used instead. |
digits |
controls number of digits printed in output. |
signif.stars |
Should significance stars be printed alongside output. |
... |
other arguments. |
summary.scam produces the same list of summary information for a fitted scam object as in the
unconstrained case summary.gam except for the last element BFGS termination condition.
p.coeff |
is an array of estimates of the strictly parametric model coefficients. |
p.t |
is an array of the |
p.pv |
is an array of p-values for the null hypothesis that the corresponding parameter is zero. Calculated with reference to the t distribution with the estimated residual degrees of freedom for the model fit if the dispersion parameter has been estimated, and the standard normal if not. |
m |
The number of smooth terms in the model. |
chi.sq |
An array of test statistics for assessing the significance of model smooth terms. See details. |
s.pv |
An array of approximate p-values for the null hypotheses that each smooth term is zero. Be warned, these are only approximate. |
se |
array of standard error estimates for all parameter estimates. |
r.sq |
The adjusted r-squared for the model. Defined as the proportion of variance explained, where original variance and residual variance are both estimated using unbiased estimators. This quantity can be negative if your model is worse than a one parameter constant model, and can be higher for the smaller of two nested models! Note that proportion null deviance explained is probably more appropriate for non-normal errors. |
dev.expl |
The proportion of the null deviance explained by the model. |
edf |
array of estimated degrees of freedom for the model terms. |
residual.df |
estimated residual degrees of freedom. |
n |
number of data. |
gcv |
minimized GCV score for the model, if GCV used. |
ubre |
minimized UBRE score for the model, if UBRE used. |
scale |
estimated (or given) scale parameter. |
family |
the family used. |
formula |
the original scam formula. |
dispersion |
the scale parameter. |
pTerms.df |
the degrees of freedom associated with each parameteric term (excluding the constant). |
pTerms.chi.sq |
a Wald statistic for testing the null hypothesis that the each parametric term is zero. |
pTerms.pv |
p-values associated with the tests that each term is zero. For penalized fits these are approximate. The reference distribution is an appropriate chi-squared when the scale parameter is known, and is based on an F when it is not. |
cov.unscaled |
The estimated covariance matrix of the parameters (or
estimators if |
cov.scaled |
The estimated covariance matrix of the parameters
(estimators if |
p.table |
significance table for parameters |
s.table |
significance table for smooths |
pTerms.table |
significance table for parametric model terms |
BFGS termination condition |
the value of the maximum component of the scaled GCV/UBRE gradient used as stopping
condition. This value is printed if
the termination code of the BFGS optimization process is not ‘1’ (not full convergence)
(see |
The p-values are approximate.
Natalya Pya <[email protected]> based on mgcv by Simon Wood
Wood S.N. (2006) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press.
Pya, N. and Wood, S.N. (2015) Shape constrained additive models. Statistics and Computing, 25(3), 543-559
Pya, N. (2010) Additive models with shape constraints. PhD thesis. University of Bath. Department of Mathematical Sciences
## Not run: ## simulating data... require(scam) n <- 200 set.seed(1) x1 <- runif(n)*6-3 f1 <- 3*exp(-x1^2) # unconstrained smooth term x2 <- runif(n)*4-1; f2 <- exp(4*x2)/(1+exp(4*x2)) # monotone increasing smooth x3 <- runif(n)*5; f3 <- -log(x3)/5 # monotone decreasing smooth f <- f1+f2+f3 y <- f + rnorm(n)*.3 dat <- data.frame(x1=x1,x2=x2,x3=x3,y=y) ## fit model ... b <- scam(y~s(x1,k=15,bs="cr",m=2)+s(x2,k=30,bs="mpi",m=2)+s(x3,k=30,bs="mpd",m=2), data=dat) summary(b) plot(b,pages=1,shade=TRUE) ## End(Not run)## Not run: ## simulating data... require(scam) n <- 200 set.seed(1) x1 <- runif(n)*6-3 f1 <- 3*exp(-x1^2) # unconstrained smooth term x2 <- runif(n)*4-1; f2 <- exp(4*x2)/(1+exp(4*x2)) # monotone increasing smooth x3 <- runif(n)*5; f3 <- -log(x3)/5 # monotone decreasing smooth f <- f1+f2+f3 y <- f + rnorm(n)*.3 dat <- data.frame(x1=x1,x2=x2,x3=x3,y=y) ## fit model ... b <- scam(y~s(x1,k=15,bs="cr",m=2)+s(x2,k=30,bs="mpi",m=2)+s(x3,k=30,bs="mpd",m=2), data=dat) summary(b) plot(b,pages=1,shade=TRUE) ## End(Not run)
Produces perspective or contour plot views of scam model
predictions.
The code is a clone of vis.gam of the mgcv package.
vis.scam(x,view=NULL,cond=list(),n.grid=30,too.far=0,col=NA, color="heat",contour.col=NULL,se=-1,type="link", plot.type="persp",zlim=NULL,nCol=50,...)vis.scam(x,view=NULL,cond=list(),n.grid=30,too.far=0,col=NA, color="heat",contour.col=NULL,se=-1,type="link", plot.type="persp",zlim=NULL,nCol=50,...)
The documentation below is the same as in documentation object vis.gam.
x |
a |
view |
an array containing the names of the two main effect terms to be displayed on the x and y dimensions of the plot. If omitted the first two suitable terms will be used. |
cond |
a named list of the values to use for the other predictor terms
(not in |
n.grid |
The number of grid nodes in each direction used for calculating the plotted surface. |
too.far |
plot grid nodes that are too far from the points defined by the variables given in |
col |
The colours for the facets of the plot. If this is |
color |
the colour scheme to use for plots when |
contour.col |
sets the colour of contours when using |
se |
if less than or equal to zero then only the predicted surface is plotted, but if greater than zero, then 3
surfaces are plotted, one at the predicted values minus |
type |
|
plot.type |
one of |
zlim |
a two item array giving the lower and upper limits for the z-axis
scale. |
nCol |
The number of colors to use in color schemes. |
... |
Simply produces a plot.
Simon Wood [email protected]
library(scam) # Example with factor variable set.seed(0) fac<-rep(1:4,20) x <- runif(80)*5; y <- fac+log(x)/5+rnorm(80)*0.1 fac <- factor(fac) b <- scam(y~fac+s(x,bs="mpi")) vis.scam(b,theta=-35,color="heat") # factor example # Example with "by" variable z<-rnorm(80)*0.4 y<-as.numeric(fac)+log(x)*z+rnorm(80)*0.1 b<-scam(y~fac+s(x,by=z)) g <- gam(y~fac+s(x,by=z)) vis.scam(b,theta=-35,color="terrain",cond=list(z=1)) # by variable example vis.scam(b,view=c("z","x"),theta= 65) # plot against by variable ## compare with gam(mgcv)... vis.gam(g,theta=-35,color="terrain",cond=list(z=1)) # by variable example vis.gam(g,view=c("z","x"),theta= 65) # plot against by variable ## all three smooths are increasing... set.seed(2) n <- 400 x <- runif(n, 0, 1) f1 <- log(x *5) f2 <- exp(2 * x) - 4 f3 <- 5* sin(x) e <- rnorm(n, 0, 2) fac <- as.factor(sample(1:3,n,replace=TRUE)) fac.1 <- as.numeric(fac==1) fac.2 <- as.numeric(fac==2) fac.3 <- as.numeric(fac==3) y <- f1*fac.1 + f2*fac.2 + f3*fac.3 + e dat <- data.frame(y=y,x=x,fac=fac,f1=f1,f2=f2,f3=f3) b1 <- scam(y ~ s(x,by=fac,bs="mpi"),data=dat,optimizer="efs") plot(b1,pages=1,scale=0,shade=TRUE) summary(b1) vis.scam(b1,theta=-40,color="terrain",cond=list(z=1)) ## note that the preceding, b1, fit is the same as.... b2 <- scam(y ~ s(x,by=as.numeric(fac==1),bs="mpi")+s(x,by=as.numeric(fac==2),bs="mpi")+ s(x,by=as.numeric(fac==3),bs="mpi"),data=dat,optimizer="efs") summary(b2) ## Note that as in gam() when using factor 'by' variables, centering ## constraints are applied to the smooths, which usually means that the 'by' ## variable should be included as a parametric term, as well. ## The difference with scam() is that here a 'zero intercept' constraint is ## applied in place of 'centering' (although scam's fitted smooths are centred for plotting). ## compare with the gam() fits.. g1 <- gam(y ~ fac+s(x,by=fac),data=dat) g2 <- gam(y ~ s(x,by=fac),data=dat) summary(g1) summary(g2) plot(g1,pages=1,scale=0,shade=TRUE)library(scam) # Example with factor variable set.seed(0) fac<-rep(1:4,20) x <- runif(80)*5; y <- fac+log(x)/5+rnorm(80)*0.1 fac <- factor(fac) b <- scam(y~fac+s(x,bs="mpi")) vis.scam(b,theta=-35,color="heat") # factor example # Example with "by" variable z<-rnorm(80)*0.4 y<-as.numeric(fac)+log(x)*z+rnorm(80)*0.1 b<-scam(y~fac+s(x,by=z)) g <- gam(y~fac+s(x,by=z)) vis.scam(b,theta=-35,color="terrain",cond=list(z=1)) # by variable example vis.scam(b,view=c("z","x"),theta= 65) # plot against by variable ## compare with gam(mgcv)... vis.gam(g,theta=-35,color="terrain",cond=list(z=1)) # by variable example vis.gam(g,view=c("z","x"),theta= 65) # plot against by variable ## all three smooths are increasing... set.seed(2) n <- 400 x <- runif(n, 0, 1) f1 <- log(x *5) f2 <- exp(2 * x) - 4 f3 <- 5* sin(x) e <- rnorm(n, 0, 2) fac <- as.factor(sample(1:3,n,replace=TRUE)) fac.1 <- as.numeric(fac==1) fac.2 <- as.numeric(fac==2) fac.3 <- as.numeric(fac==3) y <- f1*fac.1 + f2*fac.2 + f3*fac.3 + e dat <- data.frame(y=y,x=x,fac=fac,f1=f1,f2=f2,f3=f3) b1 <- scam(y ~ s(x,by=fac,bs="mpi"),data=dat,optimizer="efs") plot(b1,pages=1,scale=0,shade=TRUE) summary(b1) vis.scam(b1,theta=-40,color="terrain",cond=list(z=1)) ## note that the preceding, b1, fit is the same as.... b2 <- scam(y ~ s(x,by=as.numeric(fac==1),bs="mpi")+s(x,by=as.numeric(fac==2),bs="mpi")+ s(x,by=as.numeric(fac==3),bs="mpi"),data=dat,optimizer="efs") summary(b2) ## Note that as in gam() when using factor 'by' variables, centering ## constraints are applied to the smooths, which usually means that the 'by' ## variable should be included as a parametric term, as well. ## The difference with scam() is that here a 'zero intercept' constraint is ## applied in place of 'centering' (although scam's fitted smooths are centred for plotting). ## compare with the gam() fits.. g1 <- gam(y ~ fac+s(x,by=fac),data=dat) g2 <- gam(y ~ s(x,by=fac),data=dat) summary(g1) summary(g2) plot(g1,pages=1,scale=0,shade=TRUE)