[Rd] nmle: gnls freezes on difficult case

Nicholas Lewin-Koh nikko at hailmail.net
Wed Oct 17 20:49:52 CEST 2007


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"))



More information about the R-devel mailing list