

\section{R codes for Chapter 3}


\small\begin{verbatim}

#program for illustrating the time-dependent covariates 
#scenario in notes 3/14/2014
# revised 10/29/2014 
# note that U= midpoint of L intervals in observation 
#censoring during the interval.


lucy.dat <- function(n=5,p,L,p0=2,eff=1,diff=0.5){

prob0=prob=x=C=T=d=U=lucy.y=matrix(0,n,p)

x0=matrix(rnorm(n*p0),n,p0)   # fixed covariates x0
b0=matrix(0.1,p0,1)
b=matrix(0.2,p0,1)


a0= 1.3*seq(1:p)/eff
a= seq(1:p)/eff

#------------------------------------------------------
#	baseline info
#-------------------------------------------------------

x[,1]=4+rnorm(n)

prob0[,1]=1/(exp(-a0[1]+0.9*x[,1]+x0%*%b0)+1)
C[,1]=rbinom(n,1,prob0[,1])

prob[,1]=1/(exp(-a[1]+0.8*x[,1]+x0%*%b)+1)
T[,1]=rbinom(n,1,prob[,1])
T[,1]=ifelse(C[,1]==0,T[,1],NA)

tmp1=L[1]/2
U[,1]=ifelse(C[,1]==1,tmp1, L[1])

xt=cbind(0,x[,1])
#lucy.y[,1]=-2.5+0.5*x0[,1]+0.5*x0[,2]-5*x[,1]+2*(1-d)
#----------------------------------
#	follow-up
#---------------------------------
for (i in 2:p){
	x[,i]=diff+rnorm(n)+x[,i-1]
	x[,i]=ifelse(T[,i-1]==0,x[,i],NA)


      prob0[,i]=1/(exp(-a0[i]+0.9*x[,i]+x0%*%b0)+1)
	C[,i]=rbinom(n,1,prob0[,i])

	prob[,i]=1/(exp(-a[i]+0.8*x[,i]+x0%*%b)+1)
	T[,i]=rbinom(n,1,prob[,i])
	T[,i]=ifelse(C[,i]==0,T[,i],NA)


	tmp2=(L[i-1]+L[i])/2
	U[,i]=ifelse(C[,i]==1,tmp2,L[i])
	xtmp=cbind(L[i-1],x[,i])
	xt=rbind(xt,xtmp)
}

Lmid=(L+c(0,L[-length(L)]))/2
L1=rep(L,each=n)
L2=rep(Lmid,each=n)
U2=apply(U,1,function(x) max(x,na.rm=TRUE))



T[,p]=ifelse(U2==L[p],1,NA)
d=rowSums(T,na.rm=TRUE)
colnames(xt)=c("t",paste("xt",seq(1:(dim(xt)[2]-1)),sep=""))
id=seq(1:n)

tmpx=x  # processed below to represent most recently xt

for (j in 1:(p)){
	if(j==1){tmpx[,j]=x[,1]}
	else{tmpx[,j]=apply(x[,1:j], 1,function(x) x[length(x[!is.na(x)])])}
}


#identify the last non-missing x
xx=apply(x,1,function(x) x[length(x[!is.na(x)])])

y=1.5+x0[,1]+0.5*x0[,2]-2*d+log(U2)+0.5*xx+rnorm(n)

dat.brent=data.frame(id=id,L=L1,U=U2,d=d,x0=x0,xt=xt,Lmid=L2,y=y,xx=xx)
dat.brent$yb=as.numeric(dat.brent$U == L & dat.brent$d == 1)

return(list(U=U2,d=d,T=T,C=C,x0=x0,xt=xt,tmpx=tmpx,y=y,x=x,
			xx=xx,id=id,df=dat.brent, prob0=prob0, prob=prob))
}
#######################################################
#-----------Simulated Truth----------------------------------
p.sim <-length(L)
n.sim=500000

sim.x0 <- matrix(NA,p.sim,n.sim)
sim.lambda <- matrix(NA,p.sim,n.sim)
sim.Cens <- sim.U<-sim.Delta<-matrix(NA, p.sim, n.sim)


sim.x1 <- rnorm(n.sim,0,1)
sim.x2 <- rnorm(n.sim,0,1)

Lmid=(L+c(0,L[-p.sim]))/2

for (k in 1){

	sim.x0[k,] <- rnorm(n.sim,(3.5+0.5*k),1)
	sim.lambda[k,] <- 1/(1+exp(-1.3*k/eff+0.9*sim.x0[k,]+0.1*sim.x1+0.1*sim.x2))	
	sim.Cens[k,]<-rbinom(n.sim,1,sim.lambda[k,])
	sim.U[k,]=ifelse(sim.Cens[k,]==1,Lmid[k], L[k])
	sim.Delta[k,]= ifelse(sim.Cens[k,]==0,1,0)

}

for (k in 2:p.sim){

	sim.x0[k,] <- diff+rnorm(n.sim)+sim.x0[k-1,]
	sim.lambda[k,] <- 1/(1+exp(-1.3*k/eff+0.9*sim.x0[k,]+0.1*sim.x1+0.1*sim.x2))	

	sim.Cens[k,]<-rbinom(n.sim,1,sim.lambda[k,])
	sim.Cens[k,]=ifelse(sim.Cens[k-1,]==0,sim.Cens[k,],NA)

	sim.Delta[k,]=ifelse(sim.Cens[k,]==1|is.na(sim.Cens[k,]),0,1)
	
	sim.U[k,]=ifelse(is.na(sim.Cens[k,]),sim.U[k-1,],Lmid[k])
	sim.U[k,]=ifelse(sim.Cens[k,]==1|is.na(sim.Cens[k,]),sim.U[k,],L[k])
}
x1rep=matrix(sim.x1,p.sim,n.sim,byrow=TRUE)
x2rep=matrix(sim.x2,p.sim,n.sim,byrow=TRUE)

sim.Y <- 1.5 +x1rep+0.5*x2rep-2*sim.Delta+0.5*sim.x0 +log(sim.U)
#sim.Y <- 1.5 +x1rep+0.5*x2rep+0.5*sim.x0 
trumu<-rowMeans(sim.Y)

t1=matrix(sim.Y,ncol=1)
t2=rep(L,n.sim)
t3=data.frame(Y=t1,L=t2)

ft.t=lm(Y~L,data=t3)
trumu2=predict(ft.t,newdata=list(L=L))
trubeta=ft.t$coeff

#--------------------------------------------------------------------
# 5_11_2014
#--------------------------------------------------------------------

