[R] writing selfStart models that can deal with treatment effects

[Michael A. Gilchrist]

Hello,

I'm trying to do some non-linear regression with 2 cell types and 4 tissue 
type treatments using selfStart models

Following Ritz and Streibig (2009), I wrote the following routines:

##Selfstart
expDecayAndConstantInflowModel <- function(Tb0, time, aL, aN, T0){
   exp(-time*aL)*(T0*aL+(-1+exp(time * aL))*Tb0 * aN)/aL
}

expDecayAndConstantInflowModelInit <- function(mCall, LHS, data){
   print(paste("SelfStart mCall:", mCall));
   print(attributes(mCall));
   print(mCall[["aN"]]);
   xy <-  sortedXyData(mCall[["time"]], LHS, data);
   lmFit <- lm(log(xy[, "y"]) ~ xy[, "x"]);
   coefs <- coef(lmFit);
   T0 <- exp(coefs[1]);
   aL <- -coefs[2]*0.99;
   aN <- 0.0001;
   value <- c(aL, aN, T0);
   names(value) <- mCall[c("aL", "aN", "T0")];
   value
}

SSexpDecayAndConstantInflow <- selfStart(expDecayAndConstantInflowModel, 
expDecayAndConstantInflowModelInit, c("aL", "aN", "T0"))

Ignoring the treatment effects, the routines seem to work find:
-----------------------------------
> getInitial(Count ~ SSexpDecayAndConstantInflow(B0, Time, aL, aN, T0), data = 
tissueData)
           aL           aN           T0
4.600144e-02 1.000000e-04 1.082172e+03
> nls(Count ~ SSexpDecayAndConstantInflow(B0, Time, aL, aN, T0), data = 
tissueData)
Nonlinear regression model
   model:  Count ~ SSexpDecayAndConstantInflow(B0, Time, aL, aN, T0)
    data:  tissueData
        aL        aN        T0
    0.1131    0.1206 1348.0646
  residual sum-of-squares: 38821317

Number of iterations to convergence: 11
Achieved convergence tolerance: 6.267e-06
> 
----------------------------------

However, I need to be able to test for these treatment effects (they are quite 
heterogenous), but when I try using the standard syntax for nls(), the 
selfStart model only estimates a set of global values and these don't work 
with nls()

> getInitial(Count ~ SSexpDecayAndConstantInflow(B0, Time, aL[Type], aN[Type], 
T0[TypeTissue]), data = tissueData)
       aL[Type]       aN[Type] T0[TypeTissue]
   4.600144e-02   1.000000e-04   1.082172e+03
> nls(Count ~ SSexpDecayAndConstantInflow(B0, Time, aL[Type], aN[Type], 
T0[TypeTissue]), 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 ~ SSexpDecayAndConstantInflow(B0, Time, aL[Type], aN[Type],  :
   No starting values specified for some parameters.
Intializing 'aL', 'aN', 'T0' to '1.'.
Consider specifying 'start' or using a selfStart model
>
----------------------------------

I've tried the following syntax, but it seems that the global estimates give 
singularities when applied to each group

----------------------------------
> nls(Count ~ SSexpDecayAndConstantInflow(B0, Time, aL[Type], aN[Type], 
T0[TypeTissue]),
     data = tissueData,
     start = getInitial(Count ~ SSexpDecayAndConstantInflow(B0, Time, aL, aN, 
T0), data = tissueData))

nls(Count ~ SSexpDecayAndConstantInflow(B0, Time, aL[Type], aN[Type], 
T0[TypeTissue]),
+     data = tissueData,
+     start = getInitial(Count ~ SSexpDecayAndConstantInflow(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
>
----------------------------------

Does anyone have any suggestions for how to work around this problem or have 
an example of a selfStart model that can deal with treatment effects?  The 
dataset is quite heterogenous and so I need to explore a number of alternative 
transformations and models so I really need something that works without me 
having to continually reset the start values by hand.

Many thanks.

Mike

-----------------------------------------------------
Department of Ecology & Evolutionary Biology
569 Dabney Hall
University of Tennessee
Knoxville, TN 37996-1610

phone:(865) 974-6453
fax:  (865) 974-6042

web: http://eeb.bio.utk.edu/gilchrist.asp