R/summary.phenology.R

In phenology: Tools to Manage a Parametric Function that Describes Phenology and More

#' summary.phenology prints the information from a result object.
#' @title Print the result information from a result object.
#' @author Marc Girondot \email{marc.girondot@@gmail.com}
#' @return A list with five objects: synthesis is a data.frame with global estimate of nesting.\cr
#' details_MCMC, details_ML, details_mean are lists with day by day information for each series, and 
#' sum_synthesis is the synthesis of the sum of all analyzed time-series.
#' @param object A result file generated by fit_phenology
#' @param resultmcmc A mcmc object
#' @param season The number of season to analyze
#' @param series Names of the series to be analyzed or "all"
#' @param chain The number of chain to be used in resultmcmc
#' @param replicate.CI.mcmc Number of iterations to be used or "all"
#' @param replicate.CI Number of replicates for ML resampling
#' @param level Level to estimate confidence interval or credibility interval
#' @param print Should information be shown
#' @param ... Not used
#' @description The function summary.phenology shows result and estimates confidence interval.\cr
#' If several years are analyzed, the sum_synthesis result can be obtained only if there is 
#' not a mix of bisextile and non-bisextile years.
#' @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)
#' # Display information from the result
#' s <- summary(result_Gratiot)
#' # Using mcmc
#' s <- summary(object=result_Gratiot, resultmcmc=result_Gratiot_mcmc)
#' }
#' @method summary phenology
#' @export



