anova.manyglm fails with cor.type='shrink'

[signalbash]

## Visualise the effect of treatment on copepod abundance
tasm.cop <- mvabund(Tasmania$copepods)
#treatment <- Tasmania$treatment
# create >2 treatment variables variables
treatment = c(rep(c("Undis", "Dist"), 4),rep(c("Undis2", "Dist2"), 4))
block <- Tasmania$block
#plot(tasm.cop ~ treatment, col=as.numeric(block))

## Fitting predictive models using a negative binomial model for counts:
tasm.cop.nb <- manyglm(tasm.cop ~ treatment, family="negative.binomial")
# works fine...
anova_tasm.cop.nb = anova.manyglm(tasm.cop.nb, pairwise.comp = treatment, nBoot = 3, test="wald")

## Fitting predictive models using a negative binomial model for counts, and shrink:
tasm.cop.nb_shrink <- manyglm(tasm.cop ~ treatment, family="negative.binomial", cor.type = "shrink", test="wald")
# doesn't work
anova_pairwise_tasm.cop.nb = anova.manyglm(tasm.cop.nb_shrink, pairwise.comp = treatment, nBoot = 3, test="wald")

Time elapsed: 0 hr 0 min 0 sec
Error in dimnames(x) <- dn :
length of 'dimnames' [2] not equal to array extent
In addition: Warning messages:
1: In anova.manyglm(m, show.time = "none", keep.boot = TRUE, nBoot = anova_obj$nBoot, :
The likelihood ratio test can only be used if correlation matrix of the abundances is is assumed to be the Identity matrix. The Wald Test will be used.
2: In anova.manyglm(m, show.time = "none", keep.boot = TRUE, nBoot = anova_obj$nBoot, :
The likelihood ratio test can only be used if correlation matrix of the abundances is is assumed to be the Identity matrix. The Wald Test will be used.
3: In anova.manyglm(m, show.time = "none", keep.boot = TRUE, nBoot = anova_obj$nBoot, :
The likelihood ratio test can only be used if correlation matrix of the abundances is is assumed to be the Identity matrix. The Wald Test will be used.
4: In anova.manyglm(m, show.time = "none", keep.boot = TRUE, nBoot = anova_obj$nBoot, :
The likelihood ratio test can only be used if correlation matrix of the abundances is is assumed to be the Identity matrix. The Wald Test will be used.
5: In anova.manyglm(m, show.time = "none", keep.boot = TRUE, nBoot = anova_obj$nBoot, :
The likelihood ratio test can only be used if correlation matrix of the abundances is is assumed to be the Identity matrix. The Wald Test will be used.
6: In anova.manyglm(m, show.time = "none", keep.boot = TRUE, nBoot = anova_obj$nBoot, :
The likelihood ratio test can only be used if correlation matrix of the abundances is is assumed to be the Identity matrix. The Wald Test will be used.

I've done my best to find where it's going wrong and as far as I can tell it's in do_pairwise_comp ()

am <- anova.manyglm(m,
                          show.time = 'none',
                          keep.boot = TRUE,
                          nBoot = anova_obj$nBoot,resamp = resamp, cor.type = cor.type)

which gives a table like
Analysis of Variance Table

Model: subY ~ subWhat

Multivariate test:
Res.Df Df.diff wald Pr(>wald)
(Intercept) 5
subWhat 4 2 0.25
Arguments:
Test statistics calculated assuming correlated response via ridge regularization
P-value calculated using 3 iterations via PIT-trap resampling.

and am$table[2, 3] is NaN

whereas running it without cor.type = 'shrink' gives am$table[2, 3] a numerical value (165.5)
Analysis of Deviance Table

Model: subY ~ subWhat

Multivariate test:
Res.Df Df.diff Dev Pr(>Dev)
(Intercept) 5
subWhat 4 2 165.5 0.25
Arguments:
Test statistics calculated assuming uncorrelated response (for faster computation)
P-value calculated using 3 iterations via PIT-trap resampling.

Is there a fix for this? I'm using the latest version of mvabund from github

[dwarton]

thanks for this. The issue seems to be that the pairwise Wald test stats in this case are undefined. As a simple fix I have replaced these test stat values with zero, which is somewhat conservative, but the lesson here is that Wald stats sometimes do weird things with lots of zeros and we are safer sticking with (or ground-truthing against) LR stats where possible, as you did above.