2
$\begingroup$

Consider performing a statistical hypothesis test and wants to achieve significance $\alpha$. What is a good rule of thumb for an allowable deviation from the desired confidence level?

Say I want to perform a hypothesis test with $\alpha=0.05$, but given the actual nature of the underlying processes and test applied, the true probability of type I error is $0.01$ or $0.10$, is this generally acceptable? Surely there are at least a few widely agreed upon standards for what decision from desired $\alpha$ is acceptable.

As an example, I've been investigating one-way ANOVA, to test the interaction of unequal variance, non-normality, and unequal sample sizes on the actual probability of rejecting the hypothesis of equal means. As long as the population variances are within a factor of even 10 of each other (for normal, factor of 4 seems good even for exponential) then it does seem to give significance close to the desired level. I also find that the populations with large variance having the same or larger sample size (but not like an order of magnitude larger) as compared to the populations with smaller variance, then it is also fairly robust still. The true probability of type I error is usually from $1\%$ to $10\%$ when actually shooting for a $5\%$ level, if these factors are controlled reasonably.

My gut tells me that a deviation of $\pm 3\%$ or so from a desired level of $5\%$ would be acceptable.

So the question is: what deviation from $\alpha$ is generally acceptable for hypothesis tests?

Here is my (probably very inefficient) R code for investigating 1-way anova with deviations from normality and equal variances:

NS=10000
sstrd=vector("numeric",length=NS)
ssed=sstrd
sstd=sstrd
nv=round(c(5+runif(1)*10,
           5+runif(1)*10,
           5+runif(1)*10,
           5+runif(1)*10))
mv=c(0,0,0,0)
sv=c(1,1+runif(1)*1.5,1+runif(1)*2,1+runif(1)*3)
for (ns in 1:NS){
  x1=rnorm(nv[1],mv[1],sv[1]) #rexp(nv[1],1/sv[1])-sv[1] #
  x2=rnorm(nv[2],mv[2],sv[2]) #rexp(nv[2],1/sv[2])-sv[2]
  x3=rnorm(nv[3],mv[3],sv[3]) #rexp(nv[3],1/sv[3])-sv[3]
  x4=rnorm(nv[4],mv[4],sv[4]) #rexp(nv[4],1/sv[4])-sv[4]
  data=c(x1,x2,x3,x4)
  group=c(rep("A1",length(x1)),rep("A2",length(x2)),rep("A3",length(x3)),rep("A4",length(x4)))
  group=as.factor(group)
  xdd=mean(data)
  SST=sum((data-xdd)^2)
  g=levels(group)
  k=length(g)
  n=length(data)
  nvec=as.numeric(table(group))
  SSTr=0
  xd=vector("numeric",length=length(g))
  for (i in 1:length(g)){
    xd[i]=mean(data[group==g[i]])
    SSTr=SSTr+nvec[i]*(xd[i]-xdd)^2
  }
  xdvec=NULL
  for (i in 1:length(g)){
    xdvec=c(xdvec,rep(xd[i],nvec[i]))
  }
  SSE=sum((data-xdvec)^2)
  pv=nv/sum(nv)
  spvm=c(sum(pv[c(2,3,4)]),sum(pv[c(1,3,4)]),sum(pv[c(1,2,4)]),sum(pv[c(1,2,3)]))
  sta=sum(sv^2*((1-pv)^2+pv*spvm))
  sstrd[ns]=SSTr
  sstd[ns]=SST
  ssed[ns]=SSE
}
sta
c(mean(sstd),sd(sstd))
c(mean(sstrd),sd(sstrd))
c(mean(ssed),sd(ssed))
f=(sstrd/3)/(ssed/(sum(nv)-4))
plot(ecdf(f))
xv=seq(0,max(f),length=1000)
lines(xv,pf(xv,3,sum(nv)-4),col="red")
1-sum(f<=qf(0.95,3,sum(nv)-4))/length(f) # actual 1-alpha

