R/plot.fitRMU.R

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

#' plot.fitRMU plots the results of a fit RMU.
#' @title Plot the synthesis of RMU fit.
#' @author Marc Girondot
#' @return Return A list with result of CI.RMU()
#' @param x A result file generated by fitRMU
#' @param ... Parameters used by plot
#' @param col The function used to generate colors.
#' @param border The border of polygons used to represent the proportions.
#' @param resultMCMC MCMC result for fitRUM
#' @param chain Chain to be plotted for MCMC
#' @param replicate.CI Number of replicates to estimate CI
#' @param CI.RMU A result of CI.RMU()
#' @param what Can be proportions, numbers or total
#' @param criteria What criteria will be used for proportions or numbers: mean or 50\%
#' @param aggregate Can be model or both
#' @param order Give the order of series in plot, from bottom to top. Can be used to not show series.
#' @param control.legend Parameters send to legend
#' @param names.legend Names to show in legend.
#' @param show.legend If FALSE, does not show legend
#' @description The function plot.fitRMU plots the results of fitRMU().\cr
#' In most of the cases, replicate.CI can be set to 0 for what="proportions" or "numbers".\cr
#' The parameter CI.RMU can be used when this function is used several times with the same data.
#' @family Fill gaps in RMU
#' @examples
#' \dontrun{
#' library("phenology")
#' RMU.names.AtlanticW <- data.frame(mean=c("Yalimapo.French.Guiana", 
#'                                          "Galibi.Suriname", 
#'                                          "Irakumpapy.French.Guiana"), 
#'                                  se=c("se_Yalimapo.French.Guiana", 
#'                                       "se_Galibi.Suriname", 
#'                                       "se_Irakumpapy.French.Guiana"), stringsAsFactors = FALSE)
#' data.AtlanticW <- data.frame(Year=c(1990:2000), 
#'       Yalimapo.French.Guiana=c(2076, 2765, 2890, 2678, NA, 
#'                                6542, 5678, 1243, NA, 1566, 1566),
#'       se_Yalimapo.French.Guiana=c(123.2, 27.7, 62.5, 126, NA, 
#'                                  230, 129, 167, NA, 145, 20),
#'       Galibi.Suriname=c(276, 275, 290, NA, 267, 
#'                        542, 678, NA, 243, 156, 123),
#'       se_Galibi.Suriname=c(22.3, 34.2, 23.2, NA, 23.2, 
#'                            4.3, 2.3, NA, 10.3, 10.1, 8.9),
#'       Irakumpapy.French.Guiana=c(1076, 1765, 1390, 1678, NA, 
#'                                3542, 2678, 243, NA, 566, 566),
#'       se_Irakumpapy.French.Guiana=c(23.2, 29.7, 22.5, 226, NA, 
#'                                  130, 29, 67, NA, 15, 20), stringsAsFactors = FALSE)
#'                            
#' cst <- fitRMU(RMU.data=data.AtlanticW, RMU.names=RMU.names.AtlanticW, 
#'                colname.year="Year", model.trend="Constant", 
#'                model.SD="Zero")
#' expo <- fitRMU(RMU.data=data.AtlanticW, RMU.names=RMU.names.AtlanticW, 
#'                colname.year="Year", model.trend="Exponential", 
#'                model.SD="Zero")
#' YS <- fitRMU(RMU.data=data.AtlanticW, RMU.names=RMU.names.AtlanticW, 
#'              colname.year="Year", model.trend="Year-specific", 
#'              model.SD="Zero")
#' YS1 <- fitRMU(RMU.data=data.AtlanticW, RMU.names=RMU.names.AtlanticW, 
#'              colname.year="Year", model.trend="Year-specific", 
#'              model.SD="Zero", model.rookeries="First-order")
#' YS1_cst <- fitRMU(RMU.data=data.AtlanticW, RMU.names=RMU.names.AtlanticW, 
#'              colname.year="Year", model.trend="Year-specific", 
#'              model.SD="Constant", model.rookeries="First-order", 
#'              parameters=YS1$par)
#' YS2 <- fitRMU(RMU.data=data.AtlanticW, RMU.names=RMU.names.AtlanticW, 
#'              colname.year="Year", model.trend="Year-specific",
#'              model.SD="Zero", model.rookeries="Second-order", 
#'              parameters=YS1$par)
#' YS2_cst <- fitRMU(RMU.data=data.AtlanticW, RMU.names=RMU.names.AtlanticW, 
#'              colname.year="Year", model.trend="Year-specific",
#'              model.SD="Constant", model.rookeries="Second-order", 
#'              parameters=YS1_cst$par)
#'                
#' compare_AIC(Constant=cst, Exponential=expo, 
#' YearSpecific=YS)
#' 
#' compare_AIC(YearSpecific_ProportionsFirstOrder_Zero=YS1,
#' YearSpecific_ProportionsFirstOrder_Constant=YS1_cst)
#' 
#' compare_AIC(YearSpecific_ProportionsConstant=YS,
#'            YearSpecific_ProportionsFirstOrder=YS1,
#'            YearSpecific_ProportionsSecondOrder=YS2)
#'            
#' compare_AIC(YearSpecific_ProportionsFirstOrder=YS1_cst,
#'            YearSpecific_ProportionsSecondOrder=YS2_cst)
#' 
#' barplot_errbar(YS1_cst$proportions[1, ], y.plus = YS1_cst$proportions.CI.0.95[1, ], 
#'                y.minus = YS1_cst$proportions.CI.0.05[1, ], las=1, ylim=c(0, 0.7), 
#'                main="Proportion of the different rookeries in the region")
#' 
#' plot(cst, main="Use of different beaches along the time", what="total")
#' plot(expo, main="Use of different beaches along the time", what="total")
#' plot(YS2_cst, main="Use of different beaches along the time", what="total")
#' 
#' plot(YS1, main="Use of different beaches along the time")
#' plot(YS1_cst, main="Use of different beaches along the time")
#' plot(YS1_cst, main="Use of different beaches along the time", what="numbers")
#' 
#' parpre <- par(mar=c(4, 4, 2, 5)+0.4)
#' par(xpd=TRUE)
#' plot(YS, main="Use of different beaches along the time", 
#'      control.legend=list(x=2000, y=0.4, legend=c("Yalimapo", "Galibi", "Irakumpapy")))
#' par(mar=parpre)
#' 
#' # Example to modify order of series
#' plot(cst, order=c("Galibi.Suriname", "Irakumpapy.French.Guiana"))
#' plot(cst, order=c("Galibi.Suriname", "Irakumpapy.French.Guiana", "Yalimapo.French.Guiana"))
#' 
#' # Example to change the color
#' plot(cst, order=c("Galibi.Suriname", "Irakumpapy.French.Guiana", "Yalimapo.French.Guiana"), 
#'      col=function(n) rep(c("gray", "lightgrey"), floor(n/2)), border="black", 
#'      names.legend=c("Yalimapo", "Galibi", "Irakumpapy"))
#' }
#' @method plot fitRMU
#' @export

