gevdfit.R 30 KB
Newer Older
1
2
3
# This file contains the functions:
# - gev.d.fit, gev.d.init for fitting
# - gev.d.diag for diagnostic plots
4
# - gev.d.params for calculation of parameters
5
6
7
8
# and the documentation of the example data

#### gev.d.fit ####

9
#' @title Maximum-likelihood Fitting of the duration-dependent GEV Distribution
Laura Mack's avatar
Laura Mack committed
10
11
12
#' @description Modified \code{\link[ismev]{gev.fit}} function for Maximum-likelihood fitting
#' for the duration-dependent generalized extreme
#' value distribution, following Koutsoyiannis et al. (1998), including generalized linear
13
#' modeling of each parameter.
Laura Mack's avatar
Laura Mack committed
14
#' @param xdat A vector containing maxima for different durations.
15
#' This can be obtained from \code{\link{IDF.agg}}.
Laura Mack's avatar
Laura Mack committed
16
#' @param ds A vector of aggregation levels corresponding to the maxima in xdat.
17
#' 1/60 corresponds to 1 minute, 1 corresponds to 1 hour.
Laura Mack's avatar
Laura Mack committed
18
19
#' @param ydat A matrix of covariates for generalized linear modeling of the parameters
#' (or NULL (the default) for stationary fitting). The number of rows should be the same as the
20
#' length of xdat.
21
#' @param  mutl,sigma0l,xil,thetal,etal,taul,eta2l Numeric vectors of integers, giving the columns of ydat that contain
Laura Mack's avatar
Laura Mack committed
22
#'  covariates for generalized linear modeling of the parameters (or NULL (the default)
23
#'  if the corresponding parameter is stationary).
24
#'  Parameters are: modified location, scale offset, shape, duration offset, duration exponent, respectively.
25
#' @param mutlink,sigma0link,xilink,thetalink,etalink,eta2link,taulink Link functions for generalized linear
26
#' modeling of the parameters, created with \code{\link{make.link}}. The default is \code{make.link("identity")}.
27
28
29
30
#' @param init.vals list, giving initial values for all or some parameters
#' (order: mut, sigma0, xi, theta, eta, eta2, tau). If one of those parameters shall not be used (see theta_zero, eta2_zero, tau_zero),
#' the number of parameters has to be reduced accordingly. If some or all given values in init.vals are NA or 
#' no init.vals at all is declared (the default), initial parameters are obtained
Laura Mack's avatar
Laura Mack committed
31
#' internally by fitting the GEV separately for each duration and applying a linear model to obtain the
32
#' duration dependency of the location and shape parameter.
33
#' Initial values for covariate parameters are assumed as 0 if not given.
Felix Fauer's avatar
Felix Fauer committed
34
35
#' @param theta_zero Logical value, indicating whether theta should be estimated (FALSE, the default) or
#' should stay zero.
36
#' @param tau_zero,eta2_zero Logical values, indicating whether tau,eta2 should be estimated (TRUE, the default).
37
38
39
40
#' @param show Logical; if TRUE (the default), print details of the fit.
#' @param method The optimization method used in \code{\link{optim}}.
#' @param maxit The maximum number of iterations.
#' @param ... Other control parameters for the optimization.
Laura Mack's avatar
Laura Mack committed
41
42
43
44
45
46
47
#' @return A list containing the following components.
#' A subset of these components are printed after the fit.
#' If \code{show} is TRUE, then assuming that successful convergence is indicated,
#' the components nllh, mle and se are always printed.
#' \item{nllh}{single numeric giving the negative log-likelihood value}
#' \item{mle}{numeric vector giving the MLE's for the modified location, scale_0, shape,
#' duration offset and duration exponent, resp.}
48
#' \item{se}{numeric vector giving the standard errors for the MLE's (in the same order)}
49
#' \item{trans}{A logical indicator for a non-stationary fit.}
Laura Mack's avatar
Laura Mack committed
50
#' \item{model}{A list with components mutl, sigma0l, xil, thetal and etal.}
51
#' \item{link}{A character vector giving inverse link functions.}
Laura Mack's avatar
Laura Mack committed
52
#' \item{conv}{The convergence code, taken from the list returned by \code{\link{optim}}.
53
#' A zero indicates successful convergence.}
Laura Mack's avatar
Laura Mack committed
54
55
#' \item{data}{data is standardized to standard Gumbel.}
#' \item{cov}{The covariance matrix.}
56
#' \item{vals}{Parameter values for every data point.}
57
#' \item{init.vals}{Initial values that were used.}
58
#' \item{ds}{Durations for every data point.}
59
60
#' @details For details on the d-GEV and the parameter definitions, see \link{IDF-package}.
#' @seealso \code{\link{IDF-package}}, \code{\link{IDF.agg}}, \code{\link{gev.fit}}, \code{\link{optim}}
61
#' @export
Laura Mack's avatar
Laura Mack committed
62
63
64
65
#' @importFrom stats optim
#' @importFrom stats make.link
#'
#' @examples
66
67
#' # sampled random data from d-gev with covariates
#' # GEV parameters:
68
69
#' # mut = 4 + 0.2*cov1 +0.5*cov2
#' # sigma0 = 2+0.5*cov1
70
71
72
#' # xi = 0.5
#' # theta = 0
#' # eta = 0.5
73
74
#' # eta2 = 0.5
#' # tau = 0
Laura Mack's avatar
Laura Mack committed
75
#'
76
#' data('example',package ='IDF')
Laura Mack's avatar
Laura Mack committed
77
#'
78
#' gev.d.fit(xdat=example$dat,ds = example$d,ydat=as.matrix(example[,c('cov1','cov2')])
79
#' ,mutl=c(1,2),sigma0l=1)
80
81

