Prototyping Distributed Computations

Introduction

The distcomp package allows one to perform computations on data that is row-partitioned data among several sites. Implemented methods include stratified Cox regression and singular value decomposition (SVD). The distributed computations rely on a REST API for R, such as that provided by Opencpu for example.

In order to prototype new methods for distcomp, it is convenient to write code that does all computations locally, at least during the development phase. This allows for easy debugging and code changes. The step to using a REST call is merely one step off. For example, in the context of distributed multiple imputation, the distributed SVD is a building block. So one has to be able to be able to run the svd operation not only in a modular fashion, but also quickly to test code.

Previous versions of distcomp assumed all computations were done using opencpu remote calls. This meant that the results were time-consuming to check; for reasonably sized problems one would have to wait for hours.

Version 1.1 of distcomp addresses this by allowing local computations, in essence faking the opencpu calls. All calculations are done locally at normal speeds (rather than a https POST request) while ensuring that the sites are only privy to their own data.

An Example

available <- availableComputations()
( computationNames <- names(available) )

## [1] "QueryCount"         "StratifiedCoxModel" "RankKSVD"

We will perform an SVD on a 1500 by 20 matrix where the matrix is distributed across three sites each having a 500 by 20 matrix. We will extract the first 10 components using this distributed algorithm.

We first define the computation.

svdDef <- data.frame(compType = "RankKSVD",
                     rank = 20L,
                     ncol = 20L,
                     id = "SVD",
                     stringsAsFactors = FALSE)

We generate some synthetic data.

set.seed(12345)
## Three sites
nSites <- 3
siteData <- lapply(seq.int(nSites), function(i) matrix(rnorm(10000), ncol=20))

We now create local objects to represent the sites, a list of three items, each with a name and its specific data.

sites <- lapply(seq.int(nSites),
                function(x) list(name = paste0("site", x),
                                 worker = makeWorker(defn = svdDef, data = siteData[[x]])
                                 ))

We are now ready to create a master object and add sites to it.

master <- makeMaster(svdDef)
for (site in sites) {
  master$addSite(name = site$name, worker = site$worker)
}

When any of the sites specify a worker argument in the call to addSite, it is an indication that the computation is to be performed locally. Otherwise, a URL will have to be specified using the url argument to addSite. (Only one of worker or url can be specified.)

Finally, we are ready to run the decomposition, specifying a maximum number of iterations, so that things don’t go haywire. This takes less than 5 seconds on my laptop.

system.time(result <- master$run(max.iter = 10000))

##    user  system elapsed 
##   9.681  11.376  10.395

The result is a named list of two things, the \(d\) diagonal values and the \(v\) matrix.

We can compare the results obtained here with the actual.

full_data <- do.call(rbind, siteData)
full_svd <- svd(full_data)

Let’s print it side-by-side.

d_table <- data.frame(truth = full_svd$d, distcomp = result$d)
knitr::kable(d_table)

truth	distcomp
41.82263	41.82263
41.65242	41.65242
41.32027	41.32027
40.88629	40.88629
40.41849	40.41849
40.27631	40.27631
40.12698	40.12698
39.44775	39.44775
39.24815	39.24815
39.07074	39.07074
38.37076	38.37076
38.21350	38.21350
37.61047	37.61047
37.14286	37.14286
36.98763	36.98763
36.24778	36.24778
36.04962	36.04962
35.66050	35.66050
35.21343	35.21343
34.56329	34.56329

We can also compare the \(v\) matrix, taking into account the signs may be different.

norm(abs(result$v) - abs(full_svd$v), "F")

## [1] 1.124008e-05

Stratified Cox Model

The approach for the stratified Cox model is similar. We start with a definition.

coxDef <- data.frame(compType = "StratifiedCoxModel",
                     formula = "Surv(time, censor) ~ age + becktota + ndrugfp1 + ndrugfp2 + ivhx3 + race + treat",
                     projectName = "STCoxTest",
                     projectDesc = "STCox Project Desc",
                     stringsAsFactors = FALSE)

We now create two sites with appropriate data

## Two sites
siteDataFiles <- file.path(system.file("ex", package="distcomp"), c("uis-site1.csv", "uis-site2.csv"))
siteData <- lapply(siteDataFiles, read.csv)

