Difference in BF computed by numerial intergration vs. ABF

Last updated: 2023-09-24

Checks: 7 0

Knit directory: survival-susie/

This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20230201)

The command set.seed(20230201) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: b2ab762

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version b2ab762. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .DS_Store
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    data/.DS_Store

Unstaged changes:
    Modified:   analysis/run_ser_simple_dat.Rmd
    Modified:   analysis/ser_survival.Rmd
    Modified:   analysis/time_comparison.Rmd
    Modified:   data/dsc3/susie.lbf.rds

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/bf_diff.Rmd) and HTML (docs/bf_diff.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	b2ab762	yunqiyang0215	2023-09-24	wflow_publish("analysis/bf_diff.Rmd")
html	3f6c416	yunqiyang0215	2023-09-15	Build site.
Rmd	5281280	yunqiyang0215	2023-09-15	wflow_publish("analysis/bf_diff.Rmd")
html	74ef6f7	yunqiyang0215	2023-09-15	Build site.
Rmd	27950c6	yunqiyang0215	2023-09-15	wflow_publish("analysis/bf_diff.Rmd")
html	f92d5ae	yunqiyang0215	2023-08-01	Build site.
html	0527948	yunqiyang0215	2023-07-31	Build site.
Rmd	dbbc9ce	yunqiyang0215	2023-07-31	wflow_publish("analysis/bf_diff.Rmd")
html	11f0c28	yunqiyang0215	2023-07-31	Build site.
Rmd	f661067	yunqiyang0215	2023-07-31	wflow_publish("analysis/bf_diff.Rmd")

Description:

Assess difference in Bayes factor computed by numerical integration vs. ABF. The example of using numerical integration to compute BF is available at: https://yunqiyang0215.github.io/survival-susie/numerical_integration.html

We assess three scenarios under relatively large sample size, \(n= 10,000\).

Extremely high censoring rate, censoring rate = 0.99
Low allele frequency, MAF = 0.001
Both.

Result shows that low allele frequency makes the difference between BF and ABF larger.

# Function to simulate survival time under exponential model. That is,
# assuming survival time is exponentially distributed.
# lambda(t) = lambda*exp(b0 + Xb). S(t) = exp(-lambda*t),
# F(t) = 1- S(t) \sim Unif(0,1). Therefore, t = log(1-runif(0,1))/-exp(b0+Xb).
# For censored objects, we simulate the censoring time by rescale the actual survival time.
# @param b: vector of length (p+1) for true effect size, include intercept.
# @param X: variable matrix of size n by p.
# @param censor_lvl: a constant from [0,1], indicating the censoring level in the data.
# @return  dat: a dataframe that contains `y`, `x` and `status`.
# `status`: censoring status: 0 = censored, 1 = event observed. See Surv() in library(survival)
sim_surv_exp <- function(b, X, censor_lvl){
  n = nrow(X)
  p = ncol(X)
  dat = list()

  status <- ifelse(runif(n) > censor_lvl, 1, 0)
  lambda <- exp(cbind(rep(1,n), X) %*% b)
  surT <- log(1 - runif(n)) /(-lambda)
  # rescale censored subject to get observed time
  surT[status == 0] = surT[status == 0] * runif(sum(status == 0))

  y = cbind(surT, status)
  colnames(y) = c("time", "status")
  colnames(X) <- unlist(lapply(1:p, function(i) paste0("x", i)))
  dat[["X"]] = X
  dat[["y"]] = y
  return(dat)
}

compute_multiple_lbf <- function(x, y, o, prior_variance, zscore.spa, ...){
  fit <- coxph(y~ x + offset(o))
  bhat <- summary(fit)$coefficients[1, 1] # bhat = -alphahat
  sd <- summary(fit)$coefficients[1, 3]
  zscore <- bhat/sd
  lbf <- compute_lbf(zscore, sd, prior_variance)
  lbf.corr <- lbf - bhat^2/sd^2/2+ as.numeric(summary(fit)$logtest[1]/2)
  lbf.spa <- compute_lbf(zscore.spa, sd, prior_variance)
  return(list(zscore = zscore, sd = sd, 
              lbf=lbf, lbf.corr = lbf.corr, lbf.spa = lbf.spa))
}

# @param b: the effect size
# @param surT: a Surv() object, containing time and status
# @param x: covariate
get_partial_lik <- function(b, survT, x){
  ll.partial = logLik(coxph(survT ~ x, init = c(b), control=list('iter.max'= 0, timefix=FALSE)))
  max.ll.partial = logLik(coxph(survT ~ x))
  lik.partial.normed = exp(ll.partial - max.ll.partial) # partial likelihood / max(partial likelihood)
  return(as.numeric(lik.partial.normed))
}

integrand <- function(b, survT, x, prior_sd){
  val = get_partial_lik(b, survT, x) * dnorm(b, mean = 0, sd = prior_sd)
  return(val)
}

get_zscore_spa <- function(dat, X){
  # genome matrix
  n = nrow(X)
  Geno.mtx = X
  Phen.mtx = data.frame(ID = paste0("IID-",1:n),
                      event=dat$y[,2],
                      time=dat$y[,1])
  rownames(Geno.mtx) = paste0("IID-",1:nrow(Geno.mtx))
  colnames(Geno.mtx) = paste0("SNP-",1:ncol(Geno.mtx))
  obj.null = SPACox_Null_Model(Surv(time,event) ~ 1, data=Phen.mtx, 
                               pIDs=Phen.mtx$ID, gIDs=rownames(Geno.mtx))
  SPACox.res = SPACox(obj.null, Geno.mtx)
  zscore = SPACox.res[1, 7]
  return(zscore)
}

library(survival)
library(cubature)
library(SPACox)

Loading required package: seqminer

Loading required package: data.table

source("./code/surv_susie_helper.R")

Scenario 1:

seeds = c(1:5)
result = matrix(NA, ncol = 5, nrow = length(seeds))
colnames(result) = c("lbf.wakefeld", "lbf.laplace", "lbf.numerical", "lbf.spa", "error")
for (seed in seeds){
  set.seed(seed)
  # simulate 1 variable
  n = 1e4
  X = cbind(rbinom(n, size = 2, prob = 0.2))
  # the first element of b is for intercept
  b = c(1, 1)
  censor_lvl = 0.99
  dat <- sim_surv_exp(b, X, censor_lvl)
  survT <- Surv(dat$y[,1], dat$y[,2])

  # compute lbf using numerical integration 
  prior_sd = 1
  res = cubintegrate(f = integrand, survT = survT, x = dat$X[,1], prior_sd = 1, lower = -100, upper = 100, method = "hcubature")
  max.ll.partial = logLik(coxph(survT ~ dat$X[,1]))
  prob.h1 = res$integral*exp(max.ll.partial)
  prob.h0 = get_partial_lik(b = 0, survT, dat$X[,1])
  lbf.numerical = log(res$integral) - log(prob.h0)
  
  # compute other lbf
  zscore.spa <- get_zscore_spa(dat, X)
  res.lbf = compute_multiple_lbf(x = dat$X[,1], y = survT, o = rep(0, n), prior_variance = prior_sd^2, zscore.spa = zscore.spa)
  
  
  result[seed, 1] = res.lbf$lbf
  result[seed, 2] = res.lbf$lbf.corr
  result[seed, 3] = lbf.numerical
  result[seed, 4] = res.lbf$lbf.spa
  result[seed, 5] = res$error
}

[1] "Start calculating empirical CGF for martingale residuals..."
[1] "Complete 1000/10000."
[1] "Complete 2000/10000."
[1] "Complete 3000/10000."
[1] "Complete 4000/10000."
[1] "Complete 5000/10000."
[1] "Complete 6000/10000."
[1] "Complete 7000/10000."
[1] "Complete 8000/10000."
[1] "Complete 9000/10000."
[1] "Complete 10000/10000."
[1] "Sample size is 10000."
[1] "Number of variants is 1."
[1] "Start Analyzing..."
[1] "2023-09-24 13:48:02 CDT"
[1] "Analysis Complete."
[1] "2023-09-24 13:48:02 CDT"
[1] "Start calculating empirical CGF for martingale residuals..."
[1] "Complete 1000/10000."
[1] "Complete 2000/10000."
[1] "Complete 3000/10000."
[1] "Complete 4000/10000."
[1] "Complete 5000/10000."
[1] "Complete 6000/10000."
[1] "Complete 7000/10000."
[1] "Complete 8000/10000."
[1] "Complete 9000/10000."
[1] "Complete 10000/10000."
[1] "Sample size is 10000."
[1] "Number of variants is 1."
[1] "Start Analyzing..."
[1] "2023-09-24 13:48:27 CDT"
[1] "Analysis Complete."
[1] "2023-09-24 13:48:27 CDT"
[1] "Start calculating empirical CGF for martingale residuals..."
[1] "Complete 1000/10000."
[1] "Complete 2000/10000."
[1] "Complete 3000/10000."
[1] "Complete 4000/10000."
[1] "Complete 5000/10000."
[1] "Complete 6000/10000."
[1] "Complete 7000/10000."
[1] "Complete 8000/10000."
[1] "Complete 9000/10000."
[1] "Complete 10000/10000."
[1] "Sample size is 10000."
[1] "Number of variants is 1."
[1] "Start Analyzing..."
[1] "2023-09-24 13:48:54 CDT"
[1] "Analysis Complete."
[1] "2023-09-24 13:48:54 CDT"
[1] "Start calculating empirical CGF for martingale residuals..."
[1] "Complete 1000/10000."
[1] "Complete 2000/10000."
[1] "Complete 3000/10000."
[1] "Complete 4000/10000."
[1] "Complete 5000/10000."
[1] "Complete 6000/10000."
[1] "Complete 7000/10000."
[1] "Complete 8000/10000."
[1] "Complete 9000/10000."
[1] "Complete 10000/10000."
[1] "Sample size is 10000."
[1] "Number of variants is 1."
[1] "Start Analyzing..."
[1] "2023-09-24 13:49:16 CDT"
[1] "Analysis Complete."
[1] "2023-09-24 13:49:16 CDT"
[1] "Start calculating empirical CGF for martingale residuals..."
[1] "Complete 1000/10000."
[1] "Complete 2000/10000."
[1] "Complete 3000/10000."
[1] "Complete 4000/10000."
[1] "Complete 5000/10000."
[1] "Complete 6000/10000."
[1] "Complete 7000/10000."
[1] "Complete 8000/10000."
[1] "Complete 9000/10000."
[1] "Complete 10000/10000."
[1] "Sample size is 10000."
[1] "Number of variants is 1."
[1] "Start Analyzing..."
[1] "2023-09-24 13:49:45 CDT"
[1] "Analysis Complete."
[1] "2023-09-24 13:49:45 CDT"

result

     lbf.wakefeld lbf.laplace lbf.numerical   lbf.spa        error
[1,]     15.21520    12.33433      12.34616  5.826518 1.147181e-07
[2,]     18.73883    15.31115      15.32139  7.188096 1.543566e-07
[3,]     28.33746    23.63167      23.63849 13.066644 3.369106e-07
[4,]     20.55204    16.72506      16.73613  8.243768 1.205324e-07
[5,]     27.03120    22.22531      22.23338 11.701627 1.220578e-07

Scenario 2:

seeds = c(1:5)
result = matrix(NA, ncol = 5, nrow = length(seeds))
colnames(result) = c("lbf.wakefeld", "lbf.laplace", "lbf.numerical", "lbf.spa", "error")
for (seed in seeds){
  set.seed(seed)
  # simulate 1 variable
  n = 1e4
  X = cbind(rbinom(n, size = 2, prob = 0.001))
  # the first element of b is for intercept
  b = c(1, 1)
  censor_lvl = 0.6
  dat <- sim_surv_exp(b, X, censor_lvl)
  survT <- Surv(dat$y[,1], dat$y[,2])

  
  # compute lbf using numerical integration 
  prior_sd = 1
  res = cubintegrate(f = integrand, survT = survT, x = dat$X[,1], prior_sd = 1, lower = -100, upper = 100, method = "hcubature")
  max.ll.partial = logLik(coxph(survT ~ dat$X[,1]))
  prob.h1 = res$integral*exp(max.ll.partial)
  prob.h0 = get_partial_lik(b = 0, survT, dat$X[,1])
  lbf.numerical = log(res$integral) - log(prob.h0)
  
  # compute other lbf
  zscore.spa <- get_zscore_spa(dat, X)
  res.lbf = compute_multiple_lbf(x = dat$X[,1], y = survT, o = rep(0, n), prior_variance = prior_sd^2, zscore.spa = zscore.spa)
  
  
  result[seed, 1] = res.lbf$lbf
  result[seed, 2] = res.lbf$lbf.corr
  result[seed, 3] = lbf.numerical
  result[seed, 4] = res.lbf$lbf.spa
  result[seed, 5] = res$error
}

[1] "Start calculating empirical CGF for martingale residuals..."
[1] "Complete 1000/10000."
[1] "Complete 2000/10000."
[1] "Complete 3000/10000."
[1] "Complete 4000/10000."
[1] "Complete 5000/10000."
[1] "Complete 6000/10000."
[1] "Complete 7000/10000."
[1] "Complete 8000/10000."
[1] "Complete 9000/10000."
[1] "Complete 10000/10000."
[1] "Sample size is 10000."
[1] "Number of variants is 1."
[1] "Start Analyzing..."
[1] "2023-09-24 13:50:18 CDT"
[1] "Analysis Complete."
[1] "2023-09-24 13:50:18 CDT"
[1] "Start calculating empirical CGF for martingale residuals..."
[1] "Complete 1000/10000."
[1] "Complete 2000/10000."
[1] "Complete 3000/10000."
[1] "Complete 4000/10000."
[1] "Complete 5000/10000."
[1] "Complete 6000/10000."
[1] "Complete 7000/10000."
[1] "Complete 8000/10000."
[1] "Complete 9000/10000."
[1] "Complete 10000/10000."
[1] "Sample size is 10000."
[1] "Number of variants is 1."
[1] "Start Analyzing..."
[1] "2023-09-24 13:50:52 CDT"
[1] "Analysis Complete."
[1] "2023-09-24 13:50:52 CDT"
[1] "Start calculating empirical CGF for martingale residuals..."
[1] "Complete 1000/10000."
[1] "Complete 2000/10000."
[1] "Complete 3000/10000."
[1] "Complete 4000/10000."
[1] "Complete 5000/10000."
[1] "Complete 6000/10000."
[1] "Complete 7000/10000."
[1] "Complete 8000/10000."
[1] "Complete 9000/10000."
[1] "Complete 10000/10000."
[1] "Sample size is 10000."
[1] "Number of variants is 1."
[1] "Start Analyzing..."
[1] "2023-09-24 13:51:24 CDT"
[1] "Analysis Complete."
[1] "2023-09-24 13:51:24 CDT"
[1] "Start calculating empirical CGF for martingale residuals..."
[1] "Complete 1000/10000."
[1] "Complete 2000/10000."
[1] "Complete 3000/10000."
[1] "Complete 4000/10000."
[1] "Complete 5000/10000."
[1] "Complete 6000/10000."
[1] "Complete 7000/10000."
[1] "Complete 8000/10000."
[1] "Complete 9000/10000."
[1] "Complete 10000/10000."
[1] "Sample size is 10000."
[1] "Number of variants is 1."
[1] "Start Analyzing..."
[1] "2023-09-24 13:51:55 CDT"
[1] "Analysis Complete."
[1] "2023-09-24 13:51:55 CDT"
[1] "Start calculating empirical CGF for martingale residuals..."
[1] "Complete 1000/10000."
[1] "Complete 2000/10000."
[1] "Complete 3000/10000."
[1] "Complete 4000/10000."
[1] "Complete 5000/10000."
[1] "Complete 6000/10000."
[1] "Complete 7000/10000."
[1] "Complete 8000/10000."
[1] "Complete 9000/10000."
[1] "Complete 10000/10000."
[1] "Sample size is 10000."
[1] "Number of variants is 1."
[1] "Start Analyzing..."
[1] "2023-09-24 13:52:31 CDT"
[1] "Analysis Complete."
[1] "2023-09-24 13:52:31 CDT"

result

     lbf.wakefeld lbf.laplace lbf.numerical    lbf.spa        error
[1,]     1.107410   0.4316504     0.5382978 -0.2156671 1.524437e-06
[2,]     5.556143   3.3666361     3.4197635  1.3888268 6.036238e-08
[3,]     8.496663   4.6024045     4.6882092  1.8509981 6.913447e-08
[4,]     3.066519   1.4709938     1.5918999 -0.1197060 1.052610e-06
[5,]     1.963928   1.1080616     1.1633604  0.2814170 3.301337e-09

Scenario 3:

seeds = c(1:3)
result = matrix(NA, ncol = 5, nrow = length(seeds))
colnames(result) = c("lbf.wakefeld", "lbf.laplace", "lbf.numerical", "lbf.spa", "error")
for (seed in seeds){
  set.seed(seed)
  # simulate 1 variable
  n = 1e4
  X = cbind(rbinom(n, size = 2, prob = 0.001))
  # the first element of b is for intercept
  b = c(1, 1)
  censor_lvl = 0.99
  dat <- sim_surv_exp(b, X, censor_lvl)
  survT <- Surv(dat$y[,1], dat$y[,2])

  
  # compute lbf using numerical integration 
  prior_sd = 1
  res = cubintegrate(f = integrand, survT = survT, x = dat$X[,1], prior_sd = 1, lower = -100, upper = 100, method = "hcubature")
  max.ll.partial = logLik(coxph(survT ~ dat$X[,1]))
  prob.h1 = res$integral*exp(max.ll.partial)
  prob.h0 = get_partial_lik(b = 0, survT, dat$X[,1])
  lbf.numerical = log(res$integral) - log(prob.h0)
  
  # compute other lbf
  zscore.spa <- get_zscore_spa(dat, X)
  res.lbf = compute_multiple_lbf(x = dat$X[,1], y = survT, o = rep(0, n), prior_variance = prior_sd^2, zscore.spa = zscore.spa)
  
  
  result[seed, 1] = res.lbf$lbf
  result[seed, 2] = res.lbf$lbf.corr
  result[seed, 3] = lbf.numerical
  result[seed, 4] = res.lbf$lbf.spa
  result[seed, 5] = res$error
}

[1] "Start calculating empirical CGF for martingale residuals..."
[1] "Complete 1000/10000."
[1] "Complete 2000/10000."
[1] "Complete 3000/10000."
[1] "Complete 4000/10000."
[1] "Complete 5000/10000."
[1] "Complete 6000/10000."
[1] "Complete 7000/10000."
[1] "Complete 8000/10000."
[1] "Complete 9000/10000."
[1] "Complete 10000/10000."
[1] "Sample size is 10000."
[1] "Number of variants is 1."
[1] "Start Analyzing..."
[1] "2023-09-24 13:52:59 CDT"
[1] "Analysis Complete."
[1] "2023-09-24 13:52:59 CDT"
[1] "Start calculating empirical CGF for martingale residuals..."
[1] "Complete 1000/10000."
[1] "Complete 2000/10000."
[1] "Complete 3000/10000."
[1] "Complete 4000/10000."
[1] "Complete 5000/10000."
[1] "Complete 6000/10000."
[1] "Complete 7000/10000."
[1] "Complete 8000/10000."
[1] "Complete 9000/10000."
[1] "Complete 10000/10000."
[1] "Sample size is 10000."
[1] "Number of variants is 1."
[1] "Start Analyzing..."
[1] "2023-09-24 13:53:26 CDT"
[1] "Analysis Complete."
[1] "2023-09-24 13:53:26 CDT"
[1] "Start calculating empirical CGF for martingale residuals..."
[1] "Complete 1000/10000."
[1] "Complete 2000/10000."
[1] "Complete 3000/10000."
[1] "Complete 4000/10000."
[1] "Complete 5000/10000."
[1] "Complete 6000/10000."
[1] "Complete 7000/10000."
[1] "Complete 8000/10000."
[1] "Complete 9000/10000."
[1] "Complete 10000/10000."
[1] "Sample size is 10000."
[1] "Number of variants is 1."
[1] "Start Analyzing..."
[1] "2023-09-24 13:53:54 CDT"
[1] "Analysis Complete."
[1] "2023-09-24 13:53:54 CDT"

result

      lbf.wakefeld lbf.laplace lbf.numerical       lbf.spa        error
[1,] -1.713396e-07  0.02051070   -0.01238313 -1.706615e-07 3.969705e-08
[2,] -1.491693e-07  0.04862419   -0.02675245 -1.477403e-07 4.480898e-08
[3,] -1.712154e-07  0.05576799   -0.03009710 -1.690265e-07 4.704538e-08

sessionInfo()

R version 4.1.1 (2021-08-10)
Platform: x86_64-apple-darwin20.6.0 (64-bit)
Running under: macOS Monterey 12.0.1

Matrix products: default
BLAS:   /usr/local/Cellar/openblas/0.3.18/lib/libopenblasp-r0.3.18.dylib
LAPACK: /usr/local/Cellar/r/4.1.1_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] SPACox_0.1.2      data.table_1.14.9 seqminer_9.1      cubature_2.1.0   
[5] survival_3.2-11   workflowr_1.6.2  

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.8.3     pillar_1.6.4     compiler_4.1.1   bslib_0.4.1     
 [5] later_1.3.0      jquerylib_0.1.4  git2r_0.28.0     tools_4.1.1     
 [9] digest_0.6.28    lattice_0.20-44  jsonlite_1.7.2   evaluate_0.14   
[13] lifecycle_1.0.3  tibble_3.1.5     pkgconfig_2.0.3  rlang_1.1.1     
[17] Matrix_1.5-3     cli_3.6.1        rstudioapi_0.13  yaml_2.2.1      
[21] xfun_0.27        fastmap_1.1.0    stringr_1.4.0    knitr_1.36      
[25] fs_1.5.0         vctrs_0.6.3      sass_0.4.4       grid_4.1.1      
[29] rprojroot_2.0.2  glue_1.4.2       R6_2.5.1         fansi_0.5.0     
[33] rmarkdown_2.11   magrittr_2.0.1   whisker_0.4      splines_4.1.1   
[37] promises_1.2.0.1 ellipsis_0.3.2   htmltools_0.5.5  httpuv_1.6.3    
[41] utf8_1.2.2       stringi_1.7.5    cachem_1.0.6     crayon_1.4.1

Difference in BF computed by numerial intergration vs. ABF

Yunqi Yang

7/31/2023

Description:

Scenario 1:

Scenario 2:

Scenario 3: