Isotonic Regression

Introduction

Isotonic regression is regression with monotonic constraints. There are several packages in R to fit isotonic regression models. In this example, we consider isotone which uses a pooled-adjacent-violators algorithm (PAVA) and active set methods to perform the fit.

Pituitary Data Example

We will use data from the isotone package (see de Leeuw, Hornik, and Mair (2009)) on the size of pituitary fissures for 11 subjects between 8 and 14 years of age.

data("pituitary")
str(pituitary)

## 'data.frame':    11 obs. of  2 variables:
##  $ age : num  8 8 8 10 10 10 12 12 12 14 ...
##  $ size: num  21 23.5 23 24 21 25 21.5 22 19 23.5 ...

Since the size is expected to increase with age, an isotonic fit is suggested.

res_p <- with(pituitary, gpava(age, size))

The CVXR formulation expresses this pretty much in the mathematical form.

x_p <- with(pituitary, {
      n <- length(size)
      x <- Variable(n)
      objective <- Minimize(p_norm(size - x, 2))
      constraint <- list(diff(x) >= 0)
      problem <- Problem(objective, constraint)
      result <- solve(problem)
      result$getValue(x)
})

We define a variable x of size n, the number of observations. The objective to be minimized is the least-squares error (p_norm). The constraint specifies monotonicity and it uses the diff function familiar to R users (see ?diff), except that it now refers to elements of the x. A problem is created and called via solve. The resulting optimal value for x is retrieved using getValue(x).

As the output below shows, the results are very close.

cbind(res_p$x, x_p)

##         [,1]     [,2]
##  [1,] 21.000 20.99998
##  [2,] 22.375 22.37499
##  [3,] 22.375 22.37499
##  [4,] 22.375 22.37499
##  [5,] 22.375 22.37499
##  [6,] 22.375 22.37499
##  [7,] 22.375 22.37499
##  [8,] 22.375 22.37499
##  [9,] 22.375 22.37499
## [10,] 23.500 23.50003
## [11,] 25.000 25.00004

Handling Ties

Package isotone provides additional methods for handling tied data besides the default ties = "primary" method; ties = "secondary" enforces equality within ties, and ties = "tertiary" enforces monotonicity on the means. (The latter may cause individual fits to be non-monotonic.)

res_s <- with(pituitary, gpava(age, size, ties = "secondary"))
res_t <- with(pituitary, gpava(age, size, ties = "tertiary"))

In CVXR, the secondary method just requires an additional constraint to enforce equality within tied values; no other modification is necessary. We do that below by figuring out the tied observation indices using base::split and force those x values to be equal (i.e. diff == 0).

x_s <- with(pituitary, {
    n <- length(size)
    x <- Variable(n)
    objective <- Minimize(p_norm(size - x, 2))
    secondary_constraints <- lapply(base::split(x = seq_len(n),
                                                f = age),
                                    function(i) diff(x[i]) == 0)
    constraint <- c(diff(x) >= 0,
                    secondary_constraints)
    problem <- Problem(objective, constraint)
    solve(problem)$getValue(x)
})

Finally, the tertiary method requires computing the block means to be computed for use in enforcing monotonicity.

x_t <- with(pituitary, {
    n <- length(size)
    x <- Variable(n)
    objective <- Minimize(p_norm(size - x, 2))
    blocks <- base::split(x = seq_len(n),
                          f = pituitary$age)
    block_means <- lapply(blocks, function(i) {
        v <- numeric(n)
        v[i] <- 1.0 / length(i)
        matrix(v, nrow = 1) %*% x
    })
    block_mean_vector <- do.call(vstack, block_means)
    constraint <- list(diff(block_mean_vector) >= 0)
    problem <- Problem(objective, constraint)
    solve(problem)$getValue(x)
})

In the above, we make use of the vstack function to create a single vector of the block means.

We can now compare all three solutions.

Notes

Other arbitrary desired behavior can be specified for tied observations using the above as a guide.

Session Info

sessionInfo()

## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin21.6.0 (64-bit)
## Running under: macOS Ventura 13.0
## 
## Matrix products: default
## BLAS:   /usr/local/Cellar/openblas/0.3.21/lib/libopenblasp-r0.3.21.dylib
## LAPACK: /usr/local/Cellar/r/4.2.1_4/lib/R/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices datasets  utils     methods   base     
## 
## other attached packages:
## [1] kableExtra_1.3.4 isotone_1.1-0    CVXR_1.0-11     
## 
## loaded via a namespace (and not attached):
##  [1] xfun_0.34         bslib_0.4.0       slam_0.1-50       lattice_0.20-45  
##  [5] Rmosek_10.0.25    colorspace_2.0-3  htmltools_0.5.3   viridisLite_0.4.1
##  [9] yaml_2.3.6        gmp_0.6-6         rlang_1.0.6       jquerylib_0.1.4  
## [13] glue_1.6.2        Rmpfr_0.8-9       Rcplex_0.3-5      bit64_4.0.5      
## [17] lifecycle_1.0.3   stringr_1.4.1     munsell_0.5.0     blogdown_1.13    
## [21] rvest_1.0.3       gurobi_9.5-2      codetools_0.2-18  evaluate_0.17    
## [25] knitr_1.40        fastmap_1.1.0     cccp_0.2-9        highr_0.9        
## [29] Rcpp_1.0.9        scales_1.2.1      cachem_1.0.6      webshot_0.5.4    
## [33] jsonlite_1.8.3    systemfonts_1.0.4 bit_4.0.4         digest_0.6.30    
## [37] stringi_1.7.8     bookdown_0.29     grid_4.2.1        Rglpk_0.6-4      
## [41] ECOSolveR_0.5.4   cli_3.4.1         tools_4.2.1       magrittr_2.0.3   
## [45] sass_0.4.2        rcbc_0.1.0.9001   Matrix_1.5-1      xml2_1.3.3       
## [49] nnls_1.4          assertthat_0.2.1  rmarkdown_2.17    svglite_2.1.0    
## [53] httr_1.4.4        rstudioapi_0.14   R6_2.5.1          compiler_4.2.1

Source

R Markdown