I have found some "issues" (bugs?) with nls confidence intervals ...
some with the relatively new "port" algorithm, others more general
(but possibly in the "well, don't do that" category). I have
corresponded some with Prof. Ripley about them, but I thought I
would just report how far I've gotten in case anyone else has
thoughts. (I'm finding the code in stats/nls.R and stats/nle-profile.R
quite dense & scary ...)
All of this has been done with R-devel from 3 Jan 2006; the changes
that Prof. Ripley already made to allow confint.nls not to crash
when algorithm="port" are in R-devel, not R-patched.
a synopsis of the problems with confint():
with a 1-parameter model (is confint not appropriate for 1-parameter
models? it doesn't say so in the docs [by the way, "normality" is
misspelled as "nornality" in ?confint]):
algorithm=default or plinear: get a complaint from qr.qty ('qr' and
'y' must have the same number of rows)
port: "cannot allocate vector of size [large]" [caused by C code
looking for dims when they aren't there]
2-parameter models:
default OK
port "cannot allocate vector"
plinear "Error in xy.coords"
3-parameter models are OK
I can fix the 2-parameter port case by adding drop=FALSE in
appropriate places, but I wanted to check in just in case
there are better/more efficient ways than my slogging through
one case at a time ...
apologies for the long message, but I am temporarily cut
off from any way to post these files to the web.
cheers
Ben Bolker
code that tests various combinations of numbers of parameters
and algorithms:
-----------
resmat = array(dim=c(3,2,3),
dimnames=list(npar=1:3,c("fit","confint"),c("default","plinear","port")))
resmat.fix <- resmat
## sim. values
npts=1000
set.seed(1001)
x = runif(npts)
b = 0.7
y = x^b+rnorm(npts,sd=0.05)
a =0.5
y2 = a*x^b+rnorm(npts,sd=0.05)
c = 1.0
y3 = a*(x+c)^b+rnorm(npts,sd=0.05)
d = 0.5
y4 = a*(x^d+c)^b+rnorm(npts,sd=0.05)
testfit <- function(model,start,alg) {
tryfit <- try(fit <-
nls(model,start=start,algorithm=alg,control=list(maxiter=200)))
if (class(tryfit)!="try-error") {
fitcode="OK"
tryci <- try(confint(fit))
if (class(tryci)!="try-error") {
cicode="OK"
} else cicode = as.character(tryci)
} else {
fitcode = as.character(tryfit)
cicode="?"
}
c(fitcode,cicode)
}
m1 = c(y~x^b,y2~a*x^b,y3~a*(x+exp(logc))^b)
m2 = c(y2~x^b,y3~(x+exp(logc))^b,y4~(x^d+exp(logc))^b)
s1 = list(c(b=1),c(a=1,b=1),c(a=1,b=1,logc=0))
s2 = list(c(b=1),c(b=1,logc=0),c(b=1,logc=0,d=0.5))
for (p in 1:3) {
resmat[p,,"default"] <- testfit(m1[[p]],start=s1[[p]],alg=NULL)
resmat[p,,"port"] <- testfit(m1[[p]],start=s1[[p]],alg="port")
}
for (p in 1:3) {
resmat[p,,"plinear"] <- testfit(m2[[p]],start=s2[[p]],alg="plinear")
}
print(resmat)
set.seed(1002)
example(nls,local=TRUE)
diffs:
--
*** /usr/local/src/R/R-devel/src/library/stats/R/nls.R 2006-01-07
10:57:08.000000000 -0500
--- nlsnew.R 2006-01-07 19:18:53.000000000 -0500
***************
*** 266,277 ****
gradSetArgs[[1]] <- (~attr(ans, "gradient"))[[2]]
gradCall <-
switch(length(gradSetArgs) - 1,
! call("[", gradSetArgs[[1]], gradSetArgs[[2]]),
! call("[", gradSetArgs[[1]], gradSetArgs[[2]],
gradSetArgs[[2]]),
call("[", gradSetArgs[[1]], gradSetArgs[[2]],
gradSetArgs[[2]],
! gradSetArgs[[3]]),
call("[", gradSetArgs[[1]], gradSetArgs[[2]],
gradSetArgs[[2]],
! gradSetArgs[[3]], gradSetArgs[[4]]))
getRHS.varying <- function()
{
ans <- getRHS.noVarying()
--- 266,277 ----
gradSetArgs[[1]] <- (~attr(ans, "gradient"))[[2]]
gradCall <-
switch(length(gradSetArgs) - 1,
! call("[", gradSetArgs[[1]], gradSetArgs[[2]],drop=FALSE),
! call("[", gradSetArgs[[1]], gradSetArgs[[2]],
gradSetArgs[[2]],drop=FALSE),
call("[", gradSetArgs[[1]], gradSetArgs[[2]],
gradSetArgs[[2]],
! gradSetArgs[[3]],drop=FALSE),
call("[", gradSetArgs[[1]], gradSetArgs[[2]],
gradSetArgs[[2]],
! gradSetArgs[[3]], gradSetArgs[[4]],drop=FALSE))
getRHS.varying <- function()
{
ans <- getRHS.noVarying()
***************
*** 331,337 ****
else {
vary
}, envir = thisEnv)
! gradCall[[length(gradCall)]] <<- useParams
if(all(useParams)) {
assign("setPars", setPars.noVarying, envir = thisEnv)
assign("getPars", getPars.noVarying, envir = thisEnv)
--- 331,337 ----
else {
vary
}, envir = thisEnv)
! gradCall[[length(gradCall)-1]] <<- useParams
if(all(useParams)) {
assign("setPars", setPars.noVarying, envir = thisEnv)
assign("getPars", getPars.noVarying, envir = thisEnv)
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
