Loading required packages

# Load packages
# Load packages
library(gamlss)
library(gamlss.dist)
library(gamlss.add)
library(gamlss.data)
library(gamlss.ggplots)
library(RelDists)
library(ggplot2)

# Bucket plot
# Load external script for bucket plotting functions (ensure bp.R is in your working directory)
# Fernanda De Bastiani, Mikis Stasinopoulos and Bob Rigby
source("bp.R")

# this code can be obtained via the paper
# Distributional modelling of positively skewed data via the flexible Weibull extension distribution (2024)
# https://doi.org/10.1111/anzs.12423

Conditional Fitting Models

# Simple additive linear model
model <- AI ~ age + gender + BMI

# Nonlinear model using lo() (local regression based on loess)
model.2 <- AI ~ lo(~age) + gender + lo(~BMI)

# Alternatively, one could use p-splines (e.g., pbm(age) + gender + pbm(BMI))

Fitting the Simple Models

# Fit the simple models
NO <- gamlss(model, family = NO, data = sleep, trace=F)
GA <- gamlss(model, family = GA, data = sleep, trace=F)
FWE <- gamlss(model, family = FWE, data = sleep, trace=F)

# Compare models using AIC
GAIC(NO, GA, FWE)
##     df      AIC
## FWE  5 761.6130
## GA   5 766.7911
## NO   5 864.2635

Worm Plot

model_wp(NO, GA, FWE) +
  theme_bw(base_size = 15) +
  theme(legend.position = "bottom")

Bucketplot

moment_bucket(NO, GA, FWE, no_bootstrap = 999) +
  ylab("moment excess kurtosis") +
  xlab("moment skewness") +
  theme_bw(base_size = 15) +
  theme(legend.position = "none")

Assessment of Variable Combinations

# For nicer labels
age <- sleep$age
gender <- as.factor(sleep$gender)
levels(gender) <- c("male", "female")
bmi <- sleep$BMI

# Gender only
bp.plot(list(NO, GA, FWE), xvar = ~gender,
            show.legend = TRUE, no.bootstrap = 999,
            text.to.show = c("N", "GA", "FWE"),
            overlap = 0)

# Gender and age
bp.plot(list(NO, GA, FWE), xvar = ~gender * age,
            show.legend = TRUE, no.bootstrap = 999,
            text.to.show = c("N", "GA", "FWE"),
            overlap = 0)

Fitting Models with Nonlinear Effects

# Fit models with nonlinear effects using lo()

NO.lo  <- gamlss(model.2, family = NO(),  data = sleep, trace=F)

GA.lo  <- gamlss(model.2, family = GA(),  data = sleep, trace=F)

FWE.lo <- gamlss(model.2, family = FWE(), data = sleep, trace=F)

# Compare models using AIC
GAIC(NO.lo, GA.lo, FWE.lo)
             df         AIC
GA.lo  12.06974    765.2226
NO.lo  12.06974    864.1669
FWE.lo 10.05350 612077.1944

Exploring Associations Between AI and Other Variables

Correlation Plot: Age and AI

ggplot(sleep, aes(x = age, y = AI)) +
  geom_point() +
  geom_smooth(method = "loess", se = TRUE, color = "red", linetype = "solid") +
  ggtitle("Correlation Plot: Age and AI") +
  xlab("Age") +
  ylab("AI")

Correlation Plot: BMI and AI

ggplot(sleep, aes(x = BMI, y = AI)) +
  geom_point() +
  geom_smooth(method = "loess", se = TRUE, color = "red", linetype = "solid") +
  ggtitle("Correlation Plot: BMI and AI") +
  xlab("BMI") +
  ylab("AI")