Lucy.gcomp<- function(x0,time.x,x,u,d,y,L,Lx){

	n=length(u)
	#------------estimation of parameters
	ft3=lucy.discrete.haz.dep.cens(u=u,d=1-d,L=L,x0=x0,time.x=time.x,Lx=Lx)
	alpha=ft3$beta[1:4]
	bxt=ft3$beta[7]
	b0=ft3$beta[5:6]

	mu.x=colMeans(x0)
	s.x=apply(x0,2,sd)

	xt10=x[,1]
	mu.x0=mean(xt10)
	s.x0=sd(xt10)

	t20=!is.na(x[,2])
	xt10.20=xt10[t20]
	xt20.20=x[,2][t20]

	t30=!is.na(x[,3])
	xt20.30=x[,2][t30]
	xt30.30=x[,3][t30]

	t40=!is.na(x[,4])
	xt30.40=x[,3][t40]
	xt40.40=x[,4][t40]

	x.mk=c(xt10.20,xt20.30,xt30.40)
	y.mk=c(xt20.20,xt30.30,xt40.40)
	ft.mk=lm(y.mk~x.mk)
	mk.a=ft.mk$coeff[1]
	mk.b=ft.mk$coeff[2]

	m1=mu.x0
	m2=m1*mk.b+mk.a
	m3=m2*mk.b+mk.a
	m4=m3*mk.b+mk.a
	mt=c(m1,m2,m3,m4)

	#identify the last non-missing x
	last.x=apply(x,1,function(x) x[length(x[!is.na(x)])])

	ftlm=lm(y ~ x0+d+last.x+log(u))
	cf=ftlm$coef

#________________MC integration_________________

p.sim <-length(L)
n.sim=1000

sim.x0 <- matrix(NA,p.sim,n.sim)
sim.lambda <- matrix(NA,p.sim,n.sim)
sim.Cens <- sim.U<-sim.Delta<-matrix(NA, p.sim, n.sim)


sim.x1 <- rnorm(n.sim,mu.x[1],s.x[1])
sim.x2 <- rnorm(n.sim,mu.x[2],s.x[2])

Lmid=(L+c(0,L[-p.sim]))/2

for (k in 1){

	sim.x0[k,] <- rnorm(n.sim,mt[k],s.x0)
	sim.lambda[k,] <- 1/(1+exp(-alpha[k]-sim.x0[k,]*bxt-sim.x1*b0[1]-sim.x2*b0[2]))	
	sim.Cens[k,]<-rbinom(n.sim,1,sim.lambda[k,])
	sim.U[k,]=ifelse(sim.Cens[k,]==1,Lmid[k], L[k])
	sim.Delta[k,]= ifelse(sim.Cens[k,]==0,1,0)
}

for (k in 2:p.sim){

	sim.x0[k,] <- rnorm(n.sim,mk.a,s.x0)+sim.x0[k-1,]*mk.b
	sim.lambda[k,] <- 1/(1+exp(-alpha[k]-sim.x0[k,]*bxt-sim.x1*b0[1]
										-sim.x2*b0[2]))	

	sim.Cens[k,]<-rbinom(n.sim,1,sim.lambda[k,])
	sim.Cens[k,]=ifelse(sim.Cens[k-1,]==0,sim.Cens[k,],NA)

	sim.Delta[k,]=ifelse(sim.Cens[k,]==1|is.na(sim.Cens[k,]),0,1)
	
	sim.U[k,]=ifelse(is.na(sim.Cens[k,]),sim.U[k-1,],Lmid[k])
	sim.U[k,]=ifelse(sim.Cens[k,]==1|is.na(sim.Cens[k,]),sim.U[k,],L[k])
}
x1rep=matrix(sim.x1,p.sim,n.sim,byrow=TRUE)
x2rep=matrix(sim.x2,p.sim,n.sim,byrow=TRUE)

sim.Y <- cf[1] +cf[2]*x1rep+cf[3]*x2rep+ cf[4]*sim.Delta+cf[5]*sim.x0
					+cf[6]*log(sim.U) 


mu_hat <- rowMeans(sim.Y)

t1=matrix(sim.Y,ncol=1)
t2=rep(L,n.sim)
t3=data.frame(Y=t1,L=t2)

ft.t=lm(Y~L,data=t3)
beta_hat=ft.t$coeff
mu_hat_model=predict(ft.t,newdata=list(L=L))

return(list(mu=mu_hat,beta=ft3$beta, mu2=mu_hat_model,b=beta_hat))

}




#--------------------------------------------------------------------
# 8_15_2014
#--------------------------------------------------------------------

