我如何在r中使用catch for NLS功能

Question

我正在为四线性函数进行回归。我有两个选择是使用nlslm和nls2。但是，对于某些数据集，使用nlslm将某些问题置于：在初始参数估计的单数梯度矩阵 ，或者它们跑到Infinitie Loop中。我想使用尝试捕获来解决这个问题。谁能帮我吗？提前感谢大家。

这是完整的代码：

 # Packages needed for estimaton of Ideal trajectory - nonlinear regression
  #-------------------------------------------------------------------------------
  library("minpack.lm")
  library("nlstools")
  library("nlsMicrobio")
  library("stats") 
  library("tseries") #runs test for auto correlation
  #Use NLS2
  library(proto)
  library(nls2)
  
  ################################################################

  # Set working directory
  setwd("C:/Users/Kevin Le/PycharmProjects/Pig Data Black Box - Copy")
  
  #load dataset
  load("Data/JRPData_TTC.Rdata") #load dataset created in MissingData.step
  ID <- 5470
  #Create a new dataframe which will store Data after ITC estimation

  #Dataframe contains ITC parameters
  ITC.param.pos2 <- data.frame(ANIMAL_ID=factor(),
                               X0=double(),
                               Y1=double(),
                               Y2=double(),
                               Ylast=double(),
                               a=double(),
                               b=double(),
                               c=double(),
                               d=double(),
                               stringsAsFactors=FALSE)

  #Dataframe contains data points on the ITC
  Data.remain <- data.frame(ANIMAL_ID=character(),
                            Age=double(),
                            obs.CFI=double(),
                            tt=double(),
                            ttt=double(),
                            stringsAsFactors=FALSE)

  #===============================================================
  # For loop for automatically estimating ITC of all pigs
  #===============================================================
  IDC <-  seq_along(ID)  # 17, 23, 52, 57, 116
  for (idc in IDC){

    # idc = 1
    i <- ID[idc]
    Data <- No.NA.Data.1[No.NA.Data.1$ANIMAL_ID == i,]
    idc1 <- unique(as.numeric(Data$idc.1))
    ####### Create data frame of x (Age) and y (CFI) ########
    x <- as.numeric(Data$Age.plot)
    Y <- as.numeric(Data$CFI.plot)
    Z <- as.numeric(Data$DFI.plot)
    Data.xy <- as.data.frame(cbind(x,Y))

    #Initial parameteres for parameter estimation
    X0.0 <- x[1]
    Xlast <- x[length(x)]

    ##################################################################
    # 1. reparametrization CFI at X0 = 0
    #function used for reparametrization in MAPLE
    # solve({
    # 0=a+b*X_0+c*X_0**2,
    # DFIs=b+2*c*Xs,CFIs=a+b*Xs+c*Xs**2},
    # {a,b,c});
    # a = -X0*(2*CFIs*Xs-CFIs*X0-Xs^2*DFIs+Xs*DFIs*X0)/(Xs^2-2*X0*Xs+X0^2)
    # b = (-Xs^2*DFIs+DFIs*X0^2+2*CFIs*Xs)/(Xs^2-2*X0*Xs+X0^2)
    # c = -(CFIs-Xs*DFIs+X0*DFIs)/(Xs^2-2*X0*Xs+X0^2)

    #  2. with the source of the function abcd and pred
    ##################################################################

    #Provide set of initial parameters
    Xs.1 <- round(seq(X0.0 + 1, Xlast - 1, len = 30), digits = 0)
    X0.1 <- rep(X0.0, length(Xs.1))
    DFIs.1 <- NULL
    CFIs.1 <- NULL
    for(A in seq_along(Xs.1)){
      DFIs2 <- Data[Data$Age.plot == Xs.1[A],]$DFI.plot
      CFIs2 <- Data[Data$Age.plot == Xs.1[A],]$CFI.plot
      DFIs.1 <- c(DFIs.1, DFIs2)
      CFIs.1 <- c(CFIs.1, CFIs2)
    }

    st1 <- data.frame(cbind(X0.1, Xs.1, DFIs.1, CFIs.1))
    names(st1) <- c("X0","Xs", "DFIs","CFIs")

    #RUN NLS2 to find optimal initial parameters
    st2 <- nls2(Y ~ nls.func.2(X0, Xs, DFIs, CFIs),
                Data.xy,
                start = st1,
                # weights = weight,
                # trace = T,
               algorithm = "brute-force")

    par_init <- coef(st2); par_init

    #--------------------------------------------
    # Create empty lists to store data after loop
    #--------------------------------------------

    par <- list()
    AC.res <- list()
    AC.pvalue <- NULL
    data2 <- list()
    data3 <- list()
    param <- data.frame(rbind(par_init))
    par.abcd <- data.frame(rbind(abcd.2(as.vector(par_init))))
    param.2 <- data.frame(X0=double(),
                          Xs=double(),
                          DFIs=double(),
                          CFIs=double(),
                          a=double(),
                          b=double(),
                          c=double(),
                          stringsAsFactors=FALSE)
    j <- 2
    AC_pvalue <- 0
    AC.pvalue[1] <- AC_pvalue
    datapointsleft <- as.numeric(dim(Data)[1])
    dpl <- datapointsleft #vector of all dataponitsleft at each step

    #-------------------------------------------------------------------------------
    # Start the procedure of Non Linear Regression
    #-------------------------------------------------------------------------------

            while ((AC_pvalue<=0.05) && datapointsleft >= 20){
              weight <- 1/Y^2
              # ---------------- NON linear reg applied to log(Y) ---------------------------------
              st2 <- nls2(Y ~ nls.func.2(X0, Xs, DFIs, CFIs),
                          Data.xy,
                          start = st1,
                          weights = weight,
                          trace = F,
                          algorithm = "brute-force")
              par_init <- coef(st2)
              par_init
              # st1 <- st1[!(st1$Xs == par_init[2]),]
              nls.CFI <- nlsLM(Y ~ nls.func.2(X0, Xs, DFIs, CFIs),
                             Data.xy,
                             control = list(tol = 1e-2, printEval = TRUE, maxiter = 1024),
                             start = list(X0 = par_init[1], Xs = par_init[2],
                                          DFIs = par_init[3], CFIs = par_init[4]),
                             weights = weight,
                             algorithm = "port",
                             lower = c(-10000,X0.0+1, -10000, -10000),
                             upper = c(10000, Xlast-1, 10000, 10000),
                             trace = F)

              # nls.CFI <- nls2(Y ~ nls.func.2(X0, Xs, DFIs, CFIs),
              #                                 Data.xy,
              #                                 start = list(X0 = par_init[1], Xs = par_init[2],
              #                                              DFIs = par_init[3], CFIs = par_init[4]),
              #                                 weights = weight,
              #                                 control = nls.control(warnOnly = TRUE),
              #                                 trace = T,
              #                                 algorithm = "port",
              #                                 lower = c(-100000000,X0.0+1, -1000000000, -1000000000),
              #                                 upper = c(1000000000, Xlast-1, 1000000000, 1000000000))

              # nls.CFI <- nlsLM(Y ~ nls.func.2(X0, Xs, DFIs, CFIs),
              #                              Data.xy,
              #                              control = nls.control(warnOnly = TRUE),
              #                              start = list(X0 = par_init[1], Xs = par_init[2],
              #                                           DFIs = par_init[3], CFIs = par_init[4]),
              #                              weights = weight,
              #                              algorithm = "port",
              #                              lower = c(-1000000000,X0.0+1, -1000000000, -1000000000),
              #                              upper = c(1000000000, Xlast-1, 1000000000, 1000000000),
              #                              trace = F)


              #--------RESULTS analysis GOODNESS of fit
              #estimate params
              par[[j]] <- coef(nls.CFI)
              par.abcd[j,] <- abcd.2(as.vector(coef(nls.CFI) )) #calculation of a, b, c and d
              param[j,] <- par[[j]]
              param.2[j-1,] <- cbind(param[j,], par.abcd[j,])

              #summary
              # summ = overview((nls.CFI))  #summary
              #residuals
              res1 <- nlsResiduals(nls.CFI) #residuals
              res2 <- nlsResiduals(nls.CFI)$resi1
              res <- res2[, 2]
              AC.res <- test.nlsResiduals(res1)
              AC.pvalue[j] <- AC.res$p.value

              #---------Check for negative residuals----------

              #Add filtration step order to data
              Step <- rep(j - 1, length(x))
              #create a new dataset with predicted CFI included
              Data.new <- data.frame(cbind(x, Z, Y, pred.func.2(par[[j]],x)[[1]], res, Step))
              names(Data.new) <- c("Age", "Observed_DFI","Observed_CFI", "Predicted_CFI", "Residual", "Step")
              # plot(Data.new$Age, Data.new$Predicted_CFI, type = "l", col = "black",lwd = 2,
              #      ylim = c(0, max(Data.new$Predicted_CFI, Data.new$Observed_CFI)))
              # lines(Data.new$Age, Data.new$Observed_CFI, type = "p", cex = 1.5)

              #
              #remove negative res
              Data.pos <- Data.new[!Data.new$Residual<0,]
              # lines(Data.pos$Age, Data.pos$Predicted_CFI, type = "l", col = j-1, lwd = 2)
              # lines(Data.pos$Age, Data.pos$Observed_CFI, type = "p", col = j, cex = 1.5)
              #restart

              #Criteria to stop the loop when the estimated parameters are equal to initial parameters
              # Crite <- sum(param.2[dim(param.2)[1],c(1:4)] == par_init)

              datapointsleft <- as.numeric(dim(Data.pos)[1])
              par_init <- par[[j]]
              AC_pvalue <- AC.pvalue[j]
              j <- j+1
              x <- Data.pos$Age
              Y <- Data.pos$Observed_CFI
              Z <- Data.pos$Observed_DFI
              Data.xy <- as.data.frame(cbind(x,Y))
              dpl <- c(dpl, datapointsleft)
              dpl
              #Create again the grid
              X0.0 <- x[1]
              Xlast <- x[length(x)]
              #Xs
              if(par_init[2] -15 <= X0.0){
                Xs.1 <- round(seq(X0.0 + 5, Xlast - 5, len = 30), digits = 0)
              } else if(par_init[2] + 5 >= Xlast){
                Xs.1 <- round(seq(par_init[2]-10, par_init[2]-1, len = 6), digits = 0)
              } else{
                Xs.1 <- round(seq(par_init[2]-5, par_init[2] + 5, len = 6), digits = 0)
              }
              #
              X0.1 <- rep(X0.0, length(Xs.1))
              DFIs.1 <- NULL
              CFIs.1 <- NULL
              for(A in seq_along(Xs.1)){
                DFIs2 <- Data[Data$Age.plot == Xs.1[A],]$DFI.plot
                CFIs2 <- Data[Data$Age.plot == Xs.1[A],]$CFI.plot
                DFIs.1 <- c(DFIs.1, DFIs2)
                CFIs.1 <- c(CFIs.1, CFIs2)
              }
              st1 <- data.frame(cbind(X0.1, Xs.1, DFIs.1, CFIs.1))

              if(X0.0 <= par_init[2] && Xlast >=par_init[2]){
              st1 <- rbind(st1, par_init)
              }
              names(st1) <- c("X0","Xs", "DFIs","CFIs")
              }

      

  } # end FOR loop

