pweightsSet up data
rm(list=ls())
data(mtcars)
mtcars$wt <- mtcars$wt/sum(mtcars$wt)
# devtools::install_github('anniejw6/modmarg', ref = 'jacob_weights')
## . 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
## ------------------------------------------------------------------------------
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 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
## ------------------------------------------------------------------------------
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, 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.
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
See here
sandwich1 <- function(object, ...){
sandwich::sandwich(object) * nobs(object) / (nobs(object) - 1)
}
## . 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
## ------------------------------------------------------------------------------
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 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
## ------------------------------------------------------------------------------
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, 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.
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