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

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

## ----rs-data------------------------------------------------------------------
set.seed(3)
N <- 100; TT <- 6
f <- data.frame(id = rep(seq_len(N), each = TT), year = rep(0:(TT - 1), N))
f$smoker <- rep(rbinom(N, 1, 0.4), each = TT)
b0 <- rnorm(N, sd = 0.4)
b1 <- rnorm(N, sd = 0.04)
f$fev1 <- 3 + b0[f$id] + (-0.05 - 0.03 * f$smoker + b1[f$id]) * f$year -
  0.2 * f$smoker + rnorm(N * TT, sd = 0.1)
f$exac <- rpois(N * TT, exp(-0.3 + 0.3 * f$smoker +
                              (0.05 + 3 * b1[f$id]) * f$year + b0[f$id]))

## ----rs-----------------------------------------------------------------------
rs <- panglm(fev1 ~ year * smoker, data = f, index = c("id", "year"),
             model = "random", random = ~ year)
summary(rs)
head(rs$ranef, 3)

## ----rs-lrt-------------------------------------------------------------------
ri <- panglm(fev1 ~ year * smoker, data = f, index = c("id", "year"),
             model = "random")
lr <- 2 * (as.numeric(logLik(rs)) - as.numeric(logLik(ri)))
c(LR = lr, p = pchisq(lr, df = 2, lower.tail = FALSE))

## ----rs-pois------------------------------------------------------------------
rp <- panglm(exac ~ year + smoker, data = f, index = c("id", "year"),
             model = "random", family = "poisson", random = ~ year, R = 5)
summary(rp)

## ----fmm-data-----------------------------------------------------------------
set.seed(1)
N <- 300; TT <- 4
m <- data.frame(id = rep(seq_len(N), each = TT), year = rep(seq_len(TT), N))
income <- rnorm(N)
high <- rbinom(N, 1, plogis(-1 - 0.8 * income))
m$income <- income[m$id]
m$chronic <- rbinom(N * TT, 1, 0.4)
m$age <- rep(round(runif(N, 40, 80)), each = TT)
m$visits <- rpois(N * TT, exp(ifelse(high[m$id] == 1,
                                     1.5 + 0.6 * m$chronic,
                                     0 + 0.3 * m$chronic) +
                                0.01 * (m$age - 60)))

## ----fmm-select---------------------------------------------------------------
panfmmselect(visits ~ chronic + age, data = m, index = c("id", "year"),
             k = 1:3, nstart = 1, seed = 1, tol = 1e-4)

## ----fmm----------------------------------------------------------------------
fmm <- panfmm(visits ~ chronic + age, data = m, index = c("id", "year"),
              k = 2, family = "poisson", concomitant = ~ income,
              nstart = 3, seed = 1)
summary(fmm)

## ----fmm-post, message = FALSE------------------------------------------------
first <- !duplicated(m$id)
post <- predict(fmm, type = "posterior")[first, ]
assigned <- predict(fmm, type = "class")[first]
all(assigned == max.col(post))
table(true = high, assigned = assigned)

## ----qr-data------------------------------------------------------------------
set.seed(1)
N <- 200; TT <- 5
q <- data.frame(id = rep(seq_len(N), each = TT), year = rep(seq_len(TT), N))
alpha <- rep(rnorm(N), each = TT)
q$chronic <- rpois(N * TT, 1 + 0.5 * (alpha > 0))
q$age <- rep(sample(40:70, N, TRUE), each = TT) + q$year
q$cost <- 2 + alpha + 0.5 * q$chronic + 0.02 * q$age +
  (1 + 0.3 * q$chronic) * rexp(N * TT)

## ----qr, eval = has_qr--------------------------------------------------------
qr <- panqr(cost ~ chronic + age, data = q, index = c("id", "year"),
            tau = c(0.1, 0.25, 0.5, 0.75, 0.9), model = "mmqr")
summary(qr)

## ----qr-plot, eval = has_qr, fig.alt = "Coefficient of chronic and age across quantiles with confidence bands"----
plot(qr)

## ----qr-canay, eval = has_qr, warning = FALSE---------------------------------
cn <- panqr(cost ~ chronic + age, data = q, index = c("id", "year"),
            tau = c(0.5, 0.9), model = "canay", boot = 50, seed = 1)
coef(cn)
coef(panqr(cost ~ chronic + age, data = q, index = c("id", "year"),
           tau = c(0.5, 0.9), model = "pooling"))