Veganify dispweight

This PR suggests veganification and complete re-writing and re-structuring of `dispweight` of @EDiLD. The major changes in order of importance are:

1. If `dispweight` is supposed to be a function to transform data then it should return transformed data. The PR now returns transformed data with other original items as attributes. Alternatively, if we want to keep this as returning a list of items (including transferred data item) then the minimum is to provide a method to extract the transformed data plus `print.dispweight` for compact display of the results. No `$`-sign should be needed in normal **vegan** use. 

2. **vegan** functions should use **vegan** infrastructure, and the PR replaces original inlined multinomial simulation by identical `nullmodel(…, "c0_ind")`. This brings along a complete restructuring of the code. Instead of looping over each species, we loop over levels of factors and `simulate.nullmodel` for the whole subset of community in one sweep. Because of this, code had to re-structured in pieces usable in such a loop. The basic statistic (`dhat`) was also streamlined. The results remain identical with the original code, but in my laptop the current code also is 3x faster.

3. Because PR re-structured the code, we also want to expose arguments to users. Instead of hard-coding p-value to 0.05 and simulation to multinomial, PR makes these both user arguments (other than `"c0_ind"` model may not be sensible, but if users want to try them, we should not stop them: we assume they know what they do). Other arguments were also changed to accord other **vegan** functions: other methods use `groups` instead of `group`, and here we clearly do simulations instead of permutations, hence `nperm` changed to `nsimul` (like in `oecosimu`).

4. The function should work with only one group, when `groups` argument was not given (and then default to one group), and when there are empty factor levels.

In addition, there are some minor changes that are best interpreted from commit diffs.

Finally, the PR also suggests a new method and function `gdispweight`  or generalized dispersion weighting. Clarke's dispersion weighting is based on Pearson Chi-square with Poisson error for a model with one factor. The significance of this is assessed by multinomial simulation for each species within each factor level, omitting completely the factor levels where species is absent (reducing both the total N and residual d.f. for those cases). The `gdispweight` evaluates the very same Pearson Chi-square with (Quasi-)Poisson error but directly with `glm` using any `formula`.  So we can use several factors or continuous covariates. Instead of simulation, the significance of overdispersion (=Pearson Chi-square) is assessed directly from Chi-square distribution. The weights are similar for species that occur in every factor level when `formula` only has one factor. However, the weights diverge when species does not occur at every level, because then original `dispweight` will reduce the number of degrees of freedom used as divisor in estimating over-dispersion. The p-values will always be based on different models: `gdispweight` evaluates the overall model, but `dispweight` only studies the within-level variation by its multinomial simulation. This is an unpublished method, but the generalization is obvious.
`sapply()` always goes the wrong way. Perhaps a `do.call("rbind", lapply(split(comm, group), colMeans))` is what I was thinking of...?

The reasoning was whether or not splitting data first then applying a fast `colMeans()` to each sub-community object might scale better than repeatedly splitting single columns into `groups` and then using `mean()`.

Aa for `NA`; agree - for now fail loudly immediately if any `NA` present.