wolfCOPtest: Empirically Created Schweizer and Wolff Sigma Test Against an...
In copBasic: General Bivariate Copula Theory and Many Utility Functions
wolfCOPtest R Documentation
Empirically Created Schweizer and Wolff Sigma Test Against an Independence Copula
Description
Perform a Schweizer–Wolff Sigma test for \sigma_\mathbf{C} stemming from \mathrm{NULL} : \mathbf{C} = \mathbf{\Pi} (Independence copula, P), where \sigma_\mathbf{C} is defined under wolfCOP. There are no other known hypothesis tests for \sigma_\mathbf{C} in the literature.
This Schweizer–Wolff Sigma hypothesis test has been empirically tuned from small to large samples by extensive simulation, and, as implemented, the test is closely tied to the version of sample computation by wolfCOP using a particular form of the empirical copula (refer to that function for details). The wolfCOPtest function internally loads the wolfCOPtest_data_smlsam tables, and most users generally will not need to access that dataset themselves unless perhaps for checking veracity of the simulations.
Usage
wolfCOPtest(x, y, asuv=FALSE, aslist=TRUE, na.rm=TRUE, digits=6,
probs=c(0.90, 0.95, 0.98, 0.99, 0.995),
zmat=NULL, statf=mean, rndphi=20, usepade=FALSE,
ties.method=c("average", "first", "last", "random", "max", "min"), ...)
Arguments
x
Nonexceedance probabilities u in the X direction if asuv=TRUE, else data in the X direction for which nonexceedance probabilities will be computed by plotting-position formula. If x is a single value, it is assumed to be a previously computed \sigma_\mathbf{C} and y will be considered as the sample size. Note, statistical estimation is based on samples of n \ge 3, if the sample size is less than this, then internally n = 3 is set for the computations. The x can also be a two column data structure and the y will be taken from the second column;
y
Nonexceedance probabilities v in the Y direction if asuv=TRUE, else data in the Y direction for which nonexceedance probabilities will be computed by plotting-position formula (refer to lmomco::pp() function, but simply the
rank(...) function is used);
asuv
Logical identifying whether the x and y are already in nonexceedance probabilities or rather such need to be computed internally by plotting-position formula (refer to lmomco::pp() function, but simply the rank(...) function is used);
aslist
Logical to return a list of the test results rather than a vector;
na.rm
A logical to trigger a complete.cases operation on internal matrix of the data as probabilities;
digits
Digits for rounding before return;
probs
Quantiles by nonexceedance probability to extract from the fitted distribution—these are common upper-tail critical values. If probs is empty, then 0.95 nonexceedance probability is silently used;
zmat
An optional matrix of four columns and rows equally those of the data values with columns 1 and 2 containing, respectively, the lower and upper limits of the u and columns 3 and 4 containing, respectively, the lower and upper limits of v. This zmat feature provides for support if data censoring is involved and simulation between the limits made (refer to Examples). Note, if zmat is given, then ties.method is set "random" and asuv is set false;
statf
Statistical function for reducing simulated \sigma values (mean or median suggested), but a user could plot(density(result$rand_sigma)) and reach their own decisions on centrality measure;
rndphi
A multiplier on the number of ties or censored values to drive simulation of samples as part of populating the rand_sigma element, when applicable, on the function return;
usepade
A logical to switch, on estimation of the empirical L-moments of the logit transformed simulations, between (1) power-exponent optimized regression equations, or (2) optimized Padé Approximant forms;
ties.method
Argument of the same name for rank, but if there are ties, asuv is set false, and ties.method is set to random, then simulation is made internally on the ranking and statf used to compute the \sigma estimate;
...
Additional arguments to pass to the lmomco::pp() function.
Details
The implementation here is based on Monte Carlo simulations of \mathbf{\Pi} and computation of the L-moments of a logit transform of the computed \sigma_\mathbf{\Pi} for which nonlinear regression of the L-moments to the sample size have been created. Development scripts for the simulations are located in the directory inst/make_wolfCOPtest/ of copBasic sources.
The logit transformation is defined as
\mathrm{logit}\bigl(\sigma_\mathbf{C}\bigr) = \mathrm{log}\frac{\sigma_\mathbf{C}}{1-\sigma_\mathbf{C}} = \eta\mbox{,}
and its inverse is defined as
\sigma_\mathbf{C} = \mathrm{logit}^{(-1)}(\eta) = \frac{1}{1 + \mathrm{exp}(-\eta)}\mbox{.}
Sample-sample empirical distributions (n \in 3, \cdots, 40) will have p-values reported on output and missing for n \ge 41 (refer to wolfCOPtest_data_smlsam for more details).
The predicted L-moments are used to drive from the lmomco package the Pearson type III (PE3) (parpe3, cdfpe3, quape3) of the logit transform of \sigma_\mathbf{\Pi} to compute the p-value and supporting information.
Evaluation of an L-moment ratio diagram (refer to lmomco functions plotlmrdia() and lmrdia()) shows a change from a generalized normal distribution (GNO) to PE3 at about n = 40, but inspection of upper-tail empirical distributions (and hence p-values in the practical range for \mathrm{NULL} rejection) suggest the PE3 tail is better choice for overall application. The wolfCOPtest function, therefore, exclusively uses the PE3.
To conclude, the figure below shows the L-moment ratio diagram of the L-skew and L-kurtosis diagram for the logit transform of the L-moments of the \hat\sigma_\mathbf{C} based on EMPIRgrid_fast in the aforementioned build directory of the copBasic package. The figure is based on nearly 900 million simulations summed across a representative range of sample sizes n \in 8, \cdots 15001, and provides support for the decision by the author to use the PE3 as the distribution model for the hypothesis test.
Value
Extensive list or named vector output the test is returned. For the list, the estimate contains duplicate information as statistic[1]; this is so that there is an estimate entry to parallel the return of cor.test(...) to help users write nomenclature similar calls between the functions.
Author(s)
W.H. Asquith
See Also
wolfCOP, wolfCOPsamc, wolfCOPtest_data_smlsam
Examples
uv <- simCOP(n=62, cop=PSP, graphics=FALSE, seed=8)
# Approx p-value 0.001 PSP is not independence.
# There is no p.value(sample_le40) entry on the return because too large sample size.
wolfCOPtest(uv[,1], uv[,2], asuv=TRUE, aslist=FALSE, usepade=FALSE)
wolfCOPtest(uv[,1], uv[,2], asuv=TRUE, aslist=FALSE, usepade=TRUE )
# See how the p-values between the empirical approximations are about the same.
uv <- simCOP(n=62, cop=P, graphics=FALSE, seed=8)
wolfCOPtest(uv[,1], uv[,2], asuv=TRUE, aslist=FALSE) # approx p-value 0.200***
# Note, the small sample p-values are NA for two previous examples because the sample
# size has become "large", but in the next example, we can inspect two p-values
# (1) by the Pearson type III distribution to the predicted L-moments for the sample
# size and (2) by the small-sample empirical tables in object wolfCOPtest_data_smlsam.
uv <- simCOP(n=13, cop=P, graphics=FALSE, seed=81)
wolfCOPtest(uv[,1], uv[,2], asuv=TRUE, aslist=FALSE)
# p.value(dist_pe3) p.value(sample_le40) # four signif. figures kept here as
# 0.2550*** 0.25144** # simulations might tune up internals
## Not run:
# This example needs larger nsim but will show that the Sigma is performing close to
# Spearman Rho, which helps establish the veracity of the algorithms. If the user were
# to switch out the copula P and replace with CIRCcop (circular copula), then one will
# see Sigma performing similarly but Rho failing because the CIRCcop is an example of
# perfect dependence in a setting of independence. Use of much larger nsim for this
# example will hone the apparent rejection to about 0.05 for P, but clearly we see
# divergence in performance for CIRCcop.
alpha <- 0.05; n <- 42; nsim <- 20000; cop <- P; p <- NULL
set.seed(1)
wolf <- replicate(nsim, { wolfCOP(para=rCOP(n=n, cop, para=p),
as.sample=TRUE, nlarge=0, usefastgrid=TRUE) })
wolfpv <- sapply(wolf, function(w) {
wolfCOPtest(w, n, asuv=TRUE, aslist=FALSE)[4]})
message(" Schweizer-Wolff Sigma(n=", n, ") NULL=P rejection rate = ",
sum(wolfpv < alpha) / nsim)
set.seed(1)
rhopv <- replicate(nsim, { uv <- rCOP(n=n, cop=cop, para=p);
cor.test(uv[,1], uv[,2], method="spearman")$p.value })
message(" Spearman Rho(n=", n, ") NULL=P rejection rate = ",
sum(rhopv < alpha) / nsim)
# For nsim <- 20000
# P : Schweizer-Wolff Sigma(n=42) NULL=P rejection rate = 0.04935
# P : Spearman Rho(n=42) NULL=P rejection rate = 0.04805
# CIRCcop : Schweizer-Wolff Sigma(n=42) NULL=P rejection rate = 0.28125
# CIRCcop : Spearman Rho(n=42) NULL=P rejection rate = 0.00100
# Such example computations show that Rho does not "sense" the perfect dependence
# of CIRCcop because its "correlation" is zero but has obvious dependence when
# plotted. Therefore, Sigma has a much higher rejection rate that will go to unity
# as the sample size increase.
## End(Not run)
## Not run:
lmomco::plotlmrdia(lmomco::lmrdia(usrtrim=TRUE), nopoints=TRUE, xleg="topleft",
xlim=c(0.11,0.32), ylim=c(0.12,0.20), autolegend=TRUE,
expand.names=TRUE, autoaxes=TRUE, xaxs="i", yaxs="i", lwd.cex=1.3)
xy <- st <- NULL
ns <- sort(unique(c(10^seq(0,1,0.01), 10^seq(0, 4, by=0.02))))
for(n in ns) {
lmr <- wolfCOPtest(0, n, usepade=FALSE)$lmoms_logit_sigma; if(is.null(lmr)) next
xy <- rbind(xy, data.frame(n=n, mu=lmr[1], lam2=lmr[2], tau3=lmr[3], tau4=lmr[4]))
}
xy <- xy[par()$usr[1] <= xy$tau3 & xy$tau3 <= par()$usr[2],]
xy <- xy[par()$usr[3] <= xy$tau4 & xy$tau4 <= par()$usr[4],]
lines(xy$tau3, xy$tau4, col="grey60", lwd=7)
for(n in ns) {
lmr <- wolfCOPtest(0, n, usepade=TRUE)$lmoms_logit_sigma; if(is.null(lmr)) next
st <- rbind(st, data.frame(n=n, mu=lmr[1], lam2=lmr[2], tau3=lmr[3], tau4=lmr[4]))
}
st <- st[par()$usr[1] <= st$tau3 & st$tau3 <= par()$usr[2],]
st <- st[par()$usr[3] <= st$tau4 & st$tau4 <= par()$usr[4],]
lines(st$tau3, st$tau4, col="red1", lwd=2)
txt <- c("Logit transform sampling Sigma distribution, wolfCOPtest(usepade=FALSE)",
"Logit transform sampling Sigma distribution, wolfCOPtest(usepade=TRUE )")
legend("bottomright", txt, lwd=c(7,2), col=c("grey50", "red1"), bty="n", cex=0.9) #
## End(Not run)
## Not run:
# Here is a cross check in properties of the mathematics computing W.
# The examples here also show the pass of 2-columns in x; so, a local copy
# of a simulations is not needed before passing to wolfCOPtest().
alpha <- 0.05; n <- 42; nsim <- 2000; cop <- P; p <- NULL
set.seed(1)
wolf1 <- replicate(nsim, { wolfCOP(para=rCOP(n=n, cop, para=p),
as.sample=TRUE, nlarge=0, usefastgrid=TRUE) })
wolf1pv <- sapply(wolf1, function(w) { # here, sigma precomputed
wolfCOPtest(w, n, asuv=TRUE, aslist=FALSE)[4]})
set.seed(1)
wolf2pv <- replicate(nsim, { # here the "true" probabilities are known
wolfCOPtest(rCOP(n=n, cop, para=p), asuv=TRUE )$p.value[1]})
set.seed(1)
wolf3pv <- replicate(nsim, { # here the probabilities computed internally
wolfCOPtest(rCOP(n=n, cop, para=p), asuv=FALSE)$p.value[1]})
message(" Schweizer-Wolff Sigma(n=", n, ") NULL: P rejection rate = ",
sum(wolf1pv < alpha) / nsim)
message(" Schweizer-Wolff Sigma(n=", n, ") NULL: P rejection rate = ",
sum(wolf2pv < alpha) / nsim)
message(" Schweizer-Wolff Sigma(n=", n, ") NULL: P rejection rate = ",
sum(wolf3pv < alpha) / nsim) #
## End(Not run)
## Not run:
# Here is a cross check of sorts between wolfCOPtest and the Independence test
# available in the copula package. We use a Clayton copula with small bivariate
# association and will see for this circumstance that wolfCOPtest out performs
# copula::indepTest by just a small amount by rejecting the NULL more often.
sample_size <- 100; ALPHA <- 0.05
THECOP <- CLcop; THEPARA <- 0.20; rhoCOP(cop=THECOP, para=0.20) # 0.1358466
nsim <- 2000; digits <- as.integer( floor( log10(nsim) ) )
J <- copula::indepTestSim(sample_size, 2, N=100000)
W <- I <- NULL
for(i in seq_len(nsim)) {
if(is.null(W)) {
uv <- copBasic::simCOP( sample_size, cop=THECOP, para=THEPARA )
mtext("Copula simulation reference figure", line=0.5, font=2 )
} else {
uv <- copBasic::rCOP( sample_size, cop=THECOP, para=THEPARA)
}
I <- c(I, copula::indepTest( uv, J )$pvalues[1])
W <- c(W, copBasic::wolfCOPtest(uv[,1], uv[,2], asuV=TRUE )$p.value[1])
}
par(lend=2, las=1); hyphen <- "\u00ad"
plot( sort( rank(I) ) / ( length(I) + 1 ), sort(I), type="l", xaxs="i", yaxs="i",
xlab="Nonexceedance probability", lwd=3, col="grey50",
ylab=paste0("P", hyphen, "value of respective test of independence"))
lines(sort( rank(W) ) / ( length(W) + 1 ), sort(W), col="darkorange1")
txt <- c(paste0("Independence test (copula::indepTest) has rejection rate = ",
round( sum(I < ALPHA) / nsim, digits=digits), ""),
paste0("Schweizer-Wolff test (copBasic::wolfCOPtest) has rejection rate = ",
round( sum(W < ALPHA) / nsim, digits=digits), ""))
legend("topleft", txt, bty="o", box.lty=0, bg=NA, inset=0.02, cex=0.9,
lwd=c(3,1), col=c("grey50", "darkorange1"))
mtext(paste0("Empirical distribution of p", hyphen, "values"), line=0.5, font=2) #
## End(Not run)
## Not run:
# DATA CENSORING EXAMPLE ------------------------------------------------------------
set.seed(10) # Example using the zmat feature for support of censoring by simulation.
muslp <- 0.005; sdslp <- 0.019 # drift with exp used to keep parameter positive
yy <- 1980:2021; n <- length(yy) # Simulated log-Pearson III annual peak streamflows
qq <- 10^sapply(yy, function(y) { # put a small drift in the log10 mean and variation
para <- lmomco::vec2par(c(3.5+(y-min(yy))*muslp,
exp(log(0.5)-(y-min(yy))*sdslp), 0.3), type="pe3")
lmomco::rlmomco(1, para) }) # flow
zul <- zur <- rep(NA, n) # initialize the lower and upper limits for years
zvl <- zvr <- rep(NA, n) # initialize the lower and upper limits for flows
# NA years are not simulated and the data are used as is
# Make some highly contrived censoring for the 3, 6, 9, 13, and 35+th observations.
zul[ 3] <- yy[ 3] - 1.1; zur[ 3] <- yy[ 3] + 1.1
zul[13] <- yy[13] - 1.5; zur[13] <- yy[13] + 2.4
zvl[ 3] <- 10; zvr[3] <- qq[3]
zvl[ 6] <- qq[6]; zvr[6] <- max(qq) + sqrt(.Machine$double.eps) # nudge above max
zvl[ 9] <- 10^(log10(qq[9]) - 0.3); zvr[9] <- 10^(log10(qq[9]) + 0.2) # log-cycles
qq[ 9] <- sqrt( zvl[9] * zvr[9] ) # geometric mean
zvl[(n-6):n] <- 10; zvr[(n-6):n] <- qq[(n-6):n]
xy <- cbind(yy, qq); zmat <- cbind(zul, zur, zvl, zvr)
swolf <- wolfCOPtest(xy[,1], xy[,2], zmat=zmat, asuv=FALSE, statf=median, rndphi=500)
ylim <- log10( range( qq ) ); ylim <- 10^c( floor( ylim[1] ), ceiling( ylim[2] ) )
plot( yy, qq, type="n", log="y", xlab="Water year", ylab="Streamflow", ylim=ylim)
segments(x=zul, x1=zur, y=qq, y1=qq, col="darkgoldenrod3")
segments(x=yy, x1=yy, y=zvl, y1=zvr, col="darkgoldenrod3")
points(yy, qq, bg="darkgoldenrod1", col="darkgoldenrod4", pch=21)
naive <- wolfCOPtest(xy[,1], xy[,2], asuv=FALSE, ties.method="random")
txt <- paste0("Sigma simulation (nsim=",
prettyNum(length(swolf$rand_sigma), big.mark=","), ") ",
"p.value=", round(swolf$p.value[1], digits=3), "; ",
"Naive p.value=", round(naive$p.value[1], digits=3))
mtext(txt, line=1) # Visually the point values tell us to expect independence to be
# rejected, but when the data censoring is taken into account and simulation used,
# we do not reject independence that is heavily influenced by the last seven data
# values on the right being censored down through the bottom axis.
swolf$rand_sigma <- sort(swolf$rand_sigma)
p <- sapply(swolf$rand_sigma, function(k) { # distribution of p-values
n <- ifelse(swolf$sample_size <= 40, 2, 1)
wolfCOPtest(k, swolf$sample_size)$p.value[n] })
swolf$rand_p.values <- p
# plot(swolf$rand_sigma, swolf$rand_p.values, pch=16, col=grey(0, 0.05)) #
## End(Not run)
copBasic documentation built on July 23, 2026, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| wolfCOPtest | R Documentation |
Empirically Created Schweizer and Wolff Sigma Test Against an Independence Copula
Description
Perform a Schweizer–Wolff Sigma test for \sigma_\mathbf{C} stemming from \mathrm{NULL} : \mathbf{C} = \mathbf{\Pi} (Independence copula, P), where \sigma_\mathbf{C} is defined under wolfCOP. There are no other known hypothesis tests for \sigma_\mathbf{C} in the literature.
This Schweizer–Wolff Sigma hypothesis test has been empirically tuned from small to large samples by extensive simulation, and, as implemented, the test is closely tied to the version of sample computation by wolfCOP using a particular form of the empirical copula (refer to that function for details). The wolfCOPtest function internally loads the wolfCOPtest_data_smlsam tables, and most users generally will not need to access that dataset themselves unless perhaps for checking veracity of the simulations.
Usage
wolfCOPtest(x, y, asuv=FALSE, aslist=TRUE, na.rm=TRUE, digits=6,
probs=c(0.90, 0.95, 0.98, 0.99, 0.995),
zmat=NULL, statf=mean, rndphi=20, usepade=FALSE,
ties.method=c("average", "first", "last", "random", "max", "min"), ...)
Arguments
x |
Nonexceedance probabilities |
y |
Nonexceedance probabilities |
asuv |
Logical identifying whether the |
aslist |
Logical to return a list of the test results rather than a vector; |
na.rm |
A logical to trigger a |
digits |
Digits for rounding before return; |
probs |
Quantiles by nonexceedance probability to extract from the fitted distribution—these are common upper-tail critical values. If |
zmat |
An optional matrix of four columns and rows equally those of the data values with columns 1 and 2 containing, respectively, the lower and upper limits of the |
statf |
Statistical function for reducing simulated |
rndphi |
A multiplier on the number of ties or censored values to drive simulation of samples as part of populating the |
usepade |
A logical to switch, on estimation of the empirical L-moments of the logit transformed simulations, between (1) power-exponent optimized regression equations, or (2) optimized Padé Approximant forms; |
ties.method |
Argument of the same name for |
... |
Additional arguments to pass to the |
Details
The implementation here is based on Monte Carlo simulations of \mathbf{\Pi} and computation of the L-moments of a logit transform of the computed \sigma_\mathbf{\Pi} for which nonlinear regression of the L-moments to the sample size have been created. Development scripts for the simulations are located in the directory inst/make_wolfCOPtest/ of copBasic sources.
The logit transformation is defined as
\mathrm{logit}\bigl(\sigma_\mathbf{C}\bigr) = \mathrm{log}\frac{\sigma_\mathbf{C}}{1-\sigma_\mathbf{C}} = \eta\mbox{,}
and its inverse is defined as
\sigma_\mathbf{C} = \mathrm{logit}^{(-1)}(\eta) = \frac{1}{1 + \mathrm{exp}(-\eta)}\mbox{.}
Sample-sample empirical distributions (n \in 3, \cdots, 40) will have p-values reported on output and missing for n \ge 41 (refer to wolfCOPtest_data_smlsam for more details).
The predicted L-moments are used to drive from the lmomco package the Pearson type III (PE3) (parpe3, cdfpe3, quape3) of the logit transform of \sigma_\mathbf{\Pi} to compute the p-value and supporting information.
Evaluation of an L-moment ratio diagram (refer to lmomco functions plotlmrdia() and lmrdia()) shows a change from a generalized normal distribution (GNO) to PE3 at about n = 40, but inspection of upper-tail empirical distributions (and hence p-values in the practical range for \mathrm{NULL} rejection) suggest the PE3 tail is better choice for overall application. The wolfCOPtest function, therefore, exclusively uses the PE3.
To conclude, the figure below shows the L-moment ratio diagram of the L-skew and L-kurtosis diagram for the logit transform of the L-moments of the \hat\sigma_\mathbf{C} based on EMPIRgrid_fast in the aforementioned build directory of the copBasic package. The figure is based on nearly 900 million simulations summed across a representative range of sample sizes n \in 8, \cdots 15001, and provides support for the decision by the author to use the PE3 as the distribution model for the hypothesis test.
Value
Extensive list or named vector output the test is returned. For the list, the estimate contains duplicate information as statistic[1]; this is so that there is an estimate entry to parallel the return of cor.test(...) to help users write nomenclature similar calls between the functions.
Author(s)
W.H. Asquith
See Also
wolfCOP, wolfCOPsamc, wolfCOPtest_data_smlsam
Examples
uv <- simCOP(n=62, cop=PSP, graphics=FALSE, seed=8)
# Approx p-value 0.001 PSP is not independence.
# There is no p.value(sample_le40) entry on the return because too large sample size.
wolfCOPtest(uv[,1], uv[,2], asuv=TRUE, aslist=FALSE, usepade=FALSE)
wolfCOPtest(uv[,1], uv[,2], asuv=TRUE, aslist=FALSE, usepade=TRUE )
# See how the p-values between the empirical approximations are about the same.
uv <- simCOP(n=62, cop=P, graphics=FALSE, seed=8)
wolfCOPtest(uv[,1], uv[,2], asuv=TRUE, aslist=FALSE) # approx p-value 0.200***
# Note, the small sample p-values are NA for two previous examples because the sample
# size has become "large", but in the next example, we can inspect two p-values
# (1) by the Pearson type III distribution to the predicted L-moments for the sample
# size and (2) by the small-sample empirical tables in object wolfCOPtest_data_smlsam.
uv <- simCOP(n=13, cop=P, graphics=FALSE, seed=81)
wolfCOPtest(uv[,1], uv[,2], asuv=TRUE, aslist=FALSE)
# p.value(dist_pe3) p.value(sample_le40) # four signif. figures kept here as
# 0.2550*** 0.25144** # simulations might tune up internals
## Not run:
# This example needs larger nsim but will show that the Sigma is performing close to
# Spearman Rho, which helps establish the veracity of the algorithms. If the user were
# to switch out the copula P and replace with CIRCcop (circular copula), then one will
# see Sigma performing similarly but Rho failing because the CIRCcop is an example of
# perfect dependence in a setting of independence. Use of much larger nsim for this
# example will hone the apparent rejection to about 0.05 for P, but clearly we see
# divergence in performance for CIRCcop.
alpha <- 0.05; n <- 42; nsim <- 20000; cop <- P; p <- NULL
set.seed(1)
wolf <- replicate(nsim, { wolfCOP(para=rCOP(n=n, cop, para=p),
as.sample=TRUE, nlarge=0, usefastgrid=TRUE) })
wolfpv <- sapply(wolf, function(w) {
wolfCOPtest(w, n, asuv=TRUE, aslist=FALSE)[4]})
message(" Schweizer-Wolff Sigma(n=", n, ") NULL=P rejection rate = ",
sum(wolfpv < alpha) / nsim)
set.seed(1)
rhopv <- replicate(nsim, { uv <- rCOP(n=n, cop=cop, para=p);
cor.test(uv[,1], uv[,2], method="spearman")$p.value })
message(" Spearman Rho(n=", n, ") NULL=P rejection rate = ",
sum(rhopv < alpha) / nsim)
# For nsim <- 20000
# P : Schweizer-Wolff Sigma(n=42) NULL=P rejection rate = 0.04935
# P : Spearman Rho(n=42) NULL=P rejection rate = 0.04805
# CIRCcop : Schweizer-Wolff Sigma(n=42) NULL=P rejection rate = 0.28125
# CIRCcop : Spearman Rho(n=42) NULL=P rejection rate = 0.00100
# Such example computations show that Rho does not "sense" the perfect dependence
# of CIRCcop because its "correlation" is zero but has obvious dependence when
# plotted. Therefore, Sigma has a much higher rejection rate that will go to unity
# as the sample size increase.
## End(Not run)
## Not run:
lmomco::plotlmrdia(lmomco::lmrdia(usrtrim=TRUE), nopoints=TRUE, xleg="topleft",
xlim=c(0.11,0.32), ylim=c(0.12,0.20), autolegend=TRUE,
expand.names=TRUE, autoaxes=TRUE, xaxs="i", yaxs="i", lwd.cex=1.3)
xy <- st <- NULL
ns <- sort(unique(c(10^seq(0,1,0.01), 10^seq(0, 4, by=0.02))))
for(n in ns) {
lmr <- wolfCOPtest(0, n, usepade=FALSE)$lmoms_logit_sigma; if(is.null(lmr)) next
xy <- rbind(xy, data.frame(n=n, mu=lmr[1], lam2=lmr[2], tau3=lmr[3], tau4=lmr[4]))
}
xy <- xy[par()$usr[1] <= xy$tau3 & xy$tau3 <= par()$usr[2],]
xy <- xy[par()$usr[3] <= xy$tau4 & xy$tau4 <= par()$usr[4],]
lines(xy$tau3, xy$tau4, col="grey60", lwd=7)
for(n in ns) {
lmr <- wolfCOPtest(0, n, usepade=TRUE)$lmoms_logit_sigma; if(is.null(lmr)) next
st <- rbind(st, data.frame(n=n, mu=lmr[1], lam2=lmr[2], tau3=lmr[3], tau4=lmr[4]))
}
st <- st[par()$usr[1] <= st$tau3 & st$tau3 <= par()$usr[2],]
st <- st[par()$usr[3] <= st$tau4 & st$tau4 <= par()$usr[4],]
lines(st$tau3, st$tau4, col="red1", lwd=2)
txt <- c("Logit transform sampling Sigma distribution, wolfCOPtest(usepade=FALSE)",
"Logit transform sampling Sigma distribution, wolfCOPtest(usepade=TRUE )")
legend("bottomright", txt, lwd=c(7,2), col=c("grey50", "red1"), bty="n", cex=0.9) #
## End(Not run)
## Not run:
# Here is a cross check in properties of the mathematics computing W.
# The examples here also show the pass of 2-columns in x; so, a local copy
# of a simulations is not needed before passing to wolfCOPtest().
alpha <- 0.05; n <- 42; nsim <- 2000; cop <- P; p <- NULL
set.seed(1)
wolf1 <- replicate(nsim, { wolfCOP(para=rCOP(n=n, cop, para=p),
as.sample=TRUE, nlarge=0, usefastgrid=TRUE) })
wolf1pv <- sapply(wolf1, function(w) { # here, sigma precomputed
wolfCOPtest(w, n, asuv=TRUE, aslist=FALSE)[4]})
set.seed(1)
wolf2pv <- replicate(nsim, { # here the "true" probabilities are known
wolfCOPtest(rCOP(n=n, cop, para=p), asuv=TRUE )$p.value[1]})
set.seed(1)
wolf3pv <- replicate(nsim, { # here the probabilities computed internally
wolfCOPtest(rCOP(n=n, cop, para=p), asuv=FALSE)$p.value[1]})
message(" Schweizer-Wolff Sigma(n=", n, ") NULL: P rejection rate = ",
sum(wolf1pv < alpha) / nsim)
message(" Schweizer-Wolff Sigma(n=", n, ") NULL: P rejection rate = ",
sum(wolf2pv < alpha) / nsim)
message(" Schweizer-Wolff Sigma(n=", n, ") NULL: P rejection rate = ",
sum(wolf3pv < alpha) / nsim) #
## End(Not run)
## Not run:
# Here is a cross check of sorts between wolfCOPtest and the Independence test
# available in the copula package. We use a Clayton copula with small bivariate
# association and will see for this circumstance that wolfCOPtest out performs
# copula::indepTest by just a small amount by rejecting the NULL more often.
sample_size <- 100; ALPHA <- 0.05
THECOP <- CLcop; THEPARA <- 0.20; rhoCOP(cop=THECOP, para=0.20) # 0.1358466
nsim <- 2000; digits <- as.integer( floor( log10(nsim) ) )
J <- copula::indepTestSim(sample_size, 2, N=100000)
W <- I <- NULL
for(i in seq_len(nsim)) {
if(is.null(W)) {
uv <- copBasic::simCOP( sample_size, cop=THECOP, para=THEPARA )
mtext("Copula simulation reference figure", line=0.5, font=2 )
} else {
uv <- copBasic::rCOP( sample_size, cop=THECOP, para=THEPARA)
}
I <- c(I, copula::indepTest( uv, J )$pvalues[1])
W <- c(W, copBasic::wolfCOPtest(uv[,1], uv[,2], asuV=TRUE )$p.value[1])
}
par(lend=2, las=1); hyphen <- "\u00ad"
plot( sort( rank(I) ) / ( length(I) + 1 ), sort(I), type="l", xaxs="i", yaxs="i",
xlab="Nonexceedance probability", lwd=3, col="grey50",
ylab=paste0("P", hyphen, "value of respective test of independence"))
lines(sort( rank(W) ) / ( length(W) + 1 ), sort(W), col="darkorange1")
txt <- c(paste0("Independence test (copula::indepTest) has rejection rate = ",
round( sum(I < ALPHA) / nsim, digits=digits), ""),
paste0("Schweizer-Wolff test (copBasic::wolfCOPtest) has rejection rate = ",
round( sum(W < ALPHA) / nsim, digits=digits), ""))
legend("topleft", txt, bty="o", box.lty=0, bg=NA, inset=0.02, cex=0.9,
lwd=c(3,1), col=c("grey50", "darkorange1"))
mtext(paste0("Empirical distribution of p", hyphen, "values"), line=0.5, font=2) #
## End(Not run)
## Not run:
# DATA CENSORING EXAMPLE ------------------------------------------------------------
set.seed(10) # Example using the zmat feature for support of censoring by simulation.
muslp <- 0.005; sdslp <- 0.019 # drift with exp used to keep parameter positive
yy <- 1980:2021; n <- length(yy) # Simulated log-Pearson III annual peak streamflows
qq <- 10^sapply(yy, function(y) { # put a small drift in the log10 mean and variation
para <- lmomco::vec2par(c(3.5+(y-min(yy))*muslp,
exp(log(0.5)-(y-min(yy))*sdslp), 0.3), type="pe3")
lmomco::rlmomco(1, para) }) # flow
zul <- zur <- rep(NA, n) # initialize the lower and upper limits for years
zvl <- zvr <- rep(NA, n) # initialize the lower and upper limits for flows
# NA years are not simulated and the data are used as is
# Make some highly contrived censoring for the 3, 6, 9, 13, and 35+th observations.
zul[ 3] <- yy[ 3] - 1.1; zur[ 3] <- yy[ 3] + 1.1
zul[13] <- yy[13] - 1.5; zur[13] <- yy[13] + 2.4
zvl[ 3] <- 10; zvr[3] <- qq[3]
zvl[ 6] <- qq[6]; zvr[6] <- max(qq) + sqrt(.Machine$double.eps) # nudge above max
zvl[ 9] <- 10^(log10(qq[9]) - 0.3); zvr[9] <- 10^(log10(qq[9]) + 0.2) # log-cycles
qq[ 9] <- sqrt( zvl[9] * zvr[9] ) # geometric mean
zvl[(n-6):n] <- 10; zvr[(n-6):n] <- qq[(n-6):n]
xy <- cbind(yy, qq); zmat <- cbind(zul, zur, zvl, zvr)
swolf <- wolfCOPtest(xy[,1], xy[,2], zmat=zmat, asuv=FALSE, statf=median, rndphi=500)
ylim <- log10( range( qq ) ); ylim <- 10^c( floor( ylim[1] ), ceiling( ylim[2] ) )
plot( yy, qq, type="n", log="y", xlab="Water year", ylab="Streamflow", ylim=ylim)
segments(x=zul, x1=zur, y=qq, y1=qq, col="darkgoldenrod3")
segments(x=yy, x1=yy, y=zvl, y1=zvr, col="darkgoldenrod3")
points(yy, qq, bg="darkgoldenrod1", col="darkgoldenrod4", pch=21)
naive <- wolfCOPtest(xy[,1], xy[,2], asuv=FALSE, ties.method="random")
txt <- paste0("Sigma simulation (nsim=",
prettyNum(length(swolf$rand_sigma), big.mark=","), ") ",
"p.value=", round(swolf$p.value[1], digits=3), "; ",
"Naive p.value=", round(naive$p.value[1], digits=3))
mtext(txt, line=1) # Visually the point values tell us to expect independence to be
# rejected, but when the data censoring is taken into account and simulation used,
# we do not reject independence that is heavily influenced by the last seven data
# values on the right being censored down through the bottom axis.
swolf$rand_sigma <- sort(swolf$rand_sigma)
p <- sapply(swolf$rand_sigma, function(k) { # distribution of p-values
n <- ifelse(swolf$sample_size <= 40, 2, 1)
wolfCOPtest(k, swolf$sample_size)$p.value[n] })
swolf$rand_p.values <- p
# plot(swolf$rand_sigma, swolf$rand_p.values, pch=16, col=grey(0, 0.05)) #
## End(Not run)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.