ktweedie: Kernel-based Tweedie compound Poisson gamma model using high-dimensional covariates for the analyses of zero-inflated response variables.

Introduction

ktweedie is a package that fits nonparametric Tweedie compound Poisson gamma models in the reproducing kernel Hilbert space. The package is based on two algorithms, the ktweedie for kernel-based Tweedie model and the sktweedie for sparse kernel-based Tweedie model. The ktweedie supports a wide range of kernel functions implemented in the R package kernlab and the optimization algorithm is a Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm with bisection line search. The package includes cross-validation functions for one-dimensional tuning of the kernel regularization parameter and for two-dimensional joint tuning of one kernel parameter and . The sktweedie uses variable weights to achieve variable selection. It is a meta-algorithm that alternatively updates the kernel parameters and a set of variable weights.

The ktweedie solves the problem

where is the index parameter, is the dispersion parameter, is an kernel matrix computed according to the user-specified kernel function , whose entries are are calculated based on the -dimensional predictors from subjects . In the kernel-based Tweedie model, the mean parameter for the -th observation is modelled by

The sktweedie solves

where involves variable weights .

Installation

1. From the CRAN.

install.packages("ktweedie")

2. From the Github.

devtools::install_github("ly129/ktweedie")

Quick Start

First we load the ktweedie package:

library(ktweedie)

The package includes a toy data for demonstration purpose. The predictor matrix x is generated from standard normal distribution and y is generated according to

where . That said, only the first two predictors are associated with the response.

data(dat)
x <- dat$x
y <- dat$y

An input matrix x and an output vector y are now loaded. The ktd_estimate() function can be used to fit a basic ktweedie model. The regularization parameter lam1 can be a vector, which will be solved in a decreasing order with warm start.

fit.ktd <- ktd_estimate(x = x,
                        y = y,
                        kern = rbfdot(sigma = 0.1),
                        lam1 = c(0.01, 0.1, 1))
str(fit.ktd$estimates)
#> List of 3
#>  $ lambda 1   :List of 3
#>   ..$ fn         : num 110
#>   ..$ coefficient: num [1:30, 1] 0.5558 -0.062 -0.0381 0.0523 -0.0251 ...
#>   ..$ convergence: int 0
#>  $ lambda 0.1 :List of 3
#>   ..$ fn         : num 51
#>   ..$ coefficient: num [1:30, 1] 1.662 -0.235 -0.177 0.867 -0.143 ...
#>   ..$ convergence: int 0
#>  $ lambda 0.01:List of 3
#>   ..$ fn         : num 39.2
#>   ..$ coefficient: num [1:30, 1] 7.692 -0.49 -0.841 4.624 -0.696 ...
#>   ..$ convergence: int 0

fit.ktd$estimates stores the estimated coefficients and the final objective function value. The estimated kernel-based model coefficients for the -th lam1 can be accessed by the index l: fit.ktd$estimates[[l]]$coefficient.

The function can also be used to fit the sktweedie model by setting the argument sparsity to TRUE, and specifying the regularization coefficient in the argument lam2.

fit.sktd <- ktd_estimate(x = x,
                         y = y,
                         kern = rbfdot(sigma = 0.1),
                         lam1 = 5,
                         sparsity = TRUE,
                         lam2 = 1)

And we can access the fitted coefficients in a similar manner to the fit.ktd. Additionally, the fitted variable weights can be accessed by

fit.sktd$estimates[[1]]$weight
#>           [,1]
#> [1,] 1.0000000
#> [2,] 0.4462078
#> [3,] 0.0000000
#> [4,] 0.0000000
#> [5,] 0.0000000

Variables with weights close to 0 can be viewed as noise variables.

Recommended Data Analysis Pipeline

The ktweedie and sktweedie algorithms require careful tuning of one to multiple hyperparameters, depending on the choice of kernel functions. For the ktweedie, we recommend either a one-dimensional tuning for lam1 () or a two-dimensional random search for lam1 and the kernel parameter using cross-validation. Tuning is achieved by optimizing a user-specified criterion, including log likelihood loss = "LL", mean absolute difference loss = "MAD" and root mean squared error loss = "RMSE". Using the Laplacian kernel as an example.

laplacedot(sigma = 1)
#> Laplace kernel function. 
#>  Hyperparameter : sigma =  1

Cross-validation

The one-dimensional search for the optimal lam1, can be achieved with the ktd_cv() function from a user-specified vector of candidate values:

ktd.cv1d <- ktd_cv(x = x,
                   y = y,
                   kern = laplacedot(sigma = 0.1),
                   lambda = c(0.0001, 0.001, 0.01, 0.1, 1),
                   nfolds = 5,
                   loss = "LL")
ktd.cv1d
#> $LL
#>         1       0.1      0.01     0.001     1e-04 
#> -82.30040 -60.33054 -55.68177 -55.68835 -65.38823 
#> 
#> $Best_lambda
#> [1] 0.01

The two-dimensional joint search for the optimal lam1 and sigma requires ktd_cv2d(). In the example below, a total of ncoefs = 10 pairs of candidate lam1 and sigma values are randomly sampled (uniformly on the log scale) within the ranges lambda = c(1e-5, 1e0) and sigma = c(1e-5, 1e0), respectively. Then the nfolds = 5-fold cross-validation is performed to select the pair that generates the optimal cross-validation loss = "MAD".

ktd.cv2d <- ktd_cv2d(x = x,
                     y = y,
                     kernfunc = laplacedot,
                     lambda = c(1e-5, 1e0),
                     sigma = c(1e-5, 1e0),
                     nfolds = 5,
                     ncoefs = 10,
                     loss = "MAD")
