Package 'SpatialExtremes'

Anova Tables

Description

Computes analysis of deviance for objects of class 'maxstab' or 'spatgev'.

Usage

## S3 method for class 'maxstab'
anova(object, object2, method = "RJ", square = "chol", ...)
## S3 method for class 'spatgev'
anova(object, object2, method = "RJ", square = "chol", ...)

Arguments

object, object2	Two objects of class 'maxstab' or 'spatgev'.
method	Character string. Must be one of "CB" or "RJ" for the Chandler and Bate or the Rotnitzky and Jewell approaches respectively. See function profile.
square	The choice for the matrix square root. This is only useful for the 'CB' method. Must be one of 'chol' (Cholesky) or 'svd' (Singular Value Decomposition).
...	Other options to be passed to the anova function.

Details

As ”maxstab” objects are fitted using pairwise likelihood, the model is misspecified. As a consequence, the likelihood ratio statistic is no longer $\chi^2$ distributed. To compute the anova table, we use the methodology proposed by Rotnitzky and Jewell to adjust the distribution of the likelihood ratio statistic.

Value

This function returns an object of class anova. These objects represent analysis-of-deviance tables.

Author(s)

Mathieu Ribatet

References

Chandler, R. E. and Bate, S. (2007) Inference for clustered data using the independence loglikelihood Biometrika, 94, 167–183.

Rotnitzky, A. and Jewell, N. (1990) Hypothesis testing of regression parameters in semiparametric generalized linear models for cluster correlated data. Biometrika 77, 485–497.

See Also

fitmaxstab, fitspatgev, profile, TIC

Examples

##Define the coordinates of each location
n.site <- 30
locations <- matrix(rnorm(2*n.site, sd = sqrt(.2)), ncol = 2)
colnames(locations) <- c("lon", "lat")

##Simulate a max-stable process - with unit Frechet margins
data <- rmaxstab(50, locations, cov.mod = "gauss", cov11 = 100, cov12 =
25, cov22 = 220)

##Now define the spatial model for the GEV parameters
param.loc <- -10 + 2 * locations[,2]
param.scale <- 5 + 2 * locations[,1] + locations[,2]^2
param.shape <- rep(0.2, n.site)

##Transform the unit Frechet margins to GEV 
for (i in 1:n.site)
  data[,i] <- frech2gev(data[,i], param.loc[i], param.scale[i],
param.shape[i])

##Define three models for the GEV margins to be fitted
loc.form <- loc ~ lat
scale.form <- scale ~ lon + I(lat^2)
shape.form <- shape ~ lon

M0 <- fitspatgev(data, locations, loc.form, scale.form, shape.form)
M1 <- fitspatgev(data, locations, loc.form, scale.form, shape.form,
shapeCoeff2 = 0)

##Model selection
anova(M0, M1)
anova(M0, M1, method = "CB", square = "svd")

---

Pairwise empirical and extremal concurrence probabilities

Description

This function computes the pairwise empirical or the pairwise extremal concurrence probability estimates.

Usage

concprob(data, coord, fitted, n.bins, add = FALSE, xlim = c(0,
max(dist)), ylim = c(min(0, concProb), max(1, concProb)), col = 1:2,
which = "kendall", xlab, ylab, block.size = floor(nrow(data)^(1/3)),
plot = TRUE, compute.std.err = FALSE, ...)

Arguments

data	A matrix representing the data. Each column corresponds to one location.
coord	A matrix that gives the coordinates of each location. Each row corresponds to one location.
fitted	An object of class maxstab - usually the output of the fitmaxstab function. May be missing.
n.bins	The number of bins to be used. If missing, pairwise F-madogram estimates will be computed.
xlim, ylim	A numeric vector of length 2 specifying the x/y coordinate ranges.
col	The colors used for the points and optionnaly the fitted curve.
which	A character string specifying which estimator should be used. Should be one of "emp" (empirical), "boot" (bootstrap version) and "kendall" (kendall based).
xlab, ylab	The labels for the x/y-axis (may be missing).
add	Logical. If TRUE, the plot is added to the current figure; otherwhise (default) a new plot is computed.
block.size	Integer specifying the block size for the empirical and bootstrap estimator.
plot	Logical. If TRUE (default) a plot is produced.
compute.std.err	Logical. If TRUE, standard errors are estimated using a jackknife procedure. It is currently only available for the Kendall estimator.
...	Additional options to be passed to the plot function.

Value

This function returns invisibly a matrix containing the pairwise distances and the concurrence probability estimates.

Author(s)

Mathieu Ribatet

References

Dombry, C., Ribatet, M. and Stoev, S. (2017) Probabilities of concurrent extremes. To appear in JASA

See Also

fmadogram, lmadogram

Examples

n.site <- 25
locations <- matrix(runif(2*n.site, 0, 10), ncol = 2)
colnames(locations) <- c("lon", "lat")

##Simulate a max-stable process - with unit Frechet margins
n.obs <- 100
data <- rmaxstab(n.obs, locations, cov.mod = "whitmat", nugget = 0, range = 1,
smooth = 1.75)

##Compute the F-madogram
concprob(data, locations)

##Compare the F-madogram with a fitted max-stable process
fitted <- fitmaxstab(data, locations, "whitmat", nugget = 0)
concprob(fitted = fitted)

---

Maps of concurrence probabilities/expected concurrence cell area

Description

This function produces maps for concurrence probabilities or expected concurrence cell areas.

Usage

concurrencemap(data, coord, which = "kendall", type = "cell", n.grid =
100, col = cm.colors(64), plot = TRUE, plot.border = NULL,
compute.std.err = FALSE, ...)

Arguments

data	A matrix representing the data. Each column corresponds to one location.
coord	A matrix that gives the coordinates of each location. Each row corresponds to one location.
which	A character string specifying which estimator should be used. Should be one of "emp" (empirical), "boot" (bootstrap version) and "kendall" (kendall based).
type	Either "cell" for cell areas or a integer between 1 and the number of locations specifying which site should be used as reference location—see Details.
n.grid	Integer specifying the size of the prediction grid.
col	The colors used to produce the map.
plot	Logical. If TRUE (default), a map is produced; otherwise results are invisibly returned.
plot.border	The name of an R function that can be used to plot the border of the study region. If NULL, no border are plotted.
compute.std.err	Logical. If TRUE, a map of standard errors is also produced. It is currently only available for concurrence probability maps.
...	Additional options to be passed to the image function.

Value

This function returns invisibly a list with the x and y coordinates and the corresponding values for the estimated concurrence probabilities or expected concurrence cell area.

Author(s)

Mathieu Ribatet

References

Dombry, C., Ribatet, M. and Stoev, S. (2015) Probabilities of concurrent extremes. Submitted

See Also

concprob

Examples

##require(maps) ## <<-- to plot US borders
data(USHCNTemp)
coord <- as.matrix(metadata[,2:3])

## Subset the station to have a fast example
n.site <- 30
chosen.site <- sample(nrow(coord), n.site)
coord <- coord[chosen.site,]
maxima.summer <- maxima.summer[,chosen.site]

## Define a function to plot the border
border <- function(add = FALSE) maps::map("usa", add = add)

par(mar = rep(0, 4))
## Produce a pairwise concurrence probability map w.r.t. station number 15
concurrencemap(maxima.summer, coord, type = 15, plot.border = border,
 compute.std.err = TRUE)

## Produce the expected concurrence cell area
concurrencemap(maxima.summer, coord, plot.border = border)

---

Produces a conditional 2D map from a fitted max-stable process

Description

Produces a conditional 2D map from a fitted max-stable process.

Usage

condmap(fitted, fix.coord, x, y, covariates = NULL, ret.per1 = 100,
ret.per2 = ret.per1, col = terrain.colors(64), plot.contour = TRUE,
...)

Arguments

fitted	An object of class maxstab. Most often, it will be the output of the function fitmaxstab.
fix.coord	The spatial coordinates of the location from which the conditional quantile is computed.
x, y	Numeric vector defining the grid at which the levels are computed.
covariates	An array specifying the covariates at each grid point defined by x and y. If NULL, no covariate is needed. See map to see how to build it.
ret.per1, ret.per2	Numerics giving the return period for which the quantile map is plotted. See details.
col	A list of colors such as that generated by 'rainbow', 'heat.colors', 'topo.colors', 'terrain.colors' or similar functions.
plot.contour	Logical. If TRUE (default), contour lines are added to the plot.
...	Several arguments to be passed to the image function.

Details

The function solves the following equation:

$\Pr\left[Z(x_2) > z_2 | Z(x_1) > z_1 \right] = \frac{1}{T_2}$

where $z_1 = -1 / \log(1 - 1/T_1)$.

In other words, it computes, given that at location $x_1$ we exceed the level $z_1$, the levels which is expected to be exceeded in average every $T_2$ year.

