nri.R

expit<-function(eta) exp(eta)/(1+exp(eta))

one.nri<-function(N=1000){

x<-rnorm(50)
y<-rbinom(50,1,expit(x))
z<-rnorm(50)
m_old<-glm(y~x,family=binomial())
m_new<-glm(y~x+z,family=binomial())

x<-rnorm(N)
y<-rbinom(N,1,expit(x))
z<-rnorm(N)
p_old<-expit(cbind(1,x)%*%coef(m_old))
p_new<-expit(cbind(1,x,z)%*%coef(m_new))

mean((p_new>p_old)[y==1])-mean((p_new>p_old)[y==0])
}

r<-replicate(1000,one.nri())
summary(r)