#######################################################################
#  IPW 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


x <- cbind( const=1,
   DISTR.1=samp$DISTR.1, 
   BLACK=samp$BLACK,
   NBHISP=samp$NBHISP,
   GRADE=samp$GRADE,
   SLFHLTH=samp$SLFHLTH,
   SLFWGHT=samp$SLFWGHT,
   WORKHARD=samp$WORKHARD,
   GOODQUAL=samp$GOODQUAL,
   PHYSFIT=samp$PHYSFIT,
   PROUD=samp$PROUD,
   LIKESLF=samp$LIKESLF,
   ACCEPTED=samp$ACCEPTED,
   FEELLOVD=samp$FEELLOVD)

y <- samp$DISTR.2
t <- samp$DIET
xtwx <- t(x*pihat*(1-pihat))%*%x

# ACE and its standard error
mu1hat <- sum( t*y/pihat ) / sum(t/pihat)
mu0hat <- sum( (1-t)*y/(1-pihat) ) / sum((1-t)/(1-pihat))
ACE <- mu1hat - mu0hat
p <- ncol(x)
score <- matrix( NA, nrow(x), p+2 )
score[,1:p] <- (t-pihat)*x
score[,(p+1)] <- (1-t)*(1/(1-pihat))*(y-mu0hat)
score[,(p+2)] <- t*(1/pihat)*(y-mu1hat)
A <- matrix( 0, p+2, p+2 )
A[1:p,1:p] <- xtwx
A[(p+1),1:p] <- apply( -(1-t)*pihat*(1/(1-pihat))*(y-mu0hat)*x, 2, sum)
A[(p+2),1:p] <- apply( t*(1/pihat)*(1-pihat)*(y-mu1hat)*x, 2, sum)
A[(p+1),(p+1)] <- sum( (1-t)/(1-pihat) )
A[(p+2),(p+2)] <- sum( t/pihat )
Ainv <- solve(A)
cov.theta <- Ainv %*% t(score)%*%score %*% t(Ainv)
a <- c( rep(0,p), -1, 1)
SE <- sqrt( t(a)%*%cov.theta%*%a )

# > print(ACE)
# [1] -0.005064163
# > print(SE)
#            [,1]
# [1,] 0.02831887


# ACE.1 and its standard error
mu1hat <- sum( t*y ) / sum(t)
mu0hat <- sum( (1-t)*y*pihat/(1-pihat) ) / sum((1-t)*pihat/(1-pihat))
ACE.1 <- mu1hat - mu0hat
p <- ncol(x)
score[,(p+1)] <- (1-t)*pihat*(1/(1-pihat))*(y-mu0hat)
score[,(p+2)] <- t*(y-mu1hat)
A[(p+1),1:p] <- apply( -(1-t)*pihat*(1/(1-pihat))*(y-mu0hat)*x, 2, sum)
A[(p+2),1:p] <- 0
A[(p+1),(p+1)] <- sum( (1-t)*pihat/(1-pihat) )
A[(p+2),(p+2)] <- sum( t )
Ainv <- solve(A)
cov.theta <- Ainv %*% t(score)%*%score %*% t(Ainv)
a <- c( rep(0,p), -1, 1)
SE.1 <- sqrt( t(a)%*%cov.theta%*%a )

# > print(ACE.1)
# [1] -0.02528903
# > print(SE.1)
#            [,1]
# [1,] 0.01495289