Michael A. Gilchrist
2009-Nov-12 20:39 UTC
[R] writing selfStart models that can deal with treatment effects
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