leeper / margins Goto Github PK

Currently the package uses the delta method for variances. This should probably be the default (at least for OLS), but there are alternatives:

• Krinsky-Robb. Details: http://onlinelibrary.wiley.com/doi/10.1111/1475-6773.12122/pdf
• Bootstrapping. This is hard in practice because of use of the way atmeans is currently implemented. I suppose the bootstrapping would come first, then atmeans would be applied.

There's also the fact that Stata uses numerical derivatives, which I've avoided up to this point by using symbolic derivatives. That won't necessarily generalize well to other kinds of models.

Thanks for your work on this!

I was playing around with this cool package, and wondering why margins was so slow when applied to a small logistic glm with 100 observations. Basic profiling suggests that the function spends about 80% of its take taking the jacobian in the delta_once function.

I can't quite find my way around the code yet, but it looks like it's calculating the variance for every observation, which may not be necessary in every use case.

Any thoughts on how to speed things up? If you point me to the right direction, I may be able to PR something.

Hello,

This is a wonderful package. In Stata, you can run a regression with survey data, and then run margins, but with this package it appears that you cannot currently run margins for survey data.

So in Stata you can do something like this:

* california health interview survey 2014 teen dataset
* http://healthpolicy.ucla.edu/chis/data/public-use-data-file/Pages/2013-2014.aspx
use teen

svyset [pw=rakedw0], jkrw(rakedw1-rakedw80, multiplier(1)) vce(jack) mse

svy: regress bmi_p i.srsex##i.raceh2t_p1 srage_p

margins i.srsex i.raceh2t_p1

Perhaps I am doing something wrong, but this won't work in R:

library(margins)
library(haven)
library(survey)
# as.data.frame because tibbles break survey package
chis14 <- as.data.frame(read_dta("teen.dta")) 

chis14.svy <- svrepdesign(data=chis14,
                          weights=~rakedw0,
                          repweights="rakedw[1-9]",
                          scale=1,rscales=1,
                          type='other',combined.weights=TRUE,
                          mse=TRUE)

mylm <- svyglm(bmi_p~factor(srsex)*factor(raceh2t_p1)+srage_p,design=chis14.svy)

summary(mylm)

margins(mylm)
# Error in eval(expr, envir, enclos) : object 'srsex' not found

Issue is that prediction::find_data looks for data in the wrong place for svyglm objects.

I believe it shouldn't be too difficult to get margins working with the survey package. Would you consider a pull request?

Thanks

Allow someone to pass custom vcov value to margins() through the plotting commands:

• cplot()
• persp()
• image()

This error is uninformative. Can probably do this better:

> x2 <- lm(mpg ~ hp * wt, data = mtcars)
> m2 <- margins(x2, at = list(wt = wts, am = 0:1))
Error in setNames(est, c(names(est)[1], names(betas))) : 
  'names' attribute [5] must be the same length as the vector [4]

This will require passing information about the model type to the variance calculator.

See Greene p.734 for logit and probit formulae

This is only used in one place for jacobian() because the package is actually slower than just rewriting a vectorized numerical derivative. It shouldn't be difficult to replicate the jacobian() function using that same approach but for each coefficient.

https://cran.r-project.org/package=LinRegInteractive

I'm opening this up to follow-up on our twitter conversation. My package marfx replicates some of this package. I don't have time to go into detail now, but follow up later when I have more time to write my thoughts and how the two packages differ and how we can build on the strengths of each. But in the meantime, I wanted to open up this issue as a place for discussing this.

Add options for:

• Marginal effects (ME), the default
• Elasticities (ME * (x/y))
• Semi-elasticities (ME * (1/y))

In Stata this is:

dydx(): change in y for a change in x
eyex(): proportional change in y for a proportional change in x
eydx(): proportional change in y for a change in x
dyex(): change in y for a proportional change in x

Currently cplot() will only draw one line.

• It should have an at-type argument to allow predicted values or marginal effects across additional variables (e.g., to highlight effects among subgroups).
• It might also be nice to have it handle factor variables specified as x to draw predictions or effects as a dot plot with error bars.
• Update vignette here to use cplot().

