fit_phenology: Fit the phenology parameters to timeseries of counts.
In phenology: Tools to Manage a Parametric Function that Describes Phenology and More

fit_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

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)

phenology documentation built on Oct. 16, 2023, 9:06 a.m.