sail is a package that fits a linear model with non-linear interactions via penalized maximum likelihood. The regularization path is computed at a grid of values for the regularization parameter \(\lambda\) and a fixed value of the second regularization parameter \(\alpha\). The method enforces the strong heredity property, i.e., an interaction is selected only if its corresponding main effects are also included. The interactions are limited to a single exposure variable, i.e., \(y \sim e + x_1 + x_2 + e*x_1 + e*x_2 + \epsilon\). Furthermore, this package allows a user-defined basis expansion on the \(x\) variables to allow for non-linear effects. The default is bsplines (e.g. splines::bs(x, 5)). It currently only fits linear models (binomial models are due in the next release).

Model

Let \(Y=(Y_1, \ldots, Y_n) \in \mathbb{R}^n\) be a continuous outcome variable, a binary or continuous environment vector, a matrix of predictors, and \(\varepsilon = (\varepsilon_1, \ldots, \varepsilon_n) \in \mathbb{R}^n\) a vector of i.i.d random variables with mean 0. Furthermore let \(f_j: \mathbb{R} \rightarrow \mathbb{R}\) be a smoothing method for variable \(X_j\) by a projection on to a set of basis functions: \[\begin{equation} f_j(X_j) = \sum_{\ell = 1}^{m_j} \psi_{j\ell}(X_j) \beta_{j\ell} \label{eq:smooth} \end{equation}\] Here, the \(\left\lbrace \psi_{j\ell} \right\rbrace_1^{m_j}\) are a family of basis functions in \(X_j\)~. Let \(\boldsymbol{\Psi}_j\) be the \(n \times m_j\) matrix of evaluations of the \(\psi_{j\ell}\) and for \(j = 1, \ldots, p\), i.e., \(\boldsymbol{\theta}_j\) is a \(m_j\)-dimensional column vector of basis coefficients for the \(j\)th main effect. In this article we consider an additive interaction regression model of the form \[\begin{align} Y & = \beta_0 \cdot \boldsymbol{1} + \sum_{j=1}^p \boldsymbol{\Psi}_j \boldsymbol{\theta}_j + \beta_E X_E + \sum_{j=1}^p (X_E \circ \boldsymbol{\Psi}_j) \boldsymbol{\alpha}_{j} + \varepsilon \label{eq:linpred} \end{align}\] where \(\beta_0\) is the intercept, \(\beta_E\) is the coefficient for the environment variable, \(\boldsymbol{\alpha}_j = (\alpha_{j1}, \ldots, \alpha_{jm_j})\in \mathbb{R}^{m_j}\) are the basis coefficients for the \(j\)th interaction term and \((X_E \circ \boldsymbol{\Psi}_j)\) is the \(n \times m_j\) matrix formed by the component-wise multiplication of the column vector \(X_E\) by each column of \(\boldsymbol{\Psi}_j\). To enforce the strong heredity property, we reparametrize the coefficients for the interaction terms in~ as \(\boldsymbol{\alpha}_{j} = \gamma_{j} \beta_E \boldsymbol{\theta}_j\): \[\begin{align} Y & = \beta_0 \cdot \boldsymbol{1} + \sum_{j=1}^p \boldsymbol{\Psi}_j \boldsymbol{\theta}_j + \beta_E X_E + \sum_{j=1}^p \gamma_{j} \beta_E (X_E \circ \boldsymbol{\Psi}_j) \boldsymbol{\theta}_j + \varepsilon \label{eq:linpred2} \end{align}\] For a continuous response, we use the squared-error loss: \[\begin{equation} \mathcal{L}(Y;\boldsymbol{\theta}) = \frac{1}{2n}\lVert Y - \beta_0 \cdot \boldsymbol{1} - \sum_{j=1}^p \boldsymbol{\Psi}_j \boldsymbol{\theta}_j - \beta_E X_E - \sum_{j=1}^p \gamma_{j} \beta_E (X_E \circ \boldsymbol{\Psi}_j) \boldsymbol{\theta}_j \rVert_2^2 \end{equation}\]

where \(\boldsymbol{\theta} \equiv (\beta_0, \beta_E,\boldsymbol{\theta}_1, \ldots, \boldsymbol{\theta}_p, \gamma_1, \ldots, \gamma_p)\).

We consider the following penalized least squares criterion for this problem: \[\begin{equation} \arg\min_{\boldsymbol{\theta} } \mathcal{L}(Y;\boldsymbol{\theta}) + \lambda (1-\alpha) \left( w_E |\beta_E| + \sum_{j=1}^{p} w_j \lVert\boldsymbol{\theta}_j \rVert_2 \right) + \lambda\alpha \sum_{j=1}^{p} w_{jE} |\gamma_{j}| \label{eq:lassolikelihood3} \end{equation}\]

where \(\lambda >0\) and \(\alpha \in (0,1)\) are tuning parameters and \(w_E, w_j, w_{jE}\) are adaptive weights for \(j=1, \ldots, p\). These weights serve as a way of allowing parameters to be penalized differently.

Installation

The package can be installed from GitHub via

install.packages("pacman")
pacman::p_load_gh('sahirbhatnagar/sail')

Quick Start

We give a quick overview of the main functions and go into details in other vignettes. We will use the simulated data which ships with the package and can be loaded via:

library(sail)
data("sailsim")
names(sailsim)
#> [1] "x"        "y"        "e"        "f1"       "f2"       "f3"      
#> [7] "f4"       "f3.inter" "f4.inter"

We first define a basis expansion. In this example we use cubic bsplines with degree 5.

library(splines)
f.basis <- function(x) splines::bs(x, degree = 5)

Next we fit the model using the most basic call to sail

fit <- sail(x = sailsim$x, y = sailsim$y, e = sailsim$e, basis = f.basis)

fit is an object of class sail that contains all the relevant information of the fitted model including the estimated coefficients at each value of \(\lambda\) (by default the program chooses its own decreasing sequence of 100 \(\lambda\) values). There are print, plot, coef and predict methods of objects of class sail. The print method outputs the following:

