R/parameter_estimation_functions.R

In sheejamk/MarkovModel: Decision Analysis Markov Model

#######################################################################
#' Get the parameter values from reading a file
#' @param parameter  parameter of interest
#' @param param_value parameter value to be assigned
#' @return the paramvalue
#' @examples
#' a <- get_parameter_direct("cost_IT", param_value=100)
#'@export
get_parameter_direct <- function(parameter, param_value){
  if(!is.null(param_value)) {
    assigned_value <- assign(parameter, param_value)
    return(assigned_value)
  }else{
    stop("Need to provide a parameter value")
  }
}
#######################################################################
#' Get the parameter values from reading a file
#' @param parameter  parameter of interest
#' @param paramfile parameter file to be provided
#' @param strategycol treatment strategy
#' @param strategyname treatment strategy name in the column strategycol
#' @return the paramvalue
#' @examples
#' @description  the file should have these column names (atleast)
#' Parameter,	Description,	Strategy,	Value
#' a <- get_parameter_read("cost_IT", paramfile=system.file("extdata","table_param.csv",
#' package = "MarkovModel"), strategycol="Strategy", strategyname="Intervention")
#' @export
get_parameter_read <- function(parameter, paramfile, strategycol=NA,strategyname = NA){
  if (is.null(paramfile))
    stop("Need to provide a parameter file to lookup")
  if (IPDFileCheck::test_file_exist_read(paramfile) != 0)
    stop("File doesnt exists or notable to access")
  dataset <- data.frame(read.csv(paramfile,header = TRUE,sep = ",", stringsAsFactors = FALSE))
  result <- IPDFileCheck::check_column_exists("value",dataset)
  if (result != 0) {
    stop(paste("Parameter file should contain value in column name",sep = ""))
  }
  if (!is.na(strategycol)) {
    if (IPDFileCheck::check_column_exists(strategycol,dataset) == 0){
      dataset = dataset[dataset[[strategycol]] == strategyname,]
      answer <- dataset[dataset$Parameter == parameter,]$Value
    }else{
      stop(paste("Parameter file should contain strategy in column name",sep = ""))
    }
  }else{
    answer <- dataset[dataset$Parameter == parameter,]$Value
    if (is.na(answer)) {
      warning("Error- Parameter value is NA-Did you use the wrong function
              to get the parameter?")
    }
  }
  return(answer)
}
#######################################################################
#' Get the definition of given parameter distribution defined in a file
#' @param parameter  parameter of interest
#' @param paramfile data file to be provided
#' @param strategycol treatment strategy column name
#' @param strategyname treatment strategy name in the column strategycol
#' @return the definition of parameter from the given distribution
#' @description  the file should have these column names (atleast)
#' Parameter,	Description,	Strategy,	Value,	Distribution,	Param1_name, Param1_value,
#' Param2_name,	Param2_value
#' @examples
#' a <- get_parameter_def_distribution("rr", paramfile=
#' system.file("extdata", "table_param.csv", package = "MarkovModel"))
#' @export
get_parameter_def_distribution <- function(parameter, paramfile, strategycol = NA,
                                           strategyname = NA){
  if (is.null(paramfile))
    stop("Need to provide a parameter file to lookup")
  dataset <- data.frame(read.csv(paramfile,header = TRUE,sep = ",",stringsAsFactors = FALSE))
  if (!is.na(strategycol)) {
    if (IPDFileCheck::check_column_exists(strategycol,dataset) != 0)
      dataset = dataset[dataset[[strategycol]] == strategyname,]
  }
  result <- IPDFileCheck::check_column_exists("Distribution",dataset)
  if (result != 0) {
    stop(paste("Parameter file should contain distribution in column name",sep = ""))
  }
  if (is.na(dataset[dataset$Parameter == parameter,]$Distribution)) {
    stop("This parameter may not be estimated from a distribution- use get_parameter_read
         instead")
  }
  param_distribution = c("Param1_name","Param1_value")
  result <- unlist(lapply(param_distribution,IPDFileCheck::check_column_exists,dataset))
  if (sum(result) != 0) {
    stop(paste("Parameter file should contain parameters for the distribtution in
               column name",sep = ""))
  }
  this_param <- dataset[dataset$Parameter == parameter,]$Parameter
  this_distr_name <- dataset[dataset$Parameter == parameter,]$Distribution
  param1 <- dataset[dataset$Parameter == parameter,]$Param1_name
  param1_value <- dataset[dataset$Parameter == parameter,]$Param1_value
  if (!is.na(dataset[dataset$Parameter == parameter,]$Param2_name)) {
    param2 <- dataset[dataset$Parameter == parameter,]$Param2_name
    param2_value <- dataset[dataset$Parameter == parameter,]$Param2_value
    expression_created = paste(this_distr_name,"(",param1 ,"=", param1_value, "," ,
                               param2, "=",param2_value,")",sep = "")
  }else{
    expression_created = paste(this_distr_name,"(",param1, "=" , param1_value,")",sep = "")
  }
  expression_recreated <- check_estimate_substitute_proper_params(expression_created)
  param_with_expression <- paste(this_param, " = ", expression_recreated, sep = "")
  return(param_with_expression)
}

#######################################################################
#' Get the parameter values using the provided statistical regression methods
#' @param param_to_be_estimated  parameter of interest
#' @param dataset data set to be provided
#' @param method methd of estimation (for example, linear, logistic regression etc)
#' @param indep_var the independent variable (column name in data file)
#' @param info_get_method additional information on methods e.g Kaplan-Meier ot hazard
#' @param info_distribution distribution name  eg. for logistic regression -binomial
#' @param covariates list of covariates - calculations to be done before passing
#' @param strategycol column name containing arm details, default is NA
#' @param strategyname name of the  arm, default is NA
#' @param timevar_survival time variable for survival analysis, default is NA
#' @return the results of the regression analysis
#' @examples
#' mydata <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")
#' results_logit <- get_parameter_estimated_regression("admit", mydata,
#' "logistic regression", "gre", NA,"binomial", c("gpa", "factor(rank)"))
#' @export
get_parameter_estimated_regression <- function(param_to_be_estimated, dataset, method,
                                            indep_var, info_get_method, info_distribution,
                                            covariates = NA, strategycol = NA,
                                            strategyname = NA, timevar_survival = NA){
  if (is.null(dataset))
    stop("Need to provide a data set to lookup")
  if (is.na(method))
    stop("Please provide  a statistical method")
  if (!is.na(strategycol)) {
    if (IPDFileCheck::check_column_exists(strategycol,dataset) != 0)
      dataset = dataset[dataset[[strategycol]] == strategyname,]
  }
  method_rpack <- find_keyword_regression_method(method,info_get_method)
  caps_info_method <- toupper(info_get_method)
  caps_method = toupper(method)
  if (caps_method == "SURVIVAL" | caps_method == "SURVIVAL ANALYSIS") {
    if (is.na(timevar_survival)) {
      stop("For survival analysis, please provide the varaible to use as time ")
    }else{
      surv_object <- paste("survival::Surv(", timevar_survival, ",", param_to_be_estimated, ")", sep = "")
      object <- surv_object
    }
  }else{
    object <- param_to_be_estimated
  }
  if (sum(is.na(covariates)) == 0) {
    covariates_list = list()
    for (i in 1:length(covariates)) {
      if (i == length(covariates))
        covariates_list = paste(covariates_list,paste(covariates[i], sep = ""))
      else
        covariates_list = paste(covariates_list,paste(covariates[i], " + ", sep = ""))
    }
  }else{
    covariates_list = NA
  }
  if (caps_method == "SURVIVAL" | caps_method == "SURVIVAL ANALYSIS") {
    if (caps_info_method == "PARAMETRIC REGRESSION" | caps_info_method == "PARAMETRIC") {
      this_dist <- find_survreg_distribution(info_distribution)
      if (is.na(covariates_list))
        expression_recreated <- paste0(method_rpack,"(", object, " ~ ", indep_var, ", ",
                                     "data = dataset,  dist = \"", this_dist, "\" ) ", sep = "")
      else
        expression_recreated <- paste0(method_rpack,"(", object, " ~ ", indep_var, " + ",
                                     covariates_list, ", ","data = dataset,  dist = \"",
                                     this_dist, "\" ) ", sep = "")
    }else{
      if (caps_info_method == "KAPLAN-MEIER" | caps_info_method == "KM") {
        if (is.na(covariates_list))
              expression_recreated <- paste0(method_rpack,"(", object, " ~ ", indep_var,",type =",
                                           "\"kaplan-meier\"", ", ","data = dataset) ", sep = "")
        else
            expression_recreated <- paste0(method_rpack,"(", object, " ~ ", indep_var, " + ",covariates_list,
                                         ", type =", "\"kaplan-meier\"",", ","data = dataset) ", sep = "")
      }else{
        if (caps_info_method == "FLEMING-HARRINGTON" | caps_info_method == "FH") {
          if (is.na(covariates_list))
                expression_recreated <- paste0(method_rpack,"(", object, " ~ ", indep_var,", type =",
                                             "\"fleming-harrington\"",", ","data = dataset) ", sep = "")
          else
                expression_recreated <- paste0(method_rpack,"(", object, " ~ ", indep_var, " + ",covariates_list,
                                             ", type =", "\"fleming-harrington\"",", ","data = dataset) ", sep = "")
        }else{
          if (caps_info_method == "FH2") {
            if (is.na(covariates_list))
              expression_recreated <- paste0(method_rpack,"(", object, " ~ ", indep_var, ", type =", "\"fh2\"",
                                           ", ","data = dataset) ", sep = "")
            else
              expression_recreated <- paste0(method_rpack,"(", object, " ~ ", indep_var, " + ", covariates_list,
                                           ", type =", "\"fh2\"",", ","data = dataset) ", sep = "")
          }else{
            if (is.na(covariates_list))
              expression_recreated <- paste0(method_rpack,"(", object, " ~ ", indep_var, ", ","data = dataset) ", sep = "")
            else
              expression_recreated <- paste0(method_rpack,"(", object, " ~ ", indep_var, " + ", covariates_list,
                                           ", ", "data = dataset) ", sep = "")
          }
        }
      }
    }
  }else{
    if (is.na(covariates_list))
      expression_recreated <- paste0(method_rpack,"(", object, " ~ ", indep_var, ", ","data = dataset) ", sep = "")
    else
      expression_recreated <- paste0(method_rpack,"(", object, " ~ ", indep_var, " + ", covariates_list,
                                   ", ", "data = dataset) ", sep = "")
  }

  param_estimated <- eval(parse(text = expression_recreated))
  summary_regression_results = summary(param_estimated)
  if (caps_info_method %in% c("KAPLAN-MEIER","KM","FLEMING-HARRINGTON", "FH2","FH"))
    vcov_param_estimated <- chol_decomp_matrix <- list()
  else{
    vcov_param_estimated <- stats::vcov(param_estimated)
    chol_decomp_matrix <- chol(vcov_param_estimated)
  }
  if (caps_info_method %in% c("PARAMETRIC REGRESSION","PARAMETRIC"))
    plot_result <- graphics::plot(param_estimated)
  else{
    if (caps_info_method %in% c("COX-PROPORTIONAL-HAZARD","COX PROPORTIONAL HAZARD","COX-PH", "COX PH","COXPH"))
      plot_result <- ggplot2::autoplot(survival::survfit(param_estimated))
    else
    plot_result <- ggplot2::autoplot(param_estimated)
  }

  results =  structure(list(
    param_estimated = param_estimated,
    variance_covariance = vcov_param_estimated,
    summary_regression_results = summary_regression_results,
    cholesky_decomp_matrix = chol_decomp_matrix,
    plot = plot_result
  ))
  # if(caps_method == "SURVIVAL" | caps_method == "SURVIVAL ANALYSIS"){
  #   if(!is.na(info_distribution) & trimws(toupper(info_distribution))== "WEIBULL"){
  #     params_hr_etr = SurvRegCensCov::ConvertWeibull(param_estimated)
  #     results =  structure(list(
  #       param_estimated = param_estimated,
  #       variance_covariance = vcov_param_estimated,
  #       summary_regression_results = summary_regression_results,
  #       cholesky_decomp_matrix = chol_decomp_matrix,
  #       weibull_params_survival_analysis = params_hr_etr,
  #       plot = plot_result
  #     ))
  #   }
  #}
  return(results)
}
#' ##########################################################################################################
#' #' Function for mixed effect regression
#' #' @param data a dataframe
#' #' @param yy column name of dependent variable
#' #' @param fixed_effect column name of independent variable
#' #' @param random_effects, independent variables , an empty vector by default
#' #' @return result regression result with plot if success and -1, if failure
#' #' @examples performMixedEffect(data,"ShOWSscore","does_redn",c())
#' #' @export
#' performMixedEffect<-function(this_data, yy, fixed_effect,random_effects=c()){
#'   if(length(random_effects)!=0){
#'     expre=paste("1|",random_effects[1],sep="")
#'     i=2
#'     while(i<=length(random_effects)){
#'       this=paste("1|",random_effects[i],sep="")
#'       expre=paste(expre,this,sep="+")
#'       i=i+1
#'     }
#'     fmla <- as.formula(paste(yy, paste("~"),paste(fixed_effect,"+"),paste(expre, collapse= "+")))
#'     fit<-lmer(fmla,data=this_data)
#'     p_plt<-ggPredict(fit,se=TRUE,interactive=TRUE)
#'   }else{
#'     fmla <- as.formula(paste(yy, paste("~"),paste(fixed_effect,collapse= "+")))
#'     fit<-lm(fmla,data=this_data)
#'     p_plt<-ggPredict(fit,se=TRUE,interactive=TRUE)
#'   }
#'   return(list(summary(fit),p_plt))
#' }
#'
#' ##########################################################################################################
#' #' Function for linear regression
#' #' @param data a dataframe
#' #' @param yy column name of dependent variable
#' #' @param x1 column name of independent variable 1
#' #' @param extra_variables, independent variables , an empty vector by default
#' #' @param interaction, boolean value to indicate the interaction,  default is FALSE
#' #' @return result regression result with plot if success and -1, if failure
#' #' @examples performLinearReg(data,"EQ5D3L_from5L","does_redn",c(), FALSE)
#' performLinearReg<-function(this_data, yy,x1,extra_variables=c(),interaction=FALSE){
#'   if(length(extra_variables)==0){
#'     # no need to check for interaction
#'     fmla <- as.formula(paste(yy, paste("~"),paste(x1,collapse= "+")))
#'     fit<lm(fmla,data=this_data)
#'     p_plt<-ggPredict(fit,se=TRUE,interactive=TRUE)
#'   }else{
#'       expre=paste(extra_variables[1],sep="")
#'       i=2
#'       while(i<=length(extra_variables)){
#'         this=paste(extra_variables[i],sep="")
#'         expre=paste(expre,this,sep="+")
#'         i=i+1
#'       }
#'       if(interaction==FALSE){
#'         fmla <- as.formula(paste(yy, paste("~"),paste(x1,"+",sep=""),paste(expre, collapse= "+")))
#'         fit<-lm(fmla,data=this_data)
#'         p_plt<-ggPredict(fit,se=TRUE,interactive=TRUE)
#'       }else{
#'         expre=paste(extra_variables[1],sep="")
#'         i=2
#'         while(i<=length(extra_variables)){
#'           this=paste(extra_variables[i],sep="")
#'           expre=paste(expre,this,sep="*")
#'           i=i+1
#'         }
#'         fmla <- as_formula(paste(yy, paste("~"),paste(x1,"*",sep=""),paste(expre, collapse= "*")))
#'         fit<-lm(fmla,data=this_data)
#'         p_plt<-ggPredict(fit,se=TRUE,interactive=TRUE)
#'       }
#'   }
#'   results<-list(summary(fit),p_plt)
#'   return(results)
#' }
#'