# Reporting betareg outcome - how to compare non-nested models?

I am looking for advice how to gain and report results using beta regression for an ANCOVA-like model. My model is as follows:

Model <- betareg(plants ~ rain * cows)


The questions I want to answer is whether rain, cows, or interaction between these two factors have an influence to plant relative biomass (expressed as % in a range of (0;1)). If this was a GLM, I would use anova(Model) and report results from the likelihood ratio test using Chisq. Since ANOVA is not suitable for betareg, I read that I could use lrtest or wald test and compare it manually, like this:

 m1 <- betareg(plants~ 1 ,data = mydata)
m2 <- betareg(plants~ cows,data = mydata)
m3 <- betareg(plants~ rain ,data = mydata)
m4 <- betareg(plants~ rain:cows,data = mydata)
lrtest(m1,m2,m3,m4)


However, as I understand this is not correct, as my models are not fully nested (in particular m2 and m3). I am wondering if there is an alternative to compare these models and determine which factors could have significantly affected plant biomass?

P.S. This question is closely related to https://stackoverflow.com/questions/44183329/anova-like-object-for-betaregression however, the answer did not help in my case.

• One potential option is to use emmeans::joint_tests(Model). The emmeans package supports betareg objects, but be sure to read up on the documentation for the joint_tests function. – Sal Mangiafico Aug 16 '18 at 18:57
• Thank you for your comments! Concerning the comparisons, maybe you are right and it could help to avoid using reference model (plants~1). However I find it harder to interpret, for instance if both cows and rain had significant additive effect (but maybe I am confusing something) – Kriste M Aug 21 '18 at 14:32
• I'm deleting my previous comment. It was confused. I was trying to get to Achim Zeileis's answer, but I wrote it too quickly, with too little checking on the actual results. – Sal Mangiafico Aug 21 '18 at 17:33

The situation is completely analogous to (generalized) linear models where you could do forward or backward selection via the so-called "encompassing" model (cows + rain) in your case to the interaction model. Only the two non-nested models with one main effect cannot be compared directly with a LR or Wald test but both can be compared against the encompassing model. Or you adopt an information criterion (AIC or BIC, ...) instead. So you could do

m0   <- betareg(plants ~ 1)
mC   <- betareg(plants ~ cows)
mR   <- betareg(plants ~ rain)
mCR  <- betareg(plants ~ cows + rain)
mCxR <- betareg(plates ~ cows * rain)


Then you can test both paths:

lrtest(m0, mC, mCR, mCxR)
lrtest(m0, mR, mCR, mCxR)


and use these to do either forward or backward selection. Or you could do:

BIC(m0, mC, mR, mCR, mCxR)


or something like that.

• I guess that plants ~ cows * rain and plants ~ cows:rain is really of interest to you.
• Thank you for the answer! Would it make sense then for testing the effect of cows to run lrtest(m0,mC), for rain lrtest(m0, mR) and for the interaction lrtest(mCR, mCxR)? Just in this case I have another question. Here m0 is a reference model, I defined it as betareg(plants ~ 1) as I saw it like that in several sources. But is there a reason to choose ~1? Basically, I would assume I would like to check if mC explains my data better than a random model. – Kriste M Aug 21 '18 at 14:15
• If you want to do test-based model selection, I would recommend to either use a forward selection (adding the most significant term in each step, if any) or backward selection (dropping the most non-significant term in each step, if any). Alternatively, you can do information-criterion-based selection which is not stepwise. In models with so few terms, the strategies often (but not always) lead to the same results. As for y ~ 1: This means that all observations come from the same beta distribution without any effects of cows and/or rain. – Achim Zeileis Aug 21 '18 at 15:32