这是数据文件。我已经将数据导出到.rdata以进行更容易的导入。： https://drive.google.com/file/d/1gvmarnkwmeyz-nops-nosp1dhzkqnzkqntu2ups3r/view?usp = sharing？usp = sharing

？代码>在此部分的初始参数估计的单数梯度矩阵：

 nls.CFI <- nlsLM(Y ~ nls.func.2(X0, Xs, DFIs, CFIs),
                             Data.xy,
                             control = list(tol = 1e-2, printEval = TRUE, maxiter = 1024),
                             start = list(X0 = par_init[1], Xs = par_init[2],
                                          DFIs = par_init[3], CFIs = par_init[4]),
                             weights = weight,
                             algorithm = "port",
                             lower = c(-10000,X0.0+1, -10000, -10000),
                             upper = c(10000, Xlast-1, 10000, 10000),
                             trace = F)

互补函数（文件function.r）：

abcd.2 <- function(P){
  X0 <- P[1]
  Xs <- P[2]
  DFIs <- P[3]
  CFIs <- P[4]
  
  a <- -X0*(2*CFIs*Xs-CFIs*X0-Xs^2*DFIs+Xs*DFIs*X0)/(Xs^2-2*X0*Xs+X0^2)
  b <- (-Xs^2*DFIs+DFIs*X0^2+2*CFIs*Xs)/(Xs^2-2*X0*Xs+X0^2)
  c <- -(CFIs-Xs*DFIs+X0*DFIs)/(Xs^2-2*X0*Xs+X0^2)
  
  pp <- as.vector(c(a, b, c))
  return(pp)
}

