bound.ps<-function(y,z,s,y.bin=T,r1,r0,lambda1=1,lambda0=1,probit=T,tau=0.5,prob.se=T,beta.cov,x) { y1<-y[z==1] y0<-y[z==0] rho<-mean(z) mu11<-mean(y1) mu00<-mean(y0) phi11<-rho*mu11 phi00<-(1-rho)*mu00 eta11<-mean(y1^2) eta00<-mean(y0^2) dta11<-rho*eta11 dta00<-(1-rho)*eta00 # Point estimate by IPW theta1<-max(min(mean(y*z/s),phi11+1-rho),phi11) theta0<-max(min(mean(y*(1-z)/(1-s)),phi00+rho),phi00) att.est<-mu11-(theta0-phi00)/rho atc.est<-(theta1-phi11)/(1-rho)-mu00 ate.est<-theta1-theta0 s.var.1<-mean((1-s)/s) s.var.0<-mean(s/(1-s)) if (y.bin==F) { b.1<-phi11-dta11 b.0<-phi00-dta00 d.1<-rho*mu11^2+(1-rho)*eta11 d.0<-(1-rho)*mu00^2+rho*eta00 f.1.star<-(theta1^2-2*mu11*theta1+d.1)/(d.1-mu11^2) f.0.star<-(theta0^2-2*mu00*theta0+d.0)/(d.0-mu00^2) obs1<-c(rho=rho,mu=mu11,phi=phi11,eta=eta11,dta=dta11,theta=theta1,s.var=s.var.1,f.star=f.1.star,b=b.1,d=d.1) obs0<-c(rho=1-rho,mu=mu00,phi=phi00,eta=eta00,dta=dta00,theta=theta0,s.var=s.var.0,f.star=f.0.star,b=b.0,d=d.0) } else { obs1<-c(rho=rho,mu=mu11,phi=phi11,eta=eta11,dta=dta11,theta=theta1,s.var=s.var.1,b=0) obs0<-c(rho=1-rho,mu=mu00,phi=phi00,eta=eta00,dta=dta00,theta=theta0,s.var=s.var.0,b=0) } # Calculation of first bound for Y1 and Y0 bound.1.1<-first.bound(par=obs1,a=r1^2*lambda1*s.var.1,bin=y.bin,r=r1) bound.1.0<-first.bound(par=obs0,a=r0^2*lambda0*s.var.0,bin=y.bin,r=r0) # Calculation of second bound for Y1 and Y0 (when continous variable) bound.2.1<-bound.2.0<-NA if (y.bin==F) { bound.2.1<-second.bound(par=obs1,a=r1^2*lambda1*s.var.1,r=r1) bound.2.0<-second.bound(par=obs0,a=r0^2*lambda0*s.var.0,r=r0) } # Calculation of y limits y1.lim<-comb.bounds(first=bound.1.1,second=bound.2.1,par=obs1,r=r1,bin=y.bin) y0.lim<-comb.bounds(first=bound.1.0,second=bound.2.0,par=obs0,r=r0,bin=y.bin) # Calculation of ATT, ATC and ATE att.lim<-atc.lim<-ate.lim<-c(min=NA,max=NA) if (!is.na(y0.lim["min"])) att.lim["max"]<-mu11-(y0.lim["min"]-phi00)/rho if (!is.na(y0.lim["max"])) att.lim["min"]<-mu11-(y0.lim["max"]-phi00)/rho if (!is.na(y1.lim["max"])) atc.lim["max"]<-(y1.lim["max"]-phi11)/(1-rho)-mu00 if (!is.na(y1.lim["min"])) atc.lim["min"]<-(y1.lim["min"]-phi11)/(1-rho)-mu00 if (!is.na(y1.lim["max"]) & !is.na(y0.lim["min"])) ate.lim["max"]<-y1.lim["max"]-y0.lim["min"] if (!is.na(y1.lim["min"]) & !is.na(y0.lim["max"])) ate.lim["min"]<-y1.lim["min"]-y0.lim["max"] # Probit model if (probit) { cov.matr<-matrix(c(1+tau^2,tau^2,tau^2,1+tau^2),2,2) w<-sqrt(1+tau^2) lin.pre<-qnorm(s) if (tau>0) var.s.star<-apply(cbind(w*lin.pre,w*lin.pre),1,pmvnorm,lower=-Inf,sigma=cov.matr)-s^2 else var.s.star<-rep(0,length(s)) g<-mean(var.s.star/s^2) h<-mean(var.s.star/(1-s)^2) bound.1.1.pb<-first.bound(par=obs1,a=r1^2*g,bin=y.bin,r=r1,se=prob.se) bound.1.0.pb<-first.bound(par=obs0,a=r0^2*h,bin=y.bin,r=r0,se=prob.se) bount.2.1.pb<-bound.2.0.pb<-NA if (y.bin==F) { bound.2.1.pb<-second.bound(par=obs1,a=r1^2*g,r=r1,se=prob.se) bound.2.0.pb<-second.bound(par=obs0,a=r0^2*h,r=r0,se=prob.se) } y1.lim.pb<-comb.bounds(first=bound.1.1.pb,second=bound.2.1.pb,par=obs1,r=r1,bin=y.bin,se=prob.se) y0.lim.pb<-comb.bounds(first=bound.1.0.pb,second=bound.2.0.pb,par=obs0,r=r0,bin=y.bin,se=prob.se) att.lim.pb<-atc.lim.pb<-ate.lim.pb<-c(min=NA,max=NA) if (!is.na(y0.lim.pb["min"])) att.lim.pb["max"]<-mu11-(y0.lim.pb["min"]-phi00)/rho if (!is.na(y0.lim.pb["max"])) att.lim.pb["min"]<-mu11-(y0.lim.pb["max"]-phi00)/rho if (!is.na(y1.lim.pb["max"])) atc.lim.pb["max"]<-(y1.lim.pb["max"]-phi11)/(1-rho)-mu00 if (!is.na(y1.lim.pb["min"])) atc.lim.pb["min"]<-(y1.lim.pb["min"]-phi11)/(1-rho)-mu00 if (!is.na(y1.lim.pb["max"]) & !is.na(y0.lim.pb["min"])) ate.lim.pb["max"]<-y1.lim.pb["max"]-y0.lim.pb["min"] if (!is.na(y1.lim.pb["min"]) & !is.na(y0.lim.pb["max"])) ate.lim.pb["min"]<-y1.lim.pb["min"]-y0.lim.pb["max"] if (prob.se) { if (y.bin) { y1.max.d<-c(y1.lim.pb["d.max.d.theta"],0,y1.lim.pb["d.max.d.rho"],y1.lim.pb["d.max.d.phi"],0,y1.lim.pb["d.max.d.v"],0) y1.min.d<-c(y1.lim.pb["d.min.d.theta"],0,y1.lim.pb["d.min.d.rho"],y1.lim.pb["d.min.d.phi"],0,y1.lim.pb["d.min.d.v"],0) y0.max.d<-c(0,y0.lim.pb["d.max.d.theta"],-y0.lim.pb["d.max.d.rho"],0,y0.lim.pb["d.max.d.phi"],0,y0.lim.pb["d.max.d.v"]) y0.min.d<-c(0,y0.lim.pb["d.min.d.theta"],-y0.lim.pb["d.min.d.rho"],0,y0.lim.pb["d.min.d.phi"],0,y0.lim.pb["d.min.d.v"]) phi11.d<-c(0,0,0,1,0,0,0) phi00.d<-c(0,0,0,0,1,0,0) rho.d<-c(0,0,1,0,0,0,0) } else { y1.max.d<-c(y1.lim.pb["d.max.d.theta"],0,y1.lim.pb["d.max.d.rho"],y1.lim.pb["d.max.d.phi"],y1.lim.pb["d.max.d.dta"],0,0,y1.lim.pb["d.max.d.v"],0) y1.min.d<-c(y1.lim.pb["d.min.d.theta"],0,y1.lim.pb["d.min.d.rho"],y1.lim.pb["d.min.d.phi"],y1.lim.pb["d.min.d.dta"],0,0,y1.lim.pb["d.min.d.v"],0) y0.max.d<-c(0,y0.lim.pb["d.max.d.theta"],-y0.lim.pb["d.max.d.rho"],0,0,y0.lim.pb["d.max.d.phi"],y0.lim.pb["d.max.d.dta"],0,y0.lim.pb["d.max.d.v"]) y0.min.d<-c(0,y0.lim.pb["d.min.d.theta"],-y0.lim.pb["d.min.d.rho"],0,0,y0.lim.pb["d.min.d.phi"],y0.lim.pb["d.min.d.dta"],0,y0.lim.pb["d.min.d.v"]) phi11.d<-c(0,0,0,1,0,0,0,0,0) phi00.d<-c(0,0,0,0,0,1,0,0,0) rho.d<-c(0,0,1,0,0,0,0,0,0) } att.max.d<-((phi11.d+phi00.d-y0.min.d)*rho-rho.d*(phi11+phi00-y0.lim.pb["min"]))/rho^2 att.min.d<-((phi11.d+phi00.d-y0.max.d)*rho-rho.d*(phi11+phi00-y0.lim.pb["max"]))/rho^2 atc.max.d<-((y1.max.d-phi11.d-phi00.d)*(1-rho)+rho.d*(y1.lim.pb["max"]-phi11-phi00))/(1-rho)^2 atc.min.d<-((y1.min.d-phi11.d-phi00.d)*(1-rho)+rho.d*(y1.lim.pb["min"]-phi11-phi00))/(1-rho)^2 ate.max.d<-y1.max.d-y0.min.d ate.min.d<-y1.min.d-y0.max.d den<-dnorm(lin.pre) p.nor<-pnorm(lin.pre/sqrt(2*w^2-1)) influ<-(((z/s-(1-z)/(1-s))*den)*x)%*%(beta.cov*length(s)) alpha1<-apply(-(den*y*z/s^2)*x,2,mean) alpha0<-apply((den*y*(1-z)/(1-s)^2)*x,2,mean) gamma1<-apply((2*den*(p.nor*s-var.s.star-s^2)/s^3)*x,2,mean) gamma0<-apply((2*den*(p.nor+var.s.star+s^2-s-s*p.nor)/(1-s)^3)*x,2,mean) if (y.bin) cov.ran<-cov(cbind(y*z/s+influ%*%alpha1,y*(1-z)/(1-s)+influ%*%alpha0,z,y*z,y*(1-z), var.s.star/s^2+influ%*%gamma1,var.s.star/(1-s)^2+influ%*%gamma0)) else cov.ran<-cov(cbind(y*z/s+influ%*%alpha1,y*(1-z)/(1-s)+influ%*%alpha0,z,y*z,y^2*z,y*(1-z), y^2*(1-z),var.s.star/s^2+influ%*%gamma1,var.s.star/(1-s)^2+influ%*%gamma0)) d.matr<-rbind(y1.max.d,y1.min.d,y0.max.d,y0.min.d,att.max.d,att.min.d,atc.max.d,atc.min.d,ate.max.d,ate.min.d) se.pb<-sqrt(diag(d.matr%*%cov.ran%*%t(d.matr))/length(s)) } } out1<-rbind(c(theta1,theta0,att.est,atc.est,ate.est),cbind(y1.lim,y0.lim,att.lim,atc.lim,ate.lim)) dimnames(out1)<-list(c("est","min","max"),c("y1.mean","y0.mean","att","atc","ate")) out2<-c(obs1,obs0[-1]) names(out2)<-c(paste(names(obs1),".1",sep=""),paste(names(obs0)[-1],".0",sep="")) out3<-c(r1=r1,r0=r0,lambda1=lambda1,lambda0=lambda0) if (probit) { out1<-rbind(out1,cbind(y1.lim.pb,y0.lim.pb,att.lim.pb,atc.lim.pb,ate.lim.pb)[1:2,]) out2<-c(out2,g=g,h=h) out3<-c(out3,tau=tau) rownames(out1)<-c("est","min","max","min.pb","max.pb") if (prob.se) { out1<-rbind(out1,se.pb[seq(2,length(se.pb),2)],se.pb[seq(1,length(se.pb)-1,2)]) rownames(out1)<-c("est","min","max","min.pb","max.pb","se.min.pb","se.max.pb") } } return(list(result=out1,distri=out2,parameter=out3)) } first.bound<-function(par,a,bin,r,se=F) { b<-par["b"] theta<-par["theta"] phi<-par["phi"] rho<-par["rho"] sq<-sqrt(a^2*(1-4*b)+4*a*(theta*(1-theta)-b)) up<-(2*theta+a+sq)/(2*(a+1)) low<-(2*theta+a-sq)/(2*(a+1)) if (bin==F) { up<-min(up,phi+1-rho) low<-max(low,phi) } else { if (up>phi+1-rho) up<-NA if (low0) sq.d<-(0.5/sq)*(2*a*a.d*(1-4*b)-4*b.d*a^2+4*(a.d*(theta*(1-theta)-b)+((1-2*theta)*theta.d-b.d)*a)) else sq.d<-rep(0,length(a.d)) up.d<-((2*theta.d+a.d+sq.d)*2*(a+1)-2*a.d*(2*theta+a+sq))/(4*(a+1)^2) low.d<-((2*theta.d+a.d-sq.d)*2*(a+1)-2*a.d*(2*theta+a-sq))/(4*(a+1)^2) if (up==phi+1-rho) up.d<-c(0,-1,1,0,0) if (low==phi) low.d<-c(0,0,1,0,0) names(up.d)<-paste("d.up.d.",c("theta","rho","phi","dta","v"),sep="") names(low.d)<-paste("d.low.d.",c("theta","rho","phi","dta","v"),sep="") } else { up.d<-low.d<-rep(NA,4) a.d<-c(0,0,0,r^2) b.d<-rep(0,4) theta.d<-c(1,0,0,0) if (a>0) sq.d<-(0.5/sq)*(2*a*a.d*(1-4*b)-4*b.d*a^2+4*(a.d*(theta*(1-theta)-b)+((1-2*theta)*theta.d-b.d)*a)) else sq.d<-rep(0,length(a.d)) if (!is.na(up)) { if (upphi) low.d<-((2*theta.d+a.d-sq.d)*2*(a+1)-2*a.d*(2*theta+a-sq))/(4*(a+1)^2) else low.d<-c(0,0,1,0) } names(up.d)<-paste("d.up.d.",c("theta","rho","phi","v"),sep="") names(low.d)<-paste("d.low.d.",c("theta","rho","phi","v"),sep="") } out<-c(out,low.d,up.d) } return(out) } second.bound<-function(par,a,r,se=F) { theta<-par["theta"] mu<-par["mu"] d<-par["d"] rho<-par["rho"] phi<-par["phi"] dta<-par["dta"] f<-a*rho/(1-rho) f.star<-par["f.star"] if (f<=f.star) { D<-sqrt(f*((mu^2-d)*f+theta^2-2*mu*theta+d)) if (f>1) { right<-(f*mu-theta+D)/(f-1) left<-(f*mu-theta-D)/(f-1) } else if (f==1) { if (theta!=mu) left<-right<-(theta^2-d)/(2*(theta-mu)) else { if (d<=theta^2) left<-right<-theta else left<-right<-NA } } else if (f<1) { right<-(f*mu-theta-D)/(f-1) left<- (f*mu-theta+D)/(f-1) } } else left<-right<-NA out<-c(left,right,f) names(out)<-c("left","right","f") if (se==T) { f.d<-mu.d<-theta.d<-d.d<-D.d<-left.d<-right.d<-c(theta=NA,rho=NA,phi=NA,dta=NA,v=NA) if (!is.na(right)) { mu.d<-c(0,-mu/rho,1/rho,0,0) theta.d<-c(1,0,0,0,0) d.d<-c(0,-(phi^2+dta)/rho^2,2*phi/rho,(1-rho)/rho,0) if (f != 1) { f.d<-c(0,f/(rho*(1-rho)),0,0,r^2*rho/(1-rho)) if (a>0) D.d<-(0.5/D)*(f.d*((mu^2-d)*f+theta^2-2*mu*theta+d)+((2*mu*mu.d-d.d)*f+(mu^2-d)*f.d+2*theta*theta.d-2*mu.d*theta -2*mu*theta.d+d.d)*f) else D.d<-rep(0,length(mu.d)) right.d<-((f.d*mu+f*mu.d-theta.d+sign(f-1)*D.d)*(f-1)-f.d*(f*mu-theta+sign(f-1)*D))/(f-1)^2 left.d<-((f.d*mu+f*mu.d-theta.d-sign(f-1)*D.d)*(f-1)-f.d*(f*mu-theta-sign(f-1)*D))/(f-1)^2 } else if (f==1 & theta != mu) right.d<-left.d<-0.5*((2*theta*theta.d-d.d)*(theta-mu)-(theta.d-mu.d)*(theta^2-d))/(theta-mu)^2 else if (f==1 & theta==mu) right.d<-left.d<-theta.d } names(right.d)<-paste("d.right.d.",c("theta","rho","phi","dta","v"),sep="") names(left.d)<-paste("d.left.d.",c("theta","rho","phi","dta","v"),sep="") out<-c(out,left.d,right.d) } return(out) } comb.bounds<-function(first,second=NA,par,bin,r,se=F) { if (bin) { y.max<-y.min<-(r>=0)*first["low"]+(r<0)*first["up"] out<-c(y.min,y.max) names(out)<-c("min","max") if (se) { up.d<-paste("d.up.d.",c("theta","rho","phi","v"),sep="") low.d<-paste("d.low.d.",c("theta","rho","phi","v"),sep="") y.max.d<-y.min.d<-(r>=0)*first[low.d]+(r<0)*first[up.d] out<-c(out,y.min.d,y.max.d) names(out)<-c("min","max",paste("d.min.d.",c("theta","rho","phi","v"),sep=""),paste("d.max.d.",c("theta","rho","phi","v"),sep="")) } } else { y.min<-y.max<-NA m<-0 out<-rep(NA,2) if (!is.na(second["left"])) { if (second["f"]<1) { if (r>=0 & first["low"]<=second["left"]) { m<-1 y.min<-first["low"] y.max<-second["left"] } else if (r<0 & first["up"]>=second["right"]) { m<-2 y.min<-second["right"] y.max<-first["up"] } } else if (second["f"]==1) { if (par["mu"]<=par["theta"] & r>=0 & first["low"]<=second["left"]) { m<-3 y.min<-first["low"] y.max<-second["left"] } else if (par["mu"]>par["theta"] & r<0 & first["up"]>=second["right"]) { m<-4 y.min<-second["right"] y.max<-first["up"] } } else if (second["f"]>1) { if (par["mu"]<=par["theta"] & r>=0 & first["low"]<=second["right"]) { m<-5 y.min<-max(first["low"],second["left"]) y.max<-second["right"] } else if (par["mu"]>par["theta"] & r<0 & first["up"]>=second["left"]) { m<-6 y.min<-second["left"] y.max<-min(first["up"],second["right"]) } } out<-c(y.min,y.max) names(out)<-c("min","max") if (se) { up.d<-paste("d.up.d.",c("theta","rho","phi","dta","v"),sep="") low.d<-paste("d.low.d.",c("theta","rho","phi","dta","v"),sep="") right.d<-paste("d.right.d.",c("theta","rho","phi","dta","v"),sep="") left.d<-paste("d.left.d.",c("theta","rho","phi","dta","v"),sep="") if (m>0) { y.min.d<-(m==1 | m==3)*first[low.d]+(m==2 | m==4)*second[right.d]+(m==6)*second[left.d]+ (m==5)*((first["low"]>=second["left"])*first[low.d]+(first["low"]second["right"])*second[right.d]) } else y.min.d<-y.max.d<-rep(NA,5) out<-c(out,y.min.d,y.max.d) names(out)<-c("min","max",paste("d.min.d.",c("theta","rho","phi","dta","v"),sep=""),paste("d.max.d.",c("theta","rho","phi","dta","v"),sep="")) } } } return(out) } svar.lim<-function(s,lambda) { s<-sort(s) v1<-cumsum((1-s)/s) v0<-cumsum(s/(1-s)) v.fix<-v1[length(v1)]*lambda temp<-which(v10) { pos.min<-max(temp) resi<-v.fix-v1[pos.min] lambda.min<-(v0[pos.min]+resi*s[pos.min+1]^2/(1-s[pos.min+1])^2)/v0[length(v0)] } else { resi<-v.fix lambda.min<-(resi*s[1]^2/(1-s[1])^2)/v0[length(v0)] } temp<-which(v1[length(v1)]-v1