Power analysis

One-sample

Here the Canadian temperature dataset from the previous section will be used to demonstrate one-sample power analysis.

A full script for one-sample power analysis appears below, and immediately below that we shall consider it in more detail.

import power1d

J        = 8    # sample size
Q        = 365  # continuum size

# construct baseline geometry:
g0       = power1d.geom.GaussianPulse( Q , q=200 , fwhm=190 , amp=40 )
g1       = power1d.geom.Constant( Q , amp=23 )
baseline = g0 - g1  # subtract the geometries

# construct signal geometry:
signal0  = power1d.geom.Null( Q )
signal1  = power1d.geom.GaussianPulse( Q , q=200 , fwhm=100 , amp=5 )

# construct noise model:
noise0   = power1d.noise.Gaussian( J , Q , mu = 0 , sigma = 0.3 )
noise1   = power1d.noise.SmoothGaussian( J , Q , mu = 0 , sigma = 3 , fwhm = 70 )
noise    = power1d.noise.Additive( noise0 , noise1 )

# create data sample models:
model0   = power1d.models.DataSample( baseline, signal0, noise, J=J )
model1   = power1d.models.DataSample( baseline, signal1, noise, J=J )

# create experiment models:
teststat = power1d.stats.t_1sample
emodel0  = power1d.models.Experiment( model0 , teststat )    # null hypothesis
emodel1  = power1d.models.Experiment( model1 , teststat )    # alternative hypothesis

# simulate the experiments:
sim      = power1d.ExperimentSimulator( emodel0 , emodel1 )
results  = sim.simulate( 5000, progress_bar=True )

# visualize the power results:
results.plot()

(Source code, png, hires.png, pdf)

All of the commands from the code above up until “# create experiment models” are describe in the previous section regarding data modeling. The only difference is that here two signals and two data sample models are created: one representing the null hypothesis and one representing the alternative hypothesis. In power1d the “Null” geometry is a null continuum (i.e. a continuum with all zero values) as seen here.

Once the signals and data samples are constructed, the next step is to add them to a power1d.models.Experiment object along with a test statistic function. When there is just one data sample model the test statistic function should accept a (J x Q) data array and return a (Q,) array representing the test statistic continuum. Test statistic functions for a variety of common experimental designs are available in power1d.stats.

Next the two Experiment models, representing the null and alternative hypotheses are added to an ExperimentSimulator object. When the simulate method is called for N iterations (here N$=5000) the simulator will simualte both the null and the alternative experiments for *N iterations and store the test statistic continua for both as (N x Q) arrays in the results.Z0 and results.Z1 attributes, respectively.

Last, the simulation results are summarized graphically using the plot method.

The top left panel of the results depicts the null and alternative models and lists their omnibus powers. By definition the ‘power’ of the null model is alpha, and alpha is 0.05 by default. Alpha and all other power-relevant parameters can be changed using the the results.set_ methods as summarized in the power1d.results API.

Here the alternative model has an omnibus power of 0.843. This sugges that, if our noise models are correct, we will reject the null hypothesis (of null signal) with a probability of 0.843. However, this doesn’t imply that we will detect the modeled signal with a probability of 0.843 because occassionally we will detect ‘signal’ which randomly appears in continuum regions outside of the signal region. This point is emphasized in the point-of-interest (POI) and center-of-interest (COI) power continua in the lower right panel. The POI continuum indicates the probability of rejecting the null hypothesis at particular continuum points, and the COI continuum indicates the same but extends the search region to a certain radius (here 3 points) around the POI. The POI and COI continuum results suggest that we are only about 70% likely to detect the signal we have modeled, even though the omnibus power is greater than 80%.

If we were to broaden the COI radius to the full continuum size of $Q$=365 using results.set_coi_radius(365), the COI continuum would plateau at the omnbus power of 0.843. Thus the POI and COI continua quantify the regional specificity of the omnibus power results.

Two-sample

A slightly simpler model than the one above will be used to demonstrate power calculations for two-sample and other experiments.

In the example below signal is modeled as a GaussianPulse at point 65 in a continuum of length 101 and noise is modeled as uncorreletated Gaussian. The two-sample test statistic function (power1d.stats.t_2sample_fn) is used to calculate the test statistic continuum

Note that power1d.stats functions with “_fn” appended to the end are slightly more efficient calculators than are those without “_fn” appended. This is because the “_fn” functions store calculations results like matrix inversion results so that they don’t have to be repeated across every single simulation iteration. More details regarding this are provided in power1d.stats.

Also note that, unlike the one-sample results above, power1d separates the null and alternative models into separate panels. This is done to prevent single-panel clutter.

import numpy as np
import power1d