ktd.cv2d
#> $MAD
#>    Lambda=0.000435692, Sigma=0.174196   Lambda=0.00855899, Sigma=0.00201436 
#>                              354.1993                              431.4734 
#>  Lambda=0.00518177, Sigma=0.000749782   Lambda=7.25693e-05, Sigma=0.0620986 
#>                              469.7289                              327.0395 
#>   Lambda=0.0513091, Sigma=0.000344321   Lambda=0.0108477, Sigma=0.000277883 
#>                              626.3884                              589.4097 
#> Lambda=9.72691e-05, Sigma=2.19179e-05   Lambda=0.0682224, Sigma=0.000455657 
#>                              433.5755                              624.1514 
#>   Lambda=0.000228745, Sigma=0.0247239     Lambda=0.166265, Sigma=0.00695988 
#>                              332.0113                              544.0900 
#> 
#> $Best_lambda
#> [1] 7.25693e-05
#> 
#> $Best_sigma
#> [1] 0.0620986

Fitting

Then the model is fitted using the hyperparameter(s) selected by the ktd_cv() or ktd_cv2d(). In the example below, the selected lam1 and sigma values are accessed by ktd.cv2d$Best_lambda and ktd.cv2d$Best_sigma, which are then be fed into the ktd_estimate() to perform final model fitting.

ktd.fit <- ktd_estimate(x = x,
                        y = y,
                        kern = laplacedot(sigma = ktd.cv2d$Best_sigma),
                        lam1 = ktd.cv2d$Best_lambda)
str(ktd.fit$estimates)
#> List of 1
#>  $ lambda 7.25693e-05:List of 3
#>   ..$ fn         : num 36.6
#>   ..$ coefficient: num [1:30, 1] 24.82 -9.63 -17.4 44.79 3.7 ...
#>   ..$ convergence: int 0

For the sktweedie, only the Gaussian radial basis function (RBF) kernel rbfdot() is supported. We recommend using the same set of tuned parameters as if a ktweedie model is fitted and tuning lam2 manually:

sktd.cv2d <- ktd_cv2d(x = x,
                      y = y,
                      kernfunc = rbfdot,
                      lambda = c(1e-3, 1e0),
                      sigma = c(1e-3, 1e0),
                      nfolds = 5,
                      ncoefs = 10,
                      loss = "LL")

sktd.fit <- ktd_estimate(x = x,
                         y = y,
                         kern = rbfdot(sigma = sktd.cv2d$Best_sigma),
                         lam1 = sktd.cv2d$Best_lambda,
                         sparsity = TRUE,
                         lam2 = 1,
                         ftol = 1e-3,
                         partol = 1e-3,
                         innerpartol = 1e-5)

Prediction

The function ktd_predict() can identify necessary information stored in ktd.fit$data and sktd.fit$data to make predictions at the user-specified newdata. If the argument newdata is unspecified, the prediction will be made at the original x used in model training and fitting.

ktd.pred <- ktd_predict(ktd.fit, type = "response")
head(ktd.pred$prediction)
#>              [,1]
#> [1,] 6.448220e+02
#> [2,] 1.750695e-03
#> [3,] 9.215399e-02
#> [4,] 4.713962e+00
#> [5,] 1.678452e-01
#> [6,] 1.650646e+00

If newdata with the same dimension as x is provided, the prediction will be made at the new data points.

# Use a subset of the original x as newdata.
newdata <- x[1:6, ]
ktd.pred.new <- ktd_predict(ktd.fit,
                            newdata = newdata,
                            type = "response")
sktd.pred.new <- ktd_predict(sktd.fit,
                             newdata = newdata,
                             type = "response")
data.frame(ktweedie = ktd.pred.new$prediction,
           sktweedie = sktd.pred.new$prediction)
#>       ktweedie  sktweedie
#> 1 6.448220e+02 421.931421
#> 2 1.750695e-03  22.543092
#> 3 9.215399e-02  23.415272
#> 4 4.713962e+00   1.642355
#> 5 1.678452e-01  12.034229
#> 6 1.650646e+00 122.187222

Variable Selection

In practice, the variable selection results of the sktweedie is more meaningful. An effective way to fit the sktweedie is to start with an arbitrarily big lam2 that sets all weights to zero and gradually decrease its value. Note that a warning message is generated for the first lam2, suggesting that all weights are set to zero.

nlam2 <- 10
lam2.seq <- 20 * 0.8^(1:nlam2 - 1)
wts <- matrix(NA, nrow = nlam2, ncol = ncol(x))
for (i in 1:nlam2) {
  sktd.tmp <- ktd_estimate(x = x,
                           y = y,
                           kern = rbfdot(sigma = sktd.cv2d$Best_sigma),
                           lam1 = sktd.cv2d$Best_lambda,
                           sparsity = TRUE,
                           lam2 = lam2.seq[i],
                           ftol = 1e-3,
                           partol = 1e-3,
                           innerpartol = 1e-5)
  wt.tmp <- sktd.tmp$estimates[[1]]$weight
  if (is.null(wt.tmp)) wts[i, ] <- 0 else wts[i, ] <- wt.tmp
}
#> WARNING: All weights are zero in weight update iteration:
#> [1] 2
# plot the solution path with graphics::matplot()
matplot(y = wts,
        x = lam2.seq,
        type = "l",
        log = "x",
        ylab = "Weights",
        xlab = expression(paste(lambda)),
        lwd = 2)
legend("topright",
       title = "w index",
       legend = 1:5,
       lty = 1:5,
       col = 1:6,
       lwd = 2)