CODEWITHSUNDEEP
rnonlinear-functions
Brian Smith
Brian Smith

Reputation: 1353

Fitting non-linear equation to data using base R package

let say I have below data

Data = structure(list(col1 = c(31, 66, 88, 123, 249, 362, 488, 610, 
730, 842), col2 = c(2101.58953918969, 2103.57391509821, 2100.3292541732, 
2101.64107993765, 2100.51743895393, 2100.16708521627, 2102.1992412748, 
2101.06516854423, 2101.87929065226, 2101.25318636023)), row.names = c(NA, 
-10L), class = "data.frame")

Now I want to fit a non-linear equation as below -

library(stats)
nls(col2 ~ x1 + x2 / (1 + exp(-x3 * (col1 - x4))), data = Data, start = list(x1 = 0, x2 = 0, x3 = 0, x4 = 0), algorithm = "plinear")

However with this I am getting below error -

Error in qr.qty(QR.rhs, .swts * ddot(attr(rhs, "gradient"), lin)) : 
  NA/NaN/Inf in foreign function call (arg 5)

Can you please help me to understand what went wrong in my approach?

I want to use only base package to fit this equation as I can not download any contributed package from internet in my system.

Any pointer will be highly appreciated.

Upvotes: 1

Views: 200

Answers (2)

G. Grothendieck
G. Grothendieck

Reputation: 270298

There are a few problems:

  1. When using plinear the right hand side should be a matrix such that the ith column of that matrix multiplies the ith parameter that enters linearly and those linear parameters should not have starting values. The linear parameters will be reported as .lin1 and .lin2
  2. Better starting values for the nonlinear parameters are needed.
nls(col2 ~ cbind(1, 1 / (1 + exp(-x3 * (col1 - x4)))), data = Data, 
  start = list(x3 = sd(Data$col1), x4 = mean(Data$col1)), algorithm = "plinear")

giving:

Nonlinear regression model
  model: col2 ~ cbind(1, 1/(1 + exp(-x3 * (col1 - x4))))
   data: Data
       x3        x4     .lin1     .lin2 
 295.3813  358.9000 2101.5302   -0.2175 
 residual sum-of-squares: 9.145

Number of iterations to convergence: 0 
Achieved convergence tolerance: 0

Upvotes: 3

Ben Bolker
Ben Bolker

Reputation: 226911

If I use SSfpl with your current data I can get an answer.

n1 <- nls(col2 ~ SSfpl(col1, A, B, m, s), data=Data)

pframe <- data.frame(col1=seq(0,900,length=101))
pframe$col2 <- predict(n1, newdata=pframe)

library(ggplot2); theme_set(theme_bw())
ggplot(Data, aes(col1,col2)) + geom_point() + geom_smooth() +
  geom_line(data=pframe, colour="red")

data, loess curve + CI, and fitted four-parameter logistic

The parameterization is not quite the same as yours:

         A          B          m          s 
2001.56354 2002.06645  642.30178   20.76013

Based on x1 + x2 / (1 + exp(-x3 * (col1 - x4))),

I believe x4 = m (midpoint), x3 = s (scale), x1 = A (left asymptote), and x2 = B-A (B is the right asymptote).

Upvotes: 2

Related Questions