## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 6, fig.height = 4)

## ----setup--------------------------------------------------------------------
library(panglm)

## ----count-data---------------------------------------------------------------
set.seed(42)
N <- 200; TT <- 5
d <- data.frame(id = rep(1:N, each = TT), year = rep(1:TT, N))
a <- rnorm(N, sd = 0.6)
d$age <- rep(round(runif(N, 50, 80)), each = TT) + d$year - 1
d$female <- rep(rbinom(N, 1, 0.5), each = TT)
d$chronic <- rbinom(N * TT, 1, plogis(-0.5 + a[d$id]))
mu <- exp(0.5 * d$chronic + 0.03 * (d$age - 65) + 0.2 * d$female + a[d$id])
d$visits <- rnbinom(N * TT, mu = mu, size = 2)

## ----fe-poisson---------------------------------------------------------------
fe_pois <- panglm(visits ~ chronic + age, data = d, index = c("id", "year"),
                  model = "within", family = "poisson", vcov = "cluster")
summary(fe_pois)

## ----dispersion---------------------------------------------------------------
pandispersiontest(fe_pois)

## ----negbin-------------------------------------------------------------------
fe_nb <- panglm(visits ~ chronic + age, data = d, index = c("id", "year"),
                model = "within", family = "negbin")
coef(fe_nb)
fe_nb$theta
re_nb <- panglm(visits ~ chronic + age + female, data = d,
                index = c("id", "year"), model = "random", family = "negbin")
coef(re_nb)

## ----re-poisson---------------------------------------------------------------
re_pois <- panglm(visits ~ chronic + age + female, data = d,
                  index = c("id", "year"), model = "random",
                  family = "poisson", cre = TRUE, vcov = "cluster")
partialeffects(re_pois)

## ----hurdle-------------------------------------------------------------------
data(copd)
hur <- panhurdle(exacerbations ~ crp, data = copd, index = c("id", "visit"))
summary(hur)

## ----cost-data----------------------------------------------------------------
set.seed(11)
N <- 150; TT <- 4
cst <- data.frame(id = rep(1:N, each = TT), year = rep(1:TT, N))
a <- rnorm(N, sd = 0.5)
cst$age <- rep(round(runif(N, 40, 85)), each = TT)
cst$chronic <- rbinom(N * TT, 1, plogis(-0.5 + a[cst$id]))
p_any <- pnorm(-0.3 + 0.6 * cst$chronic + 0.01 * (cst$age - 60) + a[cst$id])
mu <- exp(7 + 0.5 * cst$chronic + 0.01 * (cst$age - 60) + 0.5 * a[cst$id])
cst$cost <- ifelse(runif(N * TT) < p_any,
                   rgamma(N * TT, shape = 2, rate = 2 / mu), 0)
mean(cst$cost == 0)

## ----twopart------------------------------------------------------------------
tp <- twopart(cost ~ chronic + age, data = cst, index = c("id", "year"),
              model = "pooling", vcov = "cluster")
summary(tp)
partialeffects(tp)

## ----beta---------------------------------------------------------------------
m <- plogis(1 - 0.4 * cst$chronic - 0.01 * (cst$age - 60) + a[cst$id])
cst$pdc <- rbeta(N * TT, m * 20, (1 - m) * 20)
cst$age10 <- (cst$age - 60) / 10
fit_beta <- panglm(pdc ~ chronic + age10, data = cst, model = "random",
                   index = c("id", "year"), family = "beta", R = 5)
summary(fit_beta)
partialeffects(fit_beta)

## ----ordered------------------------------------------------------------------
copd$severity <- cut(copd$exacerbations, c(-Inf, 0, 2, Inf),
                     labels = c("none", "mild", "severe"),
                     ordered_result = TRUE)
fit_ord <- panglm(severity ~ treatment + age + crp, data = copd,
                  index = c("id", "visit"), model = "random",
                  family = orderedfamily("logit"), R = 8)
summary(fit_ord)
fe_ord <- panglm(severity ~ crp, data = copd, index = c("id", "visit"),
                 model = "within", family = orderedfamily("logit"))
coef(fe_ord)

## ----ordered-ape--------------------------------------------------------------
partialeffects(fit_ord)
head(predict(fit_ord, type = "prob"), 3)

## ----tobit-data---------------------------------------------------------------
set.seed(7)
N <- 80; TT <- 4
eq <- data.frame(id = rep(1:N, each = TT), wave = rep(1:TT, N))
u <- rnorm(N, sd = 0.15)
eq$age <- rep(round(runif(N, 40, 85)), each = TT)
eq$chronic <- rpois(N * TT, 0.5 + 0.02 * (eq$age - 40))
ystar <- 1.15 - 0.006 * (eq$age - 40) - 0.08 * eq$chronic + u[eq$id] +
  rnorm(N * TT, sd = 0.12)
eq$eq5d <- pmin(ystar, 1)
mean(eq$eq5d == 1)

## ----tobit, warning = FALSE---------------------------------------------------
fit_tob <- panglm(eq5d ~ age + chronic, data = eq, index = c("id", "wave"),
                  model = "random", family = tobitfamily(left = -Inf, right = 1),
                  R = 5)
coef(fit_tob)
c(sigma_e = fit_tob$sigma_e, sigma_u = fit_tob$sigma)
fit_lin <- panglm(eq5d ~ age + chronic, data = eq, index = c("id", "wave"),
                  model = "random")
coef(fit_lin)