1.

Draw a DAG of confounding. Indicate the confounder as ‘C’, the exposure as ‘X’, and the outcome as ‘Y’. Furthermore, name the arrows as follows: βYX is the arrow from X to Y, βXC is the arrow from C to X, and βYC is the arrow from C to Y.

require(DiagrammeR)
dag1 <- "digraph dag1 {
  graph [layout = dot]
  
  node[shape = circle]
  X; Y; C;
  
  X -> Y [label = 'Byx']
  C -> X [label = 'Bxc']
  C -> Y [label = 'Byc']
  
  subgraph {
    rank = same; X; Y
  }
}"

DiagrammeR(dag1, type = "grViz")

2.

Let’s make some assumptions about the relations between the different variables, i.e., make assumptions about the values of βYX, βXC, and βYC.

byx = 1.5
bxc = 1.1
byc = 1.3

3.

Now, let’s simulate some data to illustrate confounding, and confounding adjustment. The illustrated data will resemble the DAG you just drew. Let’s assume all data are normally distributed.

set.seed(12345)
C <- rnorm(100000, 0, 1)

X <- bxc * C + rnorm(100000, 0, 1)

Y <- byx * X + byc * C + rnorm(100000, 0, 1)

4.

Now, you can fit a linear model, regressing Y on X (‘unadjusted’) or regressing Y on X and C (‘adjusted’):

fit0 <- lm(Y~X)
fit1 <- lm(Y~X + C)

list(unadjusted = fit0, adjusted = fit1) %>% map_df(tidy, .id = "model")

       model        term      estimate   std.error    statistic   p.value
1 unadjusted (Intercept) -0.0002552413 0.004197753  -0.06080426 0.9515152
2 unadjusted           X  2.1483268085 0.002820318 761.73217488 0.0000000
3   adjusted (Intercept)  0.0015838455 0.003159298   0.50132823 0.6161412
4   adjusted           X  1.5026955136 0.003154558 476.35689813 0.0000000
5   adjusted           C  1.2985982410 0.004693741 276.66595192 0.0000000

5.

Look at the estimates of the effect of X on Y from the two linear models. How do these relate to the DAG?

The unadjusted model has overestimated the coefficient for X, since all lines are positive.

The adjusted model found the true values

6.

Take different values for βYX, βXC, and βYC (including zero, or negative values) and evaluate the impact in terms of discrepancy between the 2 models from step 4. Can you explain what you observe?

byx = 1.5
bxc = -1.1
byc = 1.3

set.seed(12345)
C <- rnorm(100000, 0, 1)

X <- bxc * C + rnorm(100000, 0, 1)
Y <- byx * X + byc * C + rnorm(100000, 0, 1)

fit0 <- lm(Y~X)
fit1 <- lm(Y~X + C)

list(unadjusted = fit0, adjusted = fit1) %>% map_df(tidy, .id = "model")

       model        term    estimate   std.error   statistic   p.value
1 unadjusted (Intercept) 0.004397018 0.004207348   1.0450806 0.2959882
2 unadjusted           X 0.854037887 0.002828755 301.9130022 0.0000000
3   adjusted (Intercept) 0.001583845 0.003159298   0.5013282 0.6161412
4   adjusted           X 1.502695514 0.003154558 476.3568981 0.0000000
5   adjusted           C 1.304528371 0.004690452 278.1242430 0.0000000

Now the effect of X on Y was underestimated in the unadjusted model, since one of the arrows was negative (and the other ones positive, so the sign of the product is negative)

1.

Draw a DAG of an exposure (‘X’), an outcome (‘Y’), and an intermediate (‘M’) of the relation between X and Y. Name the arrows as follows: βYX is the direct arrow from X to Y, βMX is the arrow from X to M, and βYM is the arrow from M to Y.

require(DiagrammeR)
dag2 <- "digraph dag2 {
  graph [layout = dot]
  
  node[shape = circle]
  X; M; Y;
  
  X -> Y [label = 'Byx']
  X -> M [label = 'Bmx']
  M -> Y [label = 'Bym']
  
  subgraph {
    rank = same; X; Y; M
  }
}"

DiagrammeR(dag2, type = "grViz")

2.

Make some assumptions about the relations between the different variables, i.e., make assumptions about the values of βYX, βMX, and βYM.

byx <- 1.5
bym <- 1.2
bmx <- 1.1

3. Simulate data that resemble the DAG you just drew. Again, we assume all data are normally distributed.

Sample X (e.g., 100,000 observations) from a normal distribution, with mean 0 and standard deviation 1: b. Define the relation between X and M (based on your DAG), e.g., βMX = 2: c. Generate M, based on X plus some random error: d. Define the relation between X and Y, and M and Y (based on your DAG), e.g., βYX = 1.0 and βYM = 0.5: e. Generate Y, based on X and M plus some random error:

set.seed(12345)

X <- rnorm(100000, 0, 1)
M <- bmx * X + rnorm(100000, 0, 1)
Y <- byx * X + bym * M + rnorm(100000, 0, 1)

4. Fit a linear model, regressing Y on X (‘unadjusted’) or regressing Y on X and M (‘adjusted’):

fit0 <- lm(Y~X)
fit1 <- lm(Y~X + M)

list(unadjusted = fit0, adjusted = fit1) %>% map_df(tidy, .id = "model")

       model        term    estimate   std.error   statistic   p.value
1 unadjusted (Intercept) 0.005904210 0.004948668   1.1930909 0.2328366
2 unadjusted           X 2.822411056 0.004947113 570.5167726 0.0000000
3   adjusted (Intercept) 0.001583845 0.003159298   0.5013282 0.6161412
4   adjusted           X 1.498598241 0.004693741 319.2758898 0.0000000
5   adjusted           M 1.202695514 0.003154558 381.2564149 0.0000000

5.

Look at the estimates of the effect of X on Y from the two linear models. How do these relate to the DAG?

In the unadjusted model, the direct effect of X on Y is overestimated, since it is partially medieated by M. In the adjusted model we get the true effects

6.

Take different values for βYX, βMX, and βYM (including zero, or negative values) and evaluate the impact in terms of discrepancy between the 2 models from step 4. Can you explain what you observe?

byx <- 1.5
bym <- -1.2
bmx <- 1.1

set.seed(12345)

X <- rnorm(100000, 0, 1)
M <- bmx * X + rnorm(100000, 0, 1)
Y <- byx * X + bym * M + rnorm(100000, 0, 1)

fit0 <- lm(Y~X)
fit1 <- lm(Y~X + M)

list(unadjusted = fit0, adjusted = fit1) %>% map_df(tidy, .id = "model")

       model        term     estimate   std.error    statistic
1 unadjusted (Intercept) -0.002717154 0.004935538   -0.5505284
2 unadjusted           X  0.180719356 0.004933988   36.6274428
3   adjusted (Intercept)  0.001583845 0.003159298    0.5013282
4   adjusted           X  1.498598241 0.004693741  319.2758898
5   adjusted           M -1.197304486 0.003154558 -379.5474506
        p.value
1  5.819583e-01
2 9.108075e-292
3  6.161412e-01
4  0.000000e+00
5  0.000000e+00

Just as before, the effect of X on Y is underestimated in the unadjusted model

1.

Draw a DAG of an exposure (‘X’), an outcome (‘Y’), and a common effect (‘S’). Name the arrows as follows: βYX is the arrow from X to Y, βSX is the arrow from X to S, and βSY is the arrow from Y to S.

require(DiagrammeR)
dag3 <- "digraph {
  graph [layout = dot]
  
  node[shape = circle]
  X; Y; S;
  
  X -> Y [label = 'Byx']
  X -> S [label = 'Bsx']
  Y -> S [label = 'Bsy']
  
  subgraph {
    rank = same; X; Y
  }
}"

DiagrammeR(dag3, type = "grViz")

2.

Make some assumptions about the relations between the different variables, i.e., make assumptions about the values of βYX, βSX, and βSY.

byx = 1.5
bsx = 1.4
bsy = 1.6

3.

Simulate some data based on the DAG you just drew. Let’s assume all data are normally distributed.

1. Sample X (e.g., 100,000 observations) from a normal distribution, with mean 0 and standard deviation 1:
2. Define the relation between X and Y (based on your DAG), e.g., βYX = 0.67:
3. Generate Y, based on X plus some random error:
4. Define the relation between X and S, and Y and S (based on your DAG), e.g., βSY = 2.0 and βSY = 1.0:
5. Generate S, based on X and Y plus some random error:

set.seed(12345)
X <- rnorm(100000, 0, 1)
Y <- byx * X + rnorm(100000, 0, 1)
S <- bsx * X + bsy * Y + rnorm(10000, 0, 1)

4.

Fit a linear model, regressing Y on X (‘unadjusted’) or regressing Y on X and S (‘adjusted’):

fit0 <- lm(Y~X)
fit1 <- lm(Y~X + S)

list(unadjusted = fit0, adjusted = fit1) %>% map_df(tidy, .id = "model")

       model        term      estimate    std.error   statistic   p.value
1 unadjusted (Intercept)  0.0035922351 0.0031670408   1.1342560 0.2566899
2 unadjusted           X  1.5007048746 0.0031660461 473.9996976 0.0000000
3   adjusted (Intercept) -0.0007969837 0.0016715409  -0.4767958 0.6335086
4   adjusted           X -0.2099862821 0.0037539097 -55.9380223 0.0000000
5   adjusted           S  0.4500441703 0.0008843321 508.9085305 0.0000000

5.

Look at the estimates of the effect of X on Y from the two linear models. How do these relate to the DAG?

Now the unadjusted model find the true value, since S is a collider of X and Y

In the adjusted model, byx even switched sign

6.

Take different values for βYX, βSX, and βSY (including zero, or negative values) and evaluate the impact in terms of discrepancy between the 2 models from step 4. Can you explain what you observe?

byx = 1.5
bsx = -1.1
bsy = 1.6

set.seed(12345)
X <- rnorm(100000, 0, 1)
Y <- byx * X + rnorm(100000, 0, 1)
S <- bsx * X + bsy * Y + rnorm(10000, 0, 1)

fit0 <- lm(Y~X)
fit1 <- lm(Y~X + S)

list(unadjusted = fit0, adjusted = fit1) %>% map_df(tidy, .id = "model")

       model        term      estimate    std.error   statistic   p.value
1 unadjusted (Intercept)  0.0035922351 0.0031670408   1.1342560 0.2566899
2 unadjusted           X  1.5007048746 0.0031660461 473.9996976 0.0000000
3   adjusted (Intercept) -0.0007969837 0.0016715409  -0.4767958 0.6335086
4   adjusted           X  0.9151241436 0.0020288515 451.0552617 0.0000000
5   adjusted           S  0.4500441703 0.0008843321 508.9085305 0.0000000

With negative bsx, byx did not switch sign in the adjusted model

byx = 1.5
bsx = 0
bsy = 1.6

set.seed(12345)
X <- rnorm(100000, 0, 1)
Y <- byx * X + rnorm(100000, 0, 1)
S <- bsx * X + bsy * Y + rnorm(10000, 0, 1)

fit0 <- lm(Y~X)
fit1 <- lm(Y~X + S)

list(unadjusted = fit0, adjusted = fit1) %>% map_df(tidy, .id = "model")

       model        term      estimate    std.error   statistic   p.value
1 unadjusted (Intercept)  0.0035922351 0.0031670408   1.1342560 0.2566899
2 unadjusted           X  1.5007048746 0.0031660461 473.9996976 0.0000000
3   adjusted (Intercept) -0.0007969837 0.0016715409  -0.4767958 0.6335086
4   adjusted           X  0.4200755563 0.0027020650 155.4646367 0.0000000
5   adjusted           S  0.4500441703 0.0008843321 508.9085305 0.0000000

With bsx = 0, S is no longer a collider. In the adjusted model, byx is too low. Part of the variation in X is now passed to S through Y. Including S in the model somehow sucks a part of the true byx into bsy, shrinking byx

1.

Create a table of the distributions of the confounding variables by exposure status. You may use commands such as: mean(age[vacc==0]); mean(age[vacc==1]) table(vacc,sex)

df %>%
  group_by(vacc) %>%
  summarize(mean_age = mean(age), fraction_female = mean(sex_female))

# A tibble: 2 x 3
  vacc  mean_age fraction_female
  <lgl>    <dbl>           <dbl>
1 F         75.1           0.648
2 T         75.9           0.606

tableone::CreateTableOne(data = df, strata = "vacc", test = F)

                       Stratified by vacc
                        FALSE         TRUE          
  n                     12030         32388         
  vacc = TRUE (%)           0 ( 0.0)  32388 (100.0) 
  age (mean (sd))       75.07 (7.37)  75.88 (6.86)  
  sex_female = TRUE (%)  7801 (64.8)  19643 ( 60.6) 
  cvd = TRUE (%)         4595 (38.2)  17032 ( 52.6) 
  cvd_drug = TRUE (%)    4360 (36.2)  16515 ( 51.0) 
  pulm = TRUE (%)         787 ( 6.5)   4631 ( 14.3) 
  pulm_drug = TRUE (%)    708 ( 5.9)   4283 ( 13.2) 
  DM = TRUE (%)           467 ( 3.9)   2398 (  7.4) 
  contact (mean (sd))   11.56 (10.69) 15.83 (11.87) 
  death = TRUE (%)        113 ( 0.9)    266 (  0.8) 

2.

What is the odds ratio of the crude (unadjusted) association between vaccination status and mortality?
fit <- glm(death ~ vacc, family=‘binomial’) log.or <- fit$coef[2] se.log.or <- sqrt(diag(vcov(fit))[2]) exp(c(log.or, log.or - 1.96se.log.or, log.or + 1.96se.log.or))

Let’s use confint to get likelihood profile confidence intervals

fit0 <- glm(death ~ vacc, family = 'binomial', data = df)
summary(fit0)

Call:
glm(formula = death ~ vacc, family = "binomial", data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.1374  -0.1374  -0.1284  -0.1284   3.0991  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.65833    0.09452 -49.286   <2e-16 ***
vaccTRUE    -0.13547    0.11280  -1.201     0.23    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 4365.8  on 44417  degrees of freedom
Residual deviance: 4364.3  on 44416  degrees of freedom
AIC: 4368.3

Number of Fisher Scoring iterations: 7

confint(fit0)

                 2.5 %      97.5 %
(Intercept) -4.8493808 -4.47853095
vaccTRUE    -0.3535467  0.08908945

No significant effect in unadjusted analysis

3.

A propensity score can be constructed by means of a multivariable regression model, predicting exposure status given confounders. What is the dependent variable (outcome variable) of this model? What are the independent variables of this model? What kind of regression model is most convenient for this study?

Logistic, since exposure is binary

Independent variables can be all covariates that are not the outcome (Death), and possibly interactions and non-linear transformations of these

4.

Based on your answer to the previous question construct a propensity score. For example:

ps_fit1 <- glm(vacc ~ age + sex_female, family = "binomial", data = df)

You can obtain the predicted value of the model (i.e., probability of the determinant, given the potential confounders): PS <- fit$fitted.values

Let’s keep the data nicely together in the data.frame

df %<>% 
  mutate(ps_1 = ps_fit1$fitted.values)

5.

What is the mean propensity score of those who are vaccinated? And what is the mean propensity score of those who are not vaccinated? Is this like you would expect?

df %>%
  group_by(vacc) %>%
  summarize(mean_ps = mean(ps_1))

# A tibble: 2 x 2
  vacc  mean_ps
  <lgl>   <dbl>
1 F       0.726
2 T       0.730

The propensity scores are pretty close together. Based on the baseline table we see that age and sex are distributed pretty equally amongst vaccinated and unvaccinated subjects, so it makes sense that the propensity score does not work too well

Let’s look at the distributions

df %>%
  ggplot(aes(x = ps_1, fill = vacc)) + 
  geom_histogram(alpha = .5)

6.

The primary goal of a propensity score is to balance confounder characteristics among those exposed and those unexposed to the determinant. Balance can be assessed e.g. within quintiles of the propensity score:

n.cat <- 5 # no. categories to split PS PS_cat <- ceiling(rank(PS)*n.cat/length(PS)) # split PS for (i in 1:max(PS_cat)){ print(sapply(split(age[PS_cat==i],vacc[PS_cat==i]),mean))}

n_cat <- 5
df %<>%
  mutate(ps_range = cut(ps_1, 
                      breaks = quantile(ps_1, 
                                        probs = seq(0, 1, length.out = n_cat + 1)),
                      include.lowest = T, ordered_result = T),
         ps_cat = as.numeric(ps_range))

df %>%
  group_by(ps_cat, ps_range) %>%
  summarize(mean(vacc), mean(age), mean(sex_female), mean(pulm))

# A tibble: 5 x 6
# Groups:   ps_cat [?]
  ps_cat ps_range   `mean(vacc)` `mean(age)` `mean(sex_femal… `mean(pulm)`
   <dbl> <ord>             <dbl>       <dbl>            <dbl>        <dbl>
1   1.00 [0.672,0.…        0.674        68.9            1.00         0.104
2   2.00 (0.7,0.72…        0.723        74.1            0.863        0.107
3   3.00 (0.723,0.…        0.737        74.3            0.490        0.117
4   4.00 (0.737,0.…        0.757        77.8            0.390        0.137
5   5.00 (0.756,0.…        0.763        84.4            0.275        0.150

7.

Donstruct a table of the distributions of the confounding variables by exposure status using only those subject whose PS lies within the third category of the PS. Compare the balance of confounding variables between exposure groups in this table with that in the table of the confounders by exposure status that was constructed using all subjects. In which table are the exposure groups more comparable?

tableone::CreateTableOne(data = df %>% filter(ps_cat == 3), strata = "vacc", test = F)

                       Stratified by vacc
                        FALSE          TRUE          
  n                      2223           6242         
  vacc = TRUE (%)           0 (  0.0)   6242 (100.0) 
  age (mean (sd))       73.94 (6.12)   74.45 (6.06)  
  sex_female = TRUE (%)  1024 ( 46.1)   3122 ( 50.0) 
  cvd = TRUE (%)          893 ( 40.2)   3310 ( 53.0) 
  cvd_drug = TRUE (%)     852 ( 38.3)   3220 ( 51.6) 
  pulm = TRUE (%)         141 (  6.3)    849 ( 13.6) 
  pulm_drug = TRUE (%)    126 (  5.7)    790 ( 12.7) 
  DM = TRUE (%)            94 (  4.2)    502 (  8.0) 
  contact (mean (sd))   11.15 (10.69)  15.14 (11.36) 
  death = TRUE (%)         16 (  0.7)     33 (  0.5) 
  ps_1 (mean (sd))       0.73 (0.00)    0.73 (0.00)  
  ps_range (%)                                       
     [0.672,0.7]            0 (  0.0)      0 (  0.0) 
     (0.7,0.723]            0 (  0.0)      0 (  0.0) 
     (0.723,0.737]       2223 (100.0)   6242 (100.0) 
     (0.737,0.756]          0 (  0.0)      0 (  0.0) 
     (0.756,0.836]          0 (  0.0)      0 (  0.0) 
  ps_cat (mean (sd))     3.00 (0.00)    3.00 (0.00)  

Not a whole lot better when compared to earlier

1.

There are different ways of using the PS to control for confounding. Here we will include the PS as a continuous covariate and as a categorical covariate in the model regressing outcome on exposure and PS:

fit <- glm(death ~ vacc + PS, family=‘binomial’) log.or <- fit$coef[2] se.log.or <- sqrt(diag(vcov(fit))[2]) exp(c(log.or, log.or - 1.96se.log.or, log.or + 1.96se.log.or))

fit <- glm(death ~ vacc + factor(PS_cat), family=‘binomial’) log.or <- fit$coef[2] se.log.or <- sqrt(diag(vcov(fit))[2]) exp(c(log.or, log.or - 1.96se.log.or, log.or + 1.96se.log.or))

fit_1 <-glm(death ~ vacc + ps_1, family = "binomial", data = df)
summary(fit_1)

Call:
glm(formula = death ~ vacc + ps_1, family = "binomial", data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.4210  -0.1484  -0.1166  -0.0883   3.5402  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -21.5561     1.3049 -16.519   <2e-16 ***
vaccTRUE     -0.2157     0.1132  -1.906   0.0567 .  
ps_1         22.9303     1.7395  13.182   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 4365.8  on 44417  degrees of freedom
Residual deviance: 4184.9  on 44415  degrees of freedom
AIC: 4190.9

Number of Fisher Scoring iterations: 8

confint(fit_1)

                  2.5 %       97.5 %
(Intercept) -24.1293723 -19.01318274
vaccTRUE     -0.4345128   0.00958604
ps_1         19.5331364  26.35332361

fit_2 <-glm(death ~ vacc + factor(ps_cat), family = "binomial", data = df)
summary(fit_2)

Call:
glm(formula = death ~ vacc + factor(ps_cat), family = "binomial", 
    data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.2162  -0.1418  -0.1043  -0.1018   3.4778  

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -5.8154     0.2123 -27.398  < 2e-16 ***
vaccTRUE         -0.2298     0.1133  -2.029 0.042472 *  
factor(ps_cat)2   0.7860     0.2454   3.202 0.001363 ** 
factor(ps_cat)3   0.8335     0.2463   3.384 0.000715 ***
factor(ps_cat)4   1.4502     0.2252   6.439  1.2e-10 ***
factor(ps_cat)5   2.0707     0.2160   9.589  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 4365.8  on 44417  degrees of freedom
Residual deviance: 4195.9  on 44412  degrees of freedom
AIC: 4207.9

Number of Fisher Scoring iterations: 8

confint(fit_2)

                     2.5 %       97.5 %
(Intercept)     -6.2567028 -5.421433710
vaccTRUE        -0.4487724 -0.004363489
factor(ps_cat)2  0.3153528  1.282339288
factor(ps_cat)3  0.3608170  1.331347980
factor(ps_cat)4  1.0257024  1.912630090
factor(ps_cat)5  1.6671111  2.517402559

2.

As a comparison, estimate the relation between vaccination status and mortality and adjust for confounding by including all confounders as covariates in the model.

fit_3 <- glm(reformulate(orig_vars, response = "death"), data = df, family = "binomial")
summary(fit_3)

Call:
glm(formula = reformulate(orig_vars, response = "death"), family = "binomial", 
    data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.7894  -0.1278  -0.0936  -0.0724   3.6549  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -10.753394   0.554101 -19.407  < 2e-16 ***
vaccTRUE        -0.487711   0.117305  -4.158 3.22e-05 ***
age              0.071497   0.006961  10.271  < 2e-16 ***
sex_femaleTRUE  -0.685423   0.109941  -6.234 4.53e-10 ***
cvdTRUE          0.662615   0.300484   2.205   0.0274 *  
cvd_drugTRUE    -0.582017   0.294435  -1.977   0.0481 *  
pulmTRUE         0.690273   0.369852   1.866   0.0620 .  
pulm_drugTRUE   -0.373648   0.382207  -0.978   0.3283    
DMTRUE           0.175138   0.172056   1.018   0.3087    
contact          0.049913   0.002742  18.204  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 4365.8  on 44417  degrees of freedom
Residual deviance: 3816.7  on 44408  degrees of freedom
AIC: 3836.7

Number of Fisher Scoring iterations: 8

The estimated effect of vaccination on death is greater when including all covariates

3.

How do all these odds ratios compare to the odds ratio that was estimated without adjustment for confounding?

fit_0 <- glm(death ~ vacc, data = df, family = "binomial")
summary(fit_0)

Call:
glm(formula = death ~ vacc, family = "binomial", data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.1374  -0.1374  -0.1284  -0.1284   3.0991  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.65833    0.09452 -49.286   <2e-16 ***
vaccTRUE    -0.13547    0.11280  -1.201     0.23    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 4365.8  on 44417  degrees of freedom
Residual deviance: 4364.3  on 44416  degrees of freedom
AIC: 4368.3

Number of Fisher Scoring iterations: 7

In the unadjusted model, the log-odds for vaccination is closer to the null value (smaller effect on death)

1.

What are potential advantages of PS analysis compared to including all confounders separately in a multivariable logistic regression model?

Checks during model building; possible to use more covariates in a setting with rare events but frequnet exposure

2.

How could the balance of confounding variables between exposure groups within strata of the PS possibly be improved?

By correlations between the covariates and the covariates used for the PS

3.

Could you think of ways of (graphically) summarizing the balance of confounders for subjects with the same PS?

Yes, for example calculate the standardized difference and plot these

First for all data

require(tableone)
diffs_overall <- CreateTableOne(vars = setdiff(orig_vars, "vacc"), strata = "vacc", data = df, test = F, smd = T) %>%
  ExtractSmd()

diffs_overall <- data.frame(variable = rownames(diffs_overall), 
                            smd = as.numeric(diffs_overall))

The in each PS quantile

diffs_quantiles <- split(df, df$ps_cat) %>%
  map_df(function(data) CreateTableOne(vars = setdiff(orig_vars, "vacc"), strata = "vacc", data = data, test = F, smd = T) %>%
  ExtractSmd())
diffs_melted <- melt(diffs_quantiles, variable.name = "PS_category", value.name = "smd")
diffs_melted$variable <- (diffs_overall$variable)

Now to plot

diffs_overall %>%
  ggplot(aes(x = smd, y = variable)) + 
  geom_point(shape = 5, size = 3) + 
  lims(x = c(0,.5)) +
  geom_point(data = diffs_melted, aes(x = smd, y = variable, col = PS_category), 
             alpha = .4)

We observe that the SMD within the PS categories (the colored dots), are not my closer to 0 than the overall. So the PS does not seem to equal out the distributions of covariates for accross vaccination status

Maybe this was because we used a very simple model for propensity scoring.

table(df$vacc)

FALSE  TRUE 
12030 32388 

We have 12030 in the smallest group of treatment, so we can potential include 1203 terms.

Let’s go all out and use all possible confounders, including up to cubic terms, and all third order interactions. Of course, the cubic terms of the binary variables are senseless, but this way we can get all terms with 2 lines of code

# ps_fit2 <- glm(vacc ~ (age + sex_female + cvd + cvd_drug + pulm + pulm_drug + DM + contact)^3,
#                data = df, family = "binomial")
ps_fit2 <- glm(vacc ~ polym(age, sex_female, cvd, cvd_drug, pulm, pulm_drug, DM, contact, 
                            degree = 3, raw = T),
               data = df, family = "binomial")

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

This takes some time to fit. Let’s see how many terms we included:

length(coef(ps_fit2))

[1] 165

Well below the threshold. Let’s use this as the new propensity score and do the same calculations

n_cat <- 5
df %<>%
  mutate(
    ps_2 = ps_fit2$fitted.values,
    ps_2_range = cut(ps_2, 
                      breaks = quantile(ps_2, 
                                        probs = seq(0, 1, length.out = n_cat + 1)),
                      include.lowest = T, ordered_result = T),
    ps_2_cat = as.numeric(ps_2_range))

diffs_quantiles_2 <- split(df, df$ps_2_cat) %>%
  map_df(function(data) CreateTableOne(vars = setdiff(orig_vars, "vacc"), strata = "vacc", data = data, test = F, smd = T) %>%
  ExtractSmd())
diffs_melted_2 <- melt(diffs_quantiles_2, variable.name = "PS_category", value.name = "smd")
diffs_melted_2$variable <- (diffs_overall$variable)

Now to plot

diffs_overall %>%
  ggplot(aes(x = smd, y = variable)) + 
  geom_point(shape = 5, size = 3) + 
  lims(x = c(0,.5)) +
  geom_point(data = diffs_melted_2, aes(x = smd, y = variable, col = PS_category), 
             alpha = .4)

Now we observe that the confounder distributions within the propensity score strata are a lot more alike. However, most covariates are binary, which are not very informative to plot.

df %>%
  select(age, contact, ps_2_cat, vacc) %>% 
  as.data.table() %>%
  melt.data.table(id.vars = c("vacc", "ps_2_cat")) %>%
  ggplot(aes(x = value, fill = vacc)) +
  geom_histogram(alpha = 0.8) + 
  facet_grid(ps_2_cat ~ variable, scales = "free")

Relative distributions seem equal within the propensity score quantiles

Additional ideas

• random forest for PS
• nicer functions for smd

1.

Look at the variables measured in this study (i.e., the code book), and plan your analysis (draw a DAG, define exposure, outcome, confounders, and the model to relate exposure to outcome, etc.)

require(DiagrammeR)
dag_iv <- "digraph {
  graph [layout = dot]
  
  node[shape = circle]
  program; birth_weight;
  age; eduction; height; weight; 
  n_prev_preg; low_birthweight;
  smoking_rand; smoking_8m;
  
  program -> smoking_8m
  smoking_8m -> birth_weight
  smoking_rand -> smoking_8m
  low_birthweight -> smoking_rand
  low_birthweight -> smoking_8m
  low_birthweight -> birth_weight
  n_prev_preg -> smoking_rand
  n_prev_preg -> birth_weight
  age -> birth_weight
  age -> smoking_rand
  age -> smoking_8m
  eduction -> smoking_rand
  eduction -> smoking_8m
  height -> smoking_rand
  height -> smoking_8m
  weight -> smoking_rand
  weight -> smoking_8m
  weight -> brith_weight

  subgraph {
    rank = same; program; smoking_8m; birth_weight
  }
}"
DiagrammeR(dag_iv, type = "grViz")

2.

Load the data and attach: data <- read.table(“data_IV.txt”, sep=“”) attach(data)

iv <- read.table(here("data", "data_IV.txt"), sep = "\t")

3.

Have a quick look: summary(data)

str(iv)

'data.frame':   935 obs. of  10 variables:
 $ program        : int  0 0 0 0 0 0 0 0 0 0 ...
 $ age            : int  20 25 24 16 12 18 30 31 27 36 ...
 $ education      : int  13 11 9 13 12 12 13 13 15 11 ...
 $ height         : int  163 173 152 175 151 165 167 173 173 167 ...
 $ weight         : int  55 67 52 66 51 55 62 76 67 68 ...
 $ N.prev.preg    : int  0 0 0 0 3 0 1 0 1 0 ...
 $ low.birthweight: int  0 0 0 0 1 0 0 0 0 0 ...
 $ smoking_rand   : int  12 17 15 12 15 12 12 11 10 14 ...
 $ smoking_8months: int  5 18 14 7 3 4 11 6 0 10 ...
 $ birth.weight   : int  2355 3703 3297 3387 3258 3962 3781 3160 3235 3395 ...

Let’s curate a bit

iv %<>%
  mutate(
    program = as.logical(program),
    low.birthweight = as.logical(low.birthweight)
  )

1.

What is the association between the number of cigarettes smoked per day at 8 months gestational age and birth weight? Explain what this number means. Do you think this is a valid estimate? Why?

fit0 <- lm(birth.weight ~ smoking_8months, data =iv)
summary(fit0)

Call:
lm(formula = birth.weight ~ smoking_8months, data = iv)

Residuals:
     Min       1Q   Median       3Q      Max 
-1434.61  -343.13    -4.61   315.17  1395.37 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     3287.611     24.467 134.368   <2e-16 ***
smoking_8months   -3.397      2.900  -1.171    0.242    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 510 on 933 degrees of freedom
Multiple R-squared:  0.001469,  Adjusted R-squared:  0.0003984 
F-statistic: 1.372 on 1 and 933 DF,  p-value: 0.2417

Smoking more reduces the birth weight, but it is not a significant effect.

It does not take into account confounding variables that affect both smoking and birth weight

1.

Create a table of the distributions of the potential confounding variables by program status. You may use commands such as:

iv %>%
  group_by(program) %>%
  summarize_all(funs(mean)) %>%
  t()

                        [,1]         [,2]
program            0.0000000    1.0000000
age               25.0550847   25.2937365
education         12.1250000   12.0043197
height           163.6906780  164.2419006
weight            61.8220339   62.6220302
N.prev.preg        1.2754237    1.2829374
low.birthweight    0.1588983    0.1835853
smoking_rand      11.8474576   11.9632829
smoking_8months    8.9449153    3.3455724
birth.weight    3209.1186441 3325.2850972

Which IV assumption can you check using this table?

Pretty equal confounder distributions, only smoking 8 months and birth weight seem to be different. (what we wanted)

2.

What is the association between the advice to stop smoking and the number of cigarettes smoked per day at 8 months gestational age? Explain what this number means. Which IV assumption do you check by this analysis?

lm(smoking_8months ~ program, data = iv) %>% summary()

Call:
lm(formula = smoking_8months ~ program, data = iv)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.9449 -3.3456 -0.9449  3.0551 16.6544 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.9449     0.2315   38.64   <2e-16 ***
programTRUE  -5.5993     0.3290  -17.02   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.03 on 933 degrees of freedom
Multiple R-squared:  0.2369,    Adjusted R-squared:  0.2361 
F-statistic: 289.7 on 1 and 933 DF,  p-value: < 2.2e-16

Difference in cigarettes, due to being in there program. Although the R-squared is low, so program does not explian much of the variance. You can check that the IV is related to the exposure of interest

In a plot:

iv %>%
  ggplot(aes(x = smoking_8months, fill = program)) + 
  geom_histogram(alpha = 0.5)

In a boxplot:

iv %>%
  ggplot(aes(y = smoking_8months, x = program)) + 
  geom_boxplot()

1. What is the association between the advice to stop smoking and birth weight? Explain what this number means.

lm(birth.weight ~ program, data = iv) %>% summary()

Call:
lm(formula = birth.weight ~ program, data = iv)

Residuals:
     Min       1Q   Median       3Q      Max 
-1472.29  -331.12     4.88   329.30  1456.88 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3209.12      23.34 137.489  < 2e-16 ***
programTRUE   116.17      33.17   3.502 0.000483 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 507.1 on 933 degrees of freedom
Multiple R-squared:  0.01298,   Adjusted R-squared:  0.01192 
F-statistic: 12.27 on 1 and 933 DF,  p-value: 0.0004832

This is the ‘intention to treat’ effect. So the effect of assigning someone to a program on the birth weight.

2.

Given the two effects of the advice to stop smoking (estimated in the previous two questions), can you estimate the (IV) effect of smoking on birth weight. Explain what this number means.

We have a model for smoking at 8 months (model iv, instrumental variable)

\[smoke_i = \alpha_0 + \alpha_1 * program_i + \nu_i\]

And a model for birthweight based on the program (model intention to treat, itt)

\[birthweight_i = \gamma_0 + \gamma_1 * program_i + \eta_i\]

What we want to know is

\[birthweight_i =\beta_0 + \beta_1*smoke_i+\epsilon_i\]

So we can use the nested models iv and itt

\[birthweight_i = \beta_0 + \beta_1*(\alpha_0 + \alpha_1 * program_i + \nu_i)+\epsilon_i\]

\[=\beta_0 + \beta_1*\alpha_0 + \beta_1*\alpha_1*program_i + \beta_1*\nu_i + \epsilon_i\]

Now we can easily see that

\[\beta_1 = \frac{\gamma_1}{\alpha_1}\]

fit_itt <- lm(birth.weight ~ program, data = iv)
fit_iv  <- lm(smoking_8months ~ program, data = iv)

coef(fit_itt)[2] / coef(fit_iv)[2]

programTRUE 
  -20.74644 

For each extra cigarette a day, the birthweight goes down with by 20 grams

3.

A faster way of estimating the (IV) effect of smoking on birth weight is:

fit <- lm(smoking_8months ~ program, data = iv) 
lm(birth.weight ~ predict(fit), data = iv) %>% summary()

Call:
lm(formula = birth.weight ~ predict(fit), data = iv)

Residuals:
     Min       1Q   Median       3Q      Max 
-1472.29  -331.12     4.88   329.30  1456.88 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3394.694     40.148  84.555  < 2e-16 ***
predict(fit)  -20.746      5.924  -3.502 0.000483 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 507.1 on 933 degrees of freedom
Multiple R-squared:  0.01298,   Adjusted R-squared:  0.01192 
F-statistic: 12.27 on 1 and 933 DF,  p-value: 0.0004832

Can you explain what is done here?

It directly nests the IV model into the eventual model (see neste model above)

1.

What proportion of women completely stopped smoking in among those enrolled in the program? And what was this proportion among those in the control group? stop <- 1*(smoking_8months==0) mean(stop[program==1]); mean(stop[program==0])

iv %>%
  group_by(program) %>%
  summarize(non_smoking_8months =mean(smoking_8months == 0))

# A tibble: 2 x 2
  program non_smoking_8months
  <lgl>                 <dbl>
1 F                     0.150
2 T                     0.428

2.

How much weight would a baby, on average, gain if a mother stopped smoking? Use the following expression: Or, using R: (mean(birth.weight[program==1])-mean(birth.weight[program==0])) / (mean(stop[program==1]) - mean(stop[program==0]))

We can calculate this directly by doing another IV analysis for totally stopping with smoking, just like above

iv %<>%
  mutate(non_smoking_8months = smoking_8months == 0)

lm(birth.weight ~ 
     lm(non_smoking_8months ~ program, data = iv)$fitted.values, 
   data = iv) %>%
  summary()

Call:
lm(formula = birth.weight ~ lm(non_smoking_8months ~ program, 
    data = iv)$fitted.values, data = iv)

Residuals:
     Min       1Q   Median       3Q      Max 
-1472.29  -331.12     4.88   329.30  1456.88 

Coefficients:
                                                           Estimate
(Intercept)                                                 3146.09
lm(non_smoking_8months ~ program, data = iv)$fitted.values   419.04
                                                           Std. Error
(Intercept)                                                     38.21
lm(non_smoking_8months ~ program, data = iv)$fitted.values     119.65
                                                           t value
(Intercept)                                                 82.338
lm(non_smoking_8months ~ program, data = iv)$fitted.values   3.502
                                                           Pr(>|t|)    
(Intercept)                                                 < 2e-16 ***
lm(non_smoking_8months ~ program, data = iv)$fitted.values 0.000483 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 507.1 on 933 degrees of freedom
Multiple R-squared:  0.01298,   Adjusted R-squared:  0.01192 
F-statistic: 12.27 on 1 and 933 DF,  p-value: 0.0004832

So 419

1.

Look at the variables measured in this study (i.e., the code book), and plan your analysis (draw a DAG, define exposure, outcome, confounders, and the model to relate exposure to outcome, etc.). Also think of possible unmeasured confounders.

For DAG, see above

Possible unmeasured confounders are immunodeficiency, smoking

2.

Load the data and attach: data <- read.table(“data_PS.txt”, sep=“”) attach(data)

3.

Have a quick look: summary(data)

As above

df <- read.table(here("data", "data_PS.txt"), sep = "\t")
str(df)

'data.frame':   44418 obs. of  10 variables:
 $ vacc     : int  1 1 1 1 1 1 1 1 1 1 ...
 $ age      : int  66 73 75 76 77 78 80 81 66 67 ...
 $ sex      : int  0 1 1 1 1 1 1 1 0 0 ...
 $ cvd      : int  0 1 1 1 1 1 1 1 0 0 ...
 $ cvd_drug : int  0 1 1 1 1 1 1 1 0 0 ...
 $ pulm     : int  1 0 0 0 0 0 0 0 0 0 ...
 $ pulm_drug: int  0 0 0 0 0 0 0 0 0 0 ...
 $ DM       : int  0 0 0 0 0 0 0 0 0 0 ...
 $ contact  : int  27 4 8 7 7 5 9 17 10 13 ...
 $ death    : int  0 0 0 0 0 0 0 0 0 0 ...

Let’s curate a few 0 - 1 variables as logical vectors (true vs false) so that R treats them right internally in all functions

logical_vars <- c("vacc", "sex", "cvd", "cvd_drug", "pulm", "pulm_drug", "DM", "death")
df  %<>%
  mutate_at(vars(logical_vars), as.logical)
str(df)

'data.frame':   44418 obs. of  10 variables:
 $ vacc     : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...
 $ age      : int  66 73 75 76 77 78 80 81 66 67 ...
 $ sex      : logi  FALSE TRUE TRUE TRUE TRUE TRUE ...
 $ cvd      : logi  FALSE TRUE TRUE TRUE TRUE TRUE ...
 $ cvd_drug : logi  FALSE TRUE TRUE TRUE TRUE TRUE ...
 $ pulm     : logi  TRUE FALSE FALSE FALSE FALSE FALSE ...
 $ pulm_drug: logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ DM       : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ contact  : int  27 4 8 7 7 5 9 17 10 13 ...
 $ death    : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...

Rename sex to a more sensible name

df %<>%
  rename(sex_female = sex)

Store names of original covariates in a vector for later use

orig_covariates <- setdiff(names(df), c("death", "vacc"))
orig_covariates

[1] "age"        "sex_female" "cvd"        "cvd_drug"   "pulm"      
[6] "pulm_drug"  "DM"         "contact"   

1.

Estimate the effect of influenza vaccination on the risk of mortality, while adjusting for measured confounders, e.g.: fit <- glm(death ~ vacc + age + sex, family=‘binomial’) log.or <- fit$coef[2] se.log.or <- sqrt(diag(vcov(fit))[2]) exp(c(log.or, log.or - 1.96se.log.or, log.or + 1.96se.log.or))

Let’s do this with profile likelihoods

fit <- glm(death ~ vacc + age + sex_female, family = "binomial", data = df)
summary(fit)

Call:
glm(formula = death ~ vacc + age + sex_female, family = "binomial", 
    data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.5134  -0.1431  -0.1129  -0.0898   3.5127  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -11.601509   0.539953 -21.486  < 2e-16 ***
vaccTRUE        -0.191353   0.113312  -1.689   0.0913 .  
age              0.093608   0.006676  14.021  < 2e-16 ***
sex_femaleTRUE  -0.552634   0.105524  -5.237 1.63e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 4365.8  on 44417  degrees of freedom
Residual deviance: 4165.7  on 44414  degrees of freedom
AIC: 4173.7

Number of Fisher Scoring iterations: 8

exp(confint(fit))

                      2.5 %       97.5 %
(Intercept)    3.160118e-06 2.625595e-05
vaccTRUE       6.633449e-01 1.034765e+00
age            1.083838e+00 1.112588e+00
sex_femaleTRUE 4.680495e-01 7.080727e-01

No significant effect of vaccination here

Step 3. Sensitivity analysis of unmeasured confounding:

We will apply the method proposed by Lin, Psaty and Kronmal (Biometrics 1998) in order to quantify the impact of a potential unmeasured confounder.

1. Consider an unmeasured binary confounder (e.g., smoking yes/no).

1. Define the prevalence of the confounder among vaccinated subject (p1).
2. Define the prevalence of the confounder among unvaccinated subjects (p0).
3. Define the odds ratio of the relation between the confounder and mortality (ORyz).

p0 <- .3
p1 <- .5
ORyz <- 2

2. Calculate the factor (A) that quantifies the impact of the unmeasured confounder:

A <- (ORyz*p1 + 1-p1 ) / (ORyz*p0 + 1-p0)
A

[1] 1.153846

3.

Estimate the effect of influenza vaccination on the risk of mortality, while adjusting for measured confounders and the unmeasured confounder: obs.OR <- exp(c(log.or, log.or - 1.96se.log.or, log.or + 1.96se.log.or)) names(obs.OR) <- c(‘OR’,‘95%’,‘CI’) obs.OR # observed odds ratio obs.OR/A # odds ratio after adjustment for unmeasured confounder

obs_OR <- c(exp(coef(fit)[2]), exp(confint(fit))[2,])
obs_OR

 vaccTRUE     2.5 %    97.5 % 
0.8258413 0.6633449 1.0347650 

obs_OR / A

 vaccTRUE     2.5 %    97.5 % 
0.7157291 0.5748989 0.8967963 

4.

Evaluate different scenarios of unmeasured confounding. You may want to do this in an automated way, varying two parameters, while keeping the third fixed. For example:

With the proposed code:

ORyz    <- 2
p0  <- seq(0.3,.5,.05)
p1  <- seq(0.1,.3,.05)

M <- matrix(ncol=length(p0)*3, nrow=length(p1))
colnames(M) <- rep(p0, each=3)
rownames(M) <- p1
for (i in 1:length(p0)){
    for (j in 1:length(p1)){
        A <- (ORyz*p1[j] + 1-p1[j] ) / (ORyz*p0[i] + 1-p0[i])
        M[j,(3*i-2):(3*i)] <- obs_OR/A  }}  
round(M,2)

      0.3  0.3  0.3 0.35 0.35 0.35  0.4  0.4  0.4 0.45 0.45 0.45  0.5  0.5
0.1  0.98 0.78 1.22 1.01 0.81 1.27 1.05 0.84 1.32 1.09 0.87 1.36 1.13 0.90
0.15 0.93 0.75 1.17 0.97 0.78 1.21 1.01 0.81 1.26 1.04 0.84 1.30 1.08 0.87
0.2  0.89 0.72 1.12 0.93 0.75 1.16 0.96 0.77 1.21 1.00 0.80 1.25 1.03 0.83
0.25 0.86 0.69 1.08 0.89 0.72 1.12 0.92 0.74 1.16 0.96 0.77 1.20 0.99 0.80
0.3  0.83 0.66 1.03 0.86 0.69 1.07 0.89 0.71 1.11 0.92 0.74 1.15 0.95 0.77
      0.5
0.1  1.41
0.15 1.35
0.2  1.29
0.25 1.24
0.3  1.19

Can you explain what is done here?

A grid of possible values for p0 and p1 are made, and for each combination the value for the adjusted OR is reported

Can you modify this procedure such that the parameters p1 and p0 are fixed, and ORyz is varied over a certain range?

Instead of using nested for-loops, we can map a function to a data.frame of parameters using the function pmap from purrr.

This is a little more R-ey than creating (nested) loops

We can just as easily do this for ranges of three parameters. Here the function expand.grid comes in handy

ORs <- seq(1, 4, length.out = 5)
p0s <- seq(0.3, 0.5, 0.05)
p1s <- seq(0.1, 0.3, 0.05)

param_grid <- expand.grid(OR = ORs, p0 = p0s, p1 = p1s)
head(param_grid)

    OR   p0  p1
1 1.00 0.30 0.1
2 1.75 0.30 0.1
3 2.50 0.30 0.1
4 3.25 0.30 0.1
5 4.00 0.30 0.1
6 1.00 0.35 0.1

dim(param_grid)

[1] 125   3

So we have 125 combinations of the 3 parameters

Now apply a function to this grid. A downside of pmap is that it depends on the column orders (instead of using named arguments), so you need to make sure you the column order of the parameter grid.

adjusted_ORs <- param_grid %>%
  pmap_df(function(OR, p0, p1) {
    A = (OR * p1 + 1 - p1) / (OR * p0 + 1-p0)
    adj_OR = obs_OR / A
    data.frame(estimate = adj_OR[1], ci_low = adj_OR[2], ci_high = adj_OR[3])
  })
adjusted_ORs

     estimate    ci_low  ci_high
1   0.8258413 0.6633449 1.034765
2   0.9410749 0.7559046 1.179151
3   1.0412781 0.8363913 1.304704
4   1.1292115 0.9070226 1.414883
5   1.2069988 0.9695040 1.512349
6   0.8258413 0.6633449 1.034765
7   0.9698833 0.7790445 1.215247
8   1.0951373 0.8796530 1.372188
9   1.2050541 0.9679420 1.509912
10  1.3022881 1.0460438 1.631745
11  0.8258413 0.6633449 1.034765
12  0.9986918 0.8021845 1.251344
13  1.1489965 0.9229146 1.439673
14  1.2808966 1.0288614 1.604942
15  1.3975775 1.1225836 1.751141
16  0.8258413 0.6633449 1.034765
17  1.0275002 0.8253244 1.287440
18  1.2028557 0.9661762 1.507158
19  1.3567392 1.0897808 1.699971
20  1.4928669 1.1991234 1.870537
21  0.8258413 0.6633449 1.034765
22  1.0563086 0.8484643 1.323537
23  1.2567150 1.0094378 1.574642
24  1.4325818 1.1507003 1.795000
25  1.5881563 1.2756632 1.989933
26  0.8258413 0.6633449 1.034765
27  0.9093533 0.7304247 1.139404
28  0.9775264 0.7851837 1.224824
29  1.0342311 0.8307309 1.295874
30  1.0821368 0.8692105 1.355899
31  0.8258413 0.6633449 1.034765
32  0.9371906 0.7527846 1.174284
33  1.0280881 0.8257967 1.288177
34  1.1036944 0.8865263 1.382910
35  1.1675687 0.9378324 1.462944
36  0.8258413 0.6633449 1.034765
37  0.9650280 0.7751445 1.209164
38  1.0786498 0.8664096 1.351530
39  1.1731577 0.9423217 1.469947
40  1.2530005 1.0064543 1.569988
41  0.8258413 0.6633449 1.034765
42  0.9928653 0.7975045 1.244043
43  1.1292115 0.9070226 1.414883
44  1.2426210 0.9981170 1.556983
45  1.3384324 1.0750761 1.677033
46  0.8258413 0.6633449 1.034765
47  1.0207027 0.8198644 1.278923
48  1.1797732 0.9476355 1.478236
49  1.3120842 1.0539124 1.644019
50  1.4238642 1.1436980 1.784078
51  0.8258413 0.6633449 1.034765
52  0.8797005 0.7066065 1.102250
53  0.9211306 0.7398846 1.154161
54  0.9539890 0.7662777 1.195332
55  0.9806865 0.7877220 1.228783
56  0.8258413 0.6633449 1.034765
57  0.9066301 0.7282373 1.135992
58  0.9687753 0.7781545 1.213859
59  1.0180629 0.8177441 1.275615
60  1.0581091 0.8499106 1.325793
61  0.8258413 0.6633449 1.034765
62  0.9335597 0.7498681 1.169734
63  1.0164200 0.8164244 1.273557
64  1.0821368 0.8692105 1.355899
65  1.1355317 0.9120992 1.422802
66  0.8258413 0.6633449 1.034765
67  0.9604893 0.7714989 1.203477
68  1.0640647 0.8546943 1.333255
69  1.1462107 0.9206769 1.436182
70  1.2129543 0.9742878 1.519811
71  0.8258413 0.6633449 1.034765
72  0.9874189 0.7931297 1.237219
73  1.1117094 0.8929642 1.392953
74  1.2102846 0.9721433 1.516466
75  1.2903770 1.0364763 1.616820
76  0.8258413 0.6633449 1.034765
77  0.8519205 0.6842926 1.067442
78  0.8708871 0.6995273 1.091207
79  0.8853018 0.7111057 1.109268
80  0.8966277 0.7202030 1.123459
81  0.8258413 0.6633449 1.034765
82  0.8779997 0.7052403 1.100119
83  0.9159330 0.7357097 1.147648
84  0.9447624 0.7588665 1.183771
85  0.9674140 0.7770611 1.212153
86  0.8258413 0.6633449 1.034765
87  0.9040789 0.7261881 1.132795
88  0.9609789 0.7718922 1.204090
89  1.0042230 0.8066273 1.258274
90  1.0382004 0.8339192 1.300847
91  0.8258413 0.6633449 1.034765
92  0.9301580 0.7471358 1.165472
93  1.0060248 0.8080746 1.260532
94  1.0636835 0.8543882 1.332777
95  1.1089868 0.8907774 1.389542
96  0.8258413 0.6633449 1.034765
97  0.9562372 0.7680835 1.198149
98  1.0510707 0.8442571 1.316974
99  1.1231441 0.9021490 1.407280
100 1.1797732 0.9476355 1.478236
101 0.8258413 0.6633449 1.034765
102 0.8258413 0.6633449 1.034765
103 0.8258413 0.6633449 1.034765
104 0.8258413 0.6633449 1.034765
105 0.8258413 0.6633449 1.034765
106 0.8258413 0.6633449 1.034765
107 0.8511221 0.6836513 1.066441
108 0.8685572 0.6976558 1.088287
109 0.8813082 0.7078979 1.104264
110 0.8910393 0.7157142 1.116457
111 0.8258413 0.6633449 1.034765
112 0.8764030 0.7039578 1.098118
113 0.9112731 0.7319667 1.141810
114 0.9367752 0.7524509 1.173763
115 0.9562372 0.7680835 1.198149
116 0.8258413 0.6633449 1.034765
117 0.9016838 0.7242643 1.129794
118 0.9539890 0.7662777 1.195332
119 0.9922421 0.7970039 1.243262
120 1.0214352 0.8204528 1.279841
121 0.8258413 0.6633449 1.034765
122 0.9269647 0.7445708 1.161471
123 0.9967050 0.8005886 1.248854
124 1.0477091 0.8415569 1.312762
125 1.0866332 0.8728222 1.361533

Combine with paramater grid

sa <- data.frame(param_grid, adjusted_ORs)
sa

      OR   p0   p1  estimate    ci_low  ci_high
1   1.00 0.30 0.10 0.8258413 0.6633449 1.034765
2   1.75 0.30 0.10 0.9410749 0.7559046 1.179151
3   2.50 0.30 0.10 1.0412781 0.8363913 1.304704
4   3.25 0.30 0.10 1.1292115 0.9070226 1.414883
5   4.00 0.30 0.10 1.2069988 0.9695040 1.512349
6   1.00 0.35 0.10 0.8258413 0.6633449 1.034765
7   1.75 0.35 0.10 0.9698833 0.7790445 1.215247
8   2.50 0.35 0.10 1.0951373 0.8796530 1.372188
9   3.25 0.35 0.10 1.2050541 0.9679420 1.509912
10  4.00 0.35 0.10 1.3022881 1.0460438 1.631745
11  1.00 0.40 0.10 0.8258413 0.6633449 1.034765
12  1.75 0.40 0.10 0.9986918 0.8021845 1.251344
13  2.50 0.40 0.10 1.1489965 0.9229146 1.439673
14  3.25 0.40 0.10 1.2808966 1.0288614 1.604942
15  4.00 0.40 0.10 1.3975775 1.1225836 1.751141
16  1.00 0.45 0.10 0.8258413 0.6633449 1.034765
17  1.75 0.45 0.10 1.0275002 0.8253244 1.287440
18  2.50 0.45 0.10 1.2028557 0.9661762 1.507158
19  3.25 0.45 0.10 1.3567392 1.0897808 1.699971
20  4.00 0.45 0.10 1.4928669 1.1991234 1.870537
21  1.00 0.50 0.10 0.8258413 0.6633449 1.034765
22  1.75 0.50 0.10 1.0563086 0.8484643 1.323537
23  2.50 0.50 0.10 1.2567150 1.0094378 1.574642
24  3.25 0.50 0.10 1.4325818 1.1507003 1.795000
25  4.00 0.50 0.10 1.5881563 1.2756632 1.989933
26  1.00 0.30 0.15 0.8258413 0.6633449 1.034765
27  1.75 0.30 0.15 0.9093533 0.7304247 1.139404
28  2.50 0.30 0.15 0.9775264 0.7851837 1.224824
29  3.25 0.30 0.15 1.0342311 0.8307309 1.295874
30  4.00 0.30 0.15 1.0821368 0.8692105 1.355899
31  1.00 0.35 0.15 0.8258413 0.6633449 1.034765
32  1.75 0.35 0.15 0.9371906 0.7527846 1.174284
33  2.50 0.35 0.15 1.0280881 0.8257967 1.288177
34  3.25 0.35 0.15 1.1036944 0.8865263 1.382910
35  4.00 0.35 0.15 1.1675687 0.9378324 1.462944
36  1.00 0.40 0.15 0.8258413 0.6633449 1.034765
37  1.75 0.40 0.15 0.9650280 0.7751445 1.209164
38  2.50 0.40 0.15 1.0786498 0.8664096 1.351530
39  3.25 0.40 0.15 1.1731577 0.9423217 1.469947
40  4.00 0.40 0.15 1.2530005 1.0064543 1.569988
41  1.00 0.45 0.15 0.8258413 0.6633449 1.034765
42  1.75 0.45 0.15 0.9928653 0.7975045 1.244043
43  2.50 0.45 0.15 1.1292115 0.9070226 1.414883
44  3.25 0.45 0.15 1.2426210 0.9981170 1.556983
45  4.00 0.45 0.15 1.3384324 1.0750761 1.677033
46  1.00 0.50 0.15 0.8258413 0.6633449 1.034765
47  1.75 0.50 0.15 1.0207027 0.8198644 1.278923
48  2.50 0.50 0.15 1.1797732 0.9476355 1.478236
49  3.25 0.50 0.15 1.3120842 1.0539124 1.644019
50  4.00 0.50 0.15 1.4238642 1.1436980 1.784078
51  1.00 0.30 0.20 0.8258413 0.6633449 1.034765
52  1.75 0.30 0.20 0.8797005 0.7066065 1.102250
53  2.50 0.30 0.20 0.9211306 0.7398846 1.154161
54  3.25 0.30 0.20 0.9539890 0.7662777 1.195332
55  4.00 0.30 0.20 0.9806865 0.7877220 1.228783
56  1.00 0.35 0.20 0.8258413 0.6633449 1.034765
57  1.75 0.35 0.20 0.9066301 0.7282373 1.135992
58  2.50 0.35 0.20 0.9687753 0.7781545 1.213859
59  3.25 0.35 0.20 1.0180629 0.8177441 1.275615
60  4.00 0.35 0.20 1.0581091 0.8499106 1.325793
61  1.00 0.40 0.20 0.8258413 0.6633449 1.034765
62  1.75 0.40 0.20 0.9335597 0.7498681 1.169734
63  2.50 0.40 0.20 1.0164200 0.8164244 1.273557
64  3.25 0.40 0.20 1.0821368 0.8692105 1.355899
65  4.00 0.40 0.20 1.1355317 0.9120992 1.422802
66  1.00 0.45 0.20 0.8258413 0.6633449 1.034765
67  1.75 0.45 0.20 0.9604893 0.7714989 1.203477
68  2.50 0.45 0.20 1.0640647 0.8546943 1.333255
69  3.25 0.45 0.20 1.1462107 0.9206769 1.436182
70  4.00 0.45 0.20 1.2129543 0.9742878 1.519811
71  1.00 0.50 0.20 0.8258413 0.6633449 1.034765
72  1.75 0.50 0.20 0.9874189 0.7931297 1.237219
73  2.50 0.50 0.20 1.1117094 0.8929642 1.392953
74  3.25 0.50 0.20 1.2102846 0.9721433 1.516466
75  4.00 0.50 0.20 1.2903770 1.0364763 1.616820
76  1.00 0.30 0.25 0.8258413 0.6633449 1.034765
77  1.75 0.30 0.25 0.8519205 0.6842926 1.067442
78  2.50 0.30 0.25 0.8708871 0.6995273 1.091207
79  3.25 0.30 0.25 0.8853018 0.7111057 1.109268
80  4.00 0.30 0.25 0.8966277 0.7202030 1.123459
81  1.00 0.35 0.25 0.8258413 0.6633449 1.034765
82  1.75 0.35 0.25 0.8779997 0.7052403 1.100119
83  2.50 0.35 0.25 0.9159330 0.7357097 1.147648
84  3.25 0.35 0.25 0.9447624 0.7588665 1.183771
85  4.00 0.35 0.25 0.9674140 0.7770611 1.212153
86  1.00 0.40 0.25 0.8258413 0.6633449 1.034765
87  1.75 0.40 0.25 0.9040789 0.7261881 1.132795
88  2.50 0.40 0.25 0.9609789 0.7718922 1.204090
89  3.25 0.40 0.25 1.0042230 0.8066273 1.258274
90  4.00 0.40 0.25 1.0382004 0.8339192 1.300847
91  1.00 0.45 0.25 0.8258413 0.6633449 1.034765
92  1.75 0.45 0.25 0.9301580 0.7471358 1.165472
93  2.50 0.45 0.25 1.0060248 0.8080746 1.260532
94  3.25 0.45 0.25 1.0636835 0.8543882 1.332777
95  4.00 0.45 0.25 1.1089868 0.8907774 1.389542
96  1.00 0.50 0.25 0.8258413 0.6633449 1.034765
97  1.75 0.50 0.25 0.9562372 0.7680835 1.198149
98  2.50 0.50 0.25 1.0510707 0.8442571 1.316974
99  3.25 0.50 0.25 1.1231441 0.9021490 1.407280
100 4.00 0.50 0.25 1.1797732 0.9476355 1.478236
101 1.00 0.30 0.30 0.8258413 0.6633449 1.034765
102 1.75 0.30 0.30 0.8258413 0.6633449 1.034765
103 2.50 0.30 0.30 0.8258413 0.6633449 1.034765
104 3.25 0.30 0.30 0.8258413 0.6633449 1.034765
105 4.00 0.30 0.30 0.8258413 0.6633449 1.034765
106 1.00 0.35 0.30 0.8258413 0.6633449 1.034765
107 1.75 0.35 0.30 0.8511221 0.6836513 1.066441
108 2.50 0.35 0.30 0.8685572 0.6976558 1.088287
109 3.25 0.35 0.30 0.8813082 0.7078979 1.104264
110 4.00 0.35 0.30 0.8910393 0.7157142 1.116457
111 1.00 0.40 0.30 0.8258413 0.6633449 1.034765
112 1.75 0.40 0.30 0.8764030 0.7039578 1.098118
113 2.50 0.40 0.30 0.9112731 0.7319667 1.141810
114 3.25 0.40 0.30 0.9367752 0.7524509 1.173763
115 4.00 0.40 0.30 0.9562372 0.7680835 1.198149
116 1.00 0.45 0.30 0.8258413 0.6633449 1.034765
117 1.75 0.45 0.30 0.9016838 0.7242643 1.129794
118 2.50 0.45 0.30 0.9539890 0.7662777 1.195332
119 3.25 0.45 0.30 0.9922421 0.7970039 1.243262
120 4.00 0.45 0.30 1.0214352 0.8204528 1.279841
121 1.00 0.50 0.30 0.8258413 0.6633449 1.034765
122 1.75 0.50 0.30 0.9269647 0.7445708 1.161471
123 2.50 0.50 0.30 0.9967050 0.8005886 1.248854
124 3.25 0.50 0.30 1.0477091 0.8415569 1.312762
125 4.00 0.50 0.30 1.0866332 0.8728222 1.361533

Of course these ‘3 dimensional’ data are harder to visualize, so we can grab all values for some p1

sa[sa$p0 == 0.5,]

      OR  p0   p1  estimate    ci_low  ci_high
21  1.00 0.5 0.10 0.8258413 0.6633449 1.034765
22  1.75 0.5 0.10 1.0563086 0.8484643 1.323537
23  2.50 0.5 0.10 1.2567150 1.0094378 1.574642
24  3.25 0.5 0.10 1.4325818 1.1507003 1.795000
25  4.00 0.5 0.10 1.5881563 1.2756632 1.989933
46  1.00 0.5 0.15 0.8258413 0.6633449 1.034765
47  1.75 0.5 0.15 1.0207027 0.8198644 1.278923
48  2.50 0.5 0.15 1.1797732 0.9476355 1.478236
49  3.25 0.5 0.15 1.3120842 1.0539124 1.644019
50  4.00 0.5 0.15 1.4238642 1.1436980 1.784078
71  1.00 0.5 0.20 0.8258413 0.6633449 1.034765
72  1.75 0.5 0.20 0.9874189 0.7931297 1.237219
73  2.50 0.5 0.20 1.1117094 0.8929642 1.392953
74  3.25 0.5 0.20 1.2102846 0.9721433 1.516466
75  4.00 0.5 0.20 1.2903770 1.0364763 1.616820
96  1.00 0.5 0.25 0.8258413 0.6633449 1.034765
97  1.75 0.5 0.25 0.9562372 0.7680835 1.198149
98  2.50 0.5 0.25 1.0510707 0.8442571 1.316974
99  3.25 0.5 0.25 1.1231441 0.9021490 1.407280
100 4.00 0.5 0.25 1.1797732 0.9476355 1.478236
121 1.00 0.5 0.30 0.8258413 0.6633449 1.034765
122 1.75 0.5 0.30 0.9269647 0.7445708 1.161471
123 2.50 0.5 0.30 0.9967050 0.8005886 1.248854
124 3.25 0.5 0.30 1.0477091 0.8415569 1.312762
125 4.00 0.5 0.30 1.0866332 0.8728222 1.361533

1.

In this exercise, you evaluated the potential impact of a single (binary) unmeasured confounder. What are limitations to this approach, and how could you overcome these limitations?

You can repeat this for other potential unmeasured confounders, try multivariate sensitivity analysis.

Or (monte-carlo) simulations for hypothesized distributions of continous potential unmeasured confounders

4.

Estimate the risk ratio and confidence interval between overweight and diastolic hypertension. How do you interpret the result? fit <- glm(bpdias_dich~bmi_dich, family=binomial(link=“log”)) exp(fit\(coef[2]); exp(fit\)coef[2] - 1.96sqrt(diag(vcov(fit))[2])) exp(fit$coef[2] + 1.96sqrt(diag(vcov(fit))[2]))

fit1 <- glm(hypertension ~ overweight, family = binomial(link = "log"), data = lrg)
exp(coef(fit1))

   (Intercept) overweightTRUE 
    0.07185869     2.49288557 

exp(confint(fit1))

                    2.5 %     97.5 %
(Intercept)    0.06215579 0.08244068
overweightTRUE 2.11830960 2.94744657

Looks like a significant effect for overweight

Instead , we can also write a short function, that we can use later:

Preferably we would use likelihood-profile confidence intervals, but they take time to calculate, so let’s use the wald approximation

extract.RR <- function(fit,q=2){
        A1 <- exp(fit$coef[q]); 
        A2 <- exp(fit$coef[q] - 1.96*sqrt(diag(vcov(fit))[q]))
        A3 <- exp(fit$coef[2] + 1.96*sqrt(diag(vcov(fit))[q]))
        result <- c(A1,A2,A3)
        return(result)}
extract.RR(fit1,2)

overweightTRUE overweightTRUE overweightTRUE 
      2.492886       2.113636       2.940184 

5.

Estimate the risk ratio and confidence interval between overweight and diastolic hypertension for subjects younger than 40 years and for subjects of 40 years or older. What do you conclude? fit <- glm(bpdias_dich~bmi_dich, subset=age_dich==0, family=binomial(link=“log”)) extract.RR(fit,2) fit <- glm(bpdias_dich~bmi_dich, subset=age_dich==1, family=binomial(link=“log”)) extract.RR(fit,2)

split(lrg, lrg$age_over_40) %>% 
  map_df(function(data) glm(hypertension ~ overweight, family = binomial(link = "log"), data = data) %>%
      extract.RR()) %>% t()

          [,1]     [,2]     [,3]
FALSE 2.531006 1.946935 3.290295
TRUE  1.854241 1.507600 2.280585

OR is higher in younger patients

6.

Recode the variables bmi_dich and age_dich into one variable with four categories. Recode young normal weight subjects as 1, older normal weight subjects as 2, young overweight subjects as 3, and older overweight subjects as 4. recode <- factor(ifelse(age_dich==0 & bmi_dich==0,1, ifelse(age_dich==1 & bmi_dich==0,2, ifelse(age_dich==0 & bmi_dich==1,3,4)))) table(recode,age_dich,bmi_dich)

lrg %<>%
  mutate(group = factor(paste0(overweight, age_over_40),
                         levels = c("FALSEFALSE", "FALSETRUE", "TRUEFALSE", "TRUETRUE"),
                         labels = 1:4))

lrg[!duplicated(lrg$group),] %>% select(overweight, age_over_40, group)

# A tibble: 4 x 3
  overweight age_over_40 group
  <lgl>      <lgl>       <fct>
1 T          F           3    
2 F          F           1    
3 F          T           2    
4 T          T           4    

7.

Estimate the risk ratio and confidence interval for diastolic hypertension for older normal weight subjects, young overweight subjects and older overweight subjects with young normal weight subjects as the reference category. Compare your results with the results of question 5. Do you see similarities? Do you see disparities? How could you explain these similarities and disparities? fit <- glm(bpdias_dich~factor(recode), family=binomial(link=“log”)) extract.RR(fit,2)
extract.RR(fit,3) extract.RR(fit,4)

fit2 <- glm(hypertension ~ group, data = lrg, family = binomial(link = "log"))
rrs <- map(as.list(2:4), function(x) extract.RR(fit2, x))
rrs

[[1]]
  group2   group2   group2 
3.364375 2.538789 4.458432 

[[2]]
  group3   group3   group2 
2.531006 1.946935 4.373671 

[[3]]
  group4   group4   group2 
6.238361 4.918018 4.267611 

Group 3 (overweight, younger than 40) matches the result from 5

Group 2 (no overweight, older than 40) does not match results from 5, due to different reference category

Group 4 is pretty high, so being old and overweight gives a high risk of hypertension

8.

Estimate from the results of question 7 the measure of multiplicative interaction. Is there interaction on a multiplicative scale?

rrs[[3]][1] / (rrs[[1]][1] * rrs[[2]][1])

   group4 
0.7326103 

Looks like a negative interaction on multiplicative scale

9.

Calculate from the results of question 7 the measure of additive interaction RERI. Is there interaction on an additive scale?

rrs[[3]][1] - rrs[[1]][1] - rrs[[2]][1] + 1

  group4 
1.342981 

Positive interaction on additive scale

10. Include bmi_dich, age_dich and the product term of the two in a log-binomial regression model. What are the risk ratios and confidence intervals for BMI, age and their product term? Interpret the results. Explain in words what the risk ratio of the product term means. fit <- glm(bpdias_dich~age_dich*bmi_dich, family=binomial(link=“log”)) extract.RR(fit,2) extract.RR(fit,3) extract.RR(fit,4)

Let’s update the function for extracting rate ratios

extract_RR <- function(fit){
        A1 <- exp(fit$coef); 
        A2 <- exp(fit$coef - 1.96*sqrt(diag(vcov(fit))))
        A3 <- exp(fit$coef + 1.96*sqrt(diag(vcov(fit))))
        return(data.frame(term = names(A1), estimate = A1, ci_low = A2, ci_high = A3))}

fit3 <- glm(hypertension~overweight*age_over_40, family = binomial("log"), data = lrg)
extract_RR(fit3)

                                                         term   estimate
(Intercept)                                       (Intercept) 0.04364641
overweightTRUE                                 overweightTRUE 2.53100578
age_over_40TRUE                               age_over_40TRUE 3.36437477
overweightTRUE:age_over_40TRUE overweightTRUE:age_over_40TRUE 0.73261032
                                   ci_low    ci_high
(Intercept)                    0.03518029 0.05414989
overweightTRUE                 1.94693465 3.29029546
age_over_40TRUE                2.53878869 4.45843232
overweightTRUE:age_over_40TRUE 0.52450358 1.02328733

The risk ratio of the product means that for age over 40, the relative risk for overweight is this times as great as in the younger group

11.

Compare the results of question 10 with your results of question 7 and 8. What are the similarities and what are the disparities? Is there significant interaction on a multiplicative scale?

The interaction-term is the same as answer 8, the result for overweight is the same as answer 7

Negative interaction on multiplicative scale, but not significant (ci overlaps 1)

12.

Use the excel spreadsheet “delta method” to calculate a confidence interval around your estimate of additive interaction as calculated in question 9. What additional information do you need to calculate the confidence interval? Is there significant interaction on an additive scale? (You can use either the “dummy” sheet or the “product term” sheet)

Use values from fit

summary(fit3)

Call:
glm(formula = hypertension ~ overweight * age_over_40, family = binomial("log"), 
    data = lrg)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.7973  -0.5636  -0.4839  -0.2988   2.5027  

Coefficients:
                               Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     -3.1316     0.1100 -28.465   <2e-16 ***
overweightTRUE                   0.9286     0.1339   6.937    4e-12 ***
age_over_40TRUE                  1.2132     0.1437   8.446   <2e-16 ***
overweightTRUE:age_over_40TRUE  -0.3111     0.1705  -1.825    0.068 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 3681.8  on 4896  degrees of freedom
Residual deviance: 3375.7  on 4893  degrees of freedom
AIC: 3383.7

Number of Fisher Scoring iterations: 6

We need the covariance matrix

vcov(fit3)

                               (Intercept) overweightTRUE age_over_40TRUE
(Intercept)                     0.01210393    -0.01210393     -0.01210393
overweightTRUE                 -0.01210393     0.01791782      0.01210393
age_over_40TRUE                -0.01210393     0.01210393      0.02063548
overweightTRUE:age_over_40TRUE  0.01210393    -0.01791782     -0.02063548
                               overweightTRUE:age_over_40TRUE
(Intercept)                                        0.01210393
overweightTRUE                                    -0.01791782
age_over_40TRUE                                   -0.02063548
overweightTRUE:age_over_40TRUE                     0.02906706

RERI: 1.3323, (0.37797, 2.2866) AP: 0.2144, (0.0694, 0.3593) S: 1.3431, (1.0647, 1.6943)

Significant interaction on additive scale

13.

What would be your overall conclusion of examining the interaction between age and overweight on diastolic hypertension? The joint effect of age and BMI on diastolic blood pressure is more than what would be expected if we sum up the individual effects of age and BMI, but less than if we multiply the individual effects.

4.

Is there significant interaction on a multiplicative scale? Do not answer this question with only yes or no, but also quantify the interaction. (Drawing a graph with the individual and joint effects may be useful)

ors <- extract_RR(fit2)
ors[3,2] / (ors[4,2]*ors[2,2])

[1] 0.732671

Looks like a negative multiplicative interaction, though not significant, judging on the confidence intervall of the interaction term

ors[-1,] %>%
  ggplot() + 
  geom_errorbarh(aes(x = estimate, xmin = ci_low, xmax = ci_high, y = term)) + 
  lims(x = c(0, 3))

5.

Is there significant interaction on an additive scale? Do not answer this question with only yes or no, but also quantify the interaction.

ANSWER: Calculated by hand: 1.28-1.23-1.43+1=-0.38. RERI is -0.38 (-1.23; 0.49), so negative interaction on additive scale but not significant because 0 is in the confidence interval. The combined effect is 0.38 less than the sum of the individual effects.

6. Does adjustment for potential confounders change the findings of question 4 and 5? (Assume that all measured covariables are potential confounders)

ANSWER: Multiplicative: Without covariables: OR product term is 0.81 (0.39; 1.67). With covariables: OR product term is 0.81 (0.38; 1.75). No influence of confounding.

Additive Without covariables: RERI is -0.25 (-1.14; 0.63). With covariables: RERI is -0.26 (-1.12; 0.60). No influence of confounding.