# Post-hoc test to determine difference in variance

I'm testing whether 5 different groups have different variances, I am not interested in the mean. I performed a Brown-Forsythe test given that the data does not fit a normal distribution. The test yielded a result of heterogeneity with a p. Value of 0.0001. My question is, how do I test where the difference in variances is? meaning, between which groups is the difference significant? I'm looking for the equivalent of a Post-Hoc Tukey test but for testing variances. Thank you, this has been driving me crazy.

• Given that the Brown-Forsythe test is really an ANOVA of the absolute residuals (from the group medians), have you considered simply applying Tukey's HSD to those absolute residuals? What would be the objections to such an approach?
– whuber
Mar 29, 2018 at 13:34

I worked on a project not too long ago where I approached this type of problem within the framework of generalized least squares regression. In other words, I fitted a model to my data which estimated the means of the groups while allowing the standard deviations of groups 2, 3, ..., 4 to be equal to:

• the standard deviation of group 1 times a factor delta2 (for group 2);
• the standard deviation of group 1 times a factor delta3 (for group 3);
• the standard deviation of group 1 times a factor delta4 (for group 4);
• the standard deviation of group 1 times a factor delta5 (for group 5).

If I recall correctly, the software I used produced estimates and confidence intervals for each of these factors and that enabled me to get all the pairwise comparisons of interest between groups in terms of standard deviations. I had to transform the outcome variable to get well-behaved residuals and also fight against the software because it didn't let me choose what group would be treated as a reference (e.g., group 1), but would instead choose the reference group based on the data.

The R syntax for fitting this type of model would be something like:

library(nlme)
model <- gls(outcome ~ group, weights = varIdent(form = ~1|group),
data=mydata)


Then use something like:

summary(model)


to see a summary of the model fit.

Furthermore, use:

model$modelStruct$varStruct


to get estimates of delta1 (which will be 1 for the reference group 1), delta2, delta3, delta4 and delta5.

Finally, use:

intervals(model, which="var-cov")


to get 95% confidence intervals for delta2, delta3, delta4, delta5 and for the standard deviation sigma of the reference group (i.e., group 1). The latter will be listed under the heading Residual standard error.

See pages 159-162 of the book Linear Mixed-Effects Models Using R: A Step-by-Step Approach for an example. The book was written by Galecki and Burzykowski. Section 7.6.2 of the book gives the formula for a confidence interval for the logarithm of deltas, where s can be 2, 3, 4 or 5.

Maybe other people here will give you other ideas for how to proceed, but I thought I would share my idea in case it might spark further conversation.