CODEWITHSUNDEEP
rcurve-fittingnls
K. Maya
K. Maya

Reputation: 299

Fit data with multiple functions

I am trying to fit my data with multiple fitting functions consisting of multiple variables. An example is below for two variables. There are cases where for a certain variable the fit is not good and I get a singular gradient error. I would like to ignore those cases and proceed anyway and furthermore for the remaining variables take the better solution between the two fitting function by comparing the deviance. Like in this example for both type1 and type2 the sum of residuals is less with the first function sum(resid(myfitfun1)^2) < sum(resid(myfitfun2)^2) so take the first function for both of the variables.

myfun1<-function(x,a,b){1/(1+exp(-(x/a)+b))}
myfun2<-function(x,a,b){1+b*exp(-(x)/a)} 
   
mydata <- data.frame(v=c("type1","type1","type1","type1","type1","type1","type1","type1","type1","type1","type1",
"type1","type1","type1","type1","type1","type1","type1","type1","type1","type1","type1","type1","type1",
 "type1","type1","type1","type1","type1","type1","type1","type1","type1","type1","type1","type1",
 "type1","type1","type1","type1","type1","type1","type1","type1","type2","type2","type2","type2",
                     "type2","type2"), 
 m=c(1.116975672,1.38641493,1.423833959,1.482008121,1.513588607,1.527179033,
                     1.543512539,1.555874185,1.607579807,1.721182154,1.729059048,1.748226681,
                     1.774814055,1.815147988,1.835638421,1.854582642,1.861972,1.887704144,
                     1.915360975,1.948689331,1.97516491,1.985962227,2.011310496,2.043716548,
                     2.068918438,2.091184665,2.120366813,2.126865141,2.148241318,2.15871161,
                     2.193529738,2.256197915,2.302364722,2.316381935,2.31909683,2.325213451,
                     2.336299128,2.410419652,2.473160411,2.478302702,2.5238233,2.651124474,
                     2.70470831,2.927536062,-0.1736072,0.1235610,0.5848941,0.9016486,0.9744832,
                     1.2767238), 
 n=c(0.022727273,0.045454545,0.068181818,0.090909091,0.113636364,0.136363636,
                     0.159090909,0.181818182,0.204545455,0.227272727,0.25,0.272727273,0.295454545,
                     0.318181818,0.340909091,0.363636364,0.386363636,0.409090909,0.431818182,
                     0.454545455,0.477272727,0.5,0.522727273,0.545454545,0.568181818,0.590909091,
                     0.613636364,0.636363636,0.659090909,0.681818182,0.704545455,0.727272727,0.75,
                     0.772727273,0.795454545,0.818181818,0.840909091,0.863636364,0.886363636,
                     0.909090909,0.931818182,0.954545455,0.977272727,1,0.1666667,0.3333333,0.5000000,
                     0.6666667,0.8333333,1))
myfitfun1 <- nls(n~myfun1(m,a,b),mydata,start=list(a=1,b=1))
myfitfun2 <- nls(n~myfun2(m,a,b),mydata,start=list(a=1,b=1))

I would like to program it in a way that it handels automatically the better fit between the two functions for various type and ignoring in case of an error. Any help is appreciated.

Upvotes: 1

Views: 414

Answers (1)

jay.sf
jay.sf

Reputation: 72593

You could put both functions in a function and work with tryCatch. In the one tryCatches, just throw an NA to overcome the error. In another tryCatch set the value to Inf when an error occurs to ensure the "better" fit for the non-failing function. In normal cases the minimum is chosen. With `attr<-` we can give the MSE as an attribute to the output of the fit.

fun <- function(data) {
  myfitfun1 <- tryCatch(
    nls(n ~ myfun1(m, a, b), data, start=list(a=1, b=1)),
    error=function(e) NA)
  myfitfun2 <- tryCatch(
    nls(n ~ myfun2(m, a, b), data, start=list(a=1, b=1)),
    error=function(e) NA)
  L <- list(myfitfun1, myfitfun2)
  res <- sapply(L, function(x) {
    tryCatch(sum(resid(x)^2), error=function(e) Inf)
    })
  `attr<-`(L[[which.min(res)]], "MSE", min(res))
}

fun(mydata)
# Nonlinear regression model
# model: n ~ myfun1(m, a, b)
# data: data
# a      b 
# 0.3465 5.6461 
# residual sum-of-squares: 2.323
# 
# Number of iterations to convergence: 26 
# Achieved convergence tolerance: 7.675e-06

To get the MSE attribute, use:

attr(fun(mydata), "MSE")
# [1] 2.322945

Upvotes: 1

Related Questions