Lucy.gcomp.pred<- function(x0,time.x,x,u,d,y,L,Lx){

	n=length(u)
	#------------estimation of parameters
	ft3=lucy.discrete.haz.dep.cens(u=u,d=1-d,L=L,x0=x0,time.x=time.x,Lx=Lx)
	alpha=ft3$beta[1:4]
	bxt=ft3$beta[7]
	b0=ft3$beta[5:6]

	mu.x=colMeans(x0)
	s.x=apply(x0,2,sd)

	xt10=x[,1]
	mu.x0=mean(xt10)
	s.x0=sd(xt10)

	t20=!is.na(x[,2])
	xt10.20=xt10[t20]
	xt20.20=x[,2][t20]

	t30=!is.na(x[,3])
	xt20.30=x[,2][t30]
	xt30.30=x[,3][t30]

	t40=!is.na(x[,4])
	xt30.40=x[,3][t40]
	xt40.40=x[,4][t40]

	x.mk=c(xt10.20,xt20.30,xt30.40)
	y.mk=c(xt20.20,xt30.30,xt40.40)
	ft.mk=lm(y.mk~x.mk)
	mk.a=ft.mk$coeff[1]
	mk.b=ft.mk$coeff[2]

	m1=mu.x0
	m2=m1*mk.b+mk.a
	m3=m2*mk.b+mk.a
	m4=m3*mk.b+mk.a
	mt=c(m1,m2,m3,m4)

	#identify the last non-missing x
	last.x=apply(x,1,function(x) x[length(x[!is.na(x)])])

	ftlm=lm(y ~ x0+d+last.x+log(u))
	cf=ftlm$coef

#________________prediction_________________

predY=pred.lambda=pred.prob=matrix(NA,length(L),n)
tmpd=matrix(d,n,length(L))
tmpu=matrix(u,n,length(L))
tmpL=matrix(L,n,length(L),byrow=TRUE)

tmpp=tmpu<=tmpL
tmpd=ifelse(tmpp,tmpd,1)
tmpu=ifelse(tmpp,tmpu,tmpL)


myx2=x[,1:2]
tmpx2=apply(myx2,1,function(x) x[length(x[!is.na(x)])])
myx3=x[,1:3]
tmpx3=apply(myx3,1,function(x) x[length(x[!is.na(x)])])
myx4=x[,1:4]
tmpx4=apply(myx4,1,function(x) x[length(x[!is.na(x)])])

pred.lambda[1,] <- 1/(1+exp(-alpha[1]-x[,1]*bxt-x0[,1]*b0[1]-x0[,2]*b0[2]))
pred.lambda[2,] <- 1/(1+exp(-alpha[2]-tmpx2*bxt-x0[,1]*b0[1]-x0[,2]*b0[2]))
pred.lambda[3,] <- 1/(1+exp(-alpha[3]-tmpx3*bxt-x0[,1]*b0[1]-x0[,2]*b0[2]))
pred.lambda[4,] <- 1/(1+exp(-alpha[4]-tmpx4*bxt-x0[,1]*b0[1]-x0[,2]*b0[2]))


p1.1=(1-pred.lambda[1,])
p1.0=pred.lambda[1,]

part1=cf[1] +cf[2]*x0[,1]+cf[3]*x0[,2]+ cf[4]+cf[5]*x[,1] +cf[6]*log(L[1])
part2=cf[1] +cf[2]*x0[,1]+cf[3]*x0[,2]+cf[5]*x[,1] +cf[6]*log(Lmid[1])

predY[1,]=part1*p1.1+part2*p1.0


p2.1=(1-pred.lambda[1,])*(1-pred.lambda[2,])
p2.0a=pred.lambda[1,]
p2.0b=(1-pred.lambda[1,])*pred.lambda[2,]


m2.1=cf[1] +cf[2]*x0[,1]+cf[3]*x0[,2]+ cf[4]+cf[5]*tmpx2 +cf[6]*log(L[2])
m2.0a=cf[1] +cf[2]*x0[,1]+cf[3]*x0[,2]+ cf[5]*tmpx2 +cf[6]*log(Lmid[1])
m2.0b=cf[1] +cf[2]*x0[,1]+cf[3]*x0[,2]+ cf[5]*tmpx2 +cf[6]*log(Lmid[2])
	
predY[2,]=p2.1*m2.1+p2.0a*m2.0a+p2.0b*m2.0b


p3.1=(1-pred.lambda[1,])*(1-pred.lambda[2,])*(1-pred.lambda[3,])
p3.0a=pred.lambda[1,]
p3.0b=(1-pred.lambda[1,])*pred.lambda[2,]
p3.0c=(1-pred.lambda[1,])*(1-pred.lambda[2,])*pred.lambda[3,]

m3.1=cf[1] +cf[2]*x0[,1]+cf[3]*x0[,2]+ cf[4]+cf[5]*tmpx3 +cf[6]*log(L[3])
m3.0a=cf[1] +cf[2]*x0[,1]+cf[3]*x0[,2]+ cf[5]*tmpx3 +cf[6]*log(Lmid[1])
m3.0b=cf[1] +cf[2]*x0[,1]+cf[3]*x0[,2]+ cf[5]*tmpx3 +cf[6]*log(Lmid[2])
m3.0c=cf[1] +cf[2]*x0[,1]+cf[3]*x0[,2]+ cf[5]*tmpx3 +cf[6]*log(Lmid[3])	

predY[3,]=p3.1*m3.1+p3.0a*m3.0a+p3.0b*m3.0b+p3.0c*m3.0c


p4.1=(1-pred.lambda[1,])*(1-pred.lambda[2,])*(1-pred.lambda[3,])*(
				1-pred.lambda[4,])
p4.0a=pred.lambda[1,]
p4.0b=(1-pred.lambda[1,])*pred.lambda[2,]
p4.0c=(1-pred.lambda[1,])*(1-pred.lambda[2,])*pred.lambda[3,]
p4.0d=(1-pred.lambda[1,])*(1-pred.lambda[2,])*(1-pred.lambda[3,]
					)*pred.lambda[4,]

m4.1=cf[1] +cf[2]*x0[,1]+cf[3]*x0[,2]+ cf[4]+cf[5]*tmpx4 +cf[6]*log(L[4])
m4.0a=cf[1] +cf[2]*x0[,1]+cf[3]*x0[,2]+ cf[5]*tmpx4 +cf[6]*log(Lmid[1])
m4.0b=cf[1] +cf[2]*x0[,1]+cf[3]*x0[,2]+ cf[5]*tmpx4 +cf[6]*log(Lmid[2])
m4.0c=cf[1] +cf[2]*x0[,1]+cf[3]*x0[,2]+ cf[5]*tmpx4 +cf[6]*log(Lmid[3])	
m4.0d=cf[1] +cf[2]*x0[,1]+cf[3]*x0[,2]+ cf[5]*tmpx4 +cf[6]*log(Lmid[4])	

predY[4,]=p4.1*m4.1+p4.0a*m4.0a+p4.0b*m4.0b+p4.0c*m4.0c+p4.0d*m4.0d

mu_hat <- rowMeans(predY)


t1=matrix(predY,ncol=1)
t2=rep(L,n)
t3=data.frame(Y=t1,L=t2)

ft.t=lm(Y~L,data=t3)
beta_hat=ft.t$coeff
mu_hat_model=predict(ft.t,newdata=list(L=L))

return(list(mu=mu_hat,beta=ft3$beta, mu2=mu_hat_model,b=beta_hat))

}




####################################
# wrote 4_18_2014
# add time dependent time.x

#first column of time.x is Lx, it is the time when time-dep
# covariates are measured
# .cens is when considering censoring rather than event. 
#Number of time points in the models are different

lucy.discrete.haz.dep.cens <- function (x0=NULL,time.x=NULL,u,d,L,Lx){
	
	K=length(L)
	a0=rep(NA,K)
	n=length(u)

	if(is.null(x0)){
	p0=0
	}else{
	x0=as.matrix(x0)
	p0= dim(x0)[2]
	}	
	
	if(is.null(time.x)){
	pt=0
	}else{
	time.x=as.matrix(time.x)
	pt0= (dim(time.x)[2]-1)
	pt= pt0
	}
	

	#----get initial value for intercept-----------
	LL=c(0,L[1:(K-1)])
	midpt=(L+LL)/2

	ind1= matrix(u,n,K)==matrix(midpt,n,K,byrow=TRUE) & matrix(d,n,K)==1
	ind2=	matrix(u,n,K)>=matrix(LL,n,K,byrow=TRUE)

	p1= colSums(ind1)/colSums(ind2)
	p1=ifelse(p1==0,0.000001,p1)
	#print(p1)
	a0= log(p1/(1-p1))
	#print(a0)	
	xnew0=as.matrix(cbind(1,x0))

	p=p0+pt
	theta2=matrix(c(a0,rep(0,p)),p+K,1)
	

	theta=0
	iter <- 0
	itermax=150
	tolerance=0.001  

	#---Newton Ralphson-------------------------------
	while (sum(abs(theta - theta2)) >= tolerance) {	
		#print(abs(theta - theta2)[4]) 
		iter=iter+1
		#print(paste("iter =",iter))

		
		Ha=Sa=matrix(NA,1,K)
		Sb=Hab=matrix(NA,K,p)
		Hb=0
	
		theta=theta2
		xdf=NULL

	if(!is.null(x0)& !is.null(time.x)){
		
		for (j in 1:K){
		#print(paste("J =",j))
		
		bnew=theta[c(j,(K+1):(K+p))]
		
		Ltmp=max(Lx[which(Lx<L[j])])

		xt=time.x[which(time.x[,1]==Ltmp),-1]
		pos=which(!is.na(xt))
		
		xnew=as.matrix(cbind(1,x0,xt))[pos,]
		x=as.matrix(cbind(x0,xt))[pos,]

		xb=xnew%*%bnew

		lp=exp(-xb)
		Flp=lp/((1+lp)^2)
		Ha[j]=sum(ind2[pos,j]*Flp)
		Sa[j]=sum(ind1[pos,j])-sum(ind2[pos,j]/(lp+1))

		Hb=Hb+t(x)%*%(matrix(ind2[pos,j]*Flp,length(pos),p)*x)
		Sb[j,]=t(ind1[pos,j])%*%x-t(as.matrix(ind2[pos,j]/(lp+1)))%*%x
		Hab[j,]=t(ind2[pos,j]*Flp)%*%x
		xtmp=cbind(L=L[j],xnew)
		xdf=rbind(xdf,xtmp)
		}
	
		S=rbind(t(Sa),as.matrix(colSums(Sb)))
		H1=cbind(diag(as.vector(Ha)),Hab)
		H2=cbind(t(Hab),Hb)
		H=rbind(H1,H2)
	}

	theta2=theta+solve(H)%*%S
	beta=t(theta2)
	colnames(beta)=c( paste("a",seq(1:K),sep=""),
				paste("x0",seq(1:p0),sep=""),
				paste("xt",seq(1:pt0),sep="") )
	V=solve(H)
	colnames(V)=rownames(V)=colnames(beta)

	#print(beta)

	}
return(list(beta=beta,V=V,xdf=xdf))
}


