# Fit the phenology parameters to timeseries of counts. 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) ```r 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. ## Returns Return a list of with data and result ## Details fit_phenology fits parameters to timeseries. ## Examples ```r ## 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, '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) ## End(Not run) ``` ## 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()` ## Author(s) Marc Girondot marc.girondot@gmail.com

fit_phenology function