Value

A plot. Additionally, a list with the details for plotting the map is returned invisibly.

Author(s)

Mathieu Ribatet

See Also

map, filled.contour, heatmap, heat.colors, topo.colors, terrain.colors, rainbow

Examples

##Define the coordinate of each location
n.site <- 30
locations <- matrix(runif(2*n.site, 0, 10), ncol = 2)
colnames(locations) <- c("lon", "lat")

##Simulate a max-stable process - with unit Frechet margins
data <- rmaxstab(50, locations, cov.mod = "whitmat", nugget = 0, range =
2, smooth = 1)

##Now define the spatial model for the GEV parameters
param.loc <- -10 - 4 * locations[,1] + locations[,2]^2
param.scale <- 5 + locations[,2] + locations[,1]^2 / 10
param.shape <- rep(.2, n.site)

##Transform the unit Frechet margins to GEV
for (i in 1:n.site)
  data[,i] <- frech2gev(data[,i], param.loc[i], param.scale[i],
param.shape[i])

##Define a model for the GEV margins to be fitted
##shape ~ 1 stands for the GEV shape parameter is constant
##over the region
loc.form <- loc ~ lon + I(lat^2)
scale.form <- scale ~ lat + I(lon^2)
shape.form <- shape ~ 1

##  1- Fit a max-stable process
fitted <- fitmaxstab(data, locations, "whitmat", loc.form, scale.form,
                     shape.form, nugget = 0)

cond.coord <- c(5.1, 5.1)
condmap(fitted, cond.coord, seq(0, 10, length = 25), seq(0,10, length
 =25), ret.per1 = 100, ret.per2 = 1.5)
points(t(cond.coord), pch = "*", col = 2, cex = 2)

---

Conditional simulation of Gaussian random fields

Description

This function generates conditional simulation of Gaussian random fields from the simple kriging predictor.

Usage

condrgp(n, coord, data.coord, data, cov.mod = "powexp", mean = 0, sill =
1, range = 1, smooth = 1, grid = FALSE, control = list())

Arguments

n	Integer. The number of conditional simulations.
coord	A numeric vector or matrix specifying the coordinates where the process has to be generated. If coord is a matrix, each row specifies one location.
data.coord	A numeric vector or matrix specifying the coordinates where the process is conditioned.
data	A numeric vector giving the conditioning observations.
cov.mod	A character string specifying the covariance function family. Must be one of "whitmat", "powexp", "cauchy" or "bessel" for the Whittle-Mater, the powered exponential, the Cauchy or Bessel covariance families.
mean, sill, range, smooth	The mean, sill, range and smooth of the Gaussian process.
grid	Logical. Does coord specifies a grid?
control	A named list passing options to the simulation method of Gaussian processes — see rgp.

Value

A list with components:

coord	The coordinates at which the process was simulated;
cond.sim	The simulated process;
data.coord	The coordinates of the conditioning locations;
data	The conditioning observations;
cov.mod	The covariance function family;
grid	Does coord specifies a grid?

Author(s)

Mathieu Ribatet

See Also

kriging, rgp.

Examples

## Several conditional simulations
n.site <- 50
n.sim <- 512

x.obs <- runif(n.site, -100, 100)
x.sim <- seq(-100, 100, length = n.sim)

data <- rgp(1, x.obs, "whitmat", sill = 1, range = 10, smooth = 0.75)

sim <- condrgp(5, x.sim, x.obs, data, "whitmat", sill = 1, range =
10, smooth = 0.75)

matplot(x.sim, t(sim$cond.sim),  type = "l", lty = 1, xlab = "x", ylab =
expression(Y[cond](x)))
points(x.obs, data, pch = 21, bg = 1)
title("Five conditional simulations")

## Comparison between one conditional simulations and the kriging
## predictor on a grid
x.obs <- matrix(runif(2 * n.site, -100, 100), ncol = 2)
x <- y <- seq(-100, 100, length = 100)
x.sim <- cbind(x, y)

data <- rgp(1, x.obs, "whitmat", sill = 1, range = 50, smooth = 0.75)

krig <- kriging(data, x.obs, x.sim, "whitmat", sill = 1, range = 50,
smooth = 0.75, grid = TRUE)
sim <- condrgp(1, x.sim, x.obs, data, "whitmat", sill = 1, range = 50,
smooth = 0.75, grid = TRUE)

z.lim <- range(c(sim$cond.sim, data, krig$krig.est))
breaks <- seq(z.lim[1], z.lim[2], length = 65)
col <- heat.colors(64)
idx <- as.numeric(cut(data, breaks))

op <- par(mfrow = c(1,2))
image(x, y, krig$krig.est, col = col, breaks = breaks)
points(x.obs, bg = col[idx], pch = 21)
title("Kriging predictor")
image(x, y, sim$cond.sim, col = col, breaks = breaks)
points(x.obs, bg = col[idx], pch = 21)
title("Conditional simulation")
## Note how the background colors of the above points matches the ones
## returned by the image function
par(op)

---

Conditional simulation of max-linear random fields

Description

This function generates (approximate) conditional simulation of unit Frechet max-linear random fields. It can be used to get approximate conditional simulation for max-stable processes.

Usage

condrmaxlin(n, coord, data.coord, data, cov.mod = "gauss", ..., grid =
FALSE, p = 10000)

Arguments

n	Integer. The number of conditional simulations.
coord	A numeric vector or matrix specifying the coordinates where the process has to be generated. If coord is a matrix, each row specifies one locations.
data.coord	A numeric vector or matrix specifying the coordinates where the process is conditioned.
data	A numeric vector giving the conditioning observations.
cov.mod	A character string specifying the max-stable model. See section Details.
...	The parameters of the max-stable model. See section Details.
grid	Logical. Does coord defines a grid?
p	An integer. The number of unit Frechet random variables used in the max-linear approximation.

Details

Any unit Frechet max-stable processes $\{Z(x)\}$ can be approximated by a unit Frechet max-linear process, i.e.,

$Z(x) \approx \max_{j=1, \ldots, p} f_j(x) Z_j,$

where $f_j$ are non-negative deterministic functions, $p$ is a sufficiently large integer and $Z_j$ are independent unit Frechet random variables. Note that to ensure unit Frechet margins, the following condition has to be satisfied

$\sum_{j=1, \ldots, p} f_j(x) = 1,$

for all $x$.

Currently only the discretized Smith model is implemented for which $f_j(x) = c(p) \varphi(x - u_j ; \Sigma)$ where $\varphi(\cdot; \Sigma)$ is the zero mean (multivariate) normal density with covariance matrix $\Sigma$, $u_j$ is a sequence of deterministic points appropriately chosen and $c(p)$ is a constant ensuring unit Frechet margins.

Value

A matrix containing observations from the required max-stable model. Each column represents one stations. If grid = TRUE, the function returns an array of dimension nrow(coord) x nrow(coord) x n.

Warnings

It may happen that some conditional observations are not honored because the approximation of a max-stable process by a max-linear one isn't accurate enough! Sometimes taking a larger p solves the issue.

Author(s)

Mathieu Ribatet

References

Wang, Y. and Stoev, S. A. (2011) Conditional Sampling for Max-Stable Random Fields. Advances in Applied Probability.

See Also

rmaxstab, condrmaxlin

Examples

## One dimensional conditional simulations
n.cond.site <- 10
cond.coord <- runif(n.cond.site, -10, 10)
data <- rmaxlin(1, cond.coord, var = 3, p = 10000)

x <- seq(-10, 10, length = 250)
cond.sim <- condrmaxlin(5, x, cond.coord, data, var = 3)

matplot(x, t(log(cond.sim)), type = "l", lty = 1, pch = 1)
points(cond.coord, log(data))

## Two dimensional conditional simulation
cond.coord <- matrix(runif(2 * n.cond.site, -10, 10), ncol = 2)
data <- rmaxstab(1, cond.coord, "gauss", cov11 = 4, cov12 = 0, cov22 = 4)

x <- y <- seq(-10, 10, length = 75)
cond.sim <- condrmaxlin(4, cbind(x, y), cond.coord, data, cov11 = 4,
cov12 = 0, cov22 = 4, grid = TRUE, p = 2000)
## Note p is set to 2000 for CPU reasons but is likely to be too small

op <- par(mfrow = c(2, 2), mar = rep(1, 4))
for (i in 1:4){
image(x, y, log(cond.sim[,,i]), col = heat.colors(64), xaxt = "n", yaxt
= "n", bty = "n")
contour(x, y, log(cond.sim[,,i]), add = TRUE)
text(cond.coord[,1], cond.coord[,2], round(log(data), 2), col = 3)
}
par(op)

---

Conditional simulation of max-stable processes

Description

This function performs conditional simulation of various max-stable processes.

Usage