fit
#> 
#> Call:  sail(x = sailsim$x, y = sailsim$y, e = sailsim$e, basis = f.basis) 
#> 
#>      df_main df_interaction df_environment     %Dev  Lambda
#> s1         0              0              0 0.000000 1.48800
#> s2         0              0              1 0.001698 1.42100
#> s3         0              0              1 0.003357 1.35600
#> s4         0              0              1 0.004897 1.29400
#> s5         0              0              1 0.006352 1.23500
#> s6         1              0              1 0.014720 1.17900
#> s7         2              0              1 0.045810 1.12600
#> s8         2              0              1 0.080020 1.07500
#> s9         2              0              1 0.111900 1.02600
#> s10        2              0              1 0.141400 0.97910
#> s11        2              0              1 0.168900 0.93460
#> s12        2              0              1 0.194400 0.89210
#> s13        2              0              1 0.218100 0.85160
#> s14        2              0              1 0.240000 0.81290
#> s15        2              0              1 0.260400 0.77590
#> s16        2              0              1 0.279200 0.74070
#> s17        3              1              1 0.465100 0.70700
#> s18        3              1              1 0.479500 0.67490
#> s19        3              1              1 0.492900 0.64420
#> s20        3              2              1 0.508200 0.61490
#> s21        3              2              1 0.524200 0.58700
#> s22        3              2              1 0.536600 0.56030
#> s23        3              2              1 0.548000 0.53480
#> s24        4              2              1 0.558000 0.51050
#> s25        5              2              1 0.568900 0.48730
#> s26        5              2              1 0.579900 0.46520
#> s27        5              2              1 0.590100 0.44400
#> s28        5              2              1 0.599500 0.42380
#> s29        5              3              1 0.609800 0.40460
#> s30        5              3              1 0.619100 0.38620
#> s31        5              3              1 0.628200 0.36860
#> s32        6              3              1 0.637200 0.35190
#> s33        6              3              1 0.646800 0.33590
#> s34        6              3              1 0.655000 0.32060
#> s35        6              3              1 0.663300 0.30600
#> s36        7              3              1 0.671600 0.29210
#> s37        7              3              1 0.679500 0.27890
#> s38        7              3              1 0.687400 0.26620
#> s39        7              3              1 0.694600 0.25410
#> s40        7              4              1 0.702600 0.24250
#> s41        7              4              1 0.734100 0.23150
#> s42        7              4              1 0.739300 0.22100
#> s43        7              5              1 0.745200 0.21090
#> s44        7              6              1 0.754400 0.20140
#> s45        7              6              1 0.760700 0.19220
#> s46        7              6              1 0.767300 0.18350
#> s47        7              6              1 0.771800 0.17510
#> s48        7              6              1 0.776100 0.16720
#> s49        8              6              1 0.780700 0.15960
#> s50        8              6              1 0.784500 0.15230
#> s51        9              6              1 0.788700 0.14540
#> s52       10              6              1 0.792700 0.13880
#> s53       12              6              1 0.797800 0.13250
#> s54       13              6              1 0.802600 0.12650
#> s55       13              6              1 0.807500 0.12070
#> s56       13              6              1 0.812200 0.11520
#> s57       14              6              1 0.816700 0.11000
#> s58       14              6              1 0.821300 0.10500
#> s59       15              6              1 0.826100 0.10020
#> s60       15              7              1 0.831200 0.09566
#> s61       13             10              1 0.871400 0.09131
#> s62       13             10              1 0.873800 0.08716
#> s63       13             10              1 0.876200 0.08320
#> s64       14             10              1 0.878600 0.07942
#> s65       14             10              1 0.881000 0.07581
#> s66       14             11              1 0.883400 0.07236
#> s67       14             11              1 0.885900 0.06907
#> s68       14             11              1 0.888400 0.06593
#> s69       14             11              1 0.890600 0.06294
#> s70       16             13              1 0.930600 0.06008
#> s71       16             13              1 0.932300 0.05735
#> s72       16             13              1 0.934200 0.05474
#> s73       16             13              1 0.936100 0.05225
#> s74       16             13              1 0.938500 0.04988
#> s75       16             13              1 0.940000 0.04761
#> s76       17             13              1 0.942100 0.04545
#> s77       17             13              1 0.943800 0.04338
#> s78       17             14              1 0.945600 0.04141
#> s79       17             14              1 0.947300 0.03953
#> s80       18             14              1 0.949000 0.03773
#> s81       18             14              1 0.950700 0.03602
#> s82       18             14              1 0.952400 0.03438
#> s83       18             14              1 0.954400 0.03282
#> s84       18             14              1 0.955800 0.03132
#> s85       18             14              1 0.957200 0.02990
#> s86       19             14              1 0.959000 0.02854
#> s87       19             15              1 0.960300 0.02724
#> s88       19             15              1 0.967800 0.02601
#> s89       19             15              1 0.968900 0.02482
#> s90       19             15              1 0.970100 0.02370
#> s91       19             15              1 0.971300 0.02262
#> s92       19             15              1 0.972400 0.02159
#> s93       19             15              1 0.973500 0.02061
#> s94       19             15              1 0.974600 0.01967
#> s95       20             19              1 0.977900 0.01878
#> s96       20             19              1 0.978600 0.01792
#> s97       20             19              1 0.979300 0.01711
#> s98       20             19              1 0.980000 0.01633
#> s99       20             19              1 0.980700 0.01559
#> s100      20             19              1 0.981400 0.01488

When expand = TRUE (i.e. the user did not provide their own design matrix), the df_main and df_interaction columns correspond to the number of non-zero predictors present in the model before basis expansion. This does not correspond to the number of non-zero coefficients in the model, but rather the number of unique variables. In this example we expanded each column of \(\mathbf{X}\) to five columns. If df_main=4, df_interaction=2 and df_environment=1, then the total number of non-zero coefficients would be \(5 \times (4+2) + 1\).

The entire solution path can be plotted via the plot method for objects of class sail. The y-axis is the value of the coefficient and the x-axis is the \(\log(\lambda)\). Each line represents a coefficient in the model, and each color represents a variable (i.e. in this example a given variable will have 5 lines when it is non-zero). The numbers at the top of the plot represent the number of non-zero variables in the model: top panel (df_main + df_environment), bottom panel (df_interaction). The black line is the coefficient path for the environment variable.

plot(fit)

The estimated coefficients at each value of lambda is given by (matrix partially printed here for brevity)

coef(fit)[1:6,50:55]
#> 6 x 6 sparse Matrix of class "dgCMatrix"
#>                   s50        s51        s52        s53        s54
#> (Intercept)  5.290833  5.2837573  5.2803780  5.2753630  5.2717902
#> X1_1        -0.979303 -0.9604292 -0.9449933 -0.9220866 -0.9171202
#> X1_2         1.690176  1.7893506  1.8923047  1.9950535  2.1040903
#> X1_3         1.646209  1.7048994  1.7722054  1.8250656  1.8950572
#> X1_4         1.522426  1.5433270  1.5662846  1.5853967  1.6101828
#> X1_5         3.338527  3.4182125  3.4907005  3.5762497  3.6632614
#>                    s55
#> (Intercept)  5.2695429
#> X1_1        -0.9270781
#> X1_2         2.2057013
#> X1_3         1.9641879
#> X1_4         1.6321740
#> X1_5         3.7452560

The predicted response at each value of lambda:

predict(fit)[1:5,50:55]
#>            s50       s51       s52       s53       s54       s55
#> [1,]  6.244760  6.199363  6.185411  6.178014  6.156219  6.124321
#> [2,]  3.002808  2.995389  3.038672  3.079610  3.143616  3.208972
#> [3,]  2.073350  2.043516  2.016354  1.997213  1.966913  1.957284
#> [4,] 13.489022 13.490853 13.361077 13.384407 13.350715 13.324150
#> [5,]  1.225603  1.210385  1.134488  1.156340  1.156697  1.156134

The predicted response at a specific value of lambda can be specified by the s argument:

