1
$\begingroup$

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)

enter image description here

dotplot QC

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,
+                     adjust="tukey")
> 
> 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.

enter image description here

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

Quentin

EDIT:

RAW MEANS
enter image description here

MODEL ESTIMATES

enter image description here

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

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."

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Quentin Corbel Aug 24 '19 at 8:01
  • $\begingroup$ 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. $\endgroup$ – rvl Aug 24 '19 at 12:27
  • $\begingroup$ 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. $\endgroup$ – Quentin Corbel Aug 26 '19 at 8:17
  • $\begingroup$ 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. $\endgroup$ – rvl Aug 26 '19 at 16:18
  • $\begingroup$ 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. $\endgroup$ – rvl Aug 26 '19 at 16:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.