#____________________with respect to event

lucy.discrete.haz.dep <- function (x0=NULL,time.x=NULL,u,d,L,Lx){
	
	K=length(L)-1
	a0=rep(NA,K)
	n=length(u)

	if(is.null(x0)){
	p0=0
	}else{
	x0=as.matrix(x0)
	p0= dim(x0)[2]
	}	
	
	if(is.null(time.x)){
	pt=0
	}else{
	time.x=as.matrix(time.x)
	pt0= (dim(time.x)[2]-1)
	pt= pt0
	}
	

	#----get initial value for intercept-----------
	ind1= matrix(u,n,K)== matrix(L[1:K],n,K,byrow=TRUE) & matrix(d,n,K)==1
	ind2=	matrix(u,n,K)>=matrix(L[1:K],n,K,byrow=TRUE)

	p1= colSums(ind1)/colSums(ind2)
	p1=ifelse(p1==0,0.000001,p1)
	a0= log(p1/(1-p1))	
	xnew0=as.matrix(cbind(1,x0))

	p=p0+pt
	theta2=matrix(c(a0,rep(0,p)),p+K,1)
	

	theta=0
	iter <- 0
	itermax=150
	tolerance=0.00000000001  

	#---Newton Ralphson-------------------------------
	while (sum(abs(theta - theta2)) >= tolerance) {	
		iter=iter+1
#		print(paste("iter =",iter))

		
		Ha=Sa=matrix(NA,1,K)
		Sb=Hab=matrix(NA,K,p)
		Hb=0
	
		theta=theta2
		xdf=NULL

	if(!is.null(x0)& !is.null(time.x)){
		
		for (j in 1:K){
		#print(paste("J =",j))
		
		bnew=theta[c(j,(K+1):(K+p))]
		
		Ltmp=max(Lx[which(Lx<L[j])])

		xt=time.x[which(time.x[,1]==Ltmp),-1]
		pos=which(!is.na(xt))
		
		xnew=as.matrix(cbind(1,x0,xt))[pos,]
		x=as.matrix(cbind(x0,xt))[pos,]

		xb=xnew%*%bnew

		lp=exp(-xb)
		Flp=lp/((1+lp)^2)
		Ha[j]=sum(ind2[pos,j]*Flp)
		Sa[j]=sum(ind1[pos,j])-sum(ind2[pos,j]/(lp+1))

		Hb=Hb+t(x)%*%(matrix(ind2[pos,j]*Flp,length(pos),p)*x)
		Sb[j,]=t(ind1[pos,j])%*%x-t(as.matrix(ind2[pos,j]/(lp+1)))%*%x
		Hab[j,]=t(ind2[pos,j]*Flp)%*%x
		xtmp=cbind(L=L[j],xnew)
		xdf=rbind(xdf,xtmp)
		}
	
		S=rbind(t(Sa),as.matrix(colSums(Sb)))
		if (length(Ha)==1){
		H1=cbind(Ha,Hab)
		}
		if (length(Ha)>1){
		H1=cbind(diag(as.vector(Ha)),Hab)
		}
		H2=cbind(t(Hab),Hb)
		H=rbind(H1,H2)
	}

	theta2=theta+solve(H)%*%S
	beta=t(theta2)
	colnames(beta)=c( paste("a",seq(1:K),sep=""),
				paste("x0",seq(1:p0),sep=""),
				paste("xt",seq(1:pt0),sep="") )
	V=solve(H)
	colnames(V)=rownames(V)=colnames(beta)

#	print(beta)

	}
return(list(beta=beta,V=V,xdf=xdf))
}

\end{verbatim}




\section{R codes for Chapter 4}
\begin{verbatim}
#################################################################
# 2 stage  y is indep of U,  Dec,2014
######################################################################
lucy.dat2 <- function(n=5,p=2,L,p0=2,eff=1,diff=0){

prob0=prob=x=C=T=d=U=lucy.y=matrix(0,n,p)

x0=matrix(rnorm(n*p0),n,p0)   # fixed covariates x0
b0=matrix(0.1,p0,1)
b=matrix(0.2,p0,1)


a0= 1.3*seq(1:p)/eff
a= seq(1:p)/eff

#------------------------------------------------------
#	baseline info
#-------------------------------------------------------

x[,1]=4+rnorm(n)

prob0[,1]=1/(exp(-a0[1]+0.9*x[,1]+x0%*%b0)+1)
C[,1]=rbinom(n,1,prob0[,1])
y1= 1.5+x0[,1]+0.5*x0[,2]+0.5*x[,1]+rnorm(n)
y1=ifelse(C[,1]==1,y1,NA)

prob[,1]=1/(exp(-a[1]+0.8*x[,1]+x0%*%b)+1)
T[,1]=rbinom(n,1,prob[,1])
T[,1]=ifelse(C[,1]==0,T[,1],NA)
yt=1.5+x0[,1]+0.5*x0[,2]+0.5*x[,1]+0.5+rnorm(n)
y=ifelse(T[,1]==1 & !is.na(T[,1]),yt,y1)

tmp1=L[1]/2
U[,1]=ifelse(C[,1]==1,tmp1, L[1])

xt=cbind(0,x[,1])

#----------------------------------
#	follow-up
#---------------------------------
for (i in 2:p){
	x[,i]=diff+rnorm(n)+x[,i-1]
	x[,i]=ifelse(T[,i-1]==0,x[,i],NA)


      prob0[,i]=1/(exp(-a0[i]+0.9*x[,i]+x0%*%b0)+1)
	C[,i]=rbinom(n,1,prob0[,i])
	yt= 1.5+x0[,1]+0.5*x0[,2]+0.5*x[,i]+0.5+rnorm(n)
	y=ifelse(C[,i]==1& !is.na(C[,i]),yt,y)

	prob[,i]=1/(exp(-a[i]+0.8*x[,i]+x0%*%b)+1)
	T[,i]=rbinom(n,1,prob[,i])
	T[,i]=ifelse(C[,i]==0,T[,i],NA)
	yt= 1.5+x0[,1]+0.5*x0[,2]+0.5*x[,i]+0.5+0.3+rnorm(n)
	y=ifelse(T[,i]==1& !is.na(T[,i]),yt,y)

	tmp2=(L[i-1]+L[i])/2
	U[,i]=ifelse(C[,i]==1,tmp2,L[i])
	xtmp=cbind(L[i-1],x[,i])
	xt=rbind(xt,xtmp)
}

Lmid=(L+c(0,L[-length(L)]))/2
L1=rep(L,each=n)
L2=rep(Lmid,each=n)
U2=apply(U,1,function(x) max(x,na.rm=TRUE))



T[,p]=ifelse(U2==L[p],1,NA)
yt= 1.5+x0[,1]+0.5*x0[,2]+0.5*x[,p]+0.5+0.3+rnorm(n)
y=ifelse(T[,p]==1& !is.na(T[,p]),yt,y)

d=rowSums(T,na.rm=TRUE)
colnames(xt)=c("t",paste("xt",seq(1:(dim(xt)[2]-1)),sep=""))
id=seq(1:n)

tmpx=x  # processed below to represent most recently xt

for (j in 1:(p)){
	if(j==1){tmpx[,j]=x[,1]}
	else{tmpx[,j]=apply(x[,1:j], 1,function(x) x[length(x[!is.na(x)])])}
}


#identify the last non-missing x
xx=apply(x,1,function(x) x[length(x[!is.na(x)])])

dat.brent=data.frame(id=id,L=L1,U=U2,d=d,x0=x0,xt=xt,Lmid=L2,y=y,xx=xx)
dat.brent$yb=as.numeric(dat.brent$U == L & dat.brent$d == 1)

return(list(U=U2,d=d,T=T,C=C,x0=x0,xt=xt,tmpx=tmpx,y=y,x=x,xx=xx,
id=id,df=dat.brent, prob0=prob0, prob=prob))
}