predict(fit, s = 0.8)
#>               1
#>   [1,] 5.624232
#>   [2,] 4.940944
#>   [3,] 3.847965
#>   [4,] 6.687777
#>   [5,] 3.058124
#>   [6,] 3.607326
#>   [7,] 6.615916
#>   [8,] 5.032047
#>   [9,] 6.040452
#>  [10,] 5.486751
#>  [11,] 4.765237
#>  [12,] 7.238408
#>  [13,] 7.063706
#>  [14,] 4.590884
#>  [15,] 7.657128
#>  [16,] 4.904077
#>  [17,] 4.195213
#>  [18,] 4.504339
#>  [19,] 4.809399
#>  [20,] 5.615791
#>  [21,] 6.546767
#>  [22,] 5.673886
#>  [23,] 4.631122
#>  [24,] 6.849338
#>  [25,] 4.091877
#>  [26,] 7.020022
#>  [27,] 7.175685
#>  [28,] 3.550668
#>  [29,] 6.737500
#>  [30,] 6.366243
#>  [31,] 5.192395
#>  [32,] 5.384846
#>  [33,] 6.935418
#>  [34,] 4.258696
#>  [35,] 7.603471
#>  [36,] 5.303784
#>  [37,] 5.303826
#>  [38,] 4.285530
#>  [39,] 4.803300
#>  [40,] 5.503811
#>  [41,] 3.620280
#>  [42,] 3.628351
#>  [43,] 5.421528
#>  [44,] 7.492914
#>  [45,] 4.183416
#>  [46,] 4.536360
#>  [47,] 5.348853
#>  [48,] 7.518597
#>  [49,] 4.493770
#>  [50,] 5.276594
#>  [51,] 5.568093
#>  [52,] 4.473407
#>  [53,] 5.563292
#>  [54,] 5.744353
#>  [55,] 4.592868
#>  [56,] 4.307980
#>  [57,] 4.854115
#>  [58,] 5.244134
#>  [59,] 6.736890
#>  [60,] 5.070075
#>  [61,] 4.517567
#>  [62,] 3.129270
#>  [63,] 6.601207
#>  [64,] 4.547151
#>  [65,] 4.566028
#>  [66,] 6.376860
#>  [67,] 4.967343
#>  [68,] 6.138383
#>  [69,] 5.385448
#>  [70,] 6.923559
#>  [71,] 5.015826
#>  [72,] 3.524949
#>  [73,] 5.669353
#>  [74,] 7.044308
#>  [75,] 3.523070
#>  [76,] 4.674963
#>  [77,] 3.396258
#>  [78,] 5.592357
#>  [79,] 5.415423
#>  [80,] 4.810314
#>  [81,] 5.712646
#>  [82,] 5.990818
#>  [83,] 4.861425
#>  [84,] 6.761552
#>  [85,] 5.493592
#>  [86,] 3.607781
#>  [87,] 3.986110
#>  [88,] 4.259476
#>  [89,] 4.174571
#>  [90,] 3.666310
#>  [91,] 5.181711
#>  [92,] 7.086435
#>  [93,] 5.406326
#>  [94,] 6.583031
#>  [95,] 4.664773
#>  [96,] 3.695451
#>  [97,] 5.100857
#>  [98,] 4.984838
#>  [99,] 3.961504
#> [100,] 6.078087

You can specify more than one value for s:

predict(fit, s = c(0.8, 0.2))
#>               1          2
#>   [1,] 5.624232  6.5230740
#>   [2,] 4.940944  2.9750177
#>   [3,] 3.847965  2.3266980
#>   [4,] 6.687777 13.9562204
#>   [5,] 3.058124  1.5688962
#>   [6,] 3.607326  1.6348912
#>   [7,] 6.615916 13.2009227
#>   [8,] 5.032047  5.9480404
#>   [9,] 6.040452  5.8691675
#>  [10,] 5.486751  5.0987940
#>  [11,] 4.765237  3.5631269
#>  [12,] 7.238408 11.8329816
#>  [13,] 7.063706 12.4826512
#>  [14,] 4.590884  3.9524402
#>  [15,] 7.657128 14.2931482
#>  [16,] 4.904077  5.0900767
#>  [17,] 4.195213 -0.5844299
#>  [18,] 4.504339  4.4666357
#>  [19,] 4.809399  3.0374049
#>  [20,] 5.615791  3.5399188
#>  [21,] 6.546767  8.3512494
#>  [22,] 5.673886  4.7077896
#>  [23,] 4.631122  1.9943630
#>  [24,] 6.849338 10.7983805
#>  [25,] 4.091877  0.9601846
#>  [26,] 7.020022  9.3762352
#>  [27,] 7.175685  8.6775432
#>  [28,] 3.550668 -1.6523278
#>  [29,] 6.737500 13.5513436
#>  [30,] 6.366243  2.2204467
#>  [31,] 5.192395  5.8213763
#>  [32,] 5.384846  7.4511268
#>  [33,] 6.935418  9.0277752
#>  [34,] 4.258696  4.5323823
#>  [35,] 7.603471 16.3907390
#>  [36,] 5.303784  0.1377432
#>  [37,] 5.303826  5.6261096
#>  [38,] 4.285530  2.3206999
#>  [39,] 4.803300  3.8195444
#>  [40,] 5.503811  4.8408086
#>  [41,] 3.620280  2.4050059
#>  [42,] 3.628351  0.4144305
#>  [43,] 5.421528  6.4970752
#>  [44,] 7.492914 15.6998648
#>  [45,] 4.183416  3.1368088
#>  [46,] 4.536360  3.5137943
#>  [47,] 5.348853  4.4711267
#>  [48,] 7.518597 18.0122252
#>  [49,] 4.493770  2.3526655
#>  [50,] 5.276594  7.4618284
#>  [51,] 5.568093  5.9679862
#>  [52,] 4.473407  1.7229550
#>  [53,] 5.563292  7.3990744
#>  [54,] 5.744353  8.6534271
#>  [55,] 4.592868  2.0907568
#>  [56,] 4.307980  0.7927840
#>  [57,] 4.854115  5.3436578
#>  [58,] 5.244134  9.2910545
#>  [59,] 6.736890  9.7709887
#>  [60,] 5.070075  1.6897955
#>  [61,] 4.517567  1.4507809
#>  [62,] 3.129270 -0.9532221
#>  [63,] 6.601207  8.1123225
#>  [64,] 4.547151  1.3003285
#>  [65,] 4.566028  3.1206603
#>  [66,] 6.376860  3.5765440
#>  [67,] 4.967343  4.8268433
#>  [68,] 6.138383  6.1849972
#>  [69,] 5.385448 10.1092482
#>  [70,] 6.923559  8.2999571
#>  [71,] 5.015826  4.0040478
#>  [72,] 3.524949  1.3597820
#>  [73,] 5.669353  6.7846103
#>  [74,] 7.044308  4.5323126
#>  [75,] 3.523070  2.6555577
#>  [76,] 4.674963  1.1032379
#>  [77,] 3.396258  1.8733145
#>  [78,] 5.592357 13.3777530
#>  [79,] 5.415423  6.6235012
#>  [80,] 4.810314  2.9584783
#>  [81,] 5.712646  4.0235765
#>  [82,] 5.990818  4.9867095
#>  [83,] 4.861425  4.5053610
#>  [84,] 6.761552  5.9701019
#>  [85,] 5.493592  4.8052874
#>  [86,] 3.607781  1.8020652
#>  [87,] 3.986110  3.4415117
#>  [88,] 4.259476  0.7337290
#>  [89,] 4.174571  4.3173319
#>  [90,] 3.666310  5.4733594
#>  [91,] 5.181711  8.2829841
#>  [92,] 7.086435 10.2504954
#>  [93,] 5.406326  5.2240663
#>  [94,] 6.583031  7.7020652
#>  [95,] 4.664773  5.3974401
#>  [96,] 3.695451  2.1516883
#>  [97,] 5.100857  9.2733043
#>  [98,] 4.984838  6.5142889
#>  [99,] 3.961504  3.3387104
#> [100,] 6.078087  4.9794667

