mu0prev <- sum(y/(1-p0))/sum(1/(1-p0))
##this is a newton raphson for mu
##alterated gam
phiprev1 <- rep(0,n); phiprev2 <- rep(0,n);
Flag2 <- T
Count2 <- 1
aux.formula1 <- paste("yy1~s(x1,df=",DF5,")")
aux.formula2 <- paste("yy2~s(x2,df=",DF6,")")
ww <- exp((h1+h2+alpha*Q(y))); ww <- ww/sum(ww)
while(Flag2 && Count2 < Maxcount){ 
    yy1 <- y-mu0prev-phiprev2[r==1]
    gam1 <- gam(aux.formula1,w=ww)
    phi1 <- predict.gam(gam1,newdata=data.frame(x1=Age))
    phi0 <- phi1 + phiprev2
    yy2 <- y-mu0prev-phi1[r==1]
    gam2 <- gam(aux.formula2,w=ww)
    phi2 <- predict.gam(gam2,newdata=data.frame(x2=CD4))
    phi0 <-  phi1 + phi2
    mu0 <- mu0prev + mean(weights.r*(Y-mu0prev) + (1-weights.r)*phi0)
    Flag2 <- ((mu0-mu0prev)^2 + mean((phi1-phiprev1)^2+(phi2-phiprev2)^2))/
      (mu0^2+mean(phi1^2+phi2^2)) > epsilon
    phiprev1 <- phi1;phiprev2 <- phi2;mu0prev <- mu0
    Count2 <- Count2 + 1
  }
  derivative <- -1