Learn R Programming
phenology (version 10.3)

fit_phenology: 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")
)

Value

Return a list of with data and result

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.

Author

Marc Girondot marc.girondot@gmail.com

Details

fit_phenology fits parameters to timeseries.

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(), phenology(), phenology2fitRMU(), phenology_MHmcmc(), phenology_MHmcmc_p(), plot.phenology(), plot.phenologymap(), plot_delta(), plot_phi(), print.phenology(), print.phenologymap(), print.phenologyout(), remove_site(), result_Gratiot, result_Gratiot1, result_Gratiot2, result_Gratiot_Flat, summary.phenology(), summary.phenologymap(), summary.phenologyout()

Examples

Run this code
if (FALSE) {
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,
           'Phi' = 10)
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)
		
}

Run the code above in your browser using DataLab