OSCR

Double shrinkage transfer causal learning: An application to alzheimer's disease.

Code ↔ Paper

5 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 5 matches · 2 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § 2. Materials and methods › 2.3. Proposed method › 2.3.3 Estimation of causal effect. ↔ simulate_dstcl_main.R, lines 156–206 · score 0.70 · l1 regularized, lambda.min, lambda.1se, error, alpha, glmnet
  2. [2] § 2. Materials and methods › 2.3. Proposed method › 2.3.3 Estimation of causal effect. ↔ real_data_analysis_template.R, lines 334–393 · score 0.64 · lambda.min, lambda.1se, deviance, fold, binomial, alpha
  3. [3] § 3. Simulation › 3.3. Simulation results ↔ glmtrans.R, the whole file · a weak match · score 0.54 · lambda.min, lambda.1se, penalty, bias, fit
  4. [4] § 3. Simulation › 3.3. Simulation results ↔ simulate_dstcl_main.R, lines 156–206 · score 0.53 · lambda.min, lambda.1se, penalty, fit, DSTCL, Simulation
  5. [5] § 3. Simulation › 3.2. Simulation scenarios ↔ glmtrans.R, the whole file · a weak match · score 0.50 · lambda.min, lambda.1se, penalty

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 685 lines · 24 KB · no license · 2 matches

  1. required_pkgs <- c(
  2. "MASS", "glmnet", "Matrix", "parallel",
  3. "DoubleML", "mlr3", "mlr3learners", "tmle",
  4. "ranger", "e1071", "SuperLearner"
  5. )
  7. missing_pkgs <- required_pkgs[!sapply(required_pkgs, requireNamespace, quietly = TRUE)]
  8. if (length(missing_pkgs) > 0) {
  9. stop("Missing packages: ", paste(missing_pkgs, collapse = ", "),
  10. "\nPlease install them before running this simulation script.")
  11. }
  13. suppressPackageStartupMessages({
  14. library(MASS)
  15. library(glmnet)
  16. library(Matrix)
  17. library(parallel)
  18. library(DoubleML)
  19. library(mlr3)
  20. library(mlr3learners)
  21. library(tmle)
  22. library(ranger)
  23. library(e1071)
  24. library(SuperLearner)
  25. })
  27. if (requireNamespace("lgr", quietly = TRUE)) {
  28. lgr::get_logger("mlr3")$set_threshold("warn")
  29. }
  31. get_script_dir <- function() {
  32. cmd.args <- commandArgs(FALSE)
  33. file.arg <- cmd.args[grep("^--file=", cmd.args)]
  34. if (length(file.arg) > 0) dirname(normalizePath(sub("^--file=", "", file.arg[1]))) else getwd()
  35. }
  37. CODE_DIR <- get_script_dir()
  38. PROJECT_DIR <- normalizePath(file.path(CODE_DIR, ".."), mustWork = FALSE)
  39. output_dir <- file.path(PROJECT_DIR, "results")
  40. dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
  42. source(file.path(CODE_DIR, "gen_latest.R"))
  43. source(file.path(CODE_DIR, "glmtrans.R"))
  45. set.seed(321)
  47. #E(z=1)=0.5, link = 'logit'
  48. n=50 # 50,200,1000,2000
  49. m=1000
  50. p=100
  51. nsim=200 #模拟次数
  52. B=500#重抽样次数
  53. BOOT_CORES <- max(1, min(15, parallel::detectCores() - 1))
  54. BOOT_CLUSTER <- NULL
  56. ######################################
  57. ## setting
  58. alpha1 <- c(0.2, rep(0.3,5),rep(0.6,45),rep(0,50)) #source. ps
  59. alpha2 <- c(-0.2, rep(-0.3,5),rep(0.6,45),rep(0,50)) #target. ps
  60. beta1.t <- c(1, rep(0.3,3),rep(-0.3,2),rep(0.8,25),rep(-0.8,20),rep(0,50))
  61. beta1.c <- c(-1, rep(0.3,3),rep(-0.3,2),rep(0.8,25),rep(-0.8,20),rep(0,50))
  62. beta2.t <- c(-1, rep(-0.3,3),rep(0.3,2),rep(0.8,25),rep(-0.8,20),rep(0,50))
  63. beta2.c <- c(1, rep(-0.3,3),rep(0.3,2),rep(0.8,25),rep(-0.8,20),rep(0,50))
  65. sigmoid <- function(x) 1 / (1 + exp(-x))
  67. ######################################
  68. ## 最小诊断记录对象
  69. ######################################
  70. REG_DIAG <- data.frame()
  71. RUN_STATUS <- data.frame()
  72. ESTIMATOR_TIME <- rep(NA, nsim)
  73. BOOTSTRAP_TIME <- rep(NA, nsim)
  74. TOTAL_TIME <- rep(NA, nsim)
  75. CURRENT_SIM <- NA
  76. CURRENT_STAGE <- "main"
  77. SAVE_REG_DIAG <- TRUE
  78. METHOD_NAMES <- c(
  79. "naive_diff",
  80. "oracle_dr",
  81. "target_lasso_dr",
  82. "doubleml_glmnet_ranger",
  83. "doubleml_svm_ranger",
  84. "tmle_glmnet",
  85. "glmtrans2old",
  86. "dstcl_transfer"
  87. )
  88. MAIN_METHOD_TIME_SUM <- matrix(0, nsim, length(METHOD_NAMES))
  89. MAIN_METHOD_TIME_COUNT <- matrix(0, nsim, length(METHOD_NAMES))
  90. BOOT_METHOD_TIME_SUM <- matrix(0, nsim, length(METHOD_NAMES))
  91. BOOT_METHOD_TIME_COUNT <- matrix(0, nsim, length(METHOD_NAMES))
  92. colnames(MAIN_METHOD_TIME_SUM) <- METHOD_NAMES
  93. colnames(MAIN_METHOD_TIME_COUNT) <- METHOD_NAMES
  94. colnames(BOOT_METHOD_TIME_SUM) <- METHOD_NAMES
  95. colnames(BOOT_METHOD_TIME_COUNT) <- METHOD_NAMES
  97. record_method_time <- function(method_index, elapsed_sec) {
  98. if (is.na(CURRENT_SIM)) return(invisible(NULL))
  99. if (is.character(method_index)) method_index <- match(method_index, METHOD_NAMES)
  100. if (is.na(method_index)) return(invisible(NULL))
  101. if (CURRENT_STAGE == "bootstrap") {
  102. BOOT_METHOD_TIME_SUM[CURRENT_SIM, method_index] <<-
  103. BOOT_METHOD_TIME_SUM[CURRENT_SIM, method_index] + elapsed_sec
  104. BOOT_METHOD_TIME_COUNT[CURRENT_SIM, method_index] <<-
  105. BOOT_METHOD_TIME_COUNT[CURRENT_SIM, method_index] + 1
  106. } else {
  107. MAIN_METHOD_TIME_SUM[CURRENT_SIM, method_index] <<-
  108. MAIN_METHOD_TIME_SUM[CURRENT_SIM, method_index] + elapsed_sec
  109. MAIN_METHOD_TIME_COUNT[CURRENT_SIM, method_index] <<-
  110. MAIN_METHOD_TIME_COUNT[CURRENT_SIM, method_index] + 1
  111. }
  112. invisible(NULL)
  113. }
  115. add_reg_diag <- function(fit, method, component) {
  116. d <- fit$diag
  117. d$sim <- CURRENT_SIM
  118. d$method <- method
  119. d$component <- component
  120. diag_cols <- c("sim", "method", "component",
  121. "nfolds_used", "lambda_rule",
  122. "source_lambda", "correction_lambda",
  123. "nonzero_source", "nonzero_delta", "nonzero_final",
  124. "prop_nonzero_source", "prop_nonzero_delta", "prop_nonzero_final",
  125. "l1_source", "l1_delta", "l1_final",
  126. "n_try_correction", "converged")
  127. missing_cols <- setdiff(diag_cols, names(d))
  128. for (col in missing_cols) {
  129. d[[col]] <- NA
  130. }
  131. d <- d[, diag_cols]
  132. return(d)
  133. }
  135. ######################################
  136. ## 估计函数（baseline + DR + transfer）
  138. or <- function(x,z,beta.t,beta.c){
  139. mean((as.matrix(x)%*%beta.t)-(as.matrix(x)%*%beta.c))
  140. }
  142. ipw <- function(y,z,x,alpha){
  143. p <- sigmoid(as.matrix(x)%*%alpha)
  144. p <- pmin(pmax(p,1e-4),1-1e-4)
  145. sum(y*z/p)/sum(z/p) - sum(y*(1-z)/(1-p))/sum((1-z)/(1-p))
  146. }
  148. dr <- function(y,z,x,alpha,beta.t,beta.c){
  149. ps <- sigmoid(as.matrix(x)%*%alpha)
  150. ps <- pmin(pmax(ps,1e-4),1-1e-4)
  151. mu1 <- as.matrix(x)%*%beta.t
  152. mu0 <- as.matrix(x)%*%beta.c
  153. mean((mu1-mu0)+(z*(y-mu1))/ps-((1-z)*(y-mu0))/(1-ps))
  154. }
  156. glmtrans2old <- function(target=target, source = source, family = "gaussian", lambda = lambda,
  157. type.measure=type.measure, alpha = 1, nfolds = nfolds,
  158. intercept=intercept, penalty.factor=penalty.factor) {
  159. ## Purpose:
  160. ## Reproduce the original l1-TCL comparison method. The source model is
  161. ## fitted by ordinary GLM, followed by an L1-regularized target correction.
  162. beta1 <- glm(y ~ ., family = family, data = source)$coefficients
  163. offset <- as.numeric(as.matrix(target[, -1]) %*% beta1[-1]) + beta1[1]
  165. n.try <- 0
  166. while (TRUE) {
  167. cv.fit.correct <- try(
  168. cv.glmnet(
  169. x = as.matrix(target[, -1]),
  170. y = as.matrix(target[, 1]),
  171. offset = offset,
  172. alpha = 1,
  173. nfolds = nfolds,
  174. intercept = intercept,
  175. type.measure = type.measure,
  176. family = family,
  177. grouped = FALSE,
  178. penalty.factor = penalty.factor
  179. ),
  180. silent = TRUE
  181. )
  182. if (!inherits(cv.fit.correct, "try-error")) break
  183. n.try <- n.try + 1
  184. if (n.try > 10) stop("Errors occur during the debiasing step!!!")
  185. }
  187. if (lambda == "lambda.min") {
  188. lambda.fit <- cv.fit.correct$lambda.min
  189. coef.correct <- coef(cv.fit.correct, s = cv.fit.correct$lambda.min)
  190. } else if (lambda == "lambda.1se") {
  191. lambda.fit <- cv.fit.correct$lambda.1se
  192. coef.correct <- coef(cv.fit.correct, s = cv.fit.correct$lambda.1se)
  193. } else {
  194. stop("lambda must be either 'lambda.min' or 'lambda.1se'")
  195. }
  197. beta <- beta1 + coef.correct
  198. obj <- list(
  199. beta = beta,
  200. family = family,
  201. fitting.list = list(beta1 = beta1, delta_a = coef.correct),
  202. lambda = lambda.fit
  203. )
  204. class(obj) <- "glmtrans"
  205. return(obj)
  206. }
  208. estimate_doubleml_glmnet_ranger <- function(y.t, z.t, x.t) {
  209. data_dml <- DoubleMLData$new(
  210. data = data.frame(y = y.t, d = z.t, x.t),
  211. y_col = "y",
  212. d_cols = "d"
  213. )
  214. ml_g <- lrn("regr.cv_glmnet")
  215. ml_m <- lrn("classif.ranger", predict_type = "prob")
  216. dml_irm <- DoubleMLIRM$new(data_dml, ml_g = ml_g, ml_m = ml_m, n_folds = 3)
  217. dml_irm$fit()
  218. as.numeric(dml_irm$coef)
  219. }
  221. estimate_doubleml_svm_ranger <- function(y.t, z.t, x.t) {
  222. data_dml <- DoubleMLData$new(
  223. data = data.frame(y = y.t, d = z.t, x.t),
  224. y_col = "y",
  225. d_cols = "d"
  226. )
  227. ml_g <- lrn("regr.svm")
  228. ml_m <- lrn("classif.ranger", predict_type = "prob")
  229. dml_irm <- DoubleMLIRM$new(data_dml, ml_g = ml_g, ml_m = ml_m, n_folds = 3)
  230. dml_irm$fit()
  231. as.numeric(dml_irm$coef)
  232. }
  234. estimate_tmle_glmnet <- function(y.t, z.t, x.t) {
  235. tmleresult <- tmle(
  236. y.t, z.t, x.t,
  237. family = "gaussian",
  238. Q.SL.library = "SL.glmnet",
  239. g.SL.library = "SL.glmnet"
  240. )
  241. as.numeric(tmleresult$estimates$ATE$psi)
  242. }
  244. run_safe_method <- function(method_name, expr) {
  245. method.start <- proc.time()[3]
  246. ans <- tryCatch(expr, error = function(e) NA_real_)
  247. record_method_time(method_name, proc.time()[3] - method.start)
  248. ans
  249. }
  251. init_boot_cluster <- function() {
  252. if (BOOT_CORES <= 1) return(NULL)
  254. cl <- parallel::makeCluster(BOOT_CORES)
  255. boot_wd <- getwd()
  256. parallel::clusterCall(cl, setwd, boot_wd)
  257. parallel::clusterEvalQ(cl, {
  258. suppressPackageStartupMessages({
  259. library(MASS)
  260. library(glmnet)
  261. library(Matrix)
  262. library(DoubleML)
  263. library(mlr3)
  264. library(mlr3learners)
  265. library(tmle)
  266. library(ranger)
  267. library(e1071)
  268. library(SuperLearner)
  269. })
  270. if (requireNamespace("lgr", quietly = TRUE)) {
  271. lgr::get_logger("mlr3")$set_threshold("warn")
  272. }
  273. NULL
  274. })
  275. parallel::clusterExport(
  276. cl,
  277. varlist = c(
  278. "m", "n", "p", "METHOD_NAMES",
  279. "alpha1", "alpha2", "beta1.t", "beta1.c", "beta2.t", "beta2.c",
  280. "sigmoid", "dr", "record_method_time", "add_reg_diag", "estimator",
  281. "bootstrap_one", "glmtrans2", "glmtrans2old",
  282. "estimate_doubleml_glmnet_ranger", "estimate_doubleml_svm_ranger",
  283. "estimate_tmle_glmnet", "run_safe_method"
  284. ),
  285. envir = .GlobalEnv
  286. )
  287. parallel::clusterSetRNGStream(cl, 321)
  288. cat("Initialized PSOCK cluster with", BOOT_CORES, "workers\n")
  289. flush.console()
  290. cl
  291. }
  293. ######################################
  294. ## estimator 主函数
  295. ######################################
  296. estimator <- function(data_s,data_t){
  297. ATE <- rep(NA, length(METHOD_NAMES))
  298. names(ATE) <- METHOD_NAMES
  299. y.s<-data_s[,1]; z.s<-data_s[,2]; x.s<-data_s[,c(-1,-2)]
  300. data_s<-data_s[,-3]; colnames(data_s)[1:2]<-c('y','z')
  302. y.t<-data_t[,1]; z.t<-data_t[,2]; x.t<-data_t[,c(-1,-2)]
  303. data_t<-data_t[,-3]; colnames(data_t)[1:2]<-c('y','z')
  305. ATE["naive_diff"] <- run_safe_method(
  306. "naive_diff",
  307. mean(y.t[z.t==1]) - mean(y.t[z.t==0])
  308. )
  310. ATE["oracle_dr"] <- run_safe_method(
  311. "oracle_dr",
  312. dr(y.t,z.t,x.t,alpha2,beta2.t,beta2.c)
  313. )
  315. if(sum(z.t)==n||sum(z.t)==0) stop("Target treatment has only one class")
  317. # target lasso, estimate alpha beta.t beta.c
  318. x.t1=x.t[z.t==1,]
  319. y.t1=y.t[z.t==1]
  320. x.t0=x.t[z.t==0,]
  321. y.t0=y.t[z.t==0]
  322. x.s1=x.s[z.s==1,]
  323. y.s1=y.s[z.s==1]
  324. x.s0=x.s[z.s==0,]
  325. y.s0=y.s[z.s==0]
  327. ATE["target_lasso_dr"] <- run_safe_method("target_lasso_dr", {
  328. lamda<-cv.glmnet(as.matrix(x.t), as.matrix(z.t), family = "binomial", alpha = 1,
  329. type.measure="deviance", nfolds = min(5, length(z.t)), intercept = T, grouped=F)
  330. target.ps <- glmnet(x.t, z.t, family = "binomial", alpha = 1, intercept = F, lambda = lamda$lambda.1se)
  331. lamda<-cv.glmnet(as.matrix(x.t0), as.matrix(y.t0), family = "gaussian", alpha = 1,
  332. type.measure = "mse", nfolds = min(5, length(y.t0)), grouped=F)
  333. target.or.c <- glmnet(x.t0, y.t0, family = "gaussian", alpha = 1, intercept = TRUE, lambda = lamda$lambda.1se)
  334. lamda<-cv.glmnet(as.matrix(x.t1), as.matrix(y.t1), family = "gaussian", alpha = 1,
  335. type.measure = "mse", nfolds = min(5, length(y.t1)), grouped=F)
  336. target.or.t <- glmnet(x.t1, y.t1, family = "gaussian", alpha = 1, intercept = TRUE, lambda = lamda$lambda.1se)
  338. dr(y.t,z.t,x.t,target.ps$beta,target.or.t$beta,target.or.c$beta)
  339. })
  341. penalty.factor = rep(1,p)
  342. ps.target<-data_t[,-1]
  343. colnames(ps.target)[1]<-"y"
  344. ps.source<-data_s[,-1]
  345. colnames(ps.source)[1]<-"y"
  347. ATE["doubleml_glmnet_ranger"] <- run_safe_method(
  348. "doubleml_glmnet_ranger",
  349. estimate_doubleml_glmnet_ranger(y.t, z.t, x.t)
  350. )
  352. ATE["doubleml_svm_ranger"] <- run_safe_method(
  353. "doubleml_svm_ranger",
  354. estimate_doubleml_svm_ranger(y.t, z.t, x.t)
  355. )
  357. ATE["tmle_glmnet"] <- run_safe_method(
  358. "tmle_glmnet",
  359. estimate_tmle_glmnet(y.t, z.t, x.t)
  360. )
  362. ATE["glmtrans2old"] <- run_safe_method("glmtrans2old", {
  363. old_alpha<-glmtrans2old(target=ps.target, source = ps.source, family = "binomial", lambda = "lambda.1se", type.measure="deviance",
  364. alpha = 1, nfolds = min(5,nrow(ps.target)), intercept=F, penalty.factor=penalty.factor)
  366. or.t.target<-data.frame(cbind(y.t1,x.t1[,-1]))
  367. names(or.t.target)[1]<-c("y")
  368. or.t.source<-data.frame(cbind(y.s1,x.s1[,-1]))
  369. names(or.t.source)[1]<-c("y")
  370. old_beta.t <- glmtrans2old(target = or.t.target, source = or.t.source, family = "gaussian", lambda = "lambda.1se", type.measure="mse",
  371. alpha = 1, nfolds = min(5,length(y.t1)), intercept=T, penalty.factor=penalty.factor)
  373. or.c.target<-data.frame(cbind(y.t0,x.t0[,-1]))
  374. names(or.c.target)[1]<-c("y")
  375. or.c.source<-data.frame(cbind(y.s0,x.s0[,-1]))
  376. names(or.c.source)[1]<-c("y")
  377. old_beta.c <- glmtrans2old(target = or.c.target, source = or.c.source, family = "gaussian", lambda = "lambda.1se", type.measure="mse",
  378. alpha = 1, nfolds = min(5,length(y.t0)), intercept=T, penalty.factor=penalty.factor)
  380. dr(y.t,z.t,x.t,old_alpha$beta,old_beta.t$beta,old_beta.c$beta)
  381. })
  383. method.start <- proc.time()[3]
  384. tran_alpha3<-glmtrans2(target=ps.target, source = ps.source, family = "binomial", lambda = "lambda.1se", type.measure="deviance",
  385. alpha = 1, nfolds = min(5,nrow(ps.target)), intercept=F, penalty.factor=penalty.factor)
  387. or.t.target<-data.frame(cbind(y.t1,x.t1[,-1]))
  388. names(or.t.target)[1]<-c("y")
  389. or.t.source<-data.frame(cbind(y.s1,x.s1[,-1]))
  390. names(or.t.source)[1]<-c("y")
  391. tran_beta3.t <- glmtrans2(target = or.t.target, source = or.t.source, family = "gaussian", lambda = "lambda.1se", type.measure="mse",
  392. alpha = 1, nfolds = min(5,length(y.t1)), intercept=T, penalty.factor=penalty.factor)
  394. or.c.target<-data.frame(cbind(y.t0,x.t0[,-1]))
  395. names(or.c.target)[1]<-c("y")
  396. or.c.source<-data.frame(cbind(y.s0,x.s0[,-1]))
  397. names(or.c.source)[1]<-c("y")
  398. tran_beta3.c <- glmtrans2(target = or.c.target, source = or.c.source, family = "gaussian", lambda = "lambda.1se", type.measure="mse",
  399. alpha = 1, nfolds = min(5,length(y.t0)), intercept=T, penalty.factor=penalty.factor)
  401. if (SAVE_REG_DIAG) {
  402. REG_DIAG <<- rbind(
  403. REG_DIAG,
  404. add_reg_diag(tran_alpha3, "DSTCL", "PS"),
  405. add_reg_diag(tran_beta3.t, "DSTCL", "OR_treated"),
  406. add_reg_diag(tran_beta3.c, "DSTCL", "OR_control")
  407. )
  408. }
  410. ATE["dstcl_transfer"] <- dr(y.t,z.t,x.t,tran_alpha3$beta,tran_beta3.t$beta,tran_beta3.c$beta)
  411. record_method_time("dstcl_transfer", proc.time()[3] - method.start)
  412. return(ATE)
  413. }
  415. bootstrap_one <- function(i) {
  416. SAVE_REG_DIAG <<- FALSE
  417. CURRENT_STAGE <<- "bootstrap"
  418. CURRENT_SIM <<- sim_id_boot
  419. MAIN_METHOD_TIME_SUM <<- matrix(0, 1, length(METHOD_NAMES))
  420. MAIN_METHOD_TIME_COUNT <<- matrix(0, 1, length(METHOD_NAMES))
  421. BOOT_METHOD_TIME_SUM <<- matrix(0, sim_id_boot, length(METHOD_NAMES))
  422. BOOT_METHOD_TIME_COUNT <<- matrix(0, sim_id_boot, length(METHOD_NAMES))
  423. colnames(BOOT_METHOD_TIME_SUM) <<- METHOD_NAMES
  424. colnames(BOOT_METHOD_TIME_COUNT) <<- METHOD_NAMES
  426. data_ss <- data_s_boot[sample(1:m,m,replace=TRUE),]
  427. data_tt <- data_t_boot[sample(1:n,n,replace=TRUE),]
  428. fit.i <- try(estimator(data_ss,data_tt), silent = TRUE)
  430. if (inherits(fit.i, "try-error")) {
  431. list(
  432. index = i,
  433. ate = rep(NA, length(METHOD_NAMES)),
  434. time_sum = rep(0, length(METHOD_NAMES)),
  435. time_count = rep(0, length(METHOD_NAMES)),
  436. error = as.character(fit.i)
  437. )
  438. } else {
  439. list(
  440. index = i,
  441. ate = fit.i,
  442. time_sum = as.numeric(BOOT_METHOD_TIME_SUM[sim_id_boot, ]),
  443. time_count = as.numeric(BOOT_METHOD_TIME_COUNT[sim_id_boot, ]),
  444. error = NA_character_
  445. )
  446. }
  447. }
  449. ######################################
  450. ## bootstrap 函数
  451. ######################################
  452. boo <- function(data_s,data_t,B){
  453. ATE_B <- matrix(NA,B,length(METHOD_NAMES))
  454. old.save <- SAVE_REG_DIAG
  455. old.stage <- CURRENT_STAGE
  456. SAVE_REG_DIAG <<- FALSE
  457. CURRENT_STAGE <<- "bootstrap"
  458. on.exit({
  459. SAVE_REG_DIAG <<- old.save
  460. CURRENT_STAGE <<- old.stage
  461. }, add = TRUE)
  463. sim_id <- CURRENT_SIM
  464. cat(" Bootstrap 1 to", B, "of", B,
  465. "in simulation", sim_id, "using", BOOT_CORES, "workers with load balancing\n")
  466. flush.console()
  468. if (is.null(BOOT_CLUSTER)) {
  469. data_s_boot <<- data_s
  470. data_t_boot <<- data_t
  471. sim_id_boot <<- sim_id
  472. boot_res <- lapply(seq_len(B), bootstrap_one)
  473. rm(data_s_boot, data_t_boot, sim_id_boot, envir = .GlobalEnv)
  474. } else {
  475. data_s_boot <- data_s
  476. data_t_boot <- data_t
  477. sim_id_boot <- sim_id
  478. parallel::clusterExport(
  479. BOOT_CLUSTER,
  480. varlist = c("data_s_boot", "data_t_boot", "sim_id_boot"),
  481. envir = environment()
  482. )
  483. boot_res <- parallel::parLapplyLB(BOOT_CLUSTER, seq_len(B), bootstrap_one)
  484. }
  486. for (res in boot_res) {
  487. ATE_B[res$index, ] <- res$ate
  488. BOOT_METHOD_TIME_SUM[sim_id, ] <<- BOOT_METHOD_TIME_SUM[sim_id, ] + res$time_sum
  489. BOOT_METHOD_TIME_COUNT[sim_id, ] <<- BOOT_METHOD_TIME_COUNT[sim_id, ] + res$time_count
  490. }
  491. return(ATE_B)
  492. }
  494. ######################################
  495. ## 主模拟循环
  496. ######################################
  497. EST <- matrix(NA,nsim,length(METHOD_NAMES))
  498. cover<- matrix(NA,nsim,length(METHOD_NAMES))
  499. colnames(EST) <- METHOD_NAMES
  500. colnames(cover) <- METHOD_NAMES
  501. True_ate<-rep(-2,length(METHOD_NAMES))
  503. BOOT_CLUSTER <- init_boot_cluster()
  504. on.exit({
  505. if (!is.null(BOOT_CLUSTER)) {
  506. parallel::stopCluster(BOOT_CLUSTER)
  507. }
  508. }, add = TRUE)
  510. for(k in 1:nsim){
  511. sim.start <- proc.time()[3]
  512. cat("Starting simulation", k, "of", nsim, "\n")
  513. flush.console()
  515. data<-gen.transfer.train(m,n,p,alpha1,alpha2,beta1.t,beta1.c,beta2.t,beta2.c,scenario="ps_correct")
  516. data_s<-as.data.frame(data[["data.source"]])
  517. data_t<-as.data.frame(data[["data.target"]])
  519. CURRENT_SIM <<- k
  520. CURRENT_STAGE <<- "main"
  521. est.start <- proc.time()[3]
  522. fit.k <- try(estimator(data_s,data_t), silent = TRUE)
  523. ESTIMATOR_TIME[k] <- proc.time()[3] - est.start
  525. if (inherits(fit.k, "try-error")) {
  526. cat("Simulation", k, "main estimator failed:\n", as.character(fit.k), "\n")
  527. flush.console()
  528. EST[k,] <- NA
  529. RUN_STATUS <- rbind(
  530. RUN_STATUS,
  531. data.frame(sim = k, status = "failed", message = as.character(fit.k),
  532. estimator_runtime_sec = ESTIMATOR_TIME[k],
  533. bootstrap_runtime_sec = NA,
  534. total_runtime_sec = NA)
  535. )
  536. } else {
  537. EST[k,] <- fit.k
  538. RUN_STATUS <- rbind(
  539. RUN_STATUS,
  540. data.frame(sim = k, status = "ok", message = NA,
  541. estimator_runtime_sec = ESTIMATOR_TIME[k],
  542. bootstrap_runtime_sec = NA,
  543. total_runtime_sec = NA)
  544. )
  545. }
  547. boot.start <- proc.time()[3]
  548. ATE_B<-boo(data_s,data_t,B)
  549. BOOTSTRAP_TIME[k] <- proc.time()[3] - boot.start
  550. upper<-apply(ATE_B,2,function(x)quantile(x,0.975,na.rm=TRUE))
  551. lower<-apply(ATE_B,2,function(x)quantile(x,0.025,na.rm=TRUE))
  552. cover[k,]<-ifelse(((True_ate-lower)>-1e-06)&((True_ate-upper)<1e-06),1,0)
  554. TOTAL_TIME[k] <- proc.time()[3] - sim.start
  555. RUN_STATUS[RUN_STATUS$sim == k, "bootstrap_runtime_sec"] <- BOOTSTRAP_TIME[k]
  556. RUN_STATUS[RUN_STATUS$sim == k, "total_runtime_sec"] <- TOTAL_TIME[k]
  558. cat("Completed", k, "of", nsim,
  559. "simulations; estimator:", round(ESTIMATOR_TIME[k], 2), "sec;",
  560. "bootstrap:", round(BOOTSTRAP_TIME[k], 2), "sec;",
  561. "total:", round(TOTAL_TIME[k], 2), "sec\n")
  563. sim_method_total <- BOOT_METHOD_TIME_SUM[k, ]
  564. sim_method_count <- BOOT_METHOD_TIME_COUNT[k, ]
  565. sim_method_mean <- ifelse(sim_method_count > 0, sim_method_total / sim_method_count, NA)
  566. top_methods <- order(sim_method_total, decreasing = TRUE)[1:min(3, length(METHOD_NAMES))]
  567. cat(" Slowest bootstrap methods in simulation", k, ":\n")
  568. for (idx in top_methods) {
  569. if (sim_method_count[idx] > 0) {
  570. cat(" -", METHOD_NAMES[idx],
  571. "total:", round(sim_method_total[idx], 2), "sec;",
  572. "mean:", round(sim_method_mean[idx], 2), "sec;",
  573. "n:", sim_method_count[idx], "\n")
  574. }
  575. }
  576. flush.console()
  577. }
  579. ######################################
  580. ## 结果整理与保存
  581. ######################################
  582. True_ate <- -2
  583. BIAS<-colMeans(EST,na.rm=TRUE)-True_ate
  584. SD<-apply(EST,2,sd,na.rm=TRUE)
  585. MSE<-apply(EST,2,function(x)mean((x-True_ate)^2,na.rm=TRUE))
  586. Coverage<-colMeans(cover,na.rm=TRUE)
  588. result<-data.frame(
  589. "Method"=METHOD_NAMES,
  590. "Bias"=round(BIAS,2),
  591. "SD"=round(SD,2),
  592. "MSE"=round(MSE,2),
  593. "Coverage"=round(Coverage,2)
  594. )
  596. print(result)
  598. runtime_summary <- data.frame(
  599. mean_estimator_runtime_sec = mean(ESTIMATOR_TIME, na.rm = TRUE),
  600. median_estimator_runtime_sec = median(ESTIMATOR_TIME, na.rm = TRUE),
  601. mean_bootstrap_runtime_sec = mean(BOOTSTRAP_TIME, na.rm = TRUE),
  602. median_bootstrap_runtime_sec = median(BOOTSTRAP_TIME, na.rm = TRUE),
  603. mean_total_runtime_sec = mean(TOTAL_TIME, na.rm = TRUE),
  604. median_total_runtime_sec = median(TOTAL_TIME, na.rm = TRUE),
  605. sd_total_runtime_sec = sd(TOTAL_TIME, na.rm = TRUE),
  606. min_total_runtime_sec = min(TOTAL_TIME, na.rm = TRUE),
  607. max_total_runtime_sec = max(TOTAL_TIME, na.rm = TRUE),
  608. n_failed = sum(RUN_STATUS$status != "ok"),
  609. failure_rate = mean(RUN_STATUS$status != "ok")
  610. )
  612. main_method_total <- colSums(MAIN_METHOD_TIME_SUM, na.rm = TRUE)
  613. main_method_count <- colSums(MAIN_METHOD_TIME_COUNT, na.rm = TRUE)
  614. bootstrap_method_total <- colSums(BOOT_METHOD_TIME_SUM, na.rm = TRUE)
  615. bootstrap_method_count <- colSums(BOOT_METHOD_TIME_COUNT, na.rm = TRUE)
  617. method_time_summary <- data.frame(
  618. method_index = seq_along(METHOD_NAMES),
  619. method = METHOD_NAMES,
  620. main_total_sec = main_method_total,
  621. main_count = main_method_count,
  622. main_mean_sec = ifelse(main_method_count > 0, main_method_total / main_method_count, NA),
  623. bootstrap_total_sec = bootstrap_method_total,
  624. bootstrap_count = bootstrap_method_count,
  625. bootstrap_mean_sec = ifelse(bootstrap_method_count > 0, bootstrap_method_total / bootstrap_method_count, NA)
  626. )
  628. bootstrap_method_time_by_sim <- do.call(
  629. rbind,
  630. lapply(1:nsim, function(sim_id) {
  631. total <- BOOT_METHOD_TIME_SUM[sim_id, ]
  632. count <- BOOT_METHOD_TIME_COUNT[sim_id, ]
  633. data.frame(
  634. sim = sim_id,
  635. method_index = seq_along(METHOD_NAMES),
  636. method = METHOD_NAMES,
  637. bootstrap_total_sec = as.numeric(total),
  638. bootstrap_count = as.numeric(count),
  639. bootstrap_mean_sec = ifelse(count > 0, as.numeric(total) / as.numeric(count), NA)
  640. )
  641. })
  642. )
  644. if (nrow(REG_DIAG) > 0) {
  645. reg_summary <- aggregate(
  646. cbind(nfolds_used,
  647. nonzero_source, nonzero_delta, nonzero_final,
  648. prop_nonzero_source, prop_nonzero_delta, prop_nonzero_final,
  649. l1_source, l1_delta, l1_final,
  650. source_lambda, correction_lambda,
  651. n_try_correction) ~ method + component,
  652. data = REG_DIAG,
  653. FUN = function(x) mean(x, na.rm = TRUE),
  654. na.action = na.pass
  655. )
  657. reg_lambda_summary <- aggregate(
  658. lambda_rule ~ method + component,
  659. data = REG_DIAG,
  660. FUN = function(x) {
  661. x <- unique(na.omit(x))
  662. if (length(x) == 0) NA_character_ else paste(x, collapse = ";")
  663. },
  664. na.action = na.pass
  665. )
  667. reg_summary <- merge(
  668. reg_summary,
  669. reg_lambda_summary,
  670. by = c("method", "component"),
  671. all.x = TRUE
  672. )
  673. } else {
  674. reg_summary <- data.frame()
  675. }
  677. write.csv(result,file = file.path(output_dir, "result_main_ps_50_2000_200_6.csv"))
  678. write.csv(EST,file = file.path(output_dir, "EST_main_ps_50_2000_200_6.csv"))
  679. write.csv(cover,file = file.path(output_dir, "cover_main_ps_50_2000_200_6.csv"))
  680. write.csv(REG_DIAG, file = file.path(output_dir, "regularization_diagnostics_raw_main_ps_50_2000_200_6.csv"), row.names = FALSE)
  681. write.csv(reg_summary, file = file.path(output_dir, "regularization_diagnostics_summary_main_ps_50_2000_200_6.csv"), row.names = FALSE)
  682. write.csv(RUN_STATUS, file = file.path(output_dir, "run_status_main_ps_50_2000_200_6.csv"), row.names = FALSE)
  683. write.csv(runtime_summary, file = file.path(output_dir, "runtime_summary_main_ps_50_2000_200_6.csv"), row.names = FALSE)
  684. write.csv(method_time_summary, file = file.path(output_dir, "method_time_summary_main_ps_50_2000_200_6.csv"), row.names = FALSE)
  685. write.csv(bootstrap_method_time_by_sim, file = file.path(output_dir, "bootstrap_method_time_main_ps_50_2000_200_6.csv"), row.names = FALSE)

