I'm attempting to follow this example for a zero-inflated beta regression on my vegetation cover data, but I get this error message after I try to run the emmeans function in R. Is it because of an over-abundance of zeros?
Each observation is independent, and is the average of 10 replicate observations per "Site" (N=47 per season). SP = Species of seagrass, CYR = Calendar year, Month: 3 = March, 4 = April, etc.., Percent = proportion of 1m^2 quadrat covered by SP, SEAS = Season (wet,dry - this is only dry season data). The goal is to see if there are differences in vegetation cover between months, by year.
(752 obs)
head(dry, 15)
YR_SEAS_MO_Site use_for_analysis CYR Month Season Site SP Percent YR_SEAS
1 2014_DRY_3_10 1 2014 3 DRY 10 Sf 0.00 2014_DRY
2 2014_DRY_3_10 1 2014 3 DRY 10 Rm 0.00 2014_DRY
3 2014_DRY_3_10 1 2014 3 DRY 10 Tt 0.12 2014_DRY
4 2014_DRY_3_10 1 2014 3 DRY 10 Hw 0.09 2014_DRY
5 2014_DRY_3_11 1 2014 3 DRY 11 Hw 0.17 2014_DRY
6 2014_DRY_3_11 1 2014 3 DRY 11 Sf 0.00 2014_DRY
7 2014_DRY_3_11 1 2014 3 DRY 11 Rm 0.00 2014_DRY
8 2014_DRY_3_11 1 2014 3 DRY 11 Tt 0.06 2014_DRY
9 2014_DRY_3_12 1 2014 3 DRY 12 Tt 0.01 2014_DRY
10 2014_DRY_3_12 1 2014 3 DRY 12 Hw 0.34 2014_DRY
11 2014_DRY_3_12 1 2014 3 DRY 12 Sf 0.00 2014_DRY
12 2014_DRY_3_12 1 2014 3 DRY 12 Rm 0.00 2014_DRY
13 2014_DRY_3_13 1 2014 3 DRY 13 Rm 0.00 2014_DRY
14 2014_DRY_3_13 1 2014 3 DRY 13 Tt 0.01 2014_DRY
15 2014_DRY_3_13 1 2014 3 DRY 13 Hw 0.53 2014_DRY
> summary(dry$Percent)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.0000 0.0000 0.0592 0.0700 0.7600
> require(proxyC)
(colZeros(dry$Percent)/length(dry$Percent))*100
[1] 63.29787
require(gamlss)
m1 <- gamlss(Percent ~ Month*CYR, family = BEZI, data = dry, trace = F)
summary(m1)
******************************************************************
Family: c("BEZI", "Zero Inflated Beta")
Call: gamlss(formula = Percent ~ Month * CYR, family = BEZI, data = dry, trace = F)
Fitting method: RS()
------------------------------------------------------------------
Mu link function: logit
Mu Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.49687 0.22213 -6.739 3.21e-11 ***
Month2 0.12405 0.37449 0.331 0.7406
Month3 -0.11640 0.24616 -0.473 0.6364
CYR2015 -0.09874 0.27365 -0.361 0.7183
CYR2016 -0.08591 0.27175 -0.316 0.7520
CYR2022 -0.29094 0.16122 -1.805 0.0715 .
Month3:CYR2015 0.27618 0.31985 0.863 0.3882
Month3:CYR2016 -0.02327 0.32402 -0.072 0.9428
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
------------------------------------------------------------------
Sigma link function: log
Sigma Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.93327 0.08568 22.56 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
------------------------------------------------------------------
Nu link function: logit
Nu Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.54502 0.07566 7.204 1.44e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
------------------------------------------------------------------
No. of observations in the fit: 752
Degrees of Freedom for the fit: 10
Residual Deg. of Freedom: 742
at cycle: 7
Global Deviance: 513.7825
AIC: 533.7825
SBC: 580.0098
******************************************************************
require(emmeans)
emmeans(m1, "Month", type = "response")
Error in V[idx, idx, drop = FALSE] : subscript out of bounds