You can also extract a list of active variables (i.e. variables with a non-zero estimated coefficient) for each value of lambda:

fit[["active"]]
#> [[1]]
#> character(0)
#> 
#> [[2]]
#> [1] "E"
#> 
#> [[3]]
#> [1] "E"
#> 
#> [[4]]
#> [1] "E"
#> 
#> [[5]]
#> [1] "E"
#> 
#> [[6]]
#> [1] "X4" "E" 
#> 
#> [[7]]
#> [1] "X3" "X4" "E" 
#> 
#> [[8]]
#> [1] "X3" "X4" "E" 
#> 
#> [[9]]
#> [1] "X3" "X4" "E" 
#> 
#> [[10]]
#> [1] "X3" "X4" "E" 
#> 
#> [[11]]
#> [1] "X3" "X4" "E" 
#> 
#> [[12]]
#> [1] "X3" "X4" "E" 
#> 
#> [[13]]
#> [1] "X3" "X4" "E" 
#> 
#> [[14]]
#> [1] "X3" "X4" "E" 
#> 
#> [[15]]
#> [1] "X3" "X4" "E" 
#> 
#> [[16]]
#> [1] "X3" "X4" "E" 
#> 
#> [[17]]
#> [1] "X1"   "X3"   "X4"   "X4:E" "E"   
#> 
#> [[18]]
#> [1] "X1"   "X3"   "X4"   "X4:E" "E"   
#> 
#> [[19]]
#> [1] "X1"   "X3"   "X4"   "X4:E" "E"   
#> 
#> [[20]]
#> [1] "X1"   "X3"   "X4"   "X3:E" "X4:E" "E"   
#> 
#> [[21]]
#> [1] "X1"   "X3"   "X4"   "X3:E" "X4:E" "E"   
#> 
#> [[22]]
#> [1] "X1"   "X3"   "X4"   "X3:E" "X4:E" "E"   
#> 
#> [[23]]
#> [1] "X1"   "X3"   "X4"   "X3:E" "X4:E" "E"   
#> 
#> [[24]]
#> [1] "X1"   "X3"   "X4"   "X8"   "X3:E" "X4:E" "E"   
#> 
#> [[25]]
#> [1] "X1"   "X3"   "X4"   "X8"   "X11"  "X3:E" "X4:E" "E"   
#> 
#> [[26]]
#> [1] "X1"   "X3"   "X4"   "X8"   "X11"  "X3:E" "X4:E" "E"   
#> 
#> [[27]]
#> [1] "X1"   "X3"   "X4"   "X8"   "X11"  "X3:E" "X4:E" "E"   
#> 
#> [[28]]
#> [1] "X1"   "X3"   "X4"   "X8"   "X11"  "X3:E" "X4:E" "E"   
#> 
#> [[29]]
#> [1] "X1"   "X3"   "X4"   "X8"   "X11"  "X1:E" "X3:E" "X4:E" "E"   
#> 
#> [[30]]
#> [1] "X1"   "X3"   "X4"   "X8"   "X11"  "X1:E" "X3:E" "X4:E" "E"   
#> 
#> [[31]]
#> [1] "X1"   "X3"   "X4"   "X8"   "X11"  "X1:E" "X3:E" "X4:E" "E"   
#> 
#> [[32]]
#>  [1] "X1"   "X2"   "X3"   "X4"   "X8"   "X11"  "X1:E" "X3:E" "X4:E" "E"   
#> 
#> [[33]]
#>  [1] "X1"   "X2"   "X3"   "X4"   "X8"   "X11"  "X1:E" "X3:E" "X4:E" "E"   
#> 
#> [[34]]
#>  [1] "X1"   "X2"   "X3"   "X4"   "X8"   "X11"  "X1:E" "X3:E" "X4:E" "E"   
#> 
#> [[35]]
#>  [1] "X1"   "X2"   "X3"   "X4"   "X8"   "X11"  "X1:E" "X3:E" "X4:E" "E"   
#> 
#> [[36]]
#>  [1] "X1"   "X2"   "X3"   "X4"   "X8"   "X11"  "X20"  "X1:E" "X3:E" "X4:E"
#> [11] "E"   
#> 
#> [[37]]
#>  [1] "X1"   "X2"   "X3"   "X4"   "X8"   "X11"  "X20"  "X1:E" "X3:E" "X4:E"
#> [11] "E"   
#> 
#> [[38]]
#>  [1] "X1"   "X2"   "X3"   "X4"   "X8"   "X11"  "X20"  "X1:E" "X3:E" "X4:E"
#> [11] "E"   
#> 
#> [[39]]
#>  [1] "X1"   "X2"   "X3"   "X4"   "X8"   "X11"  "X20"  "X1:E" "X3:E" "X4:E"
#> [11] "E"   
#> 
#> [[40]]
#>  [1] "X1"   "X2"   "X3"   "X4"   "X8"   "X11"  "X20"  "X1:E" "X2:E" "X3:E"
#> [11] "X4:E" "E"   
#> 
#> [[41]]
#>  [1] "X1"   "X2"   "X3"   "X4"   "X8"   "X11"  "X20"  "X1:E" "X2:E" "X3:E"
#> [11] "X4:E" "E"   
#> 
#> [[42]]
#>  [1] "X1"   "X2"   "X3"   "X4"   "X8"   "X11"  "X20"  "X1:E" "X2:E" "X3:E"
#> [11] "X4:E" "E"   
#> 
#> [[43]]
#>  [1] "X1"   "X2"   "X3"   "X4"   "X8"   "X11"  "X20"  "X1:E" "X2:E" "X3:E"
#> [11] "X4:E" "X8:E" "E"   
#> 
#> [[44]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X8"    "X11"   "X20"   "X1:E" 
#>  [9] "X2:E"  "X3:E"  "X4:E"  "X8:E"  "X11:E" "E"    
#> 
#> [[45]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X8"    "X11"   "X20"   "X1:E" 
#>  [9] "X2:E"  "X3:E"  "X4:E"  "X8:E"  "X11:E" "E"    
#> 
#> [[46]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X8"    "X11"   "X20"   "X1:E" 
#>  [9] "X2:E"  "X3:E"  "X4:E"  "X8:E"  "X11:E" "E"    
#> 
#> [[47]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X8"    "X11"   "X20"   "X1:E" 
#>  [9] "X2:E"  "X3:E"  "X4:E"  "X8:E"  "X11:E" "E"    
#> 
#> [[48]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X8"    "X11"   "X20"   "X1:E" 
#>  [9] "X2:E"  "X3:E"  "X4:E"  "X8:E"  "X11:E" "E"    
#> 
#> [[49]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X8"    "X10"   "X11"   "X20"  
#>  [9] "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X8:E"  "X11:E" "E"    
#> 
#> [[50]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X8"    "X10"   "X11"   "X20"  
#>  [9] "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X8:E"  "X11:E" "E"    
#> 
#> [[51]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X6"    "X8"    "X10"   "X11"  
#>  [9] "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X8:E"  "X11:E" "E"    
#> 
#> [[52]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X6"    "X8"    "X10"   "X11"  
#>  [9] "X16"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X8:E"  "X11:E"
#> [17] "E"    
#> 
#> [[53]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X6"    "X8"    "X10"   "X11"  
#>  [9] "X15"   "X16"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E" 
#> [17] "X8:E"  "X11:E" "E"    
#> 
#> [[54]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X10"  
#>  [9] "X11"   "X15"   "X16"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E" 
#> [17] "X4:E"  "X8:E"  "X11:E" "E"    
#> 
#> [[55]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X10"  
#>  [9] "X11"   "X15"   "X16"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E" 
#> [17] "X4:E"  "X8:E"  "X11:E" "E"    
#> 
#> [[56]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X10"  
#>  [9] "X11"   "X15"   "X16"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E" 
#> [17] "X4:E"  "X8:E"  "X11:E" "E"    
#> 
#> [[57]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X10"  
#>  [9] "X11"   "X14"   "X15"   "X16"   "X19"   "X20"   "X1:E"  "X2:E" 
#> [17] "X3:E"  "X4:E"  "X8:E"  "X11:E" "E"    
#> 
#> [[58]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X10"  
#>  [9] "X11"   "X14"   "X15"   "X16"   "X19"   "X20"   "X1:E"  "X2:E" 
#> [17] "X3:E"  "X4:E"  "X8:E"  "X11:E" "E"    
#> 
#> [[59]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X10"  
#>  [9] "X11"   "X14"   "X15"   "X16"   "X17"   "X19"   "X20"   "X1:E" 
#> [17] "X2:E"  "X3:E"  "X4:E"  "X8:E"  "X11:E" "E"    
#> 
#> [[60]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X10"  
#>  [9] "X11"   "X14"   "X15"   "X16"   "X17"   "X19"   "X20"   "X1:E" 
#> [17] "X2:E"  "X3:E"  "X4:E"  "X8:E"  "X11:E" "X16:E" "E"    
#> 
#> [[61]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X9"   
#>  [9] "X10"   "X11"   "X15"   "X16"   "X20"   "X1:E"  "X2:E"  "X3:E" 
#> [17] "X4:E"  "X5:E"  "X6:E"  "X8:E"  "X9:E"  "X11:E" "X16:E" "E"    
#> 
#> [[62]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X9"   
#>  [9] "X10"   "X11"   "X15"   "X16"   "X20"   "X1:E"  "X2:E"  "X3:E" 
#> [17] "X4:E"  "X5:E"  "X6:E"  "X8:E"  "X9:E"  "X11:E" "X16:E" "E"    
#> 
#> [[63]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X9"   
#>  [9] "X10"   "X11"   "X15"   "X16"   "X20"   "X1:E"  "X2:E"  "X3:E" 
#> [17] "X4:E"  "X5:E"  "X6:E"  "X8:E"  "X9:E"  "X11:E" "X16:E" "E"    
#> 
#> [[64]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X9"   
#>  [9] "X10"   "X11"   "X12"   "X15"   "X16"   "X20"   "X1:E"  "X2:E" 
#> [17] "X3:E"  "X4:E"  "X5:E"  "X6:E"  "X8:E"  "X9:E"  "X11:E" "X16:E"
#> [25] "E"    
#> 
#> [[65]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X9"   
#>  [9] "X10"   "X11"   "X12"   "X15"   "X16"   "X20"   "X1:E"  "X2:E" 
#> [17] "X3:E"  "X4:E"  "X5:E"  "X6:E"  "X8:E"  "X9:E"  "X11:E" "X16:E"
#> [25] "E"    
#> 
#> [[66]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X9"   
#>  [9] "X10"   "X11"   "X12"   "X15"   "X16"   "X20"   "X1:E"  "X2:E" 
#> [17] "X3:E"  "X4:E"  "X5:E"  "X6:E"  "X8:E"  "X9:E"  "X11:E" "X16:E"
#> [25] "X20:E" "E"    
#> 
#> [[67]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X9"   
#>  [9] "X10"   "X11"   "X12"   "X15"   "X16"   "X20"   "X1:E"  "X2:E" 
#> [17] "X3:E"  "X4:E"  "X5:E"  "X6:E"  "X8:E"  "X9:E"  "X11:E" "X16:E"
#> [25] "X20:E" "E"    
#> 
#> [[68]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X9"   
#>  [9] "X10"   "X11"   "X12"   "X15"   "X16"   "X20"   "X1:E"  "X2:E" 
#> [17] "X3:E"  "X4:E"  "X5:E"  "X6:E"  "X8:E"  "X9:E"  "X11:E" "X16:E"
#> [25] "X20:E" "E"    
#> 
#> [[69]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X9"   
#>  [9] "X10"   "X11"   "X12"   "X15"   "X16"   "X20"   "X1:E"  "X2:E" 
#> [17] "X3:E"  "X4:E"  "X5:E"  "X6:E"  "X8:E"  "X9:E"  "X11:E" "X16:E"
#> [25] "X20:E" "E"    
#> 
#> [[70]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X9"   
#>  [9] "X10"   "X11"   "X12"   "X14"   "X15"   "X16"   "X18"   "X20"  
#> [17] "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E"  "X6:E"  "X8:E"  "X9:E" 
#> [25] "X11:E" "X12:E" "X14:E" "X15:E" "X16:E" "E"    
#> 
#> [[71]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X9"   
#>  [9] "X10"   "X11"   "X12"   "X14"   "X15"   "X16"   "X18"   "X20"  
#> [17] "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E"  "X6:E"  "X8:E"  "X9:E" 
#> [25] "X11:E" "X12:E" "X14:E" "X15:E" "X16:E" "E"    
#> 
#> [[72]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X9"   
#>  [9] "X10"   "X11"   "X12"   "X14"   "X15"   "X16"   "X18"   "X20"  
#> [17] "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E"  "X6:E"  "X8:E"  "X9:E" 
#> [25] "X11:E" "X12:E" "X14:E" "X15:E" "X16:E" "E"    
#> 
#> [[73]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X9"   
#>  [9] "X10"   "X11"   "X12"   "X14"   "X15"   "X16"   "X18"   "X20"  
#> [17] "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E"  "X6:E"  "X8:E"  "X9:E" 
#> [25] "X11:E" "X12:E" "X14:E" "X15:E" "X16:E" "E"    
#> 
#> [[74]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X9"   
#>  [9] "X10"   "X11"   "X12"   "X14"   "X15"   "X16"   "X18"   "X20"  
#> [17] "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E"  "X6:E"  "X8:E"  "X9:E" 
#> [25] "X11:E" "X12:E" "X14:E" "X15:E" "X16:E" "E"    
#> 
#> [[75]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X8"    "X9"   
#>  [9] "X10"   "X11"   "X12"   "X14"   "X15"   "X16"   "X18"   "X20"  
#> [17] "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E"  "X6:E"  "X8:E"  "X9:E" 
#> [25] "X11:E" "X12:E" "X14:E" "X15:E" "X16:E" "E"    
#> 
#> [[76]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X14"   "X15"   "X16"   "X18"  
#> [17] "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E"  "X6:E"  "X8:E" 
#> [25] "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E" "X16:E" "E"    
#> 
#> [[77]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X14"   "X15"   "X16"   "X18"  
#> [17] "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E"  "X6:E"  "X8:E" 
#> [25] "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E" "X16:E" "E"    
#> 
#> [[78]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X14"   "X15"   "X16"   "X18"  
#> [17] "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E"  "X6:E"  "X8:E" 
#> [25] "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E" "X16:E" "X20:E" "E"    
#> 
#> [[79]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X14"   "X15"   "X16"   "X18"  
#> [17] "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E"  "X6:E"  "X8:E" 
#> [25] "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E" "X16:E" "X20:E" "E"    
#> 
#> [[80]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X14"   "X15"   "X16"   "X18"  
#> [17] "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E"  "X6:E" 
#> [25] "X8:E"  "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E" "X16:E" "X20:E"
#> [33] "E"    
#> 
#> [[81]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X14"   "X15"   "X16"   "X18"  
#> [17] "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E"  "X6:E" 
#> [25] "X8:E"  "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E" "X16:E" "X20:E"
#> [33] "E"    
#> 
#> [[82]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X14"   "X15"   "X16"   "X18"  
#> [17] "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E"  "X6:E" 
#> [25] "X8:E"  "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E" "X16:E" "X20:E"
#> [33] "E"    
#> 
#> [[83]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X14"   "X15"   "X16"   "X18"  
#> [17] "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E"  "X6:E" 
#> [25] "X8:E"  "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E" "X16:E" "X20:E"
#> [33] "E"    
#> 
#> [[84]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X14"   "X15"   "X16"   "X18"  
#> [17] "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E"  "X6:E" 
#> [25] "X8:E"  "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E" "X16:E" "X20:E"
#> [33] "E"    
#> 
#> [[85]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X14"   "X15"   "X16"   "X18"  
#> [17] "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E"  "X6:E" 
#> [25] "X8:E"  "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E" "X16:E" "X20:E"
#> [33] "E"    
#> 
#> [[86]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X13"   "X14"   "X15"   "X16"  
#> [17] "X18"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E" 
#> [25] "X6:E"  "X8:E"  "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E" "X16:E"
#> [33] "X20:E" "E"    
#> 
#> [[87]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X13"   "X14"   "X15"   "X16"  
#> [17] "X18"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E" 
#> [25] "X6:E"  "X7:E"  "X8:E"  "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E"
#> [33] "X16:E" "X20:E" "E"    
#> 
#> [[88]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X13"   "X14"   "X15"   "X16"  
#> [17] "X18"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E" 
#> [25] "X6:E"  "X7:E"  "X8:E"  "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E"
#> [33] "X16:E" "X20:E" "E"    
#> 
#> [[89]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X13"   "X14"   "X15"   "X16"  
#> [17] "X18"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E" 
#> [25] "X6:E"  "X7:E"  "X8:E"  "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E"
#> [33] "X16:E" "X20:E" "E"    
#> 
#> [[90]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X13"   "X14"   "X15"   "X16"  
#> [17] "X18"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E" 
#> [25] "X6:E"  "X7:E"  "X8:E"  "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E"
#> [33] "X16:E" "X20:E" "E"    
#> 
#> [[91]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X13"   "X14"   "X15"   "X16"  
#> [17] "X18"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E" 
#> [25] "X6:E"  "X7:E"  "X8:E"  "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E"
#> [33] "X16:E" "X20:E" "E"    
#> 
#> [[92]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X13"   "X14"   "X15"   "X16"  
#> [17] "X18"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E" 
#> [25] "X6:E"  "X7:E"  "X8:E"  "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E"
#> [33] "X16:E" "X20:E" "E"    
#> 
#> [[93]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X13"   "X14"   "X15"   "X16"  
#> [17] "X18"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E" 
#> [25] "X6:E"  "X7:E"  "X8:E"  "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E"
#> [33] "X16:E" "X20:E" "E"    
#> 
#> [[94]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X13"   "X14"   "X15"   "X16"  
#> [17] "X18"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E"  "X5:E" 
#> [25] "X6:E"  "X7:E"  "X8:E"  "X9:E"  "X11:E" "X12:E" "X14:E" "X15:E"
#> [33] "X16:E" "X20:E" "E"    
#> 
#> [[95]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X13"   "X14"   "X15"   "X16"  
#> [17] "X17"   "X18"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E" 
#> [25] "X5:E"  "X6:E"  "X7:E"  "X8:E"  "X9:E"  "X10:E" "X11:E" "X12:E"
#> [33] "X13:E" "X14:E" "X15:E" "X16:E" "X18:E" "X19:E" "X20:E" "E"    
#> 
#> [[96]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X13"   "X14"   "X15"   "X16"  
#> [17] "X17"   "X18"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E" 
#> [25] "X5:E"  "X6:E"  "X7:E"  "X8:E"  "X9:E"  "X10:E" "X11:E" "X12:E"
#> [33] "X13:E" "X14:E" "X15:E" "X16:E" "X18:E" "X19:E" "X20:E" "E"    
#> 
#> [[97]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X13"   "X14"   "X15"   "X16"  
#> [17] "X17"   "X18"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E" 
#> [25] "X5:E"  "X6:E"  "X7:E"  "X8:E"  "X9:E"  "X10:E" "X11:E" "X12:E"
#> [33] "X13:E" "X14:E" "X15:E" "X16:E" "X18:E" "X19:E" "X20:E" "E"    
#> 
#> [[98]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X13"   "X14"   "X15"   "X16"  
#> [17] "X17"   "X18"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E" 
#> [25] "X5:E"  "X6:E"  "X7:E"  "X8:E"  "X9:E"  "X10:E" "X11:E" "X12:E"
#> [33] "X13:E" "X14:E" "X15:E" "X16:E" "X18:E" "X19:E" "X20:E" "E"    
#> 
#> [[99]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X13"   "X14"   "X15"   "X16"  
#> [17] "X17"   "X18"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E" 
#> [25] "X5:E"  "X6:E"  "X7:E"  "X8:E"  "X9:E"  "X10:E" "X11:E" "X12:E"
#> [33] "X13:E" "X14:E" "X15:E" "X16:E" "X18:E" "X19:E" "X20:E" "E"    
#> 
#> [[100]]
#>  [1] "X1"    "X2"    "X3"    "X4"    "X5"    "X6"    "X7"    "X8"   
#>  [9] "X9"    "X10"   "X11"   "X12"   "X13"   "X14"   "X15"   "X16"  
#> [17] "X17"   "X18"   "X19"   "X20"   "X1:E"  "X2:E"  "X3:E"  "X4:E" 
#> [25] "X5:E"  "X6:E"  "X7:E"  "X8:E"  "X9:E"  "X10:E" "X11:E" "X12:E"
#> [33] "X13:E" "X14:E" "X15:E" "X16:E" "X18:E" "X19:E" "X20:E" "E"

