Hello guys, I need assistance in applying makeCluster to run this program 1000 times and getting 1000 param_final. Note each time the program runs, a different param_final is obtained since we generate a uniform random variable at the start of the program. Many thanks
m_ps <- vector(mode = "list", length = 10)
prs_df <- vector(mode = "list", length = 10)
#for (sim in 1:10){
kyoto_1<-kyoto_f %>%
rowwise() %>%
mutate(
dumy=1,
USUBJID=UNIQUE_PATIENT_ID,
SURV100=ifelse(survtim <100,"N","Y"),
a=runif(1)#### each iteration will generate a different random variable
)
kyoto_1 <- kyoto_1 %>% filter (nrisk>1)
kyoto_1<-kyoto_1[order(kyoto_1$a,decreasing=FALSE),]
kyoto_1.1<-merge(iwant,kyoto_1,by.x="dumy",by.y="dumy")
kyoto_1.1<-kyoto_1.1 %>% mutate(b=1:nrow(kyoto_1.1))
kyoto_1.1$diag2fd <- NA
for (i in 2:21) {
kyoto_1.1$diag2fd <- ifelse(kyoto_1.1$b > kyoto_1.1[[paste0("perc", i-1)]], kyoto_1.1[[paste0("day", i)]], kyoto_1.1$diag2fd)
}
kyoto_1.1$diag2fd <- ifelse(kyoto_1.1$b <= kyoto_1.1$perc1, kyoto_1.1$day1, kyoto_1.1$diag2fd )
kyoto_1.2<-kyoto_1.1 %>% mutate(
nsurvtim=survtim-abs(diag2fd),
nsurv100=ifelse(nsurvtim<100,"N","Y"),
trt = "K"
)
kyoto_1.2 <- kyoto_1.2 %>% filter (nsurvtim>0)
kyoto_1.3=subset(kyoto_1.2 , select = c(USUBJID,dthfl,pulmdys,renadyf,gvhd, gidys,neurdys, sex,sinfect,SURV100,survtim,ORGDYS,AGECAT,LDHG2,nrisk,nsurv100,diag2fd,nsurvtim,trt,dumy))
tma<- tma %>% mutate(trt = "N",dumy=1)
tma.1<- subset(tma, select=c(USUBJID,DTHFL,GIDYS,GVHD,NEURDYS, PULMDYS, RENADYF, SEX, SINFECT, SURV100, SURVTIM,
ORGDYS,AGECAT,LDHG2,nrisk, diag2fd, nsurvtim, nsurv100, trt,dumy))
names(kyoto_1.3)[names(kyoto_1.3)=="dthfl"] <- "DTHFL"
names(kyoto_1.3)[names(kyoto_1.3)=="pulmdys"] <- "PULMDYS"
names(kyoto_1.3)[names(kyoto_1.3)=="renadyf"] <- "RENADYF"
names(kyoto_1.3)[names(kyoto_1.3)=="sex"] <- "SEX"
names(kyoto_1.3)[names(kyoto_1.3)=="sinfect"] <- "SINFECT"
names(kyoto_1.3)[names(kyoto_1.3)=="survtim"] <- "SURVTIM"
names(kyoto_1.3)[names(kyoto_1.3)=="surv100"] <- "SURV100"
names(kyoto_1.3)[names(kyoto_1.3)=="gvhd"] <- "GVHD"
names(kyoto_1.3)[names(kyoto_1.3)=="gidys"] <- "GIDYS"
names(kyoto_1.3)[names(kyoto_1.3)=="neurdys"] <- "NEURDYS"
both<-rbind(tma.1,kyoto_1.3)
both.1<-both %>% mutate(cnsr=ifelse(DTHFL == "Y",1,0), exp=ifelse(trt == "N",1,0),nsv100=ifelse(nsurv100 == "Y",1,0))
#outcome<-both.1$y
x<- subset(both.1, select=c( SEX, AGECAT,GVHD,SINFECT,NEURDYS, PULMDYS,RENADYF,LDHG2,trt,nsurvtim,cnsr,exp,nsurv100,nsv100,dumy))
mean(both.1$nsurvtim)
#####################################################################
############## PS MATCHING ##########################################
#####################################################################
max_comb_len <- 254
Combinations <- expand.grid(SEX=0:1,AGECAT=0:1,GVHD=0:1,SINFECT=0:1,NEURDYS=0:1,PULMDYS=0:1,RENADYF=0:1,LDHG2=0:1)
Combinations <- sapply(Combinations, as.logical)
Combinations <- Combinations[2:(max_comb_len+1),] #drop the first row that is all FALSE
DF_List <- vector(mode = "list", length = max_comb_len) #Make list to store results
lm_list <- vector(mode = "list", length = max_comb_len) #Make list to store results
BIC_list <- vector(mode = "list", length = max_comb_len) #Make list to store results
#min_var <- vector(mode = "list", length = max_comb_len) #Make list to store results with final variables
(DF <- tibble(x) |> mutate(x1= exp,time=nsurvtim,cens=cnsr))
min_bic <- min_bic_pos <- Inf
for(i in 1:max_comb_len) {
vars_in_subset <- names(which(Combinations[i,]))
#print(vars_in_subset)
DF_List[[i]] <- DF[, c("x1", vars_in_subset)]
lm_list[[i]] <- glm(x1 ~ ., family=binomial(), data=DF_List[[i]])
BIC_list[[i]] <- BIC(lm_list[[i]])
if(BIC_list[[i]] < min_bic){
min_bic <- BIC_list[[i]]
min_bic_pos <- i
}
}
min_bic
min_bic_pos
min_var <- as.data.frame(colnames(DF_List[[min_bic_pos]]))
#min_var <- min_var[-1,]
min_var_1 <- transpose(min_var)
m_ps[[min_bic_pos]]<-glm(x1 ~ ., family=binomial(),data=DF_List[[min_bic_pos]])
summary(m_ps[[min_bic_pos]])
prs_df[[min_bic_pos]]<-data.frame(pr_score=predict(m_ps[[min_bic_pos]],type="response"),x1=m_ps[[min_bic_pos]]$model$x1)
head(prs_df[[min_bic_pos]])
DF_List_all<-cbind(x,prs_df[[min_bic_pos]])
########Quantiles computation from PScore
pct<-quantile(DF_List_all$pr_score, probs =c(0,0.2,0.4,0.6,0.8,1) , na.rm = T)
pct<-as.data.frame(pct)
pct_t<-transpose(pct)
pct_t<-pct_t %>%
rowwise() %>%
mutate(
dumy=1
)
DF_List_all_1<-merge(pct_t,DF_List_all,by.x="dumy",by.y="dumy")
DF_List_all_1$pscat <- ifelse(DF_List_all_1$pr_score <= DF_List_all_1$V2, 1,
ifelse(DF_List_all_1$pr_score <= DF_List_all_1$V3, 2,
ifelse(DF_List_all_1$pr_score <= DF_List_all_1$V4, 3,
ifelse(DF_List_all_1$pr_score <= DF_List_all_1$V5, 4, 5))))
Cox Model
mod.allison <- coxph(Surv(nsurvtim, cnsr) ~ trt + pscat,
data=DF_List_all_1)
summary(mod.allison)
param<-as.data.frame(cbind(mod.allison$coef ,sqrt(mod.allison$var[1,1])))
param<-param %>% slice(1:1)
param<-param %>%
rename(
Estimate = V1,
SE = V2
)
param_final<-cbind(param,min_var_1)