We now create local objects to represent the two sites, each with a name and its specific data.

sites <- lapply(seq_along(siteData),
                function(x) list(name = paste0("site", x),
                                 worker = makeWorker(defn = coxDef, data = siteData[[x]])
                                 ))

We are now ready to create a master object and add sites to it.

master <- makeMaster(coxDef)
for (site in sites) {
  master$addSite(name = site$name, worker = site$worker)
}

Now we can run the model.

result <- master$run()
print(master$summary(), digits = 4)

##        coef exp(coef) se(coef)      z         p
## 1 -0.028076    0.9723 0.008131 -3.453 5.542e-04
## 2  0.009146    1.0092 0.004991  1.832 6.691e-02
## 3 -0.521973    0.5933 0.124424 -4.195 2.727e-05
## 4 -0.194178    0.8235 0.048252 -4.024 5.717e-05
## 5  0.263634    1.3017 0.108243  2.436 1.487e-02
## 6 -0.240021    0.7866 0.115632 -2.076 3.792e-02
## 7 -0.212616    0.8085 0.093747 -2.268 2.333e-02

The above result will be the same as what we obtain from a Cox regression fit.

coxFit <- survival::coxph(formula=Surv(time, censor) ~ age + becktota + ndrugfp1 + ndrugfp2 +
                              ivhx3 + race + treat + strata(site),
                          data = rbind(siteData[[1]], siteData[[2]]))
summary(coxFit)$coefficients

##                   coef exp(coef)    se(coef)         z     Pr(>|z|)
## age       -0.028075893 0.9723146 0.008130685 -3.453078 5.542280e-04
## becktota   0.009145528 1.0091875 0.004991421  1.832250 6.691425e-02
## ndrugfp1  -0.521973047 0.5933487 0.124423881 -4.195119 2.727278e-05
## ndrugfp2  -0.194177573 0.8235117 0.048252289 -4.024215 5.716573e-05
## ivhx3TRUE  0.263634281 1.3016521 0.108243388  2.435569 1.486837e-02
## race      -0.240020862 0.7866115 0.115632433 -2.075723 3.791961e-02
## treat     -0.212616368 0.8084662 0.093747124 -2.267978 2.333058e-02

Session Info

sessionInfo()

## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin21.5.0 (64-bit)
## Running under: macOS Monterey 12.5.1
## 
## Matrix products: default
## BLAS:   /usr/local/Cellar/openblas/0.3.20/lib/libopenblasp-r0.3.20.dylib
## LAPACK: /usr/local/Cellar/r/4.2.1/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 utils     datasets  methods   base     
## 
## other attached packages:
## [1] distcomp_1.3-3 survival_3.3-1
## 
## loaded via a namespace (and not attached):
##  [1] gmp_0.6-5         Rcpp_1.0.9        highr_0.9         pillar_1.8.0     
##  [5] bslib_0.4.0       compiler_4.2.1    later_1.3.0       jquerylib_0.1.4  
##  [9] tools_4.2.1       digest_0.6.29     tibble_3.1.8      jsonlite_1.8.0   
## [13] evaluate_0.15     lifecycle_1.0.1   lattice_0.20-45   pkgconfig_2.0.3  
## [17] rlang_1.0.4       Matrix_1.4-1      DBI_1.1.3         shiny_1.7.2      
## [21] cli_3.3.0         yaml_2.3.5        xfun_0.31         fastmap_1.1.0    
## [25] httr_1.4.3        stringr_1.4.0     dplyr_1.0.9       knitr_1.39       
## [29] generics_0.1.3    sass_0.4.2        vctrs_0.4.1       tidyselect_1.1.2 
## [33] grid_4.2.1        glue_1.6.2        R6_2.5.1          fansi_1.0.3      
## [37] rmarkdown_2.14    homomorpheR_0.2-2 purrr_0.3.4       magrittr_2.0.3   
## [41] promises_1.2.0.1  htmltools_0.5.3   ellipsis_0.3.2    splines_4.2.1    
## [45] assertthat_0.2.1  mime_0.12         xtable_1.8-4      httpuv_1.6.5     
## [49] utf8_1.2.2        stringi_1.7.8     sodium_1.2.1      cachem_1.0.6