#######################################################################
#  Propensity-subclassified estimates
#######################################################################


samp  <- read.table("Diet0001.dat", row.names="ID",
   col.names=c("ID", "DISTR.1", "BLACK", "NBHISP", "GRADE", "SLFHLTH",
   "SLFWGHT", "WORKHARD", "GOODQUAL", "PHYSFIT", "PROUD",
   "LIKESLF", "ACCEPTED", "FEELLOVD", "DIET", "Y.0", "Y.1", "DISTR.2",
   "PI.TRUE"))


pihat <- glm( DIET ~ DISTR.1 + BLACK + NBHISP + GRADE + 
   SLFHLTH + SLFWGHT + WORKHARD + GOODQUAL + PHYSFIT + PROUD + 
   LIKESLF + ACCEPTED + FEELLOVD, family=binomial, data=samp)$fitted


# define nine subclasses

stratum <- cut(pihat, 
   breaks=quantile(pihat, c(0,.20,.40,.60,.70,.80,.85,.90,.95,1)),
   include.lowest=T, labels=1:9)


# compute subclassified mean difference

t <- samp$DIET
y <- samp$DISTR.2
N <- nrow(samp)
N.s <- table(stratum)
Nstar <- sum(t)
Nstar.s <- table(stratum[t==1])
thetahat <- rep(NA,9)
vhat <- rep(NA,9)
for(s in 1:9){
   tmp <- lm( y[stratum==s] ~ t[stratum==s] )
   thetahat[s] <- tmp$coef[2]
   vhat[s] <- summary(tmp)$coef[2,2]^2}
ACE <- sum( (N.s/N) * thetahat )
SE <- sqrt( sum( (N.s/N)^2 * vhat ) )
ACE.1 <- sum( (Nstar.s/Nstar) * thetahat )
SE.1 <- sqrt( sum( (Nstar.s/Nstar)^2 * vhat ) )

# > print(ACE)
# [1] 0.005307301
# > print(SE)
# [1] 0.01868749
# > print(ACE.1)
# [1] -0.02375379
# > print(SE.1)
# [1] 0.01696347

# compute subclassified estimate w/ ANCOVA adjustment in classes

x <- samp$DISTR.1
for(s in 1:9){
   tmp <- lm( y[stratum==s] ~ t[stratum==s] + x[stratum==s])
   thetahat[s] <- tmp$coef[2]
   vhat[s] <- summary(tmp)$coef[2,2]^2}
ACE <- sum( (N.s/N) * thetahat )
SE <- sqrt( sum( (N.s/N)^2 * vhat ) )
ACE.1 <- sum( (Nstar.s/Nstar) * thetahat )
SE.1 <- sqrt( sum( (Nstar.s/Nstar)^2 * vhat ) )

# > print(ACE)
# [1] -0.003727728
# > print(SE)
# [1] 0.01547183
# > print(ACE.1)
# [1] -0.02349376
# > print(SE.1)
# [1] 0.01458467