############################################
#simulated truth for lucy.dat2
##########################################

my.simtru<-function(L,eff=2,diff=0){

p.sim <-length(L)
n.sim=500000

sim.x0 <- matrix(NA,p.sim,n.sim)
sim.lambda <- matrix(NA,p.sim,n.sim)
sim.Y<-sim.Cens <- sim.U<-sim.Delta<-matrix(NA, p.sim, n.sim)

sim.x1 <- rnorm(n.sim,0,1)
sim.x2 <- rnorm(n.sim,0,1)

Lmid=(L+c(0,L[-p.sim]))/2

for (k in 1){

	sim.x0[k,] <- rnorm(n.sim,(3.5+0.5*k),1)
	sim.lambda[k,] <- 1/(1+exp(-1.3*k/eff+0.9*sim.x0[k,]+0.1*sim.x1+0.1*sim.x2))	
	sim.Cens[k,]<-rbinom(n.sim,1,sim.lambda[k,])
	sim.U[k,]=ifelse(sim.Cens[k,]==1,Lmid[k], L[k])
	sim.Delta[k,]= ifelse(sim.Cens[k,]==0,1,0)
	sim.y0= 1.5+sim.x1+0.5*sim.x2+0.5*sim.x0[k,]
	sim.yt=1.5+sim.x1+0.5*sim.x2+0.5*sim.x0[k,]+0.5
	sim.Y[k,]=ifelse(sim.Cens[k,]==1,sim.y0,sim.yt)
}

for (k in 2:p.sim){

	sim.x0[k,] <- diff+rnorm(n.sim)+sim.x0[k-1,]
	sim.lambda[k,] <- 1/(1+exp(-1.3*k/eff+0.9*sim.x0[k,]+0.1*sim.x1+0.1*sim.x2))	

	sim.Cens[k,]<-rbinom(n.sim,1,sim.lambda[k,])
	sim.Cens[k,]=ifelse(sim.Cens[k-1,]==0,sim.Cens[k,],NA)

	sim.Delta[k,]=ifelse(sim.Cens[k,]==1|is.na(sim.Cens[k,]),0,1)
	
	sim.U[k,]=ifelse(is.na(sim.Cens[k,]),sim.U[k-1,],Lmid[k])
	sim.U[k,]=ifelse(sim.Cens[k,]==1|is.na(sim.Cens[k,]),sim.U[k,],L[k])

	sim.y02= 1.5+sim.x1+0.5*sim.x2+0.5*sim.x0[k,]+0.5
	sim.yt=1.5+sim.x1+0.5*sim.x2+0.5*sim.x0[k,]+0.5+0.3

	sim.Y[k,]=ifelse(sim.Cens[1,]==1,sim.y0,0)
	sim.Y[k,]=ifelse(sim.Cens[k,]==1& !is.na(sim.Cens[k,]),sim.y02,
						sim.Y[k,])
	sim.Y[k,]=ifelse(sim.Cens[k,]==0 & !is.na(sim.Cens[k,]),sim.yt,
						sim.Y[k,])


}

trumu<-rowMeans(sim.Y)
return(trumu)
}

#4_9_2015
# add the c parameter of control variance when ps is missspecified
#5_2_2015