simulate_dstcl_main.R at commit cca9e24, no license · at the source

Overview

Authors: Yunxin Shi1, Lulu Pan1, Yu Hu1, Yongfu Yu1, Guoyou Qin1
ORCID iDs: Yunxin Shi, Guoyou Qin
  1. Department of Biostatistics, Key Laboratory of Public Health Safety of Ministry of Education, NHC Key Laboratory for Health Technology Assessment, School of Public Health, Fudan University, Shanghai, China
Institutions: Fudan University (China)
Journal: PLoS computational biology, volume 22, issue 8, article e1014706
Dates: received 20 February 2026; accepted 11 August 2026; published online 24 August 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pcbi.1014706 · PMID 42636277 · PMCID PMC13533430 · OpenAlex W7204095745
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), Alzheimer's / dementia (population)
Methods: Machine learning, Statistics
MeSH: Alzheimer Disease*, Amyloid beta-Peptides, Biomarkers, Black or African American, Computational Biology, Computer Simulation, Humans, Neuroimaging, Peptide Fragments (* major topic)
Journal subjects: People and places, Population groupings, Ethnicities, African American people, Medicine and Health Sciences, Mental Health and Psychiatry, Dementia, Alzheimer's Disease, Neurology, Medical Conditions, Neurodegenerative Diseases, Research and Analysis Methods, Simulation and Modeling, Computer and Information Sciences, Artificial Intelligence, Machine Learning, Biology and Life Sciences, Anatomy, Body Fluids, Cerebrospinal Fluid, Physiology, Nervous System, Cell biology, Chromosome biology, Chromatin, Chromatin modification, DNA methylation, Genetics, Epigenetics, Gene expression, DNA, DNA modification, Biochemistry, Nucleic acids, Biomarkers, Health Care, Health Care Policy, Treatment Guidelines
Topic: Bayesian Modeling and Causal Inference (Artificial Intelligence, Computer Science), according to OpenAlex
Funding: National Natural Science Foundation of China (82473724, 82273730); Shanghai Municipal Science and Technology Major Project (ZD2021CY001); Shanghai Talent Programs (BJKJ2024050); Shanghai Education Development Foundation and Shanghai Municipal Education Commission
Citations: not cited yet (Europe PMC); 27 references in the paper