As far as I can tell, the terms object in predict is only used for getting the original x-values, which aren't always wanted. Worse, a terms object will carry around it's entire environment when generated not in .GlobalEnv, which can be terrible for memory when working with big data and/or generating many models. See here https://stat.ethz.ch/pipermail/r-devel/2016-July/072924.html

My wish is that predict methods would allow me to manually pass newdata and only use the terms object when required, so that I can delete this from the original object and not run into these problems of carrying around all this data unnecessarily.

And, your package seems like a great place to fit this in for a variety of model objects.

The internals of persp.lm() can also be used to generate the input values for a contour() method. Thus:

• Abstract these internals into stand-alone functions
• Create image.lm() method that uses the same internals

The more I think about this, it is the wrong approach in almost every use case. If someone wants to do this, they should have to it manually, so next step is to remove all atmeans references.

• Extract the plotting code from cplot.lm() and make it standalone functions
  • A plot setup function
  • A line drawing function
• Create a cplot.polr() method that calls the drawing functions repeatedly for each outcome

Here's an example structure of a polr "prediction" as a reference:

> library("MASS")
> house.plr <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
> prediction(house.plr)
Modal prediction: "High" for 42 of 72 observations with total 3 levels
> str(prediction(house.plr))
Classes ‘prediction’ and 'data.frame':  72 obs. of  5 variables:
 $ fitted    : Factor w/ 3 levels "Low","Medium",..: 1 1 1 3 3 3 3 3 3 1 ...
 $ se.fitted :Classes 'se.fit', 'numeric'  num [1:72] NA NA NA NA NA NA NA NA NA NA ...
 $ Pr(Low)   : num  0.378 0.378 0.378 0.257 0.257 ...
 $ Pr(Medium): num  0.288 0.288 0.288 0.274 0.274 ...
 $ Pr(High)  : num  0.334 0.334 0.334 0.469 0.469 ...
 - attr(*, "model.class")= chr "polr"

Vignettes should import Stata datasets for examples.

Sorry, I placed this first as a comment to commit 03bdb72, but don not know if it gets through to you this way. So I opened this issue as suggestions:

Would it be possible to pass a vcov calculation function? This is fairly standard option, see, e.g. lmtest::coeftest.

I think, in other packages/functions the argument is called vcov or vcov., see packages lmtest, sandwich, car, ... (with the dot being there for imho historical reasons (pre NAMESPACE times) and allowing for partial argument matching). Thus, maybe renaming the vc argument could be beneficial from a consistency perspective?

This call to numDeriv::grad() is necessarily slow because it iterates over rows of the data and numDeriv::grad() further iterates over terms in the model.

A faster way to do this would be to create a function factory that expresses a vector of predicted values with respect to one variable. Then we can skip the for-loop over rows by simply doing a numerical derivative a la numDeriv::grad() but using vectorized algebra rather than a loop over values. We'll still need to loop over terms but in most cases there are going to be many fewer terms than there are observations.

Some pseudo-code:

eps <- .Machine$double.eps # make new formal argument to `marginal_effects()`
for (i in seq_along(nnames)) {
    FUN <- .build_predict_fun2(data = data, model = model, term = nnames[i], type = type)
    out[, nnames[i] ] <- ( FUN(data[[nnames[i]]]) - FUN(data[[nnames[i]]] + eps) ) / eps
}

Currently, an error is issued when trying to display the conditional effect(s) of a factor variable that has more than two levels in cplot():

> cplot(lm(y~x*fac,data=df), "x", "fac", what = "effect")
Error in cplot.lm(lm(y ~ x * fac, data = df), "x", "fac", what = "effect") : 
  Displaying on a factor variable with > 2 levels is not currently supported!

This could be fixed to display effect of each level/transition for the factor. margins() already does this fine, but it's a matter of figuring out how to draw it and return the data it appropriately.

Tests:

• predicted/fitted values match predict()
• Golder tests that allow direct calculation of effects and variances
• "response" versus "link" for GLMs
• non-base prediction() methods
• margins.loess()

Hi,

is there currently any way to view the vignettes in a nice format?

Trying vignette('Introduction', 'margins') results in the warning:

Warning message:
vignette ‘Introduction’ has no PDF/HTML 

http://cran.r-project.org/web/packages/interactionTest/index.html

