s-news
[Top] [All Lists]

Sum: Problems porting R selfStart SSweibull to S

To: "'s-news@lists.biostat.wustl.edu'" <s-news@lists.biostat.wustl.edu>
Subject: Sum: Problems porting R selfStart SSweibull to S
From: "Yang, Richard" <dyang@NRCan.gc.ca>
Date: Thu, 22 Nov 2001 11:33:40 -0500
Many thanks to Drs. Jose Pinheiro and Doulas Bates for their responses and
suggestions. Both pointed out the problem has to do with different scope
rules in R and S-PLUS and and to
minor differences in syntax (Dr. Pinheiro).

The scoping problem was solved by either breaking down the original pars
assignment in whInit() into two (Dr. Pinheiro):

        xy$ytrans <- log(-log((Rasym - xy$y)/(Rasym - Lasym)))
        pars <- coef(lm(ytrans ~ log(x), data = xy, subset = x > 0))

or  without modifying the pars assigment but adding two assignments(Dr.
Bates) making Rasym and Lasym visible:

        assignment("Rasym", NLSstRtAsymptote(xy), frame=1)
      assignment("Lasym", NLSstLfAsymptote(xy), frame=1)

In either case, the "start =" in val assignment needs modified as follows:

        val <- coef(nls(y ~ cbind(1, -exp(-exp(lrc) * x^pwr)), data = xy,
      alg = "plinear", start = list(lrc = pars[1], pwr = pars[2])))[c(3,
4,1, 2)]

The selfStart function SSwb works fine in S-Plus following the
modifications. Dr. Pinheiro will include this function in the future NLME
release.

Again, thank to Drs. Pinheiro and Bates.

Richard


> -----Original Message-----
> From: Yang, Richard [mailto:dyang@nrcan.gc.ca]
> Sent: Tuesday, November 20, 2001 4:18 PM
> To: 's-news@lists.biostat.wustl.edu'
> Subject: [S] Problems porting R selfStart SSweibull to S
> 
> 
> S++ers;
> 
> I experienced difficulties porting the SSweibull function in 
> R to S (Splus 6
> or 2000). The code looks like the following:
> 
>       >wb <- deriv( ~ Asym - Drop*exp(-exp(lrc)*x^pwr), 
> c("Asym", "Drop",
> "lrc", "pwr"),
>     function(x, Asym, Drop, lrc, pwr) {})
> 
>       >wbInit <-             ## from Dr. Doug Bates
>       function (mCall, data, LHS) 
>       {
>     xy <- sortedXyData(mCall[["x"]], LHS, data)
>     if (nrow(xy) < 5) {
>         stop("Too few distinct input values to fit the Weibull growth
> model")
>     }
>     if (any(xy[["x"]] < 0)) {
>         error("All x values must be non-negative to fit the 
> Weibull growth
> model")
>     }
>     Rasym <- NLSstRtAsymptote(xy)
>     Lasym <- NLSstLfAsymptote(xy)
>     pars <- coef(lm(log(-log((Rasym - y)/(Rasym - Lasym))) ~ 
>         log(x), data = xy, subset = x > 0))
>     val <- coef(nls(y ~ cbind(1, -exp(-exp(lrc) * x^pwr)), data = xy, 
>         alg = "plinear", start = c(lrc = pars[[1]], pwr = 
> pars[[2]])))[c(3, 
>         4, 1, 2)]
>     names(val) <- mCall[c("Asym", "Drop", "lrc", "pwr")]
>     val
>       }
> 
>       >SSwb <- selfStart(wb, initial = wbInit)
> 
> Following the creation of SSwb selfStart function, a test 
> with the Orange
> data resulted in an error:
> 
>       > getInitial(circumference ~ SSwb(age, Asym, Drop, lrc, 
> pwr), data =
> 
>       Orange)
>       Error: Object "Rasym" not found
> 
> It sounded like the NLSstRtsymptote() was not working. However,
> 
>       >Orange.sortAvg <- sortedXyData("age", "circumference", Orange)
> 
>       >NLSstRtAsymptote(Orange.sortAvg)
>       >[1] 193.9
>       > NLSstLfAsymptote(Orange.sortAvg)
>       >[1] 12.9
> 
> The SSweibull function works fine with the Orange dataset in 
> R, but the same
> error occurred in both Splus 6 and Splus 2000. 
> 
>       What did I miss? Please advise. 
> 
>       TIA
> 
> Richard
> 
> ---------------------------------------------------------------------
> This message was distributed by s-news@lists.biostat.wustl.edu.  To
> unsubscribe send e-mail to s-news-request@lists.biostat.wustl.edu with
> the BODY of the message:  unsubscribe s-news
> 

<Prev in Thread] Current Thread [Next in Thread>
  • Sum: Problems porting R selfStart SSweibull to S, Yang, Richard <=