###### This files contains the constant, linear, and quadratic sensitivity functions mentioned in the paper ##### ###Written by Lingling Li, last updated on March 03, 2014 ### ### Type I sensitivity function that c(t=1,ps) = c, and c(t=0,ps) = c*cr, both c and cr are constants ### ##input parameters: #--- y: the vector of outcomes #--- t: the vector of treatment variables #--- ps: the vector of propensity scores #--- c: c(t=1,ps)=c of the sensitivity function; c=0 means no unmeasured confounding #--- cr: c(t=0,ps)/c(t=1,ps), alternatively, c(t=0,ps)=cr*c ## output: #--- out: an 3*1 matrix for the Inverse probalbity weighted (IPW) estimates of risk (binary outcomes) or mean (continuous outcomes) difference, ######### relative risk (binary outcomes only), and odds ratio (binary outcomes only) sa1 = function(y, t, ps, c, cr=1) { c1 = c; c0 = cr*c; y = y - (t*(1-ps)*c1 + (1-t)*(-ps)*c0) mu1 = mean(t*y/ps)/mean(t/ps); mu0 = mean((1-t)*y/(1-ps))/mean((1-t)/(1-ps)); #mu1 = mean(t*y/ps); #mu0 = mean((1-t)*y/(1-ps)); rd.ipw = mu1 - mu0; rr.ipw = mu1/mu0; or.ipw = mu1/(1-mu1)*(1-mu0)/mu0; out = matrix(c(rd.ipw, rr.ipw, or.ipw), ncol=1); return(out); } ### Type II sensitivity function that c(t=1,ps)=c+s*ps, and c(t=0,ps)=c*cr+s*sr*ps, both c,s,cr,and sr are constants ### ##input parameters: #--- y: the vector of outcomes #--- t: the vector of treatment variables #--- ps: the vector of propensity scores #--- c and s: c(t=1,ps)=c+s*ps, c and s are intercept and slope for c(t=1, ps) #--- cr and sr: c(t=0,ps)=c*cr+s*sr*ps; cr is the ratio between the two intercepts (for t=0 and t=1) and sr is the ratio between the two slopes (for t=0 and t=1) ## output: #--- out: an 3*1 matrix for the Inverse probability weighted (IPW) estimates of risk (binary outcomes) or mean (continuous outcomes) difference, ######### relative risk (binary outcomes only), and odds ratio (binary outcomes only) sa2 = function(y, t, ps, c, cr=1, s, sr=1) { #--- y: the vector of outcomes #--- t: the vector of treatment indicators #--- ps: the vector of propensity scores # the Sensitivity analysis with a linear SF c0 = cr*c; c1 = c; s0 = sr*s; s1 = s; y1 = y - (t*(1-ps)*(c1+s1*ps) + (1-t)*(-ps)*(c0+s0*ps) ) mu1 = mean(t*y1/ps)/mean(t/ps); mu0 = mean((1-t)*y1/(1-ps))/mean((1-t)/(1-ps)); rd.ipw = mu1 - mu0; or.ipw = mu1/(1-mu1)*(1-mu0)/mu0; rr.ipw = mu1/mu0; out = matrix(c(rd.ipw, rr.ipw, or.ipw), ncol=1); return(out) } ### Type III sensitivity function that c(t=1,ps)=c1+s1*ps+q1*ps*ps, and c(t=0,ps)=c0+s0*ps+q0*ps*ps ### ##input parameters: #--- y: the vector of outcomes #--- t: the vector of treatment variables #--- ps: the vector of propensity scores #--- (c1,s1, q1): c(t=1,ps)=c1+s1*ps+q1*ps*ps #--- (c0,s0, q0): c(t=0,ps)=c0+s0*ps+q0*ps*ps ## output: #--- out: an 3*1 matrix for the Inverse probability weighted (IPW) estimates of risk (binary outcomes) or mean (continuous outcomes) difference, ######### relative risk (binary outcomes only), and odds ratio (binary outcomes only) sa3 = function(y, t, ps, c0, c1, s0, s1, q0, q1) { #--- y: the vector of outcomes #--- t: the vector of treatment indicators #--- ps: the vector of propensity scores #--- c(1,e) = c1 + s1*e + q1*e^2 #--- c(0,e) = c0 + s0*e + q0*e^2 y = y - (t*(1-ps)*(c1+s1*ps+q1*ps*ps) + (1-t)*(-ps)*(c0+s0*ps+q0*ps*ps) ) mu1 = mean(t*y/ps)/mean(t/ps); mu0 = mean((1-t)*y/(1-ps))/mean((1-t)/(1-ps)); rd.ipw = mu1 - mu0; or.ipw = mu1/(1-mu1)*(1-mu0)/mu0; rr.ipw = mu1/mu0; out = matrix(c(rd.ipw, rr.ipw, or.ipw), ncol=1); return(out); } ### Implement two specific Type III (quadratic) sensitivity functions given the values of ### c(t=1,ps=ps.min), c(t=1,ps=ps.max), c(t=0,ps=ps.min), c(t=0,ps=ps.max), where ps.min and ps.max are ### the lower and upper bounds of the PS in your dataset. Note that ps.min>=0 and ps.max<=1. #--- y: the vector of outcomes #--- t: the vector of treatment variables #--- ps: the vector of propensity scores #--- ps.min: the lower bound of the PS in your dataset #--- ps.max: the upper bound of the PS in your dataset #--- sf1.min: c(t=1,ps=ps.min) #--- sf1.max: c(t=1,ps=ps.max) #--- sf0.min: c(t=0,ps=ps.min) #--- sf1.max: c(t=0,ps=ps.max) #--- this function implement two Type III sensitivity functions with the same starting and ending points (as defined above) #--- SF 1: c(t, ps=((ps.min+ps.max)/2)=c(t,ps=ps.min)+(c(t,ps=ps.max)-c(t,ps=ps.min))/4 #--- SF 2: c(t, ps=((ps.min+ps.max)/2)=c(t,ps=ps.min)+(c(t,ps=ps.min)-c(t,ps=ps.min))*3/4 ## output: #--- out: an 6*1 vector, #--- the first two elements for the two IPW risk or mean difference estimates for SF1 and SF2 respectively #--- the next two elements for the two RR estimates for SF1 and SF2 respectively #--- the last two elements for the two OR estimates for SF1 and SF2 respectively sa3_2 = function(y, t, ps, ps.min=0, ps.max=1, sf1.min, sf1.max, sf0.min, sf0.max) { # SF1 # c(t=0,ps)=a0+b0*ps+q0*ps*ps # c(t=1,ps)=a1+b1*ps+q1*ps*ps # (a0, b0, q0) are determined by (i) c(t=0, ps=ps.min)=sf0.min, (ii) c(t=0, ps=ps.max)=sf0.max, # and (iii) c(t=0, ps=((ps.min+ps.max)/2)=c(t,ps=ps.min)+(c(t,ps=ps.max)-c(t,ps=ps.min))/4 # (a1, b1, q1) are determined similarly tem=c(1, ps.min, ps.min*ps.min, 1, ps.max, ps.max*ps.max, 1, (ps.min+ps.max)/2, ((ps.min+ps.max)/2)^2) A=matrix(tem, nrow=3, ncol=3, byrow=T) vecb0=matrix(c(sf0.min, sf0.max, (sf0.min+(sf0.max-sf0.min)/4)), ncol=1) vecb1=matrix(c(sf1.min, sf1.max, (sf1.min+(sf1.max-sf1.min)/4)), ncol=1) tem2=solve(A, vecb0); tem3=solve(A, vecb1) a0=tem2[1]; b0=tem2[2]; q0=tem2[3]; a1=tem3[1]; b1=tem3[2]; q1=tem3[3]; out1=sa3(y, t, ps, a0, a1, b0, b1, q0, q1) # SF2 # c(t=0,ps)=a0+b0*ps+q0*ps*ps # c(t=1,ps)=a1+b1*ps+q1*ps*ps # (a0, b0, q0) are determined by (i) c(t=0, ps=ps.min)=sf0.min, (ii) c(t=0, ps=ps.max)=sf0.max, # and (iii) c(t=0, ps=((ps.min+ps.max)/2)=c(t,ps=ps.min)+(c(t,ps=ps.max)-c(t,ps=ps.min))*3/4 # (a1, b1, q1) are determined similarly vecb0=matrix(c(sf0.min, sf0.max, (sf0.min+(sf0.max-sf0.min)*3/4)), ncol=1) vecb1=matrix(c(sf1.min, sf1.max, (sf1.min+(sf1.max-sf1.min)*3/4)), ncol=1) tem2=solve(A, vecb0); tem3=solve(A, vecb1) a0=tem2[1]; b0=tem2[2]; q0=tem2[3]; a1=tem3[1]; b1=tem3[2]; q1=tem3[3]; out2=sa3(y, t, ps, a0, a1, b0, b1, q0, q1) out=matrix(c(out1[1,1], out2[1,1], out1[2,1], out2[2,1], out1[3,1], out2[3,1]), ncol=1) return(out); }