R/Fit_AFS.R

In meganmichelle/RAREsim_package: Simulation of Rare Variant Genetic Data

#' RAREsim
#'
#' Simulate rare variant genetic data
#'
#' @param prop_df data frame with 3 columns, Lower, Upper (of MAC bin)  and proportion of variants in that MAC bin
#'
#' @param N Number of individuals in the target data
#'
#' @param p_rv percent of rare variants - just the sum of the proportions
#'
#' @return Vector of parameters - alpha, beta, and b as well as  fitted values
#'
#' @author Megan M Null, \email{megan.null@ucdenver.edu}
#' 
#' @keywords RAREsim
#'
#'
#' @export
#' @importFrom nloptr slsqp
#'

Fit_afs <- function(prop_df, N, p_rv = NULL){ ### only works with N>2200

  
  prop_df$prop <- (as.numeric(as.character(prop_df$prop)))
  
  if(anyNA(prop_df$prop) == TRUE){
    stop('proportions need to be numeric with no NA values')
  }
  
  if(is.numeric(N)==FALSE){
    stop('N needs to be a number')
  }
  
  if(is.null(p_rv)){
    p_rv <- sum(prop_df$prop)
  }
  
  hin.tune <- function(x) {
    h <- numeric(1)
    h[1] <- x[1] 
    return(h)
  }

  c1 <- 1:prop_df$Upper[nrow(prop_df)] #### MACs to use in the function
  
  #suppressMessages()
  calc_prob_LS<-function(tune){
    ### calculate b
    alpha_mac <- 1/((c1+tune[2])^(tune[1]))
    b <- p_rv/sum(alpha_mac) #### now we have b!
    b_mac <- b*alpha_mac
    all <- 0
    for(i in 1:nrow(prop_df)){
      E <- sum(b_mac[prop_df$Lower[i]:prop_df$Upper[i]])
      O <- prop_df$prop[i]
      c <- (E-O)^2
      all <- all+c
    }
    return(all)
  }  
  
  tune <- c(1,0) #### start with the function 1/x
  S <- suppressMessages(slsqp(tune, fn = calc_prob_LS, hin = hin.tune ))
  b <- p_rv/sum(1/((c1+S$par[2])^(S$par[1])))
  
  ### now get the proportion for each bin?

  re <- prop_df
  re$prop <- '.'
  colnames(re)[3] <- 'Prop'
  
  fit <- as.data.frame(matrix(nrow =  1,  ncol = prop_df$Upper[nrow(prop_df)]))
  
  for(i in 1:prop_df$Upper[nrow(prop_df)]){
    fit[,i] <- b/((S$par[2]+i)^S$par[1])
  }
  
  for(i in 1:nrow(re)){
  re$Prop[i] <- sum(fit[,c(re$Lower[i]:re$Upper[i])])
  }
  
  re$Prop <- as.numeric(as.character(re$Prop))
  
  return(list(alpha = S$par[1], beta = S$par[2], b = b, Fitted_results = re))
  
}