## ----echo=FALSE--------------------------------------------------------------- ### get knitr just the way we like it knitr::opts_chunk$set( message = FALSE, warning = FALSE, error = FALSE, tidy = FALSE, cache = FALSE ) ## ----------------------------------------------------------------------------- if (!require("survival") || !requireNamespace("digest", quietly = TRUE)) { stop("this vignette requires both the survival & digest package") } library(homomorpheR) ## ----------------------------------------------------------------------------- sampleSize <- c(n1 = 1000, n2 = 500, n3 = 1500) set.seed(12345) beta.1 <- -.015; beta.2 <- .2; beta.3 <- .001; lambdaT <- c(5, 4, 3) lambdaC <- 2 coxData <- lapply(seq_along(sampleSize), function(i) { sex <- sample(c(0, 1), size = sampleSize[i], replace = TRUE) age <- sample(40:70, size = sampleSize[i], replace = TRUE) bm <- rnorm(sampleSize[i]) trueTime <- rweibull(sampleSize[i], shape = 1, scale = lambdaT[i] * exp(beta.1 * age + beta.2 * sex + beta.3 * bm )) censoringTime <- rweibull(sampleSize[i], shape = 1, scale = lambdaC) time <- pmin(trueTime, censoringTime) event <- (time == trueTime) data.frame(stratum = i, sex = sex, age = age, bm = bm, time = time, event = event) }) ## ----------------------------------------------------------------------------- str(coxData[[1]]) ## ----------------------------------------------------------------------------- str(coxData[[2]]) ## ----------------------------------------------------------------------------- str(coxData[[3]]) ## ----------------------------------------------------------------------------- aggModel <- coxph(formula = Surv(time, event) ~ sex + age + bm + strata(stratum), data = do.call(rbind, coxData)) aggModel ## ----------------------------------------------------------------------------- aggModel$loglik ## ----------------------------------------------------------------------------- Site <- R6::R6Class( "Site", private = list( ## name of the site name = NA, ## local data data = NA, ## Control variable for cox regression cph.control = NA, beta_cache = list(), local_nll = function(beta) { ## Check if value is cached beta_hash <- paste0("b", digest::digest(beta, algo = "xxhash64")) result <- private$beta_cache[[beta_hash]] if (is.null(result)) { ## We're worker, so compute local negative log likelihood nllValue <- tryCatch({ m <- coxph(formula = Surv(time, event) ~ sex + age + bm, data = private$data, init = beta, control = private$cph.control) -(m$loglik[1]) }, error = function(e) NA) if (!is.na(nllValue)) { pubkey <- self$pubkey ## Generate random offset for int and frac parts offset <- list(int = random.bigz(nBits = 256), frac = random.bigz(nBits = 256)) ## 2. Add to neg log likelihood result.int <- floor(nllValue) result.frac <- nllValue - result.int ## Approximate fractional part by a rational result.fracnum <- gmp::as.bigz(gmp::numerator(gmp::as.bigq(result.frac) * self$den)) result <- list( int1 = pubkey$encrypt(result.int - offset$int), frac1 = pubkey$encrypt(result.fracnum - offset$frac), int2 = pubkey$encrypt(result.int + offset$int), frac2 = pubkey$encrypt(result.fracnum + offset$frac) ) private$beta_cache[[beta_hash]] <- result } else { result <- list(int1 = NA, frac1 = NA, int2 = NA, frac2 = NA) } } result } ), public = list( count = NA, ## Common denominator for approximate real arithmetic den = NA, ## The master's public key; everyone has this pubkey = NA, initialize = function(name, data) { private$name <- name private$data <- data private$cph.control <- replace(coxph.control(), "iter.max", 0) }, setPublicKey = function(pubkey) { self$pubkey <- pubkey }, setDenominator = function(den) { self$den = den }, ## neg log lik, nll = function(beta, party) { result <- private$local_nll(beta) if (party == 1) { list(int = result$int1, frac = result$frac1) } else { list(int = result$int2, frac = result$frac2) } } ) ) ## ----------------------------------------------------------------------------- NCParty <- R6::R6Class( "NCParty", private = list( ## name of the site name = NA, ## NC party number number = NA, ## The master master = NA, ## The sites sites = list() ), public = list( ## The master's public key; everyone has this pubkey = NA, ## The denoinator for rational arithmetic den = NA, initialize = function(name, number) { private$name <- name private$number <- number }, setPublicKey = function(pubkey) { self$pubkey <- pubkey ## Propagate to sites for (site in sites) { site$setPublicKey(pubkey) } }, setDenominator = function(den) { self$den <- den ## Propagate to sites for (site in sites) { site$setDenominator(den) } }, addSite = function(site) { private$sites <- c(private$sites, list(site)) }, ## neg log lik nll = function(beta) { pubkey <- self$pubkey results <- lapply(sites, function(x) x$nll(beta, private$number)) ## Accumulate the integer and fractional parts n <- length(results) sumInt <- results[[1L]]$int sumFrac <- results[[1L]]$frac for (i in 2:n) { sumInt <- pubkey$add(sumInt, results[[i]]$int) sumFrac <- pubkey$add(sumFrac, results[[i]]$frac) } list(int = sumInt, frac = sumFrac) } ) ) ## ----------------------------------------------------------------------------- Master <- R6::R6Class( "Master", private = list( ## name of the site name = NA, ## Private and public keys keys = NA, ## Non cooperating party 1 nc_party_1 = NA, ## Non cooperating party 2 nc_party_2 = NA ), public = list( ## Denominator for rational arithmetic den = NA, initialize = function(name) { private$name <- name private$keys <- PaillierKeyPair$new(1024) ## Generate new public and private key. self$den <- gmp::as.bigq(2)^256 #Our denominator for rational approximations }, setNCParty1 = function(site) { private$nc_party_1 <- site private$nc_party_1$setPublicKey(private$keys$pubkey) private$nc_party_1$setDenominator(self$den) }, setNCParty2 = function(site) { private$nc_party_2 <- site private$nc_party_2$setPublicKey(private$keys$pubkey) private$nc_party_2$setDenominator(self$den) }, ## neg log lik nLL = function(beta) { pubkey <- private$keys$pubkey privkey <- private$keys$getPrivateKey() result1 <- private$nc_party_1$nll(beta) result2 <- private$nc_party_2$nll(beta) ## Accumulate the integer and fractional parts sumInt <- pubkey$add(result1$int, result2$int) sumFrac <- pubkey$add(result1$frac, result2$frac) intResult <- as.double(privkey$decrypt(sumInt)) fracResult <- as.double(gmp::as.bigq(privkey$decrypt(sumFrac)) / self$den) ## Since we 2L, we divide by 2. (intResult + fracResult) / 2.0 } ) ) ## ----------------------------------------------------------------------------- site1 <- Site$new(name = "Site 1", data = coxData[[1]]) site2 <- Site$new(name = "Site 2", data = coxData[[2]]) site3 <- Site$new(name = "Site 3", data = coxData[[3]]) sites <- list(site1 = site1, site2 = site2, site3 = site3) ## ----------------------------------------------------------------------------- ncp1 <- NCParty$new("NCP1", 1) ncp2 <- NCParty$new("NCP1", 2) ## ----------------------------------------------------------------------------- for (s in sites) { ncp1$addSite(s) ncp2$addSite(s) } ## ----------------------------------------------------------------------------- master <- Master$new("Master") ## ----------------------------------------------------------------------------- master$setNCParty1(ncp1) master$setNCParty2(ncp2) ## ----------------------------------------------------------------------------- library(stats4) nll <- function(age, sex, bm) master$nLL(c(age, sex, bm)) fit <- mle(nll, start = list(age = 0, sex = 0, bm = 0)) ## ----------------------------------------------------------------------------- summary(fit) ## ----------------------------------------------------------------------------- summary(aggModel) ## ----------------------------------------------------------------------------- cat(sprintf("logLik(MLE fit): %f, logLik(Agg. fit): %f.\n", logLik(fit), aggModel$loglik[2]))