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@(protected)
> 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@(protected)
www.pdf.com <http://www.pdf.com>
Tel: 408-938-4420
Fax: 408-280-7915
______________________________________________
R-help@(protected)
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
