Skip to content

confint/nls

2 messages · Ben Bolker, Brian Ripley

Original
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)
1 day later
Brian Ripley
# 
These are all solved, except those for "plinear", where you cannot profile 
linear parameters and so you must specify parms (or call 
confint(profile(fm)).  And the third plinear model does not converge, so 
isn't a useful test.
On Sat, 7 Jan 2006, Ben Bolker wrote: