| fit_phenology {phenology} | R Documentation |
Fit the phenology parameters to timeseries of counts.
Description
Function of the package phenology to fit parameters to timeseries.
To fit data, the syntax is :
Result <- fit_phenology(data=dataset, fitted.parameters=par, fixed.parameters=pfixed, trace=1, hessian=TRUE)
or if no parameter is fixed :
Result <- fit_phenology(data=dataset, fitted.parameters=par)
Add trace=1 [default] to have information on the fit progression or trace=0 to hide information on the fit progression.
hessian = FALSE does not estimate Hessian matrix and SE of parameters.
If the parameter Theta is fixed to +Inf, a Poissonian model of daily nest distribution is implemented.
Special section about cofactors:
cofactors must be a data.frame with a column Date and a column for each cofactor
add.cofactors are the names of the column of parameter cofactors to use as a cofactor;
The model is then: parameter[add.cofactors] * cofactor[, add.cofactors]
If the name of the parameter is paste0(add.cofactors, "multi"), then the model is:
parameter[paste0(add.cofactors, "multi")] * cofactor[, add.cofactors] *
(number of nests without cofactor)
About parallel computing:
Set options mc.cores and forking to tell what sort of parallel computing
Example:
options(mc.cores = detectCores())
options(forking = FALSE)
Usage
fit_phenology(
data = file.choose(),
fitted.parameters = NULL,
fixed.parameters = NULL,
model_before = NULL,
store.intermediate = FALSE,
file.intermediate = "Intermediate.rda",
hessian = FALSE,
silent = FALSE,
cofactors = NULL,
add.cofactors = NULL,
zero = 1e-09,
lower = 0,
upper = Inf,
stop.fit = FALSE,
method_Snbinom = "saddlepoint",
control = list(trace = 1, REPORT = 1, maxit = 1000),
method = c("Nelder-Mead", "L-BFGS-B")
)
Arguments
data |
A dataset generated by add_format |
fitted.parameters |
Set of parameters to be fitted |
fixed.parameters |
Set of fixed parameters |
model_before |
The change of parameters before to estimate daily counts. |
store.intermediate |
TRUE or FALSE to save the intermediates |
file.intermediate |
Name of the file where to save the intermediates as a list |
hessian |
If FALSE does not estimate se of parameters |
silent |
If TRUE does not show any message |
cofactors |
data.frame with a column Date and a column for each cofactor |
add.cofactors |
Names of the column of parameter cofactors to use as a cofactor |
zero |
If the theoretical nest number is under this value, this value will be used |
lower |
Lower bound for each parameter |
upper |
Upper bound for each parameter |
stop.fit |
If TRUE, will stop search for parameters even if not ML |
method_Snbinom |
Can be Furman, exact, or saddlepoint. |
control |
List for control parameters for optim |
method |
Method used by optim. Several can be setup. |
Details
fit_phenology fits parameters to timeseries.
Value
Return a list of with data and result
Author(s)
Marc Girondot marc.girondot@gmail.com
See Also
Other Phenology model:
AutoFitPhenology(),
BE_to_LBLE(),
Gratiot,
LBLE_to_BE(),
LBLE_to_L(),
L_to_LBLE(),
MarineTurtles_2002,
MinBMinE_to_Min(),
adapt_parameters(),
add_SE(),
add_phenology(),
extract_result(),
likelihood_phenology(),
logLik.phenology(),
map_Gratiot,
map_phenology(),
par_init(),
phenology2fitRMU(),
phenology_MHmcmc_p(),
phenology_MHmcmc(),
phenology(),
plot.phenologymap(),
plot.phenology(),
plot_delta(),
plot_phi(),
print.phenologymap(),
print.phenologyout(),
print.phenology(),
remove_site(),
result_Gratiot1,
result_Gratiot2,
result_Gratiot_Flat,
result_Gratiot_mcmc,
result_Gratiot,
summary.phenologymap(),
summary.phenologyout(),
summary.phenology()
Examples
## Not run:
library(phenology)
# Read a file with data
data(Gratiot)
# Generate a formatted list nammed data_Gratiot
data_Gratiot <- add_phenology(Gratiot, name="Complete",
reference=as.Date("2001-01-01"), format="%d/%m/%Y")
# Generate initial points for the optimisation
parg <- par_init(data_Gratiot, fixed.parameters=NULL)
# Run the optimisation
result_Gratiot <- fit_phenology(data=data_Gratiot,
fitted.parameters=parg, fixed.parameters=NULL)
data(result_Gratiot)
# Plot the phenology and get some stats
output <- plot(result_Gratiot)
# or
output <- summary(result_Gratiot)
# With one sinusoid
parg <- c('LengthB' = 95.187648986054285,
'Peak' = 173.86111755419921,
'LengthE' = 63.183281230994481,
'Max_Complete' = 33.052091519419136,
'MinB_Complete' = 0.21716081973776738,
'MinE_Complete' = 0.42444245288475702,
'Theta' = 3.9554976911657187,
'Alpha' = 0,
'Beta' = 0.21019683898423902,
'Delta' = 3.6798422076201724,
result_Gratiot1 <- fit_phenology(data=data_Gratiot,
fitted.parameters=parg, fixed.parameters=NULL)
plot(result_Gratiot1)
# With two sinusoids
parg <- c('LengthB' = 95.821859173220659,
'Peak' = 174.89756503758881,
'LengthE' = 60.954167489825352,
'Max_Complete' = 33.497846416498774,
'MinB_Complete' = 0.21331670808871037,
'MinE_Complete' = 0.43730327416828613,
'Theta' = 4.2569203480842814,
'Alpha1' = 0.15305562092653918,
'Beta1' = 0,
'Delta1' = 9.435730952200263,
'Phi1' = 15.7944494669335,
'Alpha' = 0,
'Beta' = 0.28212926001402688,
'Delta' = 18.651087311957518,
'Phi' = 9.7549929595313056)
result_Gratiot2 <- fit_phenology(data=data_Gratiot,
fitted.parameters=parg, fixed.parameters=NULL)
plot(result_Gratiot2)
compare_AICc(no=result_Gratiot,
one=result_Gratiot1,
two=result_Gratiot2)
# With parametrization based on Girondot 2010
parg <- c('Peak' = 173.52272236775076,
'Flat' = 0,
'LengthB' = 94.433284359804205,
'LengthE' = 64.288485646867329,
'Max_Complete' = 32.841568389778033,
'PMin' = 1.0368261725650889,
'Theta' = 3.5534818927979592)
result_Gratiot_par1 <- fit_phenology(data=data_Gratiot,
fitted.parameters=parg, fixed.parameters=NULL)
# With new parametrization based on Omeyer et al. (2022):
# Omeyer, L. C. M., McKinley, T. J., Bréheret, N., Bal, G., Balchin, G. P.,
# Bitsindou, A., Chauvet, E., Collins, T., Curran, B. K., Formia, A., Girard, A.,
# Girondot, M., Godley, B. J., Mavoungou, J.-G., Poli, L., Tilley, D.,
# VanLeeuwe, H. & Metcalfe, K. 2022. Missing data in sea turtle population
# monitoring: a Bayesian statistical framework accounting for incomplete
# sampling Front. Mar. Sci. (IF 3.661), 9, 817014.
parg <- result_Gratiot_par1$par
parg <- c(tp=unname(parg["Peak"]), tf=unname(parg["Flat"]),
s1=unname(parg["LengthB"])/4.8, s2=unname(parg["LengthE"])/4.8,
alpha=unname(parg["Max_Complete"]), Theta=0.66)
result_Gratiot_par2 <- fit_phenology(data=data_Gratiot,
fitted.parameters=parg,
fixed.parameters=NULL)
compare_AICc(GirondotModel=result_Gratiot_par1,
OmeyerModel=result_Gratiot_par2)
plot(result_Gratiot_par1)
plot(result_Gratiot_par2)
# Use fit with co-factor
# First extract tide information for that place
td <- tide.info(year=2001, latitude=4.9167, longitude=-52.3333)
# I keep only High tide level
td2 <- td[td$Phase=="High Tide", ]
# I get the date
td3 <- cbind(td2, Date=as.Date(td2$DateTime.local))
td5 <- aggregate(x=td3[, c("Date", "DateTime.local", "Tide.meter")],
by=list(Date=td3[, "Date"]), FUN=max)[, 2:4]
with(td5, plot(DateTime.local, Tide.meter, type="l"))
td6 <- td5[, c("Date", "Tide.meter")]
parg <- par_init(data_Gratiot, fixed.parameters=NULL,
add.cofactors="Tide.meter")
likelihood_phenology(data=data_Gratiot, fitted.parameters = parg,
cofactors=td6, add.cofactors="Tide.meter")
result_Gratiot_CF <- fit_phenology(data=data_Gratiot,
fitted.parameters=parg, fixed.parameters=NULL, cofactors=td6,
add.cofactors="Tide.meter")
compare_AIC(WithoutCF=result_Gratiot, WithCF=result_Gratiot_CF)
plot(result_Gratiot_CF)
# Example with two series fitted with different peaks but same Length of season
Gratiot2 <- Gratiot
Gratiot2[, 2] <- floor(Gratiot2[, 2]*runif(n=nrow(Gratiot2)))
data_Gratiot <- add_phenology(Gratiot, name="Complete1",
reference=as.Date("2001-01-01"), format="%d/%m/%Y")
data_Gratiot <- add_phenology(Gratiot2, name="Complete2",
reference=as.Date("2001-01-01"),
format="%d/%m/%Y", previous=data_Gratiot)
pfixed=c(Min=0)
p <- par_init(data_Gratiot, fixed.parameters = pfixed)
p <- c(p, Peak_Complete1=175, Peak_Complete2=175)
p <- p[-4]
p <- c(p, Length=90)
p <- p[-(3:4)]
result_Gratiot <- fit_phenology(data=data_Gratiot, fitted.parameters=p,
fixed.parameters=pfixed)
# An example with bimodality
g <- Gratiot
g[30:60, 2] <- sample(10:20, 31, replace = TRUE)
data_g <- add_phenology(g, name="Complete", reference=as.Date("2001-01-01"),
format="%d/%m/%Y")
parg <- c('Max.1_Complete' = 5.6344636692341856,
'MinB.1_Complete' = 0.15488810581002324,
'MinE.1_Complete' = 0.2,
'LengthB.1' = 22.366647176407636,
'Peak.1' = 47.902473939250036,
'LengthE.1' = 17.828495918533015,
'Max.2_Complete' = 33.053364083447434,
'MinE.2_Complete' = 0.42438173496989717,
'LengthB.2' = 96.651564706802702,
'Peak.2' = 175.3451874571835,
'LengthE.2' = 62.481968743789835,
'Theta' = 3.6423908093342572)
pfixed <- c('MinB.2_Complete' = 0,
'Flat.1' = 0,
'Flat.2' = 0)
result_g <- fit_phenology(data=data_g, fitted.parameters=parg, fixed.parameters=pfixed)
plot(result_g)
# Exemple with some minimum counts
nb <- Gratiot[, 2]
nbs <- sample(0:1, length(nb), replace=TRUE)
nb <- ifelse(nb != 0, ifelse(nbs == 1, 1, nb), 0)
nbc <- ifelse(nb != 0, ifelse(nbs == 1, "minimum", "exact"), "exact")
Gratiot_minimal <- cbind(Gratiot, CountTypes=nbc)
Gratiot_minimal[, 2] <- nb
data_Gratiot_minimal <- add_phenology(add=Gratiot_minimal,
colname.CountTypes = "CountTypes",
month_ref=1)
parg <- par_init(data_Gratiot_minimal, fixed.parameters=NULL)
result_Gratiot_minimal <- fit_phenology(data=data_Gratiot_minimal,
fitted.parameters=parg, fixed.parameters=NULL)
plot(result_Gratiot_minimal)
summary(result_Gratiot_minimal)
# Exemple with all zero counts being not recorded
Gratiot_NoZeroCounts <- cbind(Gratiot[Gratiot[, 2] != 0, ], ZeroCounts=FALSE)
data_Gratiot_NoZeroCounts <- add_phenology(add=Gratiot_NoZeroCounts,
colname.ZeroCounts = "ZeroCounts",
ZeroCounts.defaul=FALSE,
month_ref=1)
# or
data_Gratiot_NoZeroCounts <- add_phenology(add=Gratiot[Gratiot[, 2] != 0, ],
ZeroCounts.default=FALSE,
month_ref=1)
parg <- par_init(data_Gratiot_NoZeroCounts, fixed.parameters=NULL)
result_Gratiot_NoZeroCounts <- fit_phenology(data=data_Gratiot_NoZeroCounts,
fitted.parameters=parg, fixed.parameters=NULL)
plot(result_Gratiot_NoZeroCounts)
summary(result_Gratiot_NoZeroCounts)
# Exemple with data in range of date
Gratiot_rangedate <- Gratiot
Gratiot_rangedate[, 1] <- as.character(Gratiot_rangedate[, 1])
Gratiot_rangedate[148, 1] <- paste0(Gratiot_rangedate[148, 1], "-", Gratiot_rangedate[157, 1])
Gratiot_rangedate[148, 2] <- sum(Gratiot_rangedate[148:157, 2])
Gratiot_rangedate <- Gratiot_rangedate[-(149:157), ]
data_Gratiot_rangedate <- add_phenology(add=Gratiot_rangedate,
month_ref=1)
parg <- par_init(data_Gratiot_rangedate, fixed.parameters=NULL)
likelihood_phenology(data=data_Gratiot_rangedate,
fitted.parameters=parg)
result_Gratiot_rangedate <- fit_phenology(data=data_Gratiot_rangedate,
fitted.parameters=parg,
fixed.parameters=NULL)
plot(result_Gratiot_rangedate)
likelihood_phenology(result=result_Gratiot_rangedate)
# Exemple with data in range of date and CountTypes being minimum
Gratiot_rangedate <- Gratiot
Gratiot_rangedate[, 1] <- as.character(Gratiot_rangedate[, 1])
Gratiot_rangedate[148, 1] <- paste0(Gratiot_rangedate[148, 1], "-", Gratiot_rangedate[157, 1])
Gratiot_rangedate[148, 2] <- sum(Gratiot_rangedate[148:157, 2])
Gratiot_rangedate <- Gratiot_rangedate[-(149:157), ]
Gratiot_rangedate <- cbind(Gratiot_rangedate, CountTypes="exact")
Gratiot_rangedate[148, 2] <- 100
Gratiot_rangedate[148, "CountTypes"] <- "minimum"
Gratiot_rangedate[28, "CountTypes"] <- "minimum"
data_Gratiot_rangedate <- add_phenology(add=Gratiot_rangedate,
colname.CountTypes="CountTypes",
month_ref=1)
parg <- par_init(data_Gratiot_rangedate, fixed.parameters=NULL)
result_Gratiot_rangedate <- fit_phenology(data=data_Gratiot_rangedate,
fitted.parameters=parg, fixed.parameters=NULL)
likelihood_phenology(result=result_Gratiot_rangedate)
plot(result_Gratiot_rangedate)
## End(Not run)