logreg_bias.R

par(mfrow=c(1,2))
seeds = 10
x = seq(1,10,by=0.001)
all_u_coefs = seq(-5,5,by=0.1)

# lin reg
all_bs = rep(NA,length(all_u_coefs))
for ( i in 1:length(all_u_coefs) ) {
  u_coef = all_u_coefs[i]
  cur_r = rep(NA,seeds)
  for ( s in 1:seeds ) {
    u = rnorm(length(x))
    xb = -2 + 0.5*x + u_coef*u
    y = xb
    cur_r[s] = summary(lm(y ~ x))$coef[2,1]
  }
  all_bs[i] = mean(cur_r)
  #cat( r , '\n' )
}
plot(all_u_coefs,all_bs,type="l",col="royalblue",lwd=3,las=1,ylim=c(0,1),main="Linear Regression",xlab="Effect of uncorrelated confounder",ylab="Estimated effect of primary variable")
abline(h=0.5,lty=3,col="gray")

# log reg
all_bs = rep(NA,length(all_u_coefs))
for ( i in 1:length(all_u_coefs) ) {
  u_coef = all_u_coefs[i]
  cur_r = rep(NA,seeds)
  for ( s in 1:seeds ) {
    u = rnorm(length(x))
    xb = -2 + 0.5*x + u_coef*u
    p <- 1/(1 + exp(-xb))
    y = rbinom(n=length(p),size=1,prob=p)
    cur_r[s] = summary(glm(y ~ x, family="binomial" ))$coef[2,1]
  }
  all_bs[i] = mean(cur_r)
  #cat( r , '\n' )
}
plot(all_u_coefs,all_bs,type="l",col="royalblue",lwd=3,las=1,ylim=c(0,1),main="Logistic Regression",xlab="Effect of uncorrelated confounder",ylab="Estimated effect of primary variable")
abline(h=0.5,lty=3,col="gray")