R: Standardized Climate Index (SCI)

fitSCI {SCI}

R Documentation

Standardized Climate Index (SCI)

Description

fitSCI identifies parameters for the Standardized Climate Index (SCI) transformation. transformSCI applies the transformation

Usage

fitSCI(x, ...)
## Default S3 method:
fitSCI(x, first.mon, time.scale, distr, p0, 
       p0.center.mass=FALSE, scaling=c("no","max","sd"),mledist.par =  list(),
       start.fun = dist.start, start.fun.fix = FALSE, warn = TRUE, ...)

transformSCI(x, ...)
## Default S3 method:
transformSCI(x, first.mon, obj, sci.limit = Inf, warn=TRUE, ...)

Arguments

`x`	`numeric` vector, representing a monthly univariate time series.
`first.mon`	value in 1:12 indicating the month of the first element of x
`time.scale`	The time scale (`integer`) of the SCI calculation. The time scale is the window length of an backward looking running mean.
`distr`	A character string `"name"` naming a distribution for which the corresponding density function (`dname`), the corresponding distribution function (`pname`) and the quantile function (`qname`) must be defined (see for example `GammaDist`)
`p0`	if TRUE, model Probability of zero (precipitation) months is modeled with a mixed distribution as `D(x) = p_0 + (1-p_0)G(x)`, where `G(x) > 0` is the reference distribution (e.g. Gamma) `p0` is the probability of a zero (precipitation) month.
`p0.center.mass`	If TRUE, the Probability of zero (precipitation) is estimated using a "center of mass" estimate based on the Weibull plotting position function (see details). Only applies if `p0=TRUE`.
`scaling`	Indicates whether to do some scaling of `x` prior to parameter identification. `"no"` (the default) indicates no scaling. `"max"` indicates scaling by the maximum of `x`, such that `x <- x/max(x,na.rm=TRUE)`. `"sd"` stands for scaling by the standard deviation. Scaling can stabilize parameter estimation.
`mledist.par`	named `list` that can be used to pass parameters to `mledist` in package fitdistrplus.
`start.fun`	Function with arguments `x` and `distr` estimating initial parameters of the function `distr` for each month. The function should return a named list corresponding to the parameters of `distr`. (See also `dist.start`)
`start.fun.fix`	`logical` argument, indicating if parameter estimates of `start.fun` should be used if maximum likelihood estimation breaks down. This stabilizes the implementation but can introduce biases in the resulting SCI.
`obj`	an object of class `fitSCI`, output from `fitSCI`.
`sci.limit`	Truncate absolute values of SCI that are lage than `sci.limit`. See details.
`warn`	Issue warnings if problems in parameter estimation occur.
`...`	further arguments passed to methods

Details

fitSCI estimates the parameters for transforming a meteorological and environmental time series to a Standardized Climate Index (SCI). transformSCI applies the standardisation. Typical SCI are the Standardized Precipitation Index (SPI), the Standardized Runoff Index (SRI) or the Standardized Precipitation Evapotranspiration Index (SPEI).

To reduce biases in the presence of many zero (precipitation) events, the probability of these events (p_0) can be estimated using a "center of mass" estimate based on the Weibull plotting position function (p0.center.mass=TRUE). Following Stagge et al. (2014) the probability of zero events is then estimated as p_0 = \frac{n_p}{n + 1}, where n_p refers to the number of zero events and n is the sample size. The resulting mixed distribution used fro SCI transformation is then

