tests/testthat/test_pool.R

In martindurocher/floodStat: FloodNet Regional Frequency Analysis

##############################################################################
# Testing fitting pooling group function
#############################################################################

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

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

lmm0 <- c(100,30,.3)
lmom20 <- t(replicate(20,lmm0))
coord20 <- replicate(2,runif(20))
h20 <- as.matrix(dist(coord20))
colnames(coord20) <- c('lon','lat')

## 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))

nrec0 <- 10001:10021
sim20 <- RegSim(lmom20, distr = 'gev', nrec = nrec0,
                lscale = TRUE, long = TRUE)

ifit20 <- Intersite(value ~ site + time, sim20, distance = h20)

fit1 <- PoolGroup(ifit20$lmom, ifit20$nrec, distr = 'gno',
                  distance = ifit20$distance[2,], nk = 15)

## Verify that the right station were selected
## indirectly the number of site is verified
sid <- order(h20[2,])[1:15]
expect_true(all(fit1$id == sid))
expect_null(fit1$stat)
expect_true(all(dim(fit1$lmom) == c(15,4)))

## Verify that the right regional L-moment and parameter are found
expect_true(mad(fit1$rlmom[1:3], c(1,.3,.3)) <.01)
expect_true(mad(pelgno(fit1$rlmom),fit1$para) < 1e-4)

## Verify that the record length and total station-years are OK
expect_true(all(fit1$nrec == nrec0[sid]))
expect_equal(fit1$ntot, sum(nrec0[sid]))

## Verify the distance and correlation format.
expect_equal(fit1$distance,h20[2,sid])
expect_equal(fit1$inter.corr,0)
expect_equal(fit1$inter.avg,0)
expect_equal(fit1$inter,'ind')

## verify that right number of lmom is returned
ifit <- Intersite(value~site+time,sim20, distance = h20, nmom = 5)
fit <- PoolGroup(ifit, distance = 1, nk = 15, distr = 'gev')
expect_equal(ncol(fit$lmom), 5)

ifit <- Intersite(value~site+time,sim20, distance = h20, nmom = 2)
expect_error(PoolGroup(ifit, distance = 1, nk = 15, distr = 'gev'))

## Verify that the diagnostics work

nrec0 <- 101:120
sim20 <- RegSim(lmom20, distr = 'gev', corr = .4,
                nrec = nrec0, lscale = TRUE, long = TRUE)

ifit20 <- Intersite(value~site+time, sim20, distance = h20)

fit1 <- PoolGroup(ifit20, distance = 2, nk = 15)

## right format
expect_true(length(fit1$stat)==8)
expect_true(length(fit1$discord)==15)

## right selection
lstd <- c('glo','gev','gno','pe3','gpa')
expect_true(lstd[which.min(abs(fit1$stat[4:8]))] == fit1$distr)

## diagnostics done with the right arguments
fit1 <- PoolGroup(ifit20, distance = 2, nk = 15, distr = 'gev', diagnostic = TRUE)
expect_true(length(fit1$stat)==8)

fit1 <- PoolGroup(ifit20, distance = 2, nk = 15, diagnostic = FALSE)
expect_true(length(fit1$stat)==8)

fit1 <- PoolGroup(ifit20, distance = 2, nk = 15, distr = 'gev')
expect_null(fit1$stat)
expect_null(fit1$discord)
expect_null(fit1$pvalue)

fit1 <- PoolGroup(ifit20, distance = 2, nk = 15, distr = 'gev', pvalue = TRUE)

fit1 <- PoolGroup(ifit20, distance = 2, nk = 15, distr = 'gev',
                  diagnostic = TRUE, pvalue = TRUE,
                  pvalue.nrep = 5, pvalue.nsim = 5)

expect_true(length(fit1$pvalue)==8)

## Verify simulation of the pooling group return the right format
fit1 <- PoolGroup(ifit20, distance = 2, nk = 15, distr = 'gev')

sim <- RegSim(fit1, n=1, long = FALSE)
expect_equal(length(dim(sim)),2)
expect_equal(ncol(sim),15)