Cross-Validation

cv.sail is the main function to do cross-validation along with plot, predict, and coef methods for objects of class cv.sail. We run it in parallel:

set.seed(432) # to reproduce results (randomness due to CV folds)
library(doMC)
#> Loading required package: foreach
#> Loading required package: iterators
#> Loading required package: parallel
registerDoMC(cores = 2) 
cvfit <- cv.sail(x = sailsim$x, y = sailsim$y, e = sailsim$e, basis = f.basis,
                 nfolds = 5, parallel = TRUE)

We plot the cross-validated error curve which has the mean-squared error on the y-axis and \(\log(\lambda)\) on the x-axis. It includes the cross-validation curve (red dotted line), and upper and lower standard deviation curves along the \(\lambda\) sequence (error bars). Two selected \(\lambda\)’s are indicated by the vertical dotted lines (see below). The numbers at the top of the plot represent the total number of non-zero variables at that value of \(\lambda\) (df_main + df_environment + df_interaction):

plot(cvfit)

lambda.min is the value of \(\lambda\) that gives minimum mean cross-validated error. The other \(\lambda\) saved is lambda.1se, which gives the most regularized model such that error is within one standard error of the minimum. We can view the selected \(\lambda\)’s and the corresponding coefficients:

cvfit[["lambda.min"]]
#> [1] 0.2788543
cvfit[["lambda.1se"]]
#> [1] 0.4238337

