Vignette of the pengls package

Stijn Hawinkel

1 Introduction

This vignette demonstrates the use of the pengls package for high-dimensional data with spatial or temporal autocorrelation. It consists of an iterative loop around the nlme and glmnet packages. Currently, only continuous outcomes and \(R^2\) as performance measure are implemented.

2 Installation instuctions

The pengls package is available from BioConductor, and can be installed as follows:

library(BiocManager)
install("pengls")

Once installed, it can be loaded and version info printed.

suppressPackageStartupMessages(library(pengls))
cat("pengls package version", as.character(packageVersion("pengls")), "\n")
## pengls package version 1.3.0

3 Illustration

3.1 Spatial autocorrelation

We first create a toy dataset with spatial coordinates.

The pengls method requires prespecification of a functional form for the autocorrelation. This is done through the corStruct objects defined by the nlme package. We specify a correlation decaying as a Gaussian curve with distance, and with a nugget parameter. The nugget parameter is a proportion that indicates how much of the correlation structure explained by independent errors; the rest is attributed to spatial autocorrelation. The starting values are chosen as reasonable guesses; they will be overwritten in the fitting process.

Finally the model is fitted with a single outcome variable and large number of regressors, with the chosen covariance structure and for a prespecified penalty parameter \(\lambda=0.2\).

## Starting iterations...
## Iteration 1 
## Iteration 2 
## Iteration 3

Standard extraction functions like print(), coef() and predict() are defined for the new “pengls” object.

## pengls model with correlation structure: corGaus 
##  and 20 non-zero coefficients

3.2 Temporal autocorrelation

The method can also account for temporal autocorrelation by defining another correlation structure from the nlme package, e.g. autocorrelation structure of order 1:

The fitting command is similar, this time the \(\lambda\) parameter is found through cross-validation of the naive glmnet (for full cross-validation , see below). We choose \(\alpha=0.5\) this time, fitting an elastic net model.

## Fitting naieve model...
## Starting iterations...
## Iteration 1 
## Iteration 2

Show the output

## pengls model with correlation structure: corAR1 
##  and 3 non-zero coefficients

3.3 Penalty parameter and cross-validation

The pengls package also provides cross-validation for finding the optimal \(\lambda\) value. If the tuning parameter \(\lambda\) is not supplied, the optimal \(\lambda\) according to cross-validation with the naive glmnet function (the one that ignores dependence) is used. Hence we recommend to use the following function to use cross-validation. Multithreading is supported through the BiocParallel package :

The function is called similarly to cv.glmnet:

Check the result:

## Cross-validated pengls model with correlation structure: corGaus 
##  and 0 non-zero coefficients.
##  3 fold cross-validation yielded an estimated R2 of -5.274116 .

By default, the 1 standard error is used to determine the optimal value of \(\lambda\) :

## [1] 2.397772
## [1] -5.274116

Extract coefficients and fold IDs.

## [1] 0 0 0 0 0 0
##  56  80 201  37  49  99 111 121 123 119 158   8 152 192 159  28 204 133   7 113 
##   3   1   1   2   3   3   1   1   3   2   1   2   1   1   1   2   1   3   3   3 
##  55 176 162 122 100 
##   3   3   2   1   2

By default, blocked cross-validation is used, but random cross-validation is also available (but not recommended for timecourse or spatial data). First we illustrate the different ways graphically, again using the timecourse example:

To perform random cross-validation

To negate baseline differences at different timepoints, it may be useful to center or scale the outcomes in the cross validation. For instance for centering only:

## [1] 0.9959231

4 Session info

sessionInfo()
## R version 4.2.0 RC (2022-04-21 r82226)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] BiocParallel_1.31.0 nlme_3.1-157        pengls_1.3.0       
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.8.3     knitr_1.38       magrittr_2.0.3   splines_4.2.0   
##  [5] lattice_0.20-45  R6_2.5.1         rlang_1.0.2      fastmap_1.1.0   
##  [9] foreach_1.5.2    highr_0.9        stringr_1.4.0    tools_4.2.0     
## [13] parallel_4.2.0   grid_4.2.0       glmnet_4.1-4     xfun_0.30       
## [17] cli_3.3.0        jquerylib_0.1.4  htmltools_0.5.2  iterators_1.0.14
## [21] survival_3.3-1   yaml_2.3.5       digest_0.6.29    Matrix_1.4-1    
## [25] sass_0.4.1       codetools_0.2-18 shape_1.4.6      evaluate_0.15   
## [29] rmarkdown_2.14   stringi_1.7.6    compiler_4.2.0   bslib_0.3.1     
## [33] jsonlite_1.8.0