--- title: "Introduction to predictive margins with svymargins package" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Introduction} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` # Predictive margins (PM) 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, 1. A regression analysis is performed, for example with three covariates: `age` (continuous), `sex` (factor) and `education` (factor). 2. The covariate `sex` is given value `"M"` (men) for all individuals. 3. Using the parameter estimates of the regression model and the covariate values, the predicted value is calculated for each individual. 4. The (weighted) mean of the individual predictive values is calculated ("predictive margin"). 5. Steps 2-4 are repeated, but using different values for the selected covariate, in this case by setting `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). ## Other software 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). # Load R packages ```{r, message=FALSE} library(survey) library(svymargins) library(dplyr) library(svyVGAM) # for multinomial logistic regression models ``` # Linear model ## Generate simulated data Generate the outcome variable `y` from normal distribution with three covariates: `age` (continuous), `sex` (factor) and `education` (factor). ```{r} 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) ``` ## Analysis with Stata ```{r, eval=FALSE, echo=FALSE} cmd <- "svyset [pw=one] svy: reg y i.education age i.sex margins sex education, grand // vce(unconditional) margins education, at(age=(40 50 60 70)) vce(unconditional)" # options("RStata.StataPath" =here::here("/usr/bin/stata-se")); options("RStata.StataVersion" = 16); stata(cmd, data.in = d |> mutate(one=1)) ``` An additional variable `one` with constant value 1 for all observations was added to the dataframe `d`. ```{text, eval=FALSE, echo=TRUE} . svyset [pw=one] Sampling weights: one VCE: linearized Single unit: missing Strata 1: Sampling unit 1: FPC 1: . 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 ------------------------------------------------------------------------------ ``` ## Regression analysis For simplicity, we use independent sampling design with weights equal to unity. The linear regression model is a main effects model. ```{r} # 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) ``` ## Predictive margins ### Define the desired margins In 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. * The package automatically searches for all observed levels of factor variables, so it is sufficient to give only the names of the covariates: `list("sex")`. * For numeric covariates, the desired covariate values must be given: e.g. `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`. ```{r} # Define margins as a named list: target.l <- list(null=list(), sex=list("sex"), educ=list("education"), age=list(age=seq(40,70,10)), educ_age=list("education", age=seq(40,70,10))) ``` ### Calculate the predictive margins After the regression analysis has been conducted and desired margins defined, the predictive margins can be calculated with the `svymargins` function: ```{r} # Calculate predictive margins: svymargins(res, groupfactor=target.l) # Get the output table containing the covariate information from "groups" attribute: attr(svymargins(res, groupfactor=target.l), "groups") ``` ### Confidence intervals ```{r} marg <- svymargins(res, groupfactor=target.l) cbind(marg, confint(marg)) marg <- svymargins(res, groupfactor=target.l) cbind(attr(marg, "groups"), confint(marg)) ``` Note 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()`. ### Contrasts using the `survey` package ```{r} marg <- svymargins(res, groupfactor=target.l) # Difference of adjusted means between high education and population: svycontrast(marg, quote(`educ,high` - `null`)) # Difference of adjusted means between high and low education: svycontrast(marg, quote(`educ,high` - `educ,low`)) # ... and confidence interval: svycontrast(marg, quote(`educ,high` - `educ,low`)) |> confint() # Ratio of adjusted means between high and low education: svycontrast(marg, quote(`educ,high` / `educ,low`)) ``` # Multinomial logistic regression model with `svyVGAM` package The 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. ```{r} 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) # 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) ``` ## Predictive margins Here the user needs to specify, for which outcome level (category) the adjusted prevalences are calculated using the `y.lev` argument. ```{r} # 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) # 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") ``` ### Confidence intervals #### Symmetric confidence intervals The usual `confint` method works also here. ```{r} marg <- svymargins(res, groupfactor=target.l, y.lev=1) cbind(marg, confint(marg)) marg <- svymargins(res, groupfactor=target.l, y.lev="B") cbind(attr(marg, "groups"), confint(marg)) ``` #### Asymmetric confidence intervals Probability 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. ```{r} marg <- svymargins(res, groupfactor=target.l, y.lev=1) cbind(marg, confint(marg) |> asymm_ci()) marg <- svymargins(res, groupfactor=target.l, y.lev="B") cbind(attr(marg, "groups"), confint(marg) |> asymm_ci()) ``` ### Contrasts using the `survey` package ```{r} marg <- svymargins(res, groupfactor=target.l, y.lev=2) # Adjusted risk (or prevalence) difference between high education and population: svycontrast(marg, quote(`educ,high` - `null`)) # Adjusted risk (or prevalence) difference between high and low education: svycontrast(marg, quote(`educ,high` - `educ,low`)) # ... and confidence interval: svycontrast(marg, quote(`educ,high` - `educ,low`)) |> confint() # Adjusted risk (or prevalence) ratio between high and low education: svycontrast(marg, quote(`educ,high` / `educ,low`)) # Population attributable fraction (PAF) of high education: svycontrast(marg, quote(1 - `educ,high` / `null`)) ``` # More formal presentation of the variance estimation Let $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}}. $$