#(0) Create geometry and noise models:
JA,JB,Q   = 5, 7, 101
baseline  = power1d.geom.Null(Q=Q)
signal0   = power1d.geom.Null(Q=Q)
signal1   = power1d.geom.GaussianPulse(Q=101, q=65, amp=2.5, sigma=10)
noise     = power1d.noise.Gaussian(J=5, Q=101, sigma=1)

#(1) Create data sample models:
modelA0   = power1d.models.DataSample(baseline, signal0, noise, J=JA) #null A
modelB0   = power1d.models.DataSample(baseline, signal0, noise, J=JB) #null N
modelA1   = power1d.models.DataSample(baseline, signal0, noise, J=JA) #alternative A
modelB1   = power1d.models.DataSample(baseline, signal1, noise, J=JB) #alternative B

#(2) Create experiment models:
teststat  = power1d.stats.t_2sample_fn(JA, JB)
expmodel0 = power1d.models.Experiment([modelA0, modelB0], teststat) #null
expmodel1 = power1d.models.Experiment([modelA1, modelB1], teststat) #alternative

#(3) Simulate experiments:
np.random.seed(0)
sim       = power1d.ExperimentSimulator(expmodel0, expmodel1)
results   = sim.simulate(1000, progress_bar=True)
results.plot()

(Source code, png, hires.png, pdf)

Regression

Regression can be implemented using the “regressor” keyword argument when constructing a DataSample model as indicated below.

Note that regression models are plotted with a different color scheme from those above: dark and bright colors correspond to small and large regressor values, respectively. In the “Alternative” model below (top right panel) it is apparent that dependent variable values get larger in the vicinity of Q=65 as the regressor values get larger, implying positive correlation between the independent and dependent variables in this vicinity.

import numpy as np
import power1d

#(0) Create geometry and noise models:
J,Q       = 30, 101
x         = np.linspace(0, 2, J)  #regressor (must have J values)
baseline  = power1d.geom.Null(Q=Q)
signal0   = power1d.geom.Null(Q=Q)
signal1   = power1d.geom.GaussianPulse(Q=101, q=65, amp=2.0, sigma=10)
noise     = power1d.noise.Gaussian(J=5, Q=101, sigma=1)

#(1) Create data sample models:
model0    = power1d.models.DataSample(baseline, signal0, noise, J=J, regressor=x)
model1    = power1d.models.DataSample(baseline, signal1, noise, J=J, regressor=x)

#(2) Create experiment models:
teststat  = power1d.stats.t_regress_fn(x)
expmodel0 = power1d.models.Experiment(model0, teststat)
expmodel1 = power1d.models.Experiment(model1, teststat)

#(3) Simulate experiments:
np.random.seed(0)
sim       = power1d.ExperimentSimulator(expmodel0, expmodel1)
results   = sim.simulate(100, progress_bar=True)
results.plot()

(Source code, png, hires.png, pdf)

One-way ANOVA

Power analyses for one-way ANOVA can be conducted almost identically to the two-sample analysis described above. The only difference is that more DataSample models are needed to represent each of the groups.

import numpy as np
import power1d

#(0) Create geometry and noise models:
JA,JB,JC,Q = 5, 7, 12, 101
baseline   = power1d.geom.Null(Q=Q)
signal0    = power1d.geom.Null(Q=Q)
signal1    = power1d.geom.GaussianPulse(Q=101, q=65, amp=1.5, sigma=10)
noise      = power1d.noise.Gaussian(J=5, Q=101, sigma=1)

#(1) Create data sample models:
modelA0   = power1d.models.DataSample(baseline, signal0, noise, J=JA)  #null A
modelB0   = power1d.models.DataSample(baseline, signal0, noise, J=JB)  #null B
modelC0   = power1d.models.DataSample(baseline, signal0, noise, J=JC)  #null C
modelA1   = power1d.models.DataSample(baseline, signal0, noise, J=JA)  #alternative A
modelB1   = power1d.models.DataSample(baseline, signal0, noise, J=JB)  #alternative B
modelC1   = power1d.models.DataSample(baseline, signal1, noise, J=JC)  #alternative C

#(2) Create experiment models:
teststat  = power1d.stats.f_anova1_fn(JA, JB, JC)
expmodel0 = power1d.models.Experiment([modelA0, modelB0, modelC0], teststat)
expmodel1 = power1d.models.Experiment([modelA1, modelB1, modelC1], teststat)

#(3) Simulate experiments:
np.random.seed(0)
sim       = power1d.ExperimentSimulator(expmodel0, expmodel1)
results   = sim.simulate(500, progress_bar=True)
results.plot()

(Source code, png, hires.png, pdf)