Predict with pseudo-mean factors in new data

Question

R and Stata have different default behaviors when making predictions from a model that uses categorical/factor covariates. For example, if I want to predict outcomes for both levels of a two-level covariate factor (in this case foreign and domestic cars), holding all other values at their means, Stata's margins [varname], atmeans does weird stuff with factors, calculating the mean 0/1 value for each level:

. sysuse auto2

. reg price mpg i.foreign i.rep78

...

------------------------------------------------------------------------------
       price |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         mpg |  -299.6068   63.34525    -4.73   0.000    -426.2322   -172.9815
             |
     foreign |
    Foreign  |   1102.334   901.7772     1.22   0.226    -700.2928    2904.961
             |
       rep78 |
       Fair  |   841.3622   2055.452     0.41   0.684    -3267.428    4950.153
    Average  |   1285.116   1901.486     0.68   0.502    -2515.901    5086.132
       Good  |   1155.571   1984.561     0.58   0.562     -2811.51    5122.652
  Excellent  |   2353.179   2130.577     1.10   0.274    -1905.784    6612.142
             |
       _cons |   10856.24   2266.757     4.79   0.000      6325.06    15387.43
------------------------------------------------------------------------------

. margins foreign, atmeans

Adjusted predictions                            Number of obs     =         69
Model VCE    : OLS

Expression   : Linear prediction, predict()
at           : mpg             =    21.28986 (mean)
               0.foreign       =    .6956522 (mean)
               1.foreign       =    .3043478 (mean)
               1.rep78         =    .0289855 (mean)
               2.rep78         =     .115942 (mean)
               3.rep78         =    .4347826 (mean)
               4.rep78         =    .2608696 (mean)
               5.rep78         =    .1594203 (mean)

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     foreign |
   Domestic  |    5810.55    415.892    13.97   0.000     4979.194    6641.907
    Foreign  |   6912.884   700.8393     9.86   0.000     5511.927    8313.842
------------------------------------------------------------------------------

R, on the other hand, cannot calculate the mean of a factor (since that's mathematically impossible anyway), and it doesn't divide factors up into numerical proportions like Stata. Instead, when creating a new dataframe of covariates to pass into the model, I have to choose one of the factor levels:

library(haven)
auto <- read_stata("http://www.stata-press.com/data/r13/auto2.dta")

model <- lm(price ~ mpg + as.factor(foreign) + as.factor(rep78), data=auto)
summary(model)
#> Coefficients:
#>                     Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)         10856.24    2266.76   4.789 1.08e-05 ***
#> mpg                  -299.61      63.35  -4.730 1.34e-05 ***
#> as.factor(foreign)1  1102.33     901.78   1.222    0.226    
#> as.factor(rep78)2     841.36    2055.45   0.409    0.684    
#> as.factor(rep78)3    1285.12    1901.49   0.676    0.502    
#> as.factor(rep78)4    1155.57    1984.56   0.582    0.562    
#> as.factor(rep78)5    2353.18    2130.58   1.104    0.274    
#> ---

# Create new data with average values of all covariates for both foreign and
# domestic cars
newdata <- expand.grid(mpg = mean(auto$mpg, na.rm=TRUE),
                       foreign = c(0, 1),
                       rep78 = 3)  # One of the factor levels

# Not the same as Stata, obviously
predict(model, newdata=newdata)
#>        1        2 
#> 5760.544 6862.878

I'm using R to replicate a study that was originally done in Stata that used margins [varname], atmeans to generate predicted outcomes from a model with several categorical covariates. Is there a way to replicate the pseudo-mean factor value like Stata does (decomposing the factor into its individual levels, coded as dummy 0/1 values), or is there a more accurate way to use predict() with "average" categories in R (other than just arbitrarily choosing one of the levels)? Which approach (Stata's mean-of-each-level vs. R's choose-one-level) is more accurate/appropriate?

Answer 1

The correct way to use Stata's margins command in this context would be to add the over(foreign) option. Stata has a long tradition of taking the user literally rather than trying to interpret what might really be wanted. So if you add the atmeans option you ask for an evaluation at the means, and you get just that. If you want it separated by a factor, you add the over() option.

I agree that if you forget the over() option you will get results I would not want to present, but that is why Stata warns you by giving that rather "unStataish" additional output stating at which values the margins are evaluated.

Answer 2

My personal preference when calculating marginal effects at the mean (MEMs) would be to use the base value for factor variables and the mean for continuous ones. These are also sometimes called marginal effects at Representative Values (MERs).

For example,

sysuse auto2, clear
reg price mpg i.foreign i.rep78, coefl
margins foreign, at((means) mpg (base) rep78)

This is equivalent to:

sum mpg if e(sample)
di _b[_cons] + _b[mpg]*r(mean)  + _b[1.foreign]
di _b[_cons] + _b[mpg]*r(mean) 

Essentially, you are predicting once as if every car was foreign and then again as if every car was domestic, with the mpg set to the overall mean and rep78 set to "Poor" for all vehicles. This is different from predicting for foreign and domestic cars separately, with the covariates set to origin-specific means, which is what over() does.

In this age of big data, if you don't want to write out all the variables names, you can simply use this shorthand:

margins foreign, at((means) _continuous (base) _factor)

You can replicate what R is doing with:

margins foreign, at((means) _continuous rep78=3)

or by pre-specifying the base like this:

fvset base 3 rep78
reg price mpg i.foreign i.rep78, coefl
margins foreign, at((means) _continuous (base) _factor)

As far as best practices, I don't think there's a consensus, even within a particular fields. Personally, I prefer Average Marginal Effects (AMEs) to MEMs and MERs. That is what Stata calculates by default. There is an older SJ paper by Tamas Bartus that goes over the various choices in Section 2.