R/WACSestim.R

In WACS: Multivariate Weather-State Approach Conditionally Skew-Normal Generator

###################################################################
#
# This function is part of WACSgen V1.0
# Copyright © 2013,2014,2015, D. Allard, BioSP,
# and Ronan Trépos MIA-T, INRA
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details. http://www.gnu.org
#
###################################################################
  
  
#' Estimation of the parameters of a WACS model
#'
#' @importFrom mclust Mclust
#' @importFrom mclust densityMclust
#' @importFrom mclust plotDensityMclustd
#' @importFrom mclust mclustBIC
#' @import mnormt
#' @import tmvtnorm
#' 
#' @export
#' 
#' @param wacsdata   Data, as returned by WACSdata
#' @param spar       Smoothing parameter for estimating annual cycle
#' @param trend.norm Type of norm used in for computing central tendency
#'                   and variation. Must be  \code{"L1"} or \code{"L2"}.
#' @param rain.model Model for precipitation. Must be \code{"Gamma"} or \code{"None"}
#' @param method     \code{"MLE"} or \code{"MOM"}. Estimation method for the rain model.
#' @param Vsel       Variables (other than rain) on which clustering is 
#'                   performed when \code{Vsel=NULL}, all variables are considered.
#' @param Nclusters  Number of clusters to consider. When \code{Nclusters = NULL},
#'                   absolute best clustering is sought for wet and dry weather
#'                   states in each season (up to 4).
#' @param clustering Indicates whether \code{"hard"} or \code{"soft"} clustering is
#'                   considered.
#' @param plot.it    Boolean indicating whether plots are produced                  
#' @param DIR        Directory in which placing plot
#'                  
#'                  
#' @return A list containing all parameters; see the user guide for details.
#' 
#' 
#' 
#' @examples
#' \dontrun{
#'
#'  ## For an estimation with default setting 
#'  ThisPar  = WACSestim(ThisData)
#'
#'  ## For an estimation with max. 2 dry and wet weather types per season, 
#'  ## and production of plots
#'  ThisPar  = WACSestim(ThisData, Nclusters = 1:2, plot.it = TRUE) 
#'
#'  ## For an estimation with exactly 2 dry and wet weather states per season, 
#'  ## clustering on variables 3 and 5 only and no production of plots
#'  ThisPar  = WACSestim(ThisData, Nclusters = 2, Vsel = c(3,5)) 
#'  }
#'
#' 
#' @note
#' 
#' Larger values of \code{spar} produce smoother estimates. Smaller values produce less smooth estimates. \code{spar=0.7} 
#' is a good compromise
#' 
#' 
#' Soft clustering means that days have probabilities to belong to each weather state. With hard clustering, 
#' this probability is set to 1 to the most likely weather state and 0 to all others. Density parameter estimates 
#' are more robust with \code{clustering="soft"}. Clustering is done by means of the mclust package with 
#' \env{modelNames="VVV"}