nrec1 <- apply(sim,2, function(z) sum(!is.na(z)))
expect_equal(nrec1, fit1$nrec)

sim <- RegSim(fit1, n=2, long = FALSE)
expect_equal(length(dim(sim)),3)
expect_equal(ncol(sim),15)

sim <- RegSim(fit1, n=2, long = TRUE)
expect_equal(length(dim(sim)),2)
expect_equal(ncol(sim),4)

sim <- RegSim(fit1, n=1, long = TRUE)
expect_equal(length(dim(sim)),2)
expect_equal(ncol(sim),3)

## test fitting a pooling group using POT

set.seed(12)

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

para <- cbind(rep(0,5), 1.1, -.1)
para0 <- c(0,.9,-.1)

sim <- RegSim(para, distr = 'gpa', nrec = 1e5, corr = .3, long = TRUE, lmom = FALSE)
ifit <- Intersite(value ~ site + time, sim)

fit1 <- PoolGroup(ifit, distance = h5[1,], nk = 5, method = 'pot')
fit2 <- PoolGroup(ifit, distance = h5[1,], nk = 5, distr = 'gpa')

expect_true(axd(fit1$para,para0) < 0.01)
expect_true(axd(fit2$para,para0) < 0.01)

#########################################################
## Verify the predict.poolgrp function
#########################################################

set.seed(48)

coord20 <- replicate(2,runif(20))
h20 <- as.matrix(dist(coord20))
colnames(coord20) <- c('lon','lat')

para <- cbind(rep(100,20), 30, 0)
sim20 <- RegSim(para, distr = 'gev', nrec = 1e4,
                lscale = TRUE, long = TRUE, lmom = FALSE)

ifit20 <- Intersite(value~site+time, sim20, distance = h20, smooth=TRUE)

fit1 <- PoolGroup(ifit20, distr = 'gev',
                  distance = 2, nk = 15)

q0 <- c(.7,.93)
q1 <- .1
q2 <- c(.5, .8, .9, .95, .98, .99)

hat0 <- lmom::quagev(q0,c(100,30,0))
hat1 <- lmom::quagev(q1,c(100,30,0))
hat2 <- lmom::quagev(q2,c(100,30,0))

## verify output format
out <- predict(fit1, q0)
expect_equal(names(out), c('rt3.3','rt14.3'))
expect_equal(class(out), 'numeric')

out <- predict(fit1,q0, ci = TRUE, nsim = 5)
expect_equal(dim(out),c(2,5))
expect_equal(class(out),'matrix')
expect_equal(colnames(out),c('rt','pred','se','lower','upper'))
expect_null(rownames(out))


## Verify prediction
expect_true(rxd(predict(fit1, q0),hat0) < 0.01)
expect_true(rxd(predict(fit1, q1),hat1) < 0.01)
expect_true(rxd(predict(fit1),hat2) < 0.01)

out <- predict(fit1,q1, ci = TRUE, nsim = 5)
expect_equal(dim(out),c(1,5))


#########################################################
## Verify the removing sites
#########################################################

set.seed(20)

## prepare a dataset
coord20 <- replicate(2,runif(20))
h20 <- as.matrix(dist(coord20))
colnames(coord20) <- c('lon','lat')

para <- cbind(rep(100,20), 30, -.2)
distr <- c(rep('gev',18), rep('pe3',3))
sim20 <- RegSim(para, distr = distr, nrec = 300,
                lscale = TRUE, long = TRUE, lmom = FALSE)

ifit20 <- Intersite(value ~ site + time, sim20, distance = h20)


fit1 <- PoolGroup(ifit20, distance = 2, nk = 20)

## Remove one site
out1 <- PoolRemove(fit1, id = 20)

mid <- which(fit1$id == 20)
expect_equal(out1$id, fit1$id[-mid])
expect_equal(out1$nrec, fit1$nrec[-mid])
expect_equal(out1$lmom, fit1$lmom[-mid,])
expect_equal(out1$distance, fit1$distance[-mid])
expect_equal(length(out1$discord), length(fit1$discord[-mid]))

expect_false(any(out1$para == fit1$para))
expect_true(out1$stat[1] < fit1$stat[1])