CODEWITHSUNDEEP
r
jn025
jn025

Reputation: 2895

R nls function and starting values

I'm wondering how I can find/choose the starting values for the nls function as I'm getting errors with any I put in. I also want to confirm that I can actually use the nls function with my data set.

    data
 [1] 108 128  93  96 107 126 105 137  78 117 111 131 106 123 112  90  79 106 120
[20]  91 100 103 112 138  61  57  97 100  95  92  78

week = (1:31)

> data.fit = nls(data~M*(((P+Q)^2/P)*exp((P+Q)*week)/(1+(Q/P)*exp(-(P+Q)*week))^2), start=c(M=?, P=?, Q=?))

Upvotes: 1

Views: 3107

Answers (1)

G. Grothendieck
G. Grothendieck

Reputation: 269421

If we change the function a bit and use nls2 to get starting values then we can get it to converge. The model we are using is:

log(data) = .lin1 + .lin2 * log((exp((P+Q)*week)/(1+(Q/P)*exp(-(P+Q)*week))^2))) +error

In this model .lin1 = log(M*(((P+Q)^2/P)) and when .lin2=1 it reduces to the model in the question (except for the multiplicative rather than additive error and the fact that the parameterization is different but when appropriately reduced gives the same predictions). This is a 4 parameter rather than 3 parameter model.

The linear parameters are .lin1 and .lin2. We are using algorithm = "plinear" which does not require starting values for these parameters. The RHS of plinear formulas is specified as a matrix with one column for each linear parameter specifying its coefficient (which may be a nonlinear function of the nonlinear parameters).

The code is:

data <- c(108, 128, 93, 96, 107, 126, 105, 137, 78, 117, 111, 131, 106, 
123, 112, 90, 79, 106, 120, 91, 100, 103, 112, 138, 61, 57, 97, 
100, 95, 92, 78)
week <- 1:31
if (exists("fit2")) rm(fit2)
library(nls2)

fo <- log(data) ~ cbind(1, log((exp((P+Q)*week)/(1+(Q/P)*exp(-(P+Q)*week))^2)))

# try maxiter random starting values
set.seed(123)
fit2 = nls2(fo, alg = "plinear-random", 
     start = data.frame(P = c(-10, 10), Q = c(-10, 10)),
     control = nls.control(maxiter = 1000))

# use starting value found by nls2
fit = nls(fo, alg = "plinear", start = coef(fit2)[1:2])
plot(log(data) ~ week)
lines(fitted(fit) ~ week, col = "red")

giving:

> fit
Nonlinear regression model
  model: log(data) ~ cbind(1, log((exp((P + Q) * week)/(1 + (Q/P) * exp(-(P + Q) * week))^2)))
   data: parent.frame()
       P        Q    .lin1    .lin2 
 0.05974 -0.02538  5.63199 -0.87963 
 residual sum-of-squares: 1.069

Number of iterations to convergence: 16 
Achieved convergence tolerance: 9.421e-06

enter image description here

Upvotes: 3

Related Questions