[Rd] confint/nls

[Ben Bolker]

   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)