R/nls_1d_fit.R

In fzwaeustc/pcrfn: Peng's Collection of R Functions at Nielsen

nls_1d_fit <-
function(mat, picked, func, 
                       para.initials=list(a1 = 0.03, b1 = 1, a2=0.3, b2=0.6)){
  
  #     mat=FEB_5KS30_ls_1_360_3D
  #     picked <- c(60, 100, 155, 215, 250, 300, 360, 367)
  #     func <- weibull
  #     adjust <- 0
  #     para.initials <- list(a = 0.03, b = 1)
  
  cases_in_data <- as.numeric(gsub("\\D", "", names(mat)))[-1]
  cat("The cases in the dataframe range from", min(cases_in_data), "to", max(cases_in_data), ".\n")
  
  if (sum(!picked%in%cases_in_data)!=0){
    cat("The following picked cases are not in the dataframe:", picked[!picked%in%cases_in_data], ".\n")
    picked <- picked[picked%in%cases_in_data]
  }
  
  
  funmatch <- match.fun(func)
  
  length_case <- mat[,c(TRUE, cases_in_data%in%picked)]
  
  names(length_case)=c("Minutes",paste("L", picked, sep="_"))
  length_case_melt <- melt(length_case, id=c("Minutes"), variable="Length",
                           value.name="Purity")
  
  #I need a vetor to store all the fitted values of the functional fitting
  fitted_value <- rep(NA, nrow(length_case_melt))
  fitted_coe <- list()
  
  problematic_length <- c()
  
  if (length(formalArgs(funmatch))==6) {
    cat("Double Weibull Survival model is fitted.\n")
    if (length(para.initials)!=4) {
      stop("4 Initial values (a1, b1, a2, b2) should be supplied for Double Weibull Function")
    }
    x_up_bound=max(picked[1:min(8, length(picked))])+1
    
    for (i in 1:length(picked)){
      fit <- try(eval(substitute(
        nls(var ~ funmatch(Minutes, a1, b1, a2, b2, duration=picked[i]), 
            data = length_case, 
            start = para.initials, 
            control=nls.control(maxiter=50, warnOnly=FALSE)), 
        list(var=as.name(names(length_case)[i+1]))
      )), TRUE)
      
      #Record all the length with has a nls fitting problem
      if (inherits(fit, "try-error"))
        problematic_length <- c(problematic_length,picked[i] )
      else {
        fitted_value[((i-1)*nrow(length_case)+1):(i*nrow(length_case))]<-
          predict(fit, list(Minutes = length_case[,1]))
        
        eval(substitute(fitted_coe$var <- summary(fit)$coe, 
                        list(var=as.name(names(length_case)[i+1]))
        ))
      }
    }
  } else if (length(formalArgs(funmatch))==3) {
    cat("Single Weibull Survival model is fitted.\n");
    if (length(para.initials)!=2) {
      stop("2 Initial values (a, b) should be supplied for single Weibull Function")
    }
    x_up_bound=nrow(mat)
    
    for (i in 1:length(picked)){
      fit <- try(eval(substitute(
        nls(var ~ funmatch(Minutes, a, b), 
            data = length_case, 
            start = para.initials, 
            control=nls.control(maxiter=50, warnOnly=FALSE)), 
        list(var=as.name(names(length_case)[i+1]))
      )), TRUE)
      
      #Record all the length with has a nls fitting problem
      if (inherits(fit, "try-error"))
        problematic_length <- c(problematic_length,picked[i] )
      else {
        fitted_value[((i-1)*nrow(length_case)+1):(i*nrow(length_case))]<-
          predict(fit, list(Minutes = length_case[,1]))
        
        eval(substitute(fitted_coe$var <- summary(fit)$coe, 
                        list(var=as.name(names(length_case)[i+1]))
        ))
      }
    }
  } else if (length(formalArgs(funmatch))==4) {
    cat("A weibull + linear model is fitted.\n");
    if (length(para.initials)!=3) {
      stop("3 Initial values (a, b, beta) should be supplied for weib.line Function")
    }
    x_up_bound=nrow(mat)
    
    for (i in 1:length(picked)){
      fit <- try(eval(substitute(
        nls(var ~ funmatch(Minutes, a, b, beta), 
            data = length_case, 
            start = para.initials, 
            control=nls.control(maxiter=50, warnOnly=FALSE)), 
        list(var=as.name(names(length_case)[i+1]))
      )), TRUE)
      
      #Record all the length with has a nls fitting problem
      if (inherits(fit, "try-error"))
        problematic_length <- c(problematic_length,picked[i] )
      else {
        fitted_value[((i-1)*nrow(length_case)+1):(i*nrow(length_case))]<-
          predict(fit, list(Minutes = length_case[,1]))
        
        eval(substitute(fitted_coe$var <- summary(fit)$coe, 
                        list(var=as.name(names(length_case)[i+1]))
        ))
      }
    }
  } else {
    stop("This function can only handle 1 dimensional non-linear models with 2 to 4 parameters")
  }
  
  if (length(problematic_length)==length(picked))
    stop("Probably all effective picked lengths have problems in fitting")
  
  length_case_melt$fitted <- fitted_value
  
  #get rid of all the records with problematic lengths
  
  if (length(picked)>8)
    cat("You picked more than 8 lengthes, thus only first 8 fittings will be shown in the nls fitting plot.\n")
  
  #This is to plot a max of the first 8 fit, or brewer is overwhelmed
  length_case_melt <- 
    length_case_melt[1:(min(8, length(picked))*nrow(length_case)),]
  
  y_low_bound=max((1-(1-min(length_case_melt$Purity, na.rm=TRUE))*1.1), 0)
  
  length_case_plot<-
    ggplot(data=length_case_melt)+
    geom_point(aes(x=Minutes, y=Purity, color=Length))+
    geom_line(aes(x=Minutes, y=fitted, color=Length), size=1)+
    coord_cartesian(xlim=c(0, x_up_bound), ylim=c(y_low_bound,1))+
    scale_colour_brewer(palette="Dark2")
  
  return(list(effective_picked=picked, fitted_coe=fitted_coe, plot=length_case_plot, 
              problematic_length=problematic_length))
}