From 6fbd170f3b0970cc04e42771ef343c178d2c336e Mon Sep 17 00:00:00 2001 From: awamaeva Date: Mon, 29 Jan 2024 16:32:04 -0500 Subject: [PATCH 1/4] Rename trajHRMSM_gform.R to trajhrmsm_gform.R --- R/{trajHRMSM_gform.R => trajhrmsm_gform.R} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename R/{trajHRMSM_gform.R => trajhrmsm_gform.R} (100%) diff --git a/R/trajHRMSM_gform.R b/R/trajhrmsm_gform.R similarity index 100% rename from R/trajHRMSM_gform.R rename to R/trajhrmsm_gform.R From 2b57664e62f78c2ac6b0ad7b652568a238a798db Mon Sep 17 00:00:00 2001 From: awamaeva Date: Mon, 29 Jan 2024 16:32:23 -0500 Subject: [PATCH 2/4] Rename trajHRMSM_IPW.R to trajhrmsm_ipw.R --- R/{trajHRMSM_IPW.R => trajhrmsm_ipw.R} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename R/{trajHRMSM_IPW.R => trajhrmsm_ipw.R} (100%) diff --git a/R/trajHRMSM_IPW.R b/R/trajhrmsm_ipw.R similarity index 100% rename from R/trajHRMSM_IPW.R rename to R/trajhrmsm_ipw.R From e941463e6eac9d8f82ea03ed6716cd27e9d0f086 Mon Sep 17 00:00:00 2001 From: awamaeva Date: Mon, 29 Jan 2024 16:32:38 -0500 Subject: [PATCH 3/4] Rename trajHRMSM_pltmle.R to trajhrmsm_pltmle.R --- R/{trajHRMSM_pltmle.R => trajhrmsm_pltmle.R} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename R/{trajHRMSM_pltmle.R => trajhrmsm_pltmle.R} (100%) diff --git a/R/trajHRMSM_pltmle.R b/R/trajhrmsm_pltmle.R similarity index 100% rename from R/trajHRMSM_pltmle.R rename to R/trajhrmsm_pltmle.R From 430fd178ab18bf90c4d17c4ff58257cd98dced87 Mon Sep 17 00:00:00 2001 From: awamaeva Date: Mon, 29 Jan 2024 16:32:56 -0500 Subject: [PATCH 4/4] Delete R/pltmle_countermeans2.R --- R/pltmle_countermeans2.R | 153 --------------------------------------- 1 file changed, 153 deletions(-) delete mode 100644 R/pltmle_countermeans2.R diff --git a/R/pltmle_countermeans2.R b/R/pltmle_countermeans2.R deleted file mode 100644 index 7ef9434..0000000 --- a/R/pltmle_countermeans2.R +++ /dev/null @@ -1,153 +0,0 @@ -#' @title Counterfactual means for a Pooled LTMLE -#' @description function to estimate counterfactual means for a pooled LTMLE. -#' @name sub_pltmle2 -#' @param formula specification of the model for the outcome to be fitted. -#' @param id name of the column for unique identifiant. -#' @param V baseline covariates. -#' @param L time-varying covariates. -#' @param A time-varying treatment. -#' @param Y outcome variable. -#' @param s number of measuring times. -#' @param J number of trajectory groups. -#' @param time measuring times. -#' @param obsdata observed data in wide format. -#' @returns \item{list_pltmle_countermeans}{Counterfactual means and influence functions with the pooled ltmle.} -#' \item{D}{Influence functions} -#' @import e1071 -#' @examples -#' obsdata_msm <- genobsdata_msm(n=2500) -#' V <- "V" -#' L <- "L" -#' time <- 1:3 -#' A <- "A" -#' Y <- "Y" -#' @author Awa Diop, Denis Talbot - - -sub_pltmle2<-function(formula,Y,A,L,V,s,time,timevar,idvar,obsdata,Traj){ - A. = paste0(A,".") - L. = paste0(L,".") - Y. = paste0(Y,".") - varA <- paste0(A.,time[1:s]) - D = NULL - j = 0; - obsdata0.all=list() - obsdata.all=list() - nb_regimes = 2^s #number of treatment regimes - logit = qlogis; - expit = plogis; - - if(nb_regimes > 1000) { - cat("With that number of periods, there is more than 1,000 thousand potential treatment regimes:",nb_regimes); - cont <- readline("Continue? y/n: "); - if(cont == "n" | cont == "no" | cont == "NO") stop("Aborted.") - } - #Treatment regimes - dat.combn = bincombinations(s); - list.regimes=lapply(1:nb_regimes,function(x){; - regime = dat.combn[x,]; - return(regime); - }) - - j = 0 - for(regimes in list.regimes){ - j=j+1 - obsdata0 = obsdata; - obsdata0[,paste0(A.,time[1:s])] <- do.call(cbind,lapply(1:s, function(x){ - rgm <- rep(regimes[x],nrow(obsdata0)) - return(rgm)})); - - obsdata0.all[j]=list(obsdata0) - obsdata.all[j]=list(obsdata) - } - obsdata.all2 = data.frame(do.call(rbind,obsdata.all)) - obsdata0.all2 = data.frame(do.call(rbind,obsdata0.all)) - - #Model - modQs = glm(formula = formula, family = binomial, data = obsdata); - #for t=3 predict the outcome for all different regimes of treatment - Qs = lapply(1:nb_regimes,function(i)predict(modQs, newdata = obsdata0.all[[i]], type = "res")); - - #Compute the weights for all different regimes of treatment - W = IPW(numerator = "unstabilized", id = idvar , V = V,L =L.,A = A., C = FALSE, - censor = censor, K = s,obsdata = obsdata); - - Hs.all=list() - j=0; - for(regimes in list.regimes){ - j=j+1 - Hs.all[j] = list(as.matrix(1*(rowSums(sapply(1:s,function(x) - 1*(obsdata[,varA[x]] == unlist(regimes)[x])))==s)*W[,s])%*%t(as.matrix(Traj[j,1:J]))) - } - - Hs=do.call(rbind,Hs.all) - # Update the risk for each regime of treatment - ## First estimate the vector of fluctuation error - - modEs = glm(as.formula(paste0(paste0(Y.,s), "~", "-1 + offset(logit(unlist(Qs))) + Hs")), family = binomial, data = obsdata.all2); - coef_Es <- coef(modEs) - coef_Es<- ifelse(is.na(coef_Es),0, coef_Es) - Qstar = lapply(1:nb_regimes,function(x) expit(logit(Qs[[x]]) + - as.numeric(t(t(as.matrix(coef_Es))%*%t(Hs.all[[x]]))))); - - ##Influence curve for each regime of treatment - D_list <- list() - Ds = lapply(1:nb_regimes,function(x)sapply(1:J,function(y) - (Hs.all[[x]][,y]*(obsdata[,paste0(Y.,s)] - Qstar[[x]])))) - - D_list[s]<-list(Ds) - - for(i in (s-1):1){ - #For t=s-1 repeat the same process as for t=s - sf1 <- paste0(A.,time[1:i], collapse = "+"); - sf2 <- paste0(unlist(lapply(1:length(L.), function(x) paste0(L.[x],time[1:i], collapse = "+"))),collapse = "+"); - sf3 <- paste0(V,collapse = "+") - - - modQs = lapply(1:nb_regimes,function(x){ - qstar = Qstar[[x]]; - form.up <- as.formula(paste("qstar ~",paste0(sf1, "+", sf2, "+", sf3))); - mod_temp = glm(form.up, family = binomial, data = obsdata); - return(mod_temp); - }); - - Qs = lapply(1:nb_regimes,function(x)predict(modQs[[x]], newdata = obsdata0.all[[x]], type = "res")); - - Hs.all=list() - j=0 - for(regimes in list.regimes){ - j=j+1 - Hs.all[j] = list(as.matrix(1*(rowSums(sapply(1:i,function(x) 1*(obsdata[,varA[x]] == unlist(regimes)[x])))==i)*W[,i])%*%t(as.matrix(Traj[j,1:J]))) - } - - Hs = do.call(rbind,Hs.all) - - modEs = glm(unlist(Qstar) ~-1+offset(logit(unlist(Qs)))+Hs, family = binomial,data = obsdata.all2); - coef_Es <- coef(modEs) - coef_Es<- ifelse(is.na(coef_Es),0, coef_Es) - Qstarm1 = lapply(1:nb_regimes,function(x)expit(logit(Qs[[x]]) + - as.numeric(t(t(as.matrix( coef_Es))%*%t(Hs.all[[x]]))))); - - ##Influence curve for each regime of treatment - Ds = lapply(1:nb_regimes,function(x) sapply(1:J,function(y) - as.matrix(Hs.all[[x]][,y]*(Qstar[[x]]-Qstarm1[[x]])))) - D_list[i]<-list(Ds) - Qstar = Qstarm1; - } - - Q0 = unlist(Qstar); - - D = lapply(1:nb_regimes,function(x) as.matrix(Reduce('+',lapply(1:J, function(y){ - comps = as.matrix(D_list[[y]][[x]] + as.numeric(Qstar[[x]] - mean(Qstar[[x]]))) - return(comps) - })))) - - obsdata0.all2$Y = Q0 - obsdata0.all2$atraj <- apply(obsdata0.all2[,varA],1,paste0,collapse="") - obsdataT2 = aggregate(as.formula(paste0(Y,"~", "atraj")), FUN = mean,data = obsdata0.all2) - obsdataT2[,Y] = as.numeric(as.character(obsdataT2[,Y])) - - - list_pltmle_countermeans = list(counter.means = obsdataT2[,Y],D=D) - return(list_pltmle_countermeans) -}