---
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: <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
------------------------------------------------------------------------------
```

## 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}}.
$$

