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

## ----data---------------------------------------------------------------------
library(panglm)
data(copd)
head(copd)

## ----pooled-------------------------------------------------------------------
fit_pool <- panglm(fev1 ~ treatment + age + smoker + crp, data = copd,
                   index = c("id", "visit"), model = "pooling",
                   vcov = "cluster")
summary(fit_pool)

## ----within-------------------------------------------------------------------
fit_fe <- panglm(fev1 ~ crp, data = copd, index = c("id", "visit"),
                 model = "within", vcov = "cluster")
summary(fit_fe)

## ----between------------------------------------------------------------------
fit_be <- panglm(fev1 ~ crp + age, data = copd, index = c("id", "visit"),
                 model = "between")
coef(fit_be)

## ----random-------------------------------------------------------------------
fit_re <- panglm(fev1 ~ treatment + age + smoker + crp, data = copd,
                 index = c("id", "visit"), model = "random")
coef(fit_re)

## ----cre----------------------------------------------------------------------
fit_cre <- panglm(fev1 ~ treatment + age + smoker + crp, data = copd,
                  index = c("id", "visit"), model = "random", cre = TRUE,
                  vcov = "cluster")
summary(fit_cre)

## ----hausman------------------------------------------------------------------
fit_re1 <- panglm(fev1 ~ crp, data = copd, index = c("id", "visit"),
                  model = "random")
panhausman(fit_fe, fit_re1)
panhausman(fit_fe, method = "mundlak", data = copd)

## ----vcov---------------------------------------------------------------------
sqrt(diag(vcov(fit_fe, type = "classical")))
sqrt(diag(vcov(fit_fe, type = "cluster")))
confint(fit_fe)

## ----ape----------------------------------------------------------------------
fit_hosp <- panglm(hospitalized ~ treatment + age + crp, data = copd,
                   index = c("id", "visit"), model = "random",
                   family = "probit", cre = TRUE, vcov = "cluster")
partialeffects(fit_hosp)

## ----predict------------------------------------------------------------------
new <- data.frame(id = 1, visit = 6, treatment = c(0, 1), age = 65,
                  smoker = 0, crp = 3)
predict(fit_re, newdata = new, type = "response")
predict(fit_fe, newdata = data.frame(id = 1:2, visit = 6, crp = 3))

## ----tidy, eval = has_generics------------------------------------------------
generics::tidy(fit_fe, conf.int = TRUE)
generics::glance(fit_fe)

## ----modelsummary, eval = has_ms----------------------------------------------
modelsummary::modelsummary(
  list(Pooled = fit_pool, "Fixed effects" = fit_fe, "Random effects" = fit_re),
  gof_map = c("nobs", "n.groups"), output = "markdown"
)

## ----plot, fig.alt = "Coefficient estimates of the correlated random effects model with 95% confidence intervals"----
plot(fit_cre, which = "coefficients")

## ----plot-resid, fig.alt = "Residuals against fitted values of the fixed-effects model"----
plot(fit_fe, which = "residuals")