lucy.DR.m2 <- function(y,U,L,x0,d,time.x,x,tmpx,xx){

	ind1=ifelse(U<L[1],0,1) #I(C>t1)
	ind2=ifelse(U==L[2],1,0) # I(C>t2)
	ind3=ifelse(U==L[1],1,0) # I(U=t1)

	p0=dim(x0)[2]
	n=length(y)
	#------------estimation of parameters-------------------------

	ft3=lucy.discrete.haz.dep.cens(u=U,d=1-d,L=L,x0=x0,time.x=time.x,Lx=L-10)
	alpha=ft3$beta[1:p]
	bxt=ft3$beta[p+p0+1]
	b0=ft3$beta[(p+1):(p+p0)]

	ft4=lucy.discrete.haz.dep(u=U,d=d,L=L,x0=x0,time.x=time.x,Lx=L-10)
	a=ft4$beta[1:(p-1)]
	bxt4=ft4$beta[p+p0]
	b=ft4$beta[(p):(p+p0-1)]
	ftm=lm(x[,2]~x[,1])

	# predict lambda's and xt1 
	pred.xt1=ftm$coeff[1]+ftm$coeff[2]*x[,1]

	pred.lambda0.0=1/(1+exp(-alpha[1]-x[,1]*bxt-x0%*%b0))
	pred.lambda0.1=1/(1+exp(-alpha[2]-pred.xt1*bxt-x0%*%b0))
	
	#pred.lambda0.0=1/(1+exp(-alpha[1]-x0%*%b0))
	#pred.lambda0.1=1/(1+exp(-alpha[2]-x0%*%b0))

	pred.lambda.1<- 1/(1+exp(-a[1]-x[,1]*bxt4-x0[,1]*b[1]-x0[,2]*b[2]))
	

	# estimate w wt
	ft1=calc.w(u=U,d=d,L=L,x0=x0,time.x=time.x,Lx=L-10,M=1,tmpx=tmpx)
	ft2=calc.w(u=U,d=d,L=L,x0=x0,time.x=time.x,Lx=L-10,M=2,tmpx=tmpx)


############################## 
# solving alpha_g for m=1 
################################
# a%*%x=b
wt1.1=1
p1d=cbind(1,ind1,x0,x[,1])
MxA=(t(p1d)*matrix(wt1.1,ncol=n,nrow=dim(x0)[2]+3,byrow=TRUE))%*%p1d/n
MxB=t(p1d)%*%as.matrix(y*wt1.1)/n
sol1=solve(MxA,MxB)


# variance minimization
wt.i1=ind1*(1-1/pred.lambda.1)


p1d=cbind(1,ind1,x0,x[,1])
MxA.i=(t(p1d)*matrix(wt.i1,ncol=n,nrow=dim(x0)[2]+3,byrow=TRUE))%*%p1d/n
MxB.i=t(p1d)%*%as.matrix(y*wt.i1)/n
sol1.t=solve(MxA.i[2:5,2:5],MxB.i[2:5])
sol1.i=c(0,sol1.t)

#cT adjustment
#derivative of g1.1
dg11=cbind(ind1,x0,x[,1])
	
#d(lambda.1)/d(gamma)=lambda.1^2*exp(-gammaX)  .X
dlam.nox=pred.lambda.1^2*exp(-a[1]-x[,1]*bxt4-x0[,1]*b[1]-x0[,2]*b[2])
dx=cbind(1,x[,1],x0)

theta2=c(sol1.t,rep(0.1, dim(dx)[2]))
	iter <- 0
	#itermax=150
	tolerance=0.000001  
theta=0
	#---Newton Ralphson-------------------------------
	while (sum(abs(theta - theta2)) >= tolerance) {	
		#print(abs(theta - theta2)[4]) 
		iter=iter+1
		#print(paste("iter =",iter))
	
		theta=theta2
		c=theta[(length(sol1.t)+1):length(theta)]
		ag=theta[1:length(sol1.t)]
g11=dg11%*%ag
l1=(dlam.nox/(1-pred.lambda.1)*(dx%*%c))

Sa=t(ind1*(1/pred.lambda.1-1)*(y-g11-l1))%*% dg11
Sb=t(ind1/pred.lambda.1*(y-g11-l1)*dlam.nox)%*%dx
S=cbind(Sa,Sb)

#-I(theta)
Iaa=t(matrix(ind1*(1/pred.lambda.1-1),ncol=dim(dg11)[2],nrow=n)*dg11)%*%dg11
Ica=t(matrix(ind1/pred.lambda.1*dlam.nox,ncol=dim(dx)[2],nrow=n)*dx)%*%dg11
Icc=t(matrix(ind1*dlam.nox^2/((1-pred.lambda.1)*pred.lambda.1),
			ncol=dim(dx)[2],nrow=n)*(dx))%*%dx
I=rbind(cbind(Iaa,t(Ica)),cbind(Ica, Icc))

theta2=theta+solve(I)%*%t(S)

	#print(theta2)

	}




############################## 
# solving alpha_g for m=2 
################################
wt1=1
wt2=1#1-ind3

q1d=cbind(1,ind1,ind2,x0,xx) # 1st derivative for g2^(2)
prc2=(1-pred.lambda0.1)

xlast=pred.xt1*ind1+ x[,1]*(1-ind1)
q2d=cbind(1,ind1,ind1*prc2,x0,xlast)  # 1st derivative for g1^(2)
# 1st derivative for sum p(a|s)*g2^(2)-g1^(2)
q2p=cbind(alpha0=0,alpha1=0,alpha2=ind1*(ind2-prc2),alphab1=0,alphab2=0,
			alphat=ind1*(xx-pred.xt1))  

Mxb=t(q1d)%*%as.matrix(y*wt1)/n
Mxa=((t(q1d)*matrix(wt1,ncol=n,nrow=dim(x0)[2]+4,byrow=TRUE))
%*%q1d-(t(q2d)*matrix(wt2,ncol=n,nrow=dim(x0)[2]+4,byrow=TRUE))%*%q2p)/n

sol2=solve(Mxa,Mxb)



# variance minimization
wt.i2=ind1*(pred.lambda.1/(1-pred.lambda.1))


Mxb.i=t(q1d)%*%as.matrix(y*wt.i2)/n
Mxa.i=((t(q1d)*matrix(wt.i2,ncol=n,nrow=dim(x0)[2]+4,byrow=TRUE))
%*%q1d-(t(q2d)*matrix(0,ncol=n,nrow=dim(x0)[2]+4,byrow=TRUE))%*%q2p)/n

sol2t=solve(Mxa.i[2:6,2:6],Mxb.i[2:6])
sol2.i=c(0,sol2t)

#cT adjustment
#derivative of g1.2
#dg120=cbind(ind1,x0,xlast,ind1*prc2)
dg12=cbind(ind1,ind1*prc2,x0,xlast)
#	dg12=cbind(ind1,ind2,x0,xx)
#d(lambda.1)/d(gamma)=lambda.1^2*exp(-gammaX)  .X
dlam.nox2=(pred.lambda.1)^2*exp(-a[1]-x[,1]*bxt4-x0[,1]*b[1]-x0[,2]*b[2])
dx=cbind(1,x[,1],x0)

theta4=c(sol2t,rep(-0.1, dim(dx)[2]))
	iter <- 0
	#itermax=150
	tolerance=0.000001  
theta3=0
	#---Newton Ralphson-------------------------------
	while (sum(abs(theta3 - theta4)) >= tolerance) {	
		#print(abs(theta3 - theta4)[4]) 
		iter=iter+1
		#print(paste("iter =",iter))
	
		theta3=theta4
		c2=theta3[(length(sol2t)+1):length(theta3)]
		ag2=theta3[1:length(sol2t)]
g12=dg12%*%ag2
l2=(dlam.nox2/(pred.lambda.1)*(dx%*%c2))

Sa2=t(ind1*(pred.lambda.1/((1-pred.lambda.1)))*(y-g12-l2))%*% dg12
Sb2=t(ind1/((1-pred.lambda.1))*(y-g12-l2)*dlam.nox2)%*%dx
S2=cbind(Sa2,Sb2)

#-I(theta)
Iaa2=t(matrix(ind1*(pred.lambda.1/((1-pred.lambda.1))),
			ncol=dim(dg12)[2],nrow=n)*dg12)%*%dg12
Ica2=t(matrix(ind1/((1-pred.lambda.1))*dlam.nox2,ncol=dim(dx)[2],
				nrow=n)*dx)%*%dg12
Icc2=t(matrix(ind1*dlam.nox2^2/((1-pred.lambda.1)*pred.lambda.1),
			ncol=dim(dx)[2],nrow=n)*(dx))%*%dx
I2=rbind(cbind(Iaa2,t(Ica2)),cbind(Ica2, Icc2))

theta4=theta3+solve(I2)%*%t(S2)

	#print(theta4)

	}





####### fit in g funtions########################


g1.1=p1d%*%as.matrix(sol2[-3])
g2.2=q1d%*%as.matrix(sol2)
g1.2=q2d%*%as.matrix(sol2)
g1.1i=p1d%*%as.matrix(sol1.i)
g1.2i=q2d%*%as.matrix(sol2.i)

g1.1c=dg11%*%as.matrix(theta2[1:length(sol1.t)])
g1.2c=cbind(ind1,ind1*prc2,x0,xlast)%*%as.matrix(theta4[1:length(sol2t)])



mu.dr=mu.jt=mu.or=mu.or.i=mu.dr.i=mu.dr.c=rep(NA,length(L))
#------------DR m=1-------------------------
mu.dr[1]=sum(ft1$W*y-ft1$W*g1.1+g1.1)/n
mu.jt[1]=sum(ft1$W*y)/sum(ft1$W)
mu.or[1]=mean(g1.1)
mu.or.i[1]=mean(g1.1i)
mu.dr.i[1]=sum(ft1$W*y-ft1$W*g1.1i+g1.1i)/n
mu.dr.c[1]=sum(ft1$W*y-ft1$W*g1.1c+g1.1c)/n


#------------DR m=2-------------------------
t.num.d=sum((ft2$W*y-(ft2$W-1)*g1.2))
t.num.di=sum((ft2$W*y-(ft2$W-1)*g1.2i))	#(1-ind3)*
t.num.dc=sum((ft2$W*y-(ft2$W-1)*g1.2c))

t.den=n#sum(1-ind3)
mu.jt[2]=sum(ft2$W*y)/sum(ft2$W)
mu.or[2]=mean(g1.2)
mu.dr[2]=t.num.d/t.den
mu.dr.i[2]=t.num.di/t.den
mu.or.i[2]=mean(g1.2i)
mu.dr.c[2]=t.num.dc/t.den

return(list(mu.dr=mu.dr,mu.jt=mu.jt,mu.or=mu.or,
		mu.or.i=mu.or.i, mu.dr.i=mu.dr.i,mu.dr.c=mu.dr.c))
}

