[R] suggestions for nls error: false convergence

[Spencer Graves]

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[1]
m = par[2]
n = par[3]
tau = par[4]
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)
[1] 8.42
est.no2.hessEig <-
eigen(est.no2$hessian, symmetric=TRUE)

 > round(est.no2.hessEig$values, 1)
[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)
[1] 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 
reply to my comments on this example.

	  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
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

-- 
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