SSposnegRichards: Self-Starting Positive-Negative Richards Model (double-Richards)

Description

This selfStart function evaluates a range of flexible logistic functions. It also has an initial attribute that creates initial estimates of the parameters for the model specified.

Usage

SSposnegRichards(x,

Asym = NA,

K = NA,

Infl = NA,

M = NA,

RAsym = NA,

Rk = NA,

Ri = NA,

RM = NA,

modno,

pn.options)

Arguments

a numeric vector of the primary predictor variable at which to evaluate the model

Asym

a numeric value for the asymptote of the positive (increasing) curve

a numeric value for the rate parameter of the positive (increasing) curve

Infl

a numeric value for the point of inflection of the positive (increasing) curve

a numeric value for the shape parameter of the positive (increasing) curve

RAsym

a numeric value for the asymptote of the negative (decreasing) curve

a numeric value for the rate parameter of the negative (decreasing) curve

a numeric value for the point of inflection of the negative (decreasing) curve

a numeric value for the shape parameter of the negative (decreasing) curve

modno

a numeric value (currently integer only) between 1 and 36 specifying the identification number of the equation to be fitted

pn.options

character string specifying name of list object populated with starting parameter estimates, fitting options and bounds

Value

a numeric vector of the same length as x containing parameter estimates for equation specified (fixed parameters are not return but are substituted in calls to nls nlsList and nlme with the fixed parameters stored in pnmodelparams; see modpar

Details

This selfStart function evaluates a range of flexible logistic functions. It also has an initial attribute that creates initial estimates of the parameters for the model specified. Equations can fit both monotonic and non-monotonic curves (two different trajectories). These equations have also been described as double-Richards curves, or positive-negative Richards curves. **From version 1.2 onwards this function can fit curves that exhibit negative followed by positive trajectories or double positive or double negative trajectories.*** The 32 possible equations (plus custom model #17) are all based on the subtraction of one Richards curve from another, producing: $y = A / ([1+ m exp(-k (t-i))]1/m) + A' / ([1+ m' exp(-k' (t-i' ))]1/m' )$, where A=Asym, k=K, i=Infl, m=M, A'=RAsym, k'=Rk, i'=Ri, m'=RM; as described in the Arguments section above. All 32 possible equations are simply reformulations of this equation, in each case fixing a parameter or multiple parameters to (by default) the mean parameter across all individuals in the dataset (such as produced by a nls model). Thus, a model in which one parameter is fixed has a 7-parameter equation, and one in which four are fixed has a 4-parameter equation, thus reducing complexity and computation and compensatory parameter changes when a parameter does not vary across group levels (e.g individuals) [the most appropriate equation can be determined using model selection in pn.modselect.step or pn.mod.compare]. Any models that require parameter fixing (i.e. all except #1) extract appropriate values from the list object specified by pn.options for the fixed parameters. This object is most easily created by running modpar and can be adjusted manually or by using change.pnparameters to user-required specification. Each of the 32 equations is identified by an integer value for modno (1 to 36).Models 21-36 are the same as 1-16 except that in the former the first curve parameter m is fixed and not estimated. All equations (except 17 - see below) contain parameters Asym, K, and Infl. The list below summarizes which of the other 5 parameters are contained in which of the models (Y indicates that the parameter is estimated, blank indicates it is fixed). modno M RAsym Rk Ri RM NOTES 1 Y Y Y Y Y 8 parameter model 2 Y Y Y Y 3 Y Y Y 4 Y Y Y 5 Y Y 6 Y Y Y Y 7 Y Y Y Y 8 Y Y Y Y 9 Y Y Y 10 Y Y 11 Y Y 12 Y 4 parameter, standard Richards model 13 Y Y Y 14 Y Y Y 15 Y Y Y 16 Y Y 17 see below 18 see below 19 see below 20 see below 21 Y Y Y Y 7 parameter model, 4 recession params 22 Y Y Y 6 parameters (used for double-logistic, double-Gompertz or double-Von Bertalannfy) 23 Y Y 24 Y Y 25 Y 26 Y Y Y 27 Y Y Y 28 Y Y Y 29 Y Y 30 Y 31 Y 32 only 3 parameters (used for logistic, Gompertz or Von Bertalannfy - see below) 33 Y Y 34 Y Y 35 Y Y 36 Y modno 17 represents a different parameterization for a custom model: (Asym/ 1 + exp(Infl - x)/ M) - (RAsym / 1 + exp(Ri - x)/ RM), in which M and RM actually represent scale parameters not shape parameters. This model and a suite of reductions: modno 17.1: modno 17.2: modno 17.3: are designed for use in modeling migration, $sensu.$ Bunnefeld et al. 2011, Singh et al. In press. modnos 18, 19 and 20 are reserved for internal use by modpar. To access common 3 parameter sigmoidal models use modno = 32, fixing parameters (using change.pnparameters) to M = 1 for logistic, M = 0.1 for Gompertz, and M= -0.3 for von Bertalanffy. The same settings can be used with modno = 2 to fit the double-logistic, double-Gompertz or double-Von Bertalannfy. Note that to fit only 3 or 4 parameter curves, the option force4par = TRUE should be specified when running modpar. The call for SSposnegRichards only differs from conventional selfStart models in that it requires a value for modno and a list of fitting options and values from modpar to which to fix parameters in the reduced models. Depending on the model chosen, different combinations of the 8 possible parameters are required: if one is missing the routine will stop with an appropriate error, if an extra one is added, it will be ignored (provided that it is labelled, e.g. M = 1; this is good practice to prevent accidental misassignments). Here are two examples (7 parameter and 3 parameter): richardsR2.lis <- nlsList(mass ~ SSposnegRichards(age, Asym = Asym, K = K, Infl = Infl, M = M, RAsym = RAsym, Rk = Rk, Ri = Ri, , modno = 2), data = posneg.data) #correct call includes all necessary parameters richardsR20.lis <- nlsList(mass ~ SSposnegRichards(age, Asym = Asym, K = K, Infl = Infl, modno = 2), data = logist.data) #incorrect call missing required parameters, function terminates and generates an error message Examples for all models can be found in the list object posnegRichards.calls. If specified using modpar optional constraints may be placed to specify response values at the minimum value and/or maximum values of the predictor. Such constraint allows realistic fits for datasets which are missing data at either end of the curve (e.g. hatching weight for some growth curves). Estimates are produced by splitting the two curves into separate positive and negative curves and estimating parameters for each curve separately in a similar manner to SSlogis. Each curve is fit first by optim using the parameter bounds in pnmodelparamsbounds (see modpar) and a subsequent refinement is attempted using nls with more restrictive parameter bounds. Finally, both curves are annealed and parameters are again estimated using restrictive bounds and starting values already determined during separate estimates. Equations for which the positive curve was inestimable are not estimated further, but if negative curve estimation or overall curve estimation fail, partial estimates are used: either default negative parameters (RAsym = 0.05*Asym, Rk = K, Ri = Infl, RM = M) annealed to positive curves or separate estimates annealed; both with compensation for interation between asymptotes. From version 1.2 onwards, this function can now fit two component models, where the first curve is used up to the x-value (e.g. age) of intersection and the second curve is used afterwards. Confusingly, these are also called "double Richards", "double Gompertz" or "double logistic": see Murphy et al. (2009) or Ross et al. (1995) for examples. To specify such models set twocomponent.x = TRUE (this will estimate the x of intersection) when running modpar. Alternatively, if known, the x of intersection can be set directly by setting twocomponent.x = # (where # is the x of intersection). Whne {code{modpar} is run this option will be saved in pnmodelparas and can be changed at will, either manually or using change.pnparameters. From version 1.2 onwards, this function can now fit bilogistic (and more generally biRichards) curves, where the final curve is effectively either two positive curves or two negative curves. See Meyer (1994) for examples. This functionality is default and does not need to be specified.

References

## Oswald, S.A. et al. (Early view) FlexParamCurve: R package for flexible fitting of nonlinear parametric curves. Methods in Ecology and Evolution. doi: 10.1111/j.2041-210X.2012.00231.x (see also tutorial and introductory videos at: http://www.methodsinecologyandevolution.org/view/0/podcasts.html (posted September 2012 - if no longer at this link, check the archived videos at: http://www.methodsinecologyandevolution.org/view/0/VideoPodcastArchive.html#allcontent #1# Nelder, J.A. (1962) Note: an alternative form of a generalized logistic equation. Biometrics, 18, 614-616. #2# Huin, N. & Prince, P.A. (2000) Chick growth in albatrosses: curve fitting with a twist. Journal of Avian Biology, 31, 418-425. #3# Meyer, P. (1994) Bi-logistic growth. Technological Forecasting and Social Change. 47: 89-102 #4# Murphy, S. et al. (2009) Importance of biological parameters in assessing the status of Delphinus delphis. Marine Ecology Progress Series 388: 273-291. #5# Pinheiro, J. & Bates, D. (2000) Mixed-Effects Models in S and S-Plus. Springer Verlag, Berlin. #6# Ross, J.L. et al. (1994) Age, growth, mortality, and reproductive biology of red drums in North Carolina waters. Transactions of the American Fisheries Society 124: 37-54. #7# Bunnefeld et al. (2011) A model-driven approach to quantify migration patterns: individual, regional and yearly differences. Journal of Animal Ecology 80: 466-476. #8# Singh et al. (In press) From migration to nomadism: movement variability in a northern ungulate across its latitudinal range. Ecological Applications.

Examples

Run this code

set.seed(3) #for compatability issues

 require(graphics)

    # retrieve mean estimates of 8 parameters using getInitial

    # and posneg.data object

    data(posneg.data)

    modpar(posneg.data$age, posneg.data$mass,verbose=TRUE, pn.options = "myoptions", wide.bounds=TRUE)

    getInitial(mass ~ SSposnegRichards(age, Asym, K, Infl, M, 

        RAsym, Rk, Ri, RM, modno = 1, pn.options = "myoptions"), data = posneg.data)



    # retrieve mean estimates and produce plot to illustrate fit for 

    # curve with M, Ri and Rk fixed

    pars <- coef(nls(mass ~ SSposnegRichards(age, 

        Asym = Asym, K = K, Infl = Infl, RAsym = RAsym, 

        	RM = RM, modno = 24, pn.options = "myoptions"), data = posneg.data,

        	control=list(tolerance = 10)))

    plot(posneg.data$age, posneg.data$mass, pch=".")

    curve(posnegRichards.eqn(x, Asym = pars[1], K = pars[2], 

        Infl = pars[3], RAsym = pars[4],  

        RM = pars[5],  modno = 24, pn.options = "myoptions"), xlim = c(0, 

        200), add = TRUE)

    

    # retrieve mean estimates and produce plot to illustrate fit for custom model 17

    \dontrun{    pars<-as.numeric( getInitial(mass ~ SSposnegRichards(age, Asym, K, Infl,

           M, RAsym, Rk, Ri, RM, modno = 17, pn.options = "myoptions"), data = datansd) )

     plot(datansd$jday21March, datansd$moosensd)

     curve( posnegRichards.eqn(x, Asym = pars[1], K = 1, Infl = pars[2], 

            M = pars[3], RAsym = pars[4], Rk = 1, Ri = pars[5], RM = pars[6], 

            modno = 17, pn.options = "myoptions"), lty = 3, xlim = c(0, 200) , add = TRUE)}

        

    # fit nls object using 8 parameter model

    # note: ensure data object is a groupedData object

    richardsR1.nls <- nls(mass ~ SSposnegRichards(age, Asym = Asym, 

        K = K, Infl = Infl, M = M, RAsym = RAsym, Rk = Rk, Ri = Ri, 

        RM = RM, modno = 1, pn.options = "myoptions"), data = posneg.data)

        

    # fit nlsList object using 8 parameter model

    # note: ensure data object is a groupedData object

    # also note: not many datasets require all 8 parameters

    \dontrun{

    richardsR1.lis <- nlsList(mass ~ SSposnegRichards(age, Asym = Asym, 

        K = K, Infl = Infl, M = M, RAsym = RAsym, Rk = Rk, Ri = Ri, 

        RM = RM, modno = 1, pn.options = "myoptions"), data = posneg.data)

    summary(richardsR1.lis)}



    # fit nlsList object using 6 parameter model with value M and RM

    # fixed to value in pnmodelparams and then fit nlme model

    # note data is subset to provide estimates for a few individuals

    # as an example

    subdata <- subset(posneg.data,posneg.data$id == as.character(36)

   		| posneg.data$id == as.character(9) 

   		| posneg.data$id == as.character(32) 

   		| posneg.data$id == as.character(43))

    richardsR22.lis <- nlsList(mass ~ SSposnegRichards(age, Asym = Asym, 

        K = K, Infl = Infl, RAsym = RAsym, Rk = Rk, Ri = Ri, 

        modno = 22, pn.options = "myoptions"), data = subdata)

    summary(richardsR22.lis )

    richardsR22.nlme <- nlme(richardsR22.lis, random = pdDiag(Asym + Infl ~ 1) )

    summary(richardsR22.nlme)

         

    # fit nls object using simple logistic model, with 

    # M, RAsym, Rk, Ri, and RM fixed to values in pnmodelparams

    data(logist.data)

    modpar(logist.data$age, logist.data$mass ,force4par = TRUE, pn.options = "myoptions")

    change.pnparameters(M = 1, pn.options = "myoptions") #set to logistic (M =1) prior to fit

    richardsR32.nls <- nls(mass ~ SSposnegRichards(age, Asym = Asym, 

        K = K, Infl = Infl, modno = 32, pn.options = "myoptions"), data = logist.data)

    coef(richardsR32.nls)

                

    # fit a two component model - enter your own data in place of "mydata"

    \dontrun{

    modpar(mydata$x,mydata$y,twocomponent.x=TRUE, pn.options = "myoptions") # if x of intersection unknown

    modpar(mydata$x,mydata$y,twocomponent.x=75, pn.options = "myoptions") # if x of intersection = 75

    richardsR1.nls <- nls(y~ SSposnegRichards(x, Asym = Asym, K = K,

      Infl = Infl, M = M, RAsym = RAsym, Rk = Rk, Ri = Ri, RM = RM, modno = 1, pn.options = "myoptions")

                       , data = mydata)

    coef(richardsR1.nls)}

Run the code above in your browser using DataLab