Set up data

rm(list=ls())
data(mtcars)
mtcars$wt <- mtcars$wt/sum(mtcars$wt)

# devtools::install_github('anniejw6/modmarg', ref = 'jacob_weights')

Gaussian

Model Coefficents

Stata

## . reg vs i.gear mpg [pweight = wt]
## (sum of wgt is   1.0000e+00)
## 
## Linear regression                                      Number of obs =      32
##                                                        F(  3,    28) =   31.63
##                                                        Prob > F      =  0.0000
##                                                        R-squared     =  0.5525
##                                                        Root MSE      =  .34221
## 
## ------------------------------------------------------------------------------
##              |               Robust
##           vs |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
## -------------+----------------------------------------------------------------
##         gear |
##           4  |    .333215   .2398524     1.39   0.176    -.1581005    .8245304
##           5  |  -.2160992   .1424822    -1.52   0.141    -.5079608    .0757624
##              |
##          mpg |   .0419613   .0145497     2.88   0.007     .0121576    .0717651
##        _cons |  -.4972822   .1911988    -2.60   0.015    -.8889352   -.1056292
## ------------------------------------------------------------------------------

R

mod <- glm(vs ~ as.factor(gear) + mpg, data = mtcars, family = 'gaussian',
           weights = wt)

lmtest::coeftest(mod, vcov = sandwich::vcovHC(mod, "HC1"),
                 df = mod$df.residual)
## 
## t test of coefficients:
## 
##                   Estimate Std. Error t value Pr(>|t|)   
## (Intercept)      -0.497282   0.191199 -2.6009 0.014686 * 
## as.factor(gear)4  0.333215   0.239852  1.3892 0.175707   
## as.factor(gear)5 -0.216099   0.142482 -1.5167 0.140559   
## mpg               0.041961   0.014550  2.8840 0.007469 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Margins (Levels)

Stata

## . margins i.gear
## 
## Predictive margins                                Number of obs   =         32
## Model VCE    : Robust
## 
## Expression   : Linear prediction, predict()
## 
## ------------------------------------------------------------------------------
##              |            Delta-method
##              |     Margin   Std. Err.      t    P>|t|     [95% Conf. Interval]
## -------------+----------------------------------------------------------------
##         gear |
##           3  |   .2810976    .108443     2.59   0.015     .0589621    .5032332
##           4  |   .6143126   .1681143     3.65   0.001      .269946    .9586792
##           5  |   .0649984   .0890096     0.73   0.471    -.1173296    .2473264
## ------------------------------------------------------------------------------

R

modmarg::marg(mod, var_interest = 'gear', type = 'levels',
              vcov = sandwich::vcovHC(mod, "HC1"),
              dof = mod$df.residual,
              weights = mtcars$wt)[[1]]
##      Label     Margin Standard.Error Test.Stat     P.Value Lower CI (95%)
## 1 gear = 3 0.28109763     0.10844305 2.5921222 0.014989089     0.05896212
## 2 gear = 4 0.61431261     0.16811433 3.6541359 0.001053408     0.26994601
## 3 gear = 5 0.06499843     0.08900964 0.7302403 0.471309415    -0.11732955
##   Upper CI (95%)
## 1      0.5032332
## 2      0.9586792
## 3      0.2473264

Margins (Effects)

Stata

## . margins, dydx(gear)
## 
## Average marginal effects                          Number of obs   =         32
## Model VCE    : Robust
## 
## Expression   : Linear prediction, predict()
## dy/dx w.r.t. : 4.gear 5.gear
## 
## ------------------------------------------------------------------------------
##              |            Delta-method
##              |      dy/dx   Std. Err.      t    P>|t|     [95% Conf. Interval]
## -------------+----------------------------------------------------------------
##         gear |
##           4  |    .333215   .2398524     1.39   0.176    -.1581005    .8245304
##           5  |  -.2160992   .1424822    -1.52   0.141    -.5079608    .0757624
## ------------------------------------------------------------------------------
## Note: dy/dx for factor levels is the discrete change from the base level.

R

modmarg::marg(mod, var_interest = 'gear', type = 'effects',
              vcov = sandwich::vcovHC(mod, "HC1"),
              dof = mod$df.residual,
              weights = mtcars$wt)[[1]]
##      Label     Margin Standard.Error Test.Stat   P.Value Lower CI (95%)
## 1 gear = 3  0.0000000      0.0000000       NaN       NaN      0.0000000
## 2 gear = 4  0.3332150      0.2398524  1.389250 0.1757072     -0.1581005
## 3 gear = 5 -0.2160992      0.1424822 -1.516675 0.1405588     -0.5079608
##   Upper CI (95%)
## 1     0.00000000
## 2     0.82453043
## 3     0.07576236

Binomial

See here