D(x) = \left\{ \begin{array}{l l} p_0 + (1-p_0) G(x) & \quad \mbox{if } x > 0 \\ \frac{n_p + 1}{2(n+1)} & \quad \mbox{if } x = 0 \end{array} \right.

where G(x)>0 is a model (e.g. gamma) distribution.

Uncertainty in distribution parameters can cause unrealistically large (small) SCI values if values in x exceed the values used for parameter estimation (see fitSCI). Therefore transformSCI allows for a truncation of the SCI series such that abs(sci)<=sci.limit. The truncation can be disabled by setting sci.limit=Inf.

Value

fitSCI returns an object of class "fitSCI" with the following components:

`dist.para`	A column `matrix` containing the parameters of distribution `distr` for each month. Row names correspond to the distribution parameters. If `p0=TUE` an additional row named `P0` is introduced, indicating the probability of zero (precipitation) events.
`dist.para.flag`	an vector indicating possible issues occurring throughout parameter estimation. Possible values are: 0. no problems occurred; 1. starting values could not be estimated; 2. `mledist` crashed with unknown error; 3. `mledist` did not converge; 4. all values in this month are `NA`; 5. all values in this month are constant, distribution not defined.
`time.scale`	The time scale (`integer`) of the SCI calculation.
`distr`	A character string `"name"` naming a distribution used
`p0`	`logical` indicating whether probability of zero (precipitation) events is estimated separately.
`p0.center.mass`	`logical` indicating whether probability of zero (precipitation) events is estimated using the "centre of mass" estimator (see Stagge et al. (2014) for details).
`scaling`	`numeric` value that has been used to scale `x` (see argument `scaling`). A value of 1 results from `scaling="no"`, other values are the maximum value or the standard deviation of `x`, depending on the choice of the parameter `scaling`.
`call`	the function call

transformSCI returns a numeric vector containing the SCI, having values of the standard normal distribution.

Note

This function is intended to be used together with transformSCI.

Author(s)

Lukas Gudmundsson & James Stagge

References

Stagee, J.H. ; Tallaksen, L.M.; Gudmundsson, L.; van Loon, A.; Stahl, K.: Candidate Distributions for Climatological Drought Indices (SPI and SPEI), 2015, International Journal of Climatology, 35, 4027-4040, doi:10.1002/joc.4267.

Stagee, J.H. ; Tallaksen, L.M.; Gudmundsson, L.; van Loon, A.; Stahl, K.: Response to comment on "Candidate Distributions for Climatological Drought Indices (SPI and SPEI)", 2016, International Journal of Climatology, 36, 2132-2138, doi:10.1002/joc.4564.

McKee, T.; Doesken, N. & Kleist, J.: The relationship of drought frequency and duration to time scales Preprints, 8th Conference on Applied Climatology, 1993, 179-184.

Shukla, S. & Wood, A. W.: Use of a standardized runoff index for characterizing hydrologic drought Geophysical Research Letters, 2008, 35, L02405.

Vicente-Serrano, S. M.; Begueria, S. & Lopez-Moreno, J. I.: A Multiscalar Drought Index Sensitive to Global Warming: The Standardized Precipitation Evapotranspiration Index J. Climate, Journal of Climate, American Meteorological Society, 2009, 23, 1696-1718.

Examples

##
## generate artificial data
##
set.seed(101)
n.years <- 60
date <- rep(1:n.years,each=12) + 1950 + rep((0:11)/12,times=n.years)
## Precipitation
PRECIP <- (0.25*sin( 2 * pi * date) + 0.3)*rgamma(n.years*12, shape = 3, scale = 1)
PRECIP[PRECIP<0.1] <- 0
## Potential Evapotranspiration
PET <- 0.5*sin( 2 * pi * date)+1.2+rnorm(n.years*12,0,0.2)
## display test data
matplot(date,cbind(PRECIP,PET),t=c("h","l"),col=c("blue","red"),lty=1)
legend("topright",legend=c("PRECIPitation","temperature"),fill=c("blue","red"))

##
## example SPI
##
spi.para <- fitSCI(PRECIP,first.mon=1,distr="gamma",time.scale=6,p0=TRUE)
spi.para
spi <- transformSCI(PRECIP,first.mon=1,obj=spi.para)
plot(date,spi,t="l")

##
## effect of time.scale on SPI
##
spi.1.para <- fitSCI(PRECIP,first.mon=1,time.scale=1,distr="gamma",p0=TRUE)
spi.12.para <- fitSCI(PRECIP,first.mon=1,time.scale=12,distr="gamma",p0=TRUE)
spi.1 <- transformSCI(PRECIP,first.mon=1,obj=spi.1.para)
spi.12 <- transformSCI(PRECIP,first.mon=1,obj=spi.12.para)
matplot(date,cbind(spi.1,spi.12),t="l",lty=1,col=c("red","blue"),lwd=c(1,2))
legend("topright",legend=c("time.scale=1","time.scale=12"),fill=c("red","blue"))

##
## example SPEI
##
if(require(evd)){
    spei.para <- fitSCI(PRECIP-PET,first.mon=1,time.scale=6,distr="gev",p0=FALSE)
    spei <- transformSCI(PRECIP-PET,first.mon=1,obj=spei.para)
    plot(date,spei,t="l")
}

##
## effect of changing different distribution for SPEI computation
##
spei.genlog.para <- fitSCI(PRECIP-PET,first.mon=1,time.scale=6,distr="genlog",p0=FALSE)
spei.genlog <- transformSCI(PRECIP-PET,first.mon=1,obj=spei.genlog.para)
if(require(evd)){lines(date,spei.genlog, col="red")} else {plot(date,spei.genlog,t="l")}
## in this case: only limited effect.
## generally: optimal choice of distribution: user responsibility.

##
## use a 30 year reference period for SPI parameter estimation
##
sel.date <- date>=1970 & date<2000
spi.ref.para <- fitSCI(PRECIP[sel.date],first.mon=1,distr="gamma",time.scale=6,p0=TRUE)
## apply the the parameters of the reference period to all data
## also outside the reference period
spi.ref <- transformSCI(PRECIP,first.mon=1,obj=spi.ref.para)
plot(date,spi.ref,t="l",col="blue",ylim=c(-5,5),lwd=2)
lines(date[sel.date],spi.ref[sel.date],col="red",lwd=3)
legend("bottom",legend=c("reference period","extrapolation period"),fill=c("red","blue"),
       horiz=TRUE)

##
## use "start.fun.fix" in instances where maximum likelyhood estimation fails
##
## force failure of maximum likelyhood estimation by adding "strange" value
## a warning should be issued
xx <- PRECIP-PET; xx[300] <- 1000
spei.para <- fitSCI(xx,first.mon=2,time.scale=1,p0=FALSE,distr="gev")
spei.para$dist.para
## use start.fun, usually ment for estimating inital values for
## parameter optimisation if maximum likelihood estimation fails
spei.para <- fitSCI(xx,first.mon=2,time.scale=1,p0=FALSE,distr="gev",
                    start.fun.fix=TRUE)
spei.para$dist.para

##
## usage of sci.limit to truncate unrealistic SCI values
## 
PRECIP.mod <- PRECIP
PRECIP.mod[300] <- 100 ## introduce spuriously large value
spi.mod.para <- fitSCI(PRECIP.mod,first.mon=1,time.scale=3,p0=TRUE,distr="gamma")
plot(transformSCI(PRECIP.mod,first.mon=1,obj=spi.mod.para,sci.limit=Inf),
     t="l",col="blue",lwd=2)
lines(transformSCI(PRECIP.mod,first.mon=1,obj=spi.mod.para,sci.limit=4),col="red")

##
## how to modify settings of function "mledist" used for parameter identification
## 
## identify parameters with standard settings
spi.para <- fitSCI(PRECIP,first.mon=1,distr="gamma",time.scale=6,p0=TRUE)
## add lower and upper limits for parameter identification
lower.lim <- apply(spi.para$dist.para,1,min) - 0.5*apply(spi.para$dist.para,1,sd)
upper.lim <- apply(spi.para$dist.para,1,max) + 0.5*apply(spi.para$dist.para,1,sd)
spi.para.limit <- fitSCI(PRECIP,first.mon=1,distr="gamma",time.scale=6,p0=TRUE,
                         mledist.par=list(lower=lower.lim, upper=upper.lim))

##
## how to write an own start.fun
## (required if distributions not mentioned in "dist.start" are used)
##
## function with same arguments as "dist.start" 
my.start <- function(x,distr="gamma"){
### code based on "mmedist" in package "fitdistrplus"
    ppar <- try({
        n <- length(x)
        m <- mean(x)
        v <- (n - 1)/n * var(x)
        shape <- m^2/v
        rate <- m/v
        list(shape = shape, rate = rate)},TRUE)
    if (class(ppar) == "try-error") ## function has to be able to return NA parameters
        ppar <- list(shape = NA, rate = NA)
    return(ppar)
}
my.start(PRECIP)
spi.para <- fitSCI(PRECIP,first.mon=1,time.scale=6,p0=TRUE,
                   distr="gamma",start.fun=my.start)

[Package SCI version 1.0-2 Index]