This is admitted a very brief lecture on sensitivity analysis under the Rosenbaum’s sensitivity analysis framework. The purpose of this 40 minute lecture is to illustrate the key ideas of sensitivity analysis under a very particular setup (i.e. matched pairs in finite samples) and some computer code to do this. For more details, see Rosenbaum’s 2010 book on the Design of Observational Studies.
Suppose the study units are paired where each pair consists of one treated unit and a control unit and we are interested in testing the sharp null hypothesis of no treatment effect. I assume that the population is fixed and the only thing that is random is the treatment assignment within each pair. This is knowin # Paired Study of Smoking and Lung Cancer
Smoking is treatment (\(A=1\) smoking; \(A=0\) non-smoking). Outcome is lung cancer (\(Y=1\) is lung cancer; \(Y=0\) is no cancer). There were 36,975 pairs. 122 pairs were discordant where one of the pair had lung cancer; the rest of 36,853 pairs had same outcomes between the pairs.
Suppose we want to test the sharp null \(H_0: Y_{si}(1) = Y_{si}(0)\) for all \(i\) and we conduct randomization inference (i.e. potential outcomes fixed). We want to test this null using the test statistic \[ T = \sum_{s =1}^I \delta_{s} \] where \(I\) is total number of pairs and \(\delta_{s}\) is the treated-minus-control difference in outcome for each pair \(s\). Let \(d\) be the number of discordant pairs. In our observed data, \(d = 122\).
The null distribution of \(T\) was
shown in lecture to be \[
T = \sum_{s: \text{discordant pairs}} Q_i, \quad{} Q_i = \begin{cases} 1
& p_s \\ -1 & 1- p_s \end{cases}
\] where \(p_s\) is the
probability that first unit in the pair gets treated and second unit in
the pair gets control and \(1 -p_s\) is
the probability that first unit in the pair gets control and second unit
in the pair gets treated. We can transform this variable to look like
\[
\frac{T + d}{2} = \sum_{s: \text{discordant pairs}} \frac{Q_i + 1}{2}
\] where within each summand, it is a Bernounilli distribution
with probability \(p_s\) of being \(1\) and \(1-p_s\) of being \(0\). If each pair was equally likely be a
smoker, then \(p_s = 0.5\) for all
\(s\), then \((T+d)/2 \sim {\rm Binom}(d,1/2)\) and we
can compute the p-value as
\[
P\left( \frac{T + d}{2} \geq \frac{t_{\rm obs} +d}{2} \right) =
P\left({\rm Binom}(122,0.5) \geq \frac{110 - 12}{2} + \frac{122}{2}
\right) = P\left({\rm Binom(122, 0.5) } \geq 110 \right) \approx 0
\] In short, the p-value is close to zero and it’s unlikely that
\(H_0\) is true assuming \(p_s = 1/2\)
1 - pbinom(109,122,0.5) #It's 109, not 110 because pbinom counts *upto* 110## [1] 0Now, suppose we think one unit in the pair is more likely to get treatment than control due to unobserved differences between the pair. In lecture, we bounded this probability as \[ \frac{1}{1 + \Gamma} \leq p_s \leq \frac{\Gamma}{1+\Gamma} \] We shown in lecture that this corresponds to the following worst-case distributions \[ P\left({\rm Binom}\left(122,\frac{1}{1+\Gamma} \right) \geq 110 \right) \leq P\left(\frac{T+d}{2}\geq 110 \right) \leq P\left({\rm Binom}\left(122,\frac{\Gamma}{1+\Gamma}\right) \geq 110 \right) \] Here is a list of the upper and lower bounds as a function of \(\Gamma\)
GammaSeq = c(1,2,3,4,5,5.5,6,7)
boundMatrix = matrix(0,length(GammaSeq),2)
for(i in 1:length(GammaSeq)) {
boundMatrix[i,2] =1 - pbinom(109,122,GammaSeq[i]/(1 + GammaSeq[i]))
boundMatrix[i,1] =1 - pbinom(109,122,1/(1 + GammaSeq[i]))
}
out = cbind(GammaSeq,boundMatrix)
colnames(out) = c("Gamma","Lower Bound P-value","Upper Bound P-value")
round(out,5)## Gamma Lower Bound P-value Upper Bound P-value
## [1,] 1.0 0 0.00000
## [2,] 2.0 0 0.00000
## [3,] 3.0 0 0.00002
## [4,] 4.0 0 0.00196
## [5,] 5.0 0 0.02317
## [6,] 5.5 0 0.05197
## [7,] 6.0 0 0.09693
## [8,] 7.0 0 0.23031Even if \(p_s \neq 1/2\) and one pair is more likely to smoke because of an unobserved variable, this unobserved variable must have a \(\Gamma> 5.5\)-fold difference away from 50/50 randomization to not reject the null hypothesis at \(\alpha= 0.05\) that smoking does not cause lung cancer. In other words, an unmeasured variable must change \(p_s\) outside of the range \(0.15 \leq p_s \leq 0.85\) for us to retain the null hypothesis.
We can re-interpret \(\Gamma\) for pairs and binary outcomes by using the following equation \[ \Gamma = \frac{ \beta_{\rm UA} \beta_{\rm UY} + 1}{\beta_{\rm UA} + \beta_{\rm UY}} \] where \(\beta_{\rm UA}\) measures an unmeasured confounder \(U\)’s odds of increasing smoking \(A=1\) and \(\beta_{\rm UY}\) measures an unmeasured confounder \(U\)’s odds of increasing lung cancer \(Y=1\). Typically,we set \(\beta_{\rm UA} =\beta_{\rm UY} = E\) and find the ``evidence value’’ \(E\) given \(\Gamma\) where the p-value range goes just below \(\alpha = 0.05\), i.e. \(E = \Gamma \pm \sqrt{\Gamma^2 -1}\).
For example, if \(\Gamma = 5\), which was the \(\Gamma\) value that still rejected the null in favor of smoking causing lung cancer, this means even if an unmeasured confounder such as a genetic variant increased the odds of smoking 9.9 \(= 5 + \sqrt{5^2 - 1}\) times and the odds of lung cancer 9.9 times, we would still have evidence that smoking causes lung cancer.
In short, even with such an unmeasured confounder, we would still be able to conclude that smoking caused lung cancer. This observational study’s conclusion is very insensitive to unmeasured confounders.
See this paper for details on software. Rosenbaum even has a website/shiny app to do sensitivity analysis.
We’ll analyze the effect of veteran status on smoking.
## iduser pop57 edfa57q edmo57q sesp57 hsdm57 hsrankq tchevl gwiiq_bm relfml
## 1 1007 6 7 12 4 3 41 0 86 44
## 2 1008 6 7 10 3 3 17 0 86 40
## 3 1010 6 13 14 3 3 35 0 88 1
## 4 1013 6 10 10 3 3 87 0 123 96
## 5 1016 6 -3 16 3 3 84 0 107 39
## 6 1018 6 -1 12 3 3 29 0 92 96
## milty packyears outcome
## 1 0 0 0
## 2 1 70 70
## 3 1 25 25
## 4 0 40 40
## 5 0 0 0
## 6 0 16 16Here is the result for pair matching.
# PairMatching
# I realize this code is highly inefficient. Improvements are welcome!
matchvec=pairmatch(distmat.propensity,controls=nocontrols.per.match,data=datatemp)
datatemp$matchvec=matchvec
## Create a matrix saying which control units each treated unit is matched to
## Create vectors of the subject indices of the treatment units ordered by
## their matched set and corresponding control unit
treated.subject.index=rep(0,sum(treated==1))
matched.control.subject.index.mat=matrix(rep(0,nocontrols.per.match*length(treated.subject.index)),ncol=nocontrols.per.match)
matchedset.index=substr(matchvec,start=3,stop=10)
matchedset.index.numeric=as.numeric(matchedset.index)
for(i in 1:length(treated.subject.index)){
matched.set.temp=which(matchedset.index.numeric==i)
treated.temp.index=which(datatemp$treated[matched.set.temp]==1)
treated.subject.index[i]=matched.set.temp[treated.temp.index]
matched.control.subject.index.mat[i,]=matched.set.temp[-treated.temp.index]
}
matched.control.subject.index=matched.control.subject.index.mat
Yvec = propscore.model$data$outcome
Ymat = matrix(rep(0,(1+nocontrols.per.match)*length(treated.subject.index)),ncol=1+nocontrols.per.match)
for(i in 1:length(treated.subject.index)) {
Ymat[i,1] = Yvec[treated.subject.index[i]]
Ymat[i,2:(1+nocontrols.per.match)] = Yvec[matched.control.subject.index.mat[i,]]
}
colnames(Ymat) = c("Treated Outcome",paste("Control",1:nocontrols.per.match,"Outcome"))
head(Ymat)## Treated Outcome Control 1 Outcome
## [1,] 70 25
## [2,] 25 35
## [3,] 20 74
## [4,] 6 0
## [5,] 15 0
## [6,] 70 0library(sensitivitymw)
# senmw does matching with multiple controls (same number of controls)
GammaSeq = seq(1,2,0.1)
upperBound = matrix(0,length(GammaSeq),3)
for(i in 1:length(GammaSeq)) {
upperBound[i,1] =senmw(Ymat, gamma = GammaSeq[i], method = "t")$pval
upperBound[i,2] =senmw(Ymat, gamma = GammaSeq[i], method = "p")$pval
upperBound[i,3] = senmw(Ymat, gamma = GammaSeq[i], method = "w")$pval
}
out = cbind(GammaSeq,upperBound)
colnames(out) = c("Gamma","t-test","trimmed mean test","weighted trimmed mean test")
round(out,3)## Gamma t-test trimmed mean test weighted trimmed mean test
## [1,] 1.0 0.000 0.000 0.000
## [2,] 1.1 0.000 0.000 0.000
## [3,] 1.2 0.011 0.000 0.000
## [4,] 1.3 0.091 0.004 0.007
## [5,] 1.4 0.321 0.040 0.050
## [6,] 1.5 0.636 0.184 0.192
## [7,] 1.6 0.866 0.456 0.441
## [8,] 1.7 0.966 0.736 0.702
## [9,] 1.8 0.994 0.909 0.879
## [10,] 1.9 0.999 0.977 0.962
## [11,] 2.0 1.000 0.996 0.991Notice that the weighted trimmed mean is least sensitive to unobserved confounder of magnitude \(\Gamma\). The R package also has the ability to compute range of point estimates as well as range of confidence intervals.
senmwCI(Ymat, gamma = 1.35, method = "w", one.sided = TRUE)## $PointEstimate
## minimum maximum
## 2.50 10.41
##
## $Confidence.Interval
## minimum maximum
## 0.5 InfWe can also amplify the sensitivity analysis like discussed in lecture using the sensitivitymv package. The amplification also works for continuous outcomes, but the intreprtation is based purely in terms of the odds of observing a positive outcome.
library(sensitivitymv)##
## Attaching package: 'sensitivitymv'## The following objects are masked from 'package:sensitivitymw':
##
## mscorev, multrnks, newurksamplify(1.35, c(4 : 7)) #Here 4:7 is the \beta_{UA} amplitude## 4 5 6 7
## 1.660377 1.575342 1.526882 1.495575uniroot(function(x){amplify(1.35,x) - x},c(1.35+0.01,10))$root # E-value computation, with 0.01 added for uniroot numerical issues## [1] 2.256935Thus, to overturn the conclusion that veteran status causes an increase in number of cigarettes smoked, you would need an unmeasured confounder whose odds of increasing the outcome and treatment is about \(2.26\) (from E-value). Alternatively, you need a \(\Gamma = 1.35\), which roughly corresponds to the probability bound \(0.426 \leq p_s \leq 0.574\) for one of the units in the pair taking the treatment. Broadly speaking, the conclusion is sensitive to unmeasured biases.