Assuming a desired level of $\alpha=5\%$, the last line of code is where I calculate the actual probability of going beyond the threshold $f$ value.

I'm not looking for corrections to my R code, per se, but feel free to point out any errors if there are any.

$\endgroup$
  • $\begingroup$ Good idea to google 'Welch ANOVA simulation rank transformation, Kruskal'. Go to page 3 to get beyond commercial stuff with pay walls. Some people have compared Welch ANOVA (unequal variances), with standard ANOVA, and other alternatives. In the process they have found some examples such as you seek. $\endgroup$ – BruceET Jun 29 '19 at 23:36
  • 1
    $\begingroup$ You are accumulating votes to close. Your question about 'generally considered kosher' is a bit too vague. Some may think you're asking us to check your code; off-topic for this site. Maybe you can edit your question, beginning by briefly stating your purpose and method of approach? Then succinctly ask your central question. After that giving code will seem more relevant. Maybe skip CIs since testing seems your main focus. $\endgroup$ – BruceET Jun 30 '19 at 0:28
  • 1
    $\begingroup$ I'll edit the question. I can see why it may seem somewhat unclear. Mostly I'm interested in standards on what decision from significance or confidence is widely thought of as acceptable. I am using anova as a test case/example. $\endgroup$ – jdods Jun 30 '19 at 0:42
  • $\begingroup$ Since coverage is basically a binomial variable (fraction of iid outcomes of simulations that satisfy a criterion) I would think of this in terms of standard deviations of the nominal coverage (sqrt(n*p*(1-p))). $\endgroup$ – Ben Bolker Jun 30 '19 at 1:15
  • 3
    $\begingroup$ What is "acceptable" depends on who is doing the accepting. There's no universal standard of what people agree is okay. $\endgroup$ – Glen_b -Reinstate Monica Jun 30 '19 at 3:13
2
$\begingroup$

I think this kind of exploration of the robustness of one-factor ANOVA is very worthwhile. There are several reasons you might not have found interesting results. Let me explore in the direction of my top two best guesses: (a) You might not be looking at the deviations from assumption that cause real problems. (b) Your code may not be doing exactly what you expect.

Regarding (a): It is well known that the pooled two-sample t test fails when sample sizes differ and the smaller sample has the greater population variance. In particular, the true level of significance may be much larger than the intended 5% level. The Welch t test handles this situation much better.

Because the traditional one-factor ANOVA is a generalization of the pooled t test to accommodate $k > 2$ levels of the factor, it is reasonable to expect that the ANOVA will behave badly for an unbalanced design in which the smallest number of replications matches the level with the greatest population variance. So let's start by looking there.

Regarding (b): At least at the start, it is best to use built-in functions to use pre-debugged code. Then when everything is working as it should, maybe try to substitute your own code to speed up the simulation.

For example, here is R code for getting the P-value of a particular ANOVA.

# Generate random data to specifications
set.seed(2019)
x1 = rnorm(5, 100, 10)  # note small nr of reps
x2 = rnorm(10, 100, 10)
x3 = rnorm(10, 100, 10)
x = c(x1,x2,x3)
g = as.factor(rep(1:3, c(5,10,10)))
# 'Steal' P-value from built-in function
anova(lm(x ~ g))
Analysis of Variance Table

Response: x
          Df  Sum Sq Mean Sq F value Pr(>F)
g          2   37.53  18.767  0.1547 0.8576  # <-- P-value here
Residuals 22 2668.23 121.283  
anova(lm(x ~ g))[1,5]                        # Pick it off
[1] 0.8575656

All assumptions of ANOVA are met, $H_0$ is true, and P-value exceeds 5%, as expected (most of the time).

Assumptions met: Now let's replicate the above m = 10^4 times and check that the true significance level is 5%.

