Power Estimation
Discrete Data/Model
Simple Null Hypothesis
vals = 0:10
pnull = function() pbinom(0:10, 10, 0.5)
rnull =function () table(c(0:10, rbinom(100, 10, 0.5)))-1
ralt =function (p) table(c(0:10, rbinom(100, 10, p)))-1
P=gof_power_disc(pnull, rnull, vals, ralt,
param_alt=seq(0.5, 0.6, 0.01), B=c(500, 500), nbins=c(11, 5))
plot_power(P, "p", Smooth=FALSE)
In all cases the arguments vals, pnull and rnull are the same as for the testing routines. In addition we now have
ralt: a routine with one parameter that generates data under some alternative hypothesis.
param_alt: values to be passed to ralt. This allows the calculation of the power for many different values.
B=c(1000, 1000) the first number is the number of simulation runs for power estimation and the second the number of runs to be used to find the null distribution.
Composite Null Hypothesis
vals = 0:10
pnull = function(p) pbinom(0:10, 10, p)
rnull = function (p) table(c(0:10, rbinom(100, 10, p)))-1
phat = function(x) sum(0:10*x)/1000- Null Hypothesis is true
ralt =function (p) table(c(0:10, rbinom(100, 10, p)))-1
gof_power_disc(pnull, rnull, vals, ralt, c(0.5, 0.6), phat,
B=c(100, 200), nbins=c(11, 5), maxProcessors = 2)
#> KS K AD CvM W ZA ZC Wassp1 chi large Pearson
#> 0.5 0.06 0.06 0.04 0.08 0.05 0.05 0.04 0.06 0.02
#> 0.6 0.05 0.06 0.07 0.07 0.03 0.05 0.00 0.06 0.03
#> chi small Pearson chi large LR-m chi small LR-m
#> 0.5 0.06 0.02 0.06
#> 0.6 0.02 0.03 0.02Note that power estimation in the case of a composite hypothesis (aka with parameters estimated) is much slower than the simple hypothesis case.
- Null Hypothesis is false
ralt =function (p) table(c(rep(0:10, 2), rbinom(100, 10, p)))
gof_power_disc(pnull, rnull, vals, ralt, 0.5, phat,
B=c(100, 200), nbins=c(11, 5), maxProcessors = 2)
#> KS K AD CvM
#> 1.00 1.00 1.00 1.00
#> W ZA ZC Wassp1
#> 0.63 1.00 0.00 1.00
#> chi large Pearson chi small Pearson chi large LR-m chi small LR-m
#> 0.23 0.09 0.17 0.10Continuous Data/Model
Simple Null Hypothesis
pnull = function(x) pnorm(x)
qnull = function(x) qnorm(x)
rnull = function() rnorm(100)
ralt = function(mu=0) rnorm(100, mu)
gof_power_cont(pnull, rnull, qnull, ralt, c(0, 1), B=c(100, 200))
#> KS K AD CvM W ZA ZK ZC Wassp1 EP large Pearson
#> 0 0.02 0.02 0.04 0.04 0.04 0.07 0.06 0.05 0.04 0.07
#> 1 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
#> ES large Pearson EP small Pearson ES small Pearson EP large LR-m
#> 0 0.12 0.12 0.05 0.03
#> 1 1.00 1.00 1.00 1.00
#> ES large LR-m EP small LR-m ES small LR-m
#> 0 0.06 0.06 0.05
#> 1 1.00 1.00 1.00Composite Null Hypothesis
pnull = function(x, p=c(0,1)) pnorm(x, p[1], ifelse(p[2]>0, p[2], 0.01))
qnull = function(x, p=c(0,1)) qnorm(x, p[1], ifelse(p[2]>0, p[2], 0.01))
rnull = function(p=c(0,1)) rnorm(500, p[1], p[2])
ralt = function(mu=0) rnorm(100, mu)
phat = function(x) c(mean(x), sd(x))
gof_power_cont(pnull, rnull, qnull, ralt, c(0, 1), phat, B=c(100, 200),
maxProcessor=2)
#> KS K AD CvM W ZA ZK ZC Wassp1 EP large Pearson
#> 0 0.94 0.90 0.08 0.09 0.09 0.64 0.01 0.01 1 0.06
#> 1 0.97 0.89 0.09 0.06 0.06 0.60 0.02 0.03 1 0.04
#> ES large Pearson EP small Pearson ES small Pearson EP large LR-m
#> 0 0.03 0.09 0.04 0.04
#> 1 0.04 0.08 0.04 0.02
#> ES large LR-m EP small LR-m ES small LR-m
#> 0 0.04 0.06 0.05
#> 1 0.03 0.06 0.04ralt = function(df=1) {
# t distribution truncated at +- 5
x=rt(1000, df)
x=x[abs(x)<5]
x[1:100]
}
gof_power_cont(pnull, rnull, qnull, ralt, c(2, 50), phat, Range=c(-5,5), B=c(100, 200), maxProcessor=2)
#> KS K AD CvM W ZA ZK ZC Wassp1 EP large Pearson
#> 2 1.00 1.00 0.63 0.61 0.62 0.96 0.27 0.17 1 0.19
#> 50 0.96 0.94 0.04 0.06 0.04 0.69 0.00 0.02 1 0.09
#> ES large Pearson EP small Pearson ES small Pearson EP large LR-m
#> 2 0.28 0.24 0.36 0.10
#> 50 0.00 0.11 0.03 0.03
#> ES large LR-m EP small LR-m ES small LR-m
#> 2 0.24 0.21 0.37
#> 50 0.01 0.07 0.03Running other tests
Its is very easy for a user to add other goodness-of-fit tests to the package. This can be done by editing the routines TS_cont and/or TS_disc, which are located in the folder inst/examples in the Rgof library folder. Or a user can write their own version of these files.
Example
Say we wish to use a test that is a variant of the Cramer-vonMises test, using the integrated absolute difference of the empirical and the theoretical distribution function:
\[\int_{-\infty}^{\infty} \vert F(x) - \hat{F}(x) \vert dF(x)\] For continuous data we have the routine
newTS_cont = function(x, Fx, param, qnull) {
Fx=sort(Fx)
n=length(x)
out = sum(abs( (2*1:n-1)/2/n-Fx ))
names(out) = "CvM alt"
out
}This routine has to have four arguments x, Fx, param and qnull, just like TS_cont. Note that the return object has to be a named vector. Also, the arguments param and qnull are likely not needed for the calculation of the new test statistic but need to be present for consistency reasons.
Then we can run this test with
pnull = function(x) punif(x)
qnull = function(x) qunif(x)
rnull = function() runif(500)
x = rnull()
Rgof::gof_test_cont(x, pnull, rnull, qnull, TS=newTS_cont)
#> $Statistics
#> CvM alt
#> 6.9863
#>
#> $p.value
#> CvM alt
#> 0.3946Say we want to find the power of this test when the true distribution is a linear:
ralt = function(slope) {
if(slope==0) y=runif(500)
else y=(slope-1+sqrt((1-slope)^2+4*slope* runif(500)))/2/slope
}gof_power_cont(pnull, rnull, qnull, ralt, TS=newTS_cont, param_alt=round(seq(0, 0.5, length=5), 3), Range=c(0,1))
#> CvM alt
#> 0 0.049
#> 0.125 0.337
#> 0.25 0.893
#> 0.375 0.998
#> 0.5 1.000for discrete data we write the routine using Rcpp:
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector newTS_disc(IntegerVector x,
NumericVector Fx,
NumericMatrix nm,
NumericVector vals) {
Rcpp::CharacterVector methods=CharacterVector::create("CvM alt");
int const nummethods=methods.size();
int k=x.size(), n, i;
NumericVector TS(nummethods), ecdf(k);
double tmp;
TS.names() = methods;
n=0;
for(i=0;i<k;++i) n = n + x[i];
ecdf(0) = double(x(0))/double(n);
for(i=1;i<k;++i) {
ecdf(i) = ecdf(i-1) + x(i)/double(n);
}
tmp = std::abs(ecdf[0]-Fx(0))*Fx(0);
for(i=1;i<k;++i)
tmp = tmp + std::abs(ecdf(i)-Fx(i))*(Fx(i)-Fx(i-1));
TS(0) = tmp;
return TS;
}Again the routine has to have four arguments x, Fx, nm and vals, just like TS_disc, and the output vector has to have names. Aagin, the arguments nm and vals are likely not needed for the calculation of the new test statistic but need to be present for consistency reasons.
Note that one drawback of writing the routine in Rcpp is that it is then not possible to used multiple processors.
vals=1:50/51
pnull = function() (1:50)/50
rnull = function() c(rmultinom(1, 500, rep(1/50,50)))
x = rnull()
gof_test_disc(x, pnull, rnull, vals, TS=newTS_disc)
#> $Statistics
#> CvM alt
#> 0.013
#>
#> $p.value
#> CvM alt
#> 0.4648and for power calculations we can run
ralt = function(slope) {
if(slope==0) p=rep(1/50, 50)
else p=diff(slope * (0:50/50)^2 + (1 - slope) * 0:50/50)
c(rmultinom(1, 500, p))
}
gof_power_disc(pnull, rnull, vals, ralt, TS=newTS_disc, param_alt=round(seq(0, 0.5, length=5), 3))
#> Parallel Programming is not possible if custom TS is written in C++. Switching to single processor
#> CvM alt
#> 0 0.055
#> 0.125 0.367
#> 0.25 0.898
#> 0.375 0.999
#> 0.5 1.000