Stablised weights

stabwt.R

expit<-function(eta) exp(eta)/(1+exp(eta))
logit<-function(p) log(p/(1-p))

one.sim<-function(beta, N=1000){
pop<-data.frame(x=rnorm(50000))
pop$y<-with(pop,rbinom(50000,1,expit(beta*x-3-beta)))
	
rval<-replicate(N,{	
cases<-which(pop$y==1)
controls<-sample(which(pop$y==0),length(cases))

thesample<-pop[c(cases,controls),]
thesample$wt<-with(thesample,ifelse(y==1,1,(50000-length(cases))/length(cases)))

m0<-glm(y~x,data=thesample,family=binomial)
m1<-glm(y~x,data=thesample,family=quasibinomial,weights=wt)
s1<-glm(I(1/wt)~x,data=thesample,family=Gamma)
thesample$stabwt<-thesample$wt/fitted(s1)
m2<-glm(y~x,data=thesample,family=quasibinomial,weights=stabwt)
s2<-glm(I(1/wt)~ns(x,3),data=thesample,family=Gamma)
thesample$stabwt2<-thesample$wt/fitted(s2)
m3<-glm(y~x,data=thesample,family=quasibinomial,weights=stabwt2)

c(coef(m0),coef(m1),coef(m2),coef(m3))})

c(apply(rval,1,median),apply(rval,1,mad))
}