The estimated coefficients at lambda.1se and lambda.min:

cbind(coef(cvfit, s="lambda.1se"), # lambda.1se is the default
coef(cvfit, s = "lambda.min"))
#> 202 x 2 sparse Matrix of class "dgCMatrix"
#>                        1            1
#> (Intercept)   5.42162243   5.43619701
#> X1_1         -0.78253656  -1.23830516
#> X1_2          0.13429262   0.48732240
#> X1_3          0.43865579   0.88911732
#> X1_4          0.79446840   1.10810480
#> X1_5          1.92178800   2.68114985
#> X2_1          .           -0.21979939
#> X2_2          .           -0.95046682
#> X2_3          .           -0.76046594
#> X2_4          .           -0.17408785
#> X2_5          .           -0.06842359
#> X3_1          3.45105131   4.42740541
#> X3_2          1.56761411   2.28691047
#> X3_3         -1.27174989  -1.36817236
#> X3_4         -2.27784639  -2.56642476
#> X3_5         -1.23352658  -1.03097615
#> X4_1          4.47171678   5.37947206
#> X4_2         -2.87477819  -3.35562770
#> X4_3         -6.65050142  -8.07333662
#> X4_4         -3.44081136  -3.50308509
#> X4_5         -0.47029997  -0.55029116
#> X5_1          .            .         
#> X5_2          .            .         
#> X5_3          .            .         
#> X5_4          .            .         
#> X5_5          .            .         
#> X6_1          .            .         
#> X6_2          .            .         
#> X6_3          .            .         
#> X6_4          .            .         
#> X6_5          .            .         
#> X7_1          .            .         
#> X7_2          .            .         
#> X7_3          .            .         
#> X7_4          .            .         
#> X7_5          .            .         
#> X8_1          0.17501411   0.47920984
#> X8_2          0.11615654   0.33988638
#> X8_3          0.01782976   0.07650873
#> X8_4         -0.10502457  -0.28112125
#> X8_5         -0.24477578  -0.66610209
#> X9_1          .            .         
#> X9_2          .            .         
#> X9_3          .            .         
#> X9_4          .            .         
#> X9_5          .            .         
#> X10_1         .            .         
#> X10_2         .            .         
#> X10_3         .            .         
#> X10_4         .            .         
#> X10_5         .            .         
#> X11_1         0.13685648   0.34767475
#> X11_2         0.02857112   0.06273705
#> X11_3        -0.07321859  -0.20148715
#> X11_4        -0.20935470  -0.67854244
#> X11_5        -0.22638940  -0.69583368
#> X12_1         .            .         
#> X12_2         .            .         
#> X12_3         .            .         
#> X12_4         .            .         
#> X12_5         .            .         
#> X13_1         .            .         
#> X13_2         .            .         
#> X13_3         .            .         
#> X13_4         .            .         
#> X13_5         .            .         
#> X14_1         .            .         
#> X14_2         .            .         
#> X14_3         .            .         
#> X14_4         .            .         
#> X14_5         .            .         
#> X15_1         .            .         
#> X15_2         .            .         
#> X15_3         .            .         
#> X15_4         .            .         
#> X15_5         .            .         
#> X16_1         .            .         
#> X16_2         .            .         
#> X16_3         .            .         
#> X16_4         .            .         
#> X16_5         .            .         
#> X17_1         .            .         
#> X17_2         .            .         
#> X17_3         .            .         
#> X17_4         .            .         
#> X17_5         .            .         
#> X18_1         .            .         
#> X18_2         .            .         
#> X18_3         .            .         
#> X18_4         .            .         
#> X18_5         .            .         
#> X19_1         .            .         
#> X19_2         .            .         
#> X19_3         .            .         
#> X19_4         .            .         
#> X19_5         .            .         
#> X20_1         .            0.10215979
#> X20_2         .           -0.03829014
#> X20_3         .           -0.17766509
#> X20_4         .           -0.18547953
#> X20_5         .            0.09733455
#> E             1.99574559   2.06409655
#> X1_1:E        .           -0.53455653
#> X1_2:E        .            0.21036928
#> X1_3:E        .            0.38381772
#> X1_4:E        .            0.47835111
#> X1_5:E        .            1.15740950
#> X2_1:E        .            .         
#> X2_2:E        .            .         
#> X2_3:E        .            .         
#> X2_4:E        .            .         
#> X2_5:E        .            .         
#> X3_1:E        1.80698510   2.37991172
#> X3_2:E        0.82080940   1.22930803
#> X3_3:E       -0.66589363  -0.73544867
#> X3_4:E       -1.19269003  -1.37955841
#> X3_5:E       -0.64587975  -0.55419190
#> X4_1:E        8.31246506   8.87335933
#> X4_2:E       -5.34391922  -5.53505808
#> X4_3:E      -12.36260331 -13.31684893
#> X4_4:E       -6.39611711  -5.77828687
#> X4_5:E       -0.87423964  -0.90769710
#> X5_1:E        .            .         
#> X5_2:E        .            .         
#> X5_3:E        .            .         
#> X5_4:E        .            .         
#> X5_5:E        .            .         
#> X6_1:E        .            .         
#> X6_2:E        .            .         
#> X6_3:E        .            .         
#> X6_4:E        .            .         
#> X6_5:E        .            .         
#> X7_1:E        .            .         
#> X7_2:E        .            .         
#> X7_3:E        .            .         
#> X7_4:E        .            .         
#> X7_5:E        .            .         
#> X8_1:E        .            .         
#> X8_2:E        .            .         
#> X8_3:E        .            .         
#> X8_4:E        .            .         
#> X8_5:E        .            .         
#> X9_1:E        .            .         
#> X9_2:E        .            .         
#> X9_3:E        .            .         
#> X9_4:E        .            .         
#> X9_5:E        .            .         
#> X10_1:E       .            .         
#> X10_2:E       .            .         
#> X10_3:E       .            .         
#> X10_4:E       .            .         
#> X10_5:E       .            .         
#> X11_1:E       .            .         
#> X11_2:E       .            .         
#> X11_3:E       .            .         
#> X11_4:E       .            .         
#> X11_5:E       .            .         
#> X12_1:E       .            .         
#> X12_2:E       .            .         
#> X12_3:E       .            .         
#> X12_4:E       .            .         
#> X12_5:E       .            .         
#> X13_1:E       .            .         
#> X13_2:E       .            .         
#> X13_3:E       .            .         
#> X13_4:E       .            .         
#> X13_5:E       .            .         
#> X14_1:E       .            .         
#> X14_2:E       .            .         
#> X14_3:E       .            .         
#> X14_4:E       .            .         
#> X14_5:E       .            .         
#> X15_1:E       .            .         
#> X15_2:E       .            .         
#> X15_3:E       .            .         
#> X15_4:E       .            .         
#> X15_5:E       .            .         
#> X16_1:E       .            .         
#> X16_2:E       .            .         
#> X16_3:E       .            .         
#> X16_4:E       .            .         
#> X16_5:E       .            .         
#> X17_1:E       .            .         
#> X17_2:E       .            .         
#> X17_3:E       .            .         
#> X17_4:E       .            .         
#> X17_5:E       .            .         
#> X18_1:E       .            .         
#> X18_2:E       .            .         
#> X18_3:E       .            .         
#> X18_4:E       .            .         
#> X18_5:E       .            .         
#> X19_1:E       .            .         
#> X19_2:E       .            .         
#> X19_3:E       .            .         
#> X19_4:E       .            .         
#> X19_5:E       .            .         
#> X20_1:E       .            .         
#> X20_2:E       .            .         
#> X20_3:E       .            .         
#> X20_4:E       .            .         
#> X20_5:E       .            .

