Generates posterior samples from a phenology model using an option of several different likelihoods and prior configurations.
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, ...)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.
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.
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.
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 \(\sigma^2\) parameter
is the variance for the Beta, Normal, and truncated Normal likelihoods. See Details for VI variable's support for each likelihood.
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
\(\alpha\) starting values are provided, they must be within their
respective prior distributions, see the Details for additional
guidance on selecting starting values.
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.
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
\(\sigma^2\)'s inverse-Gamma prior distribution. At minimum the priors list
must include a prior for sigma.sq.IG. If \(\alpha\) priors are not specified then they are assumed to follow their default values, see the Details section below.
the number of MCMC samples to collect.
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.
if TRUE, regression fitted are returned. The
argument sub.sample controls which MCMC samples are used to
generate the fitted values.
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.
if TRUE, model specification and progress of the
sampler is printed to the screen. Otherwise, nothing is printed to the screen.
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).
The output is a list of class pheno. List components include:
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.
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.
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.
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 \(\alpha\)'s are:
| \(\alpha_1 \sim Unif(\gamma_1, \gamma_2)\), | \(\alpha_2 \sim Unif(0, \gamma_2-\alpha_1)\), | \(\alpha_3 \sim Unif(0, 1)\) |
| \(\alpha_4 \sim Unif(0, \alpha_7)\), | \(\alpha_5 \sim Unif(0, 0.001)\), | \(\alpha_6 \sim Unif(0, 1)\) |
| \(\alpha_7 \sim Unif(1, 365)\) |
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.
# NOT RUN {
##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)
# }
Run the code above in your browser using DataLab