condrmaxstab(k = 1, coord, cond.coord, cond.data, cov.mod = "powexp",
..., do.sim = TRUE, thin = n.cond, burnin = 50, parts)

Arguments

k	An integer. The number of conditional simulations to be generated.
coord	A vector or matrix that gives the coordinates of each location. Each row corresponds to one location - if any.
cond.coord	A vector or matrix that gives the coordinates of each conditional location. Each row corresponds to one location - if any.
cond.data	A vector that gives the conditional values at the corresponding conditioning locations. Each row corresponds to one location - if any.
cov.mod	A character string that gives the max-stable model. This must be one of "brown" for the Brown-Resnick model, or "whitmat", "cauchy", "powexp" and "bessel" for the Schlather model with the given correlation family.
...	The parameters of the max-stable model. See rmaxstab for more details.
do.sim	A logical value. If TRUE (the default), the conditional simulations are performed; otherwise only the simulated random partitions, i.e., the hitting scenarios, are returned.
thin	A positive integer giving by which amount the generated Markov chain should be thinned. This is only useful when the number of conditioning locations is greater than 7.
burnin	A positive integer giving the duration of the burnin period of the Markov chain.
parts	A matrix giving the hitting scenarios. Each row corresponds to one hitting scenarios. If missing then a Gibbs sampler will be used to generate such hitting scenarios.

Details

The algorithm consists in three steps:

1. Draw a random partition $\theta$ from

   $\Pr\{\theta = \tau \mid Z(x) = z\}$

2. Given the random partition, draw the extremal functions from

   $\Pr\{\varphi^+ \in \cdot \mid Z(x) = z, \theta = \tau\}$

3. Independently, draw the sub-extremal functions, i.e.,

   $\max_{i \ge 1} \varphi_i 1_{\{\varphi_i(x) < z\}}$

The distribution in Step 1 is usually intractable and in such cases a random scan Gibbs sampler will be used to sample from this distribution.

Value

This function returns a list whose components are

sim	The conditional simulations. Beware the first values corresponds to the conditioning values.
sub.ext.fct	The values of the sub-extremal functions.
ext.fct	The values of the extremal functions.
timings	The timings in seconds for each step of the algorithm.

Warning

This function can be extremely time consuming when the number of conditioning locations is large.

Author(s)

Mathieu Ribatet

References

Dombry, C. and Eyi-Minko, F. (2012) Regular conditional distributions of max infinitely divisible processes. Submitted.

Dombry, C., Eyi-Minko, F. and Ribatet, M. (2012) Conditional simulation of max-stable processes. To appear in Biometrika.

See Also

rmaxstab, condrgp

Examples

n.sim <- 50
n.cond <- 5

range <- 10
smooth <- 1.5

n.site <- 200
coord <- seq(-5, 5, length = n.site)
cond.coord <- seq(-4, 4, length = n.cond)
all.coord <- c(cond.coord, coord)

all.cond.data <- rmaxstab(1, all.coord, "powexp", nugget = 0, range = range,
                      smooth = smooth)
cond.data <- all.cond.data[1:n.cond]

ans <- condrmaxstab(n.sim, coord, cond.coord, cond.data, range = range,
                    smooth = smooth, cov.mod = "powexp")

idx <- order(all.coord)
matplot(coord, t(log(ans$sim)), type = "l", col = "grey", lty = 1,
        xlab = expression(x), ylab = expression(Z(x)))
lines(all.coord[idx], log(all.cond.data)[idx])
points(cond.coord, log(cond.data), pch = 15, col = 2)

---

Defines and computes covariance functions

Description

This function defines and computes several covariance function either from a fitted “max-stable” model; either by specifying directly the covariance parameters.

Usage

covariance(fitted, nugget, sill, range, smooth, smooth2 = NULL, cov.mod =
"whitmat", plot = TRUE, dist, xlab, ylab, col = 1, ...)

Arguments

fitted	An object of class “maxstab”. Most often this will be the output of the fitmaxstab function. May be missing if sill, range, smooth and cov.mod are given.
nugget, sill, range, smooth, smooth2	The nugget, sill, scale and smooth parameters for the covariance function. May be missing if fitted is given.
cov.mod	Character string. The name of the covariance model. Must be one of "whitmat", "cauchy", "powexp", "bessel" or "caugen" for the Whittle-Matern, Cauchy, Powered Exponential, Bessel and Generalized Cauchy models. May be missing if fitted is given.
plot	Logical. If TRUE (default) the covariance function is plotted.
dist	A numeric vector corresponding to the distance at which the covariance function should be evaluated. May be missing.
xlab, ylab	The x-axis and y-axis labels. May be missing.
col	The color to be used for the plot.
...	Several option to be passed to the plot function.

Details

Currently, four covariance functions are defined: the Whittle-Matern, powered exponential (also known as "stable"), Cauchy and Bessel models. These covariance functions are defined as follows for $h > 0$

Whittle-Matern

$\gamma(h) = \sigma \frac{2^{1-\kappa}}{\Gamma(\kappa)} \left(\frac{h}{\lambda} \right)^{\kappa} K_{\kappa}\left(\frac{h}{\lambda} \right)$

Powered Exponential

$\gamma(h) = \sigma \exp \left[- \left(\frac{h}{\lambda} \right)^{\kappa} \right]$

Cauchy

$\gamma(h) = \sigma \left[1 + \left(\frac{h}{\lambda} \right)^2 \right]^{-\kappa}$

Bessel

$\gamma(h) = \sigma \left(\frac{2 \lambda}{h}\right)^{\kappa} \Gamma(\kappa + 1) J_{\kappa}\left(\frac{h}{\lambda} \right)$

Generalized Cauchy

$\gamma(h) = \sigma \left\{1 + \left(\frac{h}{\lambda} \right)^{\kappa_2} \right\}^{-\kappa / \kappa_2}$

where $\sigma$, $\lambda$ and $\kappa$ are the sill, the range and shape parameters, $\Gamma$ is the gamma function, $K_{\kappa}$ and $J_\kappa$ are both modified Bessel functions of order $\kappa$. In addition a nugget effect can be set that is there is a jump at the origin, i.e., $\gamma(o) = \nu + \sigma$, where $\nu$ is the nugget effect.

Value

This function returns the covariance function. Eventually, if dist is given, the covariance function is computed for each distance given by dist. If plot = TRUE, the covariance function is plotted.

Author(s)

Mathieu Ribatet

Examples

## 1- Calling covariance using fixed covariance parameters
covariance(nugget = 0, sill = 1, range = 1, smooth = 0.5, cov.mod = "whitmat")
covariance(nugget = 0, sill = 0.5, range = 1, smooth = 0.5, cov.mod = "whitmat",
  dist = seq(0,5, 0.2), plot = FALSE)

## 2- Calling covariance from a fitted model
##Define the coordinate of each location
n.site <- 30
locations <- matrix(runif(2*n.site, 0, 10), ncol = 2)
colnames(locations) <- c("lon", "lat")

##Simulate a max-stable process - with unit Frechet margins
data <- rmaxstab(30, locations, cov.mod = "whitmat", nugget = 0, range =
3, smooth = 1)

##Fit a max-stable model
fitted <- fitmaxstab(data, locations, "whitmat", nugget = 0)
covariance(fitted, ylim = c(0, 1))
covariance(nugget = 0, sill = 1, range = 3, smooth = 1, cov.mod = "whitmat", add =
TRUE, col = 3)
title("Whittle-Matern covariance function")
legend("topright", c("Theo.", "Fitted"), lty = 1, col = c(3,1), inset =
.05)

---

Estimates the penalty coefficient from the cross-validation criterion

Description

Estimates the penalty coefficient from the cross-validation criterion.

Usage

cv(y, x, knots, degree, plot = TRUE, n.points = 150, ...)

Arguments

y	The response vector.
x	A vector/matrix giving the values of the predictor variable(s). If x is a matrix, each row corresponds to one observation.
knots	A vector givint the coordinates of the knots.
degree	The degree of the penalized smoothing spline.
plot	Logical. If TRUE (default), the evolution of the CV score as the penalty increases is plotted.
n.points	A numeric giving the number of CV computations needed to produce the plot.
...	Options to be passed to the nlm function.

Details

For every linear smoother e.g. $\hat{y} = S_\lambda y$, the cross-validation criterion consists in minimizing the following quantity:

$CV(\lambda) = \sum_{i=1}^n \left(\frac{y_i - \hat{y}_i}{1 - S_{\lambda,ii}} \right)^2$

where $\lambda$ is the penalty coefficient, $n$ the number of observations and $S_{\lambda,ii}$ the i-th diagonal element of the matrix $S_\lambda$.

Value

A list with components 'penalty', 'cv' and 'nlm.code' which give the location of the minimum, the value of the cross-validation criterion at that point and the code returned by the nlm function - useful to assess for convergence issues.

Author(s)

