x1 <- Age[r==1]
x2 <- CD4[r==1]
y <- Y[r==1]


if(alpha==0){
  aux <- gam((1-r)~s(Age,DF1)+s(CD4,DF2),family="binomial")
  p0 <- aux$fitted[r==1]
  aux <- predict.gam(aux,type="terms")
  h1 <- aux[r==1,1];h1 <- h1-mean(h1)
  h2 <- aux[r==1,2];h2 <- h2-mean(h2)
  STOP <- any(is.na(p0))
} else{
##starting values, etc..
  h1 <- trueh1[r==1]; h2 <- trueh2[r==1];eta <- trueeta

h11 <- gam((1-r)~s(Age,MAXDF),family="binomial")$add[r==1]
h21 <- gam((1-r)~s(CD4,MAXDF),family="binomial")$add[r==1]
if(any(abs(h11)>maxh) | any(abs(h21)>maxh)){
  STOP <- T
  cat("\n\nNO ITERATION, stopped at 1\n\n")
  Flag <- F
} else  Flag <- T
Count <- 0
while(Flag && Count < Maxcount){ 
  Count <- Count + 1

  h1prev <- h1
  aux <- h11 - (h2+eta+alpha*Q(y))
  if(any(is.na(aux))) STOP2 <- T
  else{
    h1 <- gam(aux~s(x1,df=DF1))$fitted
    auxeta <- mean(h1*(1+exp(eta+h1+h2+alpha*Q(y))))
    h1 <- h1 - auxeta
    eta <- eta + auxeta
    
    h2prev <- h2 
    aux <- h21 - (h1+eta +alpha*Q(y))
    if(any(is.na(aux))) STOP2 <- T
    else{
      h2 <- gam(aux~s(x2,df=DF2))$fitted
      ##for identifyability put mean in eta
      auxeta <- mean(h2*(1+exp(eta+h1+h2+alpha*Q(y))))
      h2 <- h2 - auxeta
      etaprev <- eta
      eta <- eta + auxeta
      
###now we can compute the p0's
      p0 <- ilogit(eta+h1+h2+alpha*Q(y))
    }
  }
 
  STOP <- any(is.na(p0)) 
  if(STOP | any(abs(h1)>maxh) | any(abs(h2)>maxh)) {
    Flag <- F ##maxh is too big
    cat("stopped at 2.")
  }
  else{
    Flag <-  (sum((h1-h1prev)^2+(h2-h2prev)^2+(eta-etaprev)^2)
              /sum(h1^2+h2^2+eta^2) > epsilon)
  }

}
if(Count>=Maxcount) {

  cat("Maxcount reached!, delta=",sum((h1-h1prev)^2+(h2-h2prev)^2+(eta-etaprev)^2)
              /sum(h1^2+h2^2+eta^2),"\n")
}

if(STOP | any(abs(h1)>maxh) | any(abs(h2)>maxh)) {
  STOP <- T;
  cat("\n\nNO ITERATION, stopped at 3\n\n")
  print(c(max(abs(h1)),max(abs(h2))))
} 






}