alternative survival bias

survival_bias.r

library(dplyr)

# X causes DEATH
# P causes DEATH
# AGE causes DEATH
# G causes X
# AGE causes Y
# - IF cases and controls are AGE matched then NO SURVIVAL BIAS

n <- 1000000
age <- runif(n)
g <- rbinom(n, 2, 0.5)
x <- g + rnorm(n)
p <- rnorm(n)
o <- age + rnorm(n)


probdeath <- pnorm(scale(p + x + age))
death <- rbinom(n, 1, probdeath)
dat <- data.frame(a=age, g=g, x=x, p=p, o=o, d=death)


summary(lm(x ~ p, filter(dat, d == 1)))$coefficients[2,4]
summary(lm(x ~ p, filter(dat)))$coefficients[2,4]

summary(lm(g ~ p, filter(dat, d == 1)))$coefficients[2,4]
summary(lm(g ~ p, filter(dat)))$coefficients[2,4]

summary(lm(x ~ o, filter(dat, d == 1)))$coefficients[2,4]
summary(lm(x ~ o, filter(dat)))$coefficients[2,4]

summary(lm(x ~ a, filter(dat, d == 1)))$coefficients[2,4]
summary(lm(x ~ a, filter(dat)))$coefficients[2,4]

summary(lm(g ~ a, filter(dat, d == 1)))$coefficients[2,4]
summary(lm(g ~ a, filter(dat)))$coefficients[2,4]

summary(lm(x ~ o, filter(dat, d == 1)))$coefficients[2,4]
summary(lm(x ~ o, filter(dat)))$coefficients[2,4]

summary(lm(g ~ o, filter(dat, d == 1)))
summary(lm(g ~ o + a, filter(dat, d == 1)))


# X causes DEATH
# P causes DEATH
# AGE causes DEATH
# G causes X
# AGE and P cause Y
# - IF cases and controls are AGE matched then SURVIVAL STILL BIAS OCCURS


n <- 1000000
age <- runif(n)
g <- rbinom(n, 2, 0.5)
x <- g + rnorm(n)
p <- rnorm(n)
o <- age + p + rnorm(n)


probdeath <- pnorm(scale(p + x + age))
death <- rbinom(n, 1, probdeath)
dat <- data.frame(a=age, g=g, x=x, p=p, o=o, d=death)


summary(lm(x ~ p, filter(dat, d == 1)))$coefficients[2,4]
summary(lm(x ~ p, filter(dat)))$coefficients[2,4]

summary(lm(g ~ p, filter(dat, d == 1)))$coefficients[2,4]
summary(lm(g ~ p, filter(dat)))$coefficients[2,4]

summary(lm(x ~ o, filter(dat, d == 1)))$coefficients[2,4]
summary(lm(x ~ o, filter(dat)))$coefficients[2,4]

summary(lm(x ~ a, filter(dat, d == 1)))$coefficients[2,4]
summary(lm(x ~ a, filter(dat)))$coefficients[2,4]

summary(lm(g ~ a, filter(dat, d == 1)))$coefficients[2,4]
summary(lm(g ~ a, filter(dat)))$coefficients[2,4]

summary(lm(x ~ o, filter(dat, d == 1)))$coefficients[2,4]
summary(lm(x ~ o, filter(dat)))$coefficients[2,4]

summary(lm(g ~ o, filter(dat, d == 1)))
summary(lm(g ~ o + a, filter(dat, d == 1)))