0

I am trying to run non-linear least squares using a custom function following this example The data set has a x, predictor y, and groups (a/b). I've only worked with MLLs code previously. I'm having a hard time wrapping my head around how to invoke the function into the nls and how to do subsequent nls fit. I have attached the code that I tried along with errors where applicable.

# generate data
set.seed(2)
x <- rep(sample(1:100, 9), 12)
set.seed(15)
r <- rnorm(x, mean = 0, sd = 200*x^-0.8)
y <- c(200, 300) * exp(c(-0.08, -0.05)*x) + c(120, 100) + r
df <-  data.frame(x = x, y = y, 
                  group = rep(letters[1:2], length.out = length(x)))

ga <- function(aa, K, Ka, q,c){
  p <- x / K
  lnqp <- if (q == 0) log(p) else (p^q - 1) / q
  y <- (aa * (p * K / Ka - 1) - 1) * lnqp + c
}

m <- nls(y ~ ga(aa, K, Ka, q ,c))

Error in q == 0 : 
  comparison (1) is possible only for atomic and list types
In addition: Warning message:
In nls(y ~ ga(aa, K, Ka, q, c)) :
  No starting values specified for some parameters.
Initializing ‘aa’, ‘K’, ‘Ka’ to '1.'.
Consider specifying 'start' or using a selfStart model

summary(m) 
plot(m) 

#ERROR: Error in q == 0 : 
comparison (1) is possible only for atomic and list types


# fit generalised nonlinear least squares
require(nlme)
mgnls <- gnls(y ~ ga(aa, K, Ka, q ,c), 
              start = list(aa = c(356.6851,356.6851), 
                           K = c(-exp(-2.9356),-exp(-2.9356)),
                           Ka = c(-exp(-2.9356),-exp(-2.9356)),
                           q = c(-exp(-2.9356),-exp(-2.9356)),
                           c = c(108.9860,108.9860)),
              params = list(aa ~ group, K ~ group, Ka ~ group, q ~ group, c ~ group),
              weights = varExp(),
              data = df)
plot(mgnls) # more homogenous
IRTFM
  • 258,963
  • 21
  • 364
  • 487
Rspacer
  • 2,369
  • 1
  • 14
  • 40
  • Your `ga` function should take `x` as an input--probably as the first argument (I'm not sure that it matters that it's first, but that would follow the example closely.) – Gregor Thomas May 10 '21 at 19:41
  • I don't think so. I tried it with and without x, the errors don't change. – Rspacer May 10 '21 at 19:52
  • 2
    If you are asking a debugging question then ALWAYS include full text of error message. and do so at the point where the error occurs. (It's rather pointless (and quite possibly confusing) to include code that you intend to execute after the debugging process is complete.) – IRTFM May 10 '21 at 20:01
  • 3
    The "comparison (1) is possible only for atomic and list types" error is because the `nls()` model doesn't know that `q` is a parameter, it is using the built-in function `q()` (which quits are), and that's the error you get when compare a function to 0, e.g. ,`q == 0` or `median == 0`. If you give it a list of starting values, that error goes away (though you get a different error instead about infinite values... you might need to do something to make sure `K` is positive so that `log(p)` is finite. – Gregor Thomas May 10 '21 at 20:01
  • same issue with the `c` – jeremycg May 10 '21 at 20:04

0 Answers0