Zero inflated beta regression using gamlss for vegetation cover data

Question

My goal is to analyse vegetation cover data. The way the data collection works is that you throw a quadrat (0.5m x 0.5m) in a sample plot and estimate the percent cover of the target species. Here is an example:

df <- structure(list(site = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 
3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L), .Label = c("A", 
"B", "C", "D", "E"), class = "factor"), cover = c(0.0013, 0, 
0.0208, 0.0038, 0, 0, 0, 0.0043, 0, 0.002, 0.0068, 0.0213, 0, 
0.0069, 0.0075, 0, 0, 0, 0.013, 0.0803, 0.0328, 0.1742, 0, 0, 
0.0179, 0, 0.3848, 0.1875, 0, 0.2775, 0.03, 0, 0.0042, 0.0429
)), .Names = c("site", "cover"), class = "data.frame", row.names = c(NA, 
-34L))

To visualize the data:

require(ggplot2)
ggplot(df, aes(site, cover)) + geom_boxplot() + labs(x = "Site", y = "Vegetation cover") +
  scale_y_continuous(labels = scales::percent)

Since the data are percentages/proportions on the interval [0,1), I figured a zero inflated beta regression would be appropriate. I do this using the gamlss package in R:

require(gamlss)

m1 <- gamlss(cover ~ site,  family = BEZI, data = df, trace = F)
summary(m1)

#******************************************************************
#  Family:  c("BEZI", "Zero Inflated Beta") 
#
#Call:  gamlss(formula = cover ~ site, family = BEZI, data = df, trace = F) 
#
#Fitting method: RS() 
#
#------------------------------------------------------------------
#Mu link function:  logit
#Mu Coefficients:
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  -3.6698     0.4265  -8.605 3.21e-09 ***
#siteB         0.4444     0.5714   0.778 0.443510    
#siteC         1.1225     0.5834   1.924 0.064948 .  
#siteD         2.2081     0.4984   4.430 0.000141 ***
#siteE         0.3819     0.7289   0.524 0.604566    
#------------------------------------------------------------------
#Sigma link function:  log
#Sigma Coefficients:
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)   2.8155     0.3639   7.738 2.54e-08 ***
#------------------------------------------------------------------
#Nu link function:  logit 
#Nu Coefficients:
#            Estimate Std. Error t value Pr(>|t|)
#(Intercept)  -0.3567     0.3485  -1.024    0.315
#
#------------------------------------------------------------------
#No. of observations in the fit:  34 
#Degrees of Freedom for the fit:  7
#      Residual Deg. of Freedom:  27 
#                      at cycle:  12 
#
#Global Deviance:     -43.22582 
#            AIC:     -29.22582 
#            SBC:     -18.54129 
#******************************************************************

So far so good. Now on to the estimates and their standard errors:

means_m1 <- lpred(m1, type='response', what='mu', se.fit=T)

df_fit <- data.frame(SITE = df$site, M = means_m1$fit, SE = means_m1$se.fit)

ggplot(df_fit, aes(SITE, M)) + geom_pointrange(aes(ymin=M-SE, ymax=M+SE)) + 
  labs(x="Site",y="Vegetation cover") + scale_y_continuous(labels=scales::percent)

My last point, and that's where I get stuck, is to figure out whether there are significant differences among those estimates.

My questions are (1) does this gamlss approach seems sound? and (2) is there a way to do mean comparisons for those models (i.e. posthoc tests). Or is this simply not possible with gamlss models?

Answer 1

I have added preliminary support for gamlss to the emmeans package...

> emmeans(m1, "site", type = "response")
 site   response         SE df  asymp.LCL  asymp.UCL
 A    0.02484787 0.01033381 NA 0.01092510 0.05551769
 B    0.03821898 0.01775760 NA 0.01518276 0.09290972
 C    0.07260824 0.03116918 NA 0.03063369 0.16245783
 D    0.18820172 0.04509320 NA 0.11504471 0.29250294
 E    0.03598930 0.02304948 NA 0.01005072 0.12070702

Confidence level used: 0.95 
Intervals are back-transformed from the logit scale 

> pairs(.Last.value)
 contrast odds.ratio         SE df z.ratio p.value
 A / B     0.6412302 0.36639164 NA  -0.778  0.9371
 A / C     0.3254574 0.18987682 NA  -1.924  0.3044
 A / D     0.1099111 0.05478393 NA  -4.430  0.0001
 A / E     0.6825356 0.49751299 NA  -0.524  0.9850
 B / C     0.5075516 0.32214846 NA  -1.068  0.8228
 B / D     0.1714066 0.09450455 NA  -3.199  0.0120
 B / E     1.0644159 0.82513475 NA   0.081  1.0000
 C / D     0.3377125 0.18217951 NA  -2.012  0.2599
 C / E     2.0971579 1.63694350 NA   0.949  0.8777
 D / E     6.2098905 4.44082214 NA   2.554  0.0792

P value adjustment: tukey method for comparing a family of 5 estimates 
Tests are performed on the log odds ratio scale 

A caution: The estimates and comparisons apply to the non-zero-inflated part of the model. You used a common zero inflation for all sites, so that doesn't mess-up the comparisons much, though the estimated means are too large because the model estimates that about 40% are zero.

This does not (yet, if ever) support models that include smoothing components It is available on the emmeans github repository soon. No guarantee it is bug-free.