sandwich1 <- function(object, ...){
  sandwich::sandwich(object) * nobs(object) / (nobs(object) - 1)
}

Model Coefficients

Stata

## . logit vs i.gear mpg [pweight = wt]
## 
## Iteration 0:   log pseudolikelihood = -.65054522  
## Iteration 1:   log pseudolikelihood = -.31888591  
## Iteration 2:   log pseudolikelihood = -.28116092  
## Iteration 3:   log pseudolikelihood = -.27451436  
## Iteration 4:   log pseudolikelihood =   -.274456  
## Iteration 5:   log pseudolikelihood = -.27445597  
## 
## Logistic regression                               Number of obs   =         32
##                                                   Wald chi2(3)    =      14.02
##                                                   Prob > chi2     =     0.0029
## Log pseudolikelihood = -.27445597                 Pseudo R2       =     0.5781
## 
## ------------------------------------------------------------------------------
##              |               Robust
##           vs |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
## -------------+----------------------------------------------------------------
##         gear |
##           4  |   .8351264   1.563955     0.53   0.593    -2.230169    3.900422
##           5  |  -5.550662   2.211583    -2.51   0.012    -9.885285   -1.216038
##              |
##          mpg |   .6466109    .194564     3.32   0.001     .2652725    1.027949
##        _cons |  -12.91761   3.592946    -3.60   0.000    -19.95966    -5.87557
## ------------------------------------------------------------------------------

R

mod <- glm(vs ~ as.factor(gear) + mpg, data = mtcars, family = 'binomial',
           weights = wt)
## Warning: non-integer #successes in a binomial glm!
lmtest::coeftest(mod, vcov = sandwich1(mod),
                 df = mod$df.residual)
## 
## t test of coefficients:
## 
##                   Estimate Std. Error t value Pr(>|t|)   
## (Intercept)      -12.91762    3.59293 -3.5953 0.001229 **
## as.factor(gear)4   0.83513    1.56396  0.5340 0.597567   
## as.factor(gear)5  -5.55066    2.21158 -2.5098 0.018137 * 
## mpg                0.64661    0.19456  3.3234 0.002487 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Margins (Levels)

Stata

## . margins i.gear
## 
## Predictive margins                                Number of obs   =         32
## Model VCE    : Robust
## 
## Expression   : Pr(vs), predict()
## 
## ------------------------------------------------------------------------------
##              |            Delta-method
##              |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
## -------------+----------------------------------------------------------------
##         gear |
##           3  |   .3615676   .0729504     4.96   0.000     .2185875    .5045477
##           4  |   .4466534   .1465856     3.05   0.002      .159351    .7339558
##           5  |   .0742723   .0167933     4.42   0.000      .041358    .1071866
## ------------------------------------------------------------------------------

R

modmarg::marg(mod, var_interest = 'gear', type = 'levels',
              vcov = sandwich1(mod), dof = mod$df.residual,
              weights = mtcars$wt)[[1]]
##      Label     Margin Standard.Error Test.Stat      P.Value Lower CI (95%)
## 1 gear = 3 0.36156762     0.07295030  4.956356 7.182762e-07     0.21858767
## 2 gear = 4 0.44665335     0.14658547  3.047051 2.310988e-03     0.15935110
## 3 gear = 5 0.07427229     0.01679331  4.422731 9.746087e-06     0.04135801
##   Upper CI (95%)
## 1      0.5045476
## 2      0.7339556
## 3      0.1071866

Margins (Effects)

Stata

## . margins, dydx(gear)
## 
## Average marginal effects                          Number of obs   =         32
## Model VCE    : Robust
## 
## Expression   : Pr(vs), predict()
## dy/dx w.r.t. : 4.gear 5.gear
## 
## ------------------------------------------------------------------------------
##              |            Delta-method
##              |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
## -------------+----------------------------------------------------------------
##         gear |
##           4  |   .0850857   .1703424     0.50   0.617    -.2487792    .4189507
##           5  |  -.2872953   .0754265    -3.81   0.000    -.4351286   -.1394621
## ------------------------------------------------------------------------------
## Note: dy/dx for factor levels is the discrete change from the base level.

R

modmarg::marg(mod, var_interest = 'gear', type = 'effects',
              vcov = sandwich1(mod), dof = mod$df.residual,
              weights = mtcars$wt)[[1]]
##      Label      Margin Standard.Error  Test.Stat      P.Value
## 1 gear = 3  0.00000000     0.00000000        NaN          NaN
## 2 gear = 4  0.08508572     0.17034236  0.4994983 0.6174283719
## 3 gear = 5 -0.28729534     0.07542642 -3.8089484 0.0001395591
##   Lower CI (95%) Upper CI (95%)
## 1      0.0000000      0.0000000
## 2     -0.2487792      0.4189506
## 3     -0.4351284     -0.1394623