gev.d.fit<-
82
  function(xdat, ds, ydat = NULL, mutl = NULL, sigma0l = NULL, xil = NULL, thetal = NULL, etal = NULL, taul = NULL, eta2l = NULL,
83
           mutlink = make.link("identity"), sigma0link = make.link("identity"), xilink = make.link("identity"),
84
           thetalink = make.link("identity"), etalink = make.link("identity"), taulink = make.link("identity"), eta2link = make.link("identity"),
85
           init.vals = NULL, theta_zero = FALSE, tau_zero=TRUE, eta2_zero=TRUE,
86
           show = TRUE, method = "Nelder-Mead", maxit = 10000, ...)
87
  {
88
89
    if (length(xdat) != length(ds)) {
      stop(paste0('The length of xdat is ',length(xdat),', but the length of ds is ',length(ds),'.'))
Laura Mack's avatar
Laura Mack committed
90
    }
91
    z <- list()
Laura Mack's avatar
Laura Mack committed
92
    # number of parameters (betas) to estimate for each parameter:
93
94
95
    npmu <- length(mutl) + 1
    npsc <- length(sigma0l) + 1
    npsh <- length(xil) + 1
96
    npth <- ifelse(!theta_zero,length(thetal) + 1,0)
97
    npet <- length(etal) + 1
98
    npta <- ifelse(!tau_zero,  length(taul)   + 1,0)
99
    npe2 <- ifelse(!eta2_zero,  length(eta2l)   + 1,0)
100
    z$trans <- FALSE  # indicates if fit is non-stationary
101
102
    z$model <- list(mutl, sigma0l, xil, thetal, etal, eta2l, taul)
    z$link <- list(mutlink=mutlink, sigma0link=sigma0link, xilink=xilink, thetalink=thetalink, etalink=etalink, eta2link=eta2link, taulink=taulink)
Laura Mack's avatar
Laura Mack committed
103

104
105
    # test for NA values:
    if(any(is.na(xdat))) stop('xdat contains NA values. NA values need to be removed first.')
106
107
    # test for finite values:
    if(any(is.infinite(xdat))) stop('xdat contains non finite values. Inf and -Inf need to be removed first.')
Laura Mack's avatar
Laura Mack committed
108

109
110
    # test if covariates matrix is given correctly
    npar <- max(sapply(z$model,function(x){return(ifelse(is.null(x),0,max(x)))}))
111
    if(any(npar>ncol(ydat),npar>0 & is.null(ydat)))stop("Not enough columns in covariates matrix 'ydat'.")
Laura Mack's avatar
Laura Mack committed
112

113
    # initial values
114
    init.necessary.length = 4 + as.numeric(!theta_zero) + as.numeric(!eta2_zero) + as.numeric(!tau_zero)  # mut, sigma0, xi, theta, eta, eta2, tau
115
116
117
118
119
    if (is.null(init.vals)) {init.vals = as.list(rep(NA,init.necessary.length))
    }else{init.vals = as.list(init.vals)}
    
    if(length(init.vals)!=init.necessary.length | !is.list(init.vals)) {
      warning(paste0('Parameter init.vals is not used, because it is no list of length ',init.necessary.length,'.'))
Felix Fauer's avatar
Felix Fauer committed
120
121
122
123
124
125
126
127
      init.vals <- gev.d.init(xdat,ds,z$link)
      
    }else{ # length of given values is correct

      # name given initial values
      names1=c('mu','sigma','xi')                 # standard set of parameters
      if (!theta_zero){names1=c(names1,'theta')}  # add theta (in case)
      names1=c(names1,'eta')                      # add eta   (always)
128
      if (!eta2_zero){names1=c(names1,'eta2')}    # add eta2  (in case)
Felix Fauer's avatar
Felix Fauer committed
129
130
131
132
      if (!tau_zero){names1=c(names1,'tau')}      # add tau   (in case)
      names(init.vals) <- names1
      # add missing initial value (fixed internal number of parameters: 7)
      if (theta_zero) init.vals$theta = 0
133
      if (eta2_zero) init.vals$eta2 = init.vals$eta
Felix Fauer's avatar
Felix Fauer committed
134
      if (tau_zero) init.vals$tau = 0
135
      
Felix Fauer's avatar
Felix Fauer committed
136
137
138
      if(!any(is.na(init.vals))){ #all initial values are given
        # do nothing
      }else if(any(!is.na(init.vals))) { #some initial values are given
139
        #if (!eta2_zero) print("autmoatic inital value setting not implemented yet for multiscaling (eta2_zero=FALSE)")
Felix Fauer's avatar
Felix Fauer committed
140
141
142
143
144
145
        init.vals.user <- init.vals
        init.vals <- gev.d.init(xdat,ds,z$link) #calculate init.vals using gev.d.init
        for (i in 1:length(init.vals)){ #overwrite the calculated initial values with the ones given by the user
          if(!is.na(init.vals.user[[names(init.vals.user)[i]]])) {
            init.vals[[names(init.vals.user)[i]]]<-init.vals.user[[names(init.vals.user)[i]]]
          } 
146
        }
Felix Fauer's avatar
Felix Fauer committed
147
      }else{ #no initial values are given
148
        #if (!eta2_zero) print("autmoatic inital value setting not implemented yet for multiscaling (eta2_zero=FALSE)")
Felix Fauer's avatar
Felix Fauer committed
149
        init.vals <- gev.d.init(xdat,ds,z$link)
150
      }
Felix Fauer's avatar
Felix Fauer committed
151
    } 
Laura Mack's avatar
Laura Mack committed
152
153

    # generate covariates matrices:
154
    if (is.null(mutl)) { #stationary
155
      mumat <- as.matrix(rep(1, length(xdat)))
156
      muinit <- init.vals$mu
157
    }else { #non stationary
158
      z$trans <- TRUE
159
160
      mumat <- cbind(rep(1, length(xdat)), ydat[, mutl])
      muinit <- c(init.vals$mu, rep(0, length(mutl)))[1:npmu] #fill with 0s to length npmu
161
    }
162
    if (is.null(sigma0l)) {
163
      sigmat <- as.matrix(rep(1, length(xdat)))
164
      siginit <- init.vals$sigma
165
166
    }else {
      z$trans <- TRUE
167
168
      sigmat <- cbind(rep(1, length(xdat)), ydat[, sigma0l])
      siginit <- c(init.vals$sigma, rep(0, length(sigma0l)))[1:npsc]
169
    }
170
    if (is.null(xil)) {
171
      shmat <- as.matrix(rep(1, length(xdat)))
Laura Mack's avatar
Laura Mack committed
172
      shinit <- init.vals$xi
173
174
    }else {
      z$trans <- TRUE
175
176
      shmat <- cbind(rep(1, length(xdat)), ydat[, xil])
      shinit <- c(init.vals$xi, rep(0, length(xil)))[1:npsh]
177
178
179
    }
    if (is.null(thetal)) {
      thmat <- as.matrix(rep(1, length(xdat)))
180
      thetainit <- init.vals$theta
181
182
183
    }else {
      z$trans <- TRUE
      thmat <- cbind(rep(1, length(xdat)), ydat[, thetal])
184
      thetainit <- c(init.vals$theta, rep(0, length(thetal)))[1:npth]
185
186
187
    }
    if (is.null(etal)) {
      etmat <- as.matrix(rep(1, length(xdat)))
188
      etainit <- init.vals$eta
189
190
191
    }else {
      z$trans <- TRUE
      etmat <- cbind(rep(1, length(xdat)), ydat[, etal])
192
      etainit <- c(init.vals$eta, rep(0, length(etal)))[1:npet]
193
    }
Felix Fauer's avatar
Felix Fauer committed
194
195
196
197
198
199
200
    if (is.null(taul)) {
      tamat <- as.matrix(rep(1, length(xdat)))
      tauinit <- init.vals$tau
    }else {
      z$trans <- TRUE
      tamat <- cbind(rep(1, length(xdat)), ydat[, taul])
      tauinit <- c(init.vals$tau, rep(0, length(taul)))[1:npta]
201
    }
202
203
204
205
206
207
208
209
210
    if (is.null(eta2l)) {
      e2mat <- as.matrix(rep(1, length(xdat)))
      eta2init <- init.vals$eta2
    }else {
      z$trans <- TRUE
      e2mat <- cbind(rep(1, length(xdat)), ydat[, eta2l])
      eta2init <- c(init.vals$eta2, rep(0, length(eta2l)))[1:npe2]
    }
    
211
    init <- c(muinit,siginit,shinit)
Felix Fauer's avatar
Felix Fauer committed
212
213
    if (!theta_zero) {init <- c(init,thetainit)} # add theta init (in case)
    init <- c(init,etainit)                      # add eta init   (always)
214
    if (!eta2_zero)  {init <- c(init,eta2init)}  # add eta2 init  (in case)
Felix Fauer's avatar
Felix Fauer committed
215
    if (!tau_zero)   {init <- c(init,tauinit)}   # add tau init   (in case)
216
     
217
218
219
    # function to calculate neg log-likelihood:
    gev.lik <- function(a) {
      # computes neg log lik of d-gev model
220
221
222
      mu <- mutlink$linkinv(mumat %*% (a[1:npmu]))
      sigma <- sigma0link$linkinv(sigmat %*% (a[seq(npmu + 1, length = npsc)]))
      xi <- xilink$linkinv(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)]))
