| Title: | Fit Density Curves to Peak Timing Data that Varies over Time |
|---|---|
| Description: | The 'phenomix' package fits time-varying density curves to run timing type data commonly encountered in fisheries and ecology. Example applications include to peak run timing curves collected for juvenile or adult Pacific salmon, though could also be applied to other kinds of data such as hydrographs, plant phenology (flowering, leaf out). |
| Authors: | Eric J. Ward [aut, cre], Samantha M. Wilson [ctb], Joseph H. Anderson [ctb] |
| Maintainer: | Eric J. Ward <[email protected]> |
| License: | GPL (>=3) |
| Version: | 1.0.4 |
| Built: | 2024-11-26 16:31:15 UTC |
| Source: | https://github.com/nwfsc-cb/phenomix |
Does minimal processing of data to use as argument to fitting function
create_data( data, min_number = 0, variable = "number", time = "year", date = "doy", asymmetric_model = TRUE, mu = ~1, sigma = ~1, covar_data = NULL, est_sigma_re = TRUE, est_mu_re = TRUE, tail_model = "student_t", family = "lognormal", max_theta = 10, share_shape = TRUE, nu_prior = c(2, 10), beta_prior = c(2, 1) )create_data( data, min_number = 0, variable = "number", time = "year", date = "doy", asymmetric_model = TRUE, mu = ~1, sigma = ~1, covar_data = NULL, est_sigma_re = TRUE, est_mu_re = TRUE, tail_model = "student_t", family = "lognormal", max_theta = 10, share_shape = TRUE, nu_prior = c(2, 10), beta_prior = c(2, 1) )
data |
A data frame |
min_number |
A minimum threshold to use, defaults to 0 |
variable |
A character string of the name of the variable in 'data' that contains the response (e.g. counts) |
time |
A character string of the name of the variable in 'data' that contains the time variable (e.g. year) |
date |
A character string of the name of the variable in 'data' that contains the response (e.g. day of year). The actual #' column should contain a numeric response – for example, the result from using lubridate::yday(x) |
asymmetric_model |
Boolean, whether or not to let model be asymmetric (e.g. run timing before peak has a different shape than run timing after peak) |
mu |
An optional formula allowing the mean to be a function of covariates. Random effects are not included in the formula
but specified with the |
sigma |
An optional formula allowing the standard deviation to be a function of covariates. For asymmetric models,
each side of the distribution is allowed a different set of covariates. Random effects are not included in the formula
but specified with the |
covar_data |
a data frame containing covariates specific to each time step. These are used in the formulas |
est_sigma_re |
Whether to estimate random effects by year in sigma parameter controlling tail of distribution. Defaults to TRUE |
est_mu_re |
Whether to estimate random effects by year in mu parameter controlling location of distribution. Defaults to TRUE |
tail_model |
Whether to fit Gaussian ("gaussian"), Student-t ("student_t") or generalized normal ("gnorm"). Defaults to Student-t |
family |
Response for observation model, options are "gaussian", "poisson", "negbin", "binomial", "lognormal". The default ("lognormal") is not a true lognormal distribution, but a normal-log in that it assumes log(y) ~ Normal() |
max_theta |
Maximum value of log(pred) when |
share_shape |
Boolean argument for whether asymmetric student-t and generalized normal distributions should share the shape parameter (nu for the student-t; beta for the generalized normal). Defaults to TRUE |
nu_prior |
Two element vector (optional) for penalized prior on student t df, defaults to a Gamma(shape=2, scale=10) distribution |
beta_prior |
Two element vector (optional) for penalized prior on generalized normal beta, defaults to a Normal(2, 1) distribution |
data(fishdist) datalist <- create_data(fishdist, min_number = 0, variable = "number", time = "year", date = "doy", asymmetric_model = TRUE, family = "gaussian" )data(fishdist) datalist <- create_data(fishdist, min_number = 0, variable = "number", time = "year", date = "doy", asymmetric_model = TRUE, family = "gaussian" )
This function extracts the means, sigmas, thetas, lower (25%) and upper (75%) quartiles, and respective sds
extract_all(fit)extract_all(fit)
fit |
A fitted object returned from fit() |
This function extracts the annual totals
extract_annual(fit, log = TRUE)extract_annual(fit, log = TRUE)
fit |
A fitted object returned from fit() |
log |
Whether to return estimates in log space, defaults to TRUE |
This function extracts the lower quartiles (25%) and respective sds
extract_lower(fit)extract_lower(fit)
fit |
A fitted object returned from fit() |
This function extracts the parameter means and respective sds
extract_means(fit)extract_means(fit)
fit |
A fitted object returned from fit() |
This function extracts the parameter sigma and respective sds
extract_sigma(fit)extract_sigma(fit)
fit |
A fitted object returned from fit() |
This function extracts the parameter theta and respective sds
extract_theta(fit)extract_theta(fit)
fit |
A fitted object returned from fit() |
This function extracts the upper quartiles (25%) and respective sds
extract_upper(fit)extract_upper(fit)
fit |
A fitted object returned from fit() |
This function creates a list of parameters, sets up TMB object and attempts to do fitting / estimation
fit( data_list, silent = FALSE, inits = NULL, control = list(eval.max = 2000, iter.max = 1000, rel.tol = 1e-10), limits = NULL, fit_model = TRUE )fit( data_list, silent = FALSE, inits = NULL, control = list(eval.max = 2000, iter.max = 1000, rel.tol = 1e-10), limits = NULL, fit_model = TRUE )
data_list |
A list of data, as output from create_data |
silent |
Boolean passed to TMB::MakeADFun, whether to be verbose or not (defaults to FALSE) |
inits |
Optional named list of parameters for starting values, defaults to NULL |
control |
Optional control list for stats::nlminb. For arguments see ?nlminb. Defaults to eval.max=2000, iter.max=1000, rel.tol=1e-10. For final model runs, the rel.tol should be even smaller |
limits |
Whether to include limits for stats::nlminb. Can be a list of (lower, upper), or TRUE to use suggested hardcoded limits. Defaults to NULL, where no limits used. |
fit_model |
Whether to fit the model. If not, returns a list including the data, parameters, and initial values. Defaults to TRUE |
data(fishdist) # example of fitting fixed effects, no trends, no random effects set.seed(1) datalist <- create_data(fishdist[which(fishdist$year > 1970), ], asymmetric_model = FALSE, est_mu_re = FALSE, est_sigma_re = FALSE ) fit <- fit(datalist) # # # example of model with random effects in means only, and symmetric distribution # set.seed(1) # datalist <- create_data(fishdist[which(fishdist$year > 1970),], asymmetric_model = FALSE, # est_sigma_re = FALSE) # fit <- fit(datalist) # # example of model with random effects in variances # set.seed(1) # datalist <- create_data(fishdist[which(fishdist$year > 1970),], asymmetric_model = TRUE, # est_mu_re = TRUE) # fit <- fit(datalist) # # # example of model with poisson response # set.seed(1) # datalist <- create_data(fishdist[which(fishdist$year > 1970),], asymmetric_model = FALSE, # est_sigma_trend=FALSE, est_mu_trend=FALSE, est_mu_re = TRUE, # family="poisson") # fit <- fit(datalist)data(fishdist) # example of fitting fixed effects, no trends, no random effects set.seed(1) datalist <- create_data(fishdist[which(fishdist$year > 1970), ], asymmetric_model = FALSE, est_mu_re = FALSE, est_sigma_re = FALSE ) fit <- fit(datalist) # # # example of model with random effects in means only, and symmetric distribution # set.seed(1) # datalist <- create_data(fishdist[which(fishdist$year > 1970),], asymmetric_model = FALSE, # est_sigma_re = FALSE) # fit <- fit(datalist) # # example of model with random effects in variances # set.seed(1) # datalist <- create_data(fishdist[which(fishdist$year > 1970),], asymmetric_model = TRUE, # est_mu_re = TRUE) # fit <- fit(datalist) # # # example of model with poisson response # set.seed(1) # datalist <- create_data(fishdist[which(fishdist$year > 1970),], asymmetric_model = FALSE, # est_sigma_trend=FALSE, est_mu_trend=FALSE, est_mu_re = TRUE, # family="poisson") # fit <- fit(datalist)
Get fixed effects parameters from model object, copying glmmTMB 'fast' implementation
## S3 method for class 'phenomix' fixef(object, ...)## S3 method for class 'phenomix' fixef(object, ...)
object |
The fitted phenomix model |
... |
Additional arguments |
Arguments can be adjusted as needed
limits(parnames, max_theta)limits(parnames, max_theta)
parnames |
A vector of character strings or names used for estimation |
max_theta |
A scalar (optional) giving the maximum value of theta, log(pred) |
These functions make some basic plots for the user
plot_diagnostics(fitted, type = "timing", logspace = TRUE)plot_diagnostics(fitted, type = "timing", logspace = TRUE)
fitted |
A fitted model object |
type |
A plot type for ggplot, either "timing" or "scatter" |
logspace |
whether to plot the space in log space, defaults to TRUE |
Get random effects parameters from model object, copying glmmTMB 'fast' implementation
## S3 method for class 'phenomix' ranef(object, ...)## S3 method for class 'phenomix' ranef(object, ...)
object |
The fitted phenomix model |
... |
Additional arguments |