---
title: "Quick Usage Guide to the 'fastglm' Package"
author: "Jared Huling"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Quick Usage Guide to the 'fastglm' Package}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning = FALSE,
  echo = TRUE,
  eval = requireNamespace("glm2", quietly = TRUE)
)
```

Generalized linear models relate an outcome mean to predictors while allowing response distributions such as binomial or Poisson. The *fastglm* package performs the fitting calculations in C++ and provides several ways to solve the weighted regressions used during fitting. *fastglm* is compatible with `R`'s `family` objects (see `?family`). The package can be installed via

```{r, eval = FALSE, echo=TRUE}
pak::pak("fastglm")
# pak::pak("jaredhuling/fastglm") #development version
```

and loaded via:

```{r, echo=TRUE}
library(fastglm)
```

## Example

The main function, `fastglm()`, accepts either a model formula and data frame or a predictor matrix and response vector. A matrix is convenient when the predictor columns have already been constructed. The following example illustrates that interface:

```{r, echo = TRUE}
data(esoph)
x <- model.matrix(~ agegp + unclass(tobgp) + unclass(alcgp),
                  data = esoph)

y <- cbind(esoph$ncases, esoph$ncontrols)

gfit1 <- fastglm(x = x, y = y, family = binomial(link = "cloglog"))

summary(gfit1)
```

A formula can be passed directly to `fastglm()`. The fitting routine can also be used inside `glm()` by supplying `method = fastglm_fit`, which is useful when existing analysis code already uses `glm()`:

```{r}
gfit2 <- glm(cbind(ncases, ncontrols) ~ agegp + unclass(tobgp) +
                 unclass(alcgp), data = esoph,
             family = binomial(link = "cloglog"),
             method = fastglm_fit)

summary(gfit2)
```

The coefficient summary uses the covariance calculated during fitting, evaluated at the reported coefficient estimates. It therefore obtains standard errors without fitting the model again.

## Computational stability

A GLM is usually fitted through a sequence of weighted regression updates, called iteratively reweighted least squares (IRLS). An update can increase the discrepancy between the fitted and observed outcomes or produce means outside the distribution's allowed range. The *fastglm* safeguards shorten such steps until an acceptable update is found. This step-halving procedure can improve convergence, although difficult models can still fail.

The following example deliberately fits an inappropriate link, the function relating the response mean to the linear predictor. It illustrates numerical difficulties that can arise when proposed means leave the allowed range. Different starting values and step safeguards can lead to different convergence behavior. Numerical convergence alone does not make the chosen response distribution or link scientifically appropriate.

```{r, fig.show='hold'}
set.seed(1)
x <- matrix(rnorm(10000 * 100), ncol = 100)
y <- exp(0.25 * x[,1] - 0.25 * x[,3] + 0.5 * x[,4] - 0.5 * x[,5] + rnorm(10000)) + 0.1


# fastglm
system.time(gfit1 <- glm(y ~ x, family = Gamma(link = "sqrt"),
                         method = fastglm_fit))

# glm
system.time(gfit2 <- glm(y ~ x, family = Gamma(link = "sqrt")))

# glm2
system.time(gfit3 <- glm(y ~ x, family = Gamma(link = "sqrt"),
                         method = glm2::glm.fit2))

## Note that fastglm() returns estimates with the
## largest likelihood
logLik(gfit1)
logLik(gfit2)
logLik(gfit3)

coef(gfit1)[1:5]
coef(gfit2)[1:5]
coef(gfit3)[1:5]

## check convergence of fastglm
gfit1$converged
## number of IRLS iterations
gfit1$iter

## now check convergence for glm()
gfit2$converged
gfit2$iter

## check convergence for glm2()
gfit3$converged
gfit3$iter


## increasing number of IRLS iterations for glm() does not help that much
system.time(gfit2 <- glm(y ~ x, family = Gamma(link = "sqrt"),
                         control = list(maxit = 100)))

gfit2$converged
gfit2$iter

logLik(gfit1)
logLik(gfit2)

```