Estimated non-zero coefficients at lambda.1se:

predict(cvfit, type = "nonzero")
#>                        1
#> (Intercept)   5.42162243
#> X1_1         -0.78253656
#> X1_2          0.13429262
#> X1_3          0.43865579
#> X1_4          0.79446840
#> X1_5          1.92178800
#> X3_1          3.45105131
#> X3_2          1.56761411
#> X3_3         -1.27174989
#> X3_4         -2.27784639
#> X3_5         -1.23352658
#> X4_1          4.47171678
#> X4_2         -2.87477819
#> X4_3         -6.65050142
#> X4_4         -3.44081136
#> X4_5         -0.47029997
#> X8_1          0.17501411
#> X8_2          0.11615654
#> X8_3          0.01782976
#> X8_4         -0.10502457
#> X8_5         -0.24477578
#> X11_1         0.13685648
#> X11_2         0.02857112
#> X11_3        -0.07321859
#> X11_4        -0.20935470
#> X11_5        -0.22638940
#> E             1.99574559
#> X3_1:E        1.80698510
#> X3_2:E        0.82080940
#> X3_3:E       -0.66589363
#> X3_4:E       -1.19269003
#> X3_5:E       -0.64587975
#> X4_1:E        8.31246506
#> X4_2:E       -5.34391922
#> X4_3:E      -12.36260331
#> X4_4:E       -6.39611711
#> X4_5:E       -0.87423964

