suggestions for nls error: false convergence

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

suggestions for nls error: false convergence

Thread (6 messages)