markmfredrickson/RItools — Issue #9: Randomization Dist of d2: Where should it live?
Author: jwbowers
Created: 2011-10-11 18:37:24+00:00
State: open

I have some rough code which generates the enumerated (or simulated) distribution of the d2 statistic. I am wondering where to put it (or perhaps Mark has some cool S4 idea for it). Here is the rough code, using the newspapers dataset (available separately):

```

covs.to.check<-c("rpre","blkpct","medhhi","medage","cands","pop")

xb1<-xBalance(as.formula(paste("z~",paste(covs.to.check,collapse="+"))),
data=news.df[,c("z","s",covs.to.check)],strata=list(s=~s),report="all")

## A very lazy way to get the enumerated distribution of the overall chisq test from xBalance
library(MASS) ##for the generalized inverse function b/d I'm lazy

Omega <- produceRandomizations(news.df$z, news.df$s, 1000) ##only gives 16 because that is the maximum.
Omega <- as.data.frame(apply(Omega, 2, function(idx) { x <- numeric(nrow(news.df)); x[idx] <- 1; return(x) }))

## For all z \in \Omega get adj.diffs. So, this is the distribution of the d statistic
d.dist<-sapply(Omega,function(thez){
thexb<-xBalance(as.formula(paste("z~",paste(covs.to.check,collapse="+"))),
data=data.frame(z=thez,news.df[,c("s",covs.to.check)]),
strata=list(s=~s),report="adj.mean.diffs")
thexb$results[,"adj.diff",]
})

## The distribution of the d2 statistic arises from the distribution of the d statistics
d2.dist<-apply(d.dist,2,function(thed){
t(thed) %*% ginv(cov(t(d.dist))) %*% thed
})

##Now simply compute the p-value
obs.d<-xb1$results[,"adj.diff",]
obs.d2<- as.numeric( t(obs.d) %*% ginv(cov(t(d.dist))) %*% obs.d)

two.p.d2<-2*min(mean(d2.dist>=obs.d2),mean(d2.dist<=obs.d2) )

print(two.p.d2)
```