WACSestim=function(wacsdata,spar=0.7,trend.norm="L2",rain.model="Gamma",method="MLE",
                     Vsel=NULL,Nclusters=NULL,clustering="soft",plot.it=FALSE,DIR="./"){

  # Checking wacsdata
  if (class(wacsdata) != "WACSdata") {
    stop ("[WACSestim] Data should be of class 'WACSdata', as generated by calling WACSdata")
  }
  
  DATA=wacsdata$data
    
  # Checking if the rain model has been defined
  if ( (rain.model != "Gamma") & (rain.model != "None")){
    stop ("[WACSestim] rain model should either be 'Gamma' or 'None'")
  }
  # Checking is season has been previously defined
  if (! ("season" %in% names(DATA))) {
    stop("[WACSestim] 'season' not found (maybe you forgot to use WACSdata)")
  }
  # Checking if type of clustering has been defined
  if ( (clustering != "hard") & (clustering != "soft")) {
    stop("[WACSestim] clustering must be 'hard' or 'soft' ")
  }
  #Collect names of the variables  
  vdn   = c("year","month","day","season") #date variables
  vvn   = setdiff(names(DATA), vdn) #var names 
  vvorn = setdiff(vvn, c("rain"))  #var names without rain
  
  ###############################################
  # Some intialisations and checking
  # Creating : Ns,Nd,Ny,Nv = number of seasons, days, years and variables in addition to rain
  #            seqS        = vector of length Nd containing the season (from 1 to Ns)
  #            OccR        = vector of booleans of length Nd. Is TRUE if rain is > 0                       
  ###############################################
  
  # Checking if the number of days is a multiple of number of years
  Nd = nrow(DATA)            # number of days
  Ny = Nd%/%365                 # number of years -- must be an integer
  
  if (Nd != Ny*365) {
    stop("[WACSestim] length of data should be multiple of 365 (maybe you forgot to use WACSdata)")
  }
  
  Nv=length(vvorn);         # mumber of variables other than rain
  Ns       = length(unique(DATA$season)); # number of seasons
  NumbDays = rep(0,Ns);
  for (s in 1:Ns) {
    NumbDays[s] = length(which(DATA$season == s))/Ny
  }



  
  ###############################################
  # Estimating and removing the annual cycle to all variables except rain
  #
  # Creating:  par.gamma = parameters of the gamma transform (scale and shape)
  #            NscTR     = vector of normal score transform. Equal to -999 if no rain
  #            AnnualTrend = to be saved
  #            Resid     = matrix of size (Nd x Nv) of residuals of variables
  #    
  # Calls :    wacs.estimcycle()
  #              
  ###############################################
  Ds        = wacs.estimcycle(DATA[vvorn],spar=spar,trend.norm=trend.norm,plot.it=plot.it,DIR=DIR)
  
  Resid     = Ds[[1]]
  Central   = Ds[[2]]   
  colnames(Central) = vvorn
  Deviation = Ds[[3]]
  colnames(Deviation) = vvorn
  par.trend = list(Param=c(spar,trend.norm),Central=Central,Deviation=Deviation)
  print("Trend removed")
  rm(Ds)
  
  ###############################################
  # Estimating rain model and transforming rain 
  #
  # Creating:  par.rain  = list containing rain model and its parameters 
  #            Ns.rain   = vector of normal score transform. Equal to -999 if no rain
  # Calls :    wacs.estimrain()         
  ###############################################

  
  Er       = wacs.estimrain(DATA,rain.model=rain.model,method=method,plot.it=plot.it,DIR=DIR)
  par.rain = list(RainModel=rain.model,RainPar=Er[[1]])
  NS.rain  = Er[[2]]
  print("Rain transformed")
  rm(Er)
  
  ###############################################
  #
  # Entering a loop on seasons
  #            
  ###############################################
  
  year.min        = min(DATA$year)
  julianday.month = c(0,31,59,90,120,151,181,212,243,273,304,334)
  par.seasons     = list()
  
  options(warn=-1)
  for (s in 1:Ns){
    season.name = paste("Season_",s,sep="") 
    print(season.name)
    Seasonlist  = list(NumbDays=NumbDays[s])
    day         = (DATA$year[DATA$season==s] - year.min)*365 + (julianday.month[DATA$month[DATA$season==s]]) + DATA$day[DATA$season==s]
    y           = cbind(NS.rain[DATA$season==s],Resid[DATA$season==s,])

    #############################################
    # Estimating the clusters for each season
    #
    # Creating:  WT.z: probability of belonging to each weather type
    #
    # Calls :    wacs.estimWT()            
    #############################################

    eWT        = wacs.estimWT(y,season=s,Vsel=Vsel,Nclusters=Nclusters,plot.it=plot.it,DIR=DIR)
    WT.z       = eWT$Prob
    NumbWT     = eWT$NumbWT
    WT         = eWT$WT
    Seasonlist = c(Seasonlist,list(NumbWT=NumbWT))
    print(paste("       Weather types estimated:",NumbWT[1],"dry WT;  ",NumbWT[2],"wet WT"),sep="")


    #############################################
    #
    # Estimating the transition matrix for each season
    #
    # Creating:  TransM: the transition matrix
    #
    # Calls :    wacs.estimMarkov()
    #            
    #############################################  
    
    TransM     = wacs.estimMarkov(WT.z,clustering=clustering)
    Seasonlist = c(Seasonlist,list(TransM=TransM))
    print("       Markov Transition estimated")

    #############################################
    #
    # Estimating the parameters of the CSN for each season
    #
    # Creating: par.csn = a list containing, for each weather state,
    #                     the estimated parameters (location, covariance, skewness) 
    #                     of the 1+Nv variables
    #                     
    #    
    # Calls:    wacs.estimCSNstar which does the actual estimation 
    #            
    #############################################

    par.csn=list()
    Nwt.dry = NumbWT[1]
    Nwt.wet = NumbWT[2] 
    DRY =    (y[,1] == -999)
    
    # Dry weather states    
    for (wt in 1:Nwt.dry){
      weight=rep(1,dim(WT.z)[1])
      if(clustering == "soft") weight = weight*WT.z[,wt]
      if(clustering == "hard") weight = weight*as.integer(WT==wt) 
      WT.name=paste("WT_",wt,sep="")
      print(paste("        ",WT.name),sep="")
      par.csn[[paste("W",wt,sep="")]] = wacs.estimCSNstar(y[DRY,-1],weight[DRY],day[DRY])
    }

    # Wet weather states  
    for (wt in (Nwt.dry+1):(Nwt.dry+Nwt.wet)){
      weight=rep(1,dim(WT.z)[1])
      if(clustering == "soft") weight = WT.z[,wt]
      if(clustering == "hard") weight = weight*as.integer(WT==wt)
      WT.name=paste("WT_",wt,sep="")
      print(paste("        ",WT.name),sep="")
      par.csn[[paste("W",wt,sep="")]] = wacs.estimCSNstar(y[!DRY,],weight[!DRY],day[!DRY])
    }
    Seasonlist      = c(Seasonlist,par.csn)    
    par.seasons[[season.name]] = Seasonlist
  }

  options(warn=0)
  res = c(list(bounds=wacsdata$bounds, mapping=wacsdata$mapping, 
               Trend=par.trend,Rain=par.rain,
               seasons=wacsdata$seasons,Trange=wacsdata$Trange,varnames=vvn),
          par.seasons)
  class(res) = "WACSpar"
  return(res)
}