I am trying to install the package on a new mac. MacOS 10.12.2, R version 3.3.2 (2016-10-31). Followed instructions on bottom of page. ghit and leeper/prediction were installed, I get the following error when trying to install margins:
``
install_github("leeper/margins")
Error in build_and_insert(p$pkgname, d, vers, build_args, verbose = verbose) :
Package build for margins failed!
''

It would be useful to have a plotting function that uses persp() to display a predicted value surface across two variables.

https://cran.r-project.org/web/packages/condvis/condvis.pdf

When using margins.glm with type = "response", a call to predict.glm is made. This is fine, except when using factor() variables because they are not included in model.matrix(x) (which extracts the relevant factor levels). The default should use data contained inside the model object, unless newdata has been specified rather than the output of model.matrix(x).

Same issue as #39

Add plot vignette that shows:

• Marginal effect dot plots
• Marginal effects conditional on a variable using base
• Marignal effects conditional on a variable using ggplot2

All in the title. Thanks in advance!

Despite resolving #8, factors are still problematic:

> x <- lm(mpg ~ factor(cyl), data = mtcars)
> margins(x, at = list(cyl = c(2,4,6)))
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
  contrasts can be applied only to factors with 2 or more levels

The issue is the use of model.matrix to build the design matrix. If multiple factor levels aren't present (and they aren't because at_builder will specify only one level for the factor variable), then model.matrix fails.

One strategy might be to rename the factor variables all the way through the process rather than just inside .margins. That will make using predict really difficult, though.

You state:

Note: While margins also implements fitted values from model estimation results, that functionality is already well-implemented by R's predict function, so that wheel is not reinvented here.

This is true, but this does not give you standard errors or confidence intervals, which stata's margins command does. So, it would still be nice to have this functionality at some point.

> x <- lm(mpg ~ factor(cyl)*am, data = mtcars)
> margins(x)
Error in `[.data.frame`(x$model, , attributes(terms(x))$term.labels[i]) : 
  undefined columns selected

Problem seems to be in handling the termorder object inside .margins

Currently everything is assumed to be continuous, which is not really a problem in OLS, but is a problem for GLMs. It should be an option (the default?, as in Stata) to do discrete differences rather than derivatives for factor variables.

x <- glm(am ~ factor(cyl) + hp + wt, data = mtcars, family = binomial)
out <- lapply(c(4,6,8), function(z) predict(x, newdata = within(mtcars, cyl <- z), type = 'response', se.fit = TRUE))
> mean(out[[3]]$fit - out[[2]]$fit)
[1] -0.4676555
> mean(out[[2]]$fit - out[[1]]$fit)
[1] 0.1197979
> mean(out[[3]]$fit - out[[1]]$fit)
[1] -0.3478576

Greene gives variance formulae. But what are the SEs when these are included in a discrete-by-continuous interaction.

Rather than return the model object as an attribute, return a data.frame consisting of:

• Original variables
• Predicted value for each case (numeric value with custom "fittedvalue" class)
• Marginal effects at each case value (numeric value with custom "marginaleffect" class)

(for atmeans = TRUE, this should be a data.frame with one row)

In attributes:

• variances
• current attributes
• model formula/call

This has the advantage of being able to easily do.call(rbind, margins(x)) to create a single data.frame with the marginal effects estimated at different specified at values, because marginal_effect() will return a data.frame rather than a list.

• loess
  • predict() method uses se rather than se.fit
• plm/pglm
  • no predict() method
• MASS::polr()
  • has predict() method with type = c("class", "probs")
• MASS::glm.nb()
  • should work out of the box
• lme4 (lmer(), glmer(), ...)
  • lme4:::predict.merMod() has no se.fit argument
• nlme
  • nlme:::predict.lme() has no se.fit or type arguments
  • nlme:::predict.gls() has no se.fit argument.
• nls
  • requires setting model = TRUE in original call; and then getting term(object[["model"]])
• lfe::felm()
• survey
  • (implements alternative modelling functions for non-SRS survey datasets)
• survival
  • coxph class has different type values for predict(): c("lp", "risk", "expected")
  • survival::survreg()
• censReg
  • has a margEff() generic but no predict() method
• nnet::multinom() or any nnet model generally
• AER::ivreg()
• quantreg::rq()
• Document which model types have prediction() and marginal_effects()/margins() methods

