# Predict with pseudo-mean factors in new data

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)

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?

• Side comment. Calculating the mean of a factor (taking "factor" here to mean a coded categorical variable) is not necessarily impossible mathematically . Sometimes it even makes sense, e.g. if a binary factor is coded numerically as 0 or 1, as is very common in much software, then its mean is defined and central. What is allowed or not allowed specifically in R is not the same thing (although clearly there is a logic to it). Commented Apr 29, 2016 at 10:43
• This might be nitpicky, but it seems strange to think of this as default behavior when you are using the atmeans option. Commented Apr 29, 2016 at 20:16
• "default" in that not specifying over() or at() or anything else results in the proportion of each category level. So, the default atmeans behavior… Commented Apr 29, 2016 at 20:29
• You can use the effects package in R, which "averages" over the levels of a factor (calculating the "mean"). Commented Aug 10, 2018 at 8:22

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.

• So it is not correct/best practice to use Stata's fake categorical means then, and instead make predictions at each factor level (or just select one level)? Commented Apr 29, 2016 at 14:12
• There is nothing fake about those means, they are just proportions. It is not what I would use in such cases, but others are allowed to have different opinions. Commented Apr 30, 2016 at 8:40

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.

• How would you implement Average Marginal Effects in R?
– RNB
Commented Sep 21, 2016 at 9:55
• @RNB There is R port of Stata's margins command called leeper. You can also do a lot of the MEs by hand. Standard errors will be a bit harder, especially for non-linear models, with this route. Commented Sep 21, 2016 at 17:57