FlexParamCurve (version 1.5-6)

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,
Envir = ".GlobalEnv")

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

Arguments

x

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

K

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

M

a numeric value for the shape parameter of

the positive (increasing) curve

RAsym

a numeric value for the asymptote of the

negative (decreasing) curve

Rk

a numeric value for the rate parameter of

the negative (decreasing) curve

Ri

a numeric value for the point of inflection

of the negative (decreasing) curve

RM

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

Envir

a character vector that represents the valid R environment in which to

find pn.options in and write any output to,

by default this is the global environment

Author

Stephen Oswald <[e-mail]>

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 parameter (double-logistic/double-Gompertz/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. 2012.

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). When 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. (2012) FlexParamCurve: R package for flexible

fitting of nonlinear parametric curves. Methods in Ecology and Evolution 3: 1073-1077.

     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. (2012) From migration to nomadism: movement variability in a northern

ungulate across its latitudinal range. Ecological Applications 22: 2007-2020.

See Also

SSlogis

SSgompertz

posnegRichards.eqn

Examples

 set.seed(3) #for compatability issues

 require(graphics)

    # retrieve mean estimates of 8 parameters using getInitial

    # and posneg.data object

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

    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)

    



    

        # following example not run as appropriate data are not available in the package

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

     if (FALSE) { 

     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

   # \donttest{ 
   
        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)# }

        



    # following example not run as it fits very few levels in these data - as noted

    # such a comprehensive equation is rarely required

    # fit nlsList object using 8 parameter model

    # note: ensure data object is a groupedData object

    # also note: not many datasets require all 8 parameters

         if (FALSE) {

         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
    
    # \donttest{ 

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

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

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

    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

    # \donttest{ 
    
    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"

    # this is not run for want of an appropriate dataset

    # if x of intersection unknown

    if (FALSE) {

    modpar(mydata$x,mydata$y,twocomponent.x=TRUE, pn.options = "myoptions")

    # if x of intersection = 75

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

    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)}