R/CL_estimation.R

In EHRmuse: Multi-Cohort Selection Bias Correction using IPW and AIPW Methods

CL_wt_multi_cohort <- function(K,intdata_list,select_var_list,marginals_list){
  # expit<-function(x){
  #   return(exp(x)/(1+exp(x)))
  # }

  
  # Initialize lists to store results for each cohort
  estweights_list <- list()
  gamma_estimates <- list()
  
  # Define expit function if not already defined
  # expit <- function(x) 1 / (1 + exp(-x))
  
  for (k in 1:K){
    intdata=intdata_list[[k]]
    
    # Check select_k_var and select_k_e to be data frames or not
    if (!is.data.frame(intdata)) stop("intdata must be a data frame.")
    
    intdata_select=intdata%>%
      dplyr::select(select_var_list[[k]])
    
    
    ##marginals for cohort K
    marg_k=marginals_list[[k]]
    
    ## marg_k should be a numeric vector
    if (!is.numeric(marg_k)) stop("marg_k must be a numeric vector.")
    
    
    prop<-function(gamma){
      y <- c(rep(0,(ncol(intdata_select)+1)))
      for(i in 1:nrow(intdata_select)){
        vec=c(1, as.numeric(intdata_select[i,]))
        y = y + 1/(as.vector(expit(gamma %*% vec))) * vec 
      }
      y= y - marg_k
      y
    }
    
    # Set starting values for gamma based on dimensions
    start <- c(rep(0,(ncol(intdata_select)+1)))
    
    # Optimization
    z <- nleqslv(x = start, fn = prop, method = "Newton", global = "dbldog",
                 control = list(trace = 1, allowSingular = TRUE))
    
    # gamma_values_for_cohort_k
    gamma_estimates[[k]] <- z$x
  }
  comb_est_weights_list=NULL
  for(i in 1:K){
    select_i_i=intdata_list[[i]]%>%
      dplyr::select(select_var_list[[i]])
    
    estweights_list[[i]]=matrix(0,nrow(select_i_i),K)
    for(j in 1:K){
      select_i_j=intdata_list[[i]]%>%
        dplyr::select(select_var_list[[j]])
      estweights_list[[i]][,j]= 1/expit(as.matrix(cbind(1,select_i_j)) %*% 
                                          as.numeric(gamma_estimates[[j]]))
    }
    comb_est_weights=rep(1,nrow(select_i_i))
    for(l in 1:K){
      comb_est_weights=comb_est_weights*(1-1/estweights_list[[i]][,l])
    }
    comb_est_weights_list[[i]]=1/(1-comb_est_weights)
  }  
  
  return(list(gamma_estimates = gamma_estimates, 
              combined_weights = comb_est_weights_list))
}

# source("weighted.R")

CL_est_multi_cohort <- function(K,intdata_list,select_var_list,
                                marginals_list,Z_names){
  wts_fun=CL_wt_multi_cohort(K,intdata_list,select_var_list,marginals_list)
  wts=wts_fun$combined_weights
  intdata_comb=NULL
  comb_est_weights=NULL
  for(i in 1:K){
    intdata=intdata_list[[i]]
    if (!is.data.frame(intdata)) stop("All Internal data should be a data frame.")
    intdata_comb=rbind.data.frame(intdata_comb,intdata)
    comb_est_weights=c(comb_est_weights,wts[[i]])
  }
  intdata_comb1=intdata_comb[!duplicated(intdata_comb$id),]
  comb_est_weights1=comb_est_weights[!duplicated(intdata_comb$id)]

  fin_est=weighted(intdata=intdata_comb1,
                   estweights=comb_est_weights1,
                   Z_names=Z_names)
  
  return(list(gamma_estimates = wts_fun$gamma_estimates, 
              combined_weights = wts,
              final_est=fin_est$final,
              variance_est=fin_est$var))
}