Abstract

Amyloid-Beta 42 (ABeta42) is a key biomarker of cerebral amyloidosis in Alzheimer’s disease, and estimating its causal effect on subsequent cognitive outcomes is important for understanding disease progression in racial minority populations, such as African Americans. However, causal effect analyses in such populations are often challenged by limited sample sizes, leading to biased or unstable estimates when only target-population data are used. Existing transfer causal learning methods borrow information from a larger related source domain, but may be unstable when source-domain estimation is high-dimensional or when irrelevant covariates introduce noise into transferred nuisance models. We propose a double shrinkage transfer causal learning (DSTCL) estimator that regularizes both the initial source domain estimation and the source-target parameter differences, and then estimates causal effects using a doubly robust framework. We evaluated DSTCL in simulations under varying source sample sizes, inter-domain similarities, covariate dimensions, and other settings, and then applied it to the Alzheimer’s Disease Neuroimaging Initiative dataset using African American participants as the target domain and non-Hispanic White participants as the source domain. In simulations, DSTCL achieved small bias, low mean squared error, and stable confidence interval coverage, generally outperforming competing methods. In the ADNI application, lower baseline ABeta42, reflecting greater amyloid burden, was associated with worse 12-month cognitive performance among African American participants, with DSTCL producing the narrowest confidence interval among the methods compared. These findings suggest that DSTCL provides a practical framework for transfer causal inference in biomedical studies with limited target population.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 5 matches between paragraphs and lines of code.