Mathieu Ribatet

References

Ruppert, D. Wand, M.P. and Carrol, R.J. (2003) Semiparametric Regression Cambridge Series in Statistical and Probabilistic Mathematics.

See Also

cv

Examples

n <- 200
x <- runif(n)
fun <- function(x) sin(3 * pi * x)
y <- fun(x) + rnorm(n, 0, sqrt(0.4))
knots <- quantile(x, prob = 1:(n/4) / (n/4 + 1))
cv(y, x, knots = knots, degree = 3)

---

Deviance Information Criterion

Description

This function computes the Deviance Information Criterion (DIC), and related quantities, which is a hierarchical modeling generalization of the Akaike Information Criterion. It is useful for Bayesian model selection.

Usage

DIC(object, ..., fun = "mean")

Arguments

object	An object of class latent — typically this will be the output of latent.
...	Optional arguments. Not implemented.
fun	A chararcter string giving the name of the function to be used to summarize the Markov chain. The default is to consider the posterior mean.

Details

The deviance is

$D(\theta) = -2 \log \pi(y \mid \theta),$

where $y$ are the data, $\theta$ are the unknown parameters of the models and $\pi(y \mid \theta)$ is the likelihood function. Thus the expected deviance, a measure of how well the model fits the data, is given by

$\overline{D} = {\rm E}_{\theta}[D(\theta)],$

while the effective number of parameters is

$p_D = \overline{D} - D(\theta^*),$

where $\theta^*$ is point estimate of the posterior distribution, e.g., the posterior mean. Finally the DIC is given by

${\rm DIC} = p_D + \overline{D}.$

In accordance with the AIC, models with smaller DIC should be preferred to models with larger DIC. Roughly speaking, differences of more than 10 might rule out the model with the higher DIC, differences between 5 and 10 are substantial.

Value

A vector of length 3 that returns the DIC, effective number of parameters eNoP and an estimate of the expected deviance Dbar.

Author(s)

Mathieu Ribatet

References

Spiegelhalter, D. J., Best, N. G., Carlin, B. P. and Van Der Linde, A. (2002) Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society: Series B 64, 583–639.

See Also

AIC

Examples

## Generate realizations from the model
n.site <- 15
n.obs <- 35
coord <- cbind(lon = runif(n.site, -10, 10), lat = runif(n.site, -10 , 10))

gp.loc <- rgp(1, coord, "powexp", sill = 4, range = 20, smooth = 1)
gp.scale <- rgp(1, coord, "powexp", sill = 0.4, range = 5, smooth = 1)
gp.shape <- rgp(1, coord, "powexp", sill = 0.01, range = 10, smooth = 1)

locs <- 26 + 0.5 * coord[,"lon"] + gp.loc
scales <- 10 + 0.2 * coord[,"lat"] + gp.scale
shapes <- 0.15 + gp.shape

data <- matrix(NA, n.obs, n.site)
for (i in 1:n.site)
  data[,i] <- rgev(n.obs, locs[i], scales[i], shapes[i])

loc.form <- y ~ lon
scale.form <- y ~ lat
shape.form <- y ~ 1

hyper <- list()
hyper$sills <- list(loc = c(1,8), scale = c(1,1), shape = c(1,0.02))
hyper$ranges <- list(loc = c(2,20), scale = c(1,5), shape = c(1, 10))
hyper$smooths <- list(loc = c(1,1/3), scale = c(1,1/3), shape = c(1, 1/3))
hyper$betaMeans <- list(loc = rep(0, 2), scale = c(9, 0), shape = 0)
hyper$betaIcov <- list(loc = solve(diag(c(400, 100))),
                       scale = solve(diag(c(400, 100))),
                       shape = solve(diag(c(10), 1, 1)))

## We will use an exponential covariance function so the jump sizes for
## the shape parameter of the covariance function are null.
prop <- list(gev = c(2.5, 1.5, 0.3), ranges = c(40, 20, 20), smooths = c(0,0,0))
start <- list(sills = c(4, .36, 0.009), ranges = c(24, 17, 16), smooths
              = c(1, 1, 1),  beta = list(loc = c(26, 0), scale = c(10, 0),
                               shape = c(0.15)))

mc <- latent(data, coord, loc.form = loc.form, scale.form = scale.form,
             shape.form = shape.form, hyper = hyper, prop = prop, start = start,
             n = 500)
DIC(mc)

---

Computes distance between pairs of locations

Description

This function computes euclidean distance or distance vector for each pair of a set of spatial locations.

Usage

distance(coord, vec = FALSE)

Arguments

coord	A matrix representing the coordinates of each locations. Each row corresponds to one location.
vec	Logical. If FALSE (default), the euclidean distance is computed. Otherwise, “distance vectors” are returned.

Value

If vec = FALSE, this function returns a vector that gives the euclidean distance for each pair of locations. Otherwise, this is a matrix where each column correspond to one dimension - i.e. longitude, latitude, ...

Author(s)

Mathieu Ribatet

See Also

dist

Examples

coords <- cbind(seq(0, 10, by = 2), seq(10, 0, by = -2))
distance(coords)
distance(coords, vec = TRUE)

---

Plots the extremal coefficient

Description

Plots the extremal coefficient evolution as the coordinates evolves.

Usage

extcoeff(fitted, cov.mod, param, n = 200, xlab, ylab, ...)

Arguments

fitted	A object of class "maxstab". Most often, it will be the output of the function fitmaxstab. If missing, then cov.mod and param should be supplied.
cov.mod	A character string corresponding the the covariance model in the max-stable representation. Must be one of "gauss" for the Smith's model; or "whitmat", "cauchy" or "powexp" for the Whittle-Matern, the Cauchy and the Powered Exponential covariance family with the Schlather's model. May be missing if fitted is given.
param	Numeric vector of length 3. The parameters for the Smith's or Schlather model - i.e. c(cov11, cov12, cov22) or c(nugget, range, smooth). Please respect this order.
n	Numeric. n^2 corresponds to the total number of estimated extremal coefficients for the contour plot.
xlab, ylab	The x-axis and y-axis labels. May be missing.
...	Several options to be passed to the contour function.

Value

A plot.

Author(s)

Mathieu Ribatet

See Also

fitmaxstab

Examples

## 1- Random field generation
n.site <- 30
locations <- matrix(runif(2*n.site, 0, 10), ncol = 2)
colnames(locations) <- c("lon", "lat")

data <- rmaxstab(60, locations, cov.mod = "whitmat", nugget = 0, range =
3, smooth = 1)

## 2- Fit a max-stable processes
schlather <- fitmaxstab(data, locations, "whitmat", nugget = 0)

## 3- Plot the extremal coefficient
extcoeff(schlather)

---

Fit a copula-based model to spatial extremes

Description

This function fits various copula-based models to spatial extremes data sets.

Usage

fitcopula(data, coord, copula = "gaussian", cov.mod = "whitmat",
loc.form, scale.form, shape.form, marg.cov = NULL, temp.cov = NULL,
temp.form.loc = NULL, temp.form.scale = NULL, temp.form.shape = NULL,
..., start, control = list(maxit = 10000), method = "Nelder", std.err =
TRUE, warn = TRUE, corr = FALSE)

Arguments

data	A matrix representing the data. Each column corresponds to one location.
coord	A matrix that gives the coordinates of each location. Each row corresponds to one location.
copula	A character string. Must be one of "gaussian" and "student" for a Gaussian and Student copula.
cov.mod	A character string corresponding to the correlation function family used in the copula. Must be one of "whitmat", "cauchy", "powexp", "bessel" or "caugen" for the Whittle-Matern, the Cauchy, the Powered Exponential, the Bessel and the Generalized Cauchy correlation families.
loc.form, scale.form, shape.form	R formulas defining the spatial linear model for the GEV parameters. May be missing. See section Details of function fitmaxstab.
marg.cov	Matrix with named columns giving additional covariates for the GEV parameters. If NULL, no extra covariates are used.
temp.cov	Matrix with names columns giving additional *temporal* covariates for the GEV parameters. If NULL, no temporal trend are assume for the GEV parameters — see section Details of function fitmaxstab.
temp.form.loc, temp.form.scale, temp.form.shape	R formulas defining the temporal trends for the GEV parameters. May be missing. See section Details of function fitmaxstab.
...	Several arguments to be passed to the optim, nlm or nlminb functions.
start	A named list giving the initial values for the parameters over which the pairwise likelihood is to be minimized. If start is omitted the routine attempts to find good starting values - but might fail.
control	A list giving the control parameters to be passed to the optim function.
method	The method used for the numerical optimisation procedure. Must be one of BFGS, Nelder-Mead, CG, L-BFGS-B, SANN, nlm or nlminb. See optim for details. Please note that passing nlm or nlminb will use the nlm or nlminb functions instead of optim.
std.err	Logical. Should the standard errors be computed ? The default is to return the standard errors, i.e., std.err = TRUE.
warn	Logical. If TRUE (default), users are warned if the log-likelihood is infinite at starting values and/or problems arised while computing the standard errors.
corr	Logical. If TRUE (non default), the asymptotic correlation matrix is computed.

