Nothing
R/fit_phenology.R
In phenology: Tools to Manage a Parametric Function that Describes Phenology and More
Defines functions
fit_phenology
Documented in
fit_phenology
#' fit_phenology fits parameters to timeseries.
#' @title Fit the phenology parameters to timeseries of counts.
#' @author Marc Girondot \email{marc.girondot@@gmail.com}
#' @return Return a list of with data and result
#' @param data A dataset generated by add_format
#' @param fixed.parameters Set of fixed parameters
#' @param fitted.parameters Set of parameters to be fitted
#' @param lower Lower bound for each parameter
#' @param upper Upper bound for each parameter
#' @param hessian If FALSE does not estimate se of parameters
#' @param model_before The change of parameters before to estimate daily counts.
#' @param cofactors data.frame with a column Date and a column for each cofactor
#' @param add.cofactors Names of the column of parameter cofactors to use as a cofactor
#' @param silent If TRUE does not show any message
#' @param zero If the theoretical nest number is under this value, this value will be used
#' @param stop.fit If TRUE, will stop search for parameters even if not ML
#' @param control List for control parameters for optim
#' @param method Method used by optim. Several can be setup.
#' @param method_Snbinom Can be Furman, exact, or saddlepoint.
#' @param store.intermediate TRUE or FALSE to save the intermediates
#' @param file.intermediate Name of the file where to save the intermediates as a list
#' @description Function of the package phenology to fit parameters to timeseries.\cr
#' To fit data, the syntax is :\cr
#' Result <- fit_phenology(data=dataset, fitted.parameters=par, fixed.parameters=pfixed, trace=1, hessian=TRUE)\cr
#' or if no parameter is fixed :\cr
#' Result <- fit_phenology(data=dataset, fitted.parameters=par)\cr
#' Add trace=1 [default] to have information on the fit progression or trace=0 to hide information on the fit progression.\cr
#' hessian = FALSE does not estimate Hessian matrix and SE of parameters.\cr
#' If the parameter Theta is fixed to +Inf, a Poissonian model of daily nest distribution is implemented.\cr
#' Special section about cofactors:\cr
#' cofactors must be a data.frame with a column Date and a column for each cofactor\cr
#' add.cofactors are the names of the column of parameter cofactors to use as a cofactor;\cr
#' The model is then: parameter[add.cofactors] * cofactor[, add.cofactors]\cr
#' If the name of the parameter is paste0(add.cofactors, "multi"), then the model is:\cr
#' parameter[paste0(add.cofactors, "multi")] * cofactor[, add.cofactors] *
#' (number of nests without cofactor)\cr
#' About parallel computing:\cr
#' Set options mc.cores and forking to tell what sort of parallel computing\cr
#' Example:\cr
#' options(mc.cores = detectCores())\cr
#' options(forking = FALSE)
#' @family Phenology model
#' @examples
#' \dontrun{
#' 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)
#'
#' }
#' @export
fit_phenology <-
function(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-9,
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")) {
# data=NULL
# lower = 0
# upper = Inf
# fitted.parameters=NULL
# fixed.parameters=NA
# cofactors=NULL
# add.cofactors=NULL
# hessian=TRUE
# silent=FALSE
# zero=1E-9
# control=list(trace=1, REPORT=1, maxit=1000)
# method = c("Nelder-Mead", "L-BFGS-B")
# stop.fit=FALSE
# store.intermediate=FALSE
# file.intermediate="Intermediate.rda"
# data_Gratiot <- add_phenology(Gratiot, name="Complete", reference=as.Date("2001-01-01"), format="%d/%m/%Y")
# data=data_Gratiot; fitted.parameters=result_Gratiot$par; fixed.parameters=result_Gratiot$fixed.parameters; trace=1
# if (is.null(fixed.parameters)) {fixed.parameters<-NA}
# if (!requireNamespace("optimx", quietly = TRUE)) {
# stop("optimx package is absent; Please install it first")
# }
# method <- c("Nelder-Mead", "L-BFGS-B")
# if (length(fitted.parameters) == 1) method <- "Brent"
# library("optimx")
Lnegbin <- getFromNamespace(".Lnegbin", ns="phenology")
if (store.intermediate) {
store <- list()
save(store, file=file.intermediate)
}
if (!inherits(data, "phenologydata")) {
message("Data should have been formated first using the function add_phenology(). I format it now.")
data <- add_phenology(data)
}
if (is.null(fitted.parameters)) {
message("No initial parameters set has been defined. I estimate one set using par_init().")
fitted.parameters <- par_init(data, fixed.parameters=fixed.parameters)
}
if ((!is.null(add.cofactors)) & (!is.null(cofactors))) {
cf1 <- cofactors
for (ff in add.cofactors)
cf1[, ff] <- cofactors[, ff] - mean(cofactors[, ff])
cofactors <- cf1
}
preLNL <- +Inf
if (control[["maxit"]] != 0) {
repeat {
for (m in seq_along(method)) {
if ((method[m] == "L-BFGS-B") | (method[m] == "Brent")) {
resul <- optim(par=fitted.parameters, fn=Lnegbin,
pt=list(data=data, fixed=fixed.parameters,
model_before=model_before,
out=TRUE, cofactors=cofactors,
add.cofactors=add.cofactors, zero=zero,
namespar=names(fitted.parameters),
method_Snbinom=method_Snbinom,
store.intermediate=store.intermediate,
file.intermediate=file.intermediate),
lower = lower,
upper = upper,
method = method[m],
control=control,
hessian=FALSE)
} else {
resul <- optim(par=fitted.parameters, fn=Lnegbin,
pt=list(data=data, fixed=fixed.parameters,
model_before=model_before,
out=TRUE, cofactors=cofactors,
add.cofactors=add.cofactors, zero=zero,
namespar=names(fitted.parameters),
method_Snbinom=method_Snbinom,
store.intermediate=store.intermediate,
file.intermediate=file.intermediate),
method = method[m],
control=control,
hessian=FALSE)
}
resfit <- resul$par
if (is.null(names(resfit))) {
names(resfit) <- names(fitted.parameters)
}
if (any(grepl("^Peak", names(resfit)))) resfit[substr(names(resfit), 1, 4)=="Peak"] <- abs(resfit[substr(names(resfit), 1, 4)=="Peak"])
if (any(grepl("^Theta", names(resfit)))) resfit[substr(names(resfit), 1, 5)=="Theta"] <- abs(resfit[substr(names(resfit), 1, 5)=="Theta"])
if (any(grepl("^PMin", names(resfit)))) resfit["PMin"] <- abs(resfit["PMin"])
# 28/3/2018. J'avais oublié ceux là
if (any(grepl("^PMinE", names(resfit)))) resfit["PMinE"] <- abs(resfit["PMinE"])
if (any(grepl("^PMinB", names(resfit)))) resfit["PMinB"] <- abs(resfit["PMinB"])
if (any(grepl("^Flat", names(resfit)))) resfit["Flat"] <- abs(resfit["Flat"])
if (any(grepl("^Length", names(resfit)))) resfit[substr(names(resfit), 1, 6)=="Length"] <- abs(resfit[substr(names(resfit), 1, 6)=="Length"])
# Ca fait en même temps MinE et MinB
if (any(grepl("^Min", names(resfit)))) resfit[substr(names(resfit), 1, 3)=="Min"]<-abs(resfit[substr(names(resfit), 1, 3)=="Min"])
if (any(grepl("^PMax", names(resfit)))) resfit[substr(names(resfit), 1, 3)=="Max"]<-abs(resfit[substr(names(resfit), 1, 3)=="Max"])
resul$par <- resfit
fitted.parameters <- resfit
# 22/8/2022
fake <- d(fitted.parameters)
}
conv <- resul$convergence
value <- resul$value
if ((conv == 0) | (abs(preLNL - value)<1e-6) | (stop.fit)) break
if (!silent) message("Convergence is not achieved. Optimization continues !")
preLNL <- value
}
if ((conv != 0) & (!stop.fit)) message("Convergence is not achieved but change in likelihood is <1e-6. Use the results with caution.")
if ((conv != 0) & (stop.fit)) message("Convergence is not achieved but maximum number of iterations is reached. Use the results with caution.")
if (!silent) {
message("Fit done!")
if (is.finite(value)) {
cat(paste("-Ln L=", format(value, digits=max(3, trunc(log10(abs(value)))+4)), "\n", sep=""))
} else {
cat("-Ln L= -Inf\n")
}
}
result_list <- resul
# result_list <- list()
# result_list$par <- resfit
# result_list$value <- value
# result_list$convergence <- conv
# resul <- result_list
} else {
# J'avais maxit == 0
resul <- list()
resul$data <- data
resul$par <- fitted.parameters
resul$value <- NA
resul$convergence <- NA
resfit <- fitted.parameters
}
if (hessian) {
if (!requireNamespace("numDeriv", quietly = TRUE)) {
stop("numDeriv package is absent; Please install it first")
}
if (!silent) message("Estimation of the standard error of parameters. Be patient please.")
mathessian <- try(getFromNamespace("hessian", ns="numDeriv")(func=Lnegbin,
x=resfit,
method="Richardson",
pt=list(data=data, fixed=fixed.parameters,
model_before=model_before,
out=TRUE, cofactors=cofactors,
add.cofactors=add.cofactors, zero=zero,
namespar=names(fitted.parameters),
method_Snbinom=method_Snbinom,
store.intermediate=store.intermediate,
file.intermediate=file.intermediate))
, silent=TRUE)
if (substr(mathessian[1], 1, 5)=="Error") {
if (!silent) {
warning("Error in the fit; probably one or more parameters are not estimable.")
warning("Standard errors of parameters cannot be estimated.")
}
res_se <- rep(NA, length(resfit))
names(res_se) <- names(resfit)
} else {
rownames(mathessian) <- colnames(mathessian) <- names(resfit)
resul$hessian <- mathessian
res_se <- SEfromHessian(mathessian, silent=silent)
}
} else {
if (!silent) warning("Standard errors are not estimated.")
res_se <- rep(NA, length(resfit))
names(res_se) <- names(resfit)
}
resul$se <- res_se
resul$fixed.parameters <- fixed.parameters
resul$method_incertitude <- "convolution"
resul$method_Snbinom <- method_Snbinom
resul$data <- data
# Lnegbin(x=resfit, pt=list(data=data, fixed=fixed.parameters,
# out=FALSE, cofactors=cofactors,
# method_Snbinom=method_Snbinom,
# add.cofactors=add.cofactors, zero=zero))
for(kl in 1:length(res_se)) {
if (is.na(res_se[kl])) {
if (!silent) cat(paste(names(resfit[kl]), "=", format(resfit[kl], digits=max(3, trunc(log10(abs(resfit[kl])))+4)), " SE= NaN\n", sep=""))
} else {
if (!silent) cat(paste(names(resfit[kl]), "=", format(resfit[kl], digits=max(3, trunc(log10(abs(resfit[kl])))+4)), " SE=", format(res_se[kl], digits=max(3, trunc(log10(abs(res_se[kl])))+4)), "\n", sep=""))
}
}
dtout <- list()
for (kl in 1:length(resul$data)) {
if (!silent) cat(paste("Series: ", names(resul$data[kl]), "\n", sep=""))
# la date de référence est resul$data[[kl]][1, "Date"]-resul$data[[kl]][1, "ordinal"]
ref <- resul$data[[kl]][1, "Date"]-resul$data[[kl]][1, "ordinal"]
intdtout <- c(reference=ref)
# save(list = ls(all.names = TRUE), file = "total.RData", envir = environment())
par <- getFromNamespace(".format_par", ns="phenology")(c(resfit, fixed.parameters), names(resul$data[kl]), model_before=model_before)
# C'est quoi ça ? Je retire 26/8/2019
# par <- par[1:(length(par)-9)]
sepfixed <- fixed.parameters[strtrim(names(fixed.parameters), 3)=="se#"]
if (!is.null(sepfixed) & (!identical(unname(sepfixed), numeric(0)))) names(sepfixed) <- substring(names(sepfixed), 4)
se <- c(res_se, sepfixed)
se <- getFromNamespace(".format_par", ns="phenology")(se, names(resul$data[kl]), model_before=model_before)
# C'est quoi ça ? Je retire 26/8/2019
# se <- se[1:(length(se)-9)]
d1 <- ref + par["Peak"]
if (!silent) cat(paste("Peak: ", d1, "\n", sep=""))
intdtout <- c(intdtout, Peak=unname(d1))
if (!is.na(se["Peak"])) {
d2 <- d1-1.96*se["Peak"]
d3 <- d1+1.96*se["Peak"]
if (!silent) cat(paste("confidence interval:", d2, " to ", d3, "\n", sep=""))
intdtout <- c(intdtout, PeakCI1=unname(d2), PeakCI2=unname(d3))
} else {
if (!silent) cat(paste("confidence interval not available\n", sep=""))
intdtout <- c(intdtout, PeakCI1=NA, PeakCI2=NA)
}
d1 <- ref+par["Begin"]
if (!silent) cat(paste("Begin: ", d1, "\n", sep=""))
intdtout <- c(intdtout, Begin=unname(d1))
# pour l'intervalle de confiance, il faut modifier soit directement
# Begin
# Peak Length
# Peak LengthB
d2 <- NULL
d3 <- NULL
if (!is.na(se["Begin"])) {
d2 <- ref+par["Begin"]-1.96*se["Begin"]
d3 <- ref+par["Begin"]+1.96*se["Begin"]
} else {
sel <- 0
l <- NA
if (!is.na(se["Length"])) {
l <- par["Length"]
sel <- se["Length"]
} else {
if (!is.na(se["LengthB"])) {
l <- par["LengthB"]
sel <- se["LengthB"]
}
}
if (!is.na(se["Peak"])) {
d2 <- ref+par["Peak"]-1.96*se["Peak"]-l-1.96*sel
d3 <- ref+par["Peak"]+1.96*se["Peak"]-l+1.96*sel
} else {
d2 <- ref+par["Peak"]-l-1.96*sel
d3 <- ref+par["Peak"]-l+1.96*sel
}
}
if (!is.null(d2) & !is.na(d2)) {
if (!silent) cat(paste("confidence interval:", d2, " to ", d3, "\n", sep=""))
intdtout <- c(intdtout, BeginCI1=unname(d2), BeginCI2=unname(d3))
} else {
if (!silent) cat(paste("confidence interval not available\n", sep=""))
intdtout <- c(intdtout, BeginCI1=NA, BeginCI2=NA)
}
d1 <- ref+par["End"]
if (!silent) cat(paste("End: ", d1, "\n", sep=""))
intdtout <- c(intdtout, End=unname(d1))
# pour l'intervalle de confiance, il faut modifier soit directement
# End
# Peak Length
# Peak LengthE
d2 <- NULL
d3 <- NULL
if (!is.na(se["End"])) {
d2 <- ref+par["End"]-1.96*se["End"]
d3 <- ref+par["End"]+1.96*se["End"]
} else {
sel <- 0
l <- NA
if (!is.na(se["Length"])) {
l <- par["Length"]
sel <- se["Length"]
} else {
if (!is.na(se["LengthE"])) {
l <- par["LengthE"]
sel <- se["LengthE"]
}
}
if (!is.na(se["Peak"])) {
d2 <- ref+par["Peak"]-2*se["Peak"]+l-1.96*sel
d3 <- ref+par["Peak"]+2*se["Peak"]+l+1.96*sel
} else {
d2 <- ref+par["Peak"]+l-1.96*sel
d3 <- ref+par["Peak"]+l+1.96*sel
}
}
if (!is.null(d2) & !is.na(d2)) {
if (!silent) cat(paste("confidence interval:", d2, " to ", d3, "\n", sep=""))
intdtout <- c(intdtout, EndCI1=unname(d2), EndCI2=unname(d3))
} else {
if (!silent) cat(paste("confidence interval not available\n", sep=""))
intdtout <- c(intdtout, EndCI1=NA, EndCI2=NA)
}
dtout <- c(dtout, list(intdtout))
}
names(dtout) <- names(resul$data)
resul$Dates <- dtout
resul <- addS3Class(resul, "phenology")
if (!silent)
if (is.finite(resul$value)) {
cat(paste("-Ln L=", format(resul$value, digits=max(3, trunc(log10(abs(resul$value)))+4)), "\n", sep=""))
cat(paste("Parameters=", format(length(resul$par), digits=max(3, trunc(log10(length(resul$par)))+4)), "\n", sep=""))
cat(paste("AIC=", format(2*resul$value+2*length(resul$par), digits=max(3, trunc(log10(abs(2*resul$value+2*length(resul$par))))+4)), "\n", sep=""))
} else {
cat("-Ln L=-Inf\n")
cat(paste("Parameters=", format(length(resul$par), digits=max(3, trunc(log10(length(resul$par)))+4)), "\n", sep=""))
cat("AIC=-Inf\n")
}
resul$cofactors <- cofactors
resul$add.cofactors <- add.cofactors
resul$model_before <- model_before
resul$zero <- zero
return(resul)
}
Try the phenology package in your browser
Any scripts or data that you put into this service are public.
phenology documentation built on Oct. 16, 2023, 9:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions fit_phenology
Documented in fit_phenology
#' fit_phenology fits parameters to timeseries.
#' @title Fit the phenology parameters to timeseries of counts.
#' @author Marc Girondot \email{marc.girondot@@gmail.com}
#' @return Return a list of with data and result
#' @param data A dataset generated by add_format
#' @param fixed.parameters Set of fixed parameters
#' @param fitted.parameters Set of parameters to be fitted
#' @param lower Lower bound for each parameter
#' @param upper Upper bound for each parameter
#' @param hessian If FALSE does not estimate se of parameters
#' @param model_before The change of parameters before to estimate daily counts.
#' @param cofactors data.frame with a column Date and a column for each cofactor
#' @param add.cofactors Names of the column of parameter cofactors to use as a cofactor
#' @param silent If TRUE does not show any message
#' @param zero If the theoretical nest number is under this value, this value will be used
#' @param stop.fit If TRUE, will stop search for parameters even if not ML
#' @param control List for control parameters for optim
#' @param method Method used by optim. Several can be setup.
#' @param method_Snbinom Can be Furman, exact, or saddlepoint.
#' @param store.intermediate TRUE or FALSE to save the intermediates
#' @param file.intermediate Name of the file where to save the intermediates as a list
#' @description Function of the package phenology to fit parameters to timeseries.\cr
#' To fit data, the syntax is :\cr
#' Result <- fit_phenology(data=dataset, fitted.parameters=par, fixed.parameters=pfixed, trace=1, hessian=TRUE)\cr
#' or if no parameter is fixed :\cr
#' Result <- fit_phenology(data=dataset, fitted.parameters=par)\cr
#' Add trace=1 [default] to have information on the fit progression or trace=0 to hide information on the fit progression.\cr
#' hessian = FALSE does not estimate Hessian matrix and SE of parameters.\cr
#' If the parameter Theta is fixed to +Inf, a Poissonian model of daily nest distribution is implemented.\cr
#' Special section about cofactors:\cr
#' cofactors must be a data.frame with a column Date and a column for each cofactor\cr
#' add.cofactors are the names of the column of parameter cofactors to use as a cofactor;\cr
#' The model is then: parameter[add.cofactors] * cofactor[, add.cofactors]\cr
#' If the name of the parameter is paste0(add.cofactors, "multi"), then the model is:\cr
#' parameter[paste0(add.cofactors, "multi")] * cofactor[, add.cofactors] *
#' (number of nests without cofactor)\cr
#' About parallel computing:\cr
#' Set options mc.cores and forking to tell what sort of parallel computing\cr
#' Example:\cr
#' options(mc.cores = detectCores())\cr
#' options(forking = FALSE)
#' @family Phenology model
#' @examples
#' \dontrun{
#' 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)
#'
#' }
#' @export
fit_phenology <-
function(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-9,
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")) {
# data=NULL
# lower = 0
# upper = Inf
# fitted.parameters=NULL
# fixed.parameters=NA
# cofactors=NULL
# add.cofactors=NULL
# hessian=TRUE
# silent=FALSE
# zero=1E-9
# control=list(trace=1, REPORT=1, maxit=1000)
# method = c("Nelder-Mead", "L-BFGS-B")
# stop.fit=FALSE
# store.intermediate=FALSE
# file.intermediate="Intermediate.rda"
# data_Gratiot <- add_phenology(Gratiot, name="Complete", reference=as.Date("2001-01-01"), format="%d/%m/%Y")
# data=data_Gratiot; fitted.parameters=result_Gratiot$par; fixed.parameters=result_Gratiot$fixed.parameters; trace=1
# if (is.null(fixed.parameters)) {fixed.parameters<-NA}
# if (!requireNamespace("optimx", quietly = TRUE)) {
# stop("optimx package is absent; Please install it first")
# }
# method <- c("Nelder-Mead", "L-BFGS-B")
# if (length(fitted.parameters) == 1) method <- "Brent"
# library("optimx")
Lnegbin <- getFromNamespace(".Lnegbin", ns="phenology")
if (store.intermediate) {
store <- list()
save(store, file=file.intermediate)
}
if (!inherits(data, "phenologydata")) {
message("Data should have been formated first using the function add_phenology(). I format it now.")
data <- add_phenology(data)
}
if (is.null(fitted.parameters)) {
message("No initial parameters set has been defined. I estimate one set using par_init().")
fitted.parameters <- par_init(data, fixed.parameters=fixed.parameters)
}
if ((!is.null(add.cofactors)) & (!is.null(cofactors))) {
cf1 <- cofactors
for (ff in add.cofactors)
cf1[, ff] <- cofactors[, ff] - mean(cofactors[, ff])
cofactors <- cf1
}
preLNL <- +Inf
if (control[["maxit"]] != 0) {
repeat {
for (m in seq_along(method)) {
if ((method[m] == "L-BFGS-B") | (method[m] == "Brent")) {
resul <- optim(par=fitted.parameters, fn=Lnegbin,
pt=list(data=data, fixed=fixed.parameters,
model_before=model_before,
out=TRUE, cofactors=cofactors,
add.cofactors=add.cofactors, zero=zero,
namespar=names(fitted.parameters),
method_Snbinom=method_Snbinom,
store.intermediate=store.intermediate,
file.intermediate=file.intermediate),
lower = lower,
upper = upper,
method = method[m],
control=control,
hessian=FALSE)
} else {
resul <- optim(par=fitted.parameters, fn=Lnegbin,
pt=list(data=data, fixed=fixed.parameters,
model_before=model_before,
out=TRUE, cofactors=cofactors,
add.cofactors=add.cofactors, zero=zero,
namespar=names(fitted.parameters),
method_Snbinom=method_Snbinom,
store.intermediate=store.intermediate,
file.intermediate=file.intermediate),
method = method[m],
control=control,
hessian=FALSE)
}
resfit <- resul$par
if (is.null(names(resfit))) {
names(resfit) <- names(fitted.parameters)
}
if (any(grepl("^Peak", names(resfit)))) resfit[substr(names(resfit), 1, 4)=="Peak"] <- abs(resfit[substr(names(resfit), 1, 4)=="Peak"])
if (any(grepl("^Theta", names(resfit)))) resfit[substr(names(resfit), 1, 5)=="Theta"] <- abs(resfit[substr(names(resfit), 1, 5)=="Theta"])
if (any(grepl("^PMin", names(resfit)))) resfit["PMin"] <- abs(resfit["PMin"])
# 28/3/2018. J'avais oublié ceux là
if (any(grepl("^PMinE", names(resfit)))) resfit["PMinE"] <- abs(resfit["PMinE"])
if (any(grepl("^PMinB", names(resfit)))) resfit["PMinB"] <- abs(resfit["PMinB"])
if (any(grepl("^Flat", names(resfit)))) resfit["Flat"] <- abs(resfit["Flat"])
if (any(grepl("^Length", names(resfit)))) resfit[substr(names(resfit), 1, 6)=="Length"] <- abs(resfit[substr(names(resfit), 1, 6)=="Length"])
# Ca fait en même temps MinE et MinB
if (any(grepl("^Min", names(resfit)))) resfit[substr(names(resfit), 1, 3)=="Min"]<-abs(resfit[substr(names(resfit), 1, 3)=="Min"])
if (any(grepl("^PMax", names(resfit)))) resfit[substr(names(resfit), 1, 3)=="Max"]<-abs(resfit[substr(names(resfit), 1, 3)=="Max"])
resul$par <- resfit
fitted.parameters <- resfit
# 22/8/2022
fake <- d(fitted.parameters)
}
conv <- resul$convergence
value <- resul$value
if ((conv == 0) | (abs(preLNL - value)<1e-6) | (stop.fit)) break
if (!silent) message("Convergence is not achieved. Optimization continues !")
preLNL <- value
}
if ((conv != 0) & (!stop.fit)) message("Convergence is not achieved but change in likelihood is <1e-6. Use the results with caution.")
if ((conv != 0) & (stop.fit)) message("Convergence is not achieved but maximum number of iterations is reached. Use the results with caution.")
if (!silent) {
message("Fit done!")
if (is.finite(value)) {
cat(paste("-Ln L=", format(value, digits=max(3, trunc(log10(abs(value)))+4)), "\n", sep=""))
} else {
cat("-Ln L= -Inf\n")
}
}
result_list <- resul
# result_list <- list()
# result_list$par <- resfit
# result_list$value <- value
# result_list$convergence <- conv
# resul <- result_list
} else {
# J'avais maxit == 0
resul <- list()
resul$data <- data
resul$par <- fitted.parameters
resul$value <- NA
resul$convergence <- NA
resfit <- fitted.parameters
}
if (hessian) {
if (!requireNamespace("numDeriv", quietly = TRUE)) {
stop("numDeriv package is absent; Please install it first")
}
if (!silent) message("Estimation of the standard error of parameters. Be patient please.")
mathessian <- try(getFromNamespace("hessian", ns="numDeriv")(func=Lnegbin,
x=resfit,
method="Richardson",
pt=list(data=data, fixed=fixed.parameters,
model_before=model_before,
out=TRUE, cofactors=cofactors,
add.cofactors=add.cofactors, zero=zero,
namespar=names(fitted.parameters),
method_Snbinom=method_Snbinom,
store.intermediate=store.intermediate,
file.intermediate=file.intermediate))
, silent=TRUE)
if (substr(mathessian[1], 1, 5)=="Error") {
if (!silent) {
warning("Error in the fit; probably one or more parameters are not estimable.")
warning("Standard errors of parameters cannot be estimated.")
}
res_se <- rep(NA, length(resfit))
names(res_se) <- names(resfit)
} else {
rownames(mathessian) <- colnames(mathessian) <- names(resfit)
resul$hessian <- mathessian
res_se <- SEfromHessian(mathessian, silent=silent)
}
} else {
if (!silent) warning("Standard errors are not estimated.")
res_se <- rep(NA, length(resfit))
names(res_se) <- names(resfit)
}
resul$se <- res_se
resul$fixed.parameters <- fixed.parameters
resul$method_incertitude <- "convolution"
resul$method_Snbinom <- method_Snbinom
resul$data <- data
# Lnegbin(x=resfit, pt=list(data=data, fixed=fixed.parameters,
# out=FALSE, cofactors=cofactors,
# method_Snbinom=method_Snbinom,
# add.cofactors=add.cofactors, zero=zero))
for(kl in 1:length(res_se)) {
if (is.na(res_se[kl])) {
if (!silent) cat(paste(names(resfit[kl]), "=", format(resfit[kl], digits=max(3, trunc(log10(abs(resfit[kl])))+4)), " SE= NaN\n", sep=""))
} else {
if (!silent) cat(paste(names(resfit[kl]), "=", format(resfit[kl], digits=max(3, trunc(log10(abs(resfit[kl])))+4)), " SE=", format(res_se[kl], digits=max(3, trunc(log10(abs(res_se[kl])))+4)), "\n", sep=""))
}
}
dtout <- list()
for (kl in 1:length(resul$data)) {
if (!silent) cat(paste("Series: ", names(resul$data[kl]), "\n", sep=""))
# la date de référence est resul$data[[kl]][1, "Date"]-resul$data[[kl]][1, "ordinal"]
ref <- resul$data[[kl]][1, "Date"]-resul$data[[kl]][1, "ordinal"]
intdtout <- c(reference=ref)
# save(list = ls(all.names = TRUE), file = "total.RData", envir = environment())
par <- getFromNamespace(".format_par", ns="phenology")(c(resfit, fixed.parameters), names(resul$data[kl]), model_before=model_before)
# C'est quoi ça ? Je retire 26/8/2019
# par <- par[1:(length(par)-9)]
sepfixed <- fixed.parameters[strtrim(names(fixed.parameters), 3)=="se#"]
if (!is.null(sepfixed) & (!identical(unname(sepfixed), numeric(0)))) names(sepfixed) <- substring(names(sepfixed), 4)
se <- c(res_se, sepfixed)
se <- getFromNamespace(".format_par", ns="phenology")(se, names(resul$data[kl]), model_before=model_before)
# C'est quoi ça ? Je retire 26/8/2019
# se <- se[1:(length(se)-9)]
d1 <- ref + par["Peak"]
if (!silent) cat(paste("Peak: ", d1, "\n", sep=""))
intdtout <- c(intdtout, Peak=unname(d1))
if (!is.na(se["Peak"])) {
d2 <- d1-1.96*se["Peak"]
d3 <- d1+1.96*se["Peak"]
if (!silent) cat(paste("confidence interval:", d2, " to ", d3, "\n", sep=""))
intdtout <- c(intdtout, PeakCI1=unname(d2), PeakCI2=unname(d3))
} else {
if (!silent) cat(paste("confidence interval not available\n", sep=""))
intdtout <- c(intdtout, PeakCI1=NA, PeakCI2=NA)
}
d1 <- ref+par["Begin"]
if (!silent) cat(paste("Begin: ", d1, "\n", sep=""))
intdtout <- c(intdtout, Begin=unname(d1))
# pour l'intervalle de confiance, il faut modifier soit directement
# Begin
# Peak Length
# Peak LengthB
d2 <- NULL
d3 <- NULL
if (!is.na(se["Begin"])) {
d2 <- ref+par["Begin"]-1.96*se["Begin"]
d3 <- ref+par["Begin"]+1.96*se["Begin"]
} else {
sel <- 0
l <- NA
if (!is.na(se["Length"])) {
l <- par["Length"]
sel <- se["Length"]
} else {
if (!is.na(se["LengthB"])) {
l <- par["LengthB"]
sel <- se["LengthB"]
}
}
if (!is.na(se["Peak"])) {
d2 <- ref+par["Peak"]-1.96*se["Peak"]-l-1.96*sel
d3 <- ref+par["Peak"]+1.96*se["Peak"]-l+1.96*sel
} else {
d2 <- ref+par["Peak"]-l-1.96*sel
d3 <- ref+par["Peak"]-l+1.96*sel
}
}
if (!is.null(d2) & !is.na(d2)) {
if (!silent) cat(paste("confidence interval:", d2, " to ", d3, "\n", sep=""))
intdtout <- c(intdtout, BeginCI1=unname(d2), BeginCI2=unname(d3))
} else {
if (!silent) cat(paste("confidence interval not available\n", sep=""))
intdtout <- c(intdtout, BeginCI1=NA, BeginCI2=NA)
}
d1 <- ref+par["End"]
if (!silent) cat(paste("End: ", d1, "\n", sep=""))
intdtout <- c(intdtout, End=unname(d1))
# pour l'intervalle de confiance, il faut modifier soit directement
# End
# Peak Length
# Peak LengthE
d2 <- NULL
d3 <- NULL
if (!is.na(se["End"])) {
d2 <- ref+par["End"]-1.96*se["End"]
d3 <- ref+par["End"]+1.96*se["End"]
} else {
sel <- 0
l <- NA
if (!is.na(se["Length"])) {
l <- par["Length"]
sel <- se["Length"]
} else {
if (!is.na(se["LengthE"])) {
l <- par["LengthE"]
sel <- se["LengthE"]
}
}
if (!is.na(se["Peak"])) {
d2 <- ref+par["Peak"]-2*se["Peak"]+l-1.96*sel
d3 <- ref+par["Peak"]+2*se["Peak"]+l+1.96*sel
} else {
d2 <- ref+par["Peak"]+l-1.96*sel
d3 <- ref+par["Peak"]+l+1.96*sel
}
}
if (!is.null(d2) & !is.na(d2)) {
if (!silent) cat(paste("confidence interval:", d2, " to ", d3, "\n", sep=""))
intdtout <- c(intdtout, EndCI1=unname(d2), EndCI2=unname(d3))
} else {
if (!silent) cat(paste("confidence interval not available\n", sep=""))
intdtout <- c(intdtout, EndCI1=NA, EndCI2=NA)
}
dtout <- c(dtout, list(intdtout))
}
names(dtout) <- names(resul$data)
resul$Dates <- dtout
resul <- addS3Class(resul, "phenology")
if (!silent)
if (is.finite(resul$value)) {
cat(paste("-Ln L=", format(resul$value, digits=max(3, trunc(log10(abs(resul$value)))+4)), "\n", sep=""))
cat(paste("Parameters=", format(length(resul$par), digits=max(3, trunc(log10(length(resul$par)))+4)), "\n", sep=""))
cat(paste("AIC=", format(2*resul$value+2*length(resul$par), digits=max(3, trunc(log10(abs(2*resul$value+2*length(resul$par))))+4)), "\n", sep=""))
} else {
cat("-Ln L=-Inf\n")
cat(paste("Parameters=", format(length(resul$par), digits=max(3, trunc(log10(length(resul$par)))+4)), "\n", sep=""))
cat("AIC=-Inf\n")
}
resul$cofactors <- cofactors
resul$add.cofactors <- add.cofactors
resul$model_before <- model_before
resul$zero <- zero
return(resul)
}
Try the phenology package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.