pheno: Generates posterior samples from a phenology model

In rsBayes: Bayesian Models for Remotely Sensed Data

Description

Generates posterior samples from a phenology model using an option of several different likelihoods and prior configurations.

Usage

pheno(formula, data = parent.frame(), gamma = c(0, 1), family = "normal",
      starting, tuning, priors, n.samples, sub.sample,
      summary.only = FALSE, fitted = FALSE, verbose = TRUE, n.report=100, ...)

Arguments

formula	a symbolic description of the regression model to be fit. The left hand size is the name of the Vegetation Index (VI) variable and the right and size is the name of the time variable.
data	an optional data frame containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which pheno is called.
gamma	a vector of length two that holds the VI variable's theoretical bounds. Lower and upper bounds are given in element 1 and 2, respectively. These must be finite values.
family	a quoted string that indicates which likelihood to use for modeling the VI variable. Options are Beta, Normal, and truncated Normal, specified using quoted argument values beta, normal, and t.normal, respectively. If family equals "beta" then The σ^2 parameter is the variance for the Beta, Normal, and truncated Normal likelihoods. See Details for VI variable's support for each likelihood.
starting	a list with each tag corresponding to a parameter name. Valid tags are alpha.1, alpha.2, ..., alpha.7, and sigma.sq. The value portion of each tag is the parameter's starting value. If starting is not provided, then each parameter's starting value is set as a random value from its prior distribution (in most cases this does not provide good starting values). If α starting values are provided, they must be within their respective prior distributions, see the Details for additional guidance on selecting starting values.
tuning	a list with each tag corresponding to a parameter name. Valid tags are alpha.1, alpha.2, ..., alpha.7, and sigma.sq. The value portion of each tag defines the variance of the Metropolis sampler Normal proposal distribution. Tuning values should be selected to keep the acceptance rate between approximately 20 and 50 percent.
priors	a list with tags alpha and sigma.sq.IG. alpha should be a list comprising alpha.1.Unif, alpha.2.Unif, ..., alpha.7.Unif with each of these tags set equal to a vector of length two that holds the lower and upper bound for the Uniform priors on the given alpha. The sigma.sq.IG is a vector of length two that holds the shape and scale parameters for the σ^2's inverse-Gamma prior distribution. At minimum the priors list must include a prior for sigma.sq.IG. If α priors are not specified then they are assumed to follow their default values, see the Details section below.
n.samples	the number of MCMC samples to collect.
sub.sample	an optional list that specifies the subset of MCMC samples should be returned. Valid tags are start, end, and thin. The default values are start=floor(0.5*n.samples), end=n.samples and thin=1.
fitted	if TRUE, regression fitted are returned. The argument sub.sample controls which MCMC samples are used to generate the fitted values.
summary.only	if TRUE, returns the mean, standard deviation, and quantiles of MCMC samples as well as fitted values if fitted=TRUE. The summaries are based on the entire chain or subsamples if sub.sample is provided. The defualt probabilities used in the quantile function are 0.025, 0.5, and 0.975. These default quantiles can be changed, see Details.
verbose	if TRUE, model specification and progress of the sampler is printed to the screen. Otherwise, nothing is printed to the screen.
n.report	the interval to report Metropolis sampler acceptance and MCMC progress.
...	optional argument summary.prob is used to pass a vector of probabilities to the quantile function's probs argument which is used to compute summaries of MCMC samples and fitted values if fitted = TRUE. Optional argument t.normal.bounds takes a vector of length two with the first and second elements defining the truncated Normal's lower and upper bounds, respectively. The default bounds for the truncated Normal are (0,1).

Details

Selection of the likelihood via the family argument should be respective of possible VI value range. The Beta likelihood assumes support for VI between 0 and 1. The Normal likelihood assumes support for VI on the whole real line. The default for the truncated Normal likelihood is support between 0 and 1; however, specifying the optional argument t.normal.bounds allows for user defined support bounds.

The default priors for the α's are:

alpha.1 ~ Unif(gamma[1], gamma[2]),	alpha.2 ~ Unif(0, gamma[2]-alpha.1),	alpha.3 ~ Unif(0, 1)
alpha.4 ~ Unif(0, alpha.7),	alpha.5 ~ Unif(0, 0.001),	alpha.6 ~ Unif(0, 1)
alpha.7 ~ Unif(0, 1)

The hyperparameters of these prior distributions can be change using the priors argument.

An error will be thrown if an alpha starting value provided in the starting list is outside the default or supplied priors. Note, alpha.2 and alpha.4 have upper bounds determined by starting values of alpha.1 and alpha.7, respectively. As noted in the description of the starting argument, starting values will be automatically generated if not provided. We strongly suggest that reasonable starting values are provided or the MCMC sampler might not perform well.

Value

The output is a list of class pheno. List components include:

p.theta.samples	is a coda object of n.samples (or number of subsamples specified in the optional argument sub.samples) with columns labeled alpha.1, alpha.2, ..., alpha.7, and sigma^2.
MH.acceptance	is a named vector with values giving the Metropolis samplers overall acceptance percent and the last batch which is the percent acceptance in the last n.report samples.
p.fitted.samples	if fitted=TRUE then this n, i.e., the number of VI observations, by n.samples (or number of subsamples specified in the optional argument sub.samples) matrix holds the posterior fitted values sampled one-for-one from the parameters posterior values.

Author(s)

Andrew O. Finley [e-mail], Chad Babcock [e-mail]

Examples

  
##select a single pixel and a few years
df <- subset(aoi_hls, pixel == 1 & year %in% c(2017,2018,2019))
plot(df$doy, df$evi)

set.seed(1)

m.beta <- pheno(evi~doy, data=df, family="beta",
                starting=list(alpha.1=0.1, alpha.2=0.5, alpha.3=0.25,
                         alpha.4=50, alpha.5=0.0001, alpha.6 = 0.25, alpha.7=200, sigma.sq=0.01),
		tuning=list(alpha.1=0.001, alpha.2=0.01, alpha.3=0.001,
			 alpha.4=0.1, alpha.5=0.00005, alpha.6=0.001, alpha.7=0.5, sigma.sq=0.1),
                n.samples=50000, n.report=25000,
                priors=list(alpha=list(alpha.5=c(-1, 1)), sigma.sq.IG=c(2, 0.01)),
                fitted = TRUE, sub.sample=list("start"=25000)
                )

summary(m.beta)

rsBayes documentation built on Sept. 14, 2020, 5:06 p.m.