fit_with() is a pipe-friendly tool that applies a list of
formulas to a fitting function such as stats::lm().
The list of formulas is typically created with formulas().
Arguments
- data
A dataset used to fit the models.
- .f
A fitting function such as
stats::lm(),lme4::lmer()orrstanarm::stan_glmer().- .formulas
A list of formulas specifying a model.
- ...
Additional arguments passed on to
.f
Details
Assumes that .f takes the formula either as first argument
or as second argument if the first argument is data. Most
fitting functions should fit these requirements.
Examples
# fit_with() is typically used with formulas().
disp_fits <- mtcars %>% fit_with(lm, formulas(~disp,
additive = ~drat + cyl,
interaction = ~drat * cyl,
full = add_predictors(interaction, ~am, ~vs)
))
# The list of fitted models is named after the names of the list of
# formulas:
disp_fits$full
#>
#> Call:
#> .f(formula = disp ~ drat * cyl + am + vs, data = data)
#>
#> Coefficients:
#> (Intercept) drat cyl am vs
#> -164.83 27.58 75.84 -38.22 -15.53
#> drat:cyl
#> -6.97
#>
# Additional arguments are passed on to .f
mtcars %>% fit_with(glm, list(am ~ disp), family = binomial)
#> [[1]]
#>
#> Call: .f(formula = am ~ disp, family = function (link = "logit")
#> {
#> linktemp <- substitute(link)
#> if (!is.character(linktemp))
#> linktemp <- deparse(linktemp)
#> okLinks <- c("logit", "probit", "cloglog", "cauchit", "log")
#> family <- "binomial"
#> if (linktemp %in% okLinks)
#> stats <- make.link(linktemp)
#> else if (is.character(link)) {
#> stats <- make.link(link)
#> linktemp <- link
#> }
#> else {
#> if (inherits(link, "link-glm")) {
#> stats <- link
#> if (!is.null(stats$name))
#> linktemp <- stats$name
#> }
#> else {
#> stop(gettextf("link \"%s\" not available for %s family; available links are %s",
#> linktemp, family, paste(sQuote(okLinks), collapse = ", ")),
#> domain = NA)
#> }
#> }
#> variance <- function(mu) mu * (1 - mu)
#> validmu <- function(mu) all(is.finite(mu)) && all(mu > 0 &
#> mu < 1)
#> dev.resids <- function(y, mu, wt) .Call(C_binomial_dev_resids,
#> y, mu, wt)
#> aic <- function(y, n, mu, wt, dev) {
#> m <- if (any(n > 1))
#> n
#> else wt
#> -2 * sum(ifelse(m > 0, (wt/m), 0) * dbinom(round(m *
#> y), round(m), mu, log = TRUE))
#> }
#> simfun <- function(object, nsim) {
#> ftd <- fitted(object)
#> n <- length(ftd)
#> ntot <- n * nsim
#> wts <- object$prior.weights
#> if (any(wts%%1 != 0))
#> stop("cannot simulate from non-integer prior.weights")
#> if (!is.null(m <- object$model)) {
#> y <- model.response(m)
#> if (is.factor(y)) {
#> yy <- factor(1 + rbinom(ntot, size = 1, prob = ftd),
#> labels = levels(y))
#> split(yy, rep(seq_len(nsim), each = n))
#> }
#> else if (is.matrix(y) && ncol(y) == 2) {
#> yy <- vector("list", nsim)
#> for (i in seq_len(nsim)) {
#> Y <- rbinom(n, size = wts, prob = ftd)
#> YY <- cbind(Y, wts - Y)
#> colnames(YY) <- colnames(y)
#> yy[[i]] <- YY
#> }
#> yy
#> }
#> else rbinom(ntot, size = wts, prob = ftd)/wts
#> }
#> else rbinom(ntot, size = wts, prob = ftd)/wts
#> }
#> structure(list(family = family, link = linktemp, linkfun = stats$linkfun,
#> linkinv = stats$linkinv, variance = variance, dev.resids = dev.resids,
#> aic = aic, mu.eta = stats$mu.eta, initialize = binomInitialize(family),
#> validmu = validmu, valideta = stats$valideta, simulate = simfun,
#> dispersion = 1), class = "family")
#> }, data = data)
#>
#> Coefficients:
#> (Intercept) disp
#> 2.6308 -0.0146
#>
#> Degrees of Freedom: 31 Total (i.e. Null); 30 Residual
#> Null Deviance: 43.23
#> Residual Deviance: 29.73 AIC: 33.73
#>