Skip to content

Instantly share code, notes, and snippets.

@juanchiem
Created April 28, 2023 18:24
Show Gist options
  • Save juanchiem/48c95634921346eb6fa511d6e150237d to your computer and use it in GitHub Desktop.
Save juanchiem/48c95634921346eb6fa511d6e150237d to your computer and use it in GitHub Desktop.
library(tidyverse)
# https://stats.stackexchange.com/questions/309047/zero-inflated-beta-regression-using-gamlss-for-vegetation-cover-data
# 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:
df <- data.frame(
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),
site = as.factor(c("A","A","A","A","A",
"A","A","A","A","A","A","A","A","B","B",
"B","B","B","B","B","C","C","C","C","C","D",
"D","D","D","D","D","E","E","E")))
df %>%
ggplot() +
aes(site, cover) +
geom_boxplot() +
labs(x = "Site", y = "Vegetation cover") +
scale_y_continuous(labels = scales::percent)
library(gamlss)
library(emmeans)
m1 <- gamlss(cover ~ site, family = BEZI, data = df, trace = F)
summary(m1)
plot(m1)
em <- emmeans(m1, "site", type = "response")
pairs(em)
library(see)
ggplot(df,
aes(x = site, y = cover)) +
geom_violinhalf(flip = TRUE) +
scale_fill_material_d()
df %>%
ggplot()+
aes(x=cover, fill = site) +
geom_density(alpha = 0.4) +
facet_wrap("site")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment