[Rd] nmle: gnls freezes on difficult case

From: Nicholas Lewin-Koh <nikko_at_hailmail.net>
Date: Wed, 17 Oct 2007 14:49:52 -0400

Hi,
I am not sure this is a bug but I can repeat it, The functions and data are below.
I know this is nasty data, and it is very questionable whether a 4pl model
is appropriate, but it is data fed to an automated tool and I would have hoped for an error. Does this repeat for anyone else? My details:
> version

               _                           
platform       i686-pc-linux-gnu           
arch           i686                        
os             linux-gnu                   
system         i686, linux-gnu             
status                                     
major          2                           
minor          6.0                         
year           2007                        
month          10                          
day            03                          
svn rev        43063                       
language       R                           
version.string R version 2.6.0 (2007-10-03)

start=c(-1.5, 9.5, 0.09, 10.25)
names(start)<-c("A","B","xmid","scal")

gnls(response~SSllogis(conc,A,B,xmid,scal),tdat,start=start,weights=varPower(),verbose=TRUE)

**Iteration 1
GLS step: Objective: NULLvarStruct parameters:

    power
0.3373199

NLS step: RSS = 0
 model parameters:-0.799941 8.99983 -0.522623 212.314  iterations: 2

Convergence:

   params varStruct
 1.172208 1.000000
### After about 10 min hit cntrl C

Warning messages:
1: In log(xmid) : NaNs produced
2: In log(xmid) : NaNs produced
3: In log(xmid) : NaNs produced




#####################################################################
#### Data
tdat<-data.frame(conc=c(0.00203,0.0061,0.0183,0.0549,0.165,0.494,1.48,4.44,13.3,40),
                 response=c(12,-4,19,11,-5,-3,1,6,0,-8))
#### Self start function
SSllogis <- selfStart(~ A + (B-A)/(1 + exp(scal*(log(x)-log(xmid)))),   function(mCall, data, LHS)
  {
    #browser()
    xy <- sortedXyData(mCall[["x"]], LHS, data)     if (nrow(xy) < 5) {
        stop("too few distinct input values to fit a four-parameter
        logistic")
      }

    rng <- range(xy$y)
    drng <- diff(rng)
    B <- rng[2] + 0.001
    A <- rng[1] - 0.001
    xy$prop <- log((B-xy$y)/(xy$y-A+0.001))     #(xy$y - rng[1] + 0.05 * drng)/(1.1 * drng)     ir <- as.vector(coef(lm(prop ~ log(x), data = xy)))     scal <- ir[2]
    xmid <- exp(-ir[1]/ir[2])
    #pars <- as.vector(coef(nls(y ~ cbind(1, 1/(1 + exp((xmid - 
    #    x)/exp(lscal)))), data = xy, start = list(xmid = ir[1], 
    #    lscal = log(abs(ir[2]))), algorithm = "plinear")))
    value <- c(A, B, xmid, scal)
    names(value) <- mCall[c("A", "B", "xmid", "scal")]     value

  }, c("A", "B", "xmid", "scal"))



R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Wed 17 Oct 2007 - 18:54:36 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 25 Oct 2007 - 11:37:11 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.