[R] non-linear regression and root finding

[Ivan Krylov]

В Mon, 6 Nov 2023 17:53:49 +0100
Troels Ring <tring using gvdnet.dk> пишет:

> Hence I wonder if I could somehow have non linear regression to find
> the 3 pK values. Below is HEPESFUNC which delivers charge in the
> fluid for known pKs, HEPTOT and SID. Is it possible to have
> root-finding in the formula with nls?

Sure. Just reformulate the problem in terms of a function that takes a
vector of predictors (your independent variable SID) and the desired
parameters (pK1, pK2, pK3) as separate arguments and then returns
predicted values of the dependent variable (to compare against pHobs):

kw <- 1e-14 # I'm assuming
pHm <- Vectorize(function(SID, pK1, pK2, pK3)
 -log10(uniroot(HEPESFUNC,c(1e-20,1),tol=1e-20,maxiter=1e4,
        HEPTOT=HEPTOT,SID = SID, pK1=pK1,pK2=pK2,pK3=pK3)$root))

(Yes, Vectorize() doesn't make the function any faster, but I don't see
an easy way to rewrite this function to make it truly vectorised.)

Unfortunately, nls() seems to take a step somewhere where crossprod()
of the Jacobian of the model function cannot be inverted and fails with
the message "Singular gradient". I wish that R could have a more
reliable built-in nonlinear least squares solver. (I could also be
holding it wrong.) Meanwhile, we have excellent CRAN packages nlsr and
minpack.lm:

minpack.lm::nlsLM(
 pHobs ~ pHm(SID, pK1, pK2, pK3),
 data.frame(pHobs = pHobs, SID = SID),
 start = c(pK1 = pK1, pK2 = pK2, pK3 = pK3),
 # the following is also needed to avoid MINPACK failing to fit
 lower = rep(-1, 3), upper = rep(9, 3)
)
# Nonlinear regression model
#   model: pHobs ~ pHm(SID, pK1, pK2, pK3)
#    data: data.frame(pHobs = pHobs, SID = SID)
#    pK1     pK2    pK3
# -1.000   2.966  7.606
#  residual sum-of-squares: 0.001195
#
# Number of iterations to convergence: 15 
# Achieved convergence tolerance: 1.49e-08

(Unfortunately, your code seemed to have a lot of non-breakable spaces,
which confuse R's parser and make it harder to copy&paste it.)

I think that your function can also be presented as a degree-5
polynomial in H, so it should also be possible to use polyroot() to
obtain your solutions in a more exact manner.

-- 
Best regards,
Ivan