LESWinnie/DSTCL

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: cca9e2446b942bf93aad7eeeecc96ec807bee1fc, 13 July 2026
Languages: R (4)
Size: 5 files, 4 scripts
Software Heritage: not archived
Found in: “Data Availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: glmnet (2 files), tidyverse (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
5 files
At the source: github.com/LESWinnie/DSTCL

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 4 scripts, each with its path and the digest of its content;
  • 5 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability

The R code for implementing DSTCL, generating the simulated data, and reproducing the analyses is publicly available on GitHub at https://github.com/LESWinnie/DSTCL. The simulated datasets can be generated using the code provided in the repository. The data used in the real-data application were obtained from the ADNI database (https://adni.loni.usc.edu/). Under the ADNI Data Use Agreement, participant-level data cannot be redistributed by the authors, but qualified researchers may apply for access at https://adni.loni.usc.edu/data-samples/adni-data/. The authors had no special access privileges.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 9 MeSH terms, 4 funders, 23 references.

Cite

This paper

Shi, Y., Pan, L., Hu, Y., Yu, Y., & Qin, G. (2026). Double shrinkage transfer causal learning: An application to alzheimer's disease. PLoS computational biology, 22(8), e1014706. https://doi.org/10.1371/journal.pcbi.1014706

BibTeX

@article{shi2026double,
author = {Shi, Yunxin and Pan, Lulu and Hu, Yu and Yu, Yongfu and Qin, Guoyou},
title = {{Double shrinkage transfer causal learning: An application to alzheimer's disease}},
journal = {PLoS computational biology},
year = {2026},
month = aug,
volume = {22},
number = {8},
pages = {e1014706},
publisher = {PLOS},
issn = {1553-734X},
doi = {10.1371/journal.pcbi.1014706},
url = {https://doi.org/10.1371/journal.pcbi.1014706},
pmid = {42636277},
pmcid = {PMC13533430}
}

RIS

TY - JOUR
AU - Shi, Yunxin
AU - Pan, Lulu
AU - Hu, Yu
AU - Yu, Yongfu
AU - Qin, Guoyou
TI - Double shrinkage transfer causal learning: An application to alzheimer's disease
T2 - PLoS computational biology
J2 - PLoS Comput Biol
PY - 2026
DA - 2026/08/24
VL - 22
IS - 8
SP - e1014706
SN - 1553-734X
PB - PLOS
DO - 10.1371/journal.pcbi.1014706
UR - https://doi.org/10.1371/journal.pcbi.1014706
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pcbi.1014706",
"type": "article-journal",
"title": "Double shrinkage transfer causal learning: An application to alzheimer's disease",
"container-title": "PLoS computational biology",
"author": [
{
"family": "Shi",
"given": "Yunxin"
},
{
"family": "Pan",
"given": "Lulu"
},
{
"family": "Hu",
"given": "Yu"
},
{
"family": "Yu",
"given": "Yongfu"
},
{
"family": "Qin",
"given": "Guoyou"
}
],
"container-title-short": "PLoS Comput Biol",
"volume": "22",
"issue": "8",
"page": "e1014706",
"DOI": "10.1371/journal.pcbi.1014706",
"PMID": "42636277",
"PMCID": "PMC13533430",
"ISSN": "1553-734X",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pcbi.1014706",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
24
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1186/s13073-026-01698-8 [code]
From aging to Alzheimer's disease: concordant brain DNA methylation changes in late life.
Journal: Genome medicine
In common: tidyverse, Alzheimer's / dementia, 3 references
[2] doi:10.3390/diagnostics16132029
FLAME: Federated Learning and Aggregated Multi-Model Ensemble for Multi-Class Alzheimer's Disease Stage Classification from Structured Clinical Data.
Journal: Diagnostics (Basel, Switzerland)
In common: Alzheimer's / dementia, 3 references
[3] doi:10.1038/s41586-026-10877-x [code]
Human brain organoids record the passage of time over multiple years.
Journal: Nature
In common: glmnet, tidyverse, 1 reference
[4] doi:10.1016/j.isci.2026.116439 [code]
Decoding the role of transcriptomic clocks in the human prefrontal cortex.
Journal: iScience
In common: glmnet, tidyverse, 1 reference
[5] doi:10.1038/s41380-026-03686-1 [code]
Early oligodendrocyte dysfunction signature in Alzheimer's disease: Insights from DNA methylomics and transcriptomics.
Journal: Molecular psychiatry
In common: tidyverse, Alzheimer's / dementia, 2 references
[6] doi:10.3389/fnagi.2026.1847611 [code]
APOE ε4-associated hippocampal atrophy trajectories across the Alzheimer's disease continuum: a systematic review, meta-analysis, and longitudinal validation.
Journal: Frontiers in aging neuroscience
In common: tidyverse, Alzheimer's / dementia, 2 references
[7] doi:10.1080/20002297.2026.2705667 [code]
Oral microbiota dysbiosis related to the cortical thinning and cognitive impairment in cerebral small vessel disease.
Journal: Journal of oral microbiology
In common: glmnet, tidyverse, Alzheimer's / dementia
[8] doi:10.1177/13872877261471049 [code]
Predicting future brain atrophy based on longitudinal MRI.
Journal: Journal of Alzheimer's disease : JAD
In common: glmnet, tidyverse, Alzheimer's / dementia
[9] doi:10.1002/advs.76205 [code]
CHCHD10 Mitigates Alzheimer's Disease-Related Phenotypes in Association With Epigenetic Remodeling in Directly Reprogrammed Neurons.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: glmnet, tidyverse, Alzheimer's / dementia
[10] doi:10.1038/s43587-026-01149-4 [code]
The blood metabolome of brain health in midlife and influences of genes, microbiome and exposome.
Journal: Nature aging
In common: glmnet, tidyverse, Alzheimer's / dementia

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.