[R] non-linear regression and root finding

[J C Nash]

Your script is missing something (in particular kw).

I presume you are trying to estimate the pK values. You may have more success
with package nlsr than nls(). nlsr::nlxb() tries to get the Jacobian of the
model specified by a formula and do so by applying symbolic or automatic
differentiation. The multi expression function would probably not work.
(That is one advantage of nls(), but it has an unstabilized solver and unless
carefully called will use a simple forward derivative approximation.)
There is also nlsr::nlfb() that can use functions for the residuals and
jacobian. Your function is messy, but likely the jacobian can be developed
with a bit of work.

Some of the symbolic derivative features of R can help here. Alternatively,
there are several numerical approximations in nlsr and the vignette "Intro
to nlsr" explains how to use these. Note that simple forward and backward
approximations are generally not worth using, but central approximation is
quite good.

The solver in the nlsr package allows of bounds on the parameters. Since you
want them ordered, you may want to transform

    pK1 < pK2 <- pK3

to  pk1, deltapk2, deltapk3 so pK2 = pk1+deltapk2
and pk3 = pk2 + deltapk3 = pk1 + deltapk2 + deltapk3
and all values bounded below by 0 or possibly some small numbers to
keep the paramters apart.

I won't pretend any of this is trivial. It is VERY easy to make small errors
that still allow for output that is disastrously incorrect. If it is possible
to get a single "formula" as one expression even if spread over multiple lines,
then nlxb() might be able to handle it.


J Nash (maintainer of nlsr and optimx)

On 2023-11-06 11:53, Troels Ring wrote:

....
> HEPESFUNC <-
>    function(H,SID,HEPTOT,pK1,pK2,pK3) {
>      XX <- (H^3/(10^-pK3*10^-pK2*10^-pK1)+H^2/(10^-pK3*10^-pK2)+H/(10^-pK3))
>      IV <- HEPTOT/(1+XX)
>      I <- IV*H^3/(10^-pK3*10^-pK2*10^-pK1)
>      II <- IV*H^2/(10^-pK3*10^-pK2)
>      III <- IV*H/10^-pK3
>      H - kw/H + SID + I*2 + II - IV
>    }
> 
> 
> HEPTOT <- 0.050
> SID <- c(-seq(10,1)*1e-4,0,seq(1,10)*1e-4)
> 
> pHobs  <- c(4.63,4.68,4.72,4.77,4.83,4.9,4.96,5.04,5.12,5.21,5.3,
>              5.39,5.48,5.55,5.63,5.69,5.74,5.8,5.85,5.89,5.93)
> 
> pK1 <- -1; pK2 <- 3; pK3 <- 7.55 # literature values
> pK3 <- 7.6; pK2 <- 2.96 # values eye-balled to be better
> 
> pH <- c()
> for (i in 1:length(SID)) {
>    pH[i] <- -log10(uniroot(HEPESFUNC,c(1e-20,1),tol=1e-20,maxiter=1e4,
>                         HEPTOT=HEPTOT,SID = SID[i],
>                         pK1=pK1,pK2=pK2,pK3=pK3)$root)
> }
> 
> plot(SID,pH)
> points(SID,pHobs,col="red")