Value

This function returns a object of class copula.

Warning

This function does not use max-stable copula and the use of non max-stable copula for modelling spatial extreme is highly questionable. This function was mainly implemented for educational purposes and not for concrete modelling purposes.

Author(s)

Mathieu Ribatet

References

Davison, A.C., Padoan, S.A., Ribatet, M. (2010) Statistical Modelling of Spatial Extremes. Submitted to Statistical Science.

See Also

fitmaxstab, latent

Examples

## Not run: 
n.site <- 30
n.obs <- 50

coord <- matrix(runif(2 * n.site, -10, 10), ncol = 2)
colnames(coord) <- c("lon", "lat")

## Generate data from a Gaussian copula model
data <- rcopula(n.obs, coord, "gaussian", "powexp", nugget = 0, range = 4, smooth = 1.2)

## Transform the margins to GEV
locs <- -5 + coord[,"lon"] / 10
scales <- 10 + coord[,"lat"] / 2
shapes <- rep(0.2, n.site)

for (i in 1:n.site)
  data[,i] <- frech2gev(data[,i], locs[i], scales[i], shapes[i])

## Fit a Gaussian copula model

## 1. Define trend surfaces
loc.form <- y ~ lon
scale.form <- y ~ lat
shape.form <- y ~ 1

## 2. Fit
M0 <- fitcopula(data, coord, "gaussian", "powexp", loc.form, scale.form,
                shape.form, nugget = 0)

## End(Not run)

---

Estimates the covariance function for the Schlather's model

Description

Estimates the covariance function for the Schlather's model using non-parametric estimates of the pairwise extremal coefficients.

Usage

fitcovariance(data, coord, cov.mod, marge = "emp", control = list(),
..., start, weighted = TRUE)

Arguments

data	A matrix representing the data. Each column corresponds to one location.
coord	A matrix that gives the coordinates of each location. Each row corresponds to one location.
cov.mod	A character string corresponding the the covariance model in the Schlather's model. Must be one of "whitmat", "cauchy", "powexp", "bessel" or "caugen" for the Whittle-Matern, the Cauchy, the Powered Exponential, the Bessel and the Generalized Cauchy correlation families.
marge	Character string specifying how margins are transformed to unit Frechet. Must be one of "emp", "frech" or "mle" - see function fitextcoeff.
control	The control arguments to be passed to the optim function.
...	Optional arguments to be passed to the optim function.
start	A named list giving the initial values for the parameters over which the weighted sum of square is to be minimized. If start is omitted the routine attempts to find good starting values.
weighted	Logical. Should weighted least squares be used?

Details

The fitting procedure is based on weighted least squares. More precisely, the fitting criteria is to minimize:

$\sum_{i,j} \left(\frac{\tilde{\theta}_{i,j} - \hat{\theta}_{i,j}}{s_{i,j}}\right)^2$

where $\tilde{\theta}_{i,j}$ is a non parametric estimate of the extremal coefficient related to location i and j, $\hat{\theta}_{i,j}$ is the fitted extremal coefficient derived from the Schlather's model and $s_{i,j}$ are the standard errors related to the estimates $\tilde{\theta}_{i,j}$.

Value

An object of class maxstab.

Author(s)

Mathieu Ribatet

References

Smith, R. L. (1990) Max-stable processes and spatial extremes. Unpublished manuscript.

See Also

fitcovmat, fitmaxstab, fitextcoeff

Examples

n.site <- 50
locations <- matrix(runif(2*n.site, 0, 10), ncol = 2)
colnames(locations) <- c("lon", "lat")

##Simulating a max-stable process using RandomFields
##This is the Schlather's approach
data <- rmaxstab(50, locations, cov.mod = "whitmat", nugget = 0, range =
30, smooth = 1)

fitcovariance(data, locations, "whitmat")

---

Estimates the covariance matrix for the Smith's model

Description

Estimates the covariance matrix for the Smith's model using non-parametric estimates of the pairwise extremal coefficients.

Usage

fitcovmat(data, coord, marge = "emp", iso = FALSE, control = list(),
..., start, weighted = TRUE)

Arguments

data	A matrix representing the data. Each column corresponds to one location.
coord	A matrix that gives the coordinates of each location. Each row corresponds to one location.
marge	Character string specifying how margins are transformed to unit Frechet. Must be one of "emp", "frech" or "mle" - see function fitextcoeff.
iso	Logical. If TRUE, isotropy is supposed. Otherwise (default), anisotropy is allowed.
control	The control arguments to be passed to the optim function.
...	Optional arguments to be passed to the optim function.
start	A named list giving the initial values for the parameters over which the weighted sum of square is to be minimized. If start is omitted the routine attempts to find good starting values.
weighted	Logical. Should weighted least squares be used?

Details

The fitting procedure is based on weighted least squares. More precisely, the fitting criteria is to minimize:

$\sum_{i,j} \left(\frac{\tilde{\theta}_{i,j} - \hat{\theta}_{i,j}}{s_{i,j}}\right)^2$

where $\tilde{\theta}_{i,j}$ is a non parametric estimate of the extremal coefficient related to location i and j, $\hat{\theta}_{i,j}$ is the fitted extremal coefficient derived from the Smith's model and $s_{i,j}$ are the standard errors related to the estimates $\tilde{\theta}_{i,j}$.

Value

An object of class maxstab.

Author(s)

Mathieu Ribatet

References

Smith, R. L. (1990) Max-stable processes and spatial extremes. Unpublished manuscript.

See Also

fitcovariance, fitmaxstab, fitextcoeff

Examples

n.site <- 50
n.obs <- 100
locations <- matrix(runif(2*n.site, 0, 40), ncol = 2)
colnames(locations) <- c("lon", "lat")

## Simulate a max-stable process - with unit Frechet margins
data <- rmaxstab(50, locations, cov.mod = "gauss", cov11 = 200, cov12 =
0, cov22 = 200)

fitcovmat(data, locations)

##Force an isotropic model
fitcovmat(data, locations, iso = TRUE)

---

Non parametric estimators of the extremal coefficient function

Description

Estimates non parametrically the extremal coefficient function.

Usage

fitextcoeff(data, coord, ..., estim = "ST", marge = "emp", prob = 0,
plot = TRUE, loess = TRUE, method = "BFGS", std.err = TRUE, xlab,
ylab, angles = NULL, identify = FALSE)

Arguments

data	A matrix representing the data. Each column corresponds to one location.
coord	A matrix that gives the coordinates of each location. Each row corresponds to one location.
...	Additional options to be passed to the plot function.
estim	Character string specifying the estimator to be used. Must be one of "ST" (Schlather and Tawn) or "Smith".
marge	Character string specifying how margins are transformed to unit Frechet. Must be one of "emp", "mle" or "none" - see Details
prob	The probability related to the threshold. Only useful with the ST estimator.
plot	Logical. If TRUE (default), the extremal coefficient function is plotted
loess	If TRUE (default), a local polynomial curve is plotted - see function loess.
method	The optimizer used when fitting the GEV distribution to data. See function gevmle.
std.err	Logical. If TRUE, standard errors are computed. Note that standard errors are not available with the "ST" estimator.
xlab, ylab	The x-axis and y-axis labels. May be missing.
angles	A numeric vector. A partition of the interval $(-\pi, \pi)$ to help detecting anisotropy.
identify	Logical. If TRUE, users can use the identify function to identify pairs of stations on the plot.

Details

During the estimation procedure, data need to be transformed to unit Frechet margins firts. This can be done in two different ways ; by using the empirical CDF or the GEV ML estimates.

If marge = "emp", then the data are transformed using the following relation:

$z_i = - \frac{1}{\log (F(y_i))}$

where $y_i$ are the observations available at location $i$, $F$ is the empirical CDF and $z_i$ are the observations transformed to unit Frechet scale.

If marge = "mle", then the data are transformed using the MLE of the GEV distribution - see function gev2frech.

Lastly, if data are already supposed to be unit Frechet, then no transformation is performed if one passed the option marge = "frech".

If data are already componentwise maxima, prob should be zero. Otherwise, users have to define a threshold $z$ (large enough) where univariate extreme value arguments are relevant. We define prob such that $\Pr[Z \leq z] = prob$.

Value

Plots the extremal coefficient function and returns the points used for the plot. If loess = TRUE, the output is a list with argument "ext.coeff" and "loess".

Author(s)

Mathieu Ribatet

References

Schlather, M. and Tawn, J. A. (2003) A dependence measure for multivariate and spatial extreme values: Properties and inference. Biometrika 90(1):139–156.

