confint/nls

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

confint/nls

Ben Bolker-2

   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)

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: confint/nls

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:

>
>   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)
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
>

--
Brian D. Ripley,                  [hidden email]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel