[R] custom selfStart model works with getInitial but not nls

Michael A. Gilchrist mikeg at utk.edu
Sat Oct 17 15:58:08 CEST 2009


On Sat, 17 Oct 2009, Douglas Bates wrote:

> On Fri, Oct 16, 2009 at 7:09 PM, Michael A. Gilchrist <mikeg at utk.edu> wrote:
>> Hello,
>>
>> I'm having problems creating and using a selfStart model with nlme.
>>  Briefly, I've defined the model, a selfStart object, and then combined them
>> to make a selfStart.default model.
>>
>> If I apply getInitial to the selfStart model, I get results.  However, if I
>> try usint it with nls or nlsList, these routines complain about a lack of
>> initial conditions.
>>
>> If someone could point out what I'm doing wrong, I'd greatly appreciate it.
>
> Thanks for providing the code.  Could you also provide some data as a test case?
>

Thanks for looking into this.  I've learned a lot from reading your book 
on nlme (and still have a lot to learn).  I've posted the data I've been 
working with at:
 	www.tiem.utk.edu/~mikeg/software/R/selfStart/tissue.rda

Mike

>>
>> Details:
>> ## Nonlinear model I want to fit to the data
>> const.PBMC.tcell.model <- function(B0, t, aL, aN, T0){
>>
>>  Tb0 = B0;
>>
>>  x = exp(-log(aL) + log(T0*aL+(-1+exp(t * aL))*Tb0 * aN) - t * aL);
>>
>>  return(x);
>> }
>>
>>
>> ##Define selfStart routine
>> const.PBMC.tcell.selfStart<- function(mCall, LHS, data){
>>
>>  t0 = 0;
>>  t1 = 24;
>>  t2 = 48;
>>
>>  ##Get B0 Value
>>  B0 =  data[1, "B0"];
>>
>>  T0 = mean(data[data$Time==t0, "Count"]);
>>  T1 = mean(data[data$Time==t1, "Count"]);
>>  T2 = mean(data[data$Time==t2, "Count"]);
>>
>>  if(T0 < T2){ ##increase -- doesn't work
>>    stop(paste("Error in const.PBMC.tcell.start: T0 < T2 for data: ", data[1,
>> ]));
>>
>>  }
>>  ##Estimate aL based on exponential decline from t=0 to t=24
>>  aLVal = -(log(T1) - log(T0))/(t1-t0);
>>
>>  ##Estimate aNVal based on final value
>>  aNVal = aLVal*T2/B0;
>>
>>  values = list(aLVal, aNVal, T0);
>>  names(values) <- mCall[c("aL", "aN", "T0")]; #mimic syntax used by P&B
>>  return(values)
>> }
>>
>>
>>
>> ##Now create new model with selfStart attributes
>> const.PBMC.tcell.modelSS<-  selfStart(model = const.PBMC.tcell.model,
>> initial=const.PBMC.tcell.selfStart)
>>
>>
>> ##Test routines using getInitial -- This works
>>>
>>> getInitial(Count ~ const.PBMC.tcell.modelSS(B0, Time,aL, aN, T0), data =
>>> tissueData)
>>
>> [1] 0.05720924
>> $aL
>> [1] 0.05720924
>>
>> $aN
>> [1] 0.1981895
>>
>> $T0
>> [1] 1360.292
>>
>> ##Now try to use the SS model -- this doesn't work
>>>
>>> nls(Count ~ const.PBMC.tcell.modelSS(B0, Time,aL, aN, T0), data =
>>> tissueData)
>>
>> Error in numericDeriv(form[[3L]], names(ind), env) :
>>  Missing value or an infinity produced when evaluating the model
>> In addition: Warning message:
>> In nls(Count ~ const.PBMC.tcell.modelSS(B0, Time, aL, aN, T0), data =
>> tissueData) :
>>  No starting values specified for some parameters.
>> Intializing 'aL', 'aN', 'T0' to '1.'.
>> Consider specifying 'start' or using a selfStart model
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>


More information about the R-help mailing list