Visualizing the Effect of the Non-linear Terms

bsplines are difficult to interpret. We provide a plotting function to visualize the effect of the non-linear function on the response.

Main Effects

Since we are using simulated data, we also plot the true curve:

plotMain(cvfit$sail.fit, x = sailsim$x, xvar = "X3",
         legend.position = "topright",
         s = cvfit$lambda.min, f.truth = sailsim$f3)

Interaction Effects

Again, since we are using simulated data, we also plot the true interaction:

plotInter(cvfit$sail.fit, x = sailsim$x, xvar = "X4",
          f.truth = sailsim$f4.inter,
          s = cvfit$lambda.min,
          title_z = "Estimated")

Linear Interactions

The basis argument in the sail function is very flexible in that it allows you to apply any basis expansion to the columns of \(\mathbf{X}\). Of course, there might be situations where you do not expect any non-linear main effects or interactions to be present in your data. You can still use the sail method to search for linear main effects and interactions. This can be accomplished by specifying an identity map:

f.identity <- function(i) i

We then pass this function to basis argument in cv.sail:

cvfit_linear <- cv.sail(x = sailsim$x, y = sailsim$y, e = sailsim$e,
                        basis = f.identity, nfolds = 5, parallel = TRUE)

Next we plot the cross-validated curve:

plot(cvfit_linear)

And extract the model at lambda.min:

coef(cvfit_linear, s = "lambda.min")
#> 42 x 1 sparse Matrix of class "dgCMatrix"
#>                       1
#> (Intercept)   5.4385233
#> X1_1          2.4568905
#> X2_1         -1.6516411
#> X3_1         -5.1515618
#> X4_1         -7.1355682
#> X5_1          .        
#> X6_1          .        
#> X7_1         -1.1055226
#> X8_1         -0.8602264
#> X9_1          .        
#> X10_1         .        
#> X11_1        -2.1990113
#> X12_1         .        
#> X13_1         .        
#> X14_1        -1.9208173
#> X15_1         .        
#> X16_1         3.2004658
#> X17_1         .        
#> X18_1         1.0234816
#> X19_1         .        
#> X20_1        -0.2161608
#> E             2.4466054
#> X1_1:E       -4.3993869
#> X2_1:E       -2.7114248
#> X3_1:E       -4.7196594
#> X4_1:E      -12.9774650
#> X5_1:E        .        
#> X6_1:E        .        
#> X7_1:E        2.3703946
#> X8_1:E        .        
#> X9_1:E        .        
#> X10_1:E       .        
#> X11_1:E      -2.1462173
#> X12_1:E       .        
#> X13_1:E       .        
#> X14_1:E      -0.9787678
#> X15_1:E       .        
#> X16_1:E       6.9265731
#> X17_1:E       .        
#> X18_1:E       1.9004167
#> X19_1:E       .        
#> X20_1:E       .

Applying a different penalty to each predictor

Recall that we consider the following penalized least squares criterion for this problem:

\[\begin{equation} \arg\min_{\boldsymbol{\theta} } \mathcal{L}(Y;\boldsymbol{\theta}) + \lambda (1-\alpha) \left( w_E |\beta_E| + \sum_{j=1}^{p} w_j \lVert\boldsymbol{\theta}_j \rVert_2 \right) + \lambda\alpha \sum_{j=1}^{p} w_{jE} |\gamma_{j}| \end{equation}\]

The weights \(w_E, w_j, w_{jE}\) are by default set to 1 as specified by the penalty.factor argument. This argument allows users to apply separate penalty factors to each coefficient. In particular, any variable with penalty.factor equal to zero is not penalized at all. This feature can be applied mainly for two reasons:

  1. Prior knowledge about the importance of certain variables is known. Larger weights will penalize the variable more, while smaller weights will penalize the variable less
  2. Allows users to apply the Adaptive sail, similar to the Adaptive Lasso

In the following example, we want the environment variable to always be included so we set the first element of p.fac to zero. We also want to apply less of a penalty to the main effects for \(X_2, X_3, X_4\):

# the weights correspond to E, X1, X2, X3, ... X_p, X1:E, X2:E, ... X_p:E
p.fac <- c(0, 1, 0.4, 0.6, 0.7, rep(1, 2*ncol(sailsim$x) - 4))
fit_pf <- sail(x = sailsim$x, y = sailsim$y, e = sailsim$e, basis = f.basis,
               penalty.factor = p.fac)
plot(fit_pf)

We see from the plot above that the black line (corresponding to the \(E\) variable with penalty.factor equal to zero) is always included in the model.