Smith, R. L. (1990) Max-stable processes and spatial extremes. Unpublished manuscript.

See Also

madogram

Examples

n.site <- 30
locations <- matrix(runif(2*n.site, 0, 10), ncol = 2)
colnames(locations) <- c("lon", "lat")

##Simulate a max-stable process - with unit Frechet margins
data <- rmaxstab(50, locations, cov.mod = "gauss", cov11 = 10, cov12 =
40, cov22 = 220)

##Plot the extremal coefficient function
op <- par(mfrow=c(1,2))
fitextcoeff(data, locations, estim = "Smith")
fitextcoeff(data, locations, angles = seq(-pi, pi, length = 4), estim = "Smith")
par(op)

---

Fits a max-stable process to data

Description

This function fits max-stable processes to data using pairwise likelihood. Two max-stable characterisations are available: the Smith and Schlather representations.

Usage

fitmaxstab(data,  coord, cov.mod, loc.form, scale.form,
shape.form, marg.cov = NULL, temp.cov = NULL, temp.form.loc = NULL,
temp.form.scale = NULL, temp.form.shape = NULL, iso = FALSE, ...,
fit.marge = FALSE, warn = TRUE, method = "Nelder", start, control =
list(), weights = NULL, corr = FALSE, check.grad
= FALSE)

Arguments

data	A matrix representing the data. Each column corresponds to one location.
coord	A matrix that gives the coordinates of each location. Each row corresponds to one location.
cov.mod	A character string corresponding to the covariance model in the max-stable representation. Must be one of "gauss" for the Smith's model; "whitmat", "cauchy", "powexp", "bessel" or "caugen" for the Whittle-Matern, the Cauchy, the Powered Exponential, the Bessel and the Generalized Cauchy correlation families with the Schlather's model; "brown" for Brown-Resnick processes. The geometric Gaussian and Extremal-t models with a Whittle-Matern correlation function can be fitted by passing respectively "gwhitmat" or "twhitmat". Other correlation function families are considered in a similar way.
loc.form, scale.form, shape.form	R formulas defining the spatial linear model for the GEV parameters. May be missing. See section Details.
marg.cov	Matrix with named columns giving additional covariates for the GEV parameters. If NULL, no extra covariates are used.
temp.cov	Matrix with names columns giving additional *temporal* covariates for the GEV parameters. If NULL, no temporal trend are assume for the GEV parameters — see section Details.
temp.form.loc, temp.form.scale, temp.form.shape	R formulas defining the temporal trends for the GEV parameters. May be missing. See section Details.
iso	Logical. If TRUE an isotropic model is fitted to data. Otherwise (default), anisotropy is allowed. Currently, this is only implemented for the Smith's model.
...	Several arguments to be passed to the optim, nlm or nlminb functions. See section details.
fit.marge	Logical. If TRUE, the GEV parameters are estimated pointwise or using the formulas given by loc.form, scale.form and shape.form. If FALSE, observations are supposed to be unit Frechet distributed. Note that when formulas are given, fit.marge is automatically set to TRUE.
warn	Logical. If TRUE (default), users are warned if the log-likelihood is infinite at starting values and/or problems arised while computing the standard errors.
method	The method used for the numerical optimisation procedure. Must be one of BFGS, Nelder-Mead, CG, L-BFGS-B, SANN, nlm or nlminb. See optim for details. Please note that passing nlm or nlminb will use the nlm or nlminb functions instead of optim.
start	A named list giving the initial values for the parameters over which the pairwise likelihood is to be minimized. If start is omitted the routine attempts to find good starting values - but might fail.
control	A list giving the control parameters to be passed to the optim function.
weights	A numeric vector specifying the weights in the pairwise likelihood - and so has length the number of pairs. If NULL (default), no weighting scheme is used.
corr	Logical. If TRUE (non default), the asymptotic correlation matrix is computed.
check.grad	Logical. If TRUE (non default), the analytic gradient of the pairwise likelihood will be compared to the numerical one. Such a checking might be usefull for ill-conditionned situation diagnosis.

Details

As spatial data often deal with a large number of locations, it is impossible to write analytically the joint distribution. Consequently, the fitting procedure substitutes the "full likelihood" for the pairwise likelihood.

Let define $L_{i,j}(x_{i,j}, \theta)$ the likelihood for site $i$ and $j$, where $i = 1, \dots, N-1$, $j = i+1, \dots, N$, $N$ is the number of site within the region and $x_{i,j}$ are the joint observations for site $i$ and $j$. Then the pairwise likelihood $PL(\theta)$ is defined by:

$\ell_P = \log PL(\theta) = \sum_{i = 1}^{N-1} \sum_{j=i+1}^{N} \log L_{i,j} (x_{i,j}, \theta)$

As pairwise likelihood is an approximation of the “full likelihood”, standard errors cannot be computed directly by the inverse of the Fisher information matrix. Instead, a sandwich estimate must be used to account for model mispecification e.g.

$\hat{\theta} \sim N(\theta, H^{-1} J H^{-1})$

where $H$ is the Fisher information matrix (computed from the misspecified model) and $J$ is the variance of the score function.

There are two different kind of covariates : "spatial" and "temporal".

A "spatial" covariate may have different values accross station but does not depend on time. For example the coordinates of the stations are obviously "spatial". These "spatial" covariates should be used with the marg.cov and loc.form, scale.form, shape.form.

A "temporal" covariates may have different values accross time but does not depend on space. For example the years where the annual maxima were recorded is "temporal". These "temporal" covariates should be used with the temp.cov and temp.form.loc, temp.form.scale, temp.form.shape.

As a consequence note that marg.cov must have K rows (K being the number of sites) while temp.cov must have n rows (n being the number of observations).

Value

This function returns a object of class maxstab. Such objects are list with components:

fitted.values	A vector containing the estimated parameters.
std.err	A vector containing the standard errors.
fixed	A vector containing the parameters of the model that have been held fixed.
param	A vector containing all parameters (optimised and fixed).
deviance	The (pairwise) deviance at the maximum pairwise likelihood estimates.
corr	The correlation matrix.
convergence, counts, message	Components taken from the list returned by optim - for the mle method.
data	The data analysed.
model	The max-stable characterisation used.
fit.marge	A logical that specifies if the GEV margins were estimated.
cov.fun	The estimated covariance function - for the Schlather model only.
extCoeff	The estimated extremal coefficient function.
cov.mod	The covariance model for the spatial structure.

Warning

When using reponse surfaces to model spatially the GEV parameters, the likelihood is pretty rough so that the general purpose optimization routines may fail. It is your responsability to check if the numerical optimization succeeded or not. I tried, as best as I can, to provide warning messages if the optimizers failed but in some cases, no warning will appear!

Author(s)

Mathieu Ribatet

References

Cox, D. R. and Reid, N. (2004) A note on pseudo-likelihood constructed from marginal densities. Biometrika 91, 729–737.

Demarta, S. and McNeil, A. (2005) The t copula and Related Copulas International Statistical Review 73, 111-129.

Gholam–Rezaee, M. (2009) Spatial extreme value: A composite likelihood. PhD Thesis. Ecole Polytechnique Federale de Lausanne.

Kabluchko, Z., Schlather, M. and de Haan, L. (2009) Stationary max-stable fields associated to negative definite functions Annals of Probability 37:5, 2042–2065.

Padoan, S. A. (2008) Computational Methods for Complex Problems in Extreme Value Theory. PhD Thesis. University of Padova.

Padoan, S. A., Ribatet, M. and Sisson, S. A. (2010) Likelihood-based inference for max-stable processes. Journal of the American Statistical Association (Theory and Methods) 105:489, 263-277.

Schlather, M. (2002) Models for Stationary Max-Stable Random Fields. Extremes 5:1, 33–44.

Smith, R. L. (1990) Max-stable processes and spatial extremes. Unpublished.

Examples

## Not run: 
##Define the coordinate of each location
n.site <- 30
locations <- matrix(runif(2*n.site, 0, 10), ncol = 2)
colnames(locations) <- c("lon", "lat")

##Simulate a max-stable process - with unit Frechet margins
data <- rmaxstab(40, locations, cov.mod = "whitmat", nugget = 0, range = 3,
smooth = 0.5)

##Now define the spatial model for the GEV parameters
param.loc <- -10 + 2 * locations[,2]
param.scale <- 5 + 2 * locations[,1] + locations[,2]^2
param.shape <- rep(0.2, n.site)

##Transform the unit Frechet margins to GEV
for (i in 1:n.site)
  data[,i] <- frech2gev(data[,i], param.loc[i], param.scale[i],
param.shape[i])

##Define a model for the GEV margins to be fitted
##shape ~ 1 stands for the GEV shape parameter is constant
##over the region
loc.form <- loc ~ lat
scale.form <- scale ~ lon + I(lat^2)
shape.form <- shape ~ 1

