Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

model output for unquoted family argument (cosmetic) #31

Open
TimTaylor opened this issue Aug 4, 2020 · 3 comments
Open

model output for unquoted family argument (cosmetic) #31

TimTaylor opened this issue Aug 4, 2020 · 3 comments
Labels
enhancement New feature or request

Comments

@TimTaylor
Copy link
Member

Currently if an unquoted argument is given for the family parameter the formula in the resulting model is the entire function call. Borrowing from the glm code we can get a better input. There may be a better but this will work:

library(trendbreaker)

# current implementation
glm_poisson = glm_model(hp ~ 1 + cyl + drat + wt + qsec + am, poisson)
res <- glm_poisson$train(mtcars)
get_model(res)
#> 
#> Call:  glm(formula = hp ~ 1 + cyl + drat + wt + qsec + am, family = function (link = "log") 
#> {
#>     linktemp <- substitute(link)
#>     if (!is.character(linktemp)) 
#>         linktemp <- deparse(linktemp)
#>     okLinks <- c("log", "identity", "sqrt")
#>     family <- "poisson"
#>     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
#>     validmu <- function(mu) all(is.finite(mu)) && all(mu > 0)
#>     dev.resids <- function(y, mu, wt) {
#>         r <- mu * wt
#>         p <- which(y > 0)
#>         r[p] <- (wt * (y * log(y/mu) - (y - mu)))[p]
#>         2 * r
#>     }
#>     aic <- function(y, n, mu, wt, dev) -2 * sum(dpois(y, mu, 
#>         log = TRUE) * wt)
#>     initialize <- expression({
#>         if (any(y < 0)) stop("negative values not allowed for the 'Poisson' family")
#>         n <- rep.int(1, nobs)
#>         mustart <- y + 0.1
#>     })
#>     simfun <- function(object, nsim) {
#>         wts <- object$prior.weights
#>         if (any(wts != 1)) 
#>             warning("ignoring prior weights")
#>         ftd <- fitted(object)
#>         rpois(nsim * length(ftd), ftd)
#>     }
#>     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 = initialize, 
#>         validmu = validmu, valideta = stats$valideta, simulate = simfun), 
#>         class = "family")
#> }, data = data)
#> 
#> Coefficients:
#> (Intercept)          cyl         drat           wt         qsec           am  
#>    6.256607     0.085193    -0.009274     0.176041    -0.136688     0.036604  
#> 
#> Degrees of Freedom: 31 Total (i.e. Null);  26 Residual
#> Null Deviance:       958.3 
#> Residual Deviance: 116.5     AIC: 343.6

# add family check / conversion from glm

glm_model <- function(formula, family, ...) {
  if (is.character(family)) 
    family <- get(family, mode = "function", envir = parent.frame())
  if (is.function(family)) 
    family <- family()
  if (is.null(family$family)) {
    print(family)
    stop("'family' not recognized")
  }
  
  structure(
    eval(bquote(list(
      model_class = "glm",
      train = function(data) {
        model <- glm(formula = .(formula), family = .(family$family), data = data, ...)
        trendbreaker:::model_fit(model, formula)
      }
    ))),
    class = c("trendbreaker_glm", "trendbreaker_model")
  )
}

# better formula in model output
glm_poisson = glm_model(hp ~ 1 + cyl + drat + wt + qsec + am, family = poisson)
res <- glm_poisson$train(mtcars)
get_model(res)
#> 
#> Call:  glm(formula = hp ~ 1 + cyl + drat + wt + qsec + am, family = "poisson", 
#>     data = data)
#> 
#> Coefficients:
#> (Intercept)          cyl         drat           wt         qsec           am  
#>    6.256607     0.085193    -0.009274     0.176041    -0.136688     0.036604  
#> 
#> Degrees of Freedom: 31 Total (i.e. Null);  26 Residual
#> Null Deviance:       958.3 
#> Residual Deviance: 116.5     AIC: 343.6

# still works with character input for family
glm_poisson = glm_model(hp ~ 1 + cyl + drat + wt + qsec + am, family = "poisson")
res <- glm_poisson$train(mtcars)
get_model(res)
#> 
#> Call:  glm(formula = hp ~ 1 + cyl + drat + wt + qsec + am, family = "poisson", 
#>     data = data)
#> 
#> Coefficients:
#> (Intercept)          cyl         drat           wt         qsec           am  
#>    6.256607     0.085193    -0.009274     0.176041    -0.136688     0.036604  
#> 
#> Degrees of Freedom: 31 Total (i.e. Null);  26 Residual
#> Null Deviance:       958.3 
#> Residual Deviance: 116.5     AIC: 343.6

Created on 2020-08-04 by the reprex package (v0.3.0)

@TimTaylor TimTaylor added the enhancement New feature or request label Aug 4, 2020
@TimTaylor
Copy link
Member Author

actually, probably better to convert to character and let glm do the failing...

@TimTaylor
Copy link
Member Author

library(trendbreaker)

glm_model <- function(formula, family, ...) {
  if (!is.character(family)) {
    family <- deparse(substitute(family))
  }
  
  structure(
    eval(bquote(list(
      model_class = "glm",
      train = function(data) {
        model <- glm(formula = .(formula), family = .(family), data = data, ...)
        trendbreaker:::model_fit(model, formula)
      }
    ))),
    class = c("trendbreaker_glm", "trendbreaker_model")
  )
}


# better formula in model output
glm_poisson = glm_model(hp ~ 1 + cyl + drat + wt + qsec + am, family = poisson)
res <- glm_poisson$train(mtcars)
get_model(res)
#> 
#> Call:  glm(formula = hp ~ 1 + cyl + drat + wt + qsec + am, family = "poisson", 
#>     data = data)
#> 
#> Coefficients:
#> (Intercept)          cyl         drat           wt         qsec           am  
#>    6.256607     0.085193    -0.009274     0.176041    -0.136688     0.036604  
#> 
#> Degrees of Freedom: 31 Total (i.e. Null);  26 Residual
#> Null Deviance:       958.3 
#> Residual Deviance: 116.5     AIC: 343.6

# still works with character input for family
glm_poisson = glm_model(hp ~ 1 + cyl + drat + wt + qsec + am, family = "poisson")
res <- glm_poisson$train(mtcars)
get_model(res)
#> 
#> Call:  glm(formula = hp ~ 1 + cyl + drat + wt + qsec + am, family = "poisson", 
#>     data = data)
#> 
#> Coefficients:
#> (Intercept)          cyl         drat           wt         qsec           am  
#>    6.256607     0.085193    -0.009274     0.176041    -0.136688     0.036604  
#> 
#> Degrees of Freedom: 31 Total (i.e. Null);  26 Residual
#> Null Deviance:       958.3 
#> Residual Deviance: 116.5     AIC: 343.6

Created on 2020-08-04 by the reprex package (v0.3.0)

@TimTaylor
Copy link
Member Author

will be closed by trending

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

No branches or pull requests

1 participant