plot.fitRMU <- function (x, ..., resultMCMC=NULL, chain=1, replicate.CI=10000, CI.RMU=NULL, 
                         what = "proportions", criteria="50%", aggregate = "both", order = NULL, 
                         control.legend = list(), show.legend = TRUE, col=rainbow, border=NA, 
                         names.legend=NULL) 
{
  
  # resultMCMC=NULL; chain=1; replicate.CI=10000; CI.RMU=NULL; what = "proportions"; aggregate = "both"; order = NULL;control.legend = list(); show.legend = TRUE; col=rainbow; names.legend=NULL; border=NA
  
  p3p <- list(...) # p3p <- list()
  
  what <- tolower(what)
  aggregate <- tolower(aggregate)
  what <- match.arg(what, choices = c("proportions", "numbers", 
                                      "total"))
  aggregate <- match.arg(aggregate, choices = c("both", "model"))
  matobs <- as.matrix(x$RMU.data[, as.character(x$RMU.names[, 
                                                            "mean"])])
  
  # Pas bonne idée ; plutôt si replicate.CI != 0, utilise la médiane
  # if ((what == "proportions") | (what=="numbers")) replicate.CI <- 0
  
  if (((what == "proportions") | (what=="numbers")) & (replicate.CI != 0)) {
    replicate.CI <- 0
    warning("Median of proportion of beach usage will be shown.") 
  }
  
  if ((what == "total") & (is.null(x$hessian) & is.null(resultMCMC)) & (replicate.CI != 0)) {
    replicate.CI <- 0
    warning("No confidence interval can be estimated") 
  }
  
  if (!is.null(CI.RMU)) {
    outCI <- CI.RMU
  } else {
    # if (what != "total") replicate.CI <- 0
    outCI <- CI.RMU(result=x, 
                    resultMCMC=resultMCMC, 
                    chain=chain, 
                    replicate.CI=replicate.CI, 
                    silent=TRUE)
  }
  
  if (what == "total") {
    year <- as.numeric(colnames(outCI$Total))
    if (aggregate == "model") {
      y.minus <- outCI$Total["2.5%", ]
      y <- outCI$Total["50%", ]
      y.plus <- outCI$Total["97.5%", ]
    } else {
      y.minus <- outCI$Total_both["2.5%", ]
      y <- outCI$Total_both["50%", ]
      y.plus <- outCI$Total_both["97.5%", ]
    }
    
    if (x$model.trend == "year-specific") {
      yeap <- names(c(x$par, x$fixed.parameters.computing))[substr(names(c(x$par, x$fixed.parameters.computing)), 1, 2)=="T_"]
      yeap <- substr(yeap, 3, 6)
      y <- ifelse(names(y) %in% yeap, y, NA)
      y.plus <- ifelse(names(y.plus) %in% yeap, y.plus, NA)
      y.minus <- ifelse(names(y.minus) %in% yeap, y.minus, NA)
    }
    
    maxy <- max(c(y.plus, y), na.rm = TRUE)
    pp <- modifyList(list(ylim = c(0, maxy), xlab = "Year", 
                          ylab = "Nest number", las = 1, bty = "n",
                          main = "Total number of nests in all the rookeries in the region"), 
                     p3p)
    pp <- modifyList(pp, list(x = year, y = y, y.plus = y.plus, 
                              y.minus = y.minus))
    do.call(plot_errbar, pp)
  } else {
    # C'est soit proportion soit numbers
    
    if (what == "proportions") {
      outCI_f <- outCI$Proportions
    }
    if ((what == "numbers") & (aggregate == "model")) {
      outCI_f <- outCI$Numbers
    }
    if ((what == "numbers") & (aggregate == "both")) {
      outCI_f <- outCI$Numbers_both
    }
    
    if (!is.null(order)) {
      outCI_f <- outCI_f[, , order]
    }
    
    year <- as.numeric(names(outCI_f[1, , 1]))
    
    pp <- modifyList(list(ylim = c(0, max(rowSums(outCI_f[criteria, , ]), na.rm = TRUE)), 
                          xlab = "Year", ylab = what, las = 1, bty = "n"), 
                     p3p)
    pp <- modifyList(pp, list(x = year, y = outCI_f[criteria, , 1], type = "n"))
    do.call(plot, pp)
    p <- cbind(rep(0, nrow(outCI_f[criteria, , ])),
               rep(0, nrow(outCI_f[criteria, , ])))
    p <- cbind(p, outCI_f[criteria, , ])
    xx <- c(year, rev(year))
    ####### 22/6/2021
    col <- col(ncol(p)-1)
    for (i in 2:(ncol(p)-1)) {
      # print(i)
      y <- c(apply(p[, 1:i], 1, sum), rev(apply(p[, 1:(i + 1)], 1, sum)))
      polygon(xx[!is.na(y)], y[!is.na(y)], col = col[i - 1], border = border)
    }
    if (show.legend) {
      if (is.null(names.legend)) names.legend <- rev(names(outCI_f[1, 1, ]))
      plegend <- modifyList(list(x = "bottomright", legend = names.legend, 
                                 col = rev(col), pch = 15), 
                            control.legend)
      par(xpd=TRUE)
      do.call(legend, args = plegend)
      par(xpd=FALSE)
    }
  }
  return(invisible(outCI))
}