#--------------------------------------------------------------
# NLS function
#--------------------------------------------------------------

nls.func.2 <- function(X0, Xs, DFIs, CFIs){
  pp <- c(X0, Xs, DFIs, CFIs)
  
  #calculation of a, b and c using these new parameters
  c <- abcd.2(pp)[3]
  b <- abcd.2(pp)[2]
  a <- abcd.2(pp)[1]
  
  ind1 <- as.numeric(x < Xs)
  return (ind1*(a+b*x+c*x^2)+(1-ind1)*((a+b*(Xs)+c*(Xs)^2)+(b+2*c*(Xs))*(x-(Xs))))
}

#--------------------------------------------------------------
# Fit new parameters to a quadratic-linear function of CFI
#--------------------------------------------------------------

pred.func.2 <- function(pr,age){
  #
  X0 <- pr[1]
  Xs <- pr[2]
  DFIs <- pr[3]
  CFIs <- pr[4]
  #
  x <- age
  #calculation of a, b and c using these new parameters
  c <- abcd.2(pr)[3]
  b <- abcd.2(pr)[2]
  a <- abcd.2(pr)[1]
  
  #
  ind1 <- as.numeric(x < Xs)
  #
  results <- list()
  cfi <- ind1*(a+b*x+c*x^2)+(1-ind1)*((a+b*(Xs)+c*(Xs)^2)+(b+2*c*(Xs))*(x-(Xs)))  #CFI
  dfi <- ind1*(b+2*c*x) + (1 - ind1)*(b+2*c*(Xs))                                 #DFI
  results[[1]] <- cfi
  results[[2]] <- dfi
  return (results)
}

#---------------------------------------------------------------------------------------------------------------
# Quadratic-linear function of CFI curve and its 1st derivative (DFI) with original parameters (only a, b and c)
#---------------------------------------------------------------------------------------------------------------

pred.abcd.2 <- function(pr,age){
  #
  a <- pr[1]
  b <- pr[2]
  c <- pr[3]
  
  x <- age
  
  #calculation of a, b and c using these new parameters
  #
  ind1 <- as.numeric(x < Xs)
  #
  results <- list()
  cfi <- ind1*(a+b*x+c*x^2)+(1-ind1)*((a+b*(Xs)+c*(Xs)^2)+(b+2*c*(Xs))*(x-(Xs)))  #CFI
  dfi <- ind1*(b+2*c*x) + (1 - ind1)*(b+2*c*(Xs))                                 #DFI
  results[[1]] <- cfi
  results[[2]] <- dfi
  return (results)
}

更新：我确实从上一步，发现我的数据因此有些混乱。我已经修复了。集合F数据流入无限循环不再存在的情况下，但是此错误仍然存​​在：在初始参数估计的奇异梯度矩阵。