summary.phenology <- function(object                    , 
                              resultmcmc = NULL         , 
                              season = NULL             , 
                              chain = 1                 , 
                              series = "all"            , 
                              replicate.CI.mcmc = "all" , 
                              replicate.CI = 10000      ,
                              level= 0.95               ,
                              print = TRUE              ,
                              ...) {
  
  
  # resultmcmc = NULL;season = NULL; chain = 1;series = "all";replicate.CI.mcmc = "all";replicate.CI = 10000;level= 0.95;print = TRUE
  # object=result_Gratiot; resultmcmc=NULL; season = NULL; level=0.95;chain=1; series="all"; replicate.CI.mcmc = "all"; replicate.CI = 10000; print=TRUE
  # # object=result_Gratiot; resultmcmc=result_Gratiot_mcmc
  
  formatpar <- getFromNamespace(".format_par", ns="phenology")
  dailycount <- getFromNamespace(".daily_count", ns="phenology")
  
  model_before <- object$model_before
  
  
  nday <- 366 * (max(unlist(lapply(object$data, FUN = function(x) max(c(x[, "ordinal"], x[, "ordinal2"]), na.rm = TRUE))) %/% 366) + 1)
  
  if (print) {
    cat(paste("Number of timeseries: ", length(object$data), "\n", sep=""))
    for (i in 1:length(object$data)) {
      cat(paste(names(object$data[i]), "\n", sep=""))
    }
    cat(paste("Date uncertainty management: ", object$method_incertitude, "\n", sep=""))
    
    cat("Fitted parameters:\n")
    for (i in 1:length(object$par)) {
      cat(paste(names(object$par[i]), "=", object$par[i], " SE ", object$se[i], "\n", sep=""))
    }
    if (length(object$fixed.parameters)>0) {
      cat("Fixed parameters:\n")
      for (i in 1:length(object$fixed.parameters)) {
        cat(paste(names(object$fixed.parameters[i]), "=", object$fixed.parameters[i], "\n", sep=""))
      }
    }
    cat(paste("Ln L: ", object$value, "\n", sep=""))
    cat(paste("Parameter number: ", length(object$par), "\n", sep=""))
    cat(paste("AIC: ", 2*object$value+2*length(object$par), "\n", sep=""))
  }
  
  # chain <- 1
  # Name of the series for which we want an estimate
  if (is.numeric(series)) series <- names(object$data)[series]
  if (any(series == "all")) series <- names(object$data)
  nseries <- length(series)
  rna <- rep(NA, nseries)
  
  probs <- c((1-level)/2, 0.5, 1-(1-level)/2)
  
  
  
  retdf <- data.frame(series=series, 
                      "without_obs_Mean"=rna,
                      "with_obs_Mean"=rna,
                      
                      "without_obs_Low_ML"=rna, 
                      "without_obs_Median_ML"=rna, 
                      "without_obs_High_ML"=rna, 
                      "without_obs_Mean_ML"=rna, 
                      "without_obs_Var_ML"=rna, 
                      
                      "with_obs_Low_ML"=rna, 
                      "with_obs_Median_ML"=rna, 
                      "with_obs_High_ML"=rna, 
                      "with_obs_Mean_ML"=rna, 
                      "with_obs_Var_ML"=rna, 
                      
                      "without_obs_Low_MCMC"=rna, 
                      "without_obs_Median_MCMC"=rna, 
                      "without_obs_High_MCMC"=rna, 
                      "without_obs_Mean_MCMC"=rna, 
                      "without_obs_Var_MCMC"=rna,
                      
                      "with_obs_Low_MCMC"=rna, 
                      "with_obs_Median_MCMC"=rna, 
                      "with_obs_High_MCMC"=rna, 
                      "with_obs_Mean_MCMC"=rna, 
                      "with_obs_Var_MCMC"=rna, 
                      
                      "NbObservations"=rna, 
                      "NbMonitoredDays"=rna,
                      stringsAsFactors = FALSE)
  rownames(retdf) <- series
  klist_mcmc <- list()
  klist_ML <- list()
  klist_Mean <- list()
  
  Sum_dailydata_mean <- matrix(0, nrow = nday, ncol=1)
  Sum_dailydata_withObs_mean <- Sum_dailydata_mean
  
  if (!is.null(object$hessian)) {
    Sum_dailydata_ML <- matrix(0, nrow = nday, ncol=replicate.CI)
    Sum_dailydata_withObs_ML <- Sum_dailydata_ML
  } else {
    Sum_dailydata_ML <- NULL
    Sum_dailydata_withObs_ML <- NULL
  }
  
  replicate.CI.mcmc.x <- NULL
  if (!is.null(resultmcmc)) {
    if (replicate.CI.mcmc == "all") {
      replicate.CI.mcmc.x <- nrow(resultmcmc$resultMCMC[[chain]])
    } else {
      replicate.CI.mcmc.x <- replicate.CI.mcmc
    }
    Sum_dailydata_MCMC <- matrix(0, nrow = nday, ncol=replicate.CI.mcmc.x)
    Sum_dailydata_withObs_MCMC <- Sum_dailydata_MCMC
  } else {
    Sum_dailydata_MCMC <- NULL
    Sum_dailydata_withObs_MCMC <- NULL
  }
  
  bisextile <- NULL
  
  # nmser <- series[1]
  
  sp <- NULL
  
  for (nmser in series) {
    
    synthesisPontes_MCMC <- NA
    synthesisPontes_ML <- NA
    
    if (print) {
      tx <- paste0("Timeseries: ", nmser)
      cat(paste0(rep("-", nchar(tx)), collapse=""), "\n")
      cat(tx, "\n")
      cat(paste0(rep("-", nchar(tx)), collapse=""), "\n")
    }
    
    dref <- object$Dates[[nmser]]["reference"]
    # dref <- attributes(object$data[[nmser]])[["reference"]]
    # nday <- ifelse(as.POSIXlt(dref+365)$mday==as.POSIXlt(dref)$mday, 365, 366)
    
    dta_ec <- object$data[[nmser]]
    dta_ec <- dta_ec[is.na(dta_ec$A), ]
    # Observed counts
    observedPontes <- data.frame(ordinal=dta_ec[dta_ec$CountTypes == "exact", "ordinal"], 
                                 observed=dta_ec[dta_ec$CountTypes == "exact", "nombre"])
    # je mets des 0 à toutes les dates supplémentaires de la série
    if (any(!is.na(dta_ec[dta_ec$CountTypes == "exact", "ordinal2"]))) {
      for (i in which((!is.na(dta_ec[, "ordinal2"]) & (dta_ec$CountTypes == "exact")))) {
        rnge <- (dta_ec[i, "ordinal"]+1):(dta_ec[i, "ordinal2"])
        observedPontes <- rbind(observedPontes, 
                                data.frame(ordinal= rnge, 
                                           observed=rep(0, length(rnge))))
      }
    }
    
    
    
    
    parg <- formatpar(c(object$par, object$fixed.parameters), nmser, model_before=model_before, season=season)
    
    cof <- NULL
    if ((!is.null(object$add.cofactors)) & (!is.null(object$cofactors))) {
      
      j <- 0:(nday-1)  
      # à analyser
      cof <- object$cofactors[object$cofactors$Date %in% (dref+j), ]
      cof <- cof[, -1, drop=FALSE]
      cof <- as.data.frame(cbind(Date=j, cof))
    }
    
    
    dc_mean <- dailycount(d=0:(nday-1), xpar=parg, 
                          cofactors=cof, 
                          add.cofactors=object$add.cofactors, 
                          print=FALSE, zero=1E-9)
    retdf[nmser, "without_obs_Mean"] <- sum(dc_mean)
    
    retdf[nmser, "NbObservations"] <- sum(observedPontes$observed)
    retdf[nmser, "NbMonitoredDays"] <- nrow(observedPontes)
    
    bisextile <- c(bisextile, ifelse(length(dc_mean)==365, 365, 366))
    Sum_dailydata_mean <- Sum_dailydata_mean + ifelse(length(dc_mean)==365, c(dc_mean, 0), dc_mean)
    
    if (print) {
      cat("Total estimate not taking into account the observations: ")
      cat(paste0("Mean=", retdf[nmser, "without_obs_Mean"], "\n"))
    }
    
    # 3/11/2018 Pourquoi Theta seulement dans fixed.parameters ?
    # Non, c'est bon, il cherche sur les 2
    
    SDMin <- NULL
    SDMax <- NULL
    for (mu in dc_mean) {
      qnb <- qnbinom(p = c(probs[1], probs[3]), 
                     size=parg["Theta"], 
                     mu=mu)
      SDMin <- c(SDMin, qnb[1])
      SDMax <- c(SDMax, qnb[2])
    }
    
    dc_mean <- data.frame(Date=dref+(0:(nday-1)), Ordinal = 0:(nday-1), Mean=NA, SD.Low=SDMin, SD.High=SDMax, Observed=NA, Modelled=dc_mean)
    dc_mean[match(observedPontes[, "ordinal"], dc_mean[, "Ordinal"]), "Observed"] <- observedPontes[, "observed"]
    dc_mean[, "Mean"] <- ifelse(is.na(dc_mean[, "Observed"]), dc_mean[, "Modelled"], dc_mean[, "Observed"])
    
    Sum_dailydata_withObs_mean <- Sum_dailydata_withObs_mean + ifelse(length(dc_mean[, "Mean"])==365, c(dc_mean[, "Mean"], 0), dc_mean[, "Mean"])
    
    if (!is.null(cof)) {
      dc_mean <- cbind(dc_mean, cof[, -1, drop=FALSE])
    }
    
    rownames(dc_mean) <- dc_mean[, "Ordinal"]
    k <- list(dc_mean)
    names(k) <- nmser
    
    klist_Mean <- c(klist_Mean, k)
    
    # 22/9/2019: Je mets + 1 sinon peut faire 0 et provoque une erreur
    retdf[nmser, "with_obs_Mean"] <- sum(dc_mean[, "Mean"])
    
    if (print) {
      cat("Total estimate taking into account the observations: ")
      cat(paste0("Mean=", retdf[nmser, "with_obs_Mean"], "\n"))
    }
    
    pfixed <- object$fixed.parameters
    sepfixed <- pfixed[strtrim(names(pfixed), 3)=="se#"]
    if (identical(unname(sepfixed), numeric(0))) sepfixed <- NULL
    pfixed <- pfixed[strtrim(names(pfixed), 3) != "se#"]
    if (!is.null(sepfixed)) names(sepfixed) <- substring(names(sepfixed), 4)
    
    # J'ai un sd sur des paramètres fixés
    # if (!is.null(sepfixed) & (!identical(unname(sepfixed), numeric(0)))) {
    
    pfixed.df <- data.frame()
    pfixed.df.mcmc <- data.frame()
    
    
    
    if (!is.null(pfixed)) {
      for (i in seq_along(pfixed)) {
        dfadd <- data.frame()
        dfadd.mcmc <- data.frame()
        
        if (!is.null(sepfixed)) {
          if (!is.null(replicate.CI)) {
            dfadd <- data.frame(rnorm(n=replicate.CI, mean=unname(pfixed[i]), sd=sepfixed[names(pfixed[i])]))
            colnames(dfadd) <- names(pfixed[i])
          }
          if (!is.null(replicate.CI.mcmc.x)) {
            dfadd.mcmc <- data.frame(rnorm(n=replicate.CI.mcmc.x, mean=unname(pfixed[i]), sd=sepfixed[names(pfixed[i])]))
            colnames(dfadd.mcmc) <- names(pfixed[i])
          }
        } else {
          if (!is.null(replicate.CI)) {
            dfadd <- data.frame(rep(unname(pfixed[i]), replicate.CI))
            colnames(dfadd) <- names(pfixed[i])
          }
          if (!is.null(replicate.CI.mcmc.x)) {
            dfadd.mcmc <- data.frame(rep(unname(pfixed[i]), replicate.CI.mcmc.x))
            colnames(dfadd.mcmc) <- names(pfixed[i])
          }
        }
        if (ncol(pfixed.df.mcmc) ==0 ) {
          pfixed.df.mcmc <- dfadd.mcmc
        } else {
          pfixed.df.mcmc <- cbind(pfixed.df.mcmc, dfadd.mcmc)
        }
        if (ncol(pfixed.df) ==0 ) {
          pfixed.df <- dfadd
        } else {
          pfixed.df <- cbind(pfixed.df, dfadd)
        }
        
      }
    }
    
    pfixed.df <- as.matrix(pfixed.df)
    pfixed.df.mcmc <- as.matrix(pfixed.df.mcmc)
    
    
    lnday <- 0:(nday-1)
    # 22/9/2019: Décalage d'un jour
    opord <- observedPontes[, "ordinal"]+1
    opnumb <- observedPontes[, "observed"]
    
    if (!is.null(resultmcmc)) {
      lmcmc <- nrow(resultmcmc$resultMCMC[[chain]])
      mcmctobeused <- 1:lmcmc
      if (replicate.CI.mcmc != "all") {
        repl <- ifelse(nrow(resultmcmc$resultMCMC[[chain]]) <= replicate.CI.mcmc, TRUE, FALSE)
        mcmctobeused <- sample(x=mcmctobeused, 
                               size=replicate.CI.mcmc, 
                               replace = repl)
      } else {
        replicate.CI.mcmc <- nrow(resultmcmc$resultMCMC[[chain]])
      }
      
      
      if (ncol(pfixed.df.mcmc) != 0) {
        dailydata <- sapply(X = 1:replicate.CI.mcmc, FUN=function(xxx) {
          px <- c(unlist(resultmcmc$resultMCMC[[chain]][mcmctobeused[xxx], ]), pfixed.df.mcmc[xxx, ])
          xparec <- formatpar(px, nmser, model_before=model_before, season=season)
          dailycount(lnday, xparec, print=FALSE, cofactors=cof, 
                     add.cofactors=object$add.cofactors)
        })
      } else {
        dailydata <- sapply(X = 1:replicate.CI.mcmc, FUN=function(xxx) {
          px <- c(unlist(resultmcmc$resultMCMC[[chain]][mcmctobeused[xxx], ]))
          xparec <- formatpar(px, nmser, model_before=model_before, season=season)
          dailycount(lnday, xparec, print=FALSE, cofactors=cof, 
                     add.cofactors=object$add.cofactors)
        })
      }
      
      if (nrow(dailydata) == 365) {
        dailydata_s <- rbind(dailydata, rep(0, ncol(dailydata)))
      } else {
        dailydata_s <- dailydata
      }
      
      Sum_dailydata_MCMC <- Sum_dailydata_MCMC + dailydata_s
      synthesisPontes <- apply(X = dailydata, MARGIN = 2, FUN=sum)
      
      dailydata_withObs <- dailydata
      dailydata_withObs[opord, ] <- opnumb
      
      if (nrow(dailydata_withObs) == 365) {
        dailydata_s <- rbind(dailydata_withObs, rep(0, ncol(dailydata_withObs)))
      } else {
        dailydata_s <- dailydata_withObs
      }
      
      Sum_dailydata_withObs_MCMC <- Sum_dailydata_withObs_MCMC + dailydata_s
      synthesisPontes_withObs <- apply(X = dailydata_withObs, MARGIN = 2, FUN=sum)
      
      # synthesisPontes_withObs <- apply(X = dailydata, MARGIN = 2, FUN=function(xxx) {
      #   xxx[opord] <- opnumb
      #   sum(xxx)
      # })
      
      k <-as.data.frame(t(apply(X = dailydata, MARGIN=1, 
                                FUN = function(x) {quantile(x, probs=probs)})))
      k <- list(cbind(Ordinal=lnday, k))
      
      names(k) <- nmser
      
      klist_mcmc <- c(klist_mcmc, k)
      
      k <- unname(quantile(synthesisPontes, probs=probs))
      retdf[nmser, c("without_obs_Low_MCMC", "without_obs_Median_MCMC", "without_obs_High_MCMC")] <- k
      retdf[nmser, c("without_obs_Mean_MCMC", "without_obs_Var_MCMC")] <- c(mean(synthesisPontes), var(synthesisPontes))
      
      if (print) {
        cat("Total estimate not taking into account the observations MCMC-based:\n")
        cat(paste0("Low=", k[1], "   Median=", k[2], "   High=", k[3], "\n"))
      }
      
      k <- unname(quantile(synthesisPontes_withObs, probs=probs, na.rm = TRUE))
      retdf[nmser, c("with_obs_Low_MCMC", "with_obs_Median_MCMC", "with_obs_High_MCMC")] <- k
      retdf[nmser, c("with_obs_Mean_MCMC", "with_obs_Var_MCMC")] <- c(mean(synthesisPontes_withObs, na.rm = TRUE), var(synthesisPontes_withObs, na.rm = TRUE))
      
      if (print) {
        cat("Total estimate taking into account the observations MCMC-based:\n")
        cat(paste0("Low=", k[1], "   Median=", k[2], "   High=", k[3], "\n"))
      }
      synthesisPontes_MCMC <- synthesisPontes
    } else {
      k <- list(NA)
      names(k) <- nmser
      klist_mcmc <- c(klist_mcmc, k)
    }
    # fin du mcmc analysis
    
    
    
    # Maintenant en prenant en compte la matrice hessienne
    
    if (!is.null(object$hessian)) {
      if (all(names(object$par) %in% colnames(object$hessian))) {
        
        par2 <- RandomFromHessianOrMCMC(
          Hessian = object$hessian,
          fitted.parameters = object$par,
          fixed.parameters = object$fixed.parameters,
          method="Hessian", 
          probs = c(0.025, 0.5, 0.975),
          replicates = replicate.CI, silent=TRUE)
        par2 <- par2$random
        
        dailydata <- sapply(1:replicate.CI, FUN=function(xxx) {
          dailycount(lnday, formatpar(unlist(par2[xxx, , drop=TRUE]), nmser, model_before=model_before, season=season), 
                     print=FALSE, cofactors=cof, 
                     add.cofactors=object$add.cofactors)
        })
        
        k <- as.data.frame(t(apply(X = dailydata, MARGIN=1, 
                                   FUN = function(x) {quantile(x, probs=probs)})))
        k <- list(cbind(Ordinal=lnday, k))
        
        names(k) <- nmser
        
        klist_ML <- c(klist_ML, k)
        
        if (nrow(dailydata) == 365) {
          dailydata_s <- rbind(dailydata, rep(0, ncol(dailydata)))
        } else {
          dailydata_s <- dailydata
        }
        
        Sum_dailydata_ML <- Sum_dailydata_ML + dailydata_s
        synthesisPontes <- apply(X = dailydata, MARGIN = 2, FUN=sum)
        
        dailydata_withObs <- dailydata
        dailydata_withObs[opord, ] <- opnumb
        
        if (nrow(dailydata_withObs) == 365) {
          dailydata_s <- rbind(dailydata_withObs, rep(0, ncol(dailydata_withObs)))
        } else {
          dailydata_s <- dailydata_withObs
        }
        
        Sum_dailydata_withObs_ML <- Sum_dailydata_withObs_ML + dailydata_s
        synthesisPontes_withObs <- apply(X = dailydata_withObs, MARGIN = 2, FUN=sum)
        
        # synthesisPontes_withObs <- apply(X = dailydata, MARGIN = 2, FUN=function(xxx) {
        #   xxx[opord] <- opnumb
        #   sum(xxx)
        # })
        
        k <- unname(quantile(synthesisPontes, probs=probs))
        retdf[nmser, c("without_obs_Low_ML", "without_obs_Median_ML", "without_obs_High_ML")] <- k
        retdf[nmser, c("without_obs_Mean_ML", "without_obs_Var_ML")] <- c(mean(synthesisPontes), var(synthesisPontes))
        
        if (print) {
          cat("Total estimate not taking into account the observations ML-based:\n")
          cat(paste0("Low=", k[1], "   Median=", k[2], "   High=", k[3], "\n"))
        }
        
        k <- unname(quantile(synthesisPontes_withObs, probs=probs, na.rm = TRUE))
        retdf[nmser, c("with_obs_Low_ML", "with_obs_Median_ML", "with_obs_High_ML")] <- k
        retdf[nmser, c("with_obs_Mean_ML", "with_obs_Var_ML")] <- c(mean(synthesisPontes_withObs, na.rm = TRUE), var(synthesisPontes_withObs, na.rm = TRUE))
        
        if (print) {
          cat("Total estimate taking into account the observations ML-based:\n")
          cat(paste0("Low=", k[1], "   Median=", k[2], "   High=", k[3], "\n"))
        }
        synthesisPontes_ML <- synthesisPontes
      } else {
        k <- list(NA)
        names(k) <- nmser
        klist_ML <- c(klist_ML, k)
      }
    } else {
      k <- list(NA)
      names(k) <- nmser
      klist_ML <- c(klist_ML, k)
    }
    

    
    
    # fin de la boucle des séries
    sp <- c(sp, k=list(c(list(ML=synthesisPontes_ML), list(MCMC=synthesisPontes_MCMC))))

  }
  
  names(sp) <- series
  
  nr <- bisextile[1]
  if (any(bisextile != nr)) {
    if (print) {
      warning("The sum of series is not estimated because bisextile years are mixed with non-bisextile years.")
    }
    Sum_final <- NULL
    
  } else {
    Sum_dailydata_mean <- Sum_dailydata_mean[1:nr]
    Sum_dailydata_withObs_mean <- Sum_dailydata_withObs_mean[1:nr]
    
    Sum_final <- data.frame(OrdinalDate=0:(nr-1), WithoutObs_Mean=Sum_dailydata_mean, WithObs_Mean=Sum_dailydata_withObs_mean)
    
    if (!is.null(Sum_dailydata_ML)) {
      Sum_dailydata_ML <- Sum_dailydata_ML[1:nr, ]
      Sum_dailydata_withObs_ML <- Sum_dailydata_withObs_ML[1:nr, ]
      essai <- t(apply(Sum_dailydata_ML, MARGIN = 1, FUN = function(xxx) quantile(xxx, probs=probs)))
      colnames(essai) <- paste0("WithoutObs_ML_", colnames(essai))
      Sum_final <- cbind(Sum_final, essai)
      essai <- t(apply(Sum_dailydata_withObs_ML, MARGIN = 1, FUN = function(xxx) quantile(xxx, probs=probs)))
      colnames(essai) <- paste0("WithObs_ML_", colnames(essai))
      Sum_final <- cbind(Sum_final, essai)
    }
    if (!is.null(Sum_dailydata_MCMC)) {
      Sum_dailydata_MCMC <- Sum_dailydata_MCMC[1:nr, ]
      Sum_dailydata_withObs_MCMC <- Sum_dailydata_withObs_MCMC[1:nr, ]
      essai <- t(apply(Sum_dailydata_MCMC, MARGIN = 1, FUN = function(xxx) quantile(xxx, probs=probs)))
      colnames(essai) <- paste0("WithoutObs_MCMC_", colnames(essai))
      Sum_final <- cbind(Sum_final, essai)
      essai <- t(apply(Sum_dailydata_withObs_MCMC, MARGIN = 1, FUN = function(xxx) quantile(xxx, probs=probs)))
      colnames(essai) <- paste0("WithObs_MCMC_", colnames(essai))
      Sum_final <- cbind(Sum_final, essai)
    }
  }
  
  
  
  
  rout <- list(synthesis=retdf, details_mcmc=klist_mcmc, 
               details_ML=klist_ML, details_Mean=klist_Mean, 
               sum_synthesis=Sum_final, sum_withoutobs_replicates=sp)
  rout <- addS3Class(rout, "phenologyout")
  return(invisible(rout))
}