223
      # Next line will set the theta likelihood as non-existent in case user requested it. (same for tau)
224
      if(!theta_zero) {theta <- thetalink$linkinv(thmat %*% (a[seq(npmu + npsc + npsh + 1, length = npth)]))}
225
      eta <- etalink$linkinv(etmat %*% (a[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
226
227
      if(!eta2_zero)  {eta2  <- eta2link$linkinv( e2mat %*% (a[seq(npmu + npsc + npsh + npth + npet + 1, length = npe2)]))}
      if(!tau_zero)   {tau   <- taulink$linkinv(  tamat %*% (a[seq(npmu + npsc + npsh + npth + npet + npe2 + 1, length = npta)]))}
228
      
229
      ifelse(!theta_zero, ds.t <- ds+theta, ds.t <- ds) #Don't use theta if user requested not to have it.
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
      #ifelse(!tau_zero,   sigma.d <- sigma/(ds.t^eta)+tau, sigma.d <- sigma/(ds.t^eta)) # don't use tau if user requested not to have it.
      if (tau_zero){ # don't use tau if user requested not to have it.
        if (eta2_zero){ # don't use eta2
          sigma.d <- sigma/(ds.t^eta)
          mu.d    <- mu*sigma.d
        }else{ # use eta2 (and no tau)
          sigma.d <-    sigma/(ds.t^eta2)
          mu.d    <- mu*sigma/(ds.t^eta)
        }
      }else{ # use tau
        if (eta2_zero){ # don't use eta2 
          sigma.d <- sigma/(ds.t^eta)+tau
          mu.d    <- mu*sigma.d
        }else{ # use eta2 (and tau)
          sigma.d <-     sigma/(ds.t^eta2)+tau
          mu.d    <- mu*(sigma/(ds.t^eta)+tau)
        }
      }
248
      #sigma.d <- sigma/(ds.t^eta) # old
249
      
250
251
252
      y = (xdat - mu.d) / sigma.d # new
      #y = (xdat - mu*sigma.d) / sigma.d # derivation
      #y <- xdat/sigma.d - mu # old
253
      
254
      y <- 1 + xi * y
Laura Mack's avatar
Laura Mack committed
255

256

Felix Fauer's avatar
Felix Fauer committed
257
258
259
      if(!theta_zero) {if(any(theta < 0)) {return(10^6)} } # check definition condition for theta
      if(any(eta <= 0) || any(sigma.d <= 0) || any(y <= 0)) return(10^6)
      if(!tau_zero)   {if(any(tau < 0))    {return(10^6)} } # check definition condition for tau.
260
      if(!eta2_zero) {if(any(eta2 < 0))    {return(10^6)} } # check definition condition for eta2.
Felix Fauer's avatar
Felix Fauer committed
261
      
262
      sum(log(sigma.d)) +                 sum(y^(-1/xi)) +              sum(log(y) * (1/xi + 1))
263
    }
Laura Mack's avatar
Laura Mack committed
264
265


266
267
268
    # finding minimum of log-likelihood:
    x <- optim(init, gev.lik, hessian = TRUE, method = method,
               control = list(maxit = maxit, ...))
Laura Mack's avatar
Laura Mack committed
269

270
271
    # saving output parameters:
    z$conv <- x$convergence
272
273
274
    mut <- mutlink$linkinv(mumat %*% (x$par[1:npmu]))
    sc0 <- sigma0link$linkinv(sigmat %*% (x$par[seq(npmu + 1, length = npsc)]))
    xi <- xilink$linkinv(shmat %*% (x$par[seq(npmu + npsc + 1, length = npsh)]))
275
276
277
278
279
    if(!theta_zero){ #When user does NOT set theta parameter to zero (default)
      theta <- thetalink$linkinv(thmat %*% (x$par[seq(npmu + npsc + npsh + 1, length = npth)]))
    }else{ #When user requests theta_parameter to be zero
      theta <- thetalink$linkinv(thmat %*% (0))
    }
280
    eta <- etalink$linkinv(etmat %*% (x$par[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
281
282
283
284
285
286
287
    
    if(!eta2_zero){  #When user wants to use eta2 parameter
      eta2 <- eta2link$linkinv(e2mat %*% (x$par[seq(npmu + npsc + npsh + npth + npet + 1,length = npe2)]))
    }else{ #When user requests not to have eta2 parameter (default)
      eta2 <- eta
    }
    
288
    if(!tau_zero){  #When user does NOT set tau parameter to zero (not default)
289
      tau <- taulink$linkinv(tamat %*% (x$par[seq(npmu + npsc + npsh + npth + npet + npe2 + 1,length = npta)]))
Felix Fauer's avatar
Felix Fauer committed
290
291
292
    }else{ #When user requests tau parameter to be zero
      tau <- taulink$linkinv(tamat %*% (0))
    }
293
      
294
    z$nllh <- x$value
295
    # normalize data to standard Gumbel:
296
297
298
    #sc.d  <-      sc0/((ds+theta)^eta)+tau  # old
    sc.d  <-      sc0/((ds+theta)^eta2)+tau  # new
    mut.d <- mut*(sc0/((ds+theta)^eta )+tau) # new
Felix Fauer's avatar
Felix Fauer committed
299

300
    #z$data <-  - log(as.vector((1 + xi * (xdat/sc.d-mut))^(-1/xi)))  # old
301
    z$data <-  - log(as.vector((1 + xi * ((xdat-mut.d)/sc.d))^(-1/xi))) # new
302
    z$mle <- x$par
Laura Mack's avatar
Laura Mack committed
303
    test <- try(              # catch error
304
    z$cov <- solve(x$hessian) # invert hessian to get estimation on var-covar-matrix
305
306
307
308
309
    ,silent = TRUE )
    if("try-error" %in% class(test)){
      warning("Hessian could not be inverted. NAs were produced.")
      z$cov <- matrix(NA,length(z$mle),length(z$mle))
        }
Laura Mack's avatar
Laura Mack committed
310
    z$se <- sqrt(diag(z$cov)) # sqrt(digonal entries) = standart error of mle's
Felix Fauer's avatar
Felix Fauer committed
311
    'if (!theta_zero) {#When theta parameter is returned (default)
312
313
314
315
316
      if (!tau_zero){ # when tau is returned
        z$vals <- cbind(mut, sc0, xi, theta, eta, tau)
      }else{ # when tau is not returned
        z$vals <- cbind(mut, sc0, xi, theta, eta)
      }
317
    } else {#When theta parameter is not returned, asked by user
318
319
320
321
322
      if (!tau_zero){ # if tau is returned
        z$vals <- cbind(mut, sc0, xi, eta, tau)
      }else{ # if tau is not returned
        z$vals <- cbind(mut, sc0, xi, eta)
      }
Felix Fauer's avatar
Felix Fauer committed
323
    }'
324
    z$vals <- cbind(mut, sc0, xi, theta, eta, eta2, tau)
325
    z$init.vals <- unlist(init.vals)
326

Felix Fauer's avatar
Felix Fauer committed
327
328
329
330
331
    'names2 = c("mut","sigma0","xi")               # fixed part for set of names
    if(!theta_zero){names2=c(names2,"theta")}     # add theta (in case)
    names2 = c(names2, "eta")                     # add eta (always)
    if(!tau_zero){names2=c(names2, "tau")}        # add tau (in case)
    colnames(z$vals) <- names2'
332
    colnames(z$vals) <- c("mut","sigma0","xi","theta","eta","eta2","tau")
333
    
334
    z$ds <- ds
Laura Mack's avatar
Laura Mack committed
335
    z$theta_zero <- theta_zero #Indicates if theta parameter was set to zero by user.
336
337
    z$tau_zero <- tau_zero     #Indicates if tau parameter was set to zero by user. (default)
    z$eta2_zero <- eta2_zero   #Indicates if eta2 parameter was set to zero by user. (default)
338
    if(show) {
Laura Mack's avatar
Laura Mack committed
339
      if(z$trans) { # for nonstationary fit
Laura Mack's avatar
Laura Mack committed
340
        print(z[c(2, 4)]) # print model, link (3) , conv
Laura Mack's avatar
Laura Mack committed
341
        # print names of link functions:
342
        cat('$link\n')
Felix Fauer's avatar
Felix Fauer committed
343
        #if(!tau_zero){
344
        print(c(z$link$mutlink$name,z$link$sigma0link$name,z$link$xilink$name,z$link$thetalink$name,z$link$etalink$name,z$link$eta2link$name,z$link$taulink$name))
Felix Fauer's avatar
Felix Fauer committed
345
346
347
        #} else{
        #  print(c(z$link$mutlink$name,z$link$sigma0link$name,z$link$xilink$name,z$link$thetalink$name,z$link$etalink$name))
        #}
348
        cat('\n')
Laura Mack's avatar
Laura Mack committed
349
      }else{print(z[4])} # for stationary fit print only conv
Laura Mack's avatar
Laura Mack committed
350
      if(!z$conv){ # if fit converged
351
        print(z[c(5, 7, 9)]) # print nll, mle, se
Laura Mack's avatar
Laura Mack committed
352
      }
353
    }
354
    class(z) <- "gev.d.fit"
355
356
357
358
359
360
361
    invisible(z)
}


#### gev.d.init ####

# function to get initial values for gev.d.fit:
Laura Mack's avatar
Laura Mack committed
362
# obtain initial values
363
# by fitting every duration separately
364
365

# possible ways to improve:
366
# take given initial values into account, if there are any
367
368
# xi -> mean vs. median ... how do we improve that?
# mu_tilde -> is not very good for small sample sizes yet
Jana Ulrich's avatar
Jana Ulrich committed
369
# improved initial value for eta, by fitting both mu~d and sigma~d in log-log scale
370
371

#' @title get initial values for gev.d.fit
Jana Ulrich's avatar
Jana Ulrich committed
372
373
#' @description obtain initial values by fitting every duration separately
#' @param xdat vector of maxima for different durations
374
#' @param ds vector of durations belonging to maxima in xdat
375
#' @param link list of 5, link functions for parameters, created with \code{\link{make.link}}
376
#' @return list of initial values for mu_tilde, sigma_0, xi, theta, eta, eta2, tau
Laura Mack's avatar
Laura Mack committed
377
378
#' @importFrom stats lm
#' @importFrom stats median
379
#' @importFrom ismev gev.fit
Laura Mack's avatar
Laura Mack committed
380
#' @keywords internal
381

382
gev.d.init <- function(xdat,ds,link){
383
384
385
  durs <- unique(ds)
  mles <- matrix(NA, nrow=length(durs), ncol= 3)
  for(i in 1:length(durs)){
386
    test <- try(fit <- ismev::gev.fit(xdat[ds==durs[i]],show = FALSE),silent = TRUE)
387
    if("try-error" %in% class(test) | fit$conv!=0){mles[i,] <- rep(NA,3)}else{mles[i,] <- fit$mle}
388
  }
389
  if(all(is.na(mles))){stop('Initial values could not be computed for this dataset.')}
390
  # get values for sig0 and eta (also mu_0) from linear model in log-log scale
391
392
  lmsig <- lm(log(mles[,2])~log(durs))
  lmmu <- lm(log(mles[,1])~log(durs))
Laura Mack's avatar
Laura Mack committed
393

394
  # sig0 <- exp Intercept
395
  siginit <- link$sigma0link$linkfun(exp(lmsig$coefficients[[1]]))
Laura Mack's avatar
Laura Mack committed
396
  # eta <- mean of negativ slopes
397
  etainit <- link$etalink$linkfun(mean(c(-lmsig$coefficients[[2]],-lmmu$coefficients[[2]])))
398
  eta2init <- etainit + 0.0
Laura Mack's avatar
Laura Mack committed
399
  # mean of mu_d/sig_d
400
  # could try:
401
  # mu0/sig0 = exp(lmmu$coefficients[[1]])/exp(lmsig$coefficients[[1]])
402
  muinit <- link$mutlink$linkfun(median(c(mles[,1]/mles[,2]),na.rm = TRUE))
Laura Mack's avatar
Laura Mack committed
403
  # mean of shape parameters
404
  shinit <- link$xilink$linkfun(median(mles[,3],na.rm = TRUE))
405
  thetainit <- link$thetalink$linkfun(0)
Felix Fauer's avatar
Felix Fauer committed
406
  tauinit <- link$taulink$linkfun(0)
Laura Mack's avatar
Laura Mack committed
407

408
  return(list(mu=muinit,sigma=siginit,xi=shinit,theta=thetainit,eta=etainit, eta2=eta2init, tau=tauinit))
409
410
}

411
412
413
#### gev.d.lik ####

#' d-GEV Likelihood
414
#'
415
#' Computes (log-) likelihood of d-GEV model
416
#' @param xdat numeric vector containing observations
417
#' @param ds numeric vector containing corresponding durations (1/60 corresponds to 1 minute, 1 corresponds to 1 hour)
418
#' @param mut,sigma0,xi,theta,eta,eta2,tau numeric vectors containing corresponding estimates for each of the parameters
419
#' @param log Logical; if TRUE, the log likelihood is returned.
420
#'
Laura Mack's avatar
Laura Mack committed
421
#' @return single value containing (log) likelihood
422
423
424
#' @export
#'
#' @examples
425
#' # compute log-likelihood of observation values not included in fit
426
427
#' train.set <- example[example$d!=2,]
#' test.set <- example[example$d==2,]
428
#' fit <- gev.d.fit(train.set$dat,train.set$d,mutl = c(1,2),sigma0l = 1
429
430
#'           ,ydat = as.matrix(train.set[c('cov1','cov2')]))
#' params <- gev.d.params(fit,ydat = as.matrix(test.set[c('cov1','cov2')]))
431
#' gev.d.lik(xdat = test.set$dat,ds = test.set$d,mut = params[,1],sigma0 = params[,2],xi = params[,3]
Felix Fauer's avatar
Felix Fauer committed
432
#'           ,theta = params[,4],eta = params[,5],log=TRUE)
433
434
gev.d.lik <- function(xdat,ds,mut,sigma0,xi,theta,eta,log=FALSE,tau=0,eta2=NULL) {
  if (is.null(eta2)){eta2=eta}
435
  if(any(xi==0)){stop('Function is not defined for shape parameter of zero.')}
436
  if(any(! c(length(ds),length(mut),length(sigma0),length(xi),length(theta),length(eta),length(eta2),length(tau)) %in%
437
         c(1,length(xdat)))){
438
    stop('Input vectors differ in length, but must have the same length.')
439
  }
Laura Mack's avatar
Laura Mack committed
440

441
  ds.t <- ds+theta
442
443
444
  sigma.d <-     sigma0/(ds.t^eta2)+tau
  mu.d    <- mut*(sigma0/(ds.t^eta)+tau)
  y = (xdat - mu.d) / sigma.d # new
445
  y <- 1 + xi * y
446
447
448
449
450
  
  #sigma.d <- sigma0/(ds.t^eta) + tau # old
  #y <- xdat/sigma.d - mut            # old
  #y <- 1 + xi * y                    # old
  
451
452
453
454
455
  if(log){
    return(sum(log(sigma.d) + y^(-1/xi) + log(y) * (1/xi + 1)))
  }else{
    return(prod(sigma.d * exp(y^(-1/xi)) * y ^ (1/xi + 1)))
  }
Laura Mack's avatar
Laura Mack committed
456

457
}
458
459
460
461
462

#### gev.d.diag ####

#' Diagnostic Plots for d-gev Models
#'
Laura Mack's avatar
Laura Mack committed
463
464
#' @description  Produces diagnostic plots for d-gev models using
#' the output of the function \code{\link{gev.d.fit}}. Values for different durations can be plotted in
465
466
467
#' different colors of with different symbols.
#' @param fit object returned by \code{\link{gev.d.fit}}
#' @param subset an optional vector specifying a subset of observations to be used in the plot
468
#' @param cols optional either one value or vector of same length as \code{unique(fit$ds)} to
Laura Mack's avatar
Laura Mack committed
469
#' specify the colors of plotting points.
470
#' The default uses the \code{rainbow} function.
471
#' @param pch optional either one value or vector of same length as \code{unique(fit$ds)} containing
472
473
#' integers or symbols to specify the plotting points.
#' @param which string containing 'both', 'pp' or 'qq' to specify, which plots should be produced.
474
#' @param mfrow vector specifying layout of plots. If both plots should be produced separately,
475
476
#' set to \code{c(1,1)}.
#' @param legend logical indicating if legends should be plotted
Jana Ulrich's avatar
Jana Ulrich committed
477
478
#' @param title character vector of length 2, giving the titles for the pp- and the qq-plot
#' @param emp.lab,mod.lab character string containing names for empirical and model axis
Laura Mack's avatar
Laura Mack committed
479
#' @param ... additional parameters passed on to the plotting function
480
481
482
483
484
485
486
#'
#' @export
#' @importFrom graphics plot abline par title
#' @importFrom grDevices rainbow
#'
#' @examples
#' data('example',package ='IDF')
Laura Mack's avatar
Laura Mack committed
487
#'
488
#' fit <- gev.d.fit(xdat=example$dat,ds = example$d,ydat=as.matrix(example[,c('cov1','cov2')])
489
#'                  ,mutl=c(1,2),sigma0l=1)
Laura Mack's avatar
Laura Mack committed
490
491
492
#' # diagnostic plots for complete data
#' gev.d.diag(fit,pch=1)
#' # diagnostic plots for subset of data (e.g. one station)
493
#' gev.d.diag(fit,subset = example$cov1==1,pch=1)
Jana Ulrich's avatar
Jana Ulrich committed
494
495
496
gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1,2),legend=TRUE,
                       title=c('Residual Probability Plot','Residual Quantile Plot'),
                       emp.lab='Empirical',mod.lab='Model',...){
497
  # check parameter:
Jana Ulrich's avatar
Jana Ulrich committed
498
  if(!is.element(which,c('both','pp','qq'))) stop("Parameter 'which'= ",which,
499
500
501
                                                 " but only 'both','pp' or 'qq' are allowed.")
  # subset data
  df <- data.frame(data=fit$data,ds=fit$ds)
502
  if(!is.null(subset)){
Laura Mack's avatar
Laura Mack committed
503
    if(dim(df)[1]!=length(subset)){stop("Length of 'subset' does not match length of data
504
505
506
                                        'xdat' used for fitting.")}
    df <- df[subset,]
  }
Jana Ulrich's avatar
Jana Ulrich committed
507
508
509
510
511
  # get single durations
  durs <- sort(unique(df$ds))
  # rescale durations to assign colors
  df$cval <- sapply(df$ds,function(d){which(durs==d)})

Laura Mack's avatar
Laura Mack committed
512
  # sort data
513
  df <- df[order(df$data),]
Laura Mack's avatar
Laura Mack committed
514

515
516
517
518
519
520
521
  # plotting position
  n <- length(df$data)
  px <- (1:n)/(n + 1)

  # create plots:
  if(which=='both') par(mfrow=mfrow) # 2 subplots
  # colors and symbols
Jana Ulrich's avatar
Jana Ulrich committed
522
  if(is.null(cols))cols <- rainbow(length(durs))
523
  if(is.null(pch))pch <- df$cval
Laura Mack's avatar
Laura Mack committed
524

525
526
527
  if(which=='both'|which=='pp'){
    # pp
    plot(px, exp( - exp( - df$data)), xlab =
Jana Ulrich's avatar
Jana Ulrich committed
528
529
530
           emp.lab, ylab = mod.lab,col=cols[df$cval],pch=pch,...)
    abline(0, 1, col = 1,lwd=1)
    title(title[1])
Jana Ulrich's avatar
Jana Ulrich committed
531
    if(legend){legend('bottomright',legend = round(durs,digits = 2),pch=pch,
532
                      col = cols[1:length(durs)],title = 'Duration [h]',ncol = 2)}
533
534
535
536
  }
  if(which=='both'|which=='qq'){
    # qq
    plot( - log( - log(px)), df$data, ylab =
Jana Ulrich's avatar
Jana Ulrich committed
537
538
539
            emp.lab, xlab = mod.lab,col=cols[df$cval],pch=pch,...)
    abline(0, 1, col = 1,lwd=1)
    title(title[2])
Jana Ulrich's avatar
Jana Ulrich committed
540
    if(legend){legend('bottomright',legend = round(durs,digits = 2),pch=pch,
541
                      col = cols[1:length(durs)],title = 'Duration [h]',ncol = 2)}
542
543
544
545
546
547
548
549
  }
  if(which=='both') par(mfrow=c(1,1)) # reset par
}

#### gev.d.params ####

#' Calculate gev(d) parameters from \code{gev.d.fit} output
#'
Laura Mack's avatar
Laura Mack committed
550
551
#' @description function to calculate mut, sigma0, xi, theta, eta
#' (modified location, scale offset, shape, duration offset, duration exponent)
Jana Ulrich's avatar
Jana Ulrich committed
552
#' from results of \code{\link{gev.d.fit}} with covariates or link functions other than identity.
553
#' @param fit fit object returned by \code{\link{gev.d.fit}} or \code{\link{gev.fit}}
554
#' @param ydat A matrix containing the covariates in the same order as used in \code{gev.d.fit}.
555
#' @seealso \code{\link{IDF-package}}
556
#' @return data.frame containing mu_tilde, sigma0, xi, theta, eta, eta2, tau (or mu, sigma, xi for gev.fit objects)
557
#' @export
Laura Mack's avatar
Laura Mack committed
558
#'
559
560
561
#' @examples
#' data('example',package = 'IDF')
#' fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")])
562
#'                  ,mutl = c(1,2),sigma0l = 1)
563
564
565
#' gev.d.params(fit = fit,ydat = cbind(c(0.9,1),c(0.5,1)))


566
gev.d.params <- function(fit,ydat=NULL){
567
  if(!class(fit)%in%c("gev.d.fit","gev.fit"))stop("'fit' must be an object returned by 'gev.d.fit' or 'gev.fit'.")
568
  if(!is.null(ydat)){
Laura Mack's avatar
Laura Mack committed
569
    # check covariates matrix
570
571
572
    if(!is.matrix(ydat))stop("'ydat' must be of class matrix.")
    n.par <- max(sapply(fit$model,function(x){return(ifelse(is.null(x),0,max(x)))}))
    if(n.par>ncol(ydat))stop("Covariates-Matrix 'ydat' has ",ncol(ydat), " columns, but ", n.par," are required.")
573
574
575
576
  }else{if(!fit$trans){# no model -> no covariates matrix
    ydat <- matrix(1)
    }else{stop("To calculate parameter estimates, covariates matrix 'ydat' must be provided.")}
  }
Laura Mack's avatar
Laura Mack committed
577

578
579
580
581
  # number of parameters
  npmu <- length(fit$model[[1]]) + 1
  npsc <- length(fit$model[[2]]) + 1
  npsh <- length(fit$model[[3]]) + 1
582
  if(class(fit)=="gev.d.fit"){
583
584
585
586
587
    if(!fit$theta_zero){
      npth <- length(fit$model[[4]]) + 1 #Including theta parameter (default)]
    }else{
      npth <- 0
    }#With no theta parameter, asked by user
588
    npet <- length(fit$model[[5]]) + 1
589
590
591
592
593
    if(!fit$eta2_zero){
      npe2 <- length(fit$model[[6]]) + 1 #Including eta2 parameter (not default)]
    }else{
      npe2 <- 0
    }
594
    if(!fit$tau_zero){
595
      npta <- length(fit$model[[7]]) + 1 #Including tau parameter (not default)]
596
597
598
    }else{
      npta <- 0
    }#With no tau parameter, asked by user
599
  }
Laura Mack's avatar
Laura Mack committed
600

601
  # inverse link functions
602
  if(class(fit)=="gev.d.fit"){
603
604
605
    mulink <- fit$link$mutlink$linkinv
    siglink <- fit$link$sigma0link$linkinv
    shlink <- fit$link$xilink$linkinv
606
    if(!fit$theta_zero) thetalink <- fit$link$thetalink$linkinv
607
    etalink <- fit$link$etalink$linkinv
608
609
    if(!fit$eta2_zero) eta2link <- fit$link$eta2link$linkinv
    if(!fit$tau_zero) taulink <- fit$link$taulink$linkinv
610
611
612
613
614
  }else{
    mulink <- eval(parse(text=fit$link))[[1]]
    siglink <- eval(parse(text=fit$link))[[2]]
    shlink <- eval(parse(text=fit$link))[[3]]
  }
Laura Mack's avatar
Laura Mack committed
615

616
617
618
619
  # covariates matrices
  mumat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[1]]],dim(ydat)[1],npmu-1))
  sigmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[2]]],dim(ydat)[1],npsc-1))
  shmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[3]]],dim(ydat)[1],npsh-1))
