# [R] suggestions for nls error: false convergence

Spencer Graves spencer.graves at pdf.com
Mon Dec 19 00:09:56 CET 2005

```Hi, Christian:

That's very interesting.  May I ask how "multdrc" converged and found
all 4 parameters statistically significant when "nls" failed and "optim"
told me the hessian was singular?

To try to study this myself, I first compared your numbers with what
I got and found that you had swapped the third and fourth parameters.
When I swapped them back, I got the same estimates but substantially
different standard errors from your function:

logi <- function(dose, parm){
parm[, 1] * (1+parm[, 2]*exp(-dose/parm[, 3])) / (1+parm[,
4]*exp(-dose/parm[, 3]))}
logi. <- function(dose, parm){
parm[, 1] * (1+parm[, 2]*exp(-dose/parm[, 4])) / (1+parm[,
3]*exp(-dose/parm[, 4]))}

## Fitting the function to the data (see ?multdrc for details)
library(drc)
model <- multdrc(y~x, fct = list(logi, NULL, c("a", "m", "tau", "n")),
startVal=c(277, 100, 10, 101))
model. <- multdrc(y~x, fct = list(logi., NULL, c("a", "m", "n", "tau")),
startVal=c(277, 100, 101, 10))
signif(summary(model)\$estimates, 2)
Estimate Std. Error t-value  p-value
a:(Intercept)    2.8e+02    9.6e-01   290.0 4.5e-143
m:(Intercept)    1.5e+03    6.0e+02     2.6  1.2e-02
tau:(Intercept)  5.6e+00    9.6e-02    59.0  2.7e-77
n:(Intercept)    2.3e+05    4.7e+04     4.8  5.1e-06
> signif(summary(model.)\$estimates, 2)
Estimate Std. Error t-value  p-value
a:(Intercept)    2.8e+02    8.5e-01   320.0 4.2e-148
m:(Intercept)    1.5e+03    4.8e+02     3.2  1.7e-03
n:(Intercept)    2.3e+05    3.3e+04     6.9  6.0e-10
tau:(Intercept)  5.6e+00    6.8e-02    83.0  1.2e-91

I then tried "optim" again, rescaling all parameters by your
estimated standard errors.  If your code does what I thought it might,
would expect that the residual variance times the main diagonal of the
inverse of the Hessian should be all 1's.  That's not what I got, so I'm
confused.  Here's what I did:

rescale <- summary(model.)\$estimates[,2]

func2 <- function( par,y, x, rescale ) {
par <- rescale*par
a = par
m = par
n = par
tau = par
y. <- a * (1+m*exp(-x/tau)) / (1+n*exp(-x/tau))
sum((y-y.)^2)
}

est.no2 <- optim(est.n\$par/rescale,
func2,  hessian=TRUE, y=y, x=x, rescale=rescale)

round(var.resid <- est.no2\$value/96, 2)
 8.42
est.no2.hessEig <-
eigen(est.no2\$hessian, symmetric=TRUE)

> round(est.no2.hessEig\$values, 1)
 6134.1   52.4   24.2    4.7
# MUCH better than before.

hessInv <-
with(est.no2.hessEig, vectors %*%(t(vectors)/values))
> var.par <- (est.no2\$value/96)*hessInv
> round(sqrt(diag(var.par)), 2)
 0.72 0.76 0.79 0.77

What am I missing?
Thanks again for your work in the "drc" package and for your earlier

Spencer Graves

Christian Ritz wrote:
> Hi.
>
> An alternative is to use the package 'drc' on CRAN to fit your data!
>
> x <- 1:100
>
> y <- c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,1,1,1,2,2,2,2,2,3,4,4,4,5,
> 5,5,5,6,6,6,6,6,8,8,9,9,10,13,14,16,19,21,
> 24,28,33,40,42,44,50,54,69,70,93,96,110,127,127,141,157,169,
> 178,187,206,216,227,236,238,244,246,250,255,255,257,260,261,262,266,268,
> 268,270,272,272,272,273,275,275,275,276)
>
>
> ## Defining the function (in a bit different notation)
> logi <- function(dose, parm){parm[, 1] * (1+parm[, 2]*exp(-dose/parm[, 3])) / (1+parm[, 4]*exp(-dose/parm[, 3]))}
>
> ## Fitting the function to the data (see ?multdrc for details)
> library(drc)
> model <- multdrc(y~x, fct = list(logi, NULL, c("a", "m", "tau", "n")), startVal=c(277, 100, 10, 101))
>
> summary(model)
> plot(model, log="")
>
>
> Hope this helps?
>
> Christian
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help

--
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

spencer.graves at pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915

```