Learn R Programming

vitality (version 1.3)

vitality.k: Fitting routine for the 3-parameter vitality model.

Description

This function provides the fitting routine for the 3-parameter vitality model. Intrinsic mortality is characterized by the mean (r) and variability (s) in vitality loss rate. Extrinsic mortality is characterized by the frequency (k) of lethal random challenges. Model is appropriate to animal mortality data (e.g. Anderson 2000).

Usage

vitality.k(time, sdata, rc.data=F, se=F, gfit=F, datatype="CUM", ttol=.000001, 
    init.params=F, lower=c(0,-1,0), upper=c(100,50,50), pplot=T, tlab="days", 
    lplot=F, cplot=F, Iplot=F, silent=F)

Arguments

time

Vector. Time component of data: Defaults to 0:(length(sdata)-1). Typically this refers to ages.

sdata

Required. Survival or mortality data. The default expects cumulative survival fraction. If providing incremental mortality fraction instead, use option: datatype = "INC". The default also expects the data to represent full mortality. Otherwise, use option: rc.data = T to indicate right censored data.

rc.data

Optional, Boolean. Specifies Right Censored data. If the data does not represent full mortality, it is probably right censored. The default is rc.data = F. A third option is rc.data = "TF". Use this case to add a near-term zero survival point to data which displays nearly full mortality ( <.01 survival at end). If rc.data = F but the data does not show full mortality, rc.data = "TF" will be invoked automatically.

se

Optional, Boolean. Calculates the standard errors for the MLE parameters. Default is FALSE. Set equal to the initial study population to compute standard errors.

datatype

Optional. Defaults to "CUM" for cumulative survival fraction data. Use "INC" - for incremental mortality fraction data.

ttol

Optional. Stopping criteria tolerance. Default is 1e-6. Specify as ttol = .0001. If one of the likelihood plots (esp. for "k") does not look optimal, try decreasing ttol. If the program crashes, try increasing ttol.

init.params

Optional. Please specify the initial param values. specify init.params = c(r, s, k) in that order (e.g.. init.params = c(.1, .02, .3)).

lower

vector of lower parameter bounds in order of c(r, s, k). see nlminb

upper

vector of upper parameter bounds in order of c(r, s, k). see nlminb

pplot

Optional, Boolean. Plots of cumulative survival for both data and fitted curves? Default: TRUE. FALSE produces no plots.Note: the incremental mortality plot is a continuous representation of the appropriately binned histogram of incremental mortalities.

Iplot

Boolean. Plot incremental survival? Must have pplot=TRUE

lplot

Boolean. Plot likelihood functions? Provides likelihood function plotting. Defaults to FALSE. Note: these plots are not "likelihood profiles" in that while one parameter is varied, the others are held fixed, rather than re-optimized. (must also have pplot=T and Iplot=F.)

cplot

Boolean. Plot likelihood contour plot? Provides a likelihood contour plot for a range of r and s values (can be slow so default is FALSE). Must also have lplot=T and pplot=T to get contour plots.

tlab

Character, label for time axis. Defaults to "days".

gfit

Provides a Pearson C type test for goodness of fit. Default is FALSE. Must provide the initial study population to compute goodness of fit.

silent

Optional, Boolean. Stops all print and plot options (still get most warning and all error messages) Default is FALSE. A third option, "verbose" also enables the trace setting in the ms (minimum sum) S-Plus routine.

Value

vector of final MLE r, s, k parameter estimates. standard errors of MLE parameter estimates (if se = <population> is specified).

References

  • Anderson, J.J. (2000). "A vitality-based model relating stressors and environmental properties to organism survival." Ecological Monographs 70(3):445-470.

Examples

Run this code
# NOT RUN {
data(daphnia)
time <- daphnia$days
survival_fraction <- daphnia$lx

results.modk <- vitality.k(time = time,
                           sdata = survival_fraction,
                           rc.data=TRUE, 
                           se=FALSE,
                           gfit=FALSE, 
                           datatype="CUM", 
                           ttol=.000001, 
                           init.params=FALSE,
                           #init.params=c(0.075, 0.15, 0.001),
                           lower=c(0,-1,0), upper=c(100,50,50), 
                           pplot=TRUE, 
                           tlab="days", 
                           lplot=TRUE, 
                           cplot=TRUE, 
                           Iplot=TRUE, 
                           silent=FALSE)
# }

Run the code above in your browser using DataLab