prev.u
. With the default method
(or.method="x:u y:u"
), $U$ is sampled so that the $X:U$
odds ratio is $a$ and the $Y:U$ odds ratio is $b$. With the
second method, $U$ is sampled according to the model $logit(U=1
| X, Y) = \alpha + \beta*Y + \gamma*X$, where $\beta=\log(b)$ and
$\gamma=\log(a)$ and $\alpha$ is determined so that the
prevalence of $U=1$ is prev.u
. This second method results in
the adjusted odds ratio for $Y:U$ given $X$ being $b$
whereas the default method forces the unconditional (marginal) $Y:U$
odds ratio to be $b$. Rosenbaum uses the default method.There is a plot
method for plotting objects created by
sensuc
. Values of $a$ are placed on the x-axis and observed
marginal odds or hazards ratios for $U$ (unadjusted ratios) appear
on the y-axis. For Cox models, the hazard ratios will not agree exactly
with $X$:event indicator odds ratios but they sometimes be made
close through judicious choice of the event
function. The
default plot uses four symbols which differentiate whether for the
$a,b$ combination the effect of $X$ adjusted for $U$ (and
for any other covariables that were in the original model fit) is
positive (usually meaning an effect ratio greater than 1) and
"significant", merely positive, not positive and non significant, or not
positive but significant. There is also an option to draw the numeric
value of the $X$ effect ratio at the $a$,$b$ combination
along with its $Z$ statistic underneath in smaller letters, and an
option to draw the effect ratio in one of four colors depending on the
significance of the $Z$ statistic.
# fit <- lrm(formula=y ~ x + other.predictors, x=TRUE, y=TRUE) #or
# fit <- cph(formula=Surv(event.time,event.indicator) ~ x + other.predictors,
# x=TRUE, y=TRUE)sensuc(fit,
or.xu=seq(1, 6, by = 0.5), or.u=or.xu,
prev.u=0.5, constrain.binary.sample=TRUE,
or.method=c("x:u y:u","u|x,y"),
event=function(y) if(is.matrix(y))y[,ncol(y)] else 1*y)
## S3 method for class 'sensuc':
plot(x, ylim=c((1+trunc(min(x$effect.u)-.01))/
ifelse(type=='numbers',2,1),
1+trunc(max(x$effect.u)-.01)),
xlab='Odds Ratio for X:U',
ylab=if(x$type=='lrm')'Odds Ratio for Y:U' else
'Hazard Ratio for Y:U',
digits=2, cex.effect=.75, cex.z=.6*cex.effect,
delta=diff(par('usr')[3:4])/40,
type=c('symbols','numbers','colors'),
pch=c(15,18,5,0), col=c(2,3,1,4), alpha=.05,
impressive.effect=function(x)x > 1,...)
lrm
or cph
with x=TRUE, y=TRUE
. The
first variable in the right hand side of the model formula must have
been the binary $X$ variable, and it may not interact with other
predictors.sensuc
or.xu
.event
is left at its default value, which
is the identify function, i.e, the orplot
type="numbers"
or type="colors"
."symbols"
(the default), "numbers"
, or "colors"
(see above)pch
TRUE
for a
positive effect. By default, a positive effect is taken to mean a
ratio exceeding one.plot
sensuc
returns an object of class "sensuc"
with the following elements: OR.xu
(vector of desired $X:U$ odds ratios or $a$ values), OOR.xu
(observed marginal $X:U$ odds ratios), OR.u
(desired $Y:U$ odds
ratios or $b$ values), effect.x
(adjusted odds or hazards ratio for
$X$ in a model adjusted for $U$ and all of the other predictors),
effect.u
(unadjusted $Y:U$ odds or hazards ratios), effect.u.adj
(adjusted $Y:U$ odds or hazards ratios), $Z$ (Z-statistics), prev.u
(input to sensuc
), cond.prev.u
(matrix with one row per $a$,$b$
combination, specifying prevalences of $U$ conditional on $Y$ and $X$
combinations), and type
("lrm"
or "cph"
).Rosenbaum P, Rubin D (1983): Assessing sensitivity to an unobserved binary covariate in an observational study with binary outcome. J Roy Statist Soc B 45:212--218.
Lee WC (2011): Bounding the bias of unmeasured factors with confounding and effect-modifying potentials. Stat in Med 30:1007-1017.
lrm
, cph
, sample
set.seed(17)
x <- sample(0:1, 500,TRUE)
y <- sample(0:1, 500,TRUE)
y[1:100] <- x[1:100] # induce an association between x and y
x2 <- rnorm(500)
f <- lrm(y ~ x + x2, x=TRUE, y=TRUE)
#Note: in absence of U odds ratio for x is exp(2nd coefficient)
g <- sensuc(f, c(1,3))
# Note: If the generated sample of U was typical, the odds ratio for
# x dropped had U been known, where U had an odds ratio
# with x of 3 and an odds ratio with y of 3
plot(g)
# Fit a Cox model and check sensitivity to an unmeasured confounder
# f <- cph(Surv(d.time,death) ~ treatment + pol(age,2)*sex, x=TRUE, y=TRUE)
# sensuc(f, event=function(y) y[,2] & y[,1] < 365.25 )
# Event = failed, with event time before 1 year
# Note: Analysis uses f$y which is a 2-column Surv object
Run the code above in your browser using DataLab