From bc9622de8ecc956f71c200cfcb81cbfc1c0c09d5 Mon Sep 17 00:00:00 2001 From: "George G. Vega Yon" Date: Tue, 10 Feb 2026 16:49:44 -0700 Subject: [PATCH 01/19] Adding the README.md file --- README.md | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 README.md diff --git a/README.md b/README.md new file mode 100644 index 0000000..4ca5063 --- /dev/null +++ b/README.md @@ -0,0 +1,3 @@ +## Comprehensive Cohort Project + +This is an example of a commit From a27f064d9adfa7ee6388db16d90955dee76394ca Mon Sep 17 00:00:00 2001 From: Jennifer Xu <48976310+jenniferxu1101@users.noreply.github.com> Date: Wed, 11 Feb 2026 10:59:08 -0500 Subject: [PATCH 02/19] Add files via upload --- SIDR_Ravinew.R | 303 +++++++++++++++++++++ SIDRnew.R | 405 +++++++++++++++++++++++++++++ SensIAT_sim_outcome_modeler_mave.R | 174 +++++++++++++ 3 files changed, 882 insertions(+) create mode 100644 SIDR_Ravinew.R create mode 100644 SIDRnew.R create mode 100644 SensIAT_sim_outcome_modeler_mave.R diff --git a/SIDR_Ravinew.R b/SIDR_Ravinew.R new file mode 100644 index 0000000..b71b5b2 --- /dev/null +++ b/SIDR_Ravinew.R @@ -0,0 +1,303 @@ +SIDR_Ravinew <- function(X, Y, + Y.CP = NULL, + initial = NULL, + kernel = "dnorm", + method = "optim", + optim_method = "BFGS", + abs.tol = 1e-4, + bandwidth = NULL, + wi.boot = NULL, + index_ID) +{ + X <- as.matrix(X) + Y <- as.matrix(Y) + + number_n <- dim(X)[1] + number_p <- dim(X)[2] + + if (is.null(initial)) + { + initial <- c(1, rep(0, number_p-1)) + }else + { + initial <- as.vector(initial) + if(initial[1] != 0) + initial <- initial/initial[1] + } + + if (is.null(bandwidth)) + { + if (kernel=="K2_Biweight") + { + if (is.null(wi.boot)) + { + Eij3 <- function(parameter){ + K <- function(x, h) 15/16*(1-(x/h)^2)^2 * (abs(x) <= h) + + b <- c(1, parameter[1:(number_p-1)]) + h <- exp(parameter[number_p]) + + x <- c(X%*%b) + y <- Y + + n <- length(y) + yo <- order(y) + ys <- y[yo] + uy <- rle(ys)[[1]] + cols <- cumsum(uy) # Ravi + ei <- rep(0, n) + for (i in 1:n){ + Kih <- K(x-x[i],h=h) + + # remove all obs of ith obs's patient + index_remove <- which(index_ID == index_ID[i]) + Kih[index_remove] <- 0 + + #Kih[i] <- 0 # the fix + denom <- sum(Kih) + ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + } + return(sum(ei)/n^2) + } + }else + { + stop("There's no weighted version of the K2_Biweight kernel.") + } + }else if (kernel == "dnorm") + { + if (is.null(wi.boot)) + { + Eij3 <- function(parameter){ + K <- function(x, h) dnorm(x/h, 0, 1) + + b <- c(1, parameter[1:(number_p-1)]) + h <- exp(parameter[number_p]) + + x <- c(X%*%b) + y <- Y + + n <- length(y) + yo <- order(y) + ys <- y[yo] + uy <- rle(ys)[[1]] + cols <- cumsum(uy) + ei <- rep(0, n) + for (i in 1:n){ + Kih <- K(x-x[i],h=h) + + # remove all obs of ith obs's patient + index_remove <- which(index_ID == index_ID[i]) + Kih[index_remove] <- 0 + + # Kih[i] <- 0 # the fix + denom <- sum(Kih) + ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + } + return(sum(ei)/n^2) + } + }else + { + stop("There's no weighted version of the dnorm kernel.") + } + # }else if (kernel=="K4_Biweight") + # { + # if (is.null(wi.boot)) + # { + # Eij3 <- function(parameter){ + # K <- function(x, h) 105/64*(1-3*((x/h)^2))*(1-(x/h)^2)^2 * (abs(x) <= h) + # + # b <- c(1, parameter[1:(number_p-1)]) + # h <- exp(parameter[number_p]) + # + # x <- c(X%*%b) + # y <- Y + # + # n <- length(y) + # yo <- order(y) + # ys <- y[yo] + # uy <- rle(ys)[[1]] + # cols <- cumsum(uy) + # ei <- rep(0, n) + # for (i in 1:n){ + # Kih <- K(x-x[i],h=h) + # + # # remove all obs of ith obs's patient + # index_remove <- which(index_ID == index_ID[i]) + # Kih[index_remove] <- 0 + # + # # Kih[i] <- 0 # the fix + # denom <- sum(Kih) + # ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + # } + # return(sum(ei)/n^2) + # } + # }else + # { + # stop("There's no weighted version of the K4_Biweight kernel.") + # } + } else + rlang::abort("The kernel is not supported.") + + if(method == "nlminb") + { + esti <- nlminb(start = c(initial[-1], 0), + objective = Eij3, + control = list(abs.tol = abs.tol)) + }else if (method == "optim") + { + # the new optimize function using optim, you can change the lower and upper + esti <- optim(par = c(initial[-1], 0), + fn = Eij3, + method = optim_method, + control = list(abstol = abs.tol)) + }else if (method == "nmk") + { + assertthat::assert_that(requireNamespace("dfoptim", quietly = TRUE)) + esti <- dfoptim::nmk(par = c(initial[-1], 0), + fn = Eij3, + control = list(tol = abs.tol)) + } + + results <- list(coef = c(1, esti$par[1:(number_p-1)]), + bandwidth = exp(esti$par[number_p]), + details = esti) + }else + { + if (kernel=="K2_Biweight") + { + if (is.null(wi.boot)) + { + Eij3 <- function(parameter){ + K <- function(x, h) 15/16*(1-(x/h)^2)^2 * (abs(x) <= h) + + b <- c(1, parameter[1:(number_p-1)]) + h <- bandwidth + + x <- c(X%*%b) + y <- Y + + n <- length(y) + yo <- order(y) + ys <- y[yo] + uy <- rle(ys)[[1]] + cols <- cumsum(uy) # Ravi + ei <- rep(0, n) + for (i in 1:n){ + Kih <- K(x-x[i],h=h) + + # remove all obs of ith obs's patient + index_remove <- which(index_ID == index_ID[i]) + Kih[index_remove] <- 0 + + # Kih[i] <- 0 # the fix + denom <- sum(Kih) + ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + } + return(sum(ei)/n^2) + } + + }else + { + stop("There's no weighted version of the K2_Biweight kernel.") + } + }else if (kernel=="dnorm") + { + if (is.null(wi.boot)) + { + Eij3 <- function(parameter){ + K <- function(x, h) dnorm(x/h,0,1) + + b <- c(1, parameter[1:(number_p-1)]) + h <- bandwidth + + x <- c(X%*%b) + y <- Y + + n <- length(y) + yo <- order(y) + ys <- y[yo] + uy <- rle(ys)[[1]] + cols <- cumsum(uy) # Ravi + ei <- rep(0, n) + for (i in 1:n){ + Kih <- K(x-x[i],h=h) + + # remove all obs of ith obs's patient + index_remove <- which(index_ID == index_ID[i]) + Kih[index_remove] <- 0 + + # Kih[i] <- 0 # the fix + denom <- sum(Kih) + ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + } + return(sum(ei)/n^2) + } + }else + { + stop("There's no weighted version of the dnorm kernel.") + } + # }else if (kernel=="K4_Biweight") + # { + # if (is.null(wi.boot)) + # { + # Eij3 <- function(parameter){ + # K <- function(x, h) 105/64*(1-3*((x/h)^2))*(1-(x/h)^2)^2 * (abs(x) <= h) + # + # b <- c(1, parameter[1:(number_p-1)]) + # h <- bandwidth + # + # x <- c(X%*%b) + # y <- Y + # + # n <- length(y) + # yo <- order(y) + # ys <- y[yo] + # uy <- rle(ys)[[1]] + # cols <- cumsum(uy) # Ravi + # ei <- rep(0, n) + # for (i in 1:n){ + # Kih <- K(x-x[i],h=h) + # + # # remove all obs of ith obs's patient + # index_remove <- which(index_ID == index_ID[i]) + # Kih[index_remove] <- 0 + # + # # Kih[i] <- 0 # the fix + # denom <- sum(Kih) + # ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + # } + # return(sum(ei)/n^2) + # } + # }else + # { + # stop("There's no weighted version of the K4_Biweight kernel.") + # } + } else + rlang::abort("The kernel is not supported.") + + if(method == "nlminb") + { + esti <- nlminb(start = initial[-1], + objective = Eij3, + control = list(abs.tol = abs.tol)) + }else if (method == "optim") + { + # the new optimize function using optim, you can change the lower and upper + esti <- optim(par = initial[-1], + fn = Eij3, + method = optim_method, + control = list(abstol = abs.tol)) + }else if (method == "nmk") + { + assertthat::assert_that(requireNamespace("dfoptim", quietly = TRUE)) + esti <- dfoptim::nmk(par = initial[-1], + fn = Eij3, + control = list(tol = abs.tol)) + } + results <- list(coef = c(1, esti$par[1:(number_p-1)]), + bandwidth = bandwidth, + details = esti) + } + + return(results) +} diff --git a/SIDRnew.R b/SIDRnew.R new file mode 100644 index 0000000..b01d143 --- /dev/null +++ b/SIDRnew.R @@ -0,0 +1,405 @@ +##### This function defined by MingYueh +SIDRnew <- function(X, Y, + Y.CP = NULL, + initial = NULL, + kernel = "K2_Biweight", + method = c("optim", "nlminb", "nmk"), + optim_method = "BFGS", + abs.tol = 1e-4, + bandwidth = NULL, + wi.boot = NULL) +{ + X <- as.matrix(X) + Y <- as.matrix(Y) + + method <- match.arg(method) + + number_n <- dim(X)[1] + number_p <- dim(X)[2] + + if (is.null(initial)) + { + initial <- c(1, rep(0, number_p-1)) + }else + { + initial <- as.vector(initial) + initial <- initial/initial[1] + } + + if (is.null(bandwidth)) + { + if (kernel=="K2_Biweight") + { + if (is.null(wi.boot)) + { + Eij3 <- function(parameter){ + K <- function(x, h) 15/16*(1-(x/h)^2)^2 * (abs(x) <= h) * (x!=0) + + b <- c(1, parameter[1:(number_p-1)]) + h <- exp(parameter[number_p]) + + x <- c(X%*%b) + y <- Y + + n <- length(y) + yo <- order(y) + ys <- y[yo] + uy <- rle(ys)[[1]] + cols <- cumsum(uy) + ei <- rep(0, n) + for (i in 1:n){ + Kih <- K(x-x[i],h=h) + denom <- sum(Kih) + ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + } + return(sum(ei)/n^2) + } + # cv.bh <- function(parameter) + # { + # b <- c(1, parameter[1:(number_p-1)]) + # h <- exp(parameter[number_p]) + # cv <- mean((Y.CP-NWcv_K2B_rcpp(X = X %*% b, Y = Y.CP, + # h = h))^2) + # return(cv) + # } + }else + { + stop("There's no weighted version of the K2_Biweight kernel.") + } + }else if (kernel == "dnorm") + { + if (is.null(wi.boot)) + { + Eij3 <- function(parameter){ + K <- function(x, h) dnorm(x/h, 0, 1)* (x!=0) + + b <- c(1, parameter[1:(number_p-1)]) + h <- exp(parameter[number_p]) + + x <- c(X%*%b) + y <- Y + + n <- length(y) + yo <- order(y) + ys <- y[yo] + uy <- rle(ys)[[1]] + cols <- cumsum(uy) + ei <- rep(0, n) + for (i in 1:n){ + Kih <- K(x-x[i],h=h) + denom <- sum(Kih) + ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + } + return(sum(ei)/n^2) + } + # cv.bh <- function(parameter) + # { + # b <- c(1, parameter[1:(number_p-1)]) + # h <- exp(parameter[number_p]) + # cv <- mean((Y.CP-NWcv_dnorm_rcpp(X = X %*% b, Y = Y.CP, + # h = h))^2) + # return(cv) + # } + }else + { + # wi.boot <- as.vector(wi.boot) + # cv.bh <- function(parameter) + # { + # b <- c(1, parameter[1:(number_p-1)]) + # h <- exp(parameter[number_p]) + # cv <- mean((Y.CP-NWcv_K2B_w_rcpp(X = X %*% b, Y = Y.CP, + # h = h, w = wi.boot))^2) + # return(cv) + # } + stop("There's no weighted version of the dnorm kernel.") + } + }else if (kernel=="K4_Biweight") + { + if (is.null(wi.boot)) + { + Eij3 <- function(parameter){ + K <- function(x, h) 105/64*(1-3*((x/h)^2))*(1-(x/h)^2)^2 * (abs(x) <= h) * (x!=0) + + b <- c(1, parameter[1:(number_p-1)]) + h <- exp(parameter[number_p]) + + x <- c(X%*%b) + y <- Y + + n <- length(y) + yo <- order(y) + ys <- y[yo] + uy <- rle(ys)[[1]] + cols <- cumsum(uy) + ei <- rep(0, n) + for (i in 1:n){ + Kih <- K(x-x[i],h=h) + denom <- sum(Kih) + ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + } + return(sum(ei)/n^2) + } + # cv.bh <- function(parameter) + # { + # b <- c(1, parameter[1:(number_p-1)]) + # h <- exp(parameter[number_p]) + # cv <- mean((Y.CP-pmin(pmax(NWcv_K4B_rcpp(X = X %*% b, Y = Y.CP, + # h = h), 0), 1))^2) + # return(cv) + # } + }else + { + stop("There's no weighted version of the K4_Biweight kernel.") + } + } + + if(method == "nlminb") + { + esti <- nlminb(start = c(initial[-1], 0), + objective = Eij3, + control = list(abs.tol = abs.tol)) + }else if (method == "optim") + { + # the new optimize function using optim, you can change the lower and upper + esti <- optim(par = c(initial[-1], 0), + fn = Eij3, + method = optim_method, + control = list(abstol = abs.tol)) + }else if (method == "nmk") + { + assertthat::assert_that(requireNamespace("dfoptim", quietly = TRUE)) + esti <- dfoptim::nmk(par = c(initial[-1], 0), + fn = Eij3, + control = list(tol = abs.tol)) + } + + results <- list(coef = c(1, esti$par[1:(number_p-1)]), + bandwidth = exp(esti$par[number_p]), + details = esti) + }else + { + if (kernel=="K2_Biweight") + { + if (is.null(wi.boot)) + { + Eij3 <- function(parameter){ + K <- function(x, h) 15/16*(1-(x/h)^2)^2 * (abs(x) <= h) * (x!=0) + + b <- c(1, parameter[1:(number_p-1)]) + h <- bandwidth + + x <- c(X%*%b) + y <- Y + + n <- length(y) + yo <- order(y) + ys <- y[yo] + uy <- rle(ys)[[1]] + cols <- cumsum(uy) + ei <- rep(0, n) + for (i in 1:n){ + Kih <- K(x-x[i],h=h) + denom <- sum(Kih) + ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + } + return(sum(ei)/n^2) + } + + # cv.b <- function(parameter) + # { + # b <- c(1, parameter[1:(number_p-1)]) + # cv <- mean((Y.CP-NWcv_K2B_rcpp(X = X %*% b, Y = Y.CP, + # h = bandwidth))^2) + # return(cv) + # } + }else + { + stop("There's no weighted version of the K2_Biweight kernel.") + } + }else if (kernel=="dnorm") + { + if (is.null(wi.boot)) + { + Eij3 <- function(parameter){ + K <- function(x, h) dnorm(x/h,0,1) * (x!=0) + + b <- c(1, parameter[1:(number_p-1)]) + h <- bandwidth + + x <- c(X%*%b) + y <- Y + + n <- length(y) + yo <- order(y) + ys <- y[yo] + uy <- rle(ys)[[1]] + cols <- cumsum(uy) + ei <- rep(0, n) + for (i in 1:n){ + Kih <- K(x-x[i],h=h) + denom <- sum(Kih) + ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + } + return(sum(ei)/n^2) + } + + # cv.b <- function(parameter) + # { + # b <- c(1, parameter[1:(number_p-1)]) + # cv <- mean((Y.CP-NWcv_dnorm_rcpp(X = X %*% b, Y = Y.CP, + # h = bandwidth))^2) + # return(cv) + # } + }else + { + # wi.boot <- as.vector(wi.boot) + # cv.b <- function(parameter) + # { + # b <- c(1, parameter[1:(number_p-1)]) + # cv <- mean((Y.CP-NWcv_K2B_w_rcpp(X = X %*% b, Y = Y.CP, + # h = bandwidth, w = wi.boot))^2) + # return(cv) + # } + stop("There's no weighted version of the dnorm kernel.") + } + }else if (kernel=="K4_Biweight") + { + if (is.null(wi.boot)) + { + Eij3 <- function(parameter){ + K <- function(x, h) 105/64*(1-3*((x/h)^2))*(1-(x/h)^2)^2 * (abs(x) <= h) * (x!=0) + + b <- c(1, parameter[1:(number_p-1)]) + h <- bandwidth + + x <- c(X%*%b) + y <- Y + + n <- length(y) + yo <- order(y) + ys <- y[yo] + uy <- rle(ys)[[1]] + cols <- cumsum(uy) + ei <- rep(0, n) + for (i in 1:n){ + Kih <- K(x-x[i],h=h) + denom <- sum(Kih) + ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + } + return(sum(ei)/n^2) + } + # cv.b <- function(parameter) + # { + # b <- c(1, parameter[1:(number_p-1)]) + # cv <- mean((Y.CP-pmin(pmax(NWcv_K4B_rcpp(X = X %*% b, Y = Y.CP, + # h = bandwidth), 0), 1))^2) + # return(cv) + # } + }else + { + stop("There's no weighted version of the K4_Biweight kernel.") + } + } + + if(method == "nlminb") + { + esti <- nlminb(start = initial[-1], + objective = Eij3, + control = list(abs.tol = abs.tol)) + }else if (method == "optim") + { + # the new optimize function using optim, you can change the lower and upper + esti <- optim(par = initial[-1], + fn = Eij3, + method = optim_method, + control = list(abstol = abs.tol)) + }else if (method == "nmk") + { + assertthat::assert_that(requireNamespace("dfoptim", quietly = TRUE)) + esti <- dfoptim::nmk(par = initial[-1], + fn = Eij3, + control = list(tol = abs.tol)) + }else rlang::abort("Invalid optimization method") + results <- list(coef = c(1, esti$par[1:(number_p-1)]), + bandwidth = bandwidth, + details = esti) + } + + return(results) +} + +SIDRnew_fixed_bandwidth <- +function(X, Y, ids, + Y.CP = NULL, + initial = NULL, + kernel = "K2_Biweight", + method = c("optim", "nlminb", "nmk"), + optim_method = "BFGS", + abs.tol = 1e-4 + ){ + X <- as.matrix(X) + Y <- as.matrix(Y) + + method <- match.arg(method) + + number_n <- dim(X)[1] + number_p <- dim(X)[2] + + if (is.null(initial)) + { + initial <- rep(0, number_p) + }else + { + initial <- as.vector(initial) + } + + if (kernel=="K2_Biweight") + K <- function(x) 15/16*(1-(x)^2)^2 * (abs(x) <= 1) + else if (kernel=="dnorm") + K <- function(x) dnorm(x,0,1) + else if (kernel=="K4_Biweight") + K <- function(x) 105/64*(1-3*((x)^2))*(1-(x)^2)^2 * (abs(x) <= 1) + else rlang::abort("Bad Kernel") + + Eij3 <- function(parameter){ + x <- c(X %*% parameter) + y <- Y + + n <- length(y) + yo <- order(y) + ys <- y[yo] + uy <- rle(ys)[[1]] + cols <- cumsum(uy) + ei <- rep(0, n) + for (i in 1:n){ + Kih <- ifelse(ids==ids[i], 0, K(x-x[i])) + denom <- sum(Kih) + ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + } + return(sum(ei)/n^2) + } + + + if(method == "nlminb") + { + esti <- nlminb(start = initial, + objective = Eij3, + control = list(abs.tol = abs.tol)) + }else if (method == "optim") + { + # the new optimize function using optim, you can change the lower and upper + esti <- optim(par = initial, + fn = Eij3, + method = optim_method, + control = list(abstol = abs.tol)) + }else if (method == "nmk") + { + assertthat::assert_that(requireNamespace("dfoptim", quietly = TRUE)) + esti <- dfoptim::nmk(par = initial, + fn = Eij3, + control = list(tol = abs.tol)) + }else rlang::abort("Invalid optimization method") + results <- list(coef = esti$par, + bandwidth = 1, + details = esti) +} diff --git a/SensIAT_sim_outcome_modeler_mave.R b/SensIAT_sim_outcome_modeler_mave.R new file mode 100644 index 0000000..ef3dad8 --- /dev/null +++ b/SensIAT_sim_outcome_modeler_mave.R @@ -0,0 +1,174 @@ +# Single Index Model using MAVE and Optimizing Bandwidth. +# +# Single index model estimation using minimum average variance estimation (MAVE). +# A direction is estimated using MAVE, and then the bandwidth is selected by +# minimization of the cross-validated pseudo-integrated squared error. +# Optionally, the initial coefficients of the outcome model can be re-estimated +# by optimization on a spherical manifold. This option requires the +# [ManifoldOptim][ManifoldOptim::manifold.optim] package. + +fit_SensIAT_single_index_norm1coef_model <- +function(X, Y, ids, + kernel = "K2_Biweight", + mave.method = "meanMAVE", + id = ..id.., + bw.selection = c('ise', 'mse'), + bw.method = c('optim', 'grid', 'optimize'), + bw.range = c(0.01, 1.5), + reestimate.coef = 0, + ... +){ + unique_y <- sort(unique(Y)) + + if (kernel=="K2_Biweight"){ + K <- function(x) 15/16*(1-(x)^2)^2 * (abs(x) <= 1) + K1<- function(x) 3/4*(1-(x)^2) * (abs(x) <= 1) + } else if (kernel=="dnorm") { + K <- function(x) dnorm(x,0,1) + K1 <- function(x) -x*dnorm(x, 0, 1) + } else if (kernel=="K4_Biweight"){ + K <- function(x) 105/64*(1-3*((x)^2)) * (1-(x)^2)^2 * (abs(x) <= 1) + } else{ + stop("Unknown kernel type. Please use either 'K2_Biweight', 'dnorm', or 'K4_Biweight'.") + } + + assertthat::assert_that(requireNamespace("MAVE", quietly = TRUE)) + mave_fit <- MAVE::mave.compute(X, Y, max.dim = 1, method = mave.method) + beta_hat <- mave_fit$dir[[1]][,1,drop=TRUE] + + + bw_old <- 1 + delta_beta = NA_real_ + res <- NULL + + Id_neq <- outer(ids, ids, \(x,y)as.numeric(x!=y)) + Imat <- outer(Y, unique_y, `<=`) + repeat{ + Xbeta <- (X %*% beta_hat)[,] + sigma <- sd(Xbeta) + D <- outer(Xbeta, Xbeta, FUN = `-`) + + bw.selection <- match.arg(bw.selection) + if(bw.selection == 'ise'){ + # Use Integrated Squared Error (ISE) for bandwidth selection + err <- function(log_bandwidth){ + W <- K(D/exp(log_bandwidth)) * Id_neq + denom <- rowSums(W) + Fhat <- sweep(W %*% Imat, 1, denom, "/") + Fhat[is.nan(Fhat)] <- 0 + return (sum((Imat - Fhat)^2)/length(Fhat)^2) + } + } else if(bw.selection == 'mse'){ + # Use Mean Squared Error (MSE) for bandwidth selection + err <- function(log_bandwidth){ + W <- K(D/exp(log_bandwidth)) * Id_neq + denom <- rowSums(W) + mean_Y <- c(W %*% Y) / denom + return (mean((Y - mean_Y)^2, na.rm = TRUE)) + } + } else { + stop("Unknown bw.selection type. Please use either 'ise' or 'mse'.") + } + + bw.method <- match.arg(bw.method) + if(bw.method == 'grid'){ + # Use grid search for bandwidth selection + log_bw_seq <- log(seq(min(bw.range), max(bw.range), length.out = 100) * sigma) + err_values <- purrr::map_dbl(log_bw_seq, err) + bw_opt <- log_bw_seq[which.min(err_values)] + bw.details <- list(par = bw_opt, + value = min(err_values), + convergence = 0, + message = "Grid search completed successfully.", + bw.method = bw.method, + log_bw_seq = log_bw_seq, + err_values = err_values + ) + } else if(bw.method == 'optim'){ + # Use optimization for bandwidth selection + initial <- log(sigma * 0.30) + bw.details <- optim(initial, err, method = "L-BFGS-B", lower = log(sigma * min(bw.range)), upper = log(sigma * max(bw.range)), ...) + bw.details$initial = initial + bw.details$bw.method = bw.method + bw_opt <- bw.details$par + } else if(bw.method == 'optimize'){ + # Use optimize for bandwidth selection + result <- stats::optimize(err, + interval = c(log(sigma * min(bw.range)), log(sigma * max(bw.range))), + ...) + bw_opt <- result$minimum + bw.details <- list(minimum = result$minimum, + value = result$objective, + bw.method = bw.method, + interval = c(log(sigma * 0.05), log(sigma * 1.5))) + } else{ + stop("Unknown bw.method type. Please use either 'optim' or 'grid'.") + } + + + delta_bw = abs(bw_opt - bw_old) + + if(!is.null(res) && res$fval <= bw.details$value) break + last.optimized = 'bandwidth' + + if(reestimate.coef <= 0) break + reestimate.coef <- reestimate.coef - 1 + + + bw_old <- bw_opt + rlang::check_installed("ManifoldOptim") + requireNamespace("ManifoldOptim", quietly = TRUE) + + beta_hat_old <- beta_hat + if(bw.selection == "ise") { + objFun <- function(beta){ + Xbeta <- (X %*% beta)[,] + D <- outer(Xbeta, Xbeta, FUN = `-`) + W <- K(D/exp(bw_opt)) * Id_neq + denom <- rowSums(W) + Fhat <- sweep(W %*% Imat, 1, denom, "/") + Fhat[is.nan(Fhat)] <- 0 + return (sum((Imat - Fhat)^2)/length(Fhat)^2) + } + } else if(bw.selection == "mse") { + objFun <- function(beta){ + Xbeta <- (X %*% beta)[,] + D <- outer(Xbeta, Xbeta, FUN = `-`) + W <- K1(D/exp(bw_opt)) * Id_neq + denom <- rowSums(W) + mean_Y <- c(W %*% Y) / denom + return (mean((Y - mean_Y)^2, na.rm = TRUE)) + } + } else { + stop("Unknown bw.selection type. Please use either 'ise' or 'mse'.") + } + mod <- Rcpp::Module("ManifoldOptim_module", PACKAGE = "ManifoldOptim") + prob <- new(mod$RProblem, objFun) + mani.params <- ManifoldOptim::get.manifold.params(IsCheckParams = FALSE) + solver.params <- ManifoldOptim::get.solver.params(IsCheckParams = FALSE) + maniDef <- ManifoldOptim::get.sphere.defn(length(beta_hat)) + + res <- ManifoldOptim::manifold.optim(prob, maniDef, + mani.params = mani.params, + solver.params = solver.params, x0 = beta_hat) + if(res$fval >= bw.details$value) break + last.optimized = 'beta' + + beta_hat <- as.vector(res$xopt) + delta_beta = abs(beta_hat - beta_hat_old) + } + return( + list( + coef = beta_hat, + bandwidth = exp(bw_opt), + details = list(bw.last = bw.details, + bw.delta = delta_bw, + bw.unscaled = bw_opt/sigma, + beta.delta = delta_beta, + beta.last = res, + last.optimized = last.optimized + ) + ) + ) +} + From 0384af82674ef1d056d8cdb31d973d3e537c09d5 Mon Sep 17 00:00:00 2001 From: Jennifer Xu <48976310+jenniferxu1101@users.noreply.github.com> Date: Wed, 11 Feb 2026 11:01:02 -0500 Subject: [PATCH 03/19] Add files via upload --- singleindexmodelfunctions.R | 399 ++++++++++++++++++++++++++++++++++++ 1 file changed, 399 insertions(+) create mode 100644 singleindexmodelfunctions.R diff --git a/singleindexmodelfunctions.R b/singleindexmodelfunctions.R new file mode 100644 index 0000000..176314c --- /dev/null +++ b/singleindexmodelfunctions.R @@ -0,0 +1,399 @@ +# single index model's function + +# return the initial value for beta +cumuSIR <- function(X, Y, eps = 1e-7){ + X <- as.matrix(X) + Y <- as.matrix(Y) + + number_n <- dim(X)[1] + number_p <- dim(X)[2] + + Y.CP <- matrix( + Y[rep(1:number_n, times = number_n), ]<= + Y[rep(1:number_n, each = number_n), ], + nrow = number_n, ncol = number_n + ) + + # centralizing covariates + X.cs <- t(t(X)-colMeans(X)) + + # calculating m(y)=\E[X_i 1(Y_i\leq y)] + m.y <- t(X.cs) %*% Y.CP/number_n + # calculating K=\E[m(Y_i)m(Y_i)^T] + Km <- m.y %*% t(m.y)/number_n + + Bhat <- eigen(solve(var(X) + eps*diag(number_p), Km))$vectors[, 1] + + return(Bhat) +} + +# find the optimal bandwidth, default set for kernel is +SIDRnew_h2 <- function(X, Y, + Y.CP = NULL, + beta = NULL, + kernel = "K2_Biweight", + method = "optim", + optim_method = "BFGS", + abs.tol = 1e-4, + initial_bandwidth = NULL, + wi.boot = NULL, penalty=0) +{ + X <- as.matrix(X) + Y <- as.matrix(Y) + + number_n <- dim(X)[1] + number_p <- dim(X)[2] + + if (is.null(beta)) + { + beta <- c(1, rep(0, number_p-1)) + }else + { + beta <- as.vector(beta) + beta <- beta/beta[1] + } + + if (is.null(initial_bandwidth)) + { + if (kernel=="K2_Biweight") + { + if (is.null(wi.boot)) + { + Eij3 <- function(parameter){ + K <- function(x, h) 15/16*(1-(x/h)^2)^2 * (abs(x) <= h) + b <- beta + h <- exp(parameter) + xb <- c(X%*%b) + y <- Y + n <- length(y) + yo <- order(y) + ys <- y[yo] + uy <- rle(ys)[[1]] + cols <- cumsum(uy) + ei <- rep(0, n) + for (i in 1:n){ + Kih <- K(xb-xb[i],h=h) + Kih[i] <- 0 + denom <- sum(Kih) + ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + } + return(sum(ei)/n^2+penalty*h) + } + + # cv.bh <- function(parameter) + # { + # b <- c(1, parameter[1:(number_p-1)]) + # h <- exp(parameter[number_p]) + # cv <- mean((Y.CP-NWcv_K2B_rcpp(X = X %*% b, Y = Y.CP, + # h = h))^2) + # return(cv) + # } + }else + { + stop("There's no weighted version of the K2_Biweight kernel.") + } + }else if (kernel == "dnorm") + { + if (is.null(wi.boot)) + { + Eij3 <- function(parameter){ + K <- function(x, h) dnorm(x/h, 0, 1) + + b <- beta + h <- exp(parameter) + + x <- c(X%*%b) + y <- Y + + n <- length(y) + yo <- order(y) + ys <- y[yo] + uy <- rle(ys)[[1]] + cols <- cumsum(uy) + ei <- rep(0, n) + for (i in 1:n){ + Kih <- K(x-x[i],h=h) + Kih[i] <- 0 + denom <- sum(Kih) + ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + } + return(sum(ei)/n^2) + } + # cv.bh <- function(parameter) + # { + # b <- c(1, parameter[1:(number_p-1)]) + # h <- exp(parameter[number_p]) + # cv <- mean((Y.CP-NWcv_dnorm_rcpp(X = X %*% b, Y = Y.CP, + # h = h))^2) + # return(cv) + # } + }else + { + # wi.boot <- as.vector(wi.boot) + # cv.bh <- function(parameter) + # { + # b <- c(1, parameter[1:(number_p-1)]) + # h <- exp(parameter[number_p]) + # cv <- mean((Y.CP-NWcv_K2B_w_rcpp(X = X %*% b, Y = Y.CP, + # h = h, w = wi.boot))^2) + # return(cv) + # } + stop("There's no weighted version of the dnorm kernel.") + } + }else if (kernel=="K4_Biweight") + { + if (is.null(wi.boot)) + { + Eij3 <- function(parameter){ + K <- function(x, h) 105/64*(1-3*((x/h)^2))*(1-(x/h)^2)^2 * (abs(x) <= h) + + b <- beta + h <- exp(parameter) + + x <- c(X%*%b) + y <- Y + + n <- length(y) + yo <- order(y) + ys <- y[yo] + uy <- rle(ys)[[1]] + cols <- cumsum(uy) + ei <- rep(0, n) + for (i in 1:n){ + Kih <- K(x-x[i],h=h) + Kih[i] <- 0 + denom <- sum(Kih) + ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + } + return(sum(ei)/n^2) + } + # cv.bh <- function(parameter) + # { + # b <- c(1, parameter[1:(number_p-1)]) + # h <- exp(parameter[number_p]) + # cv <- mean((Y.CP-pmin(pmax(NWcv_K4B_rcpp(X = X %*% b, Y = Y.CP, + # h = h), 0), 1))^2) + # return(cv) + # } + }else + { + stop("There's no weighted version of the K4_Biweight kernel.") + } + } + + if(method == "nlminb") + { + esti <- nlminb(start = 0, + objective = Eij3, + control = list(abs.tol = abs.tol)) + results <- list(coef = beta, + bandwidth = exp(esti$par), + details = esti) + }else if (method == "optim") + { + # the new optimize function using optim, you can change the lower and upper + esti <- optim(par = 0, + fn = Eij3, + method = optim_method, + control = list(abstol = abs.tol)) + results <- list(coef = beta, + bandwidth = exp(esti$par), + details = esti) + } + else if (method == "optimize") + { + # # esti <- nmk(par = 0, + # # fn = Eij3, + # # control = list(tol = abs.tol)) + esti <- optimize(Eij3, + interval = c(-10, 10)) + results <- list(coef = beta, + bandwidth = exp(esti$minimum), + details = esti) + } + + }else + { + if (kernel=="K2_Biweight") + { + if (is.null(wi.boot)) + { + Eij3 <- function(parameter){ + K <- function(x, h) 15/16*(1-(x/h)^2)^2 * (abs(x) <= h) + b <- beta + h <- exp(parameter) + xb <- c(X%*%b) + y <- Y + n <- length(y) + yo <- order(y) + ys <- y[yo] + uy <- rle(ys)[[1]] + cols <- cumsum(uy) + ei <- rep(0, n) + for (i in 1:n){ + Kih <- K(xb-xb[i],h=h) + Kih[i] <- 0 + denom <- sum(Kih) + ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + } + return(sum(ei)/n^2+penalty*h) + } + + # cv.b <- function(parameter) + # { + # b <- c(1, parameter[1:(number_p-1)]) + # cv <- mean((Y.CP-NWcv_K2B_rcpp(X = X %*% b, Y = Y.CP, + # h = bandwidth))^2) + # return(cv) + # } + }else + { + stop("There's no weighted version of the K2_Biweight kernel.") + } + }else if (kernel=="dnorm") + { + if (is.null(wi.boot)) + { + Eij3 <- function(parameter){ + K <- function(x, h) dnorm(x/h,0,1) + + b <- beta + h <- exp(parameter) + + x <- c(X%*%b) + y <- Y + + n <- length(y) + yo <- order(y) + ys <- y[yo] + uy <- rle(ys)[[1]] + cols <- cumsum(uy) + ei <- rep(0, n) + for (i in 1:n){ + Kih <- K(x-x[i],h=h) + Kih[i] <- 0 + denom <- sum(Kih) + ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + } + return(sum(ei)/n^2) + } + + # cv.b <- function(parameter) + # { + # b <- c(1, parameter[1:(number_p-1)]) + # cv <- mean((Y.CP-NWcv_dnorm_rcpp(X = X %*% b, Y = Y.CP, + # h = bandwidth))^2) + # return(cv) + # } + }else + { + # wi.boot <- as.vector(wi.boot) + # cv.b <- function(parameter) + # { + # b <- c(1, parameter[1:(number_p-1)]) + # cv <- mean((Y.CP-NWcv_K2B_w_rcpp(X = X %*% b, Y = Y.CP, + # h = bandwidth, w = wi.boot))^2) + # return(cv) + # } + stop("There's no weighted version of the dnorm kernel.") + } + }else if (kernel=="K4_Biweight") + { + if (is.null(wi.boot)) + { + Eij3 <- function(parameter){ + K <- function(x, h) 105/64*(1-3*((x/h)^2))*(1-(x/h)^2)^2 * (abs(x) <= h) + + b <- beta + h <- exp(parameter) + + x <- c(X%*%b) + y <- Y + + n <- length(y) + yo <- order(y) + ys <- y[yo] + uy <- rle(ys)[[1]] + cols <- cumsum(uy) + ei <- rep(0, n) + for (i in 1:n){ + Kih <- K(x-x[i],h=h) + Kih[i] <- 0 + denom <- sum(Kih) + ei[i] <- sum(uy*(1*(y[i] <= ys)[cols] - (denom != 0)* cumsum(Kih[yo])[cols] / (denom + (denom == 0)))^2) + } + return(sum(ei)/n^2) + } + # cv.b <- function(parameter) + # { + # b <- c(1, parameter[1:(number_p-1)]) + # cv <- mean((Y.CP-pmin(pmax(NWcv_K4B_rcpp(X = X %*% b, Y = Y.CP, + # h = bandwidth), 0), 1))^2) + # return(cv) + # } + }else + { + stop("There's no weighted version of the K4_Biweight kernel.") + } + } + + if(method == "nlminb") + { + esti <- nlminb(start = initial_bandwidth, + objective = Eij3, + control = list(abs.tol = abs.tol)) + results <- list(coef = beta, + bandwidth = exp(esti$par), + details = esti) + }else if (method == "optim") + { + # the new optimize function using optim, you can change the lower and upper + esti <- optim(par = initial_bandwidth, + fn = Eij3, + method = optim_method, + control = list(abstol = abs.tol)) + results <- list(coef = beta, + bandwidth = exp(esti$par), + details = esti) + } + else if (method == "optimize") + { + # # esti <- nmk(par = initial_bandwidth, + # # fn = Eij3, + # # control = list(tol = abs.tol)) + esti <- optimize(Eij3, + interval = c(-10, 10)) + results <- list(coef = beta, + bandwidth = exp(esti$minimum), + details = esti) + + } + } + + return(results) +} + + +# return the prediction using "K2_Biweight" kernel +NW_new <- function(Xb, Y, xb, y, h, kernel = "K2_Biweight"){ + + if(kernel == "dnorm"){ + K <- function(x, h){dnorm(x/h, 0, 1)} # Gaussian + } else if(kernel == "K2_Biweight"){ + K <- function(x, h){15/16*(1-(x/h)^2)^2 * (abs(x) <= h)} # K2_Biweight + } else if(kernel=="K4_Biweight"){ + K <- function(x, h){105/64*(1-3*((x/h)^2))*(1-(x/h)^2)^2 * (abs(x) <= h) }# K4_Biweight + } + + Kxb <- sapply(xb, function(x, Xb) K(Xb-x, h), Xb=Xb) + + Ylty <- sapply(y, function(x, Y) 1*(Y <= x), Y=Y) + + denom <- colSums(Kxb) + + fyxb <- (denom!=0)*crossprod(Kxb, Ylty)/(denom + (denom==0)) + + return(fyxb) + +} + From 2ac1dd87960e8240992f7941d070e5540b38c51b Mon Sep 17 00:00:00 2001 From: jenniferxu1101 Date: Wed, 11 Feb 2026 11:14:35 -0500 Subject: [PATCH 04/19] change directory --- simulation_est_slurm_double_pboot.R | 36 +++++------------------------ 1 file changed, 6 insertions(+), 30 deletions(-) diff --git a/simulation_est_slurm_double_pboot.R b/simulation_est_slurm_double_pboot.R index 2db2c7a..4f9f832 100644 --- a/simulation_est_slurm_double_pboot.R +++ b/simulation_est_slurm_double_pboot.R @@ -1,18 +1,3 @@ -indx <- as.numeric(commandArgs(trailingOnly=TRUE)) - -scenarios <- expand.grid(sim.size_list=500, index=1:50, imputed=1, - kernel="dnorm", method="optim", - single_index_method="norm1coef") -scenarios <- scenarios[order(scenarios$sim.size_list),] - - -sim.size <- scenarios$sim.size_list[indx] -index <- scenarios$index[indx] -imputed <- scenarios$imputed[indx] -kernel <- as.character(scenarios$kernel[indx]) -method <- as.character(scenarios$method[indx]) -single_index_method <- as.character(scenarios$single_index_method[indx]) - ################### ## load packages ## ################### @@ -23,21 +8,12 @@ library(parallel) library(foreach) library(doParallel) -source("/uufs/chpc.utah.edu/common/home/u6070035/CCS/code/HelperFunction.R") -source("/uufs/chpc.utah.edu/common/home/u6070035/CCS/code/singleindexmodelfunctions.R") -source("/uufs/chpc.utah.edu/common/home/u6070035/CCS/code/est_s_t_y.R") -source("/uufs/chpc.utah.edu/common/home/u6070035/CCS/code/SIDR_Ravinew.R") -source("/uufs/chpc.utah.edu/common/home/u6070035/CCS/code/SIDRnew.R") -source("/uufs/chpc.utah.edu/common/home/u6070035/CCS/code/SensIAT_sim_outcome_modeler_mave.R") - -ntasks <- Sys.getenv("SLURM_NTASKS") -if (ntasks == '') { - ntasks <- 4 -} else { - ntasks <- strtoi(ntasks) } -cat("This script use ", ntasks, " cores\n") - -package_list <- c("splines", "tidyverse", "gridExtra") +source(paste0(getwd(), "/HelperFunction.R")) +source(paste0(getwd(), "/singleindexmodelfunctions.R")) +source(paste0(getwd(), "/est_s_t_y.R")) +source(paste0(getwd(), "/SIDR_Ravinew.R")) +source(paste0(getwd(), "/SIDRnew.R")) +source(paste0(getwd(), "/SensIAT_sim_outcome_modeler_mave.R")) nboot_B <- 500 nboot_C <- 100 From ed14902badb6ea20932bc79116754c228cdf6f48 Mon Sep 17 00:00:00 2001 From: jenniferxu1101 Date: Wed, 11 Feb 2026 12:01:22 -0500 Subject: [PATCH 05/19] For loop instead of foreach --- simulation_est_slurm_double_pboot.R | 81 ++++++++--------------------- 1 file changed, 23 insertions(+), 58 deletions(-) diff --git a/simulation_est_slurm_double_pboot.R b/simulation_est_slurm_double_pboot.R index 4f9f832..797ddd3 100644 --- a/simulation_est_slurm_double_pboot.R +++ b/simulation_est_slurm_double_pboot.R @@ -4,9 +4,6 @@ library(tidyverse) library(splines) library(gridExtra) -library(parallel) -library(foreach) -library(doParallel) source(paste0(getwd(), "/HelperFunction.R")) source(paste0(getwd(), "/singleindexmodelfunctions.R")) @@ -15,8 +12,8 @@ source(paste0(getwd(), "/SIDR_Ravinew.R")) source(paste0(getwd(), "/SIDRnew.R")) source(paste0(getwd(), "/SensIAT_sim_outcome_modeler_mave.R")) -nboot_B <- 500 -nboot_C <- 100 +nboot_B <- 10 +nboot_C <- 5 grid <- expand.grid(b = 1:nboot_B, c = 1:nboot_C) gamma_length <- length(seq(-5, 5, by=1)) @@ -25,45 +22,34 @@ gamma_length <- length(seq(-5, 5, by=1)) ## simulation ## ################ -## parametric bootstrap +## empty containers +Boot.Est <- matrix(NA, nrow=nrow(grid), ncol=264) -cl <- makeCluster(ntasks) -registerDoParallel(cl) +## parametric bootstrap +load("/Users/shiyaoxu/Documents/2025_2026/research/comprehensive_cohort/simData/imputed1/n500/sim1.rds") +X_sim = data.frame(dplyr::select(data, c(age, pain_bq, expectationb, ChronicPainb))) -for (j in (40*(index-1)+1):(40*index)){ - - set.seed(1000090+j) - seed_B <- sample(1:10000000, size=nboot_B) - - load(paste0("/uufs/chpc.utah.edu/common/home/u6070035/CCS/simData/imputed", imputed, "/n", sim.size, "/sim", j, ".rds")) - X_sim = data.frame(dplyr::select(data, c(age, pain_bq, expectationb, ChronicPainb))) +for(i in 1:nboot_B){ -Boot.Est <- foreach(i=1:nrow(grid), .combine=rbind, .packages=package_list, .errorhandling = "pass") %dopar% { - - ## simulate outer bootstrap sample - data.new <- pboot(data=data, X_sim=X_sim, sim.size=sim.size, seed=seed_B[grid$b[i]]) + data.new <- pboot(data=data, X_sim=X_sim, sim.size=500, seed=NULL) X.new <- data.frame(dplyr::select(data.new, c(age, pain_bq, expectationb, ChronicPainb))) - if(grid$c[i]==1){ - ## outer estimator - topical_vals <- fit_one(data=data.new, X=X.new, trt_val=1) - oral_vals <- fit_one(data=data.new, X=X.new, trt_val=0) - }else{ - topical_vals <- rep(NA_real_, 6*gamma_length) - oral_vals <- rep(NA_real_, 6*gamma_length) - } - - ## simulate inner bootstrap sample - data.new2 <- pboot(data=data.new, X_sim=X.new, sim.size=sim.size, seed=NULL) - X.new2 <- data.frame(dplyr::select(data.new2, c(age, pain_bq, expectationb, ChronicPainb))) - ## inner estimator - temp_topical <- fit_one(data=data.new2, X=X.new2, trt_val=1) - temp_oral <- fit_one(data=data.new2, X=X.new2, trt_val=0) + topical_vals <- fit_one(data=data.new, X=X.new, trt_val=1) + oral_vals <- fit_one(data=data.new, X=X.new, trt_val=0) + Boot.Est[which(grid$b==i&grid$c==1), ] <- c(topical_vals, oral_vals) - ## return result - c(topical_vals, oral_vals, temp_topical, temp_oral) + for(j in 1:nboot_C){ + + data.new2 <- pboot(data=data.new, X_sim=X.new, sim.size=sim.size, seed=NULL) + X.new2 <- data.frame(dplyr::select(data.new2, c(age, pain_bq, expectationb, ChronicPainb))) + ## inner estimator + temp_topical <- fit_one(data=data.new2, X=X.new2, trt_val=1) + temp_oral <- fit_one(data=data.new2, X=X.new2, trt_val=0) + Boot.Est[which(grid$b==i&grid$c==j), 133:264] <- c(temp_topical, temp_oral) + + } } - + load(paste0("/uufs/chpc.utah.edu/common/home/u6070035/CCS/simResult/imputed", imputed, "/n", sim.size, "/", kernel, "_", single_index_method, "_fold5/sim", j, ".RData")) Q_b_topical <- sapply(1:nboot_B, function(x){ loc <- which(grid$b==x) @@ -302,27 +288,6 @@ var_R0_oral_boot <- sapply(1:nboot_B, function(x){ Boot.Est[loc_b, 122:132] }) -## bootstrap result -save(est_topical_boot, est_R1_topical_boot, est_R0_topical_boot, - est_oral_boot, est_R1_oral_boot, est_R0_oral_boot, - var_topical_boot, var_R1_topical_boot, var_R0_topical_boot, - var_oral_boot, var_R1_oral_boot, var_R0_oral_boot, - Q_b_topical, Q_b_R1_topical, Q_b_R0_topical, - Q_b_oral, Q_b_R1_oral, Q_b_R0_oral, - t_b_topical, t_b_R1_topical, t_b_R0_topical, - t_b_oral, t_b_R1_oral, t_b_R0_oral, - Q_t_b_topical, Q_t_b_R1_topical, Q_t_b_R0_topical, - Q_t_b_oral, Q_t_b_R1_oral, Q_t_b_R0_oral, - abs_Q_t_b_topical, abs_Q_t_b_R1_topical, abs_Q_t_b_R0_topical, - abs_Q_t_b_oral, abs_Q_t_b_R1_oral, abs_Q_t_b_R0_oral, - t_b_sd_topical, t_b_R1_sd_topical, t_b_R0_sd_topical, - t_b_sd_oral, t_b_R1_sd_oral, t_b_R0_sd_oral, - file=paste0("/uufs/chpc.utah.edu/common/home/u6070035/CCS/simResult/imputed", imputed, "/double_pboot_n", sim.size, "/", kernel, "_", single_index_method, "/sim", j, ".RData")) - print(j) - -} - -stopCluster(cl) From 3cea46ba78ec171de78e5120596b29d2914d9ecb Mon Sep 17 00:00:00 2001 From: jenniferxu1101 Date: Wed, 11 Feb 2026 12:09:24 -0500 Subject: [PATCH 06/19] print --- simulation_est_slurm_double_pboot.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/simulation_est_slurm_double_pboot.R b/simulation_est_slurm_double_pboot.R index 797ddd3..86bffa3 100644 --- a/simulation_est_slurm_double_pboot.R +++ b/simulation_est_slurm_double_pboot.R @@ -48,6 +48,8 @@ for(i in 1:nboot_B){ } + print(i) + } load(paste0("/uufs/chpc.utah.edu/common/home/u6070035/CCS/simResult/imputed", imputed, "/n", sim.size, "/", kernel, "_", single_index_method, "_fold5/sim", j, ".RData")) From 4b8675a28bf0e078bef0b004df2fef9a737dfe6b Mon Sep 17 00:00:00 2001 From: Jennifer Xu <48976310+jenniferxu1101@users.noreply.github.com> Date: Wed, 11 Feb 2026 12:13:39 -0500 Subject: [PATCH 07/19] Add files via upload --- SensIAT_sim_outcome_modeler_mave.R | 12 ++++++++---- singleindexmodelfunctions.R | 1 + 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/SensIAT_sim_outcome_modeler_mave.R b/SensIAT_sim_outcome_modeler_mave.R index ef3dad8..187e028 100644 --- a/SensIAT_sim_outcome_modeler_mave.R +++ b/SensIAT_sim_outcome_modeler_mave.R @@ -11,6 +11,7 @@ fit_SensIAT_single_index_norm1coef_model <- function(X, Y, ids, kernel = "K2_Biweight", mave.method = "meanMAVE", + use_mave, id = ..id.., bw.selection = c('ise', 'mse'), bw.method = c('optim', 'grid', 'optimize'), @@ -32,10 +33,13 @@ function(X, Y, ids, stop("Unknown kernel type. Please use either 'K2_Biweight', 'dnorm', or 'K4_Biweight'.") } - assertthat::assert_that(requireNamespace("MAVE", quietly = TRUE)) - mave_fit <- MAVE::mave.compute(X, Y, max.dim = 1, method = mave.method) - beta_hat <- mave_fit$dir[[1]][,1,drop=TRUE] - + if(use_mave){ + assertthat::assert_that(requireNamespace("MAVE", quietly = TRUE)) + mave_fit <- MAVE::mave.compute(X, Y, max.dim = 1, method = mave.method) + beta_hat <- mave_fit$dir[[1]][,1,drop=TRUE] + }else{ + beta_hat <- cumuSIR(X = X, Y = Y) + } bw_old <- 1 delta_beta = NA_real_ diff --git a/singleindexmodelfunctions.R b/singleindexmodelfunctions.R index 176314c..cd4555c 100644 --- a/singleindexmodelfunctions.R +++ b/singleindexmodelfunctions.R @@ -27,6 +27,7 @@ cumuSIR <- function(X, Y, eps = 1e-7){ return(Bhat) } + # find the optimal bandwidth, default set for kernel is SIDRnew_h2 <- function(X, Y, Y.CP = NULL, From e1c66e646d4d5c2987e008a51a31fd53e1700a73 Mon Sep 17 00:00:00 2001 From: jenniferxu1101 Date: Wed, 11 Feb 2026 13:14:40 -0500 Subject: [PATCH 08/19] profvis --- profvis.html | 4637 +++++++++++++++++++++++++++ simulation_est_slurm_double_pboot.R | 11 +- 2 files changed, 4646 insertions(+), 2 deletions(-) create mode 100644 profvis.html diff --git a/profvis.html b/profvis.html new file mode 100644 index 0000000..cf77c96 --- /dev/null +++ b/profvis.html @@ -0,0 +1,4637 @@ + + + + +profvis + + + + + + + + + + + + + +
+
+
+ + + + diff --git a/simulation_est_slurm_double_pboot.R b/simulation_est_slurm_double_pboot.R index 86bffa3..4dc220f 100644 --- a/simulation_est_slurm_double_pboot.R +++ b/simulation_est_slurm_double_pboot.R @@ -4,6 +4,7 @@ library(tidyverse) library(splines) library(gridExtra) +library(profvis) source(paste0(getwd(), "/HelperFunction.R")) source(paste0(getwd(), "/singleindexmodelfunctions.R")) @@ -12,7 +13,7 @@ source(paste0(getwd(), "/SIDR_Ravinew.R")) source(paste0(getwd(), "/SIDRnew.R")) source(paste0(getwd(), "/SensIAT_sim_outcome_modeler_mave.R")) -nboot_B <- 10 +nboot_B <- 2 nboot_C <- 5 grid <- expand.grid(b = 1:nboot_B, c = 1:nboot_C) @@ -22,6 +23,8 @@ gamma_length <- length(seq(-5, 5, by=1)) ## simulation ## ################ +pv <- profvis({ + ## empty containers Boot.Est <- matrix(NA, nrow=nrow(grid), ncol=264) @@ -39,7 +42,7 @@ for(i in 1:nboot_B){ for(j in 1:nboot_C){ - data.new2 <- pboot(data=data.new, X_sim=X.new, sim.size=sim.size, seed=NULL) + data.new2 <- pboot(data=data.new, X_sim=X.new, sim.size=500, seed=NULL) X.new2 <- data.frame(dplyr::select(data.new2, c(age, pain_bq, expectationb, ChronicPainb))) ## inner estimator temp_topical <- fit_one(data=data.new2, X=X.new2, trt_val=1) @@ -51,6 +54,10 @@ for(i in 1:nboot_B){ print(i) } + +}) + +htmlwidgets::saveWidget(pv, "profvis.html") load(paste0("/uufs/chpc.utah.edu/common/home/u6070035/CCS/simResult/imputed", imputed, "/n", sim.size, "/", kernel, "_", single_index_method, "_fold5/sim", j, ".RData")) Q_b_topical <- sapply(1:nboot_B, function(x){ From 1f07b516409f4050f3f71013f989c17e49d44ff9 Mon Sep 17 00:00:00 2001 From: jenniferxu1101 Date: Thu, 12 Feb 2026 14:08:09 -0500 Subject: [PATCH 09/19] add initial coef for gam --- HelperFunction.R | 247 +++++++++++++++++++++++++++- est_s_t_y.R | 28 ++-- simulation_est_slurm_double_pboot.R | 27 ++- 3 files changed, 279 insertions(+), 23 deletions(-) diff --git a/HelperFunction.R b/HelperFunction.R index 73d4448..dae5a10 100644 --- a/HelperFunction.R +++ b/HelperFunction.R @@ -221,6 +221,8 @@ pboot <- function(data, X_sim, sim.size, seed){ X_R0 <- X_sim[which(data$R==0), ] t_R0 <- data$t[which(data$R==0)] + X_R1 <- X_sim[which(data$R==1), ] + t_R1 <- data$t[which(data$R==1)] X_with_t_R0 <- X_with_t[which(data$R==0), ] X_with_t_R1 <- X_with_t[which(data$R==1), ] @@ -229,6 +231,7 @@ pboot <- function(data, X_sim, sim.size, seed){ g.fit <- mgcv::gam(as.formula(paste("data$R ~", gam.var)), data=X_sim, family=binomial) t_R0.fit <- mgcv::gam(as.formula(paste("t_R0 ~", gam.var)), data=X_R0, family=binomial) ## treatment model + t_R1.fit <- mgcv::gam(as.formula(paste("t_R1 ~", gam.var)), data=X_R1, family=binomial) ## treatment model M_R0.fit <- mgcv::gam(as.formula(paste("M_R0 ~", gam.var.M)), data=X_with_t_R0, family=binomial) ## missing data model M_R1.fit <- mgcv::gam(as.formula(paste("M_R1 ~", gam.var.M)), data=X_with_t_R1, family=binomial) ## missing data model @@ -379,29 +382,259 @@ pboot <- function(data, X_sim, sim.size, seed){ data.new <- data.frame(cbind(Y.new, M.new, R.new, T.new, X.new)) colnames(data.new) <- c("Y", "M", "R", "t", colnames(X.new)) + coef_g.fit <- coef(g.fit) + coef_t_R0.fit <- coef(t_R0.fit) + coef_t_R1.fit <- coef(t_R1.fit) + coef_M_R0.fit <- coef(M_R0.fit) + coef_M_R1.fit <- coef(M_R1.fit) + + return(list(data.new=data.new, coef_g.fit=coef_g.fit, coef_t_R0.fit=coef_t_R0.fit, + coef_t_R1.fit=coef_t_R1.fit, coef_M_R0.fit=coef_M_R0.fit, coef_M_R1.fit=coef_M_R1.fit)) + +} + +## 1st step of parametric bootstrap: model fitting using original data +pboot_model <- function(data, X_sim){ + + ## fit models + X_with_t <- cbind(as.factor(data$t), X_sim) + colnames(X_with_t)[1] <- "treatment" + gam.var <- paste(gam.variables(X_sim), collapse = "+") ## gam variables for treatment assignment model + gam.var.M <- paste(gam.variables(X_with_t), collapse = "+") ## gam variables for missing data model + index.var.Y <- paste(single.index.variables(X_sim), collapse = "+") + X_adjust <- model.matrix(as.formula(paste("~", index.var.Y)), data = X_sim)[,-1] + X_adjust_scale <- scale(X_adjust) + + X_adjust_scale_t_R0 <- X_adjust_scale[which(data$t==1 & data$R==0), ] + X_adjust_scale_t_R1 <- X_adjust_scale[which(data$t==1 & data$R==1), ] + X_adjust_scale_t0_R0 <- X_adjust_scale[which(data$t==0 & data$R==0), ] + X_adjust_scale_t0_R1 <- X_adjust_scale[which(data$t==0 & data$R==1), ] + M_t_R0 <- data$M[which(data$t==1 & data$R==0)] + M_t_R1 <- data$M[which(data$t==1 & data$R==1)] + M_t0_R0 <- data$M[which(data$t==0 & data$R==0)] + M_t0_R1 <- data$M[which(data$t==0 & data$R==1)] + Y_t_R0 <- data$Y[which(data$t==1 & data$R==0)] + Y_t_R1 <- data$Y[which(data$t==1 & data$R==1)] + Y_t0_R0 <- data$Y[which(data$t==0 & data$R==0)] + Y_t0_R1 <- data$Y[which(data$t==0 & data$R==1)] + + X_R0 <- X_sim[which(data$R==0), ] + t_R0 <- data$t[which(data$R==0)] + X_R1 <- X_sim[which(data$R==1), ] + t_R1 <- data$t[which(data$R==1)] + + X_with_t_R0 <- X_with_t[which(data$R==0), ] + X_with_t_R1 <- X_with_t[which(data$R==1), ] + M_R0 <- data$M[which(data$R==0)] + M_R1 <- data$M[which(data$R==1)] + + g.fit <- mgcv::gam(as.formula(paste("data$R ~", gam.var)), data=X_sim, family=binomial) + t_R0.fit <- mgcv::gam(as.formula(paste("t_R0 ~", gam.var)), data=X_R0, family=binomial) ## treatment model + t_R1.fit <- mgcv::gam(as.formula(paste("t_R1 ~", gam.var)), data=X_R1, family=binomial) ## treatment model + M_R0.fit <- mgcv::gam(as.formula(paste("M_R0 ~", gam.var.M)), data=X_with_t_R0, family=binomial) ## missing data model + M_R1.fit <- mgcv::gam(as.formula(paste("M_R1 ~", gam.var.M)), data=X_with_t_R1, family=binomial) ## missing data model + + fit_t_R1_h <- try(fit_SensIAT_single_index_norm1coef_model(X = X_adjust_scale_t_R1[which(M_t_R1==1), ], + Y = Y_t_R1[which(M_t_R1==1)], + ids=1:length(Y_t_R1[which(M_t_R1==1)]), + kernel="dnorm", bw.selection="ise", bw.method="optim", + use_mave=TRUE), silent = TRUE) + if (inherits(fit_t_R1_h, "try-error")) { + fit_t_R1_h <- try(fit_SensIAT_single_index_norm1coef_model(X = X_adjust_scale_t_R1[which(M_t_R1==1), ], + Y = Y_t_R1[which(M_t_R1==1)], + ids=1:length(Y_t_R1[which(M_t_R1==1)]), + kernel="dnorm", bw.selection="ise", bw.method="optim", + use_mave=FALSE), silent = TRUE) + } + + fit_t_R0_h <- try(fit_SensIAT_single_index_norm1coef_model(X = X_adjust_scale_t_R0[which(M_t_R0==1), ], + Y = Y_t_R0[which(M_t_R0==1)], + ids=1:length(Y_t_R0[which(M_t_R0==1)]), + kernel="dnorm", bw.selection="ise", bw.method="optim", + use_mave=TRUE), silent = TRUE) + if (inherits(fit_t_R0_h, "try-error")) { + fit_t_R0_h <- try(fit_SensIAT_single_index_norm1coef_model(X = X_adjust_scale_t_R0[which(M_t_R0==1), ], + Y = Y_t_R0[which(M_t_R0==1)], + ids=1:length(Y_t_R0[which(M_t_R0==1)]), + kernel="dnorm", bw.selection="ise", bw.method="optim", + use_mave=FALSE), silent = TRUE) + } + + fit_t0_R1_h <- try(fit_SensIAT_single_index_norm1coef_model(X = X_adjust_scale_t0_R1[which(M_t0_R1==1), ], + Y = Y_t0_R1[which(M_t0_R1==1)], + ids=1:length(Y_t0_R1[which(M_t0_R1==1)]), + kernel="dnorm", bw.selection="ise", bw.method="optim", + use_mave=TRUE), silent = TRUE) + if (inherits(fit_t0_R1_h, "try-error")) { + fit_t0_R1_h <- try(fit_SensIAT_single_index_norm1coef_model(X = X_adjust_scale_t0_R1[which(M_t0_R1==1), ], + Y = Y_t0_R1[which(M_t0_R1==1)], + ids=1:length(Y_t0_R1[which(M_t0_R1==1)]), + kernel="dnorm", bw.selection="ise", bw.method="optim", + use_mave=FALSE), silent = TRUE) + } + + fit_t0_R0_h <- try(fit_SensIAT_single_index_norm1coef_model(X = X_adjust_scale_t0_R0[which(M_t0_R0==1), ], + Y = Y_t0_R0[which(M_t0_R0==1)], + ids=1:length(Y_t0_R0[which(M_t0_R0==1)]), + kernel="dnorm", bw.selection="ise", bw.method="optim", + use_mave=TRUE), silent = TRUE) + if (inherits(fit_t0_R0_h, "try-error")) { + fit_t0_R0_h <- try(fit_SensIAT_single_index_norm1coef_model(X = X_adjust_scale_t0_R0[which(M_t0_R0==1), ], + Y = Y_t0_R0[which(M_t0_R0==1)], + ids=1:length(Y_t0_R0[which(M_t0_R0==1)]), + kernel="dnorm", bw.selection="ise", bw.method="optim", + use_mave=FALSE), silent = TRUE) + } + + return(list(fit_t_R0_h=fit_t_R0_h, fit_t_R1_h=fit_t_R1_h, fit_t0_R0_h=fit_t0_R0_h, + fit_t0_R1_h=fit_t0_R1_h, g.fit=g.fit, t_R0.fit=t_R0.fit, t_R1.fit=t_R1.fit, + M_R0.fit=M_R0.fit, M_R1.fit=M_R1.fit)) + +} + +## 2nd step of parametric bootstrap: simulate new data +pboot_sim <- function(data, X_sim, sim.size, seed, + fit_t_R0_h, fit_t_R1_h, fit_t0_R0_h, + fit_t0_R1_h, g.fit, t_R0.fit, M_R0.fit, M_R1.fit){ + + ## fit models + index.var.Y <- paste(single.index.variables(X_sim), collapse = "+") + X_adjust <- model.matrix(as.formula(paste("~", index.var.Y)), data = X_sim)[,-1] + X_adjust_scale <- scale(X_adjust) + + X_adjust_scale_t_R0 <- X_adjust_scale[which(data$t==1 & data$R==0), ] + X_adjust_scale_t_R1 <- X_adjust_scale[which(data$t==1 & data$R==1), ] + X_adjust_scale_t0_R0 <- X_adjust_scale[which(data$t==0 & data$R==0), ] + X_adjust_scale_t0_R1 <- X_adjust_scale[which(data$t==0 & data$R==1), ] + M_t_R0 <- data$M[which(data$t==1 & data$R==0)] + M_t_R1 <- data$M[which(data$t==1 & data$R==1)] + M_t0_R0 <- data$M[which(data$t==0 & data$R==0)] + M_t0_R1 <- data$M[which(data$t==0 & data$R==1)] + Y_t_R0 <- data$Y[which(data$t==1 & data$R==0)] + Y_t_R1 <- data$Y[which(data$t==1 & data$R==1)] + Y_t0_R0 <- data$Y[which(data$t==0 & data$R==0)] + Y_t0_R1 <- data$Y[which(data$t==0 & data$R==1)] + + X_t_R0_beta_t_R0 <- as.vector(X_adjust_scale_t_R0[which(M_t_R0==1), ] %*% fit_t_R0_h$coef) + X_t_R1_beta_t_R1 <- as.vector(X_adjust_scale_t_R1[which(M_t_R1==1), ] %*% fit_t_R1_h$coef) + X_t0_R0_beta_t0_R0 <- as.vector(X_adjust_scale_t0_R0[which(M_t0_R0==1), ] %*% fit_t0_R0_h$coef) + X_t0_R1_beta_t0_R1 <- as.vector(X_adjust_scale_t0_R1[which(M_t0_R1==1), ] %*% fit_t0_R1_h$coef) + y_t_R0 = sort(unique(Y_t_R0[which(M_t_R0==1)])) + ny_t_R0 = length(y_t_R0) + y_t_R1 = sort(unique(Y_t_R1[which(M_t_R1==1)])) + ny_t_R1 = length(y_t_R1) + y_t0_R0 = sort(unique(Y_t0_R0[which(M_t0_R0==1)])) + ny_t0_R0 = length(y_t0_R0) + y_t0_R1 = sort(unique(Y_t0_R1[which(M_t0_R1==1)])) + ny_t0_R1 = length(y_t0_R1) + + ## simulate data + if(!is.null(seed)){set.seed(seed)} + boot.indices <- sample(1:sim.size, size=sim.size, replace = T) + X.new <- X_sim[boot.indices,] + X_adjust_scale.new <- X_adjust_scale[boot.indices, ] + ## new R + prob.R <- predict(g.fit, newdata=X.new, type="response") + R.new <- rbinom(n=sim.size, size=1, prob=prob.R) + ## new T + T.new_R1 <- rbinom(n=sim.size,size=1,prob=0.5) + prob.T_R0 <- predict(t_R0.fit, newdata=X.new, type="response") + T.new_R0 <- rbinom(n=sim.size,size=1,prob=prob.T_R0) + T.new <- T.new_R1*R.new+T.new_R0*(1-R.new) + + ## new Y + X.new_beta_t_R0 <- as.vector(X_adjust_scale.new %*% fit_t_R0_h$coef) + X.new_beta_t_R1 <- as.vector(X_adjust_scale.new %*% fit_t_R1_h$coef) + X.new_beta_t0_R0 <- as.vector(X_adjust_scale.new %*% fit_t0_R0_h$coef) + X.new_beta_t0_R1 <- as.vector(X_adjust_scale.new %*% fit_t0_R1_h$coef) + + F_X_t_R0 <- NW_new(Xb=X_t_R0_beta_t_R0, Y=Y_t_R0[which(M_t_R0==1)], + xb=X.new_beta_t_R0, y=y_t_R0, h=fit_t_R0_h$bandwidth, + kernel = "dnorm") + F_X_t_R1 <- NW_new(Xb=X_t_R1_beta_t_R1, Y=Y_t_R1[which(M_t_R1==1)], + xb=X.new_beta_t_R1, y=y_t_R1, h=fit_t_R1_h$bandwidth, + kernel = "dnorm") + F_X_t0_R0 <- NW_new(Xb=X_t0_R0_beta_t0_R0, Y=Y_t0_R0[which(M_t0_R0==1)], + xb=X.new_beta_t0_R0, y=y_t0_R0, h=fit_t0_R0_h$bandwidth, + kernel = "dnorm") + F_X_t0_R1 <- NW_new(Xb=X_t0_R1_beta_t0_R1, Y=Y_t0_R1[which(M_t0_R1==1)], + xb=X.new_beta_t0_R1, y=y_t0_R1, h=fit_t0_R1_h$bandwidth, + kernel = "dnorm") + + i1 = which(apply(F_X_t_R0==0,1,prod)==1) + i1.closest <- apply(abs(outer(X.new_beta_t_R0[i1], X.new_beta_t_R0[-i1], FUN = "-")), 1, which.min) + F_X_t_R0[i1, ] <- F_X_t_R0[-i1, ][i1.closest, ] + + i1 = which(apply(F_X_t_R1==0,1,prod)==1) + i1.closest <- apply(abs(outer(X.new_beta_t_R1[i1], X.new_beta_t_R1[-i1], FUN = "-")), 1, which.min) + F_X_t_R1[i1, ] <- F_X_t_R1[-i1, ][i1.closest, ] + + i1 = which(apply(F_X_t0_R0==0,1,prod)==1) + i1.closest <- apply(abs(outer(X.new_beta_t0_R0[i1], X.new_beta_t0_R0[-i1], FUN = "-")), 1, which.min) + F_X_t0_R0[i1, ] <- F_X_t0_R0[-i1, ][i1.closest, ] + + i1 = which(apply(F_X_t0_R1==0,1,prod)==1) + i1.closest <- apply(abs(outer(X.new_beta_t0_R1[i1], X.new_beta_t0_R1[-i1], FUN = "-")), 1, which.min) + F_X_t0_R1[i1, ] <- F_X_t0_R1[-i1, ][i1.closest, ] + + Y.new_t_R0 <- apply(F_X_t_R0, 1, function(x){ + ecdf <- runif(1,0,1) + y_t_R0[min(which(ecdf<=x))] + }) + Y.new_t_R1 <- apply(F_X_t_R1, 1, function(x){ + ecdf <- runif(1,0,1) + y_t_R1[min(which(ecdf<=x))] + }) + Y.new_t0_R0 <- apply(F_X_t0_R0, 1, function(x){ + ecdf <- runif(1,0,1) + y_t0_R0[min(which(ecdf<=x))] + }) + Y.new_t0_R1 <- apply(F_X_t0_R1, 1, function(x){ + ecdf <- runif(1,0,1) + y_t0_R1[min(which(ecdf<=x))] + }) + Y.new <- Y.new_t_R0*T.new*(1-R.new)+Y.new_t_R1*T.new*R.new+ + Y.new_t0_R0*(1-T.new)*(1-R.new)+Y.new_t0_R1*(1-T.new)*R.new + + ## new M + X_with_t.new <- cbind(T.new, X.new) + colnames(X_with_t.new)[1] <- "treatment" + prob.M_R0 <- predict(M_R0.fit, newdata=X_with_t.new, type="response") + prob.M_R1 <- predict(M_R1.fit, newdata=X_with_t.new, type="response") + M.new_R0 <- rbinom(n=sim.size,size=1,prob=prob.M_R0) + M.new_R1 <- rbinom(n=sim.size,size=1,prob=prob.M_R1) + M.new <- M.new_R0*(1-R.new)+M.new_R1*R.new + + Y.new[M.new==0] <- NA + + data.new <- data.frame(cbind(Y.new, M.new, R.new, T.new, X.new)) + colnames(data.new) <- c("Y", "M", "R", "t", colnames(X.new)) + return(data.new) } + ## run algorithm -fit_one <- function(data, X, trt_val) { +fit_one <- function(data, X, trt_val, coef_g.fit=NULL, coef_t_R0.fit=NULL, + coef_t_R1.fit=NULL, coef_M_R0.fit=NULL, coef_M_R1.fit=NULL) { # try mave first - out <- try(est_psi(Y=data$Y, M=data$M, R=data$R, - X=X, + out <- try(est_psi(Y=data$Y, M=data$M, R=data$R, X=X, t=data$t, trt=trt_val, gamma=seq(-5, 5, by=1), fold=5, IF_output=FALSE, simple_trunc=FALSE, quant=NULL, kernel="dnorm", method="optim", single_index_method="norm1coef", - use_mave=TRUE, seed=NULL), + use_mave=TRUE, seed=NULL, coef_g.fit=coef_g.fit, coef_t_R0.fit=coef_t_R0.fit, + coef_t_R1.fit=coef_t_R1.fit, coef_M_R0.fit=coef_M_R0.fit, coef_M_R1.fit=coef_M_R1.fit), silent = TRUE) if (inherits(out, "try-error")) { # fallback cumSIR - out <- try(est_psi(Y=data$Y, M=data$M, R=data$R, - X=X, + out <- try(est_psi(Y=data$Y, M=data$M, R=data$R, X=X, t=data$t, trt=trt_val, gamma=seq(-5, 5, by=1), fold=5, IF_output=FALSE, simple_trunc=FALSE, quant=NULL, kernel="dnorm", method="optim", single_index_method="norm1coef", - use_mave=FALSE, seed=NULL), + use_mave=FALSE, seed=NULL, coef_g.fit=coef_g.fit, coef_t_R0.fit=coef_t_R0.fit, + coef_t_R1.fit=coef_t_R1.fit, coef_M_R0.fit=coef_M_R0.fit, coef_M_R1.fit=coef_M_R1.fit), silent = TRUE) } diff --git a/est_s_t_y.R b/est_s_t_y.R index b6bd976..b7bbe53 100644 --- a/est_s_t_y.R +++ b/est_s_t_y.R @@ -1,7 +1,7 @@ est_psi <- function(Y, M, R, X, t, trt, gamma, fold, seed, IF_output, simple_trunc, quant, kernel, method, single_index_method, - use_mave=TRUE){ + use_mave=TRUE, coef_g.fit=NULL, coef_t_R0.fit=NULL, coef_M_R0.fit=NULL, coef_M_R1.fit=NULL){ n <- length(t) Y[is.na(Y)] <- 0 @@ -97,9 +97,9 @@ est_psi <- function(Y, M, R, X, t, trt, gamma, fold, seed, IF_output, X_out_fold <- X[out_fold_id_list, ] X_with_T_out_fold <- X_with_T[out_fold_id_list, ] - trt.ind_out_fold_R0 <- trt.ind_out_fold[which(R_out_fold==0)] + t_out_fold_R0 <- t_out_fold[which(R_out_fold==0)] X_out_fold_R0 <- X_out_fold[which(R_out_fold==0), ] - trt.ind_out_fold_R1 <- trt.ind_out_fold[which(R_out_fold==1)] + t_out_fold_R1 <- t_out_fold[which(R_out_fold==1)] X_out_fold_R1 <- X_out_fold[which(R_out_fold==1), ] M_out_fold_R0 <- M_out_fold[which(R_out_fold==0)] @@ -124,17 +124,25 @@ est_psi <- function(Y, M, R, X, t, trt, gamma, fold, seed, IF_output, X_with_T_in_fold_R1 <- X_with_T_in_fold[which(R_in_fold==1), ] ## fit models - t_R0.fit <- mgcv::gam(as.formula(paste("trt.ind_out_fold_R0 ~", gam.var)), data=X_out_fold_R0, family=binomial) ## treatment model - t_R1.fit <- mgcv::gam(as.formula(paste("trt.ind_out_fold_R1 ~", gam.var)), data=X_out_fold_R1, family=binomial) ## treatment model - - M_R0.fit <- mgcv::gam(as.formula(paste("M_out_fold_R0 ~", gam.var.M)), data=X_with_T_out_fold_R0, family=binomial) ## missing data model - M_R1.fit <- mgcv::gam(as.formula(paste("M_out_fold_R1 ~", gam.var.M)), data=X_with_T_out_fold_R1, family=binomial) ## missing data model + t_R0.fit <- mgcv::gam(as.formula(paste("t_out_fold_R0 ~", gam.var)), data=X_out_fold_R0, + family=binomial, start=coef_t_R0.fit) ## treatment model + t_R1.fit <- mgcv::gam(as.formula(paste("t_out_fold_R1 ~", gam.var)), data=X_out_fold_R1, + family=binomial, start=coef_t_R1.fit) ## treatment model + M_R0.fit <- mgcv::gam(as.formula(paste("M_out_fold_R0 ~", gam.var.M)), data=X_with_T_out_fold_R0, + family=binomial, start=coef_M_R0.fit) ## missing data model + M_R1.fit <- mgcv::gam(as.formula(paste("M_out_fold_R1 ~", gam.var.M)), data=X_with_T_out_fold_R1, + family=binomial, start=coef_M_R1.fit) ## missing data model prop.R1 <- mean(R_in_fold) ## get predictions for pi - pi_R0 <- predict(t_R0.fit, newdata=X_in_fold_t_R0, type="response") - pi_R1 <- predict(t_R1.fit, newdata=X_in_fold_R1, type="response") + if(trt==1){ + pi_R0 <- predict(t_R0.fit, newdata=X_in_fold_t_R0, type="response") + pi_R1 <- predict(t_R1.fit, newdata=X_in_fold_R1, type="response") + }else{ + pi_R0 <- 1-predict(t_R0.fit, newdata=X_in_fold_t_R0, type="response") + pi_R1 <- 1-predict(t_R1.fit, newdata=X_in_fold_R1, type="response") + } pi_R0_l <- c(pi_R0_l, pi_R0) pi_R1_l <- c(pi_R1_l, pi_R1) fold_index_pi_R0_l <- c(fold_index_pi_R0_l, rep(k, length(pi_R0))) diff --git a/simulation_est_slurm_double_pboot.R b/simulation_est_slurm_double_pboot.R index 4dc220f..39667cd 100644 --- a/simulation_est_slurm_double_pboot.R +++ b/simulation_est_slurm_double_pboot.R @@ -34,19 +34,34 @@ X_sim = data.frame(dplyr::select(data, c(age, pain_bq, expectationb, ChronicPain for(i in 1:nboot_B){ - data.new <- pboot(data=data, X_sim=X_sim, sim.size=500, seed=NULL) + boot_b_temp <- pboot(data=data, X_sim=X_sim, sim.size=500, seed=NULL) + data.new <- boot_b_temp$data.new X.new <- data.frame(dplyr::select(data.new, c(age, pain_bq, expectationb, ChronicPainb))) - topical_vals <- fit_one(data=data.new, X=X.new, trt_val=1) - oral_vals <- fit_one(data=data.new, X=X.new, trt_val=0) + topical_vals <- fit_one(data=data.new, X=X.new, trt_val=1, coef_g.fit=boot_b_temp$coef_g.fit, + coef_t_R0.fit=boot_b_temp$coef_t_R0.fit, coef_t_R1.fit=boot_b_temp$coef_t_R1.fit, + coef_M_R0.fit=boot_b_temp$coef_M_R0.fit, coef_M_R1.fit=boot_b_temp$coef_M_R1.fit) + oral_vals <- fit_one(data=data.new, X=X.new, trt_val=0, coef_g.fit=boot_b_temp$coef_g.fit, + coef_t_R0.fit=boot_b_temp$coef_t_R0.fit, coef_t_R1.fit=boot_b_temp$coef_t_R1.fit, + coef_M_R0.fit=boot_b_temp$coef_M_R0.fit, coef_M_R1.fit=boot_b_temp$coef_M_R1.fit) Boot.Est[which(grid$b==i&grid$c==1), ] <- c(topical_vals, oral_vals) + boot_b_fit <- pboot_model(data=data.new, X_sim=X.new) + for(j in 1:nboot_C){ - data.new2 <- pboot(data=data.new, X_sim=X.new, sim.size=500, seed=NULL) + data.new2 <- pboot_sim(data=data.new, X_sim=X.new, sim.size=500, seed=NULL, + fit_t_R0_h=boot_b_fit$fit_t_R0_h, fit_t_R1_h=boot_b_fit$fit_t_R1_h, + fit_t0_R0_h=boot_b_fit$fit_t0_R0_h, fit_t0_R1_h=boot_b_fit$fit_t0_R1_h, + g.fit=boot_b_fit$g.fit, t_R0.fit=boot_b_fit$t_R0.fit, M_R0.fit=boot_b_fit$M_R0.fit, + M_R1.fit=boot_b_fit$M_R1.fit) X.new2 <- data.frame(dplyr::select(data.new2, c(age, pain_bq, expectationb, ChronicPainb))) ## inner estimator - temp_topical <- fit_one(data=data.new2, X=X.new2, trt_val=1) - temp_oral <- fit_one(data=data.new2, X=X.new2, trt_val=0) + temp_topical <- fit_one(data=data.new2, X=X.new2, trt_val=1, coef_g.fit=coef(boot_b_fit$g.fit), + coef_t_R0.fit=coef(boot_b_fit$t_R0.fit), coef_t_R1.fit=coef(boot_b_fit$t_R1.fit), + coef_M_R0.fit=coef(boot_b_fit$M_R0.fit), coef_M_R1.fit=coef(boot_b_fit$M_R1.fit)) + temp_oral <- fit_one(data=data.new2, X=X.new2, trt_val=0, coef_g.fit=coef(boot_b_fit$g.fit), + coef_t_R0.fit=coef(boot_b_fit$t_R0.fit), coef_t_R1.fit=coef(boot_b_fit$t_R1.fit), + coef_M_R0.fit=coef(boot_b_fit$M_R0.fit), coef_M_R1.fit=coef(boot_b_fit$M_R1.fit)) Boot.Est[which(grid$b==i&grid$c==j), 133:264] <- c(temp_topical, temp_oral) } From c122edeab3f457377e69f83372a1412128cf0ec1 Mon Sep 17 00:00:00 2001 From: jenniferxu1101 Date: Thu, 12 Feb 2026 14:26:49 -0500 Subject: [PATCH 10/19] profvis with initial coef for gam --- profvis_initial_coef.html | 4637 +++++++++++++++++++++++++++++++++++++ 1 file changed, 4637 insertions(+) create mode 100644 profvis_initial_coef.html diff --git a/profvis_initial_coef.html b/profvis_initial_coef.html new file mode 100644 index 0000000..b1b1400 --- /dev/null +++ b/profvis_initial_coef.html @@ -0,0 +1,4637 @@ + + + + +profvis + + + + + + + + + + + + + +
+
+
+ + + + From 4c8a1f6870e7a20d166bdad9dbec1987d995c3e7 Mon Sep 17 00:00:00 2001 From: jenniferxu1101 Date: Mon, 16 Feb 2026 10:44:02 -0500 Subject: [PATCH 11/19] Add initial coef --- est_s_t_y.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/est_s_t_y.R b/est_s_t_y.R index b7bbe53..a71baff 100644 --- a/est_s_t_y.R +++ b/est_s_t_y.R @@ -1,7 +1,8 @@ est_psi <- function(Y, M, R, X, t, trt, gamma, fold, seed, IF_output, simple_trunc, quant, kernel, method, single_index_method, - use_mave=TRUE, coef_g.fit=NULL, coef_t_R0.fit=NULL, coef_M_R0.fit=NULL, coef_M_R1.fit=NULL){ + use_mave=TRUE, coef_g.fit=NULL, coef_t_R0.fit=NULL, coef_t_R1.fit=NULL, + coef_M_R0.fit=NULL, coef_M_R1.fit=NULL){ n <- length(t) Y[is.na(Y)] <- 0 From fe28837e0277433a1839f507a00f03e1b4cb3cf6 Mon Sep 17 00:00:00 2001 From: jenniferxu1101 Date: Mon, 16 Feb 2026 10:44:29 -0500 Subject: [PATCH 12/19] Run double bootstrap with for-loop --- simulation_est_slurm_double_pboot2.R | 239 +++++++++++++++++++++++++++ 1 file changed, 239 insertions(+) create mode 100644 simulation_est_slurm_double_pboot2.R diff --git a/simulation_est_slurm_double_pboot2.R b/simulation_est_slurm_double_pboot2.R new file mode 100644 index 0000000..d8bca2c --- /dev/null +++ b/simulation_est_slurm_double_pboot2.R @@ -0,0 +1,239 @@ +################### +## load packages ## +################### +library(tidyverse) +library(splines) +library(gridExtra) +library(parallel) +library(foreach) +library(doParallel) + +source("/uufs/chpc.utah.edu/common/home/u6070035/CCS/code/HelperFunction.R") +source("/uufs/chpc.utah.edu/common/home/u6070035/CCS/code/singleindexmodelfunctions.R") +source("/uufs/chpc.utah.edu/common/home/u6070035/CCS/code/est_s_t_y.R") +source("/uufs/chpc.utah.edu/common/home/u6070035/CCS/code/SIDR_Ravinew.R") +source("/uufs/chpc.utah.edu/common/home/u6070035/CCS/code/SIDRnew.R") +source("/uufs/chpc.utah.edu/common/home/u6070035/CCS/code/SensIAT_sim_outcome_modeler_mave.R") + +ntasks <- Sys.getenv("SLURM_NTASKS") +if (ntasks == '') { + ntasks <- 4 +} else { + ntasks <- strtoi(ntasks) } +cat("This script use ", ntasks, " cores\n") + +package_list <- c("splines", "tidyverse", "gridExtra") + +nboot_B <- 500 + +################ +## simulation ## +################ + +cl <- makeCluster(ntasks) +registerDoParallel(cl) + +for (j in (40*(index-1)+1):(40*index)){ + +Boot.Est <- foreach(i=1:nboot_B, .combine=rbind, .packages=package_list, .errorhandling = "pass") %dopar% { + + gamma_length <- length(seq(-5, 5, by=1)) + nboot_C <- 100 + + load(paste0("/uufs/chpc.utah.edu/common/home/u6070035/CCS/simData/imputed", imputed, "/n", sim.size, "/sim", j, ".rds")) + X_sim = data.frame(dplyr::select(data, c(age, pain_bq, expectationb, ChronicPainb))) + + ## simulate outer bootstrap sample + boot_b_temp <- pboot(data=data, X_sim=X_sim, sim.size=sim.size, seed=NULL) + data.new <- boot_b_temp$data.new + + X.new <- data.frame(dplyr::select(data.new, c(age, pain_bq, expectationb, ChronicPainb))) + topical_vals <- fit_one(data=data.new, X=X.new, trt_val=1, coef_g.fit=boot_b_temp$coef_g.fit, + coef_t_R0.fit=boot_b_temp$coef_t_R0.fit, coef_t_R1.fit=boot_b_temp$coef_t_R1.fit, + coef_M_R0.fit=boot_b_temp$coef_M_R0.fit, coef_M_R1.fit=boot_b_temp$coef_M_R1.fit) + oral_vals <- fit_one(data=data.new, X=X.new, trt_val=0, coef_g.fit=boot_b_temp$coef_g.fit, + coef_t_R0.fit=boot_b_temp$coef_t_R0.fit, coef_t_R1.fit=boot_b_temp$coef_t_R1.fit, + coef_M_R0.fit=boot_b_temp$coef_M_R0.fit, coef_M_R1.fit=boot_b_temp$coef_M_R1.fit) + boot_b_fit <- pboot_model(data=data.new, X_sim=X.new) + + load(paste0("/uufs/chpc.utah.edu/common/home/u6070035/CCS/simResult/imputed", imputed, "/n", sim.size, "/", kernel, "_", single_index_method, "_fold5/sim", j, ".RData")) + + t_b_topical <- (topical_vals[1:11]-topical_WOMAC_12m$est_trunc)/sqrt(topical_vals[34:44]) + t_b_R1_topical <- (topical_vals[12:22]-topical_WOMAC_12m$est_trunc_R1)/sqrt(topical_vals[45:55]) + t_b_R0_topical <- (topical_vals[23:33]-topical_WOMAC_12m$est_trunc_R0)/sqrt(topical_vals[56:66]) + t_b_oral <- (oral_vals[1:11]-oral_WOMAC_12m$est_trunc)/sqrt(oral_vals[34:44]) + t_b_R1_oral <- (oral_vals[12:22]-oral_WOMAC_12m$est_trunc_R1)/sqrt(oral_vals[45:55]) + t_b_R0_oral <- (oral_vals[23:33]-oral_WOMAC_12m$est_trunc_R0)/sqrt(oral_vals[56:66]) + + ## empty containers + Q_b_topical <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + Q_b_R1_topical <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + Q_b_R0_topical <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + Q_b_oral <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + Q_b_R1_oral <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + t_bc_topical <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + t_bc_R1_topical <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + t_bc_R0_topical <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + t_bc_oral <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + t_bc_R1_oral <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + t_bc_R0_oral <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + est_topical <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + est_R1_topical <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + est_R0_topical <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + est_oral <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + est_R1_oral <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + est_R0_oral <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + + for(j in 1:nboot_C){ + + data.new2 <- pboot_sim(data=data.new, X_sim=X.new, sim.size=sim.size, seed=NULL, + fit_t_R0_h=boot_b_fit$fit_t_R0_h, fit_t_R1_h=boot_b_fit$fit_t_R1_h, + fit_t0_R0_h=boot_b_fit$fit_t0_R0_h, fit_t0_R1_h=boot_b_fit$fit_t0_R1_h, + g.fit=boot_b_fit$g.fit, t_R0.fit=boot_b_fit$t_R0.fit, M_R0.fit=boot_b_fit$M_R0.fit, + M_R1.fit=boot_b_fit$M_R1.fit) + X.new2 <- data.frame(dplyr::select(data.new2, c(age, pain_bq, expectationb, ChronicPainb))) + ## inner estimator + temp_topical <- fit_one(data=data.new2, X=X.new2, trt_val=1, coef_g.fit=coef(boot_b_fit$g.fit), + coef_t_R0.fit=coef(boot_b_fit$t_R0.fit), coef_t_R1.fit=coef(boot_b_fit$t_R1.fit), + coef_M_R0.fit=coef(boot_b_fit$M_R0.fit), coef_M_R1.fit=coef(boot_b_fit$M_R1.fit)) + temp_oral <- fit_one(data=data.new2, X=X.new2, trt_val=0, coef_g.fit=coef(boot_b_fit$g.fit), + coef_t_R0.fit=coef(boot_b_fit$t_R0.fit), coef_t_R1.fit=coef(boot_b_fit$t_R1.fit), + coef_M_R0.fit=coef(boot_b_fit$M_R0.fit), coef_M_R1.fit=coef(boot_b_fit$M_R1.fit)) + + Q_b_topical[j, ] <- as.numeric(temp_topical[1:11]<=topical_WOMAC_12m$est_trunc) + Q_b_R1_topical[j, ] <- as.numeric(temp_topical[12:22]<=topical_WOMAC_12m$est_trunc_R1) + Q_b_R0_topical[j, ] <- as.numeric(temp_topical[23:33]<=topical_WOMAC_12m$est_trunc_R0) + Q_b_oral[j, ] <- as.numeric(temp_oral[1:11]<=oral_WOMAC_12m$est_trunc) + Q_b_R1_oral[j, ] <- as.numeric(temp_oral[12:22]<=oral_WOMAC_12m$est_trunc_R1) + Q_b_R0_oral[j, ] <- as.numeric(temp_oral[23:33]<=oral_WOMAC_12m$est_trunc_R0) + + t_bc_topical[j, ] <- (temp_topical[1:11]-topical_vals[1:11])/sqrt(temp_topical[34:44]) + t_bc_R1_topical[j, ] <- (temp_topical[12:22]-topical_vals[12:22])/sqrt(temp_topical[45:55]) + t_bc_R0_topical[j, ] <- (temp_topical[23:33]-topical_vals[23:33])/sqrt(temp_topical[56:66]) + t_bc_oral[j, ] <- (temp_oral[1:11]-oral_vals[1:11])/sqrt(temp_topical[34:44]) + t_bc_R1_oral[j, ] <- (temp_oral[12:22]-oral_vals[12:22])/sqrt(temp_oral[45:55]) + t_bc_R0_oral[j, ] <- (temp_oral[23:33]-oral_vals[23:33])/sqrt(temp_oral[56:66]) + + est_topical[j, ] <- temp_topical[1:11] + est_R1_topical[j, ] <- temp_topical[12:22] + est_R0_topical[j, ] <- temp_topical[23:33] + est_oral[j, ] <- temp_oral[1:11] + est_R1_oral[j, ] <- temp_oral[12:22] + est_R0_oral[j, ] <- temp_oral[23:33] + + } + + r_Q_b_topical <- colMeans(Q_b_topical, na.rm=T) + r_Q_b_R1_topical <- colMeans(Q_b_R1_topical, na.rm=T) + r_Q_b_R0_topical <- colMeans(Q_b_R0_topical, na.rm=T) + r_Q_b_oral <- colMeans(Q_b_oral, na.rm=T) + r_Q_b_R1_oral <- colMeans(Q_b_R1_oral, na.rm=T) + r_Q_b_R0_oral <- colMeans(Q_b_R0_oral, na.rm=T) + Q_t_b_topical <- colMeans(t_bc_topical<=matrix(t_b_topical, nrow=nboot_C, ncol=gamma_length, + byrow=T), na.rm=T) + abs_Q_t_b_topical <- colMeans(abs(t_bc_topical)<=matrix(abs(t_b_topical), nrow=nboot_C, ncol=gamma_length, + byrow=T), na.rm=T) + Q_t_b_R1_topical <- colMeans(t_bc_R1_topical<=matrix(t_b_R1_topical, nrow=nboot_C, ncol=gamma_length, + byrow=T), na.rm=T) + abs_Q_t_b_R1_topical <- colMeans(abs(t_bc_R1_topical)<=matrix(abs(t_b_R1_topical), nrow=nboot_C, ncol=gamma_length, + byrow=T), na.rm=T) + Q_t_b_R0_topical <- colMeans(t_bc_R0_topical<=matrix(t_b_R0_topical, nrow=nboot_C, ncol=gamma_length, + byrow=T), na.rm=T) + abs_Q_t_b_R0_topical <- colMeans(abs(t_bc_R0_topical)<=matrix(abs(t_b_R0_topical), nrow=nboot_C, ncol=gamma_length, + byrow=T), na.rm=T) + Q_t_b_oral <- colMeans(t_bc_oral<=matrix(t_b_oral, nrow=nboot_C, ncol=gamma_length, + byrow=T), na.rm=T) + abs_Q_t_b_oral <- colMeans(abs(t_bc_oral)<=matrix(abs(t_b_oral), nrow=nboot_C, ncol=gamma_length, + byrow=T), na.rm=T) + Q_t_b_R1_oral <- colMeans(t_bc_R1_oral<=matrix(t_b_R1_oral, nrow=nboot_C, ncol=gamma_length, + byrow=T), na.rm=T) + abs_Q_t_b_R1_oral <- colMeans(abs(t_bc_R1_oral)<=matrix(abs(t_b_R1_oral), nrow=nboot_C, ncol=gamma_length, + byrow=T), na.rm=T) + Q_t_b_R0_oral <- colMeans(t_bc_R0_oral<=matrix(t_b_R0_oral, nrow=nboot_C, ncol=gamma_length, + byrow=T), na.rm=T) + abs_Q_t_b_R0_oral <- colMeans(abs(t_bc_R0_oral)<=matrix(abs(t_b_R0_oral), nrow=nboot_C, ncol=gamma_length, + byrow=T), na.rm=T) + t_b_sd_topical <- (topical_vals[1:11]-topical_WOMAC_12m$est_trunc)/apply(est_topical, 2, function(x){sd(x, na.rm=TRUE)}) + t_b_R1_sd_topical <- (topical_vals[12:22]-topical_WOMAC_12m$est_trunc_R1)/apply(est_R1_topical, 2, function(x){sd(x, na.rm=TRUE)}) + t_b_R0_sd_topical <- (topical_vals[23:33]-topical_WOMAC_12m$est_trunc_R0)/apply(est_R0_topical, 2, function(x){sd(x, na.rm=TRUE)}) + t_b_sd_oral <- (oral_vals[1:11]-oral_WOMAC_12m$est_trunc)/apply(est_oral, 2, function(x){sd(x, na.rm=TRUE)}) + t_b_R1_sd_oral <- (oral_vals[12:22]-oral_WOMAC_12m$est_trunc_R1)/apply(est_R1_oral, 2, function(x){sd(x, na.rm=TRUE)}) + t_b_R0_sd_oral <- (oral_vals[23:33]-oral_WOMAC_12m$est_trunc_R0)/apply(est_R0_oral, 2, function(x){sd(x, na.rm=TRUE)}) + +return(c(topical_vals, oral_vals, r_Q_b_topical, r_Q_b_R1_topical, r_Q_b_R0_topical, + r_Q_b_oral, r_Q_b_R1_oral, r_Q_b_R0_oral, Q_t_b_topical, Q_t_b_R1_topical, Q_t_b_R0_topical, + Q_t_b_oral, Q_t_b_R1_oral, Q_t_b_R0_oral, abs_Q_t_b_topical, abs_Q_t_b_R1_topical, abs_Q_t_b_R0_topical, + abs_Q_t_b_oral, abs_Q_t_b_R1_oral, abs_Q_t_b_R0_oral, t_b_sd_topical, t_b_R1_sd_topical, t_b_R0_sd_topical, + t_b_sd_oral, t_b_R1_sd_oral, t_b_R0_sd_oral)) + +} + +est_topical_boot <- Boot.Est[, 1:11] +est_R1_topical_boot <- Boot.Est[, 12:22] +est_R0_topical_boot <- Boot.Est[, 23:33] +var_topical_boot <- Boot.Est[, 34:44] +var_R1_topical_boot <- Boot.Est[, 45:55] +var_R0_topical_boot <- Boot.Est[, 56:66] + +est_oral_boot <- Boot.Est[, 67:77] +est_R1_oral_boot <- Boot.Est[, 78:88] +est_R0_oral_boot <- Boot.Est[, 89:99] +var_oral_boot <- Boot.Est[, 100:110] +var_R1_oral_boot <- Boot.Est[, 111:121] +var_R0_oral_boot <- Boot.Est[, 122:132] + +Q_b_topical <- Boot.Est[, 133:143] +Q_b_R1_topical <- Boot.Est[, 144:154] +Q_b_R0_topical <- Boot.Est[, 155:165] +Q_b_oral <- Boot.Est[, 166:176] +Q_b_R1_oral <- Boot.Est[, 177:187] +Q_b_R0_oral <- Boot.Est[, 188:198] + +Q_t_b_topical <- Boot.Est[, 199:209] +Q_t_b_R1_topical <- Boot.Est[, 210:220] +Q_t_b_R0_topical <- Boot.Est[, 221:231] +Q_t_b_oral <- Boot.Est[, 232:242] +Q_t_b_R1_oral <- Boot.Est[, 243:253] +Q_t_b_R0_oral <- Boot.Est[, 254:264] + +abs_Q_t_b_topical <- Boot.Est[, 265:275] +abs_Q_t_b_R1_topical <- Boot.Est[, 276:286] +abs_Q_t_b_R0_topical <- Boot.Est[, 287:297] +abs_Q_t_b_oral <- Boot.Est[, 298:308] +abs_Q_t_b_R1_oral <- Boot.Est[, 309:319] +abs_Q_t_b_R0_oral <- Boot.Est[, 320:330] + +t_b_sd_topical <- Boot.Est[, 331:341] +t_b_R1_sd_topical <- Boot.Est[, 342:352] +t_b_R0_sd_topical <- Boot.Est[, 353:363] +t_b_sd_oral <- Boot.Est[, 364:374] +t_b_R1_sd_oral <- Boot.Est[, 375:385] +t_b_R0_sd_oral <- Boot.Est[, 386:396] + +## bootstrap result +save(est_topical_boot, est_R1_topical_boot, est_R0_topical_boot, + est_oral_boot, est_R1_oral_boot, est_R0_oral_boot, + var_topical_boot, var_R1_topical_boot, var_R0_topical_boot, + var_oral_boot, var_R1_oral_boot, var_R0_oral_boot, + Q_b_topical, Q_b_R1_topical, Q_b_R0_topical, + Q_b_oral, Q_b_R1_oral, Q_b_R0_oral, + t_b_topical, t_b_R1_topical, t_b_R0_topical, + t_b_oral, t_b_R1_oral, t_b_R0_oral, + Q_t_b_topical, Q_t_b_R1_topical, Q_t_b_R0_topical, + Q_t_b_oral, Q_t_b_R1_oral, Q_t_b_R0_oral, + abs_Q_t_b_topical, abs_Q_t_b_R1_topical, abs_Q_t_b_R0_topical, + abs_Q_t_b_oral, abs_Q_t_b_R1_oral, abs_Q_t_b_R0_oral, + t_b_sd_topical, t_b_R1_sd_topical, t_b_R0_sd_topical, + t_b_sd_oral, t_b_R1_sd_oral, t_b_R0_sd_oral, + file=paste0("/uufs/chpc.utah.edu/common/home/u6070035/CCS/simResult/imputed", imputed, "/double_pboot_n", sim.size, "/", kernel, "_", single_index_method, "/sim", j, ".RData")) +print(j) + +} + +stopCluster(cl) + + + + + + From 44248608454f2133c19714eb3fe42a2bed89f3a8 Mon Sep 17 00:00:00 2001 From: jenniferxu1101 Date: Mon, 16 Feb 2026 10:57:31 -0500 Subject: [PATCH 13/19] Add initial coef for gam --- profvis_initial_coef.html | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/profvis_initial_coef.html b/profvis_initial_coef.html index b1b1400..20acb38 100644 --- a/profvis_initial_coef.html +++ b/profvis_initial_coef.html @@ -4629,9 +4629,9 @@
-
+
- - + + From 5d39a415bf6abf38ac0965313cb67aa5ecfd8910 Mon Sep 17 00:00:00 2001 From: jenniferxu1101 Date: Wed, 18 Feb 2026 09:55:28 -0500 Subject: [PATCH 14/19] smoother for gam --- HelperFunction.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/HelperFunction.R b/HelperFunction.R index dae5a10..c59dd2e 100644 --- a/HelperFunction.R +++ b/HelperFunction.R @@ -1,6 +1,6 @@ gam.variables <- function(data){ sapply(1:ncol(data), function(i){ - if (is.numeric(data[,i]) | is.integer(data[, i])){ paste("s", "(",names(data)[i], ")", sep="")} + if (is.numeric(data[,i]) | is.integer(data[, i])){ paste("s", "(",names(data)[i], ", k=5)", sep="")} else{names(data[i]) } }) } From cb421e84eb3ec3849effb19739fb169a7f4e26a2 Mon Sep 17 00:00:00 2001 From: jenniferxu1101 Date: Wed, 18 Feb 2026 10:27:48 -0500 Subject: [PATCH 15/19] coef --- simulation_est_slurm_double_pboot.R | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/simulation_est_slurm_double_pboot.R b/simulation_est_slurm_double_pboot.R index 39667cd..71e380f 100644 --- a/simulation_est_slurm_double_pboot.R +++ b/simulation_est_slurm_double_pboot.R @@ -46,6 +46,11 @@ for(i in 1:nboot_B){ Boot.Est[which(grid$b==i&grid$c==1), ] <- c(topical_vals, oral_vals) boot_b_fit <- pboot_model(data=data.new, X_sim=X.new) + coef_g.fit <- coef(boot_b_fit$g.fit) + coef_t_R0.fit <- coef(boot_b_fit$t_R0.fit) + coef_t_R1.fit <- coef(boot_b_fit$t_R1.fit) + coef_M_R0.fit <- coef(boot_b_fit$M_R0.fit) + coef_M_R1.fit <- coef(boot_b_fit$M_R1.fit) for(j in 1:nboot_C){ @@ -56,12 +61,12 @@ for(i in 1:nboot_B){ M_R1.fit=boot_b_fit$M_R1.fit) X.new2 <- data.frame(dplyr::select(data.new2, c(age, pain_bq, expectationb, ChronicPainb))) ## inner estimator - temp_topical <- fit_one(data=data.new2, X=X.new2, trt_val=1, coef_g.fit=coef(boot_b_fit$g.fit), - coef_t_R0.fit=coef(boot_b_fit$t_R0.fit), coef_t_R1.fit=coef(boot_b_fit$t_R1.fit), - coef_M_R0.fit=coef(boot_b_fit$M_R0.fit), coef_M_R1.fit=coef(boot_b_fit$M_R1.fit)) - temp_oral <- fit_one(data=data.new2, X=X.new2, trt_val=0, coef_g.fit=coef(boot_b_fit$g.fit), - coef_t_R0.fit=coef(boot_b_fit$t_R0.fit), coef_t_R1.fit=coef(boot_b_fit$t_R1.fit), - coef_M_R0.fit=coef(boot_b_fit$M_R0.fit), coef_M_R1.fit=coef(boot_b_fit$M_R1.fit)) + temp_topical <- fit_one(data=data.new2, X=X.new2, trt_val=1, coef_g.fit=coef_g.fit, + coef_t_R0.fit=coef_t_R0.fit, coef_t_R1.fit=coef_t_R1.fit, + coef_M_R0.fit=coef_M_R0.fit, coef_M_R1.fit=coef_M_R1.fit) + temp_oral <- fit_one(data=data.new2, X=X.new2, trt_val=0, coef_g.fit=coef_g.fit, + coef_t_R0.fit=coef_t_R0.fit, coef_t_R1.fit=coef_t_R1.fit, + coef_M_R0.fit=coef_M_R0.fit, coef_M_R1.fit=coef_M_R1.fit) Boot.Est[which(grid$b==i&grid$c==j), 133:264] <- c(temp_topical, temp_oral) } @@ -72,7 +77,7 @@ for(i in 1:nboot_B){ }) -htmlwidgets::saveWidget(pv, "profvis.html") +htmlwidgets::saveWidget(pv, "profvis_initial_coef.html") load(paste0("/uufs/chpc.utah.edu/common/home/u6070035/CCS/simResult/imputed", imputed, "/n", sim.size, "/", kernel, "_", single_index_method, "_fold5/sim", j, ".RData")) Q_b_topical <- sapply(1:nboot_B, function(x){ From b14a259618f101c5fbc50bdc37cc4617838427fc Mon Sep 17 00:00:00 2001 From: jenniferxu1101 Date: Wed, 18 Feb 2026 10:28:06 -0500 Subject: [PATCH 16/19] coef --- simulation_est_slurm_double_pboot2.R | 29 ++++++++++++++++------------ 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/simulation_est_slurm_double_pboot2.R b/simulation_est_slurm_double_pboot2.R index d8bca2c..28ea71d 100644 --- a/simulation_est_slurm_double_pboot2.R +++ b/simulation_est_slurm_double_pboot2.R @@ -35,14 +35,16 @@ registerDoParallel(cl) for (j in (40*(index-1)+1):(40*index)){ +load(paste0("/uufs/chpc.utah.edu/common/home/u6070035/CCS/simData/imputed", imputed, "/n", sim.size, "/sim", j, ".rds")) +X_sim = data.frame(dplyr::select(data, c(age, pain_bq, expectationb, ChronicPainb))) + +load(paste0("/uufs/chpc.utah.edu/common/home/u6070035/CCS/simResult/imputed", imputed, "/n", sim.size, "/", kernel, "_", single_index_method, "_fold5/sim", j, ".RData")) + Boot.Est <- foreach(i=1:nboot_B, .combine=rbind, .packages=package_list, .errorhandling = "pass") %dopar% { gamma_length <- length(seq(-5, 5, by=1)) nboot_C <- 100 - load(paste0("/uufs/chpc.utah.edu/common/home/u6070035/CCS/simData/imputed", imputed, "/n", sim.size, "/sim", j, ".rds")) - X_sim = data.frame(dplyr::select(data, c(age, pain_bq, expectationb, ChronicPainb))) - ## simulate outer bootstrap sample boot_b_temp <- pboot(data=data, X_sim=X_sim, sim.size=sim.size, seed=NULL) data.new <- boot_b_temp$data.new @@ -55,9 +57,12 @@ Boot.Est <- foreach(i=1:nboot_B, .combine=rbind, .packages=package_list, .errorh coef_t_R0.fit=boot_b_temp$coef_t_R0.fit, coef_t_R1.fit=boot_b_temp$coef_t_R1.fit, coef_M_R0.fit=boot_b_temp$coef_M_R0.fit, coef_M_R1.fit=boot_b_temp$coef_M_R1.fit) boot_b_fit <- pboot_model(data=data.new, X_sim=X.new) - - load(paste0("/uufs/chpc.utah.edu/common/home/u6070035/CCS/simResult/imputed", imputed, "/n", sim.size, "/", kernel, "_", single_index_method, "_fold5/sim", j, ".RData")) - + coef_g.fit <- coef(boot_b_fit$g.fit) + coef_t_R0.fit <- coef(boot_b_fit$t_R0.fit) + coef_t_R1.fit <- coef(boot_b_fit$t_R1.fit) + coef_M_R0.fit <- coef(boot_b_fit$M_R0.fit) + coef_M_R1.fit <- coef(boot_b_fit$M_R1.fit) + t_b_topical <- (topical_vals[1:11]-topical_WOMAC_12m$est_trunc)/sqrt(topical_vals[34:44]) t_b_R1_topical <- (topical_vals[12:22]-topical_WOMAC_12m$est_trunc_R1)/sqrt(topical_vals[45:55]) t_b_R0_topical <- (topical_vals[23:33]-topical_WOMAC_12m$est_trunc_R0)/sqrt(topical_vals[56:66]) @@ -93,12 +98,12 @@ Boot.Est <- foreach(i=1:nboot_B, .combine=rbind, .packages=package_list, .errorh M_R1.fit=boot_b_fit$M_R1.fit) X.new2 <- data.frame(dplyr::select(data.new2, c(age, pain_bq, expectationb, ChronicPainb))) ## inner estimator - temp_topical <- fit_one(data=data.new2, X=X.new2, trt_val=1, coef_g.fit=coef(boot_b_fit$g.fit), - coef_t_R0.fit=coef(boot_b_fit$t_R0.fit), coef_t_R1.fit=coef(boot_b_fit$t_R1.fit), - coef_M_R0.fit=coef(boot_b_fit$M_R0.fit), coef_M_R1.fit=coef(boot_b_fit$M_R1.fit)) - temp_oral <- fit_one(data=data.new2, X=X.new2, trt_val=0, coef_g.fit=coef(boot_b_fit$g.fit), - coef_t_R0.fit=coef(boot_b_fit$t_R0.fit), coef_t_R1.fit=coef(boot_b_fit$t_R1.fit), - coef_M_R0.fit=coef(boot_b_fit$M_R0.fit), coef_M_R1.fit=coef(boot_b_fit$M_R1.fit)) + temp_topical <- fit_one(data=data.new2, X=X.new2, trt_val=1, coef_g.fit=coef_g.fit, + coef_t_R0.fit=coef_t_R0.fit, coef_t_R1.fit=coef_t_R1.fit, + coef_M_R0.fit=coef_M_R0.fit, coef_M_R1.fit=coef_M_R1.fit) + temp_oral <- fit_one(data=data.new2, X=X.new2, trt_val=0, coef_g.fit=coef_g.fit, + coef_t_R0.fit=coef_t_R0.fit, coef_t_R1.fit=coef_t_R1.fit, + coef_M_R0.fit=coef_M_R0.fit, coef_M_R1.fit=coef_M_R1.fit) Q_b_topical[j, ] <- as.numeric(temp_topical[1:11]<=topical_WOMAC_12m$est_trunc) Q_b_R1_topical[j, ] <- as.numeric(temp_topical[12:22]<=topical_WOMAC_12m$est_trunc_R1) From 411dd52f25ee5cf50b9e29d6d2595360aab6c7f3 Mon Sep 17 00:00:00 2001 From: jenniferxu1101 Date: Wed, 18 Feb 2026 12:20:43 -0500 Subject: [PATCH 17/19] NA --- profvis_initial_coef.html | 4637 ------------------------------------- 1 file changed, 4637 deletions(-) delete mode 100644 profvis_initial_coef.html diff --git a/profvis_initial_coef.html b/profvis_initial_coef.html deleted file mode 100644 index 20acb38..0000000 --- a/profvis_initial_coef.html +++ /dev/null @@ -1,4637 +0,0 @@ - - - - -profvis - - - - - - - - - - - - - -
-
-
- - - - From f62baf63b5ef371289ffb0e00ac083179c96aa84 Mon Sep 17 00:00:00 2001 From: jenniferxu1101 Date: Wed, 18 Feb 2026 12:21:02 -0500 Subject: [PATCH 18/19] smooth term --- profvis_initial_coef_k5.html | 4637 ++++++++++++++++++++++++++++++++++ 1 file changed, 4637 insertions(+) create mode 100644 profvis_initial_coef_k5.html diff --git a/profvis_initial_coef_k5.html b/profvis_initial_coef_k5.html new file mode 100644 index 0000000..b79bf3d --- /dev/null +++ b/profvis_initial_coef_k5.html @@ -0,0 +1,4637 @@ + + + + +profvis + + + + + + + + + + + + + +
+
+
+ + + + From 650eed266d104fbd13c2ba746aaff46fba3370c2 Mon Sep 17 00:00:00 2001 From: jenniferxu1101 Date: Fri, 20 Feb 2026 14:44:38 -0500 Subject: [PATCH 19/19] add t statistics for outer bootstrap --- simulation_est_slurm_double_pboot2.R | 58 ++++++++++++++++++---------- 1 file changed, 37 insertions(+), 21 deletions(-) diff --git a/simulation_est_slurm_double_pboot2.R b/simulation_est_slurm_double_pboot2.R index 28ea71d..6beb837 100644 --- a/simulation_est_slurm_double_pboot2.R +++ b/simulation_est_slurm_double_pboot2.R @@ -88,8 +88,10 @@ Boot.Est <- foreach(i=1:nboot_B, .combine=rbind, .packages=package_list, .errorh est_oral <- matrix(NA, nrow=nboot_C, ncol=gamma_length) est_R1_oral <- matrix(NA, nrow=nboot_C, ncol=gamma_length) est_R0_oral <- matrix(NA, nrow=nboot_C, ncol=gamma_length) + num_NA_topical <- 0 + num_NA_oral <- 0 - for(j in 1:nboot_C){ + for(z in 1:nboot_C){ data.new2 <- pboot_sim(data=data.new, X_sim=X.new, sim.size=sim.size, seed=NULL, fit_t_R0_h=boot_b_fit$fit_t_R0_h, fit_t_R1_h=boot_b_fit$fit_t_R1_h, @@ -105,26 +107,29 @@ Boot.Est <- foreach(i=1:nboot_B, .combine=rbind, .packages=package_list, .errorh coef_t_R0.fit=coef_t_R0.fit, coef_t_R1.fit=coef_t_R1.fit, coef_M_R0.fit=coef_M_R0.fit, coef_M_R1.fit=coef_M_R1.fit) - Q_b_topical[j, ] <- as.numeric(temp_topical[1:11]<=topical_WOMAC_12m$est_trunc) - Q_b_R1_topical[j, ] <- as.numeric(temp_topical[12:22]<=topical_WOMAC_12m$est_trunc_R1) - Q_b_R0_topical[j, ] <- as.numeric(temp_topical[23:33]<=topical_WOMAC_12m$est_trunc_R0) - Q_b_oral[j, ] <- as.numeric(temp_oral[1:11]<=oral_WOMAC_12m$est_trunc) - Q_b_R1_oral[j, ] <- as.numeric(temp_oral[12:22]<=oral_WOMAC_12m$est_trunc_R1) - Q_b_R0_oral[j, ] <- as.numeric(temp_oral[23:33]<=oral_WOMAC_12m$est_trunc_R0) + Q_b_topical[z, ] <- as.numeric(temp_topical[1:11]<=topical_WOMAC_12m$est_trunc) + Q_b_R1_topical[z, ] <- as.numeric(temp_topical[12:22]<=topical_WOMAC_12m$est_trunc_R1) + Q_b_R0_topical[z, ] <- as.numeric(temp_topical[23:33]<=topical_WOMAC_12m$est_trunc_R0) + Q_b_oral[z, ] <- as.numeric(temp_oral[1:11]<=oral_WOMAC_12m$est_trunc) + Q_b_R1_oral[z, ] <- as.numeric(temp_oral[12:22]<=oral_WOMAC_12m$est_trunc_R1) + Q_b_R0_oral[z, ] <- as.numeric(temp_oral[23:33]<=oral_WOMAC_12m$est_trunc_R0) - t_bc_topical[j, ] <- (temp_topical[1:11]-topical_vals[1:11])/sqrt(temp_topical[34:44]) - t_bc_R1_topical[j, ] <- (temp_topical[12:22]-topical_vals[12:22])/sqrt(temp_topical[45:55]) - t_bc_R0_topical[j, ] <- (temp_topical[23:33]-topical_vals[23:33])/sqrt(temp_topical[56:66]) - t_bc_oral[j, ] <- (temp_oral[1:11]-oral_vals[1:11])/sqrt(temp_topical[34:44]) - t_bc_R1_oral[j, ] <- (temp_oral[12:22]-oral_vals[12:22])/sqrt(temp_oral[45:55]) - t_bc_R0_oral[j, ] <- (temp_oral[23:33]-oral_vals[23:33])/sqrt(temp_oral[56:66]) + t_bc_topical[z, ] <- (temp_topical[1:11]-topical_vals[1:11])/sqrt(temp_topical[34:44]) + t_bc_R1_topical[z, ] <- (temp_topical[12:22]-topical_vals[12:22])/sqrt(temp_topical[45:55]) + t_bc_R0_topical[z, ] <- (temp_topical[23:33]-topical_vals[23:33])/sqrt(temp_topical[56:66]) + t_bc_oral[z, ] <- (temp_oral[1:11]-oral_vals[1:11])/sqrt(temp_oral[34:44]) + t_bc_R1_oral[z, ] <- (temp_oral[12:22]-oral_vals[12:22])/sqrt(temp_oral[45:55]) + t_bc_R0_oral[z, ] <- (temp_oral[23:33]-oral_vals[23:33])/sqrt(temp_oral[56:66]) - est_topical[j, ] <- temp_topical[1:11] - est_R1_topical[j, ] <- temp_topical[12:22] - est_R0_topical[j, ] <- temp_topical[23:33] - est_oral[j, ] <- temp_oral[1:11] - est_R1_oral[j, ] <- temp_oral[12:22] - est_R0_oral[j, ] <- temp_oral[23:33] + est_topical[z, ] <- temp_topical[1:11] + est_R1_topical[z, ] <- temp_topical[12:22] + est_R0_topical[z, ] <- temp_topical[23:33] + est_oral[z, ] <- temp_oral[1:11] + est_R1_oral[z, ] <- temp_oral[12:22] + est_R0_oral[z, ] <- temp_oral[23:33] + + num_NA_topical <- num_NA_topical+as.numeric(any(is.na(temp_topical[1:33]))) + num_NA_oral <- num_NA_oral+as.numeric(any(is.na(temp_oral[1:33]))) } @@ -169,7 +174,8 @@ return(c(topical_vals, oral_vals, r_Q_b_topical, r_Q_b_R1_topical, r_Q_b_R0_topi r_Q_b_oral, r_Q_b_R1_oral, r_Q_b_R0_oral, Q_t_b_topical, Q_t_b_R1_topical, Q_t_b_R0_topical, Q_t_b_oral, Q_t_b_R1_oral, Q_t_b_R0_oral, abs_Q_t_b_topical, abs_Q_t_b_R1_topical, abs_Q_t_b_R0_topical, abs_Q_t_b_oral, abs_Q_t_b_R1_oral, abs_Q_t_b_R0_oral, t_b_sd_topical, t_b_R1_sd_topical, t_b_R0_sd_topical, - t_b_sd_oral, t_b_R1_sd_oral, t_b_R0_sd_oral)) + t_b_sd_oral, t_b_R1_sd_oral, t_b_R0_sd_oral, num_NA_topical, num_NA_oral, + t_b_topical, t_b_R1_topical, t_b_R0_topical, t_b_oral, t_b_R1_oral, t_b_R0_oral)) } @@ -215,6 +221,16 @@ t_b_sd_oral <- Boot.Est[, 364:374] t_b_R1_sd_oral <- Boot.Est[, 375:385] t_b_R0_sd_oral <- Boot.Est[, 386:396] +num_NA_topical <- Boot.Est[, 397] +num_NA_oral <- Boot.Est[, 398] + +t_b_topical <- Boot.Est[, 399:409] +t_b_R1_topical <- Boot.Est[, 410:420] +t_b_R0_topical <- Boot.Est[, 421:431] +t_b_oral <- Boot.Est[, 432:442] +t_b_R1_oral <- Boot.Est[, 443:453] +t_b_R0_oral <- Boot.Est[, 454:464] + ## bootstrap result save(est_topical_boot, est_R1_topical_boot, est_R0_topical_boot, est_oral_boot, est_R1_oral_boot, est_R0_oral_boot, @@ -229,7 +245,7 @@ save(est_topical_boot, est_R1_topical_boot, est_R0_topical_boot, abs_Q_t_b_topical, abs_Q_t_b_R1_topical, abs_Q_t_b_R0_topical, abs_Q_t_b_oral, abs_Q_t_b_R1_oral, abs_Q_t_b_R0_oral, t_b_sd_topical, t_b_R1_sd_topical, t_b_R0_sd_topical, - t_b_sd_oral, t_b_R1_sd_oral, t_b_R0_sd_oral, + t_b_sd_oral, t_b_R1_sd_oral, t_b_R0_sd_oral, num_NA_topical, num_NA_oral, file=paste0("/uufs/chpc.utah.edu/common/home/u6070035/CCS/simResult/imputed", imputed, "/double_pboot_n", sim.size, "/", kernel, "_", single_index_method, "/sim", j, ".RData")) print(j)