SSposnegRichards: Self-Starting Positive-Negative Richards Model

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)

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

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 and their gradients. It also has an initial attribute that creates initial estimates of the parameters for the model specified. All 36 possible equations are biomodal, in that they have a positive (increasing) trajectory followed by a negative (decreasing) trajectory. 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 36 possible equations 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 36 possible equations are simply reformulations of this equation, in each case fixing a parameter or multiple parameters to 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 when a parameter does not vary across groups (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 object pnmodelparams for the fixed parameters. This object is created by running modpar and can be adjusted manually or by using change.pnparameters to user required specification. Each of the 36 equations is identified by an integer value for modno (1 to 36). 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 only 3 parameters (used for logistic, Gompertz or Von Bertalannfy - see below) 21 Y Y Y Y 7 parameter model, 4 recession params 22 Y Y Y 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 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. modnos 18 and 19 are identical to 1 and 12, respectively, but are used in the initial estimation of pnmodelparams and are not intended for use in themselves. To access common 3 parameter sigmoidal models use modno = 20, fixing parameters (using change.pnparameters) to M = 1 for logistic, M = 0.1 for Gompertz, and M= -0.3 for von Bertalanffy. 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. 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 = M; 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) richardsR20.lis <- nlsList(mass ~ SSposnegRichards(age, Asym = Asym, K = K, Infl = Infl, modno = 2), data = logist_data) 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) anealed 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 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_age = TRUE (this will estimate the age of intersection) when running modpar. Alternatively, if known, the age of intersection can be set directly by setting twocomponent_age = # (where # is the age 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

#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.

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)

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

        RAsym, Rk, Ri, RM, modno = 1), data = posneg_data)



    # retrieve mean estimates and produce plot to illustrate fit for 

    # curve with RM fixed

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

        Asym = Asym, K = K, Infl = Infl, RAsym = RAsym, Rk = Rk, M = M,

        	Ri = Ri, modno = 2), data = posneg_data))

    plot(posneg_data$age, posneg_data$mass, pch=".")

    curve(posnegRichards_eqn(x, Asym = pars[1], K = pars[2], 

        Infl = pars[3], M = pars[4], RAsym = pars[5], Rk = pars[6], 

        Ri = pars[7],  modno = 2), 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), 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), 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), 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), data = posneg_data)}



    # 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), data = subdata)

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

         

    # 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) # create pnmodelparams

    change.pnparameters(M = 1) #set to logistic (M =1) prior to fit

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

        K = K, Infl = Infl, modno = 20), data = logist_data)

                

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

    \dontrun{

    modpar<-(mydata$x,mydata$y,twocomponent_age=TRUE) # if age of intersection unknown

    modpar<-(mydata$x,mydata$y,twocomponent_age=75) # if age 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)

                       , data = mydata)}

Run the code above in your browser using DataLab

Last chance! 50% off unlimited learning