# 
## G-estimation for general K-interval problem with linear blip
##      All-At-Once (Simultaneous not recursive) method 
## 	       (Sequentially randomized treatment) 
## 	   Written by Erica E M Moodie, 3 February, 2006 		 
#  
library(MASS)
expit <- function(x) { exp(x)/(1+exp(x)) }

#AAOg_estSeq.fcn <- function(PSI) {
AAOg_estSeq.fcn <- function(PSI,Y,status,action,return.all=F,tx.seq=F) {
   PSI0 <- PSI[1]
   PSI1 <- PSI[2]
   PSI22 <- PSI[3]
   PSI23 <- PSI[4]

   n <- length(Y)
   K <- 2
   Y.dash <- Y
   A <- action
   L <- status

   if(tx.seq) {
     p.coef <- coef(glm(A[,1]~L[,1],family=binomial)) 
     p.A1.given.L <- expit(cbind(1,L[,1])%*%p.coef)  
     p.coef <- coef(glm(A[,2]~L[,2]+A[,1],family=binomial)) 
     p.A2.given.L <- expit(cbind(1,L[,2],A[,1])%*%p.coef)  
   } else {
     p.A1.given.L <- mean(A[,1])
     p.A2.given.L <- mean(A[,2])
   }

   aa1 <- as.vector(A[,1] - p.A1.given.L)
   W1 <- cbind(1,L[,1])
   W1.2.inv <- ginv(t(W1)%*%W1/n)
   aa2 <- as.vector(A[,2] - p.A2.given.L)
   W2 <- cbind(1,L[,2],A[,1],L[,2]*A[,1])
   W2.2.inv <- ginv(t(W2)%*%W2/n)

   H1 <- Y.dash - A[,1]*(W1 %*% c(PSI0,PSI1)) +
       (as.numeric(cbind(1,L[,2],A[,1],L[,2]*A[,1]) %*% c(PSI0,PSI1,PSI22,PSI23) > 0) - 
          A[,2])*(W2 %*% c(PSI0,PSI1,PSI22,PSI23))
   H1.ctr <- H1 - W1 %*% (W1.2.inv %*% (t(W1) %*% H1/n)) 
   U1 <- cbind(as.vector(H1.ctr)*aa1*W1/n,0,0)

   H2 <- Y.dash - A[,2]*(W2 %*% c(PSI0,PSI1,PSI22,PSI23))
   H2.ctr <- H2 - W2 %*% (W2.2.inv %*% (t(W2) %*% H2/n)) 
   U2 <- as.vector(H2.ctr)*aa2*W2/n

   U <- (U1 + U2)/2
   if(return.all) return(apply(U,2,mean)^2)
   return(max(apply(U,2,mean)^2))
}


# nlm(AAOg_estSeq.fcn,p=c(-3,2,1,0),Y=Y,action=A,status=L)
# nlm(f=AAOg_estSeq.fcn,p=c(-3,2,1,0),iterlim=1000)