#################### omit the time dependent x#############



lucy.DR.m2.or <- function(y=dat$y,U=dat$U,L=L,x0=dat$x0,d=dat$d,
							time.x=dat$xt,x=dat$x,tmpx=dat$tmpx,xx=dat$xx){

	ind1=ifelse(U<L[1],0,1) #I(C>t1)
	ind2=ifelse(U==L[2],1,0) # I(C>t2)
	ind3=ifelse(U==L[1],1,0) # I(U=t1)

	p0=dim(x0)[2]
	n=length(y)
	#------------estimation of parameters-------------------------

	ft3=lucy.discrete.haz.dep.cens(u=U,d=1-d,L=L,x0=x0,time.x=time.x,Lx=L-10)
	alpha=ft3$beta[1:p]
	bxt=ft3$beta[p+p0+1]
	b0=ft3$beta[(p+1):(p+p0)]

	ft4=lucy.discrete.haz.dep(u=U,d=d,L=L,x0=x0,time.x=time.x,Lx=L-10)
	a=ft4$beta[1:(p-1)]
	bxt4=ft4$beta[p+p0]
	b=ft4$beta[(p):(p+p0-1)]
	ftm=lm(x[,2]~x[,1])

	# predict lambda's and xt1 
	pred.xt1=ftm$coeff[1]+ftm$coeff[2]*x[,1]

	pred.lambda0.0=1/(1+exp(-alpha[1]-x[,1]*bxt-x0%*%b0))
	pred.lambda0.1=1/(1+exp(-alpha[2]-pred.xt1*bxt-x0%*%b0))
	pred.lambda.1<- 1/(1+exp(-a[1]-x[,1]*bxt4-x0[,1]*b[1]-x0[,2]*b[2]))
	

	# estimate w wt
	ft1=calc.w(u=U,d=d,L=L,x0=x0,time.x=time.x,Lx=L-10,M=1,tmpx=tmpx)
	ft2=calc.w(u=U,d=d,L=L,x0=x0,time.x=time.x,Lx=L-10,M=2,tmpx=tmpx)


############################## 
# solving alpha_g for m=1 
################################
# a%*%x=b
wt1.1=1
p1d=cbind(1,ind1,x0)
MxA=(t(p1d)*matrix(wt1.1,ncol=n,nrow=dim(x0)[2]+2,byrow=TRUE))%*%p1d/n
MxB=t(p1d)%*%as.matrix(y*wt1.1)/n
sol1=solve(MxA,MxB)


# variance minimization
wt.i1=ind1*(1-1/pred.lambda.1)


p1d=cbind(1,ind1,x0)
MxA.i=(t(p1d)*matrix(wt.i1,ncol=n,nrow=dim(x0)[2]+2,byrow=TRUE))%*%p1d/n
MxB.i=t(p1d)%*%as.matrix(y*wt.i1)/n
sol1.t=solve(MxA.i[2:4,2:4],MxB.i[2:4])
sol1.i=c(0,sol1.t)

#cT adjustment
#derivative of g1.1
dg11=cbind(ind1,x0)
	
#d(lambda.1)/d(gamma)=lambda.1^2*exp(-gammaX)  .X
dlam.nox=pred.lambda.1^2*exp(-a[1]-x[,1]*bxt4-x0[,1]*b[1]-x0[,2]*b[2])
dx=cbind(1,x[,1],x0)

theta2=c(sol1.t,rep(0.1, dim(dx)[2]))
	iter <- 0
	#itermax=150
	tolerance=0.000001  
theta=0
	#---Newton Ralphson-------------------------------
	while (sum(abs(theta - theta2)) >= tolerance) {	
		#print(abs(theta - theta2)[4]) 
		iter=iter+1
		#print(paste("iter =",iter))
	
		theta=theta2
		c=theta[(length(sol1.t)+1):length(theta)]
		ag=theta[1:length(sol1.t)]
g11=dg11%*%ag
l1=(dlam.nox/(1-pred.lambda.1)*(dx%*%c))

Sa=t(ind1*(1/pred.lambda.1-1)*(y-g11-l1))%*% dg11
Sb=t(ind1/pred.lambda.1*(y-g11-l1)*dlam.nox)%*%dx
S=cbind(Sa,Sb)

#-I(theta)
Iaa=t(matrix(ind1*(1/pred.lambda.1-1),ncol=dim(dg11)[2],nrow=n)*dg11)%*%dg11
Ica=t(matrix(ind1/pred.lambda.1*dlam.nox,ncol=dim(dx)[2],nrow=n)*dx)%*%dg11
Icc=t(matrix(ind1*dlam.nox^2/((1-pred.lambda.1)*pred.lambda.1),
		ncol=dim(dx)[2],nrow=n)*(dx))%*%dx
I=rbind(cbind(Iaa,t(Ica)),cbind(Ica, Icc))

theta2=theta+solve(I)%*%t(S)

	#print(theta2)

	}




############################## 
# solving alpha_g for m=2 
################################
wt1=1
wt2=1#1-ind3

q1d=cbind(1,ind1,ind2,x0) # 1st derivative for g2^(2)
prc2=(1-pred.lambda0.1)

xlast=pred.xt1*ind1+ x[,1]*(1-ind1)
q2d=cbind(1,ind1,ind1*prc2,x0)  # 1st derivative for g1^(2)
# 1st derivative for sum p(a|s)*g2^(2)-g1^(2)
q2p=cbind(alpha0=0,alpha1=0,alpha2=ind1*(ind2-prc2),alphab1=0,alphab2=0)  

Mxb=t(q1d)%*%as.matrix(y*wt1)/n
Mxa=((t(q1d)*matrix(wt1,ncol=n,nrow=dim(x0)[2]+3,byrow=TRUE))
%*%q1d-(t(q2d)*matrix(wt2,ncol=n,nrow=dim(x0)[2]+3,byrow=TRUE))%*%q2p)/n

sol2=solve(Mxa,Mxb)



# variance minimization
wt.i2=ind1*(pred.lambda.1/(1-pred.lambda.1))


Mxb.i=t(q1d)%*%as.matrix(y*wt.i2)/n
Mxa.i=((t(q1d)*matrix(wt.i2,ncol=n,nrow=dim(x0)[2]+3,byrow=TRUE))
%*%q1d-(t(q2d)*matrix(0,ncol=n,nrow=dim(x0)[2]+3,byrow=TRUE))%*%q2p)/n

sol2t=solve(Mxa.i[2:5,2:5],Mxb.i[2:5])
sol2.i=c(0,sol2t)

#cT adjustment
#derivae of g1.2
#dg120=cbind(ind1,x0,xlast,ind1*prc2)
dg12=cbind(ind1,ind1*prc2,x0)
#	dg12=cbind(ind1,ind2,x0,xx)
#d(lambda.1)/d(gamma)=lambda.1^2*exp(-gammaX)  .X
dlam.nox2=(pred.lambda.1)^2*exp(-a[1]-x[,1]*bxt4-x0[,1]*b[1]-x0[,2]*b[2])
dx=cbind(1,x[,1],x0)

theta4=c(sol2t,rep(-0.1, dim(dx)[2]))
	iter <- 0
	#itermax=150
	tolerance=0.000001  
theta3=0
	#---Newton Ralphson-------------------------------
	while (sum(abs(theta3 - theta4)) >= tolerance) {	
		#print(abs(theta3 - theta4)[4]) 
		iter=iter+1
		#print(paste("iter =",iter))
	
		theta3=theta4
		c2=theta3[(length(sol2t)+1):length(theta3)]
		ag2=theta3[1:length(sol2t)]
g12=dg12%*%ag2
l2=(dlam.nox2/(pred.lambda.1)*(dx%*%c2))

Sa2=t(ind1*(pred.lambda.1/((1-pred.lambda.1)))*(y-g12-l2))%*% dg12
Sb2=t(ind1/((1-pred.lambda.1))*(y-g12-l2)*dlam.nox2)%*%dx
S2=cbind(Sa2,Sb2)

#-I(theta)
Iaa2=t(matrix(ind1*(pred.lambda.1/((1-pred.lambda.1))),
			ncol=dim(dg12)[2],nrow=n)*dg12)%*%dg12
Ica2=t(matrix(ind1/((1-pred.lambda.1))*dlam.nox2,ncol=dim(dx)[2],
			nrow=n)*dx)%*%dg12
Icc2=t(matrix(ind1*dlam.nox2^2/((1-pred.lambda.1)*pred.lambda.1),
			ncol=dim(dx)[2],nrow=n)*(dx))%*%dx
I2=rbind(cbind(Iaa2,t(Ica2)),cbind(Ica2, Icc2))

theta4=theta3+solve(I2)%*%t(S2)

	#print(theta4)

	}





####### fit in g funtions########################


g1.1=p1d%*%as.matrix(sol2[-3])
g2.2=q1d%*%as.matrix(sol2)
g1.2=q2d%*%as.matrix(sol2)
g1.1i=p1d%*%as.matrix(sol1.i)
g1.2i=q2d%*%as.matrix(sol2.i)

g1.1c=dg11%*%as.matrix(theta2[1:length(sol1.t)])
g1.2c=cbind(ind1,ind1*prc2,x0)%*%as.matrix(theta4[1:length(sol2t)])



mu.dr=mu.jt=mu.or=mu.or.i=mu.dr.i=mu.dr.c=rep(NA,length(L))
#------------DR m=1-------------------------
mu.dr[1]=sum(ft1$W*y-ft1$W*g1.1+g1.1)/n
mu.jt[1]=sum(ft1$W*y)/sum(ft1$W)
mu.or[1]=mean(g1.1)
mu.or.i[1]=mean(g1.1i)
mu.dr.i[1]=sum(ft1$W*y-ft1$W*g1.1i+g1.1i)/n
mu.dr.c[1]=sum(ft1$W*y-ft1$W*g1.1c+g1.1c)/n


#------------DR m=2-------------------------
t.num.d=sum((ft2$W*y-(ft2$W-1)*g1.2))
t.num.di=sum((ft2$W*y-(ft2$W-1)*g1.2i))	#(1-ind3)*
t.num.dc=sum((ft2$W*y-(ft2$W-1)*g1.2c))

t.den=n#sum(1-ind3)
mu.jt[2]=sum(ft2$W*y)/sum(ft2$W)
mu.or[2]=mean(g1.2)
mu.dr[2]=t.num.d/t.den
mu.dr.i[2]=t.num.di/t.den
mu.or.i[2]=mean(g1.2i)
mu.dr.c[2]=t.num.dc/t.den

return(list(mu.dr=mu.dr,mu.jt=mu.jt,mu.or=mu.or,
mu.or.i=mu.or.i, mu.dr.i=mu.dr.i,mu.dr.c=mu.dr.c))
}

