Usually the distribution of confounding variables like
age differ in the groups we wish to compare, thus an
adjustment is needed. In a linear regression model the regression
coefficient might reveal the adjusted difference of the groups directly.
However, in generalized linear models (GLM) like in logistic regression
models the parameters represent (log) odds ratios (OR), but an easier
interpretation can be based on risk ratios (RR) or risk differences
(RD), in which case the adjusted comparisons cannot be based on a single
regression parameter any more.
In predictive margins (Graubard and Korn 1999, Biometrics), some covariates are given a fixed value for all observations. For example,
age (continuous), sex (factor) and
education (factor).sex is given value "M" (men)
for all individuals.sex to value
"F" (women).The values of other covariates are retained at their observed values,
thus the (weighted) means are adjusted means, and we can
compare, in this example men and women, assuming that the
age and education distributions were the same
(as in the full data).
The SUDAAN software implemented the PREDMARG statement
in their statistical analysis procedures, accounting for the survey
sampling design and weights.
In the R statistical software, the survey package implements the same
functionality via the svypredmean function.
The limitations in these approaches were: (i) It was not possible to
choose the covariates to be fixed freely. E.g. if categorical age and
sex were main effects in the regression model, it was not possible to
define the margins for all combinations of the age and sex categories.
The researcher had to include the interaction term in
the regression model as well. (ii) It was not possible to restrict the
predictive margins to a subset of the full dataset. For
example, if the analysis was based on repeated cross-sectional surveys,
each representing the population in different years, then the predictive
margins represented the average of the populations at different years. A
more reasonable solution would be to restrict the predictive margins to
the latest year, so that the adjusted means (or prevalences) would
represent the population at that year. (iii) The R package
survey function svypredmean does not support
multinomial logistic regression models, which is
commonly used in analysing nominal outcomes.
In Stata, the procedure margins solved these
problems.
There is a need for a flexible solution also for R users, and this package aims to provide a solution to the limitations (i)–(iii).
Generate the outcome variable y from normal distribution
with three covariates: age (continuous), sex
(factor) and education (factor).
n <- 1000
# Generate data:
set.seed(1234)
d <- data.frame(sex=factor(sample(c("M", "F"), n, replace=TRUE)),
education=factor(sample(c("low", "middle", "high"), n, replace=TRUE)))
d <- d |> mutate(age=runif(n, 0, 40) + as.numeric(education) * 20,
y=rnorm(n, sd=5) + as.numeric(education) + 0.05 * age)
summary(d)
#> sex education age y
#> F:486 high :329 Min. :20.10 Min. :-11.170
#> M:514 low :338 1st Qu.:46.09 1st Qu.: 1.400
#> middle:333 Median :59.41 Median : 5.082
#> Mean :60.32 Mean : 4.952
#> 3rd Qu.:75.46 3rd Qu.: 8.495
#> Max. :99.90 Max. : 20.009An additional variable one with constant value 1 for all
observations was added to the dataframe d.
. svyset [pw=one]
Sampling weights: one
VCE: linearized
Single unit: missing
Strata 1: <one>
Sampling unit 1: <observations>
FPC 1: <zero>
. svy: reg y i.education age i.sex
(running regress on estimation sample)
Survey: Linear regression
Number of strata = 1 Number of obs = 1,000
Number of PSUs = 1,000 Population size = 1,000
Design df = 999
F(4, 996) = 40.25
Prob > F = 0.0000
R-squared = 0.1355
------------------------------------------------------------------------------
| Linearized
y | Coefficient std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
education |
low | 1.19902 .4801791 2.50 0.013 .2567451 2.141296
middle | 2.226229 .6605355 3.37 0.001 .9300326 3.522425
|
age | .0576193 .0136955 4.21 0.000 .030744 .0844946
|
sex |
M | .0579228 .3125908 0.19 0.853 -.5554872 .6713327
_cons | .3005586 .6491026 0.46 0.643 -.9732024 1.57432
------------------------------------------------------------------------------
. margins sex education, grand vce(unconditional)
Predictive margins
Number of strata = 1 Number of obs = 1,000
Number of PSUs = 1,000 Population size = 1,000
Design df = 999
Expression: Linear prediction, predict()
------------------------------------------------------------------------------
| Linearized
| Margin std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
sex |
F | 4.922536 .2305945 21.35 0.000 4.470031 5.375041
M | 4.980459 .2281417 21.83 0.000 4.532767 5.428151
|
education |
high | 3.805706 .3909045 9.74 0.000 3.038617 4.572794
low | 5.004726 .2794313 17.91 0.000 4.456386 5.553066
middle | 6.031934 .3733541 16.16 0.000 5.299286 6.764582
|
_cons | 4.952309 .1678844 29.50 0.000 4.622862 5.281755
------------------------------------------------------------------------------
. margins education, at(age=(40 50 60 70)) vce(unconditional)
Predictive margins
Number of strata = 1 Number of obs = 1,000
Number of PSUs = 1,000 Population size = 1,000
Design df = 999
Expression: Linear prediction, predict()
1._at: age = 40
2._at: age = 50
3._at: age = 60
4._at: age = 70
------------------------------------------------------------------------------
| Linearized
| Margin std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
_at#|
education |
1#high | 2.635102 .2769989 9.51 0.000 2.091535 3.178668
1#low | 3.834122 .396328 9.67 0.000 3.056391 4.611853
1#middle | 4.86133 .6045408 8.04 0.000 3.675015 6.047646
2#high | 3.211294 .3065697 10.47 0.000 2.6097 3.812889
2#low | 4.410315 .3139029 14.05 0.000 3.79433 5.026299
2#middle | 5.437523 .4837841 11.24 0.000 4.488173 6.386873
3#high | 3.787487 .3856877 9.82 0.000 3.030636 4.544338
3#low | 4.986507 .278402 17.91 0.000 4.440188 5.532827
3#middle | 6.013716 .3743501 16.06 0.000 5.279113 6.748318
4#high | 4.36368 .4909568 8.89 0.000 3.400255 5.327104
4#low | 5.5627 .306584 18.14 0.000 4.961078 6.164322
4#middle | 6.589908 .2893828 22.77 0.000 6.02204 7.157776
------------------------------------------------------------------------------For simplicity, we use independent sampling design with weights equal to unity. The linear regression model is a main effects model.
# Create survey design:
my.svy <- svydesign(~ 1, weights=~ 1, data=d)
# Run regression analysis:
res <- svyglm(y ~ education + age + sex, design=my.svy)
summary(res)
#>
#> Call:
#> svyglm(formula = y ~ education + age + sex, design = my.svy)
#>
#> Survey design:
#> svydesign(~1, weights = ~1, data = d)
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 0.30056 0.64910 0.463 0.64344
#> educationlow 1.19902 0.48018 2.497 0.01268 *
#> educationmiddle 2.22623 0.66054 3.370 0.00078 ***
#> age 0.05762 0.01370 4.207 2.82e-05 ***
#> sexM 0.05792 0.31259 0.185 0.85303
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for gaussian family taken to be 24.3666)
#>
#> Number of Fisher Scoring iterations: 2In this package, the margins are defined as a (names) list. Each
element produces a set of predictive margins, for example, the value of
the covariate sex is given value “M” (men) for all
individuals.
list("sex").list(age=seq(40,70,10)).The null margin is an empty list, and in this case all
covariates retain their observed values. Two (or more) covariates can be
set at the same time by including them in one entry of the list, for
example, list("education", age=seq(40,70,10)) calculates
the predictive margins for all combinations of education
and age.
After the regression analysis has been conducted and desired margins
defined, the predictive margins can be calculated with the
svymargins function:
# Calculate predictive margins:
svymargins(res, groupfactor=target.l)
#> mean SE
#> null 4.9523 0.1679
#> sex,F 4.9225 0.2321
#> sex,M 4.9805 0.2267
#> educ,high 3.8057 0.3904
#> educ,low 5.0047 0.2806
#> educ,middle 6.0319 0.3729
#> age,40 3.7817 0.3216
#> age,50 4.3579 0.2135
#> age,60 4.9341 0.1588
#> age,70 5.5103 0.2059
#> educ_age,high,40 2.6351 0.2770
#> educ_age,high,50 3.2113 0.3066
#> educ_age,high,60 3.7875 0.3857
#> educ_age,high,70 4.3637 0.4910
#> educ_age,low,40 3.8341 0.3963
#> educ_age,low,50 4.4103 0.3139
#> educ_age,low,60 4.9865 0.2784
#> educ_age,low,70 5.5627 0.3066
#> educ_age,middle,40 4.8613 0.6045
#> educ_age,middle,50 5.4375 0.4838
#> educ_age,middle,60 6.0137 0.3743
#> educ_age,middle,70 6.5899 0.2894
# Get the output table containing the covariate information from "groups" attribute:
attr(svymargins(res, groupfactor=target.l), "groups")
#> group_id sex education age mean SE
#> null null <NA> <NA> NA 4.952309 0.1678844
#> sex,F sex F <NA> NA 4.922536 0.2320666
#> sex,M sex M <NA> NA 4.980459 0.2267259
#> educ,high educ <NA> high NA 3.805706 0.3904234
#> educ,low educ <NA> low NA 5.004726 0.2806190
#> educ,middle educ <NA> middle NA 6.031934 0.3729471
#> age,40 age <NA> <NA> 40 3.781705 0.3216441
#> age,50 age <NA> <NA> 50 4.357897 0.2134977
#> age,60 age <NA> <NA> 60 4.934090 0.1588113
#> age,70 age <NA> <NA> 70 5.510283 0.2058501
#> educ_age,high,40 educ_age <NA> high 40 2.635102 0.2770388
#> educ_age,high,50 educ_age <NA> high 50 3.211294 0.3066048
#> educ_age,high,60 educ_age <NA> high 60 3.787487 0.3857149
#> educ_age,high,70 educ_age <NA> high 70 4.363680 0.4909774
#> educ_age,low,40 educ_age <NA> low 40 3.834122 0.3963062
#> educ_age,low,50 educ_age <NA> low 50 4.410315 0.3138745
#> educ_age,low,60 educ_age <NA> low 60 4.986507 0.2783688
#> educ_age,low,70 educ_age <NA> low 70 5.562700 0.3065529
#> educ_age,middle,40 educ_age <NA> middle 40 4.861330 0.6045403
#> educ_age,middle,50 educ_age <NA> middle 50 5.437523 0.4837828
#> educ_age,middle,60 educ_age <NA> middle 60 6.013716 0.3743476
#> educ_age,middle,70 educ_age <NA> middle 70 6.589908 0.2893785marg <- svymargins(res, groupfactor=target.l)
cbind(marg, confint(marg))
#> marg 2.5 % 97.5 %
#> null 4.952309 4.623261 5.281356
#> sex,F 4.922536 4.467694 5.377379
#> sex,M 4.980459 4.536084 5.424834
#> educ,high 3.805706 3.040490 4.570921
#> educ,low 5.004726 4.454723 5.554729
#> educ,middle 6.031934 5.300971 6.762897
#> age,40 3.781705 3.151294 4.412116
#> age,50 4.357897 3.939449 4.776345
#> age,60 4.934090 4.622825 5.245354
#> age,70 5.510283 5.106824 5.913741
#> educ_age,high,40 2.635102 2.092115 3.178088
#> educ_age,high,50 3.211294 2.610360 3.812229
#> educ_age,high,60 3.787487 3.031500 4.543474
#> educ_age,high,70 4.363680 3.401381 5.325978
#> educ_age,low,40 3.834122 3.057376 4.610868
#> educ_age,low,50 4.410315 3.795132 5.025497
#> educ_age,low,60 4.986507 4.440914 5.532100
#> educ_age,low,70 5.562700 4.961867 6.163533
#> educ_age,middle,40 4.861330 3.676453 6.046207
#> educ_age,middle,50 5.437523 4.489326 6.385720
#> educ_age,middle,60 6.013716 5.280008 6.747423
#> educ_age,middle,70 6.589908 6.022737 7.157080
marg <- svymargins(res, groupfactor=target.l)
cbind(attr(marg, "groups"), confint(marg))
#> group_id sex education age mean SE 2.5 %
#> null null <NA> <NA> NA 4.952309 0.1678844 4.623261
#> sex,F sex F <NA> NA 4.922536 0.2320666 4.467694
#> sex,M sex M <NA> NA 4.980459 0.2267259 4.536084
#> educ,high educ <NA> high NA 3.805706 0.3904234 3.040490
#> educ,low educ <NA> low NA 5.004726 0.2806190 4.454723
#> educ,middle educ <NA> middle NA 6.031934 0.3729471 5.300971
#> age,40 age <NA> <NA> 40 3.781705 0.3216441 3.151294
#> age,50 age <NA> <NA> 50 4.357897 0.2134977 3.939449
#> age,60 age <NA> <NA> 60 4.934090 0.1588113 4.622825
#> age,70 age <NA> <NA> 70 5.510283 0.2058501 5.106824
#> educ_age,high,40 educ_age <NA> high 40 2.635102 0.2770388 2.092115
#> educ_age,high,50 educ_age <NA> high 50 3.211294 0.3066048 2.610360
#> educ_age,high,60 educ_age <NA> high 60 3.787487 0.3857149 3.031500
#> educ_age,high,70 educ_age <NA> high 70 4.363680 0.4909774 3.401381
#> educ_age,low,40 educ_age <NA> low 40 3.834122 0.3963062 3.057376
#> educ_age,low,50 educ_age <NA> low 50 4.410315 0.3138745 3.795132
#> educ_age,low,60 educ_age <NA> low 60 4.986507 0.2783688 4.440914
#> educ_age,low,70 educ_age <NA> low 70 5.562700 0.3065529 4.961867
#> educ_age,middle,40 educ_age <NA> middle 40 4.861330 0.6045403 3.676453
#> educ_age,middle,50 educ_age <NA> middle 50 5.437523 0.4837828 4.489326
#> educ_age,middle,60 educ_age <NA> middle 60 6.013716 0.3743476 5.280008
#> educ_age,middle,70 educ_age <NA> middle 70 6.589908 0.2893785 6.022737
#> 97.5 %
#> null 5.281356
#> sex,F 5.377379
#> sex,M 5.424834
#> educ,high 4.570921
#> educ,low 5.554729
#> educ,middle 6.762897
#> age,40 4.412116
#> age,50 4.776345
#> age,60 5.245354
#> age,70 5.913741
#> educ_age,high,40 3.178088
#> educ_age,high,50 3.812229
#> educ_age,high,60 4.543474
#> educ_age,high,70 5.325978
#> educ_age,low,40 4.610868
#> educ_age,low,50 5.025497
#> educ_age,low,60 5.532100
#> educ_age,low,70 6.163533
#> educ_age,middle,40 6.046207
#> educ_age,middle,50 6.385720
#> educ_age,middle,60 6.747423
#> educ_age,middle,70 7.157080Note that for binary outcomes (values 0 or 1), the
asymm_ci function can be used to obtain asymmetric
confidence intervals, which stay within 0 and 1. For example,
confint(marg) |> asymm_ci().
survey packagemarg <- svymargins(res, groupfactor=target.l)
# Difference of adjusted means between high education and population:
svycontrast(marg, quote(`educ,high` - `null`))
#> nlcon SE
#> contrast -1.1466 0.3559
# Difference of adjusted means between high and low education:
svycontrast(marg, quote(`educ,high` - `educ,low`))
#> nlcon SE
#> contrast -1.199 0.4802
# ... and confidence interval:
svycontrast(marg, quote(`educ,high` - `educ,low`)) |> confint()
#> 2.5 % 97.5 %
#> contrast -2.140154 -0.2578867
# Ratio of adjusted means between high and low education:
svycontrast(marg, quote(`educ,high` / `educ,low`))
#> nlcon SE
#> contrast 0.76042 0.0888svyVGAM
packageThe simulated dataset has a three-category outcome variable
y with levels A, B and
C. The first level A is used as the reference
outcome level.
n <- 10000
# Generate data:
set.seed(1234)
d <- data.frame(sex=factor(sample(c("M", "F"), n, replace=TRUE)),
education=factor(sample(c("low", "middle", "high"), n, replace=TRUE)))
d <- d |>
mutate(age=runif(n, 0, 40) + as.numeric(education) * 20,
pr1=1,
pr2=exp(-1 + 0.5 * (as.numeric(education)-1) + 0.02 * age),
pr3=exp(1 + -0.5 * (as.numeric(education)-1) + 0.02 * age)) |>
mutate(across(matches("pr[0-9]"), ~ .x / (pr1 + pr2 + pr3))) |>
rowwise() |>
mutate(y=which(rmultinom(1, 1, c(pr1, pr2, pr3))[,1] == 1)) |>
ungroup() |>
mutate(y=factor(y, levels=1:3, labels=LETTERS[1:3]))
summary(d)
#> sex education age pr1 pr2
#> F:5037 high :3324 Min. :20.00 Min. :0.06339 Min. :0.09793
#> M:4963 low :3358 1st Qu.:44.89 1st Qu.:0.09323 1st Qu.:0.10647
#> middle:3318 Median :60.06 Median :0.11134 Median :0.23733
#> Mean :59.89 Mean :0.11447 Mean :0.26431
#> 3rd Qu.:74.83 3rd Qu.:0.13388 3rd Qu.:0.44491
#> Max. :99.99 Max. :0.17843 Max. :0.46831
#> pr3 y
#> Min. :0.4346 A:1134
#> 1st Qu.:0.4619 B:2610
#> Median :0.6452 C:6256
#> Mean :0.6212
#> 3rd Qu.:0.7477
#> Max. :0.8025
# Create survey design:
my.svy <- svydesign(~ 1, weights=~ 1, data=d)
# Run regression analysis:
res <- svy_vglm(y ~ education + age + sex,
family=multinomial(refLevel=1), design=my.svy)
summary(res)
#> svy_vglm.survey.design(y ~ education + age + sex, family = multinomial(refLevel = 1),
#> design = my.svy)
#> Independent Sampling design (with replacement)
#> svydesign(~1, weights = ~1, data = d)
#> Coef SE z p
#> (Intercept):1 -1.1176191 0.1444215 -7.7386 1.005e-14
#> (Intercept):2 0.9156547 0.1213080 7.5482 4.414e-14
#> educationlow:1 0.4637916 0.1139856 4.0689 4.724e-05
#> educationlow:2 -0.4635464 0.0952348 -4.8674 1.131e-06
#> educationmiddle:1 0.7548945 0.1553385 4.8597 1.176e-06
#> educationmiddle:2 -1.0799278 0.1369905 -7.8832 3.190e-15
#> age:1 0.0234133 0.0031483 7.4369 1.031e-13
#> age:2 0.0219073 0.0028058 7.8078 5.822e-15
#> sexM:1 0.0434120 0.0725970 0.5980 0.5498
#> sexM:2 0.0126961 0.0647859 0.1960 0.8446Here the user needs to specify, for which outcome level (category)
the adjusted prevalences are calculated using the y.lev
argument.
# Define margins as a named list:
target.l <- list(null=list(),
educ=list("education"),
age=list(age=seq(40,70,10)),
educ_age=list("education", age=seq(40,70,10)))
# Calculate predictive margins
# for the 1st outcome level "A":
svymargins(res, groupfactor=target.l, y.lev=1)
#> mean SE
#> null 0.113400 0.0032
#> educ,high 0.091826 0.0055
#> educ,low 0.118254 0.0057
#> educ,middle 0.150369 0.0109
#> age,40 0.164884 0.0088
#> age,50 0.136518 0.0051
#> age,60 0.112343 0.0033
#> age,70 0.091969 0.0035
#> educ_age,high,40 0.126719 0.0058
#> educ_age,high,50 0.104216 0.0058
#> educ_age,high,60 0.085317 0.0063
#> educ_age,high,70 0.069577 0.0067
#> educ_age,low,40 0.162486 0.0100
#> educ_age,low,50 0.134365 0.0068
#> educ_age,low,60 0.110465 0.0054
#> educ_age,low,70 0.090369 0.0052
#> educ_age,middle,40 0.205569 0.0193
#> educ_age,middle,50 0.171078 0.0135
#> educ_age,middle,60 0.141337 0.0092
#> educ_age,middle,70 0.116037 0.0064
# Get the output table containing the covariate information from "groups" attribute
# for the 2nd outcome level "B":
attr(svymargins(res, groupfactor=target.l, y.lev="B"), "groups")
#> y group_id education age mean SE
#> null B null <NA> NA 0.2610000 0.004392512
#> educ,high B educ high NA 0.1154941 0.007253262
#> educ,low B educ low NA 0.2373699 0.007335719
#> educ,middle B educ middle NA 0.4055084 0.013056277
#> age,40 B age <NA> 40 0.2336034 0.009937456
#> age,50 B age <NA> 50 0.2447548 0.007030187
#> age,60 B age <NA> 60 0.2547999 0.004806895
#> age,70 B age <NA> 70 0.2638098 0.004808076
#> educ_age,high,40 B educ_age high 40 0.1080408 0.005393824
#> educ_age,high,50 B educ_age high 50 0.1122954 0.005937567
#> educ_age,high,60 B educ_age high 60 0.1161833 0.007150073
#> educ_age,high,70 B educ_age high 70 0.1197449 0.008849008
#> educ_age,low,40 B educ_age low 40 0.2202772 0.010047307
#> educ_age,low,50 B educ_age low 50 0.2302085 0.008103387
#> educ_age,low,60 B educ_age low 60 0.2391883 0.007370292
#> educ_age,low,70 B educ_age low 70 0.2472939 0.008459925
#> educ_age,middle,40 B educ_age middle 40 0.3728347 0.021566596
#> educ_age,middle,50 B educ_age middle 50 0.3921307 0.017247027
#> educ_age,middle,60 B educ_age middle 60 0.4094232 0.013213021
#> educ_age,middle,70 B educ_age middle 70 0.4248075 0.009972701The usual confint method works also here.
marg <- svymargins(res, groupfactor=target.l, y.lev=1)
cbind(marg, confint(marg))
#> marg 2.5 % 97.5 %
#> null 0.11340000 0.10718469 0.11961531
#> educ,high 0.09182557 0.08109644 0.10255470
#> educ,low 0.11825412 0.10698483 0.12952340
#> educ,middle 0.15036914 0.12902514 0.17171315
#> age,40 0.16488358 0.14760961 0.18215755
#> age,50 0.13651763 0.12647416 0.14656109
#> age,60 0.11234293 0.10585138 0.11883449
#> age,70 0.09196903 0.08509722 0.09884083
#> educ_age,high,40 0.12671872 0.11540394 0.13803351
#> educ_age,high,50 0.10421594 0.09290416 0.11552772
#> educ_age,high,60 0.08531684 0.07303613 0.09759756
#> educ_age,high,70 0.06957717 0.05651751 0.08263682
#> educ_age,low,40 0.16248554 0.14285770 0.18211338
#> educ_age,low,50 0.13436518 0.12095138 0.14777898
#> educ_age,low,60 0.11046513 0.09989286 0.12103741
#> educ_age,low,70 0.09036875 0.08009188 0.10064563
#> educ_age,middle,40 0.20556851 0.16764929 0.24348774
#> educ_age,middle,50 0.17107759 0.14455597 0.19759922
#> educ_age,middle,60 0.14133710 0.12336725 0.15930695
#> educ_age,middle,70 0.11603709 0.10357141 0.12850278
marg <- svymargins(res, groupfactor=target.l, y.lev="B")
cbind(attr(marg, "groups"), confint(marg))
#> y group_id education age mean SE 2.5 %
#> null B null <NA> NA 0.2610000 0.004392512 0.25239083
#> educ,high B educ high NA 0.1154941 0.007253262 0.10127793
#> educ,low B educ low NA 0.2373699 0.007335719 0.22299214
#> educ,middle B educ middle NA 0.4055084 0.013056277 0.37991854
#> age,40 B age <NA> 40 0.2336034 0.009937456 0.21412636
#> age,50 B age <NA> 50 0.2447548 0.007030187 0.23097584
#> age,60 B age <NA> 60 0.2547999 0.004806895 0.24537852
#> age,70 B age <NA> 70 0.2638098 0.004808076 0.25438614
#> educ_age,high,40 B educ_age high 40 0.1080408 0.005393824 0.09746907
#> educ_age,high,50 B educ_age high 50 0.1122954 0.005937567 0.10065796
#> educ_age,high,60 B educ_age high 60 0.1161833 0.007150073 0.10216939
#> educ_age,high,70 B educ_age high 70 0.1197449 0.008849008 0.10240113
#> educ_age,low,40 B educ_age low 40 0.2202772 0.010047307 0.20058488
#> educ_age,low,50 B educ_age low 50 0.2302085 0.008103387 0.21432614
#> educ_age,low,60 B educ_age low 60 0.2391883 0.007370292 0.22474277
#> educ_age,low,70 B educ_age low 70 0.2472939 0.008459925 0.23071275
#> educ_age,middle,40 B educ_age middle 40 0.3728347 0.021566596 0.33056499
#> educ_age,middle,50 B educ_age middle 50 0.3921307 0.017247027 0.35832717
#> educ_age,middle,60 B educ_age middle 60 0.4094232 0.013213021 0.38352615
#> educ_age,middle,70 B educ_age middle 70 0.4248075 0.009972701 0.40526137
#> 97.5 %
#> null 0.2696092
#> educ,high 0.1297102
#> educ,low 0.2517476
#> educ,middle 0.4310982
#> age,40 0.2530805
#> age,50 0.2585337
#> age,60 0.2642212
#> age,70 0.2732335
#> educ_age,high,40 0.1186125
#> educ_age,high,50 0.1239328
#> educ_age,high,60 0.1301972
#> educ_age,high,70 0.1370886
#> educ_age,low,40 0.2399696
#> educ_age,low,50 0.2460908
#> educ_age,low,60 0.2536338
#> educ_age,low,70 0.2638750
#> educ_age,middle,40 0.4151045
#> educ_age,middle,50 0.4259343
#> educ_age,middle,60 0.4353202
#> educ_age,middle,70 0.4443536Probability scale parameters (i.e. with values between 0 and 1) can
be made asymmetric using the asymm_ci helper function. This
way the confidence intervals always remain between 0 and 1.
marg <- svymargins(res, groupfactor=target.l, y.lev=1)
cbind(marg, confint(marg) |> asymm_ci())
#> marg 2.5 % 97.5 %
#> null 0.11340000 0.10733166 0.11976541
#> educ,high 0.09182557 0.08164510 0.10313290
#> educ,low 0.11825412 0.10744147 0.12999644
#> educ,middle 0.15036914 0.13024726 0.17298140
#> age,40 0.16488358 0.14832734 0.18289098
#> age,50 0.13651763 0.12678156 0.14687560
#> age,60 0.11234293 0.10601334 0.11900013
#> age,70 0.09196903 0.08532408 0.09907543
#> educ_age,high,40 0.12671872 0.11582905 0.13847185
#> educ_age,high,50 0.10421594 0.09343440 0.11608228
#> educ_age,high,60 0.08531684 0.07381073 0.09842600
#> educ_age,high,70 0.06957717 0.05759826 0.08382577
#> educ_age,low,40 0.16248554 0.14379973 0.18308026
#> educ_age,low,50 0.13436518 0.12150779 0.14835333
#> educ_age,low,60 0.11046513 0.10032751 0.12148879
#> educ_age,low,70 0.09036875 0.08060463 0.10118549
#> educ_age,middle,40 0.20556851 0.17022423 0.24607513
#> educ_age,middle,50 0.17107759 0.14616130 0.19925021
#> educ_age,middle,60 0.14133710 0.12430299 0.16027827
#> educ_age,middle,70 0.11603709 0.10414115 0.12909609
marg <- svymargins(res, groupfactor=target.l, y.lev="B")
cbind(attr(marg, "groups"), confint(marg) |> asymm_ci())
#> y group_id education age mean SE 2.5 %
#> null B null <NA> NA 0.2610000 0.004392512 0.25248311
#> educ,high B educ high NA 0.1154941 0.007253262 0.10202061
#> educ,low B educ low NA 0.2373699 0.007335719 0.22329317
#> educ,middle B educ middle NA 0.4055084 0.013056277 0.38019620
#> age,40 B age <NA> 40 0.2336034 0.009937456 0.21469306
#> age,50 B age <NA> 50 0.2447548 0.007030187 0.23123925
#> age,60 B age <NA> 60 0.2547999 0.004806895 0.24549366
#> age,70 B age <NA> 70 0.2638098 0.004808076 0.25449472
#> educ_age,high,40 B educ_age high 40 0.1080408 0.005393824 0.09791463
#> educ_age,high,50 B educ_age high 50 0.1122954 0.005937567 0.10117396
#> educ_age,high,60 B educ_age high 60 0.1161833 0.007150073 0.10288652
#> educ_age,high,70 B educ_age high 70 0.1197449 0.008849008 0.10345697
#> educ_age,low,40 B educ_age low 40 0.2202772 0.010047307 0.20121704
#> educ_age,low,50 B educ_age low 50 0.2302085 0.008103387 0.21471123
#> educ_age,low,60 B educ_age low 60 0.2391883 0.007370292 0.22504305
#> educ_age,low,70 B educ_age low 70 0.2472939 0.008459925 0.23108827
#> educ_age,middle,40 B educ_age middle 40 0.3728347 0.021566596 0.33162450
#> educ_age,middle,50 B educ_age middle 50 0.3921307 0.017247027 0.35889134
#> educ_age,middle,60 B educ_age middle 60 0.4094232 0.013213021 0.38379922
#> educ_age,middle,70 B educ_age middle 70 0.4248075 0.009972701 0.40538852
#> 97.5 %
#> null 0.2697005
#> educ,high 0.1304883
#> educ,low 0.2520460
#> educ,middle 0.4313330
#> age,40 0.2536414
#> age,50 0.2587943
#> age,60 0.2643353
#> age,70 0.2733408
#> educ_age,high,40 0.1190759
#> educ_age,high,50 0.1244700
#> educ_age,high,60 0.1309476
#> educ_age,high,70 0.1382018
#> educ_age,low,40 0.2405991
#> educ_age,low,50 0.2464732
#> educ_age,low,60 0.2539313
#> educ_age,low,70 0.2642454
#> educ_age,middle,40 0.4159789
#> educ_age,middle,50 0.4264011
#> educ_age,middle,60 0.4355487
#> educ_age,middle,70 0.4444614survey packagemarg <- svymargins(res, groupfactor=target.l, y.lev=2)
# Adjusted risk (or prevalence) difference between high education and population:
svycontrast(marg, quote(`educ,high` - `null`))
#> nlcon SE
#> contrast -0.14551 0.0072
# Adjusted risk (or prevalence) difference between high and low education:
svycontrast(marg, quote(`educ,high` - `educ,low`))
#> nlcon SE
#> contrast -0.12188 0.0102
# ... and confidence interval:
svycontrast(marg, quote(`educ,high` - `educ,low`)) |> confint()
#> 2.5 % 97.5 %
#> contrast -0.1418093 -0.1019423
# Adjusted risk (or prevalence) ratio between high and low education:
svycontrast(marg, quote(`educ,high` / `educ,low`))
#> nlcon SE
#> contrast 0.48656 0.0337
# Population attributable fraction (PAF) of high education:
svycontrast(marg, quote(1 - `educ,high` / `null`))
#> nlcon SE
#> contrast 0.55749 0.0263Let \(y_i\) denote the outcome and \(X_i\) the vector of covariates of individual \(i=1,\ldots,n\). The corresponding vector of regression parameters of a generalized linear model (GLM) are \(\beta\). The linear predictor is \(X_i\beta\). The inverse link function of the GLM is \(\phi(\cdot)\), which maps the linear predictor into the expected value of the outcome for individual \(i\): \[ \mathbf{E}[Y_i\,|\,X_i,\,\beta]=\phi(X_i\beta). \] The predictive margin is defined as the (weighted) mean of the individual expected values: \[ \text{PM}(X,\,\beta):=\frac{\sum_{i=1}^nw_i\phi(X_i\beta)}{\sum_{i=1}^nw_i}, \] where \(w_i\ge0\) are the (often poststratification or inverse probability) weights. The covariate values can be modified as described above to yield the desired margins.
Let \(\hat{\beta}\) denote the point estimates and \(\hat{V}\) denote the estimated covariance matrix of the regression parameters of the GLM. The point estimate of the PM is calculated by replacing the parameters \(\beta\) with the estimates \(\hat{\beta}\).
The variance estimate is based on the total variance, following the R
code in the svypredmeans function of the
survey package: \[
\text{Var}\bigl\{\text{PM}(X,\,\hat{\beta})\bigr\} =
\text{Var}\left\{\mathbf{E}\bigl[\text{PM}(X,\,\hat{\beta})\,|\,X\bigr]\right\}
+
\mathbf{E}\left[\text{Var}\bigl\{\text{PM}(X,\,\hat{\beta})\,|\,X\bigr\}\right].
\] The first term can be calculated using the
svymeans function, which accounts for the weights and the
sampling design in calculating the variance of the expected values of
the outcome, which were first calculated using the predict
function.
The second term representing the uncertainty in the regression parameter estimates is calculated using the delta method, based on the gradient of \(\nabla\text{PM}=\partial\text{PM}(X_i,\,\beta)/\partial\beta\) (evaluated at the point estimates \(\hat{\beta}\)) and the covariance matrix \(\hat{V}\): \[ \nabla\text{PM}^T|_{\beta=\hat{\beta}}\,\hat{V}\,\nabla\text{PM}|_{\beta=\hat{\beta}}. \]