--- title: "Introduction to cmaRs" output: rmarkdown::pdf_document vignette: > %\VignetteIndexEntry{Introduction to cmaRs} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` This package is designed to construct Conic Multivariate Adaptive Regression Splines (CMARS) in R. CMARS model is a linear combination of basis functions that is taken from the forward part of the MARS algorithm. CMARS model parameters are obtained after solving Conic Quadratic Programming (CQP) which includes two cone constraints to handle accuracy and complexity of the model. Please check Weber, et al. (2011) CMARS: a new contribution to nonparametric regression with multivariate adaptive regression splines supported by continuous optimization. Inverse Problems in Science and Engineering, 2012, 20.3: 371-400. doi:10.1080/17415977.2011.62477, for more details. Prediction and binary classification models can be constructed with cmaRs. Note that **in order to construct CMARS models, both MOSEK software and Rmosek package needed. Please follow carefully the steps available in https://docs.mosek.com/latest/rmosek/install-interface.html for successful installation.** In order to construct CMARS models, a user can define the following arguments * **formula**: description of the model * **degree**: maximum degree of interaction * **nk**: maximum number of model terms before pruning * **data**: data frame * **classification**: binary variable indicating whether the model is for prediction or binary classification * **threshold.class**: the threshold value that is used to convert probabilities to classes in binary classification models. ```{r setup} library(cmaRs) ``` This package includes three data sets; preddata.std, classdata.std and table.b6. The first two data sets are taken from UCI: Machine Learning Repository (available at http://archive.ics.uci.edu/ml/), preprocessed and standardized. The first one is used for prediction and the other one is used for binary classification. Moreover, there is one more data set called as table.b6 which is directly taken from the "MPV" package (version 1.58) (Braun, and MacQueen, 2023). * The preddata.std includes `r nrow(preddata.std)` observations and `r ncol(preddata.std)` variables. ```{r} data(preddata.std) head(preddata.std) ``` * The classdata.std includes `r nrow(classdata.std)` observations and `r ncol(classdata.std)` variables. ```{r} data(classdata.std) head(classdata.std) ``` * The table.b6 includes `r nrow(table.b6)` observations and `r ncol(table.b6)` variables. ```{r} data(table.b6) head(table.b6) ``` # The Prediction Modeling As an example, the following prediction CMARS model for the data set trees can be constructed with the following snippets. ```{r echo=TRUE, results='hide', eval=FALSE} prediction.model <- cmaRs(Volume ~ ., degree = 2, nk = 20, data = trees) ``` In order to study this model in detail, the **summary** function can be used. Note here that the CMARS model is determined with **nk** = 20 and **degree** = 2. Here, the final model contains six basis function in total with one interaction term. ```{r, eval=FALSE} summary(prediction.model) #> Call: #> cmaRs(formula = Volume ~ ., data = trees, degree = 2, nk = 20) #> Volume = +29.1634 #> +4.9278 * pmax(0,Girth-14.2) #> -3.2309 * pmax(0,14.2-Girth) #> +0.7313 * pmax(0,Height-75) #> -0.1684 * pmax(0,75-Height) #> +0.1312 * pmax(0,Girth-8.3)*pmax(0,Height-75) #> -1.2977 * pmax(0,Height-78) #> R2 0.9793 r 0.9896 RSS 168.0029 ``` Some performance measures are also printed at the end of the output. For instance, the $R^2$, r and RSS values are given for the prediction models. It is also possible to construct several graphs of a prediction model. Here, **plot** function helps for this. ```{r eval=FALSE} fig1 <- plot(prediction.model) ``` ```{r fig1, echo=FALSE, fig.height=12} knitr::include_graphics("Figure1.png") ``` Another function of this package is **predict** which calculates the fitted values of a CMARS model. An example is given below. ```{r eval=FALSE} predict(prediction.model) #> [1] 9.259021 9.386200 9.695543 16.703820 20.239753 20.164926 17.308760 #> [8] 18.824519 22.046053 19.470698 22.996184 21.255036 21.255036 20.075638 #> [15] 22.055412 24.794797 29.229259 31.136370 26.874258 26.018423 32.955379 #> [22] 34.096100 30.473296 37.528153 43.074254 52.021298 53.805352 54.756953 #> [29] 55.315354 55.315354 77.168803 ``` # The Classification Modeling cmaRs can also construct binary classification models. An example is given below with its snippets. ```{r results='hide', warning=FALSE, eval=FALSE, message=FALSE} library(earth) library(cmaRs) classification.model <- cmaRs(survived ~ age, nk = 35, classification = TRUE, data = etitanic) ``` Note that, the **classification** argument is set as TRUE here which indicates a binary classification model. Moreover, the degree argument is used as 1 which is its default value indicating a main effect model. This model is constructed by using only continuous variable in the data set, age. Similar to the previous example, the **summary** function presents the details of the model. ```{r eval=FALSE} summary(classification.model) #> Call: #> cmaRs(formula = survived ~ age, data = etitanic, classification = TRUE, #> nk = 35) #> #> survived = -4.9489 #> -0.3084 * pmax(0,age-18) #> +0.3321 * pmax(0,18-age) #> -0.427 * pmax(0,age-53) #> +0.3564 * pmax(0,age-67) #> -0.3291 * pmax(0,age-64) #> +1.1664 * pmax(0,age-46) #> +0.7742 * pmax(0,age-57) #> -0.1451 * pmax(0,age-35) #> -0.7469 * pmax(0,age-58) #> +0.2288 * pmax(0,age-61) #> +0.0725 * pmax(0,age-41) #> -1.0824 * pmax(0,age-45) #> -0.5147 * pmax(0,age-48) #> +0.3658 * pmax(0,age-51) #> +0.281 * pmax(0,age-44) #> +0.0994 * pmax(0,age-34) #> +0.6316 * pmax(0,age-2) #> -0.32 * pmax(0,age-3) #> #> AUC 0.6221 MCR 0.3681 PCC 0.6319 precision 0.6339 recall 0.895 specificity 0.2506 ``` The model, here, includes terms with no interactions. Moreover, some performance measures for this binary classification are given in the last line. Some of them such as "Misclassification Rate" is calculated by taking the threshold as 0.5 as the default value. Another value of the threshold can also be set by using the argument, **threshold.class**, available in the cmaRs function. The following example is given to examplify this. ```{r echo=TRUE, results='hide', warning=FALSE, eval=FALSE} classification.model.threshold <- cmaRs(survived ~ age, nk = 35, classification = TRUE, data = etitanic, threshold.class = 0.1, Auto.linpreds = FALSE) ``` The following binary classification model follows the same construction steps but with different performance values. ```{r eval=FALSE} summary(classification.model.threshold) #> Call: #> cmaRs(formula = survived ~ age, data = etitanic, classification = TRUE, #> threshold.class = 0.1, nk = 35, Auto.linpreds = FALSE) #> #> survived = -4.9489 #> -0.3084 * pmax(0,age-18) #> +0.3321 * pmax(0,18-age) #> -0.427 * pmax(0,age-53) #> +0.3564 * pmax(0,age-67) #> -0.3291 * pmax(0,age-64) #> +1.1664 * pmax(0,age-46) #> +0.7742 * pmax(0,age-57) #> -0.1451 * pmax(0,age-35) #> -0.7469 * pmax(0,age-58) #> +0.2288 * pmax(0,age-61) #> +0.0725 * pmax(0,age-41) #> -1.0824 * pmax(0,age-45) #> -0.5147 * pmax(0,age-48) #> +0.3658 * pmax(0,age-51) #> +0.281 * pmax(0,age-44) #> +0.0994 * pmax(0,age-34) #> +0.6316 * pmax(0,age-2) #> -0.32 * pmax(0,age-3) #> #> AUC 0.6221 MCR 0.5746 PCC 0.4254 precision 1 recall 0.0291 specificity 1 ``` In addition to the performance values of the models, it is also possible to construct the ROC curve for the classification model using the **plot** function as follows. ```{r eval=FALSE} fig2 <- plot(classification.model) ``` ```{r fig2, echo=FALSE, fig.height=6} knitr::include_graphics("Figure2.png") ``` # Components of cmaRs Object One of the data sets, named as table.b6, inherited from the package MPV (Braun, 2018), with 28 observations and four predictor variables as x1, x2, x3 and x4 with the response variable; y. ```{r} data(table.b6) head(table.b6) ``` The following snippets demonstrate how to fit and obtain a CMARS prediction model (and cmaRs object) including all predictor variables available in the data set. ```{r echo=TRUE, results='hide', eval=FALSE} model.ex1 <- cmaRs(y~., degree = 2, nk = 20, classification = FALSE, Auto.linpreds = FALSE, data = table.b6) summary(model.ex1) # y= +0.003 # -0.003 * pmax(0,x4-0.0177) # -0.0798 * pmax(0,0.0177-x4) # 0 * pmax(0,x2-436.9) # 0 * pmax(0,436.9-x2) # -7e-04 * pmax(0,x1-0.0044)*pmax(0,436.9-x2) # -2.6911 * pmax(0,x1-0.01) # +0.5728 * pmax(0,x3-0.0186) # +7e-04 * pmax(0,x1-0.0044)*pmax(0,x2-436.9) # +2.6131 * pmax(0,x1-0.0106) # -2e-04 * pmax(0,436.9-x2)*pmax(0,x4-0.0062) ``` From this point on, implementation details of cmaRs package are examplified. Ten BFs (i.e. nk = 10) built in the forward CMARS algorithm (see open source code; Step 1: Forward step) are listed below. ```{r echo=TRUE, results='hide', eval=FALSE} model.ex1$bf.cmars # [1] "pmax(0,x4-0.0177)" "pmax(0,0.0177-x4)" # [3] "pmax(0,x2-436.9)" "pmax(0,436.9-x2)" # [5] "pmax(0,x1-0.0044)*pmax(0,436.9-x2)" "pmax(0,x1-0.01)" # [7] "pmax(0,x3-0.0186)" "pmax(0,x1-0.0044)*pmax(0,x2-436.9)" # [9] "pmax(0,x1-0.0106)" "pmax(0,436.9-x2)*pmax(0,x4-0.0062)" ``` Here, the BFs `pmax(x4-0.0177)`, `pmax(0.0177-x4)`, `pmax(x2-436.9)`, `pmax(436.9-x2)`, `pmax(x1-0.01)`, `pmax(x3-0.0186)`, and `pmax(x1-0.0106)` represent the main effects and the remaining BFs (the 5th BF, the 8th BF, the 10th BF) `pmax(x1-0.0044)*pmax(436.9-x2)`, `pmax(x1-0.0044)*pmax(x2-436.9)`, and `pmax(436.9-x2)*pmax(x4-0.0062)` show the interaction effects of two predictor variables. The L matrix can be formed as an (11 × 11 in this model) dimensional matrix with the first column and first row having zero values, and the diagonal elements contain the $\sf{L_{m}}$ values for the corresponding BFs. ```{r echo=TRUE, results='hide', eval=FALSE} model.ex1$L # [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] # [1,] 0 0.0000000 0.0000000 0.00000 0.00000 0.0000 0.0000000 0.0000000 # [2,] 0 0.3978693 0.0000000 0.00000 0.00000 0.0000 0.0000000 0.0000000 # [3,] 0 0.0000000 0.3339162 0.00000 0.00000 0.0000 0.0000000 0.0000000 # [4,] 0 0.0000000 0.0000000 10.15874 0.00000 0.0000 0.0000000 0.0000000 # [5,] 0 0.0000000 0.0000000 0.00000 19.45508 0.0000 0.0000000 0.0000000 # [6,] 0 0.0000000 0.0000000 0.00000 0.00000 134.2702 0.0000000 0.0000000 # [7,] 0 0.0000000 0.0000000 0.00000 0.00000 0.0000 0.3196873 0.0000000 # [8,] 0 0.0000000 0.0000000 0.00000 0.00000 0.0000 0.0000000 0.3170173 # [9,] 0 0.0000000 0.0000000 0.00000 0.00000 0.0000 0.0000000 0.0000000 # [10,] 0 0.0000000 0.0000000 0.00000 0.00000 0.0000 0.0000000 0.0000000 # [11,] 0 0.0000000 0.0000000 0.00000 0.00000 0.0000 0.0000000 0.0000000 # [,9] [,10] [,11] # [1,] 0.00000 0.0000000 0.0000 # [2,] 0.00000 0.0000000 0.0000 # [3,] 0.00000 0.0000000 0.0000 # [4,] 0.00000 0.0000000 0.0000 # [5,] 0.00000 0.0000000 0.0000 # [6,] 0.00000 0.0000000 0.0000 # [7,] 0.00000 0.0000000 0.0000 # [8,] 0.00000 0.0000000 0.0000 # [9,] 11.53763 0.0000000 0.0000 # [10,] 0.00000 0.3182766 0.0000 # [11,] 0.00000 0.0000000 78.2987 ``` To illustrate some computational details, two BFs, the first (i.e. `pmax(x4-0.0177)`) and the tenth (i.e. `pmax(436.9-x2)*pmax(x4-0.0062)`), which represent main and interaction effects in the model, respectively, are selected. As stated beforehand, the cmaRs package automatically provides the derivatives matrix (DMS) and the variable matrix (VARMS) whose derivative exists in DMS matrix only for the last BF. In order to obtain these matrices for the other BFs, lets say, for the first one (i.e. `pmax(x4-0.0177)`), the code, ```{r echo=TRUE, results='hide', eval=FALSE} bfs_forward_names <- bfs_forward_names[1:2] ``` should be inserted into the middle of the following two lines of R codes available in the `cmaRs.fit` function under the "#Step 1: Forward step" section as follows: ```{r echo=TRUE, results='hide', eval=FALSE} bfs_forward_names <- rownames(model_mars$dirs) bfs_forward_names <- bfs_forward_names[1:2] bfs_forward_names.orig <- bfs_forward_names ``` After constructing the model again, one can also obtain the derivatives for the first BF given as a symbolic matrix form (i.e. DMS) as follows. ```{r echo=TRUE, results='hide', eval=FALSE} model.ex1$DMS # [,1] [,2] [,3] [,4] # [1,] "0" "0" "0" "0" # [2,] "0" "0" "0" "0" # [3,] "0" "0" "0" "0" # [4,] "0" "0" "0" "matrix(1, nrow <- (n+1), ncol <- 1)" ``` In addition, since this BF contains no interaction effect (i.e. just has one variable, x4), the first-order derivative of the function produces ones, and the second-order derivatives produces zeros with respect to the variable x4. The variable name whose derivative exists in the DMS matrix is obtained by the following matrix form, called VARMS. ```{r echo=TRUE, results='hide', eval=FALSE} model.ex1$VARMS # [,1] [,2] [,3] [,4] # [1,] "1" "1" "1" "1" # [2,] "1" "1" "1" "1" # [3,] "1" "1" "1" "1" # [4,] "1" "1" "1" "xplus4" ``` On the other hand, the tenth BF contains an interaction effect between the predictors x2 and x4. The first-order derivatives of the tenth BF are `(0.0062-xfirst4)` with respect to the variable x2 and `(436.9-xfirst2)` with respect to the variable x4. The second-order derivative of the tenth BF is the ```{r echo=TRUE, results='hide', eval=FALSE} # matrix(-1, nrow <- (n+1), ncol <- 1) ``` with respect to both of the variables x2 and x4. Note here that in any case where the order of the derivative is greater than two, it vanishes and thus, the full derivative matrix form of the tenth BF turns out to be the following. ```{r echo=TRUE, results='hide', eval=FALSE} model.ex1$DMS # [,1] [,2] [,3] [,4] # [1,] "0" "0" "0" "0" # [2,] "0" "0.0062 - xfirst4" "0" "matrix(-1, nrow <- (n+1), ncol <- 1)" # [3,] "0" "0" "0" "0" # [4,] "0" "0" "0" "436.9 - xfirst2" ``` The corresponding variable name matrix (i.e. VARMS matrix), simply provided by calling the model name as `model.ex1`, unless a change is made in the `cmaRs.fit` function. ```{r echo=TRUE, results='hide', eval=FALSE} model.ex1$VARMS # [,1] [,2] [,3] [,4] # [1,] "1" "1" "1" "1" # [2,] "1" "xplus2" "1" "1" # [3,] "1" "1" "1" "1" # [4,] "1" "1" "1" "xplus4" ``` The $\sf{L_{m}}$ (m = 1, 10) values corresponding to the first and the tenth BFs are $\sf{L_{1}}$ = 0.3978693 and $\sf{L_{10}}$ = 78.2987, respectively. After discritization is applied, which produces complete L matrix, the optimization part which solves the CQP takes place. For solving it, the package Rmosek (MOSEK-ApS, 2017) is utilized, and as a result, the model parameter estimates are obtained as follows. ```{r echo=TRUE, results='hide', eval=FALSE} model.ex1$coefficients # [1] 3.038673e-03 -2.996649e-03 -7.977544e-02 -6.809132e-06 3.517929e-06 # [6] -6.991181e-04 -2.691055e+00 5.728424e-01 6.722162e-04 2.613055e+00 # [11] -1.701914e-04 ``` Note here again that proper installation of MOSEK has a critical importance for successfully run of the package cmaRs, since Rmosek is used as the interface of it. Besides model parameter estimates, the cmaRs package also provides the fitted values belonging to the CMARS model built if needed. ```{r echo=TRUE, results='hide', eval=FALSE} > model.ex1$fitted.values # [1] 0.0007575982 0.0003071999 0.0004176579 0.0005434535 0.0004821474 0.0004517915 # [7] 0.0004376478 0.0004160878 0.0013674034 0.0012119509 0.0012961460 0.0011833885 # [13] 0.0009870342 0.0012386077 0.0013057215 0.0011887701 0.0012214500 0.0017378222 # [19] 0.0017998247 0.0018951991 0.0026293088 0.0031847542 0.0033385298 0.0033184915 # [25] 0.0013479697 0.0024916551 0.0027357523 0.0026836362 ```