##Fit a max-stable process using the Schlather's model
fitmaxstab(data, locations, "whitmat", loc.form, scale.form,
           shape.form)

## Model without any spatial structure for the GEV parameters
## Be careful this could be *REALLY* time consuming
fitmaxstab(data, locations, "whitmat")

##  Fixing the smooth parameter of the Whittle-Matern family
##  to 0.5 - e.g. considering exponential family. We suppose the data
##  are unit Frechet here.
fitmaxstab(data, locations, "whitmat", smooth = 0.5, fit.marge = FALSE)

##  Fitting a penalized smoothing splines for the margins with the
##     Smith's model
data <- rmaxstab(40, locations, cov.mod = "gauss", cov11 = 100, cov12 =
                 25, cov22 = 220)

##     And transform it to ordinary GEV margins with a non-linear
##     function
fun <- function(x)
  2 * sin(pi * x / 4) + 10
fun2 <- function(x)
  (fun(x) - 7 ) / 15

param.loc <- fun(locations[,2])
param.scale <- fun(locations[,2])
param.shape <- fun2(locations[,1])

##Transformation from unit Frechet to common GEV margins
for (i in 1:n.site)
  data[,i] <- frech2gev(data[,i], param.loc[i], param.scale[i],
param.shape[i])

##Defining the knots, penalty, degree for the splines
n.knots <- 5
knots <- quantile(locations[,2], prob = 1:n.knots/(n.knots+1))
knots2 <- quantile(locations[,1], prob = 1:n.knots/(n.knots+1))

##Be careful the choice of the penalty (i.e. the smoothing parameter)
##may strongly affect the result Here we use p-splines for each GEV
##parameter - so it's really CPU demanding but one can use 1 p-spline
##and 2 linear models.
##A simple linear model will be clearly faster...
loc.form <- y ~ rb(lat, knots = knots, degree = 3, penalty = .5)
scale.form <- y ~ rb(lat, knots = knots, degree = 3, penalty = .5)
shape.form <- y ~ rb(lon, knots = knots2, degree = 3, penalty = .5)

fitted <- fitmaxstab(data, locations, "gauss", loc.form, scale.form, shape.form,
                     control = list(ndeps = rep(1e-6, 24), trace = 10),
                     method = "BFGS")
fitted
op <- par(mfrow=c(1,3))
plot(locations[,2], param.loc, col = 2, ylim = c(7, 14),
     ylab = "location parameter", xlab = "latitude")
plot(fun, from = 0, to = 10, add = TRUE, col = 2)
points(locations[,2], predict(fitted)[,"loc"], col = "blue", pch = 5)
new.data <- cbind(lon = seq(0, 10, length = 100), lat = seq(0, 10, length = 100))
lines(new.data[,1], predict(fitted, new.data)[,"loc"], col = "blue")
legend("topleft", c("true values", "predict. values", "true curve", "predict. curve"),
       col = c("red", "blue", "red", "blue"), pch = c(1, 5, NA, NA), inset = 0.05,
       lty = c(0, 0, 1, 1), ncol = 2)

plot(locations[,2], param.scale, col = 2, ylim = c(7, 14),
     ylab = "scale parameter", xlab = "latitude")
plot(fun, from = 0, to = 10, add = TRUE, col = 2)
points(locations[,2], predict(fitted)[,"scale"], col = "blue", pch = 5)
lines(new.data[,1], predict(fitted, new.data)[,"scale"], col = "blue")
legend("topleft", c("true values", "predict. values", "true curve", "predict. curve"),
       col = c("red", "blue", "red", "blue"), pch = c(1, 5, NA, NA), inset = 0.05,
       lty = c(0, 0, 1, 1), ncol = 2)

plot(locations[,1], param.shape, col = 2,
     ylab = "shape parameter", xlab = "longitude")
plot(fun2, from = 0, to = 10, add = TRUE, col = 2)
points(locations[,1], predict(fitted)[,"shape"], col = "blue", pch = 5)
lines(new.data[,1], predict(fitted, new.data)[,"shape"], col = "blue")
legend("topleft", c("true values", "predict. values", "true curve", "predict. curve"),
       col = c("red", "blue", "red", "blue"), pch = c(1, 5, NA, NA), inset = 0.05,
       lty = c(0, 0, 1, 1), ncol = 2)
par(op)

## End(Not run)

---

MLE for a spatial GEV model

Description

This function derives the MLE of a spatial GEV model.

Usage

fitspatgev(data, covariables, loc.form, scale.form, shape.form,
temp.cov = NULL, temp.form.loc = NULL, temp.form.scale = NULL,
temp.form.shape = NULL, ..., start, control = list(maxit = 10000),
method = "Nelder", warn = TRUE, corr = FALSE)

Arguments

data	A matrix representing the data. Each column corresponds to one location.
covariables	Matrix with named columns giving required covariates for the GEV parameter models.
loc.form, scale.form, shape.form	R formulas defining the spatial models for the GEV parameters. See section Details.
temp.cov	Matrix with names columns giving additional *temporal* covariates for the GEV parameters. If NULL, no temporal trend are assume for the GEV parameters — see section Details.
temp.form.loc, temp.form.scale, temp.form.shape	R formulas defining the temporal trends for the GEV parameters. May be missing. See section Details.
start	A named list giving the initial values for the parameters over which the pairwise likelihood is to be minimized. If start is omitted the routine attempts to find good starting values - but might fail.
...	Several arguments to be passed to the optim functions. See section details.
control	The control argument to be passed to the optim function.
method	The method used for the numerical optimisation procedure. Must be one of BFGS, Nelder-Mead, CG, L-BFGS-B or SANN. See optim for details.
warn	Logical. If TRUE (default), users will be warned if the starting values lie in a zero density region.
corr	Logical. If TRUE, the asymptotic correlation matrix is computed.

Details

A kind of "spatial" GEV model can be defined by using response surfaces for the GEV parameters. For instance, the GEV location parameters are defined through the following equation:

$\mu = X_\mu \beta_\mu$

where $X_\mu$ is the design matrix and $\beta_\mu$ is the vector parameter to be estimated. The GEV scale and shape parameters are defined accordingly to the above equation.

The log-likelihood for the GEV spatial model is consequently defined as follows:

$\ell(\beta) = \sum_{i=1}^{n.site} \sum_{j=1}^{n.obs} \log f(y_{i,j}; \theta_i)$

where $\theta_i$ is the vector of the GEV parameters for the $i$-th site.

Most often, there will be some dependence between stations. However, it can be seen from the log-likelihood definition that we supposed that the stations are mutually independent. Consequently, to get reliable standard error estimates, these standard errors are estimated with their sandwich estimates.

There are two different kind of covariates : "spatial" and "temporal".

A "spatial" covariate may have different values accross station but does not depend on time. For example the coordinates of the stations are obviously "spatial". These "spatial" covariates should be used with the marg.cov and loc.form, scale.form, shape.form.

A "temporal" covariates may have different values accross time but does not depend on space. For example the years where the annual maxima were recorded is "temporal". These "temporal" covariates should be used with the temp.cov and temp.form.loc, temp.form.scale, temp.form.shape.

As a consequence note that marg.cov must have K rows (K being the number of sites) while temp.cov must have n rows (n being the number of observations).

Value

An object of class spatgev. Namely, this is a list with the following arguments:

fitted.values	The parameter estimates.
param	All the parameters e.g. parameter estimates and fixed parameters.
std.err	The standard errors.
var.cov	The asymptotic MLE variance covariance matrix.
counts, message, convergence	Some information about the optimization procedure.
logLik, deviance	The log-likelihood and deviance values.
loc.form, scale.form, shape.form	The formulas defining the spatial models for the GEV parameters.
covariables	The covariables used for the spatial models.
ihessian	The inverse of the Hessian matrix of the negative log-likelihood.
jacobian	The variance covariance matrix of the score.

Author(s)

Mathieu Ribatet

Examples

## 1- Simulate a max-stable random field
n.site <- 35
locations <- matrix(runif(2*n.site, 0, 10), ncol = 2)
colnames(locations) <- c("lon", "lat")
data <- rmaxstab(50, locations, cov.mod = "whitmat", nugget = 0, range =
  3, smooth = 0.5)

## 2- Transformation to non unit Frechet margins
param.loc <- -10 + 2 * locations[,2]
param.scale <- 5 + 2 * locations[,1]
param.shape <- rep(0.2, n.site)
for (i in 1:n.site)
  data[,i] <- frech2gev(data[,i], param.loc[i], param.scale[i],
param.shape[i])

## 3- Fit a ''spatial GEV'' mdoel to data with the following models for
##    the GEV parameters
form.loc <- loc ~ lat
form.scale <- scale ~ lon
form.shape <- shape ~ 1

fitspatgev(data, locations, form.loc, form.scale, form.shape)