620
621
622
  if(class(fit)=="gev.d.fit"){
    if(!fit$theta_zero){thmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[4]]],dim(ydat)[1],npth-1))}
    etmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[5]]],dim(ydat)[1],npet-1))
623
624
    if(!fit$eta2_zero) {e2mat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[6]]],dim(ydat)[1],npe2-1))}
    if(!fit$tau_zero)  {tamat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[7]]],dim(ydat)[1],npta-1))}
625
  }
Laura Mack's avatar
Laura Mack committed
626

627
628
629
630
  # calculate parameters
  mut <- mulink(mumat %*% (fit$mle[1:npmu]))
  sc0 <- siglink(sigmat %*% (fit$mle[seq(npmu + 1, length = npsc)]))
  xi <- shlink(shmat %*% (fit$mle[seq(npmu + npsc + 1, length = npsh)]))
631
  if(class(fit)=="gev.d.fit"){
632
633
634
635
636
637
    if(!fit$theta_zero){
      theta <- thetalink(thmat %*% (fit$mle[seq(npmu + npsc + npsh + 1, length = npth)]))
    }else{
      theta <- rep(0,dim(ydat)[1])
    }
    eta <- etalink(etmat %*% (fit$mle[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
638
639
640
641
642
    if(!fit$eta2_zero){
      eta2 <- eta2link(e2mat %*% (fit$mle[seq(npmu + npsc + npsh + npth + npet + 1, length = npe2)]))
    }else{
      eta2 <- eta #rep(0,dim(ydat)[1])
    }
643
    if(!fit$tau_zero){
644
      tau <- taulink(tamat %*% (fit$mle[seq(npmu + npsc + npsh + npth + npet + npe2 + 1, length = npta)]))
645
646
647
    }else{
      tau <- rep(0,dim(ydat)[1])
    }
648
    return(data.frame(mut=mut,sigma0=sc0,xi=xi,theta=theta,eta=eta, eta2=eta2, tau=tau))
649
650
651
  }else{
    return(data.frame(mu=mut,sig=sc0,xi=xi))
  }
652
653
654
655
}


#### example data ####
656
#' Sampled data for duration-dependent GEV
657
#' 
Laura Mack's avatar
Laura Mack committed
658
#' @description
659
#' Randomly sampled data set used for running the example code, containing:
660
661
662
663
#' \itemize{
#'   \item \code{$xdat}: 'annual' maxima values
#'   \item \code{$ds}: corresponding durations
#'   \item \code{$cov1}, \code{$cov2}: covariates}
664
#' d-GEV parameters used for sampling:
665
#' \itemize{
666
667
668
669
#'   \item \eqn{\tilde{\mu} = 4 + 0.2 cov_1 +0.5 cov_2}
#'   \item \eqn{\sigma_0 = 2+0.5 cov_1}
#'   \item \eqn{\xi = 0.5}
#'   \item \eqn{\theta = 0}
670
#'   \item \eqn{\eta = 0.5}
671
#'   \item \eqn{\eta_2 = 0.5}
672
#'   \item \eqn{\tau = 0}}
Laura Mack's avatar
Laura Mack committed
673
#'
674
675
676
677
678
#'
#' @docType data
#' @keywords datasets
#' @name example
#' @usage data('example',package ='IDF')
679
#' @format A data frame with 330 rows and 4 variables
680
NULL