# Fits Smoothing Splines with Parametric Effects Given a real-valued response vector $\mathbf{y}=\{y_{i}\}_{n\times1}$, a semiparametric regression model has the form [REMOVE_ME]$$ y_{i}= \eta(\mathbf{x}_{i}) + \sum_{j=1}^{t}b_{j}z_{ij} + e_{i} [REMOVE_ME_2]$$ where $y_{i}$ is the $i$-th observation's respone, $\mathbf{x}_{i}=(x_{i1},\ldots,x_{ip})$ is the $i$-th observation's nonparametric predictor vector, $\eta$ is an unknown smooth function relating the response and nonparametric predictors, $\mathbf{z}_{i}=(z_{i1},\ldots,z_{it})$ is the $i$-th observation's parametric predictor vector, and $e_{i}\sim\mathrm{N}(0,\sigma^{2})$ is iid Gaussian error. Function can fit both additive and interactive non/parametric effects, and allows for 2-way and 3-way interactions between nonparametric and parametric effects (see Details and Examples). ## Description Given a real-valued response vector $\mathbf{y}=\{y_{i}\}_{n\times1}$, a semiparametric regression model has the form $$ y_{i}= \eta(\mathbf{x}_{i}) + \sum_{j=1}^{t}b_{j}z_{ij} + e_{i} $$ where $y_{i}$ is the $i$-th observation's respone, $\mathbf{x}_{i}=(x_{i1},\ldots,x_{ip})$ is the $i$-th observation's nonparametric predictor vector, $\eta$ is an unknown smooth function relating the response and nonparametric predictors, $\mathbf{z}_{i}=(z_{i1},\ldots,z_{it})$ is the $i$-th observation's parametric predictor vector, and $e_{i}\sim\mathrm{N}(0,\sigma^{2})$ is iid Gaussian error. Function can fit both additive and interactive non/parametric effects, and allows for 2-way and 3-way interactions between nonparametric and parametric effects (see Details and Examples). ```r bigssp(formula,data=NULL,type=NULL,nknots=NULL,rparm=NA, lambdas=NULL,skip.iter=TRUE,se.fit=FALSE,rseed=1234, gcvopts=NULL,knotcheck=TRUE,thetas=NULL,weights=NULL, random=NULL,remlalg=c("FS","NR","EM","none"),remliter=500, remltol=10^-4,remltau=NULL) ``` ## Arguments - `formula`: An object of class "`formula`": a symbolic description of the model to be fitted (see Details and Examples for more information). - `data`: Optional data frame, list, or environment containing the variables in `formula`. Or an object of class "makessp", which is output from `makessp`. - `type`: List of smoothing spline types for predictors in `formula` (see Details). Options include `type="cub"` for cubic, `type="cub0"` for another cubic, `type="per"` for cubic periodic, `type="tps"` for cubic thin-plate, `type="ord"` for ordinal, and `type="nom"` for nominal. Use `type="prm"` for parametric effect. - `nknots`: Two possible options: (a) scalar giving total number of random knots to sample, or (b) vector indexing which rows of `data` to use as knots. - `rparm`: List of rounding parameters for each predictor. See Details. - `lambdas`: Vector of global smoothing parameters to try. Default `lambdas=10^-c(9:0)`. - `skip.iter`: Logical indicating whether to skip the iterative smoothing parameter update. Using `skip.iter=FALSE` should provide a more optimal solution, but the fitting time may be substantially longer. See Skip Iteration section. - `se.fit`: Logical indicating if the standard errors of the fitted values should be estimated. - `rseed`: Random seed for knot sampling. Input is ignored if `nknots` is an input vector of knot indices. Set `rseed=NULL` to obtain a different knot sample each time, or set `rseed` to any positive integer to use a different seed than the default. - `gcvopts`: Control parameters for optimization. List with 3 elements: (a) `maxit`: maximum number of algorithm iterations, (b) `gcvtol`: covergence tolerance for iterative GCV update, and (c) `alpha`: tuning parameter for GCV minimization. Default: `gcvopts=list(maxit=5,gcvtol=10^-5,alpha=1)` - `knotcheck`: If `TRUE`, only unique knots are used (for stability). - `thetas`: List of initial smoothing parameters for each predictor subspace. See Details. - `weights`: Vector of positive weights for fitting (default is vector of ones). - `random`: Adds random effects to model (see Random Effects section). - `remlalg`: REML algorithm for estimating variance components (see Random Effects section). Input is ignored if `random=NULL`. - `remliter`: Maximum number of iterations for REML estimation of variance components. Input is ignored if `random=NULL`. - `remltol`: Convergence tolerance for REML estimation of variance components. Input is ignored if `random=NULL`. - `remltau`: Initial estimate of variance parameters for REML estimation of variance components. Input is ignored if `random=NULL`. ## Details The `formula` syntax is similar to that used in `lm` and many other R regression functions. Use `y~x` to predict the response `y` from the predictor `x`. Use `y~x1+x2` to fit an additive model of the predictors `x1` and `x2`, and use `y~x1*x2` to fit an interaction model. The syntax `y~x1*x2` includes the interaction and main effects, whereas the syntax `y~x1:x2` only includes the interaction. See Computational Details for specifics about how non/parametric effects are estimated. See `bigspline` for definitions of `type="cub"`, `type="cub0"`, and `type="per"` splines, which can handle one-dimensional predictors. See Appendix of Helwig and Ma (2015) for information about `type="tps"` and `type="nom"` splines. Note that `type="tps"` can handle one-, two-, or three-dimensional predictors. I recommend using `type="cub"` if the predictor scores have no extreme outliers; when outliers are present, `type="tps"` may produce a better result. Using the rounding parameter input `rparm` can greatly speed-up and stabilize the fitting for large samples. For typical cases, I recommend using `rparm=0.01` for cubic and periodic splines, but smaller rounding parameters may be needed for particularly jagged functions. For thin-plate splines, the data are NOT transformed to the interval [0,1] before fitting, so the rounding parameter should be on the raw data scale. Also, for `type="tps"` you can enter one rounding parameter for each predictor dimension. Use `rparm=1` for ordinal and nominal splines. ## Returns - **fitted.values**: Vector of fitted values corresponding to the original data points in `xvars` (if `rparm=NA`) or the rounded data points in `xunique` (if `rparm` is used). - **se.fit**: Vector of standard errors of `fitted.values` (if input `se.fit=TRUE)`. - **yvar**: Response vector. - **xvars**: List of predictors. - **type**: Type of smoothing spline that was used for each predictor. - **yunique**: Mean of `yvar` for unique points after rounding (if `rparm` is used). - **xunique**: Unique rows of `xvars` after rounding (if `rparm` is used). - **sigma**: Estimated error standard deviation, i.e., $\hat{\sigma}$. - **ndf**: Data frame with two elements: `n` is total sample size, and `df` is effective degrees of freedom of fit model (trace of smoothing matrix). - **info**: Model fit information: vector containing the GCV, multiple R-squared, AIC, and BIC of fit model (assuming Gaussian error). - **modelspec**: List containing specifics of fit model (needed for prediction). - **converged**: Convergence status: `converged=TRUE` if iterative update converged, `converged=FALSE` if iterative update failed to converge, and `converged=NA` if option `skip.iter=TRUE` was used. - **tnames**: Names of the terms in model. - **random**: Random effects formula (same as input). - **tau**: Variance parameters such that `sigma*sqrt(tau)` gives standard deviation of random effects (if `!is.null(random)`). - **blup**: Best linear unbiased predictors (if `!is.null(random)`). - **call**: Called model in input `formula`. ## References Gu, C. (2013). **Smoothing spline ANOVA models, 2nd edition**. New York: Springer. Gu, C. and Wahba, G. (1991). Minimizing GCV/GML scores with multiple smoothing parameters via the Newton method. **SIAM Journal on Scientific and Statistical Computing, 12**, 383-398. Helwig, N. E. (2013). **Fast and stable smoothing spline analysis of variance models for large samples with applications to electroencephalography data analysis**. Unpublished doctoral dissertation. University of Illinois at Urbana-Champaign. Helwig, N. E. (2016). Efficient estimation of variance components in nonparametric mixed-effects models with large samples. **Statistics and Computing, 26**, 1319-1336. Helwig, N. E. (2017). [Regression with ordered predictors via ordinal smoothing splines](http://dx.doi.org/10.3389/fams.2017.00015). Frontiers in Applied Mathematics and Statistics, 3(15), 1-13. Helwig, N. E. and Ma, P. (2015). Fast and stable multiple smoothing parameter selection in smoothing spline analysis of variance models with large samples. **Journal of Computational and Graphical Statistics, 24**, 715-732. Helwig, N. E. and Ma, P. (2016). Smoothing spline ANOVA for super-large samples: Scalable computation via rounding parameters. **Statistics and Its Interface, 9**, 433-444. ## Author(s) Nathaniel E. Helwig ## Note The spline is estimated using penalized least-squares, which does not require the Gaussian error assumption. However, the spline inference information (e.g., standard errors and fit information) requires the Gaussian error assumption. ## Warnings Cubic and cubic periodic splines transform the predictor to the interval [0,1] before fitting. When using rounding parameters, output `fitted.values` corresponds to unique rounded predictor scores in output `xunique`. Use `predict.bigssp` function to get fitted values for full `yvar` vector. ## Computational Details To estimate $\eta$ I minimize the penalized least-squares functional $$ \frac{1}{n}\sum_{i=1}^{n}\left(y_{i} - \eta(\mathbf{x}_{i}) - \textstyle\sum_{j=1}^{t}b_{j}z_{ij} \right)^{2} + \lambda J(\eta) $$ where $J(\cdot)$ is a nonnegative penalty functional quantifying the roughness of $\eta$ and $\lambda>0$ is a smoothing parameter controlling the trade-off between fitting and smoothing the data. Note that for $p>1$ nonparametric predictors, there are additional $\theta_{k}$ smoothing parameters embedded in $J$. The penalized least squares functioncal can be rewritten as $$ \|\mathbf{y} - \mathbf{K}\mathbf{d} - \mathbf{J}_{\theta}\mathbf{c}\|^{2} + n\lambda\mathbf{c}'\mathbf{Q}_{\theta}\mathbf{c} $$ where $\mathbf{K}=\{\phi(x_{i}),\mathbf{z}_{i}\}_{n \times m}$ is the parametric space basis function matrix, $\mathbf{J}_{\theta}=\sum_{k=1}^{s}\theta_{k}\mathbf{J}_{k}$ with $\mathbf{J}_{k}=\{\rho_{k}(\mathbf{x}_{i},\mathbf{x}_{h}^{*})\}_{n \times q}$ denoting the $k$-th contrast space basis funciton matrix, $\mathbf{Q}_{\theta}=\sum_{k=1}^{s}\theta_{k}\mathbf{Q}_{k}$ with $\mathbf{Q}_{k}=\{\rho_{k}(\mathbf{x}_{g}^{*},\mathbf{x}_{h}^{*})\}_{q \times q}$ denoting the $k$-th penalty matrix, and $\mathbf{d}=(d_{0},\ldots,d_{m})'$ and $\mathbf{c}=(c_{1},\ldots,c_{q})'$ are the unknown basis function coefficients. The optimal smoothing parameters are chosen by minimizing the GCV score (see `bigspline`). Note that this function uses the classic smoothing spline parameterization (see Gu, 2013), so there is more than one smoothing parameter per predictor (if interactions are included in the model). To evaluate the GCV score, this function uses the improved (scalable) SSA algorithm discussed in Helwig (2013) and Helwig and Ma (2015). ## Skip Iteration For $p>1$ predictors, initial values for the $\theta_{k}$ parameters are estimated using Algorithm 3.2 described in Gu and Wahba (1991). Default use of this function (`skip.iter=TRUE`) fixes the $\theta_{k}$ parameters afer the smart start, and then finds the global smoothing parameter $\lambda$ (among the input `lambdas`) that minimizes the GCV score. This approach typically produces a solution very similar to the more optimal solution using `skip.iter=FALSE`. Setting `skip.iter=FALSE` uses the same smart starting algorithm as setting `skip.iter=TRUE`. However, instead of fixing the $\theta_{k}$ parameters afer the smart start, using `skip.iter=FALSE` iterates between estimating the optimal $\lambda$ and the optimal $\theta_{k}$ parameters. The R function `nlm` is used to minimize the GCV score with respect to the $\theta_{k}$ parameters, which can be time consuming for models with many predictors. ## Random Effects The input `random` adds random effects to the model assuming a variance components structure. Both nested and crossed random effects are supported. In all cases, the random effects are assumed to be indepedent zero-mean Gaussian variables with the variance depending on group membership. Random effects are distinguished by vertical bars ("|"), which separate expressions for design matrices (left) from group factors (right). For example, the syntax `~1|group` includes a random intercept for each level of `group`, whereas the syntax `~1+x|group` includes both a random intercept and a random slope for each level of `group`. For crossed random effects, parentheses are needed to distinguish different terms, e.g., `~(1|group1)+(1|group2)` includes a random intercept for each level of `group1` and a random intercept for each level of `group2`, where both `group1` and `group2` are factors. For nested random effects, the syntax `~group|subject` can be used, where both `group` and `subject` are factors such that the levels of `subject` are nested within those of `group`. The input `remlalg` determines the REML algorithm used to estimate the variance components. Setting `remlalg="FS"` uses a Fisher Scoring algorithm (default). Setting `remlalg="NR"` uses a Newton-Raphson algorithm. Setting `remlalg="EM"` uses an Expectation Maximization algorithm. Use `remlalg="none"` to fit a model with known variance components (entered through `remltau`). The input `remliter` sets the maximum number of iterations for the REML estimation. The input `remltol` sets the convergence tolerance for the REML estimation, which is determined via relative change in the REML log-likelihood. The input `remltau` sets the initial estimates of variance parameters; default is `remltau = rep(1,ntau)` where `ntau` is the number of variance components. ## Examples ```r ########## EXAMPLE ########## # function with four continuous predictors set.seed(773) myfun <- function(x1v,x2v,x3v,x4v){ sin(2*pi*x1v) + log(x2v+.1) + x3v*cos(pi*(x4v)) } x1v <- runif(500) x2v <- runif(500) x3v <- runif(500) x4v <- runif(500) y <- myfun(x1v,x2v,x3v,x4v) + rnorm(500) # fit cubic spline model with x3v*x4v interaction and x3v as "cub" # (includes x3v and x4v main effects) cubmod <- bigssp(y~x1v+x2v+x3v*x4v,type=list(x1v="cub",x2v="cub",x3v="cub",x4v="cub"),nknots=50) crossprod( myfun(x1v,x2v,x3v,x4v) - cubmod$fitted.values )/500 # fit cubic spline model with x3v*x4v interaction and x3v as "cub0" # (includes x3v and x4v main effects) cubmod <- bigssp(y~x1v+x2v+x3v*x4v,type=list(x1v="cub",x2v="cub",x3v="cub0",x4v="cub"),nknots=50) crossprod( myfun(x1v,x2v,x3v,x4v) - cubmod$fitted.values )/500 # fit model with x3v*x4v interaction treating x3v as parametric effect # (includes x3v and x4v main effects) cubmod <- bigssp(y~x1v+x2v+x3v*x4v,type=list(x1v="cub",x2v="cub",x3v="prm",x4v="cub"),nknots=50) crossprod( myfun(x1v,x2v,x3v,x4v) - cubmod$fitted.values )/500 # fit cubic spline model with x3v:x4v interaction and x3v as "cub" # (excludes x3v and x4v main effects) cubmod <- bigssp(y~x1v+x2v+x3v:x4v,type=list(x1v="cub",x2v="cub",x3v="cub",x4v="cub"),nknots=50) crossprod( myfun(x1v,x2v,x3v,x4v) - cubmod$fitted.values )/500 # fit cubic spline model with x3v:x4v interaction and x3v as "cub0" # (excludes x3v and x4v main effects) cubmod <- bigssp(y~x1v+x2v+x3v:x4v,type=list(x1v="cub",x2v="cub",x3v="cub0",x4v="cub"),nknots=50) crossprod( myfun(x1v,x2v,x3v,x4v) - cubmod$fitted.values )/500 # fit model with x3v:x4v interaction treating x3v as parametric effect # (excludes x3v and x4v main effects) cubmod <- bigssp(y~x1v+x2v+x3v:x4v,type=list(x1v="cub",x2v="cub",x3v="prm",x4v="cub"),nknots=50) crossprod( myfun(x1v,x2v,x3v,x4v) - cubmod$fitted.values )/500 ```

bigssp function