So I am interested in calculating covariate-adjusted prevalence using survey data. My question pertains to a specific proprietary dataset, but below is an example using data(nhanes) data. I apologize in advance for the length of the post...
Step #1
I fit two regression models: (1) a crude model and (2) an adjusted model
#Setup
if(!require (survey)) install.packages("survey")
library(survey)
#Read in data
data("nhanes")
#Create survey design
d.obj = svydesign(id = ~SDMVPSU,
strata = ~SDMVSTRA,
weights = ~WTMEC2YR,
nest = TRUE,
data = nhanes)
#Regression models
mod.1 = svyglm(formula = HI_CHOL ~ as.factor(agecat),
design = d.obj,
family = quasibinomial(link = "log"))
mod.2 = svyglm(formula = HI_CHOL ~ as.factor(agecat) + as.factor(race) + as.factor(RIAGENDR),
design = d.obj,
family = quasibinomial(link = "log"))
Step #2
I exponentiated the coefficients to obtain the prevalence ratios and 95% CIs. The below summary shows the PRs for the crude model above the "NA" line and the PRs for the adjusted model below the line.
#Regression Summary
mod.sum = rbind(cbind("Estimate" = exp(coef(mod.1)), exp(confint(mod.1))),
c(NA, NA, NA),
cbind("Estimate" = exp(coef(mod.2)), exp(confint(mod.2))))
> mod.sum
Estimate 2.5 % 97.5 %
(Intercept) 0.008660272 0.004735950 0.01583638
as.factor(agecat)(19,39] 9.109573963 4.818424917 17.22229550
as.factor(agecat)(39,59] 20.610647336 10.441625642 40.68320376
as.factor(agecat)(59,Inf] 17.932147454 9.227318351 34.84890193
NA NA NA
(Intercept) 0.008597533 0.004623207 0.01598838
as.factor(agecat)(19,39] 9.051050139 4.816409068 17.00883531
as.factor(agecat)(39,59] 20.472879940 10.415264252 40.24274401
as.factor(agecat)(59,Inf] 17.641599353 9.050357148 34.38825923
as.factor(race)2 0.939905542 0.824583833 1.07135550
as.factor(race)3 0.693581683 0.536698646 0.89632339
as.factor(race)4 0.876065546 0.494248225 1.55284491
as.factor(RIAGENDR)2 1.210035838 1.048810926 1.39604451
Step #3
Prevalence estimation of HI_CHOL by age category is fairly simple...
#Calculting prevalence of HI_CHOL by agecat
prv = svyby(formula = ~HI_CHOL, by = ~agecat, design = d.obj, FUN = svymean, vartype = "ci", na.rm = TRUE)
> prv
agecat HI_CHOL ci_l ci_u
(0,19] (0,19] 0.008660267 0.003433241 0.01388729
(19,39] (19,39] 0.078891392 0.061116023 0.09666676
(39,59] (39,59] 0.178493821 0.156964219 0.20002342
(59,Inf] (59,Inf] 0.155297283 0.130664250 0.17993032
Using the point estimates from prv I can calculate by hand the point estimates from mod.1 if I use agecat == (0,19] as the reference group. However, I wish to calculate the covariate-adjusted prevalence for each age category (such that I could similarly calculate by hand the point estimates from mod.2 if I again use agecat == (0,19] as the reference group). I do not know how to do this.
In Stata, I can run the margins command to produce such estimates (or so I believe, but when I run the command on my proprietary data the lower bound of the confidence intervals are negative and I'm not sure how to interpret them).
. svyset SDMVPSU [pweight = WTMEC2YR], strata(SDMVSTRA) singleunit(centered)
pweight: WTMEC2YR
VCE: linearized
Single unit: centered
Strata 1: SDMVSTRA
SU 1: SDMVPSU
FPC 1: <zero>
. svy: glm HI_CHOL i.agecat i.race i.RIAGENDR, family(binomial) link(log) eform allbase
(running glm on estimation sample)
Survey: Generalized linear models
Number of strata = 15 Number of obs = 7,846
Number of PSUs = 31 Population size = 255,345,910
Design df = 16
------------------------------------------------------------------------------
| Linearized
HI_CHOL | exp(b) Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
agecat |
(0,19] | 1 (base)
(19,39] | 9.051065 2.917803 6.83 0.000 4.569875 17.92648
(39,59] | 20.47291 7.054992 8.76 0.000 9.86093 42.50513
(59,Inf] | 17.64163 5.994874 8.45 0.000 8.583895 36.25709
|
race |
1 | 1 (base)
2 | .9399055 .0616344 -0.95 0.359 .8179215 1.080082
3 | .6935817 .0908676 -2.79 0.013 .5253873 .9156208
4 | .8760655 .2516647 -0.46 0.651 .4764973 1.610693
|
RIAGENDR |
1 | 1 (base)
2 | 1.210036 .0872042 2.65 0.018 1.0386 1.409769
|
_cons | .0085975 .0027195 -15.04 0.000 .0043969 .0168111
------------------------------------------------------------------------------
. margins agecat, expression(exp(predict(xb))) vce(unconditional)
Predictive margins
Number of strata = 15 Number of obs = 7,846
Number of PSUs = 31 Population size = 255,345,910
Design df = 16
Expression : exp(predict(xb))
------------------------------------------------------------------------------
| Linearized
| Margin Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
agecat |
(0,19] | .0087239 .002677 3.26 0.005 .0030489 .0143989
(19,39] | .0789607 .009204 8.58 0.000 .059449 .0984724
(39,59] | .1786039 .0107848 16.56 0.000 .1557412 .2014666
(59,Inf] | .153904 .0129656 11.87 0.000 .1264183 .1813898
------------------------------------------------------------------------------
Additional Context
In my above nhanes example the point-estimates between the crude and adjusted models do not meaningfully differ. However, with my proprietary data, the analogous "crude" model is more of a "minimally-adjusted model" that takes the form outcome ~ exposure + covar where covar is a structural covariate that is not accounted for when setting the survey design. If I were to calculate outcome prevalence by exposure group traditionally without accounting for covar the resulting estimate would not be meaningful. Given this, the only way I know how to calculate prevalence would be based off a regression model (like I did in Stata), but I don't know how to do this in R and would really like to learn.
UPDATE 3/21/2021
So I believe I have solved the issue code-wise, but interpreting the results has me scratching my head. In short, I can group by exposure and sum the exponentiated model-fitted values for each permutation of covariates. I was able to do this using the nhanes as follows:
#Setup
library(survey)
library(magrittr)
library(dplyr)
library(tidyverse)
library(broom)
#Read in data
data("nhanes")
#Factoring covariates
nhanes$race %<>% factor()
nhanes$RIAGENDR %<>% factor()
#Create survey design
d.obj = svydesign(id = ~SDMVPSU,
strata = ~SDMVSTRA,
weights = ~WTMEC2YR,
nest = TRUE,
data = nhanes)
#Regression models
mod.1 = svyglm(formula = HI_CHOL ~ agecat,
design = d.obj,
family = quasibinomial(link = "log"))
mod.2 = svyglm(formula = HI_CHOL ~ agecat + race + RIAGENDR,
design = d.obj,
family = quasibinomial(link = "log"))
#Regression Summary
mod.sum =
list(Model.1 = data.frame(tidy(mod.1, exponentiate = TRUE, conf.int = TRUE)),
Model.2 = data.frame(tidy(mod.2, exponentiate = TRUE, conf.int = TRUE)))
#Appending regression fitted values to data frame
df.1 = augment(mod.1)
df.2 = augment(mod.2)
#Estimating probabilities for each permutation of regression formula
x.prob.1 = aggregate(cbind(est= exp(.fitted)) ~ agecat, data = df.1, FUN = mean)
x.prob.2 = aggregate(cbind(est= exp(.fitted)) ~ race + RIAGENDR + agecat, data = df.2, FUN = mean)
#Estimating probability of outcome given exposure, independant of covars
tab.x.1 = aggregate(est~ agecat, data = x.prob.1, FUN = sum)
tab.x.1 %<>% mutate(rr = est/est[1])
tab.x.2 = aggregate(est~ agecat, data = x.prob.2, FUN = sum)
tab.x.2 %<>% mutate(rr = est/est[1])
#Summarizing tab.x
tab.x.sum = list(Model.1 = tab.x.1, Model.2 = tab.x.2)
#Estimating prevalence through traditional method
prev = svyby(formula = ~HI_CHOL, by = ~agecat, design = d.obj, FUN = svymean, na.rm = TRUE)
The results of the summary objects are below...
> mod.sum
$Model.1
term estimate std.error statistic p.value conf.low conf.high
1 (Intercept) 0.008660272 0.3079463 -15.421547 9.802627e-10 0.00473595 0.01583638
2 agecat(19,39] 9.109573963 0.3249442 6.799094 1.264216e-05 4.81842492 17.22229550
3 agecat(39,59] 20.610647336 0.3469490 8.721189 8.585730e-07 10.44162564 40.68320376
4 agecat(59,Inf] 17.932147454 0.3389994 8.515045 1.122087e-06 9.22731835 34.84890193
$Model.2
term estimate std.error statistic p.value conf.low conf.high
1 (Intercept) 0.008597533 0.31652962 -15.0263347 1.111045e-07 0.004623207 0.01598838
2 agecat(19,39] 9.051050139 0.32186925 6.8440237 7.523606e-05 4.816409068 17.00883531
3 agecat(39,59] 20.472879940 0.34481686 8.7556654 1.068803e-05 10.415264252 40.24274401
4 agecat(59,Inf] 17.641599353 0.34054477 8.4284357 1.455588e-05 9.050357148 34.38825923
5 race2 0.939905542 0.06678723 -0.9279602 3.776433e-01 0.824583833 1.07135550
6 race3 0.693581683 0.13083519 -2.7965432 2.082994e-02 0.536698646 0.89632339
7 race4 0.876065546 0.29204774 -0.4530573 6.612323e-01 0.494248225 1.55284491
8 RIAGENDR2 1.210035838 0.07295691 2.6131862 2.812343e-02 1.048810926 1.39604451
> tab.x.sum
$Model.1
agecat est rr
1 (0,19] 0.008660272 1.000000
2 (19,39] 0.078891392 9.109574
3 (39,59] 0.178493821 20.610647
4 (59,Inf] 0.155297283 17.932147
$Model.2
agecat est rr
1 (0,19] 0.0666845 1.00000
2 (19,39] 0.6035648 9.05105
3 (39,59] 1.3652238 20.47288
4 (59,Inf] 1.1764213 17.64160
> prev
agecat HI_CHOL se
(0,19] (0,19] 0.008660267 0.002666899
(19,39] (19,39] 0.078891392 0.009069233
(39,59] (39,59] 0.178493821 0.010984693
(59,Inf] (59,Inf] 0.155297283 0.012568105
There are a few (if not many) problems with this approach. First, the log-binomial point estimates from R don't align with the ones from STATA, though this can probably be attributed to underling procedural differences. Second for R, while the prevalence estimates and their corresponding prevalence ratios add up, the estimates for the adjusted model don't make sense. For instance, tab.x.sum$Model.2 has prevalence estimates above 1.00! This leads me to believe that, while I have mathematically/code-wise solved my issue, the estimands I've produced are something entirely different than I thought they were.
Any insight into what I've actually calculated would be helpful...!