set.seed(2019)
m = 10^4;  pv=numeric(m)   # wrap above code in a for-loop
for(i in 1:m) {
  x1=rnorm(5,100,10); x2=rnorm(10,100,10); x3=rnorm(10,100,10)
  x = c(x1,x2,x3); g = as.factor(rep(1:3, c(5,10,10)))
  pv[i]=anova(lm(x~g))[1,5] }
mean(pv<=.05)
[1] 0.0512
2*sd(pv<=.05)/sqrt(m)   # aprx 95% margin of simulation error
[1] 0.004408329

The rejection rate is close to 5% $(0.05 \pm 0.004).$ Using m=10^6 iterations would give almost exactly 5%, but that would require a faster computer or more patience than I have right now.

Also, the distribution of the P-values is very nearly uniform, as one expects when $H_0$ is true and assumptions are met.

enter image description here

Unbalanced, heteroscedastic: Finally, lets try looking at the result when variances are unequal, as described earlier with $n_1 = 5; \sigma_1 = 25:$

set.seed(629)
m = 10^4;  pv=numeric(m)
for(i in 1:m) {
 x1 = rnorm(5,100,25); x2 = rnorm(10,100,10); x3 = rnorm(10,100,10)
 x = c(x1,x2,x3); g = as.factor(rep(1:3, c(5,10,10)))
 pv[i]=anova(lm(x~g))[1,5] }
mean(pv<=.05)
[1] 0.1602
2*sd(pv<=.05)/sqrt(m)
[1] 0.007336196

Now, the true significance level is about 16% $(0.160 \pm 0.007).$

For this simulation, the distribution of the P-values is clearly not uniformly distributed. Because $H_0$ is true, this is an indication that assumptions are not met.

enter image description here

Notes: The method of 'stealing' P-values from built-in functions runs relatively slowly. I guess the main reason is that R takes the trouble to format the entire ANOVA table, from which I pick off only the P-value.

If you're doing a lot of these simulations, it's worthwhile writing your own code. But make sure it gives the same P-value as the built-in R procedure. Starting with the same seed, you should get the same answer.

You can check power against various alternatives using essentially the same code. Perhaps compare results with true and false $H_0$'s using anova(lm) and the Welch version of the one-factor ANOVA using oneway.test.

If you're exploring effects of non-normal data I think it may be fruitful to explore exponential data with small numbers of replications. (Poisson data is also heteroscedastic, if $H_0$ is false, but I have not found interesting difficulties with Poisson data. To the point of wondering if the square-root transformation is worthwhile in a one-factor ANOVA.)

$\endgroup$
  • $\begingroup$ Ok, this at least confirms some of my suspicions, but also makes me think I might be doing something wrong too. Would you agree that 16% instead of the target 5% is pretty bad? What I'd it was 7% instead of the target 5%? I feel like that would be really good. I've done similar calculations with the Welch t test and found similar results. $\endgroup$ – jdods Jun 30 '19 at 0:56
  • $\begingroup$ If the same parameters and sample sizes led to my 16% and your 7%, then obviously one of us is wrong. // Welch's t should give close to 5% across the board when $H_0$ is true and $\alpha = 0.05.$ $\endgroup$ – BruceET Jun 30 '19 at 2:36
  • $\begingroup$ I mean if for a given population, when performing a test and wanting a 5% level, if instead the deviation from assumptions caused your true level to be anywhere in 2% to 8% say, would you agree that's good enough? Surely at some point it is too high or to low. 16% seems like the violation of test assumptions is too severe. At what point dues it become to small? 1%, 0.5%? (When going for 5%) $\endgroup$ – jdods Jun 30 '19 at 3:02
  • 1
    $\begingroup$ My code does match yours when I use your parameters. $\endgroup$ – jdods Jun 30 '19 at 4:29
  • $\begingroup$ OK, then. Seems you're on the right track. $\endgroup$ – BruceET Jun 30 '19 at 5:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.