---

Computes the F-madogram

Description

Computes the F-madogram for max-stable processes.

Usage

fmadogram(data, coord, fitted, n.bins, which = c("mado", "ext"), xlab,
ylab, col = c(1, 2), angles = NULL, marge = "emp", add = FALSE, xlim =
c(0, max(dist)), ...)

Arguments

data	A matrix representing the data. Each column corresponds to one location.
coord	A matrix that gives the coordinates of each location. Each row corresponds to one location.
fitted	An object of class maxstab - usually the output of the fitmaxstab function. May be missing.
n.bins	The number of bins to be used. If missing, pairwise F-madogram estimates will be computed.
which	A character vector of maximum size 2. It specifies if the madogram and/or the extremal coefficient functions have to be plotted.
xlab, ylab	The x-axis and y-axis labels. May be missing. Note that ylab must have the same length has which.
col	The colors used for the points and optionnaly the fitted curve.
angles	A numeric vector. A partition of the interval $(-\pi, \pi)$ to help detecting anisotropy.
marge	Character string. If 'emp', the probabilities of non exceedances are estimated using the empirical CDF. If 'mle' (default), maximum likelihood estimates are used.
add	Logical. If TRUE, the plot is added to the current figure; otherwhise (default) a new plot is computed.
xlim	A numeric vector of length 2 specifying the x coordinate range.
...	Additional options to be passed to the plot function.

Details

Let $Z(x)$ be a stationary process. The F-madogram is defined as follows:

$\nu(h) = \frac{1}{2}\mbox{E}\left[|F(Z(x+h)) - F(Z(x))| \right]$

The extremal coefficient $\theta(h)$ satisfies:

$\theta(h) = \frac{1 + 2 \nu(h)}{1 - 2 \nu(h)}$

Value

A graphic and (invisibly) a matrix with the lag distances, the F-madogram and extremal coefficient estimates.

Author(s)

Mathieu Ribatet

References

Cooley, D., Naveau, P. and Poncet, P. (2006) Variograms for spatial max-stable random fields. Dependence in Probability and Statistics, 373–390.

See Also

madogram, lmadogram

Examples

n.site <- 15
locations <- matrix(runif(2*n.site, 0, 10), ncol = 2)
colnames(locations) <- c("lon", "lat")

##Simulate a max-stable process - with unit Frechet margins
data <- rmaxstab(40, locations, cov.mod = "whitmat", nugget = 0, range = 1,
smooth = 2)

##Compute the F-madogram
fmadogram(data, locations)

##Compare the F-madogram with a fitted max-stable process
fitted <- fitmaxstab(data, locations, "whitmat", nugget = 0)
fmadogram(fitted = fitted, which = "ext")

---

Estimates the penalty coefficient from the generalized cross-validation criterion

Description

Estimates the penalty coefficient from the generalized cross-validation criterion.

Usage

gcv(y, x, knots, degree, plot = TRUE, n.points = 150, ...)

Arguments

y	The response vector.
x	A vector/matrix giving the values of the predictor variable(s). If x is a matrix, each row corresponds to one observation.
knots	A vector givint the coordinates of the knots.
degree	The degree of the penalized smoothing spline.
plot	Logical. If TRUE (default), the evolution of the CV score as the penalty increases is plotted.
n.points	A numeric giving the number of CV computations needed to produce the plot.
...	Options to be passed to the nlm function.

Details

For every linear smoother e.g. $\hat{y} = S_\lambda y$, the cross-validation criterion consists in minimizing the following quantity:

$GCV(\lambda) = \frac{n ||y - \hat{y}||^2}{(n - tr(S_\lambda))^2}$

where $\lambda$ is the penalty coefficient, $n$ the number of observations and $tr(S_\lambda)$ is the trace of the matrix $S_\lambda$.

Value

A list with components 'penalty', 'gcv' and 'nlm.code' which give the location of the minimum, the value of the cross-validation criterion at that point and the code returned by the link{nlm} function - useful to assess for convergence issues.

Author(s)

Mathieu Ribatet

References

Ruppert, D. Wand, M.P. and Carrol, R.J. (2003) Semiparametric Regression Cambridge Series in Statistical and Probabilistic Mathematics.

See Also

cv

Examples

n <- 200
x <- runif(n)
fun <- function(x) sin(3 * pi * x)
y <- fun(x) + rnorm(n, 0, sqrt(0.4))
knots <- quantile(x, prob = 1:(n/4) / (n/4 + 1))
gcv(y, x, knots = knots, degree = 3)

---

The Generalized Extreme Value Distribution

Description

Density, distribution function, quantile function and random generation for the GP distribution with location equal to 'loc', scale equal to 'scale' and shape equal to 'shape'.

Usage

rgev(n, loc = 0, scale = 1, shape = 0)
pgev(q, loc = 0, scale = 1, shape = 0, lower.tail = TRUE)
qgev(p, loc = 0, scale = 1, shape = 0, lower.tail = TRUE)
dgev(x, loc = 0, scale = 1, shape = 0, log = FALSE)

Arguments

x, q	vector of quantiles.
p	vector of probabilities.
n	number of observations.
loc	vector of the location parameters.
scale	vector of the scale parameters.
shape	a numeric of the shape parameter.
lower.tail	logical; if TRUE (default), probabilities are $\Pr[ X \le x]$, otherwise, $\Pr[X > x]$.
log	logical; if TRUE, probabilities p are given as log(p).

Value

If 'loc', 'scale' and 'shape' are not specified they assume the default values of '0', '1' and '0', respectively.

The GEV distribution function for loc = $u$, scale = $\sigma$ and shape = $\xi$ is

$G(x) = \exp\left[-\left\{1 + \xi \frac{x - u}{\sigma} \right\}^{-1 / \xi} \right]$

for $1 + \xi ( x - u ) / \sigma > 0$ and $x > u$, where $\sigma > 0$. If $\xi = 0$, the distribution is defined by continuity corresponding to the Gumbel distribution.

Examples

dgev(0.1)
rgev(100, 1, 2, 0.2)
qgev(seq(0.1, 0.9, 0.1), 1, 0.5, -0.2)
pgev(12.6, 2, 0.5, 0.1)

---

The Generalized Pareto Distribution

Description

Density, distribution function, quantile function and random generation for the GP distribution with location equal to 'loc', scale equal to 'scale' and shape equal to 'shape'.

Usage

rgpd(n, loc = 0, scale = 1, shape = 0)
pgpd(q, loc = 0, scale = 1, shape = 0, lower.tail = TRUE, lambda = 0)
qgpd(p, loc = 0, scale = 1, shape = 0, lower.tail = TRUE, lambda = 0)
dgpd(x, loc = 0, scale = 1, shape = 0, log = FALSE)

Arguments

x, q	vector of quantiles.
p	vector of probabilities.
n	number of observations.
loc	vector of the location parameters.
scale	vector of the scale parameters.
shape	a numeric of the shape parameter.
lower.tail	logical; if TRUE (default), probabilities are $\Pr[ X \le x]$, otherwise, $\Pr[X > x]$.
log	logical; if TRUE, probabilities p are given as log(p).
lambda	a single probability - see the "value" section.

Value

If 'loc', 'scale' and 'shape' are not specified they assume the default values of '0', '1' and '0', respectively.

The GP distribution function for loc = $u$, scale = $\sigma$ and shape = $\xi$ is

$G(x) = 1 - \left[ 1 + \frac{\xi (x - u )}{ \sigma } \right] ^ { - 1 / \xi}$

for $1 + \xi ( x - u ) / \sigma > 0$ and $x > u$, where $\sigma > 0$. If $\xi = 0$, the distribution is defined by continuity corresponding to the exponential distribution.

By definition, the GP distribution models exceedances above a threshold. In particular, the $G$ function is a suited candidate to model

$\Pr\left[ X \geq x | X > u \right] = 1 - G(x)$

for $u$ large enough.

However, it may be usefull to model the "non conditional" quantiles, that is the ones related to $\Pr[ X \leq x]$. Using the conditional probability definition, one have :

$\Pr\left[ X \geq x \right] = \left(1 - \lambda\right) \left( 1 + \xi \frac{x - u}{\sigma}\right)^{-1/\xi}$

where $\lambda = \Pr[ X \leq u]$.

When $\lambda = 0$, the "conditional" distribution is equivalent to the "non conditional" distribution.

Examples

dgpd(0.1)
rgpd(100, 1, 2, 0.2)
qgpd(seq(0.1, 0.9, 0.1), 1, 0.5, -0.2)
pgpd(12.6, 2, 0.5, 0.1)
##for non conditional quantiles
qgpd(seq(0.9, 0.99, 0.01), 1, 0.5, -0.2, lambda = 0.9)
pgpd(2.6, 2, 2.5, 0.25, lambda = 0.5)

---