Skip to content

Commit

Permalink
version 0.5
Browse files Browse the repository at this point in the history
  • Loading branch information
@gaborcsardi
Luke J. Keele authored and gaborcsardi committed Aug 16, 2010
1 parent 00f29d4 commit fe1b588
Show file tree
Hide file tree
Showing 18 changed files with 780 additions and 400 deletions.
11 changes: 6 additions & 5 deletions DESCRIPTION
View file
@@ -1,17 +1,18 @@
Package: rbounds
Type: Package
Title: Perform Rosenbaum bounds sensitivity tests for matched data.
Version: 0.4
Date: 2009-10-21
Version: 0.5
Date: 2010-08-16
Author: Luke J. Keele
Maintainer: Luke J. Keele <keele.4@osu.edu>
Depends: R (>= 2.7.1), Matching (>= 4.6-2)
Depends: R (>= 2.7.1), Matching (>= 4.6-2), lattice
Imports: lattice
Description: Takes matched data and calculates Rosenbaum bounds for the
treatment effect. Calculates bounds for binary data,
Hodges-Lehmann point estimates, Wilcoxon signed-rank test, and
for data with multiple matched controls. Is designed to work
with the Matching package and operate on Match() objects.
License: GPL (>= 2)
Packaged: Wed Oct 21 21:05:14 2009; Luke
Packaged: Tue Aug 17 14:49:11 2010; Luke
Repository: CRAN
Date/Publication: 2009-10-23 06:41:12
Date/Publication: 2010-08-17 18:47:28
10 changes: 10 additions & 0 deletions NAMESPACE
View file
@@ -0,0 +1,10 @@
import(lattice)

export(binarysens, data.prep, hlsens, mcontrol,
iv_sens, psens, iv_pval, iv_gssearch, maim_dist)

S3method(plot, ivsens)
S3method(print, ivsens)
S3method(print, rbounds)
S3method(summary, ivsens)
S3method(summary, rbounds)
78 changes: 42 additions & 36 deletions R/binarysens.R
View file
@@ -1,45 +1,51 @@
binarysens <- function(x, y=NULL, Gamma=6, GammaInc=1)
{
if (length(x)==1){
ctrl <- x
trt <- y
}
else
{
y.c <- x$mdata$Y[x$mdata$Tr==0]
y.t <- x$mdata$Y[x$mdata$Tr==1]
table(y.t, y.c)
y.tmp1 <- table(y.t, y.c)[2]
y.tmp2 <- table(y.t, y.c)[3]
if(y.tmp1 >= y.tmp2){
trt <- y.tmp1
ctrl <- y.tmp2
} else {
trt <- y.tmp2
ctrl <- y.tmp1
}
binarysens <- function(x, y=NULL, Gamma=6, GammaInc=1) {
if (length(x)==1) {
ctrl <- x
trt <- y
}
else {
y.c <- x$mdata$Y[x$mdata$Tr==0]
y.t <- x$mdata$Y[x$mdata$Tr==1]
table(y.t, y.c)
y.tmp1 <- table(y.t, y.c)[2]
y.tmp2 <- table(y.t, y.c)[3]

if(y.tmp1 >= y.tmp2) {
trt <- y.tmp1
ctrl <- y.tmp2
}
else {
trt <- y.tmp2
ctrl <- y.tmp1
}
}
gamma <- seq(1, Gamma, by=GammaInc)
mx <- ctrl + trt
up <- c()
lo <- c()
series <- seq(trt, mx, by=1)
n.it <- length(gamma)
for(i in 1:n.it)
{
p.plus <- gamma[i]/(1 + gamma[i])
p.minus <- 1/(1 + gamma[i])
up.tmp <- sum(dbinom(series, mx, p=p.plus))
lo.tmp <- sum(dbinom(series, mx, p=p.minus))
up <- c(up, up.tmp)
lo <- c(lo, lo.tmp)
}
for(i in 1:n.it) {
p.plus <- gamma[i]/(1 + gamma[i])
p.minus <- 1/(1 + gamma[i])
up.tmp <- sum(dbinom(series, mx, p=p.plus))
lo.tmp <- sum(dbinom(series, mx, p=p.minus))
up <- c(up, up.tmp)
lo <- c(lo, lo.tmp)
}

pval <- lo[1]
bounds <- data.frame(gamma, round(lo, 5), round(up, 5))
colnames(bounds) <- c("Gamma", "Lower bound", "Upper bound")

msg <- "Rosenbaum Sensitivity Test \n"
note <- "Note: Gamma is Odds of Differential Assignment To
Treatment Due to Unobserved Factors \n"

Obj <- list(Gamma = Gamma, GammaInc = GammaInc, pval = pval,
msg = msg, bounds = bounds, note = note)
class(Obj) <- c("rbounds", class(Obj))

out <- cbind(gamma, round(lo, 5), round(up, 5))
colnames(out) <- c("Gamma", "Lower Bound P-Value", "Upper Bound P-Value")
cat("Rosenbaum Sensitivity Test\n")
print(out, scientific=FALSE)
cat("\n")
cat("Note: Gamma is Log Odds of Differential Assignment To Treatment Due To Unobserved Factors \n")
} #end of binary.sens function
Obj
}

93 changes: 47 additions & 46 deletions R/data.prep.R
View file
@@ -1,51 +1,52 @@

data.prep <- function(obj, Y=NULL, group.size=3)
{
if(is.null(Y))
{
data.prep <- function(obj, Y=NULL, group.size=3)
{
if(is.null(Y)) {
ctrl <- obj$mdata$Y[obj$mdata$Tr==0]
trt <- obj$mdata$Y[obj$mdata$Tr==1]
} else {
}
else {
ctrl <- Y[obj$index.control]
trt <- Y[obj$index.treated]
}
if(group.size==3){
n.i <- length(trt)
idx <- seq(2,n.i,2)
trt <- trt[-idx]
n.t <- length(trt)
grp.dta <- c()
j <- 1
for(i in 1:n.t)
{
tmp <- c(trt[i], ctrl[j:(j+1)])
grp.dta <- c(grp.dta, tmp)
j <- j + 2
rm(tmp)
}
} else {
n.i <- length(trt)
idx <- seq(2,n.i,2)
trt <- trt[-idx]
trt <- trt[-idx]
n.t <- length(trt)
grp.dta <- c()
j <- 1
for(i in 1:n.t)
{
tmp <- c(trt[i], ctrl[j:(j+2)])
grp.dta <- c(grp.dta, tmp)
j <- j + 3
rm(tmp)
}
}
id <- rep(1:n.t, each=group.size)
zeros <- rep(0,group.size-1)
trt.ind <- c(1,zeros)
treat <- rep(trt.ind, times=n.t)
mctrl <- list()
mctrl$Y <- grp.dta
mctrl$id <- id
mctrl$treat <- treat
mctrl
}

if(group.size==3){
n.i <- length(trt)
idx <- seq(2,n.i,2)
trt <- trt[-idx]
n.t <- length(trt)
grp.dta <- c()
j <- 1
for(i in 1:n.t) {
tmp <- c(trt[i], ctrl[j:(j+1)])
grp.dta <- c(grp.dta, tmp)
j <- j + 2
rm(tmp)
}
}
else {
n.i <- length(trt)
idx <- seq(2,n.i,2)
trt <- trt[-idx]
trt <- trt[-idx]
n.t <- length(trt)
grp.dta <- c()
j <- 1

for(i in 1:n.t) {
tmp <- c(trt[i], ctrl[j:(j+2)])
grp.dta <- c(grp.dta, tmp)
j <- j + 3
rm(tmp)
}
}

id <- rep(1:n.t, each=group.size)
zeros <- rep(0,group.size-1)
trt.ind <- c(1,zeros)
treat <- rep(trt.ind, times=n.t)
mctrl <- list()
mctrl$Y <- grp.dta
mctrl$id <- id
mctrl$treat <- treat
mctrl
}
63 changes: 33 additions & 30 deletions R/hlsens.R
View file
@@ -1,18 +1,16 @@

hlsens <- function(x, y=NULL, pr=.1, Gamma=6, GammaInc=1)
{
hlsens <- function(x, y = NULL, pr = 0.1, Gamma = 6, GammaInc = 1) {

if (is.numeric(x)){
trt <- x
ctrl <- y
if (is.numeric(x)) {
trt <- x
ctrl <- y
}
else {
ctrl <- x$mdata$Y[x$mdata$Tr==0]
trt <- x$mdata$Y[x$mdata$Tr==1]
}
gamma <- seq(1, Gamma, by=GammaInc)
k <- length(gamma)
trt <- x$mdata$Y[x$mdata$Tr==1]
ctrl <- x$mdata$Y[x$mdata$Tr==0]
}

gamma <- seq(1, Gamma, by=GammaInc)
k <- length(gamma)

ttau <- function(x) {
tau <- x
Expand All @@ -23,21 +21,22 @@ hlsens <- function(x, y=NULL, pr=.1, Gamma=6, GammaInc=1)
sum(psi * ranks)
}

tau.up <- tau.l <- wilcox.test(trt, ctrl, paired=TRUE, conf.int=TRUE, exact=FALSE)$estimate
tau.up <- tau.l <- wilcox.test(trt, ctrl, paired=TRUE,
conf.int=TRUE, exact=FALSE)$estimate
#base <- c(1,tau.up,tau.l)
eps <- 1.e-8
c.int <- matrix(0,k,2)
eps <- 1.0e-8
c.int <- matrix(0, k, 2)
s <- length(trt)

for(i in 1:k){
for(i in 1:k) {
p.minus = 1/(1+gamma[i])
p.plus = gamma[i]/(gamma[i]+1)
t.min <- p.minus*(s*(s+1)/2)
t.max <- p.plus*(s*(s+1)/2)
lb <- t.min
ub <- t.max

while(abs(ub - lb) > eps){
while(abs(ub - lb) > eps) {
if (lb < ub){
tau.old <- tau.up
tau.up <- tau.old + pr
Expand All @@ -46,29 +45,33 @@ hlsens <- function(x, y=NULL, pr=.1, Gamma=6, GammaInc=1)
else break
}
c.int[i,2] <- tau.up

ub <- t.max
# lb <- ttau(tau.l)
# lb <- ttau(tau.l)
lb <- t.min

while (abs(ub - lb) > eps){
if (lb <= ub){
while (abs(ub - lb) > eps) {
if (lb <= ub) {
tau.old <- tau.l
tau.l <- tau.old - pr
tau.l <- tau.old - pr
lb <- ttau(tau.l)
}
else break
}
c.int[i,1] <- tau.l
}

pval <- c.int[1,1]
bounds <- data.frame(gamma, signif(c.int, digits=5))
colnames(bounds) <- c("Gamma", "Lower bound", "Upper bound")

out <- cbind(gamma, signif(c.int, digits=5))
#out <- rbind(base, out)
colnames(out) <- c("Gamma", "L. Bound HL Est.", "U. Bound HL Est.")
cat("Rosenbaum Sensitivity Test for Hodges-Lehmann Point Estimate \n")
print(out, scientific=FALSE)
cat("\n")
cat("Note: Gamma is Log Odds of Differential Assignment To Treatment Due To Unobserved Factors \n")
msg <- "Rosenbaum Sensitivity Test for Hodges-Lehmann Point Estimate \n"
note <- "Note: Gamma is Odds of Differential Assignment To
Treatment Due to Unobserved Factors \n"

Obj <- list(Gamma = Gamma, GammaInc = GammaInc, pval = pval,
msg = msg, bounds = bounds, note = note)
class(Obj) <- c("rbounds", class(Obj))

Obj
}

0 comments on commit fe1b588

Please sign in to comment.