diff --git a/HelperFunction.R b/HelperFunction.R index 73d4448..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]) } }) } @@ -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/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 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..187e028 --- /dev/null +++ b/SensIAT_sim_outcome_modeler_mave.R @@ -0,0 +1,178 @@ +# 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", + use_mave, + 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'.") + } + + 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_ + 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 + ) + ) + ) +} + diff --git a/est_s_t_y.R b/est_s_t_y.R index b6bd976..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){ + 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 @@ -97,9 +98,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 +125,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/profvis.html b/profvis.html new file mode 100644 index 0000000..cf77c96 --- /dev/null +++ b/profvis.html @@ -0,0 +1,4637 @@ + + +
+ +