#revised 5/2/2015
# calculate weight Wpk for the IPW estimator as well as 
#the intermediate weight Wpt, t<k=m

calc.w <- function (u,d,L,M,x0,time.x,Lx,tmpx){
	
	p=length(L)
	n=length(u)
	p0=dim(x0)[2]
	#------------estimation of parameters-------------------------
	ft4=lucy.discrete.haz.dep(u=u,d=d,L=L,x0=x0,time.x=time.x,Lx=Lx)
	a=ft4$beta[1:(p-1)]
	bxt4=ft4$beta[p+p0]
	b=ft4$beta[(p):(p+p0-1)]

	# ---calculate Lambda from discrete hazard function with 
	#respect to physician stop-----
	pred.lambda=matrix(NA,p,n)
	pred.lambda2=matrix(0,p,n)
	

	for (j in 1:(p)){
		pred.lambda[j,]<- 1/(1+exp(-a[j]-tmpx[,j]*bxt4-x0[,1]*b[1]-x0[,2]*b[2]))
	}

	# lambda4=1 since it is the end of the study.
	pred.lambda[p,]<-1

	K=p
	#---calculate denominator of W-------------------------------------------
	ind1= matrix(u,n,K)==matrix(L,n,K,byrow=TRUE) & matrix(d,n,K)==1 
	#I(U=tj, d=1)
	ind2=	matrix(u,n,K)> matrix(L,n,K,byrow=TRUE) 				#I(U>tj)

		# denominator for each element by individual
	den=(t(1-pred.lambda))^ind2*(t(pred.lambda)^ind1)
	
	prodden=apply(den,1,prod) #denominator for W by individual

	#-------calculate numerator of W-----------------------
	num=matrix(NA,n,K)
	num= (matrix(L,n,K,byrow=TRUE)==L[M])^ind1*(matrix(u,n,K)<=L[M])^ind2
		# tj=tm or u<=tm

	prodnum=apply(num,1,prod) #numerator for W by individual

	##########-W-#############
	W=prodnum/prodden

	#-----------------calc Wt---------------------
	tmpden=den
	tmpnum=num

	for (j in 2:p){
		tmpden[,j]=apply(den[,1:j],1,prod)
		tmpnum[,j]=apply(num[,1:j],1, prod)
	}
	Wt=tmpnum/tmpden

	######calc W{t-1}##########
	Wt1=matrix(1,n,p)
	Wt1[,-1]=Wt[,-p]


return(list(W=W,Wt=Wt,Wt1=Wt1,num=num, den=den, tmpnum=tmpnum, tmpden=tmpden))

}



\end{verbatim}
