Thanks very much. I said I knew there were better ways, but none that I could develop in the time available this morning. You've helped introduce me to other options.
Best Wishes, Spencer Graves
Peter Dalgaard BSA wrote:
Spencer Graves <[EMAIL PROTECTED]> writes:
The problem with nls is that it is NOT maximum likelihood for the Poisson distribution. For the Poisson, the standard deviation is the square root of the mean, while nls assumes constant standard deviation. That's why I stayed with glm. The answers should be comparable, but I would prefer the glm answers. spencer graves
That's why I was suggesting either a variance stabilizing
transformation or using gnls() with a weights function. In both cases
you need to watch out for scaling of SE's stemming from the fact that
the variance is supposed to be known in the Poisson distribution, but
the fitting routines estimate a sigma from the residuals anyway.
I.e. for instance:
Phyto.nls2 <- nls(sqrt(y+.5) ~ sqrt(Ymax/(1 + x/x50)+.5),data=Phytopath,start=list(Ymax=20.0,x50=0.01),trace=T)
9.400602 : 20.00 0.01 0.9064723 : 27.21511020 0.03593643 0.06338235 : 28.21654101 0.05995647 0.02616589 : 28.50221759 0.06815302 0.02608587 : 28.54871243 0.06835046 0.02608586 : 28.54960242 0.06834083 0.02608586 : 28.5495605 0.0683413
summary(Phyto.nls2)
Formula: sqrt(y + 0.5) ~ sqrt(Ymax/(1 + x/x50) + 0.5)
Parameters: Estimate Std. Error t value Pr(>|t|) Ymax 28.54956 1.65113 17.291 0.0368 * x50 0.06834 0.01264 5.407 0.1164 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 0.1615 on 1 degrees of freedom
Correlation of Parameter Estimates: Ymax x50 -0.6833
but here you need to know that the theoretical sd is 0.5, so the Std. errors all need multiplication by 0.5/0.1615.
Using gnls() we'd get (if I got the calling sequence right...)
Phyto.gnls <- gnls(y ~ Ymax/(1 + x/x50),
+ data=Phytopath,weights=varPower(fixed=.5),start=list(Ymax=20.0,x50=0.01))
summary(Phyto.gnls)
Generalized nonlinear least squares fit Model: y ~ Ymax/(1 + x/x50) Data: Phytopath AIC BIC logLik 13.71292 11.00875 -3.856458
Variance function: Structure: Power of variance covariate Formula: ~fitted(.) Parameter estimates: power 0.5
Coefficients: Value Std.Error t-value p-value Ymax 29.393796 2.4828723 11.838626 0.0536 x50 0.063028 0.0125244 5.032376 0.1249
Correlation: Ymax x50 -0.84
Standardized residuals: [1] -0.4894266 0.7621749 -0.4237346 attr(,"std") [1] 2.8478142 1.4239071 0.8586483 attr(,"label") [1] "Standardized residuals"
Residual standard error: 0.6367906 Degrees of freedom: 3 total; 1 residual
and again, it is necessary to adjust the SE's by multiplying with 1/0.6368
It is in fact rather easy to do direct MLE for this kind of model:
with(Phytopath, optim(c(28,.7),fn=function(p){Ymax<-p[1];x50<-p[2];
+ -sum(dpois(y,lambda=Ymax/(1 + + x/x50),log=TRUE))}, method="BFGS", hessian=T)) $par [1] 28.55963487 0.06833524
$value [1] 7.21251
$counts function gradient 42 13
$convergence [1] 0
$message NULL
$hessian [,1] [,2] [1,] 0.07356072 6.631868 [2,] 6.63186764 1294.792539
Warning message: NaNs produced in: dpois(x, lambda, log)
and we can get SE's from the inverse hessian with
sqrt(diag(solve(with(Phytopath,
+ optim(c(28,.7),fn=function(p){Ymax<-p[1];x50<-p[2]; + -sum(dpois(y,lambda=Ymax/(1 + + x/x50),log=TRUE))}, method="BFGS", hessian=T))$hessian))) [1] 5.02565893 0.03788052
Notice that the variance stabilizing transform seems to do rather better than gnls() compared to true MLE. I'm slightly puzzled by this, but presumably it has to do with the fact that MLE for a Gaussian model with a mean/variance relationship is *not* identical to the iteratively reweighted least squares that glm() et al. are doing. (I wouldn't quite rule out that I've simply made a mistake somewhere, though...)
______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