The current approach works for "lm" and "glm" objects, but other classes have slightly different behavior.

Weighted survey data effects in Stata: http://www.stata.com/products/stb/journals/stb31.pdf

"Predictive Margins with Survey Data" https://www.amstat.org/sections/SRMS/Proceedings/papers/1997_111.pdf

Hi,

I'm super happy that one of Stata's biggest strengths will soon be available in R for free. What I would suggest for future versions of the package though is to add some convenience functions/parameters that make it easier to draw more sophisticated graphs without having to write a lot of code.

As an example I imported the diamonds dataset from R, created some variables for a logit model and have a lot of trouble replicating the following Stata graph:

* Importing the data:
import delimited "diamonds.csv", delimiter(comma) varnames(1) clear 

* Generate dep. variable for logistic regression
recode carat (min/1 = 0 "low") (else= 1 "high"), gen(carat_dummy)

* Generate categorical (factor variable)
gen cut_cat = 1
replace cut_cat = 2 if cut == "Very Good" | cut == "Premium"
replace cut_cat = 3 if cut == "Ideal"

* Generate dummy for color
gen col_g = color == "D"

* Fitting the model
logit carat_dummy i.cut_cat c.price c.depth  i.col_g

* Calculating margins for different depths and only "D" colored diamonds for each cut category
margins, at(depth=(40(5)80) col_g= 1) over(cut_cat)
marginsplot

I assume that it is currently not possible to replicate the outcome of these 2 lines of Stata margins without a large code chunk. It would really be awesome if you could write additional functions/parameters to make this process a little easier.

Anyway, many thanks for putting effort in this project :)

Need to account for poly in formulae.

> margins(lm(mpg ~ cyl + poly(hp, 2), data = mtcars))
Error in parse(text = termtext) : 
  <text>:1:57: unexpected numeric constant
1:  ~ -1.23647789531276*cyl + -15.8106255224777*poly(hp, 2)1
                                                            ^

• Issue warning when weights are present about them being ignored
• Provide options for how to handle
  • Stata discussion is here: http://www.stata.com/manuals13/rmargins.pdf
  • Stata's default is a weighted mean based upon those weights, with an option to noweight (ignore weights used in estimation)

http://cran.r-project.org/web/packages/interplot/index.html

hsbdemo <- read.csv("http://www.ats.ucla.edu/stat/data/hsbdemo.csv")
logit1  <- glm("honors ~ I(-1*I(female=='female')) + I(-1*read)", data = hsbdemo, family="binomial")
summary(logit1) # model looks good

marginal_effects(logit1)

Error in [[<-.data.frame(*tmp*, variable, value = numeric(0)) :
replacement has 0 rows, data has 200

Install fail, plus debug trail:

> install_github("leeper/margins")
Error in build_and_insert(p$pkgname, d, vers, build_args, verbose = verbose) : 
  Package build for margins failed!
Called from: build_and_insert(p$pkgname, d, vers, build_args, verbose = verbose)
Browse[1]> arg
[1] "CMD build"                                                                       
[2] ""                                                                                
[3] "/var/folders/gl/ds8cxcyj07x0f553zt__04mm0000gn/T//RtmpH05rkw/margins2f4b1edc38b0"
Browse[1]> rpath
[1] "/Library/Frameworks/R.framework/Resources/bin/R"
Browse[1]> ## try the offending line again, because reasons...
Browse[1]> system2(rpath, arg, stdout = if (verbose) 
+     ""
+   else FALSE)
debug: [1] FALSE

Session info:

R version 3.3.1 (2016-06-21)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.11.6 (El Capitan)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
[1] rsconnect_0.4.3 tools_3.3.1     git2r_0.15.0    ghit_0.2.12    

• Manual calculation of examples from Matt Golder's multiplicative interaction terms page
• Vignette comparing output to output from other packages
  • car
  • effects
  • mfx
  • erer
  • visreg
  • DAmisc
  • interplot
  • interactionTest
  • plotMElm
• Technical details vignette
• Example showing integration with non-standard variance-covariance matrices from sandwich

Should be straightforward, but this has a prerequisite that I return unit-level SEs from marginal_effects().