##starting values for gam
h1 <- trueh1[r==1]; h2 <- trueh2[r==1];eta <- trueeta
x1 <- Age[r==1]
x2 <- CD4[r==1]
y <- Y[r==1]

##we predict for these
px1 <- pAge[pr==1]
px2 <-pCD4[pr==1]
py <- pY[pr==1]

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)) STOP2 <- T 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{
    aux.formula <- paste("aux~s(x1,df=",DF1,")")
    h1 <- gam(aux.formula)
    h1predict <- predict.gam(h1,newdata=data.frame(x1=px1))
    h1 <- h1$fitted
    auxeta <- mean(h1*(1+exp(eta+h1+h2+alpha*Q(y))))
    h1 <- h1 - auxeta
    h1predict <- h1predict - auxeta
    eta <- eta + auxeta
  
    h2prev <- h2
    aux <- h21 - (h1+eta+alpha*Q(y))
    if(any(is.na(aux))) STOP2 <- T
    else{
      aux.formula <- paste("aux~s(x2,df=",DF2,")")
      h2 <- gam(aux.formula)
      h2predict <- predict.gam(h2,newdata=data.frame(x2=px2))
      h2 <- h2$fitted
      ##for identifyability put mean in eta
      auxeta <- mean(h2*(1+exp(eta+h1+h2+alpha*Q(y))))
      h2 <- h2 - auxeta
      h2predict <- h2predict- auxeta
      etaprev <- eta
      eta <- eta + auxeta
      
      ##now we can compute the p0's
      p0 <- ilogit(eta+h1+h2+alpha*Q(y))
      p0predict <- ilogit(eta+h1predict+h2predict+alpha*Q(py))
      STOP2 <- any(is.na(p0)) |  any(is.na(p0predict))
    }
  }
  if(STOP2 | any(abs(h1)>maxh) | any(abs(h2)>maxh)) Flag <- F ##maxh is too big
  else{
    Flag <-  (sum((h1-h1prev)^2+(h2-h2prev)^2)/sum(h1^2+h2^2) > epsilon)
  }
}