## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
 collapse = TRUE,
 comment = "#>"
)

## ----TWexample----------------------------------------------------------------
library(tweedie)
library(statmod)   # For the tweedie() family function, for use in glm()

set.seed(96)
N <- 25
# Poisson mean (lambda): the gamma summands have shape 2 and scale 1
lambda <- 1.5
# Generating Compound Poisson-Gamma data manually
y <- replicate(N, {
 n_events <- rpois(1, lambda = lambda)
 if (n_events == 0) 0 else sum(rgamma(n_events, shape = 2, scale = 1))
})

mod.tw <- glm(y ~ 1, 
             family = statmod::tweedie(var.power  = 1.5, 
                                       link.power = 0) )
             # link.power = 0  means the log-link

## ----TWqqplot-----------------------------------------------------------------
qqnorm( qr <- statmod::qresid(mod.tw),
       las = 1)
qqline(qr,
      col = "grey")

## ----TWrandom-----------------------------------------------------------------
tweedie::rtweedie(10, xi = 1.4, mu = 2, phi = 1)

## ----TWplotsPDF---------------------------------------------------------------
y_grid <- seq(0, 2, 
              length = 100)
xi     <- 1.05
mu     <- 0.5
phi    <- 0.4

twden <- tweedie::dtweedie(y_grid, xi = xi, mu = mu, phi = phi)  
twdtn <- tweedie::ptweedie(y_grid, xi = xi, mu = mu, phi = phi)

plot( twden[y_grid > 0] ~ y_grid[y_grid > 0], 
     xlab = expression(italic(y)),
     ylab = "Density function",
     type ="l",
     lwd  = 2,
     las  = 1)
points(twden[y_grid==0] ~ y_grid[y_grid == 0],
     lwd = 2,
     pch = 19)

## ----TWplotsCDF---------------------------------------------------------------
plot(twdtn ~ y_grid,
    xlab = expression(italic(y)),
    ylab = "Distribution function",
    type = "l",
    las  = 1,
    lwd  = 2,
    ylim = c(0, 1) )
points(twdtn[y_grid==0] ~ y_grid[y_grid == 0],
     lwd = 2,
     pch = 19)

## ----TWplotsPDF2--------------------------------------------------------------
tweedie::tweedie_plot(y_grid, xi = xi, mu = mu, phi = phi,
                      plot_args = list(
                        ylab = "Density function",
                        xlab = expression(italic(y)),
                        las  = 1,
                        lwd  = 2
                        ) 
)

## ----TWplotsCDF2--------------------------------------------------------------
tweedie::tweedie_plot(y_grid, xi = xi, mu = mu, phi = phi, 
                      type      = "cdf",
                      plot_args = list(
                        ylab = "Distribution function",
                        xlab = expression(italic(y)),
                        las  = 1, 
                        lwd  = 2, 
                        ylim = c(0, 1)) 
                        )

## ----TWprofile----------------------------------------------------------------
y <- rtweedie(50, xi = 1.5, mu = 1, phi = 1)
out <- tweedie_profile(y ~ 1,
                       xi.vec  = seq(1.3, 1.8, length = 5), 
                       do.plot = TRUE)

# The estimated power index:
out$xi.max

# The 95% CI for the power index:
out$ci


## ----PSNLoad------------------------------------------------------------------
poison <- read.csv(system.file("extdata", 
                               "poison.csv", 
                               package = "tweedie"))

# Convert to factors:
poison$Psn  <- factor(poison$Psn)
poison$Trmt <- factor(poison$Trmt)

head(poison)

## ----PSNPlotPsn---------------------------------------------------------------
plot(Time ~ Psn,
    data = poison,
    xlab = "Poison type",
    ylab = "Survival time (tens of hours)",
    las  = 1)

## ----PSNPlotTmt---------------------------------------------------------------
plot(Time ~ Trmt,
    data = poison,
    xlab = "Treatment type",
    ylab = "Survival time (tens of hours)",
    las  = 1)

## ----PSNProfile---------------------------------------------------------------
PSNPrf <- tweedie_profile(Time ~ Trmt + Psn, 
                          data    = poison, 
                          do.plot = TRUE, 
                          xi.vec  = seq(2.5, 5.5, length = 7) )

## ----PSNXi--------------------------------------------------------------------
PSNPrf$xi.max

## ----PSNModel-----------------------------------------------------------------
PSNMod <- glm(Time ~ Trmt + Psn, 
             data   = poison,
             family = statmod::tweedie(var.power  = round(PSNPrf$xi.max, 1), 
                                       link.power = 0))
anova(PSNMod, 
      test       = "F")

## ----PSNRes-------------------------------------------------------------------
qr <- statmod::qresid(PSNMod)
qqnorm(qr,
      las = 1)
qqline(qr,
      col = "grey")

## ----QRainLoad----------------------------------------------------------------
quilpie <- read.csv(system.file("extdata", 
                                "quilpie.csv", 
                                package = "tweedie"))

# Convert to factors:
quilpie$Phase  <- factor(quilpie$Phase)

head(quilpie)

## ----QRainPlot----------------------------------------------------------------
plot(Rain ~ Phase,
    data = quilpie,
    xlab = "SOI Phase",
    ylab = "July rainfall (mm)",
    las  = 1)

## ----QRainZero----------------------------------------------------------------
cat("Percentage of Julys with zero rain:", 
    round( sum(quilpie$Rain == 0) / length(quilpie$Rain) * 100, 1),
    "\n")

## ----QRainProfile-------------------------------------------------------------
QRainPrf <- tweedie_profile(Rain ~ Phase, 
                          data    = quilpie, 
                          do.plot = TRUE, 
                          xi.vec  = seq(1.2, 1.7, length = 7) )

## ----QRainXi------------------------------------------------------------------
QRainPrf$xi.max

## ----QRainModel---------------------------------------------------------------
QRainMod <- glm(Rain ~ Phase, 
             data   = quilpie,
             family = statmod::tweedie(var.power  = round(QRainPrf$xi.max, 1), 
                                       link.power = 0))
anova(QRainMod, 
      test       = "F")

## ----QRainRes-----------------------------------------------------------------
qr <- statmod::qresid(QRainMod)
qqnorm(qr,
      las = 1)
qqline(qr,
      col = "grey")

## ----QrainConvert-------------------------------------------------------------
out <- tweedie_convert(xi = QRainPrf$xi.max,
                       mu = fitted(QRainMod),
                       phi = summary(QRainMod)$dispersion )
out2 <- cbind( 
  Year   = quilpie$Year,
  lambda = round(out$poisson.lambda, 2), # Poisson mean number of rain events
  shape  = round(out$gamma.shape, 2),    # Gamma shape
  p0     = round(out$p0, 2),             # Pr(Rainfall = 0) = exp( -lambda)
  scale  = round(out$gamma.scale, 2)     # Gamma scale
)
head(out2, 10)