Nothing
tests/testthat/test-effect-mod.R
In EValue: Sensitivity Analyses for Unmeasured Confounding and Other Biases in Observational Studies and Meta-Analyses
# no need to source helper_testthat.R explicitly to pass CRAN checks
# because helper files that start with "helper" automatically
# get sourced first: https://testthat.r-lib.org/reference/test_dir.html
# # for local testing:
# library(testthat)
# library(devtools)
# library(dplyr)
# library(ICC)
# library(msm)
# library(MetaUtility)
# library(here())
#
# setwd(here())
# setwd("tests")
# source("helper_testthat.R")
# load_all()
# # end part for local testing
# enter example datasets (Letenneur)
# Y: dementia
# X: low education
# n: sample size
# data for women
nw_1 = 2988
nw_0 = 364
dw = data.frame( Y = c(1, 1, 0, 0),
X = c(1, 0, 1, 0),
n = c( 158, 6, nw_1-158, nw_0-6 ) )
# data for men
nm_1 = 1790
nm_0 = 605
dm = data.frame( Y = c(1, 1, 0, 0),
X = c(1, 0, 1, 0),
n = c( 64, 17, nm_1-64, nm_0-17 ) )
# P(Y = 1 | X = 1) and P(Y = 1 | X = 0) for women and for men
( pw_1 = dw$n[ dw$X == 1 & dw$Y == 1 ] / sum(dw$n[ dw$X == 1 ]) )
( pw_0 = dw$n[ dw$X == 0 & dw$Y == 1 ] / sum(dw$n[ dw$X == 0 ]) )
( pm_1 = dm$n[ dm$X == 1 & dm$Y == 1 ] / sum(dm$n[ dm$X == 1 ]) )
( pm_0 = dm$n[ dm$X == 0 & dm$Y == 1 ] / sum(dm$n[ dm$X == 0 ]) )
# prevalence of low education among women and among men
fw = nw_1 / (nw_1 + nw_0)
fm = nm_1 / (nm_1 + nm_0)
# confounded estimates
RDw = pw_1 - pw_0
RDm = pm_1 - pm_0
test_that("RDt_var, symmetry when reversing sign of RD", {
# test for symmetry
v1 = RDt_var( f = .25,
p1 = 0.3,
p0 = 0.1,
n1 = 50,
n0 = 100,
.maxB = 2 )
# here I've changed argument order to reverse sign of RD
v2 = RDt_var( f = .25,
p0 = 0.3,
p1 = 0.1,
n0 = 50,
n1 = 100,
.maxB = 2 )
expect_equal(v1, v2)
})
test_that( "RDt_bound, symmetry when shifting strata in opposite directions", {
x1 = RDt_bound( p1_1 = 0.6,
p1_0 = 0.4,
n1_1 = 100,
n1_0 = 10,
f1 = 0.25,
maxB_1 = 2,
biasDir_1 = "positive",
p0_1 = 0.4,
p0_0 = 0.6,
n0_1 = 10,
n0_0 = 100,
f0 = .25,
biasDir_0 = "negative",
alpha = 0.05 )
expect_equal( x1$RD[1], -x1$RD[2] )
expect_equal( x1$se[1], x1$se[2] )
expect_equal( x1$lo[1], -x1$hi[2] )
expect_equal( x1$lo[2], -x1$hi[1] )
expect_equal( x1$pval[1], x1$pval[2] )
# sanity check: equality when both strata are the same and are positively biased
x1 = RDt_bound( p1_1 = 0.6,
p1_0 = 0.4,
n1_1 = 100,
n1_0 = 10,
f1 = 0.25,
maxB_1 = 1.6,
biasDir_1 = "positive",
p0_1 = 0.6,
p0_0 = 0.4,
n0_1 = 100,
n0_0 = 10,
f0 = .25,
maxB_0 = 1.6,
biasDir_0 = "positive",
alpha = 0.05 )
expect_equal( x1$RD[1], x1$RD[2] )
expect_equal( x1$se[1], x1$se[2] )
expect_equal( x1$lo[1], x1$lo[2] )
expect_equal( x1$hi[1], x1$hi[2] )
expect_equal( x1$pval[1], x1$pval[2] )
# sanity check: equality when both strata are the same and are NEGATIVELY biased
x1 = RDt_bound( p1_1 = 0.6,
p1_0 = 0.4,
n1_1 = 100,
n1_0 = 10,
f1 = 0.25,
maxB_1 = 1.6,
biasDir_1 = "negative",
p0_1 = 0.6,
p0_0 = 0.4,
n0_1 = 100,
n0_0 = 10,
f0 = .25,
maxB_0 = 1.6,
biasDir_0 = "negative",
alpha = 0.05 )
expect_equal( x1$RD[1], x1$RD[2] )
expect_equal( x1$se[1], x1$se[2] )
expect_equal( x1$lo[1], x1$lo[2] )
expect_equal( x1$hi[1], x1$hi[2] )
expect_equal( x1$pval[1], x1$pval[2] )
})
test_that( "Bound from RDt_bound should agree with closed form in paper", {
B = 2.3
term1 = (pw_1 - pw_0 * B) * ( fw + (1 - fw) / B )
term2 = (pm_1 * B - pm_0) * ( fm + (1 - fm) / B )
mine = term1 - term2
x = RDt_bound( p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
biasDir_1 = "positive",
maxB_1 = B,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
biasDir_0 = "negative",
maxB_0 = B )
expect_equal( x$RD[ x$stratum == "effectMod" ], mine, tol = 0.0001 )
} )
test_that( "E-value for 1 stratum from IC_evalue should match R package", {
( evalueEst = IC_evalue_inner( stratum = "1",
varName = "RD",
true = 0,
unidirBias = FALSE,
p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
alpha = 0.05 )$evalues )
( evalueCI = IC_evalue_inner( stratum = "1",
varName = "lo",
true = 0,
unidirBias = FALSE,
p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
alpha = 0.05 )$evalues )
# now try against package:
evalueOld = EValue::evalues.RD( n11 = dw$n[ dw$X == 1 & dw$Y == 1 ],
n10 = dw$n[ dw$X == 1 & dw$Y == 0 ],
n01 = dw$n[ dw$X == 0 & dw$Y == 1 ],
n00 = dw$n[ dw$X == 0 & dw$Y == 0 ],
true = 0,
alpha = 0.05)
# they agree! :D
# woohoo!!!!!
expect_equal( evalueOld$est.Evalue, evalueEst$evalue, tol = 0.001 )
expect_equal( evalueOld$lower.Evalue, evalueCI$evalue, tol = 0.001 )
})
test_that( "Bound from RDt_bound should be symmetric after flipping signs of both RDs", {
# shift each stratum by different amount
B1 = 1.5
B0 = 1.2
# both RDs > 0
( x = RDt_bound( p1_1 = .8,
p1_0 = .6,
n1_1 = 20,
n1_0 = 200,
f1 = 0.5,
biasDir_1 = "positive",
maxB_1 = B1,
p0_1 = .8,
p0_0 = .75,
n0_1 = 30,
n0_0 = 40,
f0 = 0.5,
biasDir_0 = "positive",
maxB_0 = B0 ) )
# both RDs < 0
# strata labels reversed so that the interaction contrast itself remains positive
( x2 = RDt_bound( p0_1 = .6,
p0_0 = .8,
n0_1 = 200,
n0_0 = 20,
f0 = 0.5,
biasDir_0 = "negative",
maxB_0 = B1,
p1_1 = .75,
p1_0 = .8,
n1_1 = 40,
n1_0 = 30,
f1 = 0.5,
biasDir_1 = "negative",
maxB_1 = B0 ) )
# interaction contrast should be exactly the same
expect_equal( x[3, 3:6], x2[3, 3:6] )
# strata are reversed and signs are reversed
# inference stays same
# but CIs flip
expect_equal( as.numeric( x[1, c("se","pval")] ), as.numeric( x2[2, c("se", "pval")] ) )
expect_equal( as.numeric( x[2, c("se","pval")] ), as.numeric( x2[1, c("se", "pval")] ) )
expect_equal( x[1, "RD"], -x2[2, "RD"] )
expect_equal( x[2, "RD"], -x2[1, "RD"] )
})
test_that( "E-value should be the same when flipping strata signs", {
x = evalues.IC( stat = "est",
true = 0.1,
unidirBias = FALSE,
p1_1 = .5,
p1_0 = .3,
n1_1 = 100,
n1_0 = 20,
f1 = .6,
p0_1 = .4,
p0_0 = .35,
n0_1 = 40,
n0_0 = 60,
f0 = .2,
alpha = 0.05 )$evalues
x2 = evalues.IC( stat = "est",
true = 0.1,
unidirBias = FALSE,
p0_1 = .3,
p0_0 = .5,
n0_1 = 20,
n0_0 = 100,
f0 = .6,
p1_1 = .35,
p1_0 = .4,
n1_1 = 60,
n1_0 = 40,
f1 = .2,
alpha = 0.05 )$evalues
expect_equal( x$evalue, x2$evalue, tol = 0.001 )
expect_equal( x$biasFactor, x2$biasFactor, tol = 0.001 )
} )
test_that("evalues.IC should warn if E-value is 1", {
# E-value for estimate is 1
expect_message( evalues.IC( stat = "est",
true = .5,
unidirBias = FALSE,
p1_1 = .5,
p1_0 = .4,
n1_1 = 100,
n1_0 = 100,
f1 = .5,
p0_1 = .5,
p0_0 = .42,
n0_1 = 100,
n0_0 = 100,
f0 = .5,
alpha = 0.05 ) )
# E-value for CI is 1
# first check the CI limit
RDt_bound( p1_1 = .5,
p1_0 = .4,
n1_1 = 100,
n1_0 = 100,
f1 = .5,
p0_1 = .5,
p0_0 = .42,
n0_1 = 100,
n0_0 = 100,
f0 = .5,
# no bias
maxB_1 = 1,
maxB_0 = 1,
biasDir_1 = "positive",
biasDir_0 = "positive" )
expect_message( evalues.IC( stat = "CI",
true = .5,
unidirBias = FALSE,
p1_1 = .5,
p1_0 = .4,
n1_1 = 100,
n1_0 = 100,
f1 = .5,
p0_1 = .5,
p0_0 = .42,
n0_1 = 100,
n0_0 = 100,
f0 = .5,
alpha = 0.05 ) )
# also when trying both candidates (non-unidir, unknown)
expect_message( evalues.IC( stat = "CI",
true = 0.0,
unidirBias = TRUE,
unidirBiasDirection = "unknown",
p1_1 = .9,
p1_0 = .8,
n1_1 = 100,
n1_0 = 100,
f1 = .9,
p0_1 = .5,
p0_0 = .42,
n0_1 = 100,
n0_0 = 100,
f0 = .1,
alpha = 0.05 ) )
})
test_that("evalues.IC should reject bad input", {
# if unidirBias is TRUE, must provide unidirBiasDirection
expect_error( evalues.IC( stat = "est",
true = 0,
unidirBias = TRUE,
p1_1 = .4,
p1_0 = .4,
n1_1 = 100,
n1_0 = 100,
f1 = .5,
p0_1 = .5,
p0_0 = .42,
n0_1 = 100,
n0_0 = 100,
f0 = .5,
alpha = 0.05 ) )
# specified unidirBias = FALSE, so the argument unidirBiasDirection will be ignored
expect_warning( evalues.IC( stat = "est",
true = 0,
unidirBias = FALSE,
unidirBiasDirection = "positive",
p1_1 = .5,
p1_0 = .4,
n1_1 = 100,
n1_0 = 100,
f1 = .5,
p0_1 = .5,
p0_0 = .42,
n0_1 = 100,
n0_0 = 100,
f0 = .5,
alpha = 0.05 ) )
expect_warning( evalues.IC( stat = "est",
true = 0,
unidirBias = FALSE,
unidirBiasDirection = "negative",
p1_1 = .5,
p1_0 = .4,
n1_1 = 100,
n1_0 = 100,
f1 = .5,
p0_1 = .5,
p0_0 = .42,
n0_1 = 100,
n0_0 = 100,
f0 = .5,
alpha = 0.05 ) )
} )
test_that( "E-value from evalues.IC should be the solution to RDt_bound and should agree with theory in paper", {
Eadd.est = evalues.IC( stat = "est",
true = 0,
unidirBias = FALSE,
p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
alpha = 0.05 )$evalues
# pass the resulting bias factor (i.e., equivalent to alleged E-value) to each stratum
x = RDt_bound( p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
biasDir_1 = "positive",
# from above test
maxB_1 = Eadd.est$biasFactor,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
biasDir_0 = "negative",
maxB_0 = Eadd.est$biasFactor )
expect_equal( x$RD[ x$stratum == "effectMod" ], 0, tol = 0.0001 )
# should agree with closed form
gamma = pw_0 * (1 - fw) - pw_1 * fw + pm_1 * (1 - fm) - pm_0 * fm
term1 = 1 / ( 2 * (fw * pw_0 + fm * pm_1) )
term2 = 4 * (pw_0 * fw + pm_1 * fm) * ( pw_1 * (1 - fw) + pm_0 * (1-fm) )
# bias factor
( my.Badd.est = term1 * ( sqrt( gamma^2 + term2 ) - gamma ) )
expect_equal( Eadd.est$biasFactor, my.Badd.est, tol = 0.001 )
# E-value
expect_equal( Eadd.est$evalue, g(my.Badd.est), tol = 0.001 )
})
test_that("evalues.IC solution for one stratum should agree with existing fn, evalues.RD", {
( resNonMono = evalues.IC( stat = "est",
true = 0,
unidirBias = FALSE,
p1_1 = .6,
p1_0 = .4,
n1_1 = 100,
n1_0 = 100,
f1 = .5,
p0_1 = 0.3,
p0_0 = 0.20,
n0_1 = 100,
n0_0 = 100,
f0 = .5,
alpha = 0.05 )$evalues )
# RDt_bound
RDs2 = RDt_bound( p1_1 = .6,
p1_0 = .4,
n1_1 = 100,
n1_0 = 100,
f1 = .5,
biasDir_1 = "positive",
maxB_1 = resNonMono$biasFactor,
p0_1 = 0.3,
p0_0 = 0.20,
n0_1 = 100,
n0_0 = 100,
f0 = .5,
biasDir_0 = "negative",
maxB_0 = resNonMono$biasFactor )
# stratum 1 only
evalueOld2 = evalues.RD( n11 = 100 * (0.6),
n10 = 100 * (1-0.6),
n01 = 100 * 0.4,
n00 = 100 * (1-0.4),
true = RDs2$RD[1],
alpha = 0.05)
expect_equal( evalueOld2$est.Evalue, resNonMono$evalue, tol = 0.001 )
# # stratum M only
# evalueOld2 = EValue::evalues.RD( n11 = 100 * (0.3),
# n10 = 100 * (1-0.3),
# n01 = 100 * 0.2,
# n00 = 100 * (1-0.2),
# true = RDs2$RD[2],
# alpha = 0.05)
#
#
# # they agree! :D
# # woohoo!!!!!
# expect_equal( evalueOld2$est.Evalue, resNonMono$evalue, tol = 0.001 )
})
test_that("Evalue candidate (negative unidir bias) should successfully move RDm up to RDw", {
# superassign because we'll use for subsequent tests
( Eadd.est.mono <<- evalues.IC( stat = "est",
true = 0,
unidirBias = TRUE,
unidirBiasDirection = "unknown",
p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
alpha = 0.05 ) )
x = RDt_bound( p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
biasDir_1 = "negative",
maxB_1 = 1,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
biasDir_0 = "negative",
maxB_0 = Eadd.est.mono$candidates$biasFactor[ Eadd.est.mono$candidates$biasDir == "negative" ] )
expect_equal( x$RD[ x$stratum == "0" ], RDw, tol = 0.0001 )
# and likewise for CI limit
( Eadd.CI.mono <<- evalues.IC( stat = "CI",
true = 0,
unidirBias = TRUE,
unidirBiasDirection = "unknown",
p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
alpha = 0.05 ) )
x = RDt_bound( p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
biasDir_1 = "negative",
maxB_1 = 1,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
biasDir_0 = "negative",
maxB_0 = Eadd.CI.mono$candidates$biasFactor[ Eadd.CI.mono$candidates$biasDir == "negative" ] )
expect_equal( x$lo[ x$stratum == "effectMod" ], 0, tol = 0.0001 )
} )
test_that( "Evalue candidate (positive unidir bias) should successfully move RDw down to RDm", {
( x = RDt_bound( p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
biasDir_1 = "positive",
maxB_1 = Eadd.est.mono$candidates$biasFactor[ Eadd.est.mono$candidates$biasDir == "positive" ],
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
biasDir_0 = "positive",
maxB_0 = 1 ) )
expect_equal( x$RD[ x$stratum == "1" ], RDm, tol = 0.0001 )
# and likewise for CI limit
( x = RDt_bound( p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
biasDir_1 = "positive",
maxB_1 = Eadd.CI.mono$candidates$biasFactor[ Eadd.CI.mono$candidates$biasDir == "positive" ] ,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
biasDir_0 = "positive",
maxB_0 = 1 ) )
expect_equal( x$lo[ x$stratum == "effectMod" ], 0, tol = 0.0001 )
} )
test_that( "E-values from IC_evalue (grid search) should match closed form in paper", {
B = 4
true = ( pm_1 * B - pm_0 ) * ( fm + (1-fm) / B )
# suggestively name terms as in quadratic formula
termA = fm * pm_1
termB = pm_1 * ( 1 - fm ) - fm * pm_0 - true
termC = -pm_0 * (1 - fm)
# this is the polynomial that is to be solved for E-value
( mine = termA * B^2 + termB * B + termC )
expect_equal( mine, 0, tol = 0.0001 )
# check E-value for negatively biased stratum in paper
# this is the one that arises from reversing roles and signs in the existing E-value on
# Ding Appendix, pg 18
#bm
### check E-value for positive bias (shift only stratum W)
# check against Ding Appendix, pg 18 (Prop A.11)
targetB = Eadd.est.mono$candidates$biasFactor[ Eadd.est.mono$candidates$biasDir == "positive" ]
lambda = pw_0 * (1 - fw) - pw_1 * fw
term1 = 1 / ( 2 * (fw * pw_0) )
term2 = 4 * pw_1 * pw_0 * fw * (1 - fw)
true = RDm
term3 = true + lambda
# bias factor
( myB = term1 * ( sqrt( term3^2 + term2 ) - term3 ) )
expect_equal( targetB, myB, tol = 0.001 )
### check E-value for negative bias (shif1 only stratum M)
targetB = Eadd.est.mono$candidates$biasFactor[ Eadd.est.mono$candidates$biasDir == "negative" ]
lambda = pm_1 * (1 - fm) - pm_0 * fm
term1 = 1 / ( 2 * (fm * pm_1) )
term2 = 4 * pm_1 * pm_0 * fm * (1 - fm)
true = -RDw
term3 = true + lambda
# bias factor
( myB = term1 * ( sqrt( term3^2 + term2 ) - term3 ) )
expect_equal( targetB, myB, tol = 0.001 )
})
test_that( "IC_evalue_inner's grid search should work when needing to increase search space upper bound", {
# give it a huge IC_c so that it will have to increase the search space (i.e., B will have to be > 2)
x = evalues.IC( stat = "est",
true = 0,
unidirBias = FALSE,
p1_1 = .9,
p1_0 = .1,
n1_1 = 100,
n1_0 = 100,
f1 = 0.2,
p0_1 = 0.1,
p0_0 = 0.4,
n0_1 = 100,
n0_0 = 100,
f0 = 0.3 )
# check that eventual bias factor had to be increased above the default search limit of 2
expect_equal( x$evalues$biasFactor > 2, TRUE )
})
test_that( "RDt for each stratum should stay in [-1,1] even if passed an extreme bias factor", {
# correction should increase stratum 1 to RD=1 and should decrease
# stratum 0 to -1
( x = RDt_bound( p1_1 = .9,
p1_0 = .1,
n1_1 = 100,
n1_0 = 100,
f1 = 0.2,
maxB_1 = 50,
biasDir_1 = "negative",
p0_1 = 0.1,
p0_0 = 0.4,
n0_1 = 100,
n0_0 = 100,
f0 = 0.3,
maxB_0 = 50,
biasDir_0 = "positive" ) )
expect_equal( x$RD[ x$stratum == "1" ], 1 )
expect_equal( x$RD[ x$stratum == "0" ], -1 )
# and with strata labels flipped
( x = RDt_bound( p0_1 = .9,
p0_0 = .1,
n0_1 = 100,
n0_0 = 100,
f0 = 0.2,
maxB_0 = 50,
biasDir_0 = "negative",
p1_1 = 0.1,
p1_0 = 0.4,
n1_1 = 100,
n1_0 = 100,
f1 = 0.3,
maxB_1 = 50,
biasDir_1 = "positive" ) )
expect_equal( x$RD[ x$stratum == "1" ], -1 )
expect_equal( x$RD[ x$stratum == "0" ], 1 )
})
Try the EValue package in your browser
Any scripts or data that you put into this service are public.
EValue documentation built on Oct. 28, 2021, 9:10 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# no need to source helper_testthat.R explicitly to pass CRAN checks
# because helper files that start with "helper" automatically
# get sourced first: https://testthat.r-lib.org/reference/test_dir.html
# # for local testing:
# library(testthat)
# library(devtools)
# library(dplyr)
# library(ICC)
# library(msm)
# library(MetaUtility)
# library(here())
#
# setwd(here())
# setwd("tests")
# source("helper_testthat.R")
# load_all()
# # end part for local testing
# enter example datasets (Letenneur)
# Y: dementia
# X: low education
# n: sample size
# data for women
nw_1 = 2988
nw_0 = 364
dw = data.frame( Y = c(1, 1, 0, 0),
X = c(1, 0, 1, 0),
n = c( 158, 6, nw_1-158, nw_0-6 ) )
# data for men
nm_1 = 1790
nm_0 = 605
dm = data.frame( Y = c(1, 1, 0, 0),
X = c(1, 0, 1, 0),
n = c( 64, 17, nm_1-64, nm_0-17 ) )
# P(Y = 1 | X = 1) and P(Y = 1 | X = 0) for women and for men
( pw_1 = dw$n[ dw$X == 1 & dw$Y == 1 ] / sum(dw$n[ dw$X == 1 ]) )
( pw_0 = dw$n[ dw$X == 0 & dw$Y == 1 ] / sum(dw$n[ dw$X == 0 ]) )
( pm_1 = dm$n[ dm$X == 1 & dm$Y == 1 ] / sum(dm$n[ dm$X == 1 ]) )
( pm_0 = dm$n[ dm$X == 0 & dm$Y == 1 ] / sum(dm$n[ dm$X == 0 ]) )
# prevalence of low education among women and among men
fw = nw_1 / (nw_1 + nw_0)
fm = nm_1 / (nm_1 + nm_0)
# confounded estimates
RDw = pw_1 - pw_0
RDm = pm_1 - pm_0
test_that("RDt_var, symmetry when reversing sign of RD", {
# test for symmetry
v1 = RDt_var( f = .25,
p1 = 0.3,
p0 = 0.1,
n1 = 50,
n0 = 100,
.maxB = 2 )
# here I've changed argument order to reverse sign of RD
v2 = RDt_var( f = .25,
p0 = 0.3,
p1 = 0.1,
n0 = 50,
n1 = 100,
.maxB = 2 )
expect_equal(v1, v2)
})
test_that( "RDt_bound, symmetry when shifting strata in opposite directions", {
x1 = RDt_bound( p1_1 = 0.6,
p1_0 = 0.4,
n1_1 = 100,
n1_0 = 10,
f1 = 0.25,
maxB_1 = 2,
biasDir_1 = "positive",
p0_1 = 0.4,
p0_0 = 0.6,
n0_1 = 10,
n0_0 = 100,
f0 = .25,
biasDir_0 = "negative",
alpha = 0.05 )
expect_equal( x1$RD[1], -x1$RD[2] )
expect_equal( x1$se[1], x1$se[2] )
expect_equal( x1$lo[1], -x1$hi[2] )
expect_equal( x1$lo[2], -x1$hi[1] )
expect_equal( x1$pval[1], x1$pval[2] )
# sanity check: equality when both strata are the same and are positively biased
x1 = RDt_bound( p1_1 = 0.6,
p1_0 = 0.4,
n1_1 = 100,
n1_0 = 10,
f1 = 0.25,
maxB_1 = 1.6,
biasDir_1 = "positive",
p0_1 = 0.6,
p0_0 = 0.4,
n0_1 = 100,
n0_0 = 10,
f0 = .25,
maxB_0 = 1.6,
biasDir_0 = "positive",
alpha = 0.05 )
expect_equal( x1$RD[1], x1$RD[2] )
expect_equal( x1$se[1], x1$se[2] )
expect_equal( x1$lo[1], x1$lo[2] )
expect_equal( x1$hi[1], x1$hi[2] )
expect_equal( x1$pval[1], x1$pval[2] )
# sanity check: equality when both strata are the same and are NEGATIVELY biased
x1 = RDt_bound( p1_1 = 0.6,
p1_0 = 0.4,
n1_1 = 100,
n1_0 = 10,
f1 = 0.25,
maxB_1 = 1.6,
biasDir_1 = "negative",
p0_1 = 0.6,
p0_0 = 0.4,
n0_1 = 100,
n0_0 = 10,
f0 = .25,
maxB_0 = 1.6,
biasDir_0 = "negative",
alpha = 0.05 )
expect_equal( x1$RD[1], x1$RD[2] )
expect_equal( x1$se[1], x1$se[2] )
expect_equal( x1$lo[1], x1$lo[2] )
expect_equal( x1$hi[1], x1$hi[2] )
expect_equal( x1$pval[1], x1$pval[2] )
})
test_that( "Bound from RDt_bound should agree with closed form in paper", {
B = 2.3
term1 = (pw_1 - pw_0 * B) * ( fw + (1 - fw) / B )
term2 = (pm_1 * B - pm_0) * ( fm + (1 - fm) / B )
mine = term1 - term2
x = RDt_bound( p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
biasDir_1 = "positive",
maxB_1 = B,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
biasDir_0 = "negative",
maxB_0 = B )
expect_equal( x$RD[ x$stratum == "effectMod" ], mine, tol = 0.0001 )
} )
test_that( "E-value for 1 stratum from IC_evalue should match R package", {
( evalueEst = IC_evalue_inner( stratum = "1",
varName = "RD",
true = 0,
unidirBias = FALSE,
p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
alpha = 0.05 )$evalues )
( evalueCI = IC_evalue_inner( stratum = "1",
varName = "lo",
true = 0,
unidirBias = FALSE,
p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
alpha = 0.05 )$evalues )
# now try against package:
evalueOld = EValue::evalues.RD( n11 = dw$n[ dw$X == 1 & dw$Y == 1 ],
n10 = dw$n[ dw$X == 1 & dw$Y == 0 ],
n01 = dw$n[ dw$X == 0 & dw$Y == 1 ],
n00 = dw$n[ dw$X == 0 & dw$Y == 0 ],
true = 0,
alpha = 0.05)
# they agree! :D
# woohoo!!!!!
expect_equal( evalueOld$est.Evalue, evalueEst$evalue, tol = 0.001 )
expect_equal( evalueOld$lower.Evalue, evalueCI$evalue, tol = 0.001 )
})
test_that( "Bound from RDt_bound should be symmetric after flipping signs of both RDs", {
# shift each stratum by different amount
B1 = 1.5
B0 = 1.2
# both RDs > 0
( x = RDt_bound( p1_1 = .8,
p1_0 = .6,
n1_1 = 20,
n1_0 = 200,
f1 = 0.5,
biasDir_1 = "positive",
maxB_1 = B1,
p0_1 = .8,
p0_0 = .75,
n0_1 = 30,
n0_0 = 40,
f0 = 0.5,
biasDir_0 = "positive",
maxB_0 = B0 ) )
# both RDs < 0
# strata labels reversed so that the interaction contrast itself remains positive
( x2 = RDt_bound( p0_1 = .6,
p0_0 = .8,
n0_1 = 200,
n0_0 = 20,
f0 = 0.5,
biasDir_0 = "negative",
maxB_0 = B1,
p1_1 = .75,
p1_0 = .8,
n1_1 = 40,
n1_0 = 30,
f1 = 0.5,
biasDir_1 = "negative",
maxB_1 = B0 ) )
# interaction contrast should be exactly the same
expect_equal( x[3, 3:6], x2[3, 3:6] )
# strata are reversed and signs are reversed
# inference stays same
# but CIs flip
expect_equal( as.numeric( x[1, c("se","pval")] ), as.numeric( x2[2, c("se", "pval")] ) )
expect_equal( as.numeric( x[2, c("se","pval")] ), as.numeric( x2[1, c("se", "pval")] ) )
expect_equal( x[1, "RD"], -x2[2, "RD"] )
expect_equal( x[2, "RD"], -x2[1, "RD"] )
})
test_that( "E-value should be the same when flipping strata signs", {
x = evalues.IC( stat = "est",
true = 0.1,
unidirBias = FALSE,
p1_1 = .5,
p1_0 = .3,
n1_1 = 100,
n1_0 = 20,
f1 = .6,
p0_1 = .4,
p0_0 = .35,
n0_1 = 40,
n0_0 = 60,
f0 = .2,
alpha = 0.05 )$evalues
x2 = evalues.IC( stat = "est",
true = 0.1,
unidirBias = FALSE,
p0_1 = .3,
p0_0 = .5,
n0_1 = 20,
n0_0 = 100,
f0 = .6,
p1_1 = .35,
p1_0 = .4,
n1_1 = 60,
n1_0 = 40,
f1 = .2,
alpha = 0.05 )$evalues
expect_equal( x$evalue, x2$evalue, tol = 0.001 )
expect_equal( x$biasFactor, x2$biasFactor, tol = 0.001 )
} )
test_that("evalues.IC should warn if E-value is 1", {
# E-value for estimate is 1
expect_message( evalues.IC( stat = "est",
true = .5,
unidirBias = FALSE,
p1_1 = .5,
p1_0 = .4,
n1_1 = 100,
n1_0 = 100,
f1 = .5,
p0_1 = .5,
p0_0 = .42,
n0_1 = 100,
n0_0 = 100,
f0 = .5,
alpha = 0.05 ) )
# E-value for CI is 1
# first check the CI limit
RDt_bound( p1_1 = .5,
p1_0 = .4,
n1_1 = 100,
n1_0 = 100,
f1 = .5,
p0_1 = .5,
p0_0 = .42,
n0_1 = 100,
n0_0 = 100,
f0 = .5,
# no bias
maxB_1 = 1,
maxB_0 = 1,
biasDir_1 = "positive",
biasDir_0 = "positive" )
expect_message( evalues.IC( stat = "CI",
true = .5,
unidirBias = FALSE,
p1_1 = .5,
p1_0 = .4,
n1_1 = 100,
n1_0 = 100,
f1 = .5,
p0_1 = .5,
p0_0 = .42,
n0_1 = 100,
n0_0 = 100,
f0 = .5,
alpha = 0.05 ) )
# also when trying both candidates (non-unidir, unknown)
expect_message( evalues.IC( stat = "CI",
true = 0.0,
unidirBias = TRUE,
unidirBiasDirection = "unknown",
p1_1 = .9,
p1_0 = .8,
n1_1 = 100,
n1_0 = 100,
f1 = .9,
p0_1 = .5,
p0_0 = .42,
n0_1 = 100,
n0_0 = 100,
f0 = .1,
alpha = 0.05 ) )
})
test_that("evalues.IC should reject bad input", {
# if unidirBias is TRUE, must provide unidirBiasDirection
expect_error( evalues.IC( stat = "est",
true = 0,
unidirBias = TRUE,
p1_1 = .4,
p1_0 = .4,
n1_1 = 100,
n1_0 = 100,
f1 = .5,
p0_1 = .5,
p0_0 = .42,
n0_1 = 100,
n0_0 = 100,
f0 = .5,
alpha = 0.05 ) )
# specified unidirBias = FALSE, so the argument unidirBiasDirection will be ignored
expect_warning( evalues.IC( stat = "est",
true = 0,
unidirBias = FALSE,
unidirBiasDirection = "positive",
p1_1 = .5,
p1_0 = .4,
n1_1 = 100,
n1_0 = 100,
f1 = .5,
p0_1 = .5,
p0_0 = .42,
n0_1 = 100,
n0_0 = 100,
f0 = .5,
alpha = 0.05 ) )
expect_warning( evalues.IC( stat = "est",
true = 0,
unidirBias = FALSE,
unidirBiasDirection = "negative",
p1_1 = .5,
p1_0 = .4,
n1_1 = 100,
n1_0 = 100,
f1 = .5,
p0_1 = .5,
p0_0 = .42,
n0_1 = 100,
n0_0 = 100,
f0 = .5,
alpha = 0.05 ) )
} )
test_that( "E-value from evalues.IC should be the solution to RDt_bound and should agree with theory in paper", {
Eadd.est = evalues.IC( stat = "est",
true = 0,
unidirBias = FALSE,
p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
alpha = 0.05 )$evalues
# pass the resulting bias factor (i.e., equivalent to alleged E-value) to each stratum
x = RDt_bound( p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
biasDir_1 = "positive",
# from above test
maxB_1 = Eadd.est$biasFactor,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
biasDir_0 = "negative",
maxB_0 = Eadd.est$biasFactor )
expect_equal( x$RD[ x$stratum == "effectMod" ], 0, tol = 0.0001 )
# should agree with closed form
gamma = pw_0 * (1 - fw) - pw_1 * fw + pm_1 * (1 - fm) - pm_0 * fm
term1 = 1 / ( 2 * (fw * pw_0 + fm * pm_1) )
term2 = 4 * (pw_0 * fw + pm_1 * fm) * ( pw_1 * (1 - fw) + pm_0 * (1-fm) )
# bias factor
( my.Badd.est = term1 * ( sqrt( gamma^2 + term2 ) - gamma ) )
expect_equal( Eadd.est$biasFactor, my.Badd.est, tol = 0.001 )
# E-value
expect_equal( Eadd.est$evalue, g(my.Badd.est), tol = 0.001 )
})
test_that("evalues.IC solution for one stratum should agree with existing fn, evalues.RD", {
( resNonMono = evalues.IC( stat = "est",
true = 0,
unidirBias = FALSE,
p1_1 = .6,
p1_0 = .4,
n1_1 = 100,
n1_0 = 100,
f1 = .5,
p0_1 = 0.3,
p0_0 = 0.20,
n0_1 = 100,
n0_0 = 100,
f0 = .5,
alpha = 0.05 )$evalues )
# RDt_bound
RDs2 = RDt_bound( p1_1 = .6,
p1_0 = .4,
n1_1 = 100,
n1_0 = 100,
f1 = .5,
biasDir_1 = "positive",
maxB_1 = resNonMono$biasFactor,
p0_1 = 0.3,
p0_0 = 0.20,
n0_1 = 100,
n0_0 = 100,
f0 = .5,
biasDir_0 = "negative",
maxB_0 = resNonMono$biasFactor )
# stratum 1 only
evalueOld2 = evalues.RD( n11 = 100 * (0.6),
n10 = 100 * (1-0.6),
n01 = 100 * 0.4,
n00 = 100 * (1-0.4),
true = RDs2$RD[1],
alpha = 0.05)
expect_equal( evalueOld2$est.Evalue, resNonMono$evalue, tol = 0.001 )
# # stratum M only
# evalueOld2 = EValue::evalues.RD( n11 = 100 * (0.3),
# n10 = 100 * (1-0.3),
# n01 = 100 * 0.2,
# n00 = 100 * (1-0.2),
# true = RDs2$RD[2],
# alpha = 0.05)
#
#
# # they agree! :D
# # woohoo!!!!!
# expect_equal( evalueOld2$est.Evalue, resNonMono$evalue, tol = 0.001 )
})
test_that("Evalue candidate (negative unidir bias) should successfully move RDm up to RDw", {
# superassign because we'll use for subsequent tests
( Eadd.est.mono <<- evalues.IC( stat = "est",
true = 0,
unidirBias = TRUE,
unidirBiasDirection = "unknown",
p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
alpha = 0.05 ) )
x = RDt_bound( p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
biasDir_1 = "negative",
maxB_1 = 1,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
biasDir_0 = "negative",
maxB_0 = Eadd.est.mono$candidates$biasFactor[ Eadd.est.mono$candidates$biasDir == "negative" ] )
expect_equal( x$RD[ x$stratum == "0" ], RDw, tol = 0.0001 )
# and likewise for CI limit
( Eadd.CI.mono <<- evalues.IC( stat = "CI",
true = 0,
unidirBias = TRUE,
unidirBiasDirection = "unknown",
p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
alpha = 0.05 ) )
x = RDt_bound( p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
biasDir_1 = "negative",
maxB_1 = 1,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
biasDir_0 = "negative",
maxB_0 = Eadd.CI.mono$candidates$biasFactor[ Eadd.CI.mono$candidates$biasDir == "negative" ] )
expect_equal( x$lo[ x$stratum == "effectMod" ], 0, tol = 0.0001 )
} )
test_that( "Evalue candidate (positive unidir bias) should successfully move RDw down to RDm", {
( x = RDt_bound( p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
biasDir_1 = "positive",
maxB_1 = Eadd.est.mono$candidates$biasFactor[ Eadd.est.mono$candidates$biasDir == "positive" ],
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
biasDir_0 = "positive",
maxB_0 = 1 ) )
expect_equal( x$RD[ x$stratum == "1" ], RDm, tol = 0.0001 )
# and likewise for CI limit
( x = RDt_bound( p1_1 = pw_1,
p1_0 = pw_0,
n1_1 = nw_1,
n1_0 = nw_0,
f1 = fw,
biasDir_1 = "positive",
maxB_1 = Eadd.CI.mono$candidates$biasFactor[ Eadd.CI.mono$candidates$biasDir == "positive" ] ,
p0_1 = pm_1,
p0_0 = pm_0,
n0_1 = nm_1,
n0_0 = nm_0,
f0 = fm,
biasDir_0 = "positive",
maxB_0 = 1 ) )
expect_equal( x$lo[ x$stratum == "effectMod" ], 0, tol = 0.0001 )
} )
test_that( "E-values from IC_evalue (grid search) should match closed form in paper", {
B = 4
true = ( pm_1 * B - pm_0 ) * ( fm + (1-fm) / B )
# suggestively name terms as in quadratic formula
termA = fm * pm_1
termB = pm_1 * ( 1 - fm ) - fm * pm_0 - true
termC = -pm_0 * (1 - fm)
# this is the polynomial that is to be solved for E-value
( mine = termA * B^2 + termB * B + termC )
expect_equal( mine, 0, tol = 0.0001 )
# check E-value for negatively biased stratum in paper
# this is the one that arises from reversing roles and signs in the existing E-value on
# Ding Appendix, pg 18
#bm
### check E-value for positive bias (shift only stratum W)
# check against Ding Appendix, pg 18 (Prop A.11)
targetB = Eadd.est.mono$candidates$biasFactor[ Eadd.est.mono$candidates$biasDir == "positive" ]
lambda = pw_0 * (1 - fw) - pw_1 * fw
term1 = 1 / ( 2 * (fw * pw_0) )
term2 = 4 * pw_1 * pw_0 * fw * (1 - fw)
true = RDm
term3 = true + lambda
# bias factor
( myB = term1 * ( sqrt( term3^2 + term2 ) - term3 ) )
expect_equal( targetB, myB, tol = 0.001 )
### check E-value for negative bias (shif1 only stratum M)
targetB = Eadd.est.mono$candidates$biasFactor[ Eadd.est.mono$candidates$biasDir == "negative" ]
lambda = pm_1 * (1 - fm) - pm_0 * fm
term1 = 1 / ( 2 * (fm * pm_1) )
term2 = 4 * pm_1 * pm_0 * fm * (1 - fm)
true = -RDw
term3 = true + lambda
# bias factor
( myB = term1 * ( sqrt( term3^2 + term2 ) - term3 ) )
expect_equal( targetB, myB, tol = 0.001 )
})
test_that( "IC_evalue_inner's grid search should work when needing to increase search space upper bound", {
# give it a huge IC_c so that it will have to increase the search space (i.e., B will have to be > 2)
x = evalues.IC( stat = "est",
true = 0,
unidirBias = FALSE,
p1_1 = .9,
p1_0 = .1,
n1_1 = 100,
n1_0 = 100,
f1 = 0.2,
p0_1 = 0.1,
p0_0 = 0.4,
n0_1 = 100,
n0_0 = 100,
f0 = 0.3 )
# check that eventual bias factor had to be increased above the default search limit of 2
expect_equal( x$evalues$biasFactor > 2, TRUE )
})
test_that( "RDt for each stratum should stay in [-1,1] even if passed an extreme bias factor", {
# correction should increase stratum 1 to RD=1 and should decrease
# stratum 0 to -1
( x = RDt_bound( p1_1 = .9,
p1_0 = .1,
n1_1 = 100,
n1_0 = 100,
f1 = 0.2,
maxB_1 = 50,
biasDir_1 = "negative",
p0_1 = 0.1,
p0_0 = 0.4,
n0_1 = 100,
n0_0 = 100,
f0 = 0.3,
maxB_0 = 50,
biasDir_0 = "positive" ) )
expect_equal( x$RD[ x$stratum == "1" ], 1 )
expect_equal( x$RD[ x$stratum == "0" ], -1 )
# and with strata labels flipped
( x = RDt_bound( p0_1 = .9,
p0_0 = .1,
n0_1 = 100,
n0_0 = 100,
f0 = 0.2,
maxB_0 = 50,
biasDir_0 = "negative",
p1_1 = 0.1,
p1_0 = 0.4,
n1_1 = 100,
n1_0 = 100,
f1 = 0.3,
maxB_1 = 50,
biasDir_1 = "positive" ) )
expect_equal( x$RD[ x$stratum == "1" ], -1 )
expect_equal( x$RD[ x$stratum == "0" ], 1 )
})
Try the EValue package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.