qqmath.permustats and reference distribution

[jarioksa]

qqmath function accepts argument distribution to find the expected quantiles in the Q-Q plot. The default is to plot permutations against Gaussian distribution (qnorm). However, this is a bad choice in several cases. For instance, example(permutest.betadisper) gives: One panel shows F-statistic and another shows t-statistic which have different scales and shapes. It is possible to have separate scales for panels (scales = list(relation = "free")), but the shapes are still different. It is possible to give the distribution used for the theoretical quantile (horizontal axis), but it is more tedious to do this separately for each panel (see ?panel.qqmath), and getting the scaling right is even more cumbersome (I think this would need writing a prepanel function). However, actually we would like to get something like this: Here F is plotted against F-distribution and t against t-distribution, both with correct degrees of freedom and using separate scaling. This is much better than the current plot for diagnostic purposes.

I didn't find a nice way to do this with qqmath, and the quick graph above is based on the idea xyplot(permuted ~ expected | statistic, scales = list(relation = "free"), abline = c(0,1)).

To do this in general, we should

1. Add an item defining reference distribution and its parameters (degrees of freedom etc) in the permustats method when this makes sense. The default would be to have nothing and relying on qqmath using Gaussian distribution. I think the new item should be given for betadisper, adonis and anova.cca methods.
2. Rewrite qqmath.permustats to use the reference distribution for each panel when such is given. This is the harder part to get right.