# Re: [R] suggestions for nls error: false convergence

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Mon 19 Dec 2005 - 10:09:56 EST

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@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@pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help