```
library(magrittr)dat <- read.table(header = TRUE, text = "
Y Z X count
1 1 1 16
0 1 1 25
1 0 1 230
0 0 1 5959
1 1 2 23
0 1 2 179
1 0 2 175
0 0 2 1367
1 1 3 32
0 1 3 225
1 0 3 49
0 0 3 220
")
xtabs(count ~ Y + Z + X, data = dat)[2:1,2:1,] %>% addmargins(., c(1))## , , X = 1
##
## Z
## Y 1 0
## 1 16 230
## 0 25 5959
## Sum 41 6189
##
## , , X = 2
##
## Z
## Y 1 0
## 1 23 175
## 0 179 1367
## Sum 202 1542
##
## , , X = 3
##
## Z
## Y 1 0
## 1 32 49
## 0 225 220
## Sum 257 269## Individual level data
dat2 <- epitools::expand.table(xtabs(count ~ Y + Z + X, data = dat))
xtabs( ~ Y + Z + X, data = dat2)## , , X = 1
##
## Z
## Y 0 1
## 0 5959 25
## 1 230 16
##
## , , X = 2
##
## Z
## Y 0 1
## 0 1367 179
## 1 175 23
##
## , , X = 3
##
## Z
## Y 0 1
## 0 220 225
## 1 49 32glm1 <- glm(formula = (Y == "1") ~ Z*X,
family = poisson(link = "log"),
data = dat2)
summary(glm1)##
## Call:
## glm(formula = (Y == "1") ~ Z * X, family = poisson(link = "log"),
## data = dat2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.8834 -0.2726 -0.2726 -0.2726 2.1585
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.29245 0.06594 -49.932 < 2e-16 ***
## Z1 2.35147 0.25855 9.095 < 2e-16 ***
## X2 1.11640 0.10031 11.129 < 2e-16 ***
## X3 1.58956 0.15734 10.103 < 2e-16 ***
## Z1:X2 -2.34819 0.34065 -6.893 5.45e-12 ***
## Z1:X3 -2.73192 0.34425 -7.936 2.09e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2923.6 on 8499 degrees of freedom
## Residual deviance: 2706.4 on 8494 degrees of freedom
## AIC: 3768.4
##
## Number of Fisher Scoring iterations: 6## exp(beta1)
exp(coef(glm1)[2])## Z1
## 10.50095## exp(beta1 + beta5)
exp(sum(coef(glm1)[c(2,5)]))## [1] 1.003281## exp(beta1 + beta6)
exp(sum(coef(glm1)[c(2,6)]))## [1] 0.6835544dat2Z0 <- dat2
dat2Z0$Z <- "0"
dat2Z1 <- dat2
dat2Z1$Z <- "1"
## Parametric g-formula
sum(exp(predict(glm1, newdata = dat2Z1))) / sum(exp(predict(glm1, newdata = dat2Z0)))## [1] 5.131315## Single estimate without the model
## Standardization to the population
R1 <- c(0.390, 0.114, 0.125)
R0 <- c(0.037, 0.113, 0.182)
wPopulation <- c(6230/8500, 1744/8500, 526/8500)
sum(R1 * wPopulation) / sum(R0 * wPopulation)## [1] 5.148469(0.390*(6230/8500) + 0.114*(1744/8500) + 0.125*(526/8500)) / (0.037*(6230/8500) + 0.113*(1744/8500) + 0.182*(526/8500))## [1] 5.148469## Create relative weights of exposed
wExposed <- c(41/500, 202/500, 257/500)
sum(R1 * wExposed) / sum(R0 * wExposed)## [1] 1.000366## Treatment model
txModel <- glm(formula = (Z == 1) ~ X,
family = binomial(link = "logit"),
data = dat2)
## PS
dat2$ps <- predict(txModel, type = "response")
## Weights
dat2$iptw <- 1 / ((dat2$ps * (dat2$Z == 1)) + ((1 - dat2$ps) * (1 - (dat2$Z == 1))))
dat2$smrw <- dat2$ps * dat2$iptw
dat2$srrw <- (1 - dat2$ps) * dat2$iptw
## SMR model
glm2 <- glm(formula = (Y == "1") ~ Z,
family = poisson(link = "log"),
data = dat2,
weight = smrw)
tableone::ShowRegTable(glm2)## exp(beta) [confint] p
## (Intercept) 0.14 [0.11, 0.18] <0.001
## Z1 1.00 [0.72, 1.38] 0.982## Create relative weights of exposed
wUnexposed <- c(6189/8000, 1542/8000, 269/8000)
sum(R1 * wUnexposed) / sum(R0 * wUnexposed)## [1] 5.800841## SRR model
glm3 <- glm(formula = (Y == "1") ~ Z,
family = poisson(link = "log"),
data = dat2,
weight = srrw)
tableone::ShowRegTable(glm3)## exp(beta) [confint] p
## (Intercept) 0.06 [0.05, 0.06] <0.001
## Z1 5.78 [5.24, 6.39] <0.001## IPTW model
glm4 <- glm(formula = (Y == "1") ~ Z,
family = poisson(link = "log"),
data = dat2,
weight = iptw)
tableone::ShowRegTable(glm4)## exp(beta) [confint] p
## (Intercept) 0.06 [0.06, 0.07] <0.001
## Z1 5.13 [4.68, 5.64] <0.001