Incoherent post hoc pairwise Tukey test using beta binomial distribution and betareg()

I am fitting a beta binomial regression to my data, which is proportion data (lifetime reproductive success ratio) ranging [0,1]. I transformed the 0s and 1s following Smithson and Verkuilen 2006 (y.(N-1)+0,5)/N) to fit the beta regression, as I have quite some zero and they are important datapoints which can not be discarded (See response variable (RAT) distribution after transformation and dotplot of datapoints across treatments below)

I fitted the model below:

mod1 <- betareg(RAT ~ treatment*exposure, data = lrs2)
hist(lrs2$RAT, breaks=20) summary(mod1) Anova(mod1, type=3) par(mfrow=c(2,2)) plot(mod1) AIC(mod1) BIC(mod1)  And here is what came out of it: > mod1 <- betareg(RAT ~ treatment*exposure, data = lrs2) > hist(lrs2$RAT, breaks=20)
> summary(mod1)

Call:
betareg(formula = RAT ~ treatment * exposure, data = lrs2)

Standardized weighted residuals 2:
Min      1Q  Median      3Q     Max
-5.8206 -0.3005  0.2031  0.6408  6.3900

Coefficients (mean model with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept)             -0.8153     0.1228  -6.642  3.1e-11 ***
treatmentpc             -0.3720     0.1774  -2.097  0.03602 *
exposure1d               0.4516     0.1773   2.546  0.01088 *
exposure3d               0.4482     0.1711   2.619  0.00881 **
exposure7d               0.1035     0.1753   0.591  0.55481
treatmentpc:exposure1d   0.4941     0.2489   1.986  0.04707 *
treatmentpc:exposure3d   0.4812     0.2435   1.977  0.04809 *
treatmentpc:exposure7d   0.8268     0.2483   3.330  0.00087 ***

Phi coefficients (precision model with identity link):
Estimate Std. Error z value Pr(>|z|)
(phi)   4.2043     0.2667   15.77   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Type of estimator: ML (maximum likelihood)
Log-likelihood: 97.24 on 9 Df
Pseudo R-squared: 0.08414
Number of iterations: 32 (BFGS) + 2 (Fisher scoring)
> Anova(mod1, type=3)
Analysis of Deviance Table (Type III tests)

Response: RAT
Df   Chisq Pr(>Chisq)
(Intercept)         1 44.1108  3.103e-11 ***
treatment           1  4.3964    0.03602 *
exposure            3 10.7253    0.01331 *
treatment:exposure  3 11.2268    0.01056 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> par(mfrow=c(2,2))
> plot(mod1)
> AIC(mod1)
[1] -176.4816
> BIC(mod1)
[1] -140.1622


What I am most interested about is whether effect of treatment (pc/npc) is consistant or not across exposure times (1d, 3d, 7d, 14d).

Anova says it is quite so, with P=0.01, and this fits what the means + SE in red on the dotplot say, however when I run post hoc pairwise Tukey comparison, I get model estimates (obtained extracting LS.EM.MEANS of my betareg model) that really do not match what can be observed on the dotplot:

> lsmod1 <- lsmeans(mod1,
+                     pairwise ~ treatment*exposure,
>
> lsmod1
$lsmeans treatment exposure lsmean SE df asymp.LCL asymp.UCL npc 14d 0.307 0.0261 Inf 0.256 0.358 pc 14d 0.234 0.0233 Inf 0.188 0.279 npc 1d 0.410 0.0311 Inf 0.349 0.471 pc 1d 0.440 0.0291 Inf 0.383 0.497 npc 3d 0.409 0.0290 Inf 0.352 0.466 pc 3d 0.436 0.0285 Inf 0.380 0.492 npc 7d 0.329 0.0279 Inf 0.275 0.384 pc 7d 0.436 0.0293 Inf 0.379 0.494 Confidence level used: 0.95$contrasts
contrast          estimate     SE  df z.ratio p.value
npc,14d - pc,14d  0.073023 0.0347 Inf  2.104  0.4123
npc,1d - pc,1d   -0.029837 0.0426 Inf -0.701  0.9970
npc,3d - pc,3d   -0.026634 0.0406 Inf -0.656  0.9980
npc,7d - pc,7d   -0.106898 0.0404 Inf -2.645  0.1399

P value adjustment: tukey method for comparing a family of 8 estimates


> x =emmeans(mod1, spec='treatment','exposure')
> summary(x)
exposure = 14d:
treatment emmean     SE  df asymp.LCL asymp.UCL
npc        0.307 0.0261 Inf     0.256     0.358
pc         0.234 0.0233 Inf     0.188     0.279

exposure = 1d:
treatment emmean     SE  df asymp.LCL asymp.UCL
npc        0.410 0.0311 Inf     0.349     0.471
pc         0.440 0.0291 Inf     0.383     0.497

exposure = 3d:
treatment emmean     SE  df asymp.LCL asymp.UCL
npc        0.409 0.0290 Inf     0.352     0.466
pc         0.436 0.0285 Inf     0.380     0.492

exposure = 7d:
treatment emmean     SE  df asymp.LCL asymp.UCL
npc        0.329 0.0279 Inf     0.275     0.384
pc         0.436 0.0293 Inf     0.379     0.494

Confidence level used: 0.95
> contrast(x, method='pairwise')
exposure = 14d:
contrast estimate     SE  df z.ratio p.value
npc - pc   0.0730 0.0347 Inf  2.104  0.0354

exposure = 1d:
contrast estimate     SE  df z.ratio p.value
npc - pc  -0.0298 0.0426 Inf -0.701  0.4832

exposure = 3d:
contrast estimate     SE  df z.ratio p.value
npc - pc  -0.0266 0.0406 Inf -0.656  0.5120

exposure = 7d:
contrast estimate     SE  df z.ratio p.value
npc - pc  -0.1069 0.0404 Inf -2.645  0.0082


If any within-exposure pairwise comparison were to be significant, I would expect it to be within 1d and within 14d, but here LS/EM means are telling me otherwise (and I clearly don't trust what they are telling me, it's nonsense)... Does it mean my beta regression is not ideal? Here comes the R plots of the betareg.

Thanks a lot for any help/suggestion/correction you could provide.

Quentin

EDIT:

RAW MEANS

MODEL ESTIMATES

• I seem to recall that car::Anova isn't reliable with betareg. – Sal Mangiafico Dec 13 '19 at 19:35

Note that in both sets of comparisons, the estimates and standard errors for the comparisons within the same exposure match up. It is only the $$P$$ values that differ. The reason is that the $$P$$ values in the first set are adjusted for multiplicity, and the $$P$$ values in the second set are not.

The first set of EMM comparisons consists of $${8\choose2}=28$$ pairwise comparisons, of which only 4 are shown in the posting. The Tukey correction adjusts the $$P$$ values for simultaneously testing these 28 comparisons, making them considerably higher than the unadjusted $$P$$ values.

In the second set of comparisons, we treat the means as 4 families of 2 comparisons, and hence no multiplicity adjustment. The $$P$$ values shown are thus the unadjusted $$P$$ values, which are those you would have obtained in the first set if you had specified adjust = "none".

I wonder if you noticed that the plot has the exposures in numerical order: 1d, 3d, 7d, 14d -- whereas the factor levels are alphabetized by default: 14d, 1d, 3d, 7d -- is that's the source of the confusion? It's hard to see the intervals in the plot, and I wonder how they were computed. If computed using some other method than beta regression, then that implies a different model, so comparing those intervals with the intervals obtained using beta regression is like comparing apples and oranges. Remember, in this kind of data, not all observations have the same variance; those near .5 are more variable than those approaching 0 or 1.

Anyway, I disagree about these results being "nonsense."

• Thank you for your answer and clarifying why P values are different. I am not sure which makes more sense in my case though; I would say without multicomparison adjustment as other comparisons are irrelevant to the ecological hypothesis I am interested in, right? The plot gives the factors in a different order because I re-ordered them, for the reader's convenience. The means and SE on the dotplot were calculated manually, and added the the dotlot afterwards. Shouldn't the model give me estimates that are somehow ressembling what the raw data says, if the model is a good fit? – Quentin Corbel Aug 24 '19 at 8:01
• I think the 4 focused comparisons are the way to go. The estimates DO resemble your manual calculations, don’t they? But they won’t be the same. I think for reporting purposes you should draw the dotplots and superpose the EMMs and their SEs. – rvl Aug 24 '19 at 12:27
• I edited my original post and added the plots of raw means +-SE and model estimates +- SE. I am not sure how defendable the post hoc tests on model estimates will be, in regard to what the raw means +- SE show, it's a clearly different story. I am not even sure what conclusion to draw anymore. – Quentin Corbel Aug 26 '19 at 8:17
• These do look more different than one might expect. In the new plots, I think you mean to label the last group "14d", not "15d". Do you agree that the biggest discrepancies are associated with the situations having the most zero and 1 observations? That's the pattern I think I see. You say you transformed those 0s and 1s? In which computations? I'm thinking you probably should not transform any of the data you use in betareg. If you send me the data (maintainer e-mail for emmeans package), I promise not to distribute it to anyone, but could use it to check what emmeans is doing. – rvl Aug 26 '19 at 16:18
• Try this: Set grid = expand.grid(treatment = factor("npc","pc"), exposure = factor(c("14d","1d","3d","7d"))). Then do predict(mod1, newdata = grid, type = "response") -- and I am pretty sure you will get exactly the estimates that emmeans shows. So if there is a problem, it is in the model, not in emmeans. – rvl Aug 26 '19 at 16:52