tests/testthat/test_intersite.R

In martindurocher/floodStat: FloodNet Regional Frequency Analysis

##############################################################################
# Testing simulation function RegSim
#############################################################################

library(lmomRFA)
rm(list = ls())

## Import data for tests
set.seed(0)

data(canFlood)
lmom5 <- canFlood[1:5, c('L1','LCV','LSK')]
lmom20 <- canFlood[1:20, c('L1','LCV','LSK')]

coord5 <- replicate(2,runif(5))
h5 <- as.matrix(dist(coord5))

coord20 <- replicate(2,runif(20))
h20 <- as.matrix(dist(coord20))

colnames(coord5) <- colnames(coord20) <- c('lon','lat')

## -------------------------------- ##
## Section should work correctly    ##
## -------------------------------- ##

## Define evaluation functions
rxd <- function(x,y) max(abs(1-y/x))
axd <- function(x,y) max(abs(y-x))
rmad <- function(x,y) mean(abs(1-y/x))
mad <- function(x,y) mean(abs(y-x))

## Verify that the right dimension is output
sim5 <- RegSim(lmom5, distr = 'gev', nrec = 3)
sim20 <- RegSim(lmom20, distr = 'pe3', nrec = 7)

expect_equal(dim(sim5), c(3,5))
expect_equal(dim(sim20), c(7,20))
expect_equal(class(sim5), 'matrix')
expect_null(colnames(sim5))
expect_null(rownames(sim5))

## Verify that record length are correct
sim <- RegSim(lmom5, distr = 'gev', nrec = 11:15)
sim.n <- apply(!apply(sim,2,is.na),2, sum)

expect_equal(sim.n, 11:15)

## Verify that the long format work

sim <- RegSim(lmom5, distr = 'gno', nrec = 11:15, long = TRUE)
sim.agg <- aggregate(value ~ site, sim, length)

expect_equal(dim(sim.agg),c(5,2))
expect_equal(sim.agg[,2],11:15)
expect_equal(names(sim.agg), c('site','value'))
expect_equal(class(sim.agg),'data.frame')

## Verify that using constant correlation works
set.seed(1)

fun <- function(rho){
  sim <- RegSim(lmom5, distr = 'gev', nrec = 1e5, corr = rho)
  corr <- cor(sim, method = 'spearman')
  corr <- 2*sin(pi/6*corr)
  axd(corr[lower.tri(corr)], rho)
}

expect_true(fun(0)  < .01)
expect_true(fun(.4) < .01)
expect_true(fun(.9) < .01)


## Verify that correlation coefficients are of correct
set.seed(2)

mat.corr <- exp(-h5)

fun <- function(rho, L = FALSE){
  sim <- RegSim(lmom5, distr = 'gev', nrec = 1e5, corr = rho, corr.sqrt = L)
  corr <- cor(sim, method = 'spearman')
  corr <- 2*sin(pi/6*corr)
  axd(corr,rho)
}

expect_true(fun(mat.corr) <.01)

## Verify that correlation coefficients are correct using corr.sqrt

set.seed(2)
mat.corr.sqrt <- chol(mat.corr)
sim <- RegSim(lmom5, distr = 'gev', nrec = 1e5, corr = mat.corr.sqrt, corr.sqrt = TRUE)
corr <- cor(sim, method = 'spearman')
corr <- 2*sin(pi/6*corr)

expect_true(axd(mat.corr,corr) <.01)

## Verify that the L-moments are correct

set.seed(3)

sim <- RegSim(lmom5, distr = 'pe3', nrec = 1e5)

sim.lmom <- t(apply(sim, 2, samlmu))

expect_true(rmad(lmom5[,1],sim.lmom[,1]) < 0.01)
expect_true(rmad(lmom5[,2],sim.lmom[,2]/sim.lmom[,1]) < 0.01)
expect_true(rmad(lmom5[,3],sim.lmom[,3]) < 0.01)

## Verify that using lscale = TRUE works
set.seed(4)

lmom5.mod <- lmom5
lmom5.mod[,2] <- lmom5[,1] * lmom5[,2]

sim.mod <- RegSim(lmom5.mod, distr = 'gno', nrec = 1e5, lscale = TRUE)

sim.mod.lmom <- t(apply(sim.mod, 2, samlmu))

expect_true(rmad(lmom5[,1],sim.mod.lmom[,1]) < 0.01)
expect_true(rmad(lmom5[,2],sim.mod.lmom[,2]/sim.lmom[,1]) < 0.01)
expect_true(rmad(lmom5[,3],sim.mod.lmom[,3]) < 0.01)

## Verify that passing parameter directly works
set.seed(5)

para5 <- t(apply(lmom5.mod, 1, lmom::pelgno))

sim <- RegSim(para5, distr = 'gno', nrec = 1e5, lmom = FALSE)

sim.lmom <- t(apply(sim, 2, samlmu))

expect_true(rmad(lmom5[,1],sim.lmom[,1]) < 0.01)
expect_true(rmad(lmom5[,2],sim.lmom[,2]/sim.lmom[,1]) < 0.01)
expect_true(rmad(lmom5[,3],sim.lmom[,3]) < 0.01)

## Verify the right distribution is selected
set.seed(6)

lmm <- c(100, 30, .2)
distr <- c('gev','gno','glo','pe3')
para <- rbind(pelgev(lmm),
              pelgno(lmm),
              pelglo(lmm),
              pelpe3(lmm))

lmom <- rbind(lmrgev(para[1,],4),
              lmrgno(para[2,],4),
              lmrglo(para[3,],4),
              lmrpe3(para[4,],4))

sim <- RegSim(para, distr = distr, nrec = 1e5, lmom = FALSE)

sim.lmom <- t(apply(sim, 2, samlmu))
expect_true(rmad(lmom[,1],sim.lmom[,1]) < 0.01)
expect_true(rmad(lmom[,2],sim.lmom[,2]) < 0.01)
expect_true(rmad(lmom[,3],sim.lmom[,3]) < 0.01)
expect_true(rmad(lmom[,4],sim.lmom[,4]) < 0.01)

sim <- RegSim(lmom, distr = distr, nrec = 1e5, lscale = TRUE)
sim.lmom <- t(apply(sim, 2, samlmu))

expect_true(rmad(lmom[,1],sim.lmom[,1]) < 0.01)
expect_true(rmad(lmom[,2],sim.lmom[,2]) < 0.01)
expect_true(rmad(lmom[,3],sim.lmom[,3]) < 0.01)
expect_true(rmad(lmom[,4],sim.lmom[,4]) < 0.01)


## Verify behavior with kappa distribution
set.seed(7)

sim <- RegSim(lmom[-3,], distr = 'kap', nrec = 1e5, lscale = TRUE)

sim.lmom <- t(apply(sim, 2, samlmu))
expect_true(rmad(lmom[-3,1],sim.lmom[,1]) < 0.01)
expect_true(rmad(lmom[-3,2],sim.lmom[,2]) < 0.01)
expect_true(rmad(lmom[-3,3],sim.lmom[,3]) < 0.01)
expect_true(rmad(lmom[-3,4],sim.lmom[,4]) < 0.01)

sim <- RegSim(lmom, distr = c('nor','exp','gum','kap'),
              nrec = 10, lscale = TRUE)
expect_equal(dim(sim), c(10,4))

## -------------------------------- ##
## Section Must raised an error     ##
## -------------------------------- ##

## Verify when length of arguments are not the right format
expect_error(RegSim(lmom5, distr = 'gev', nrec = c(10,10)))
expect_error(RegSim(lmom5, distr = c('gev','glo'), nrec = 10))
expect_error(RegSim(lmom5, distr = 'kap', nrec = 10))

## Verify robustess to missing values
lmom5.mod <- lmom5
lmom5.mod[3,3] <- NA

expect_error(RegSim(lmom5.mod, distr = 'gev', nrec = 10))

expect_error(RegSim(lmom5, distr = c('gum',NA,'gev','pe3',NA), nrec = 10))
expect_error(RegSim(lmom5, distr = NA, nrec = 10))
expect_error(RegSim(lmom5, distr = 'gev', nrec = NA))
expect_error(RegSim(lmom5, distr = 'gev', nrec = c(1:4,NA)))

corr.mod <- matrix(.4,5,5)
diag(corr.mod) <- 1
corr.mod[2,1] <- NA

expect_error(RegSim(lmom5, distr = 'gev', nrec = 10, corr = NA))
expect_error(RegSim(lmom5, distr = 'gev', nrec = 10, corr = corr.mod))

## Correlation must be positive definite
expect_error(RegSim(lmom5.mod, distr = 'gev', nrec = 10, corr = -1))

##############################################################################
# Testing DataWide function
#############################################################################

## basic test
xini <- x0 <- data.frame(1:12, expand.grid(1:3,1:4)[,c(2,1)])
names(xini) <- names(x0) <- c('value','site','time')
x0[,2] <- factor(x0[,2], labels = paste0('s',c(2,1,4,3)))
x0[,3] <- factor(x0[,3], labels = paste0('t',c(2,3,1)))
x0

xmat <- DataWide(x0)
xmat0 <- matrix(c(3,1,2,
                  6,4,5,
                  9,7,8,
                  12,10,11), 3,4)

expect_true( all(xmat == xmat0))

## using formula

x0 <- x0[,c(2,3,1)]
xmat <- DataWide(x0, value~site+time)
expect_true( all(xmat == xmat0))

## Create the missing value when necessary

x1 <- x0
cnew <- xini[,3]
cnew[4] <- 4
x1[,2] <- factor(cnew, labels = paste0('t',c(2,3,1,4)))

xmat <- DataWide(x1, value~site+time)

xmat1 <- matrix(c( 3,  1,  2, NA,
                   6, NA,  5,  4,
                   9,  7,  8, NA,
                  12, 10, 11, NA), 4,4)
id1 <- c(4,6,12,16)
expect_equal(xmat[-id1], xmat1[-id1])
expect_equal(xmat[id1], xmat1[id1])

## if time is a date

x0[,2] <- rep(c(as.Date('2000/1/2'),
                as.Date('2000/1/3'),
                as.Date('2000/1/1')),4)

xmat <- DataWide(x0, value~site+time)
expect_true( all(xmat == xmat0))

## If there is a missing value

x1 <- x0
x1$value[1] <- NA

xmat <- DataWide(x1, value~site+time)
expect_equal(xmat[-2], xmat0[-2])
expect_true(is.na(xmat[2]))

## Raise error for duplicate

x1$site[5] <- 's2'
expect_error(DataWide(x1, value~site+time))


##############################################################################
# Testing intersite function
#############################################################################

set.seed(8)

corr.mat <- exp(-h20)
nrec0 <- sample(10001:10020)

sim <- RegSim(lmom20, distr = 'gev', nrec = nrec0,
              corr = exp(-h20), long = TRUE)

sim$site <- factor(sim$site)
lvl <- paste0('site_',1:20)
levels(sim$site) <- lvl

ifit <- Intersite(value ~ site + time, sim, nmom = 6)


## verify the structure of the ouput
expect_equal(class(ifit),'intersite')
expect_equal(names(ifit), c('lmom','nrec','distance','corr','para'))

expect_equal(class(ifit$lmom),'matrix')
expect_equal(dim(ifit$lmom),c(20,6))
expect_equal(colnames(ifit$lmom), c('L1','LCV','LSK','LKUR','TAU5','TAU6'))
expect_equal(rownames(ifit$lmom), lvl)

expect_equal(class(ifit$nrec),'integer')
expect_equal(names(ifit$nrec),lvl)

expect_equal(colnames(ifit$corr), lvl)
expect_equal(rownames(ifit$corr), lvl)

expect_null(ifit$para)
expect_null(ifit$distance)
expect_true(mad(ifit$nrec, nrec0) < 1e-4)

## verify the L-moments and the correlation estimated
expect_true(rmad(lmom20[,1], ifit$lmom[,1]) < 0.01)
expect_true(mad(lmom20[,2],  ifit$lmom[,2]) < 0.01)
expect_true(mad(lmom20[,3],  ifit$lmom[,3])  < 0.02)
expect_true(mad(corr.mat, ifit$corr) < 0.01)

## verify that the smoothing based on an exponential model works
ifit <- Intersite(value ~ site + time, sim, distance = h20, smooth = TRUE)

expect_true(mad(ifit$corr,(1-ifit$para[1])*exp(-3*h20/ifit$para[2])) < 1e-4)
expect_equal(names(ifit$para),c('nugget', 'range', 'smooth'))
expect_equal(ifit$distance,h20)

## Robustess to missing value

set.seed(10)
sim <- RegSim(lmom5, distr = 'gev', nrec = 20, corr = .5, long = TRUE)

## record with missing value should be removed
sim1 <- sim
sim1[1:2,1] <- NA
expect_true(mad(Intersite(sim1[,c(3,2,1)])$nrec,c(18,20,20,20,20)) == 0)

## Set correlation to zero when pairs does not have enough data
expect_true(mad(Intersite(value ~ site + time, sim, nmin = 21)$corr,
             diag(rep(1,5))) < 1e-4)

## Verify the nmom parameter
set.seed(11)
sim <- RegSim(lmom5, distr = 'gev', nrec = 20, corr = .5, long = TRUE)

for(ii in 0:5){
  ifit <- Intersite(sim, nmom = ii)
  expect_true(ncol(ifit$lmom) == max(2,ii))
}