tests/testthat/test.Design.R
In markmfredrickson/RItools: Randomization Inference Tools
################################################################################
# Tests for DesignOptions objects
################################################################################
if (!exists('nreps_')) nreps_ <- 10L
library("testthat")
context("ModelMatrixPlus S4 class carries per-covariate missingness info")
test_that("If no missing data, then NotMissing is a matrix w n rows and 0 cols",{
d <- data.frame(
id = 1:500,
x = rnorm(500),
cluster = rep(1:100, 5),
strata.good = rep(c(1:4, NA), 100),
strata.bad = sample(1:100, size = 500, replace = T),
z.good = rep(c(0,1), 250),
z.bad = sample(c(0,1), size = 500, replace = T))
d$'(weights)' = 1 # meet expectation of a weights column
simple <- RItools:::model_matrix(z.good ~ x, data = d)
expect_equivalent(dim(simple@NotMissing), c(500,1))
})
test_that("Missingness gets passed through in Covariates, recorded in NotMissing",{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=rep(c(NA, 1), c(3,2)),
x2=c(1:5),
fac=factor(c(rep(1:2,2), NA))
)
datmf <- model.frame(z ~ x1 + x2 + fac, dat, na.action = na.pass)
datmf$'(weights)' <- 1
simple2 <- RItools:::model_matrix(z ~ x1 + x2 + fac, data = datmf)
expect_equivalent(ncol(simple2@NotMissing), 3)
expect_equivalent(colnames(simple2@NotMissing), c("_any Xs recorded_", "x1", "fac"))
})
test_that("model_matrix offered missing or negative weights",{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
w_with_NAs=rep(c(NA, 1), c(3,2)),
w_with_negatives=-2:2,
x=c(1:5),
fac=factor(c(rep(1:2,2), NA))
)
datmf_with_NAs <-
model.frame(z ~ x + fac, dat, na.action = na.pass,
weights=w_with_NAs)
datmf_with_negatives <-
model.frame(z ~ x + fac, dat, na.action = na.pass,
weights=w_with_negatives)
expect_error(RItools:::model_matrix(z ~ x + fac,
data = datmf_with_NAs),
"uweights")
expect_error(RItools:::model_matrix(z ~ x + fac,
data = datmf_with_negatives),
"uweights")
})
test_that("lookup tables OK, even w/ complex & multi-column terms",{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=rep(c(NA, 1), c(3,2)),
x2=c(1:5),
fac=factor(c(rep(1:2,2), NA))
)
datmf <- model.frame(z ~ x1 + x2 + fac, dat, na.action = na.pass)
datmf$'(weights)' <- 1
simple2 <- RItools:::model_matrix(z ~ x1 + x2 + fac, data=datmf)
expect_equal(simple2@OriginalVariables, 1:3)
expect_equal(simple2@TermLabels, c( "x1", "x2", "fac"))
expect_equal(simple2@NM.Covariates, c(2,0,3))
expect_equal(simple2@NM.terms, c(2,0,3))
## check that complex term don't spell trouble in themselves
datmf <- model.frame(z ~ x1 + cut(x2, c(0,3,6)) + fac, data = dat,
na.action = na.pass)
datmf$'(weights)' <- 1
simple3 <- RItools:::model_matrix(z ~ x1 + cut(x2, c(0,3,6)) + fac, data = datmf)
expect_equal(simple3@OriginalVariables, 1:3)
expect_equal(simple3@TermLabels, c("x1", "cut(x2, c(0, 3, 6))", "fac"))
expect_equal(simple3@NM.Covariates, c(2,0,3))
expect_equal(simple3@NM.terms, c(2,0,3))
## now try a complex term that actually expands to multiple columns
datmf <- model.frame(z ~ x1 + cut(x2, c(0,3,6)) + fac, data = dat,
na.action = na.pass)
datmf$'(weights)' <- 1
simple4 <- RItools:::model_matrix(z ~ x1 + cut(x2, c(0,3,6)) + fac, data = datmf,
contrasts=list("cut(x2, c(0, 3, 6))"=diag(2)))
expect_equal(simple4@OriginalVariables, c(1,2,2,3))
expect_equal(simple4@TermLabels, c("x1", "cut(x2, c(0, 3, 6))", "fac"))
expect_equal(simple4@NM.Covariates, c(2,0,0,3))
expect_equal(simple4@NM.terms, c(2,0,3))
## now put NAs in the multi-column complex term
datmf <- model.frame(z ~ x2 + cut(x1, c(0,3,6)) + fac, data = dat,
na.action = na.pass)
datmf$'(weights)' <- 1
simple5 <- RItools:::model_matrix(z ~ x2 + cut(x1, c(0,3,6)) + fac, data = datmf,
contrasts=list("cut(x1, c(0, 3, 6))"=diag(2)))
expect_equal(simple5@OriginalVariables, c(1,2,2,3))
expect_equal(simple5@TermLabels, c("x2", "cut(x1, c(0, 3, 6))", "fac"))
expect_equal(simple5@NM.Covariates, c(0,2,2,3))
expect_equal(simple5@NM.terms, c(0,2,3))
})
test_that("Issue #76: Using I() in formulas", {
x <- data.frame(z=c(1,1)) # have to exclude
while (all(x$z==x[1L,'z'])) # degenerate case
x <- data.frame(x = rnorm(10), y = rnorm(10), z = rbinom(10, size = 1, p = 1/3))
x$"(weights)" <- 1
d <- makeDesigns(z ~ I(x * sin(y)), data = x)
expect_s4_class(d, "DesignOptions")
## While we're at it, confirm that the non-stratification
## encoded in this DesignOptions bears the column name "--".
})
test_that("Null stratification is encoded by '--'",{
data(nuclearplants)
nuclearplants$"(weights)" <- 1
d0 <- makeDesigns(pr ~ cost, data=nuclearplants)
expect_equal(colnames(d0@StrataFrame), "--")
foo <- nuclearplants$pt
d1 <- makeDesigns(pr ~ cost + strata(foo), data=nuclearplants)
expect_true("--" %in% colnames(d1@StrataFrame))
})
test_that("Issue #86: makeDesigns finds variables outside data arg",{
data(nuclearplants)
foo <- nuclearplants$pt
nuclearplants$"(weights)" <- 1
d <- makeDesigns(pr ~ cost + strata(foo), data=nuclearplants)
expect_s4_class(d, "DesignOptions")
})
test_that("Duplicated missingness patterns handled appropriately",{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=c(rep(NA, 3), 1:2),
x2=c(1:5),
fac=factor(c(rep(1:2,2), NA))
)
datmf <- model.frame(z ~ x1 + I(x1^2) + fac, data=dat, na.action = na.pass)
datmf$'(weights)' <- 1
simple2 <- RItools:::model_matrix(z ~ x1 + I(x1^2) + fac, data = datmf)
expect_equal(simple2@OriginalVariables, 1:3)
expect_equal(simple2@TermLabels, c("x1", "I(x1^2)", "fac"))
expect_equal(simple2@NM.Covariates, 1+c(1,1,2))
expect_equal(simple2@NM.terms, 1+c(1,1,2))
## If exactly two terms have missing data but in the same pattern, then
## NotMissing is a matrix w/ n rows and 1 col.
datmf <- model.frame(z ~ x1 + I(x1^2), data=dat, na.action = na.pass)
datmf$'(weights)' <- 1
simple3 <- RItools:::model_matrix(z ~ x1 + I(x1^2), data = datmf)
expect_equal(simple3@OriginalVariables, 1:2)
expect_equal(simple3@TermLabels, c("x1", "I(x1^2)"))
expect_equal(ncol(simple3@NotMissing), 1)
expect_equal(simple3@NM.Covariates, c(1,1))
expect_equal(simple3@NM.terms, c(1,1))
})
test_that("All-Xes missingness |-> NotMissing col '_any Xs recorded_'",{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=c(rep(NA, 4), 2),
x2=c(1:3, NA, 5),
fac=factor(c(NA, 1:2, NA, 2))
)
dat$'(weights)' <- 1
datmf <- model.frame(z ~ x1 + I(x1^2) + fac, data=dat, na.action = na.pass)
datmf$'(weights)' <- 1
simple6 <- RItools:::model_matrix(z ~ x1 + I(x1^2) + fac, data = datmf)
expect_equal(simple6@OriginalVariables, 1:3)
expect_equal(simple6@TermLabels, c("x1", "I(x1^2)", "fac"))
expect_equal(colnames(simple6@NotMissing)[1], "_any Xs recorded_")
expect_equal(simple6@NM.Covariates, c(2,2,1))
expect_equal(simple6@NM.terms, c(2,2,1))
simple7 <- makeDesigns(z ~ x1 + I(x1^2) + fac + 0 + strata(strat) + cluster(clus),
dat)
simple7 <- as(simple7, "StratumWeightedDesignOptions")
simple7@Sweights <-
RItools:::DesignWeights(simple7,
RItools:::effectOfTreatmentOnTreated)
## As writing of this test, DesignWeights() expects only pre-aggregated designs,
## and infers treatment:control ratios from the numbers of elements in in each
## of condition (by strata), not the number of clusters. Thus in this example
## it should believe that the ETT weights are proportional to 2 for stratum a,
## 1 for stratum b:
expect_equivalent(simple7@Sweights$strat$sweights, (2:1)/3)
## key point: had missingness of each of unit 4's covariates tricked it into ignoring
## that observation, then the ETT weight for stratum b would have been proportional to 0,
## not 1.
} )
test_that("All-Xes missingness noted in descriptives", {
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=c(rep(NA, 4), 2),
x2=c(1:3, NA, 5),
fac=factor(c(NA, 1:2, NA, 2))
)
dat$'(weights)' <- 1
simple5 <- makeDesigns(z ~ x1 + I(x1^2) + fac + cluster(clus), dat)
descr <- RItools:::designToDescriptives(simple5)
expect_equivalent(descr['(_any Xs recorded_)' ,1:2,1], c(1, 1/3)) # not just 1s all around
})
context("DesignOptions objects")
test_that("Creating DesignOptions objects", {
set.seed(20130801)
replicate(nreps_,{
d <- data.frame(
id = 1:500,
x = rnorm(500),
cluster = rep(1:100, 5),
strata.good = rep(c(1:4, NA), 100),
strata.bad = sample(1:100, size = 500, replace = T),
z.good = rep(c(0,1), 250),
z.bad = sample(c(0,1), size = 500, replace = T))
d$'(weights)' = 1 # meet expectation of a weights column
# checking input
expect_error(RItools:::makeDesigns(~ x, data = d), "treatment")
expect_error(RItools:::makeDesigns(z.bad ~ x + strata(strata.good) + cluster(cluster), data = d), "cluster")
expect_error(RItools:::makeDesigns(z.good ~ x + strata(strata.bad) + cluster(cluster), data = d), "strata")
expect_error(RItools:::makeDesigns(z.good ~ x + cluster(id) + cluster(cluster), data = d), "cluster")
expect_error(RItools:::makeDesigns(z.good ~ x - 1, data = d, "stratification"))
# actually testing that the output is as expected
simple <- RItools:::makeDesigns(z.good ~ x, data = d)
expect_equal(dim(simple@StrataFrame)[2], 1)
expect_equivalent(simple@Covariates[, "x"], d$x)
expect_equivalent(simple@Z, as.logical(d$z.good))
expect_equal(nlevels(simple@Cluster), 500) # a cluster per individual
clustered <- RItools:::makeDesigns(z.good ~ x + cluster(cluster), data = d)
expect_equal(dim(clustered@StrataFrame)[2], 1)
expect_true(nlevels(clustered@Cluster) > 1)
clustStrata <- RItools:::makeDesigns(z.good ~ x + cluster(cluster) + strata(strata.good), data = d)
expect_equal(dim(clustStrata@StrataFrame)[2], 2)
clustStrata.c <- RItools:::makeDesigns(z.good ~ x + cluster(cluster) + strata(strata.good, strata.good), data = d)
expect_equivalent(clustStrata, clustStrata.c)
# dropping the overall comparison
expect_equal(dim(RItools:::makeDesigns(z.good ~ x + cluster(cluster) + strata(strata.good) - 1, data = d)@StrataFrame)[2], 1)
## More tests to write:
# - All NA strata variables
# - Missing z or cluster
# - strata with extra levels but no observations (which can be safely dropped)
# (NB: extra levels tested upstream, in xBalance, as of commit 34861515;
# see ./test.clusters.R )
})
})
test_that("NotMissing vars correctly generated",
{
replicate(nreps_,{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=rep(c(NA, TRUE), c(3,2)),
x2 = c(1:5),
x3=c(TRUE, FALSE, NA, TRUE, FALSE),
fac=factor(c(rep(1:2,2), NA))
)
dat$'(weights)' <- 1
simple <- RItools:::makeDesigns(z ~ x1 + x2 + strata(strat) + cluster(clus), data = dat)
expect_match(colnames(simple@NotMissing), "x1", all=FALSE)
expect_false(any(grepl("x2", colnames(simple@NotMissing))))
expect_false(any(grepl("TRUE", colnames(simple@Covariates))))
expect_false(any(grepl("FALSE", colnames(simple@Covariates))))
simple2 <- RItools:::makeDesigns(z ~ x1 + x2 + fac+ strata(strat) + cluster(clus), data = dat)
expect_match(colnames(simple2@NotMissing), "x1", all=FALSE)
expect_false(any(grepl("x2", colnames(simple2@NotMissing))))
expect_match(colnames(simple2@NotMissing), "fac", all=FALSE)
expect_false(any(grepl("TRUE", colnames(simple2@Covariates))))
expect_false(any(grepl("FALSE", colnames(simple2@Covariates))))
simple3 <- RItools:::makeDesigns(z ~ x1 + x3 + fac+ strata(strat) + cluster(clus), data = dat)
expect_match(colnames(simple3@NotMissing), "x1", all=FALSE)
expect_match(colnames(simple3@NotMissing), "x3", all=FALSE)
expect_match(colnames(simple3@NotMissing), "fac", all=FALSE)
expect_false(any(grepl("TRUE", colnames(simple3@Covariates))))
expect_false(any(grepl("FALSE", colnames(simple3@Covariates))))
})
})
test_that("Issue 88: logical Covariates correctly generated",
{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=rep(c(NA, TRUE), c(3,2)),
x2 = c(1:5),
x3=c(TRUE, FALSE, NA, TRUE, FALSE),
fac=factor(c(rep(1:2,2), NA))
)
dat$'(weights)' <- 1
simple1 <- RItools:::makeDesigns(z ~ x1 + x2, data = dat)
expect_false(any(grepl("TRUE", colnames(simple1@Covariates))))
expect_false(any(grepl("FALSE", colnames(simple1@Covariates))))
simple2 <- RItools:::makeDesigns(z ~ x1 + x2 + strata(strat), data = dat)
expect_false(any(grepl("TRUE", colnames(simple2@Covariates))))
expect_false(any(grepl("FALSE", colnames(simple2@Covariates))))
simple3 <- RItools:::makeDesigns(z ~ x1 + x2 + strata(strat) - 1, data = dat)
expect_false(any(grepl("TRUE", colnames(simple3@Covariates))))
expect_false(any(grepl("FALSE", colnames(simple3@Covariates))))
})
test_that("DesignOptions to descriptive statistics", {
set.seed(20130801)
replicate(nreps_,{
d <- data.frame(
x = rnorm(500),
f = factor(sample(c("A", "B", "C"), size = 500, replace = T)),
c = rep(1:100, 5),
s = rep(c(1:4, NA), 100),
z = rep(c(0,1), 250))
d$'(weights)' = 1 # meet expectation of a weights column
simple <- RItools:::makeDesigns(z ~ x + f + strata(s) + cluster(c), data = d)
descriptives <- RItools:::designToDescriptives(simple)
expect_is(descriptives, "array")
expect_equal(dim(descriptives), c(5, 5, 2))
# descriptives ignore clustering
design.noclus <- RItools:::makeDesigns(z ~ x + f + strata(s), data = d)
expect_equal(descriptives, RItools:::designToDescriptives(design.noclus))
# the strata should imply different stats
expect_false(identical(descriptives[,,1], descriptives[,,2]))
# however, the pooled s.d.s should be the same.
expect_identical(descriptives[,"pooled.sd","--"], descriptives[,"pooled.sd","s"])
# ok, now checking that values are correct.
expect_equal(mean(d$x[d$z == 1]), descriptives["x", "Treatment", "--"])
expect_equal(mean(d$x[d$z == 0]), descriptives["x", "Control", "--"])
# with equal sized strata, the the control/treatment means are the means of the the strata means
expect_equal(mean(tapply(d$x[d$z == 1], d$s[d$z == 1], mean)), descriptives["x", "Treatment", "s"])
expect_equal(mean(tapply(d$x[d$z == 0], d$s[d$z == 0], mean)), descriptives["x", "Control", "s"])
})
})
test_that("designToDescriptives uses provided covariate scales",{
d <- data.frame(x=rnorm(50), z=rep(c(0,1), 25))
d$'(weights)' <- 1
simple <- RItools:::makeDesigns(z ~ x, data=d)
sd_x <- sd(resid(lm(x ~z, data=d)))
descriptives <- designToDescriptives(simple, covariate.scales=c(x=sd_x*10))
expect_equal(descriptives["x","pooled.sd","--"], sd_x*10)
descriptives <- designToDescriptives(simple,
covariate.scales=c(x=sd_x*10, y=Inf))
expect_equal(descriptives["x","pooled.sd","--"], sd_x*10)
expect_warning(designToDescriptives(simple, covariate.scales=sd_x),
"name")
expect_warning(designToDescriptives(simple, covariate.scales=c(x="foo")),
"numeric")
})
test_that("descriptives for NotMissing variables",
{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=rep(c(NA, 1), c(3,2)),
x2=c(1:5),
fac=factor(c(rep(1:2,2), NA))
)
dat$'(weights)' <- 1
simple <- RItools:::makeDesigns(z ~ x1 + x2 + strata(strat) + cluster(clus), data = dat)
expect_false(any(grepl("NA", colnames(simple@Covariates))))
dsimple <- RItools:::designToDescriptives(simple)
expect_match(dimnames(dsimple)[[1]], "(x1)", all=FALSE)
expect_false(any(grepl("(x2)", dimnames(dsimple)[[1]], fixed=TRUE)))
simple2 <- RItools:::makeDesigns(z ~ x1 + x2 + fac+ strata(strat) + cluster(clus), data = dat)
expect_false(any(grepl("NA", colnames(simple2@Covariates))))
dsimple2 <- RItools:::designToDescriptives(simple2)
expect_match(dimnames(dsimple2)[[1]], "(x1)", all=FALSE)
expect_match(dimnames(dsimple2)[[1]], "(fac)", all=FALSE)
expect_false(any(grepl("(x2)", dimnames(dsimple2)[[1]], fixed=TRUE)))
}
)
test_that("Issue 36: Descriptives with NAs, appropriate weighting", {
### Descriptives with missing covariates ###
set.seed(20130801)
replicate(nreps_,{
d <- data.frame(
x = rnorm(500),
f = factor(sample(c("A", "B", "C"), size = 500, replace = T)),
c = rep(1:100, 5),
s = rep(c(1:4, NA), 100),
paired = rep(c(0,1), each = 250),
z = rep(c(0,1), 250))
d$'(weights)' = 1
d.missing <- d
d.missing$x[d.missing$x < -1] <- NA
simple.all <- RItools:::makeDesigns(z ~ x + f + strata(s) + cluster(c), data = d)
descriptives.all <- RItools:::designToDescriptives(simple.all)
expect_equal(descriptives.all["x", "Treatment", "--"], mean(d$x[d$z == 1]))
expect_equal(descriptives.all["x", "Treatment", "s"], mean(d$x[d$z == 1 & !is.na(d$s)]))
simple.missing <- RItools:::makeDesigns(z ~ x + f + strata(s) + cluster(c), data = d.missing)
descriptives.missing <- RItools:::designToDescriptives(simple.missing)
with(d.missing,
expect_equal(descriptives.missing["x", "Treatment", "--"],
mean(x[z == 1], na.rm = TRUE)))
with(d.missing,
expect_equal(descriptives.missing["x", "Treatment", "s"],
mean(x[z == 1 & !is.na(s)], na.rm = TRUE)))
expect_false(identical(descriptives.all, descriptives.missing))
# ETT weighting
design.paired <- RItools:::makeDesigns(z ~ x + f + strata(paired) + strata(s), data = d)
descriptives.paired <- RItools:::designToDescriptives(design.paired)
with(d, expect_equal(descriptives.paired["x", "Control", "paired"], mean(x[z == 0])))
with(d, expect_equal(descriptives.paired["x", "Control", "--"], mean(x[z == 0])))
with(d, expect_false(identical(descriptives.paired["x", "Control", "s"], mean(x[z == 0]))))
})
})
test_that("Aggegating designs by clusters", {
set.seed(20130801)
replicate(nreps_,{
d <- data.frame(
x = rnorm(500),
f = factor(sample(c("A", "B", "C"), size = 500, replace = T)),
c = rep(1:100, 5),
s = rep(c(1:4, NA), 100),
z = rep(c(0,1), 250))
d$'(weights)' <- 1
# grab a bunch of rows and duplicate them
d <- rbind(d, d[sample(1:dim(d)[1], size = 100), ])
# swapping around the data to make sure order doesn't mask bugs
d <- d[sample(1:dim(d)[1]), ]
design <- RItools:::makeDesigns(z ~ x + f + strata(s) + cluster(c), data = d)
aggDesign <- RItools:::aggregateDesigns(design)
# one row per cluster, with columns x, fa, fb, fc
expect_equal(dim(aggDesign@Covariates), c(100, 4))
# now spot check some cluster totals of totals
expect_equal(aggDesign@Covariates[1, ], colMeans(design@Covariates[design@Cluster == 1,]))
# Z's roll up as they should
Zs <- tapply(design@Z, design@Cluster, mean)
dim(Zs) <- NULL
Zs <- as.logical(Zs)
expect_equivalent(aggDesign@Z, Zs)
# extraneous levels in the Cluster slot are ignored
design2 <- design
levels(design2@Cluster) <- c(levels(design2@Cluster), letters)
aggDesign2 <- RItools:::aggregateDesigns(design2)
expect_equal(dim(aggDesign2@Covariates), c(100, 4))
})
})
test_that("aggregateDesigns treats NA covariates as 0's" ,{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x=rep(c(NA, 1), c(3,2))
)
dat$'(weights)' <- 1
design <- RItools:::makeDesigns(z~x+strata(strat)+cluster(clus)-1, dat)
aggDesign <- RItools:::aggregateDesigns(design)
expect_equivalent(aggDesign@Covariates[,'x'], c(0, 0, 1, 1) )
nm.column.for.x <- aggDesign@NM.Covariates[match( 'x', colnames(aggDesign@Covariates))]
expect_equal(aggDesign@NotMissing[,nm.column.for.x], c(0,0,1,1) )
})
test_that("Aggregation of unit weights to cluster level",{
set.seed(20130801)
replicate(nreps_,{
d.short <- data.frame(
x = rnorm(500),
f = factor(sample(c("A", "B", "C"), size = 500, replace = T)),
c = rep(1:100, 5),
s = rep(c(1:4, NA), 100),
z = rep(c(0,1), 250))
## grab a bunch of rows and duplicate them
newrows <- c(1:nrow(d.short), sample(1:nrow(d.short), size=100, replace=T) )
d.tall <- d.short[newrows,]
d.tall$'(weights)' <- 1
d.short$'(weights)' <- as.vector(table(newrows))
design.tall <- RItools:::makeDesigns(z ~ x + f + strata(s) + cluster(c), data = d.tall)
aggDesign.tall <- RItools:::aggregateDesigns(design.tall)
design.short <- RItools:::makeDesigns(z ~ x + f + strata(s) + cluster(c), data = d.short)
aggDesign.short <- RItools:::aggregateDesigns(design.short)
## we should wind up in the same place.
expect_equal(aggDesign.tall, aggDesign.short)
d2 <- d.tall
d2$'(weights)' <- 2
design.d2 <- RItools:::makeDesigns(z ~ x + f + strata(s) + cluster(c), data = d2)
aggDesign.d2 <- RItools:::aggregateDesigns(design.d2)
expect_equal(2*aggDesign.tall@UnitWeights,aggDesign.d2@UnitWeights)
expect_equal(aggDesign.tall@NotMissing, aggDesign.d2@NotMissing)
expect_equal(aggDesign.tall@Covariates, aggDesign.d2@Covariates)
})
})
context("ModelMatrixPlus S4 class enriches model.matrix-value 'class' ")
test_that("R core hasn't revised conventions we may depend on",
{
ff <- log(Volume) ~ log(Height) + log(Girth)
m <- model.frame(ff, trees)
mm <- model.matrix(ff, m)
expect_equal("(Intercept)", colnames(mm)[1])
expect_equivalent(mm[,1,drop=TRUE], rep(1, nrow(mm)))
expect_false("(Intercept)" %in% colnames(model.matrix(update(ff, .~.-1), m)))
expect_equal(dim(mm[,integer(0)]), c(nrow(mm), 0))
expect_equal(complete.cases(mm[,integer(0)]), rep(TRUE, nrow(mm)))
expect_equivalent(weighted.mean(c(1:4, NA), c(rep(1,4), 0)), 2.5)
})
test_that("Model matrix material is properly formed",
{
ff <- log(Volume) ~ log(Height) + log(Girth)
trees1 <- trees
trees1$'(weights)' <- 1
DM0 <- model_matrix(ff, trees1, remove.intercept=FALSE)
expect_is(DM0, "ModelMatrixPlus")
m <- model.frame(ff, trees)
expect_equivalent(model.matrix(ff, m), as.matrix(DM0))
fff <- update(ff, .~.-1)
trees2 <- trees
trees2[1, "Volume"] <- NA # LHS variable, shouldn't affect anything
m2a <- model.frame(fff, trees2, na.action = na.pass)
m2a$'(weights)' <- 1
expect_equivalent(model.matrix(fff, m), as.matrix(model_matrix(fff, m2a)))
trees2[1, "Height"] <- NA # RHS variable, but still shouldn't cause rows to be dropped
m2b <- model.frame(fff, trees2, na.action = na.pass)
m2b$'(weights)' <- 1
expect_equal(dim(model.matrix(fff, m)),
dim(as.matrix(model_matrix(fff,m2b) ) )
)
## specified contrasts
dd <- data.frame(a = gl(3,4), b = gl(4,1,12)) # balanced 2-way
dd$'(weights)' <- 1
expect_equal(model.matrix(~ a + b, dd),
as.matrix(model_matrix(~ a + b, dd, remove.intercept=FALSE)))
expect_equal(model.matrix(~ a + b-1, dd),
as.matrix(model_matrix(~ a + b-1, dd)))
expect_equal(model.matrix(~ a + b, dd, contrasts = list(a = "contr.sum")),
as.matrix(model_matrix(~ a + b, dd,
contrasts = list(a = "contr.sum"), remove.intercept=FALSE))
)
expect_equal(model.matrix(~ a + b, dd,
contrasts = list(a = "contr.sum", b = "contr.poly")),
as.matrix(model_matrix(~ a + b, dd,
contrasts =
list(a = "contr.sum", b = "contr.poly"),
remove.intercept=FALSE))
)
}
)
context("alignDesignsByStrata")
test_that("alignDesigns, designToDescriptives output alignment", {
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=rep(c(NA, 1), c(3,2)),
x2=c(1:5),
fac=factor(c(rep(1:2,2), NA))
)
dat$'(weights)' <- 1
simple2 <- RItools:::makeDesigns(z ~ x1 + x2 + fac+ strata(strat) + cluster(clus), data = dat)
expect_equal(colnames(simple2@StrataFrame), c("strat", "--"))
simple2 <- as(simple2, "StratumWeightedDesignOptions")
simple2@Sweights <- RItools:::DesignWeights(simple2, # Have to aggregate 1st to figure stratum weights
RItools:::effectOfTreatmentOnTreated)
expect_true(setequal(names(simple2@Sweights), c("strat", "--")))
dsimple2 <- RItools:::designToDescriptives(simple2)
asimple2a <- RItools:::alignDesignsByStrata("--", simple2)
expect_equivalent(dimnames(dsimple2)[[1]], colnames(asimple2a@Covariates))
asimple2b <- RItools:::alignDesignsByStrata("strat", simple2)
expect_equivalent(dimnames(dsimple2)[[1]], colnames(asimple2b@Covariates))
})
test_that("alignDesigns centers covars by stratum", {
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=rep(c(NA, 1), c(3,2)),
x2=c(1:5),
fac=factor(c(rep(1:2,2), NA))
)
dat$'(weights)' <- 1
## first unweighted case
simple0 <- RItools:::makeDesigns(z ~ x1 + x2 + fac+ strata(strat) + cluster(clus), data = dat)
simple0 <- as(simple0, "StratumWeightedDesignOptions")
simple0@Sweights <- RItools:::DesignWeights(simple0, # Placeholder strat weights, shouldn't affect
RItools:::effectOfTreatmentOnTreated) # this test
asimple0u <- RItools:::alignDesignsByStrata("--",simple0)
expect_equivalent(colSums(asimple0u@Covariates),
rep(0,ncol(asimple0u@Covariates)))
asimple0s <- RItools:::alignDesignsByStrata("strat",simple0)
expect_equivalent(colSums(asimple0s@Covariates[simple0@StrataFrame[["strat"]]=="a",]),
rep(0,ncol(asimple0s@Covariates)))
expect_equivalent(as.matrix(t(asimple0s@StrataMatrix) %*% asimple0s@Covariates),
matrix(0,2,ncol(asimple0s@Covariates)))
## now with weights
dat1 <- dat
dat1$'(weights)' <- rpois(nrow(dat1), lambda=10)
while (any(dat1$'(weights)'==0)) dat1$'(weights)' <- rpois(nrow(dat1), lambda=10)
simple1 <- RItools:::makeDesigns(z ~ x1 + x2 + fac+ strata(strat) + cluster(clus), data = dat1)
simple1 <- as(simple1, "StratumWeightedDesignOptions")
simple1@Sweights <- RItools:::DesignWeights(simple1, # Placeholder strat weights, shouldn't affect
RItools:::effectOfTreatmentOnTreated) # this test
asimple1u <- RItools:::alignDesignsByStrata("--",simple1)
expect_equivalent(colSums(asimple1u@Covariates),
rep(0,ncol(asimple1u@Covariates)))
asimple1s <- RItools:::alignDesignsByStrata("strat",simple1)
tmp1 <- asimple1s@Covariates
expect_equivalent(colSums(tmp1[simple1@StrataFrame[["strat"]]=="a",]),
rep(0,ncol(asimple1s@Covariates)))
expect_equivalent(as.matrix(t(asimple1s@StrataMatrix) %*% tmp1),
matrix(0,2, ncol(asimple1s@Covariates)))
myrank <- function(x, weights) rank(x)
## now with weights, post alignment transform
asimple2u <- RItools:::alignDesignsByStrata("--", simple1,
post.align.transform = myrank)
expect_equivalent(colSums(asimple2u@Covariates ),
rep(0,ncol(asimple2u@Covariates)))
asimple2s <- RItools:::alignDesignsByStrata("strat", simple1,
post.align.transform = myrank)
tmp2 <- asimple2s@Covariates
expect_equivalent(colSums(tmp2[simple1@StrataFrame[["strat"]]=="a",]),
rep(0,ncol(asimple2s@Covariates)))
expect_equivalent(as.matrix(t(asimple2s@StrataMatrix) %*% tmp2),
matrix(0,2, ncol(asimple2s@Covariates)))
} )
test_that("scale() method wrapping to alignDesignsByStrata()",{
# at first pass, we're testing form but not content here.
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=rep(c(NA, 1), c(3,2)),
x2=c(1:5),
fac=factor(c(rep(1:2,2), NA))
)
dat$'(weights)' <- 1:5
simple2 <- RItools:::makeDesigns(z ~ x1 + x2 + fac+ strata(strat) + cluster(clus), data = dat)
scl2_scaleF <- scale(simple2, center=TRUE, scale=FALSE)
simple2b <- as(simple2, "StratumWeightedDesignOptions")
simple2b@Sweights <- RItools:::DesignWeights(simple2b, # need stratum weights to be present, even if ignored
RItools:::effectOfTreatmentOnTreated)
asimple2u <- RItools:::alignDesignsByStrata("--", simple2b,
post.align.transform = NULL)
expect_identical(scl2_scaleF, asimple2u@Covariates)
asimple2s <- RItools:::alignDesignsByStrata("strat", simple2b,
post.align.transform = NULL)
simple2c <- RItools:::makeDesigns(z ~ x1 + x2 + fac+ strata(strat) + cluster(clus) - 1, data = dat)
scl2c_scaleF <- scale(simple2c, center=TRUE, scale=FALSE)
expect_identical(scl2c_scaleF, asimple2s@Covariates)
scl2_scaleF_centerF <- scale(simple2, center=FALSE, scale=FALSE) # if it's a logical,
expect_identical(scl2_scaleF, scl2_scaleF_centerF) # `center` param is ignored
scl2_scaleT <- scale(simple2, center=TRUE, scale=TRUE)
expect_equal(length(dim(scl2_scaleT)), 2L)
expect_equivalent(is.na(scl2_scaleT),
matrix(FALSE, nrow(scl2_scaleT), ncol(scl2_scaleT)))
myrank <- function(x, weights) rank(x)
scl2_scaleF_centerrank <- scale(simple2, center = myrank, scale=FALSE)
expect_identical(dim(scl2_scaleF_centerrank), dim(scl2_scaleF))
expect_false(isTRUE(all.equal(scl2_scaleF, scl2_scaleF_centerrank, check.attributes=FALSE)),
"post alignment transform ignored")
wcent <- function(x, w) {x[order(x * w)]}
expect_silent(scl2_scaleF_wcent <-
scale(simple2, center=wcent, scale=FALSE)
)
expect_identical(dim(scl2_scaleF_wcent), dim(scl2_scaleF))
})
test_that("Issue #89: Proper strata weights", {
set.seed(20180208)
replicate(nreps_,{
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- 0.5 + 0.25 * x1 - 0.25 * x2 + rnorm(n)
idx <- 0.25 + 0.1 * x1 + 0.2 * x2 - 0.5 * x3 + rnorm(n)
y <- sample(rep(c(1,0), n/2), prob = exp(idx) / (1 + exp(idx)))
xy <- data.frame(x1, x2, x3, idx, y)
xy$m[y == 1] <- order(idx[y == 1])
xy$m[y == 0] <- order(idx[y == 0])
xy$"(weights)" <- 1
xy.wts <- xy
xy.wts$"(weights)" <- (1 + exp(idx)) / exp(idx) # inverse propensity score weights
design.nowts <- RItools:::makeDesigns(y ~ x1 + x2 + x3 + strata(m), data = xy)
design.wts <- RItools:::makeDesigns(y ~ x1 + x2 + x3 + strata(m), data = xy.wts)
## ETT weights are determined by assignment probabilities, not element counts
## or cluster masses. Accordingly presence/absence of unit weights shouldn't matter
## for their sweights. (But since they do affect h_b * m-bar_b, the corresponding
## wtratio's will be affected.)
ett.nowts <- RItools:::DesignWeights(design.nowts, stratum.weights=effectOfTreatmentOnTreated)
ett.wts <- RItools:::DesignWeights(design.wts, stratum.weights=effectOfTreatmentOnTreated)
expect_equal(ett.wts$m$sweights, ett.nowts$m$sweights)
## With a single stratum, sweights has to be 1, since it's normalized.
## wtratio is its ratio with h_b * m-bar_b, thus will generally be much
## less than 1. Check this:
h <- with(xy, 1/(1/sum(y) + 1/sum(!y)))
expect_equal(ett.nowts[['--']][,'wtratio'], 1/h)
h <- with(xy.wts, 1/(1/sum(y) + 1/sum(!y)))
expect_equal(ett.wts[['--']][,'wtratio'], 1/(h*mean(xy.wts$"(weights)")))
## split up into strata, use harmonic strata weights
## again unit weights shouldn't enter into this, although they
## would affect harmonic_times_mean_weight
dw.nowts <- RItools:::DesignWeights(design.nowts, stratum.weights=harmonic)
dw.wts <- RItools:::DesignWeights(design.wts, stratum.weights=harmonic)
expect_equal(dw.wts$m$sweights, dw.nowts$m$sweights)
## in this example by-stratum harmonic mean cluster counts are always 1 --
expect_equivalent(as.vector(dw.wts$m$sweights),
rep(1, nlevels(survival::strata(xy.wts$m)))/
nlevels(survival::strata(xy.wts$m))
)
## -- so we can check the calculation of the mean cluster mass factor as
## follows.
dw.wts2 <- RItools:::DesignWeights(design.wts)
clus_mean_weights <- tapply(xy.wts$"(weights)", xy.wts$m, mean)
expect_equivalent(dw.wts2$m$sweights, clus_mean_weights/sum(clus_mean_weights))
})
})
test_that("In inferentials, NAs imputed to stratum means",{
dat <- data.frame(z=rep(0:1, 2), x=c(0, NA, 2, 10),
x_median_imputed=c(0,2,2,10),
x_grandmean_imputed=c(0,4,2,10),
m=factor(rep(c("a", "b"), each=2)),
clus=factor(c("a", "b", "c", "b")),
w=1)
mf <- model.frame(z~x+x_median_imputed+x_grandmean_imputed+
m+clus, dat, weights=w, na.action=na.pass)
simple1 <- RItools:::makeDesigns(z ~ x+x_median_imputed+x_grandmean_imputed +
strata(m), data = mf)
simple1 <- as(simple1, "StratumWeightedDesignOptions")
simple1@Sweights <- RItools:::DesignWeights(simple1)
asimple1 <- sapply(colnames(simple1@StrataFrame),
RItools:::alignDesignsByStrata, design=simple1,
simplify=FALSE, USE.NAMES=TRUE)
expect_equivalent(asimple1[['--']]@Covariates[2,"x"],
mean(asimple1[['--']]@Covariates[c(1,3,4),"x"])
)
expect_lt(var(asimple1[['--']]@Covariates[,"x"]),
var(asimple1[['--']]@Covariates[,"x_median_imputed"])
)
expect_equivalent(var(asimple1[['--']]@Covariates[,"x"]),
var(asimple1[['--']]@Covariates[,"x_grandmean_imputed"])
)
expect_lt(var(asimple1[['m']]@Covariates[,"x"]),
var(asimple1[['m']]@Covariates[,"x_median_imputed"])
)
expect_lt(var(asimple1[['m']]@Covariates[,"x"]),
var(asimple1[['m']]@Covariates[,"x_grandmean_imputed"])
)
expect_equivalent(asimple1[['m']]@Covariates[2,"x"],
asimple1[['m']]@Covariates[1,"x"])
simple2 <- RItools:::makeDesigns(z ~ x+x_median_imputed+x_grandmean_imputed + cluster(clus), data = mf)
simple2 <- aggregateDesigns(simple2)
simple2 <- as(simple2, "StratumWeightedDesignOptions")
simple2@Sweights <- RItools:::DesignWeights(simple2)
asimple2 <- alignDesignsByStrata("--", design=simple2)
expect_equivalent(var(asimple2@Covariates[,"x"]),
var(asimple2@Covariates[,"x_grandmean_imputed"])
)
})
context("HB08*")
test_that("HB08 agreement w/ xBal()", {
set.seed(20180605)
replicate(nreps_,{
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- 0.5 + 0.25 * x1 - 0.25 * x2 + rnorm(n)
idx <- 0.25 + 0.1 * x1 + 0.2 * x2 - 0.5 * x3 + rnorm(n)
y <- sample(rep(c(1,0), n/2), prob = exp(idx) / (1 + exp(idx)))
xy <- data.frame(x1, x2, x3, idx, y)
xy$m[y == 1] <- order(idx[y == 1])
xy$m[y == 0] <- order(idx[y == 0])
## this mimics matched pairs:
expect_true(all(table(xy$y, xy$m)==1))
xy$'(weights)' <- rep(1L, n)
## first unweighted case
simple0 <- RItools:::makeDesigns(y ~ x1 + x2 + x3 + strata(m), data = xy)
simple0 <- as(simple0, "StratumWeightedDesignOptions")
simple0@Sweights <- RItools:::DesignWeights(simple0) # this test
asimple0 <- sapply(colnames(simple0@StrataFrame), RItools:::alignDesignsByStrata,
design=simple0, simplify=FALSE, USE.NAMES=TRUE)
btis0 <- lapply(asimple0, HB08)
xb0 <- xBalance(y ~ x1 + x2 + x3, data = xy,
strata = list(unmatched = NULL, matched = ~ m), report = 'all')
expect_equivalent(btis0[['--']]$adj.diff.of.totals[-4], # remove '(_non-null record_)' entry
xb0$results[,'adj.diff',"unmatched"])
expect_equivalent(btis0[['--']]$tcov[1:3,1:3], # remove '(_non-null record_)' entries
attr(xb0$overall, 'tcov')$unmatched)
expect_equivalent(btis0[['--']][c('Msq', 'DF')],
xb0[['overall']]["unmatched",c('chisquare', 'df'), drop=TRUE])
expect_equivalent(btis0[['m']]$adj.diff.of.totals[-4], # remove '(_non-null record_)' entry
xb0$results[,'adj.diff',"matched"])
expect_equivalent(btis0[['m']]$tcov[1:3,1:3], # remove '(_non-null record_)' entries
attr(xb0$overall, 'tcov')$matched)
expect_equivalent(btis0[['m']][c('Msq', 'DF')],
xb0[['overall']]["matched",c('chisquare', 'df'), drop=TRUE])
## now with weights.
## first, an example with weights that don't vary within the strata().
xy_wted1 <- xy; mwts <- 0
while (any(mwts==0)) mwts <- rpois(n/2, lambda=10)
xy_wted1$'(weights)' <- unsplit(mwts, xy$m)
simple1 <- RItools:::makeDesigns(y ~ x1 + x2 + x3 + strata(m), data = xy_wted1)
simple1 <- as(simple1, "StratumWeightedDesignOptions")
simple1@Sweights <- RItools:::DesignWeights(simple1) # this test
asimple1 <- sapply(colnames(simple1@StrataFrame),
RItools:::alignDesignsByStrata, design=simple1,
simplify=FALSE, USE.NAMES=TRUE)
btis1 <- lapply(asimple1, HB08)
wts.scaled <- xy_wted1$'(weights)' / mean(xy_wted1$'(weights)')
xy_xbwts1 <- transform(xy_wted1, x1=x1*wts.scaled,
x2=x2*wts.scaled, x3=x3*wts.scaled,
w=wts.scaled)
xb1u <- xBalance(y ~ x1 + x2 + x3 + w, data = xy_xbwts1,
strata = list(unmatched = NULL), report = 'all')
expect_equivalent(btis1[['--']]$adj.diff.of.totals,
xb1u$results[,'adj.diff',"unmatched"])
expect_equivalent(btis1[['--']]$tcov,
attr(xb1u$overall, 'tcov')$unmatched)
expect_equivalent(btis1[['--']][c('Msq', 'DF')],
xb1u[['overall']]["unmatched",c('chisquare', 'df'), drop=TRUE])
wts.scaled <- xy_wted1$'(weights)' / mean( mwts )
xy_xbwts1 <- transform(xy_wted1, x1=x1*wts.scaled,
x2=x2*wts.scaled, x3=x3*wts.scaled)
xb1m <- xBalance(y ~ x1 + x2 + x3, data = xy_xbwts1,
strata = list(matched = ~ m), report = 'all')
expect_equivalent(btis1[['m']]$adj.diff.of.totals[-4], # remove '(_non-null record_)' entry
xb1m$results[,'adj.diff',"matched"])
expect_equivalent(btis1[['m']]$tcov[1:3,1:3], # remove '(_non-null record_)' entries
attr(xb1m$overall, 'tcov')$matched)
expect_equivalent(btis1[['m']][c('Msq', 'DF')],
xb1m[['overall']]["matched",c('chisquare', 'df'), drop=TRUE])
## second and example where weights vary arbitrarily (over positive numbers)
xy_wted2 <- xy; wts <- 0
while (any(wts==0)) wts <- rpois(n, lambda=10)
xy_wted2$'(weights)' <- wts
simple2 <- RItools:::makeDesigns(y ~ x1 + x2 + x3 + strata(m), data = xy_wted2)
simple2 <- as(simple2, "StratumWeightedDesignOptions")
simple2@Sweights <- RItools:::DesignWeights(simple2) # this test
asimple2 <- sapply(colnames(simple2@StrataFrame),
RItools:::alignDesignsByStrata, design=simple2,
simplify=FALSE, USE.NAMES=TRUE)
btis1 <- lapply(asimple2, HB08)
wts.scaled <- xy_wted2$'(weights)' / mean(xy_wted2$'(weights)')
xy_xbwts2 <- transform(xy_wted2, x1=x1*wts.scaled,
x2=x2*wts.scaled, x3=x3*wts.scaled,
w=wts.scaled)
xb1u <- xBalance(y ~ x1 + x2 + x3 + w, data = xy_xbwts2,
strata = list(unmatched = NULL), report = 'all')
expect_equivalent(btis1[['--']]$adj.diff.of.totals,
xb1u$results[,'adj.diff',"unmatched"])
expect_equivalent(btis1[['--']]$tcov,
attr(xb1u$overall, 'tcov')$unmatched)
expect_equivalent(btis1[['--']][c('Msq', 'DF')],
xb1u[['overall']]["unmatched",c('chisquare', 'df'), drop=TRUE])
wt2.scaled <- xy_wted2$'(weights)' /
mean( tapply(xy_wted2$'(weights)', xy_wted2$m, mean) )
xy_xbwts2 <- transform(xy_wted2, x1=x1*wts.scaled,
x2=x2*wts.scaled, x3=x3*wts.scaled,
w=wts.scaled)
xb1m <- xBalance(y ~ x1 + x2 + x3 + w, data = xy_xbwts2,
strata = list(matched = ~ m), report = 'all')
expect_equivalent(btis1[['m']]$adj.diff.of.totals,
xb1m$results[,'adj.diff',"matched"])
expect_equivalent(btis1[['m']]$tcov,
attr(xb1m$overall, 'tcov')$matched)
expect_equivalent(btis1[['m']][c('Msq', 'DF')],
xb1m[['overall']]["matched",c('chisquare', 'df'), drop=TRUE])
})
} )
test_that("HB08_2016 agreement w/ xBal()", {
set.seed(20180605)
replicate(nreps_,{
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- 0.5 + 0.25 * x1 - 0.25 * x2 + rnorm(n)
idx <- 0.25 + 0.1 * x1 + 0.2 * x2 - 0.5 * x3 + rnorm(n)
y <- sample(rep(c(1,0), n/2), prob = exp(idx) / (1 + exp(idx)))
xy <- data.frame(x1, x2, x3, idx, y)
xy$m[y == 1] <- order(idx[y == 1])
xy$m[y == 0] <- order(idx[y == 0])
## this mimics matched pairs:
expect_true(all(table(xy$y, xy$m)==1))
xy$'(weights)' <- rep(1L, n)
## first unweighted case
simple0 <- RItools:::makeDesigns(y ~ x1 + x2 + x3 + strata(m), data = xy)
simple0 <- as(simple0, "StratumWeightedDesignOptions")
simple0@Sweights <- RItools:::DesignWeights(simple0) # this test
asimple0 <- sapply(colnames(simple0@StrataFrame),
RItools:::alignDesignsByStrata, simple0,
simplify=TRUE, USE.NAMES=TRUE)
btis0 <- lapply(asimple0, HB08_2016)
xb0 <- xBalance(y ~ x1 + x2 + x3, data = xy,
strata = list(unmatched = NULL, matched = ~ m), report = 'all')
expect_equivalent(btis0[['--']]$adj.diff.of.totals[-4], # remove '(_any Xs recorded_)' entry
xb0$results[,'adj.diff',"unmatched"])
expect_equivalent(btis0[['--']]$tcov[1:3,1:3], # remove '(_any Xs recorded_)' entries
attr(xb0$overall, 'tcov')$unmatched)
expect_equivalent(btis0[['--']][c('Msq', 'DF')],
xb0[['overall']]["unmatched",c('chisquare', 'df'), drop=TRUE])
expect_equivalent(btis0[['m']]$adj.diff.of.totals[-4], # remove '(_any Xs recorded_)' entry
xb0$results[,'adj.diff',"matched"])
expect_equivalent(btis0[['m']]$tcov[1:3,1:3], # remove '(_any Xs recorded_)' entries
attr(xb0$overall, 'tcov')$matched)
expect_equivalent(btis0[['m']][c('Msq', 'DF')],
xb0[['overall']]["matched",c('chisquare', 'df'), drop=TRUE])
## now with weights. Comparing adjusted diffs based on totals will only work
## if the weights don't vary with by stratum, at least in stratified case.
xy_wted <- xy; mwts <- 0
while (any(mwts==0)) mwts <- rpois(n/2, lambda=10)
## centering of variables is needed for unstratified mean diffs comparison.
xy_wted <- transform(xy_wted, x1=x1-weighted.mean(x1,unsplit(mwts, xy$m)),
x2=x2-weighted.mean(x2,unsplit(mwts, xy$m)),
x3=x3-weighted.mean(x3,unsplit(mwts, xy$m)))
xy_wted$'(weights)' <- unsplit(mwts, xy$m)
simple1 <- RItools:::makeDesigns(y ~ x1 + x2 + x3 + strata(m), data = xy_wted)
simple1 <- as(simple1, "StratumWeightedDesignOptions")
simple1@Sweights <- RItools:::DesignWeights(simple1) # this test
asimple1 <- sapply(colnames(simple1@StrataFrame),
RItools:::alignDesignsByStrata, simple1,
simplify=TRUE, USE.NAMES=TRUE)
btis1 <- lapply(asimple1, HB08_2016)
wts.scaled <- xy_wted$'(weights)' / mean(xy_wted$'(weights)')
xy_xbwts <- transform(xy_wted, x1=x1*wts.scaled,
x2=x2*wts.scaled, x3=x3*wts.scaled,
w=wts.scaled)
xb1u <- xBalance(y ~ x1 + x2 + x3 + w, data = xy_xbwts,
strata = list(unmatched = NULL), report = 'all')
expect_equivalent(btis1[['--']]$adj.diff.of.totals, # remove '(_any Xs recorded_)' entry
xb1u$results[,'adj.diff',"unmatched"])
expect_equivalent(btis1[['--']]$tcov, # remove '(_any Xs recorded_)' entries
attr(xb1u$overall, 'tcov')$unmatched)
expect_equivalent(btis1[['--']][c('Msq', 'DF')],
xb1u[['overall']]["unmatched",c('chisquare', 'df'), drop=TRUE])
wts.scaled <- xy_wted$'(weights)' / mean( mwts )
xy_xbwts <- transform(xy_wted, x1=x1*wts.scaled,
x2=x2*wts.scaled, x3=x3*wts.scaled,
w= wts.scaled)
xb1m <- xBalance(y ~ x1 + x2 + x3 + w, data = xy_xbwts,
strata = list(matched = ~ m), report = 'all')
expect_equivalent(btis1[['m']]$adj.diff.of.totals, # remove '(_any Xs recorded_)' entry
xb1m$results[1:3,'adj.diff',"matched"])
expect_equivalent(btis1[['m']]$tcov, # remove '(_any Xs recorded_)' entries
attr(xb1m$overall, 'tcov')$matched)
expect_equivalent(btis1[['m']][c('Msq', 'DF')],
xb1m[['overall']]["matched",c('chisquare', 'df'), drop=TRUE])
})
} )
test_that("HB08 and HB08_2016 flag degenerate statistics", {
set.seed(0303022134)
# pairs
s <- 20
n <- 2 * s
z <- rep(c(0,1), s)
b <- rep(1:s, each = 2)
## theory indicates that the statistic will be degenerate when r - (n - s) = 0
## where r is the rank of the matrix in the quadratic form of the test statistic
x <- replicate(n - s, runif(n, 0, 100))
colnames(x) <- paste0("x", 1:(n-s))
df_bad <- data.frame(z = z, x = x, b = b, '(weights)' = 1, check.names = FALSE)
expect_warning(balanceTest(z ~ . + strata(b), data = df_bad, inferentials.calculator = RItools:::HB08), "degenerate")
expect_warning(balanceTest(z ~ . + strata(b), data = df_bad, inferentials.calculator = RItools:::HB08_2016), "degenerate")
})
markmfredrickson/RItools documentation built on Jan. 15, 2024, 10:53 p.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.
################################################################################
# Tests for DesignOptions objects
################################################################################
if (!exists('nreps_')) nreps_ <- 10L
library("testthat")
context("ModelMatrixPlus S4 class carries per-covariate missingness info")
test_that("If no missing data, then NotMissing is a matrix w n rows and 0 cols",{
d <- data.frame(
id = 1:500,
x = rnorm(500),
cluster = rep(1:100, 5),
strata.good = rep(c(1:4, NA), 100),
strata.bad = sample(1:100, size = 500, replace = T),
z.good = rep(c(0,1), 250),
z.bad = sample(c(0,1), size = 500, replace = T))
d$'(weights)' = 1 # meet expectation of a weights column
simple <- RItools:::model_matrix(z.good ~ x, data = d)
expect_equivalent(dim(simple@NotMissing), c(500,1))
})
test_that("Missingness gets passed through in Covariates, recorded in NotMissing",{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=rep(c(NA, 1), c(3,2)),
x2=c(1:5),
fac=factor(c(rep(1:2,2), NA))
)
datmf <- model.frame(z ~ x1 + x2 + fac, dat, na.action = na.pass)
datmf$'(weights)' <- 1
simple2 <- RItools:::model_matrix(z ~ x1 + x2 + fac, data = datmf)
expect_equivalent(ncol(simple2@NotMissing), 3)
expect_equivalent(colnames(simple2@NotMissing), c("_any Xs recorded_", "x1", "fac"))
})
test_that("model_matrix offered missing or negative weights",{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
w_with_NAs=rep(c(NA, 1), c(3,2)),
w_with_negatives=-2:2,
x=c(1:5),
fac=factor(c(rep(1:2,2), NA))
)
datmf_with_NAs <-
model.frame(z ~ x + fac, dat, na.action = na.pass,
weights=w_with_NAs)
datmf_with_negatives <-
model.frame(z ~ x + fac, dat, na.action = na.pass,
weights=w_with_negatives)
expect_error(RItools:::model_matrix(z ~ x + fac,
data = datmf_with_NAs),
"uweights")
expect_error(RItools:::model_matrix(z ~ x + fac,
data = datmf_with_negatives),
"uweights")
})
test_that("lookup tables OK, even w/ complex & multi-column terms",{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=rep(c(NA, 1), c(3,2)),
x2=c(1:5),
fac=factor(c(rep(1:2,2), NA))
)
datmf <- model.frame(z ~ x1 + x2 + fac, dat, na.action = na.pass)
datmf$'(weights)' <- 1
simple2 <- RItools:::model_matrix(z ~ x1 + x2 + fac, data=datmf)
expect_equal(simple2@OriginalVariables, 1:3)
expect_equal(simple2@TermLabels, c( "x1", "x2", "fac"))
expect_equal(simple2@NM.Covariates, c(2,0,3))
expect_equal(simple2@NM.terms, c(2,0,3))
## check that complex term don't spell trouble in themselves
datmf <- model.frame(z ~ x1 + cut(x2, c(0,3,6)) + fac, data = dat,
na.action = na.pass)
datmf$'(weights)' <- 1
simple3 <- RItools:::model_matrix(z ~ x1 + cut(x2, c(0,3,6)) + fac, data = datmf)
expect_equal(simple3@OriginalVariables, 1:3)
expect_equal(simple3@TermLabels, c("x1", "cut(x2, c(0, 3, 6))", "fac"))
expect_equal(simple3@NM.Covariates, c(2,0,3))
expect_equal(simple3@NM.terms, c(2,0,3))
## now try a complex term that actually expands to multiple columns
datmf <- model.frame(z ~ x1 + cut(x2, c(0,3,6)) + fac, data = dat,
na.action = na.pass)
datmf$'(weights)' <- 1
simple4 <- RItools:::model_matrix(z ~ x1 + cut(x2, c(0,3,6)) + fac, data = datmf,
contrasts=list("cut(x2, c(0, 3, 6))"=diag(2)))
expect_equal(simple4@OriginalVariables, c(1,2,2,3))
expect_equal(simple4@TermLabels, c("x1", "cut(x2, c(0, 3, 6))", "fac"))
expect_equal(simple4@NM.Covariates, c(2,0,0,3))
expect_equal(simple4@NM.terms, c(2,0,3))
## now put NAs in the multi-column complex term
datmf <- model.frame(z ~ x2 + cut(x1, c(0,3,6)) + fac, data = dat,
na.action = na.pass)
datmf$'(weights)' <- 1
simple5 <- RItools:::model_matrix(z ~ x2 + cut(x1, c(0,3,6)) + fac, data = datmf,
contrasts=list("cut(x1, c(0, 3, 6))"=diag(2)))
expect_equal(simple5@OriginalVariables, c(1,2,2,3))
expect_equal(simple5@TermLabels, c("x2", "cut(x1, c(0, 3, 6))", "fac"))
expect_equal(simple5@NM.Covariates, c(0,2,2,3))
expect_equal(simple5@NM.terms, c(0,2,3))
})
test_that("Issue #76: Using I() in formulas", {
x <- data.frame(z=c(1,1)) # have to exclude
while (all(x$z==x[1L,'z'])) # degenerate case
x <- data.frame(x = rnorm(10), y = rnorm(10), z = rbinom(10, size = 1, p = 1/3))
x$"(weights)" <- 1
d <- makeDesigns(z ~ I(x * sin(y)), data = x)
expect_s4_class(d, "DesignOptions")
## While we're at it, confirm that the non-stratification
## encoded in this DesignOptions bears the column name "--".
})
test_that("Null stratification is encoded by '--'",{
data(nuclearplants)
nuclearplants$"(weights)" <- 1
d0 <- makeDesigns(pr ~ cost, data=nuclearplants)
expect_equal(colnames(d0@StrataFrame), "--")
foo <- nuclearplants$pt
d1 <- makeDesigns(pr ~ cost + strata(foo), data=nuclearplants)
expect_true("--" %in% colnames(d1@StrataFrame))
})
test_that("Issue #86: makeDesigns finds variables outside data arg",{
data(nuclearplants)
foo <- nuclearplants$pt
nuclearplants$"(weights)" <- 1
d <- makeDesigns(pr ~ cost + strata(foo), data=nuclearplants)
expect_s4_class(d, "DesignOptions")
})
test_that("Duplicated missingness patterns handled appropriately",{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=c(rep(NA, 3), 1:2),
x2=c(1:5),
fac=factor(c(rep(1:2,2), NA))
)
datmf <- model.frame(z ~ x1 + I(x1^2) + fac, data=dat, na.action = na.pass)
datmf$'(weights)' <- 1
simple2 <- RItools:::model_matrix(z ~ x1 + I(x1^2) + fac, data = datmf)
expect_equal(simple2@OriginalVariables, 1:3)
expect_equal(simple2@TermLabels, c("x1", "I(x1^2)", "fac"))
expect_equal(simple2@NM.Covariates, 1+c(1,1,2))
expect_equal(simple2@NM.terms, 1+c(1,1,2))
## If exactly two terms have missing data but in the same pattern, then
## NotMissing is a matrix w/ n rows and 1 col.
datmf <- model.frame(z ~ x1 + I(x1^2), data=dat, na.action = na.pass)
datmf$'(weights)' <- 1
simple3 <- RItools:::model_matrix(z ~ x1 + I(x1^2), data = datmf)
expect_equal(simple3@OriginalVariables, 1:2)
expect_equal(simple3@TermLabels, c("x1", "I(x1^2)"))
expect_equal(ncol(simple3@NotMissing), 1)
expect_equal(simple3@NM.Covariates, c(1,1))
expect_equal(simple3@NM.terms, c(1,1))
})
test_that("All-Xes missingness |-> NotMissing col '_any Xs recorded_'",{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=c(rep(NA, 4), 2),
x2=c(1:3, NA, 5),
fac=factor(c(NA, 1:2, NA, 2))
)
dat$'(weights)' <- 1
datmf <- model.frame(z ~ x1 + I(x1^2) + fac, data=dat, na.action = na.pass)
datmf$'(weights)' <- 1
simple6 <- RItools:::model_matrix(z ~ x1 + I(x1^2) + fac, data = datmf)
expect_equal(simple6@OriginalVariables, 1:3)
expect_equal(simple6@TermLabels, c("x1", "I(x1^2)", "fac"))
expect_equal(colnames(simple6@NotMissing)[1], "_any Xs recorded_")
expect_equal(simple6@NM.Covariates, c(2,2,1))
expect_equal(simple6@NM.terms, c(2,2,1))
simple7 <- makeDesigns(z ~ x1 + I(x1^2) + fac + 0 + strata(strat) + cluster(clus),
dat)
simple7 <- as(simple7, "StratumWeightedDesignOptions")
simple7@Sweights <-
RItools:::DesignWeights(simple7,
RItools:::effectOfTreatmentOnTreated)
## As writing of this test, DesignWeights() expects only pre-aggregated designs,
## and infers treatment:control ratios from the numbers of elements in in each
## of condition (by strata), not the number of clusters. Thus in this example
## it should believe that the ETT weights are proportional to 2 for stratum a,
## 1 for stratum b:
expect_equivalent(simple7@Sweights$strat$sweights, (2:1)/3)
## key point: had missingness of each of unit 4's covariates tricked it into ignoring
## that observation, then the ETT weight for stratum b would have been proportional to 0,
## not 1.
} )
test_that("All-Xes missingness noted in descriptives", {
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=c(rep(NA, 4), 2),
x2=c(1:3, NA, 5),
fac=factor(c(NA, 1:2, NA, 2))
)
dat$'(weights)' <- 1
simple5 <- makeDesigns(z ~ x1 + I(x1^2) + fac + cluster(clus), dat)
descr <- RItools:::designToDescriptives(simple5)
expect_equivalent(descr['(_any Xs recorded_)' ,1:2,1], c(1, 1/3)) # not just 1s all around
})
context("DesignOptions objects")
test_that("Creating DesignOptions objects", {
set.seed(20130801)
replicate(nreps_,{
d <- data.frame(
id = 1:500,
x = rnorm(500),
cluster = rep(1:100, 5),
strata.good = rep(c(1:4, NA), 100),
strata.bad = sample(1:100, size = 500, replace = T),
z.good = rep(c(0,1), 250),
z.bad = sample(c(0,1), size = 500, replace = T))
d$'(weights)' = 1 # meet expectation of a weights column
# checking input
expect_error(RItools:::makeDesigns(~ x, data = d), "treatment")
expect_error(RItools:::makeDesigns(z.bad ~ x + strata(strata.good) + cluster(cluster), data = d), "cluster")
expect_error(RItools:::makeDesigns(z.good ~ x + strata(strata.bad) + cluster(cluster), data = d), "strata")
expect_error(RItools:::makeDesigns(z.good ~ x + cluster(id) + cluster(cluster), data = d), "cluster")
expect_error(RItools:::makeDesigns(z.good ~ x - 1, data = d, "stratification"))
# actually testing that the output is as expected
simple <- RItools:::makeDesigns(z.good ~ x, data = d)
expect_equal(dim(simple@StrataFrame)[2], 1)
expect_equivalent(simple@Covariates[, "x"], d$x)
expect_equivalent(simple@Z, as.logical(d$z.good))
expect_equal(nlevels(simple@Cluster), 500) # a cluster per individual
clustered <- RItools:::makeDesigns(z.good ~ x + cluster(cluster), data = d)
expect_equal(dim(clustered@StrataFrame)[2], 1)
expect_true(nlevels(clustered@Cluster) > 1)
clustStrata <- RItools:::makeDesigns(z.good ~ x + cluster(cluster) + strata(strata.good), data = d)
expect_equal(dim(clustStrata@StrataFrame)[2], 2)
clustStrata.c <- RItools:::makeDesigns(z.good ~ x + cluster(cluster) + strata(strata.good, strata.good), data = d)
expect_equivalent(clustStrata, clustStrata.c)
# dropping the overall comparison
expect_equal(dim(RItools:::makeDesigns(z.good ~ x + cluster(cluster) + strata(strata.good) - 1, data = d)@StrataFrame)[2], 1)
## More tests to write:
# - All NA strata variables
# - Missing z or cluster
# - strata with extra levels but no observations (which can be safely dropped)
# (NB: extra levels tested upstream, in xBalance, as of commit 34861515;
# see ./test.clusters.R )
})
})
test_that("NotMissing vars correctly generated",
{
replicate(nreps_,{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=rep(c(NA, TRUE), c(3,2)),
x2 = c(1:5),
x3=c(TRUE, FALSE, NA, TRUE, FALSE),
fac=factor(c(rep(1:2,2), NA))
)
dat$'(weights)' <- 1
simple <- RItools:::makeDesigns(z ~ x1 + x2 + strata(strat) + cluster(clus), data = dat)
expect_match(colnames(simple@NotMissing), "x1", all=FALSE)
expect_false(any(grepl("x2", colnames(simple@NotMissing))))
expect_false(any(grepl("TRUE", colnames(simple@Covariates))))
expect_false(any(grepl("FALSE", colnames(simple@Covariates))))
simple2 <- RItools:::makeDesigns(z ~ x1 + x2 + fac+ strata(strat) + cluster(clus), data = dat)
expect_match(colnames(simple2@NotMissing), "x1", all=FALSE)
expect_false(any(grepl("x2", colnames(simple2@NotMissing))))
expect_match(colnames(simple2@NotMissing), "fac", all=FALSE)
expect_false(any(grepl("TRUE", colnames(simple2@Covariates))))
expect_false(any(grepl("FALSE", colnames(simple2@Covariates))))
simple3 <- RItools:::makeDesigns(z ~ x1 + x3 + fac+ strata(strat) + cluster(clus), data = dat)
expect_match(colnames(simple3@NotMissing), "x1", all=FALSE)
expect_match(colnames(simple3@NotMissing), "x3", all=FALSE)
expect_match(colnames(simple3@NotMissing), "fac", all=FALSE)
expect_false(any(grepl("TRUE", colnames(simple3@Covariates))))
expect_false(any(grepl("FALSE", colnames(simple3@Covariates))))
})
})
test_that("Issue 88: logical Covariates correctly generated",
{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=rep(c(NA, TRUE), c(3,2)),
x2 = c(1:5),
x3=c(TRUE, FALSE, NA, TRUE, FALSE),
fac=factor(c(rep(1:2,2), NA))
)
dat$'(weights)' <- 1
simple1 <- RItools:::makeDesigns(z ~ x1 + x2, data = dat)
expect_false(any(grepl("TRUE", colnames(simple1@Covariates))))
expect_false(any(grepl("FALSE", colnames(simple1@Covariates))))
simple2 <- RItools:::makeDesigns(z ~ x1 + x2 + strata(strat), data = dat)
expect_false(any(grepl("TRUE", colnames(simple2@Covariates))))
expect_false(any(grepl("FALSE", colnames(simple2@Covariates))))
simple3 <- RItools:::makeDesigns(z ~ x1 + x2 + strata(strat) - 1, data = dat)
expect_false(any(grepl("TRUE", colnames(simple3@Covariates))))
expect_false(any(grepl("FALSE", colnames(simple3@Covariates))))
})
test_that("DesignOptions to descriptive statistics", {
set.seed(20130801)
replicate(nreps_,{
d <- data.frame(
x = rnorm(500),
f = factor(sample(c("A", "B", "C"), size = 500, replace = T)),
c = rep(1:100, 5),
s = rep(c(1:4, NA), 100),
z = rep(c(0,1), 250))
d$'(weights)' = 1 # meet expectation of a weights column
simple <- RItools:::makeDesigns(z ~ x + f + strata(s) + cluster(c), data = d)
descriptives <- RItools:::designToDescriptives(simple)
expect_is(descriptives, "array")
expect_equal(dim(descriptives), c(5, 5, 2))
# descriptives ignore clustering
design.noclus <- RItools:::makeDesigns(z ~ x + f + strata(s), data = d)
expect_equal(descriptives, RItools:::designToDescriptives(design.noclus))
# the strata should imply different stats
expect_false(identical(descriptives[,,1], descriptives[,,2]))
# however, the pooled s.d.s should be the same.
expect_identical(descriptives[,"pooled.sd","--"], descriptives[,"pooled.sd","s"])
# ok, now checking that values are correct.
expect_equal(mean(d$x[d$z == 1]), descriptives["x", "Treatment", "--"])
expect_equal(mean(d$x[d$z == 0]), descriptives["x", "Control", "--"])
# with equal sized strata, the the control/treatment means are the means of the the strata means
expect_equal(mean(tapply(d$x[d$z == 1], d$s[d$z == 1], mean)), descriptives["x", "Treatment", "s"])
expect_equal(mean(tapply(d$x[d$z == 0], d$s[d$z == 0], mean)), descriptives["x", "Control", "s"])
})
})
test_that("designToDescriptives uses provided covariate scales",{
d <- data.frame(x=rnorm(50), z=rep(c(0,1), 25))
d$'(weights)' <- 1
simple <- RItools:::makeDesigns(z ~ x, data=d)
sd_x <- sd(resid(lm(x ~z, data=d)))
descriptives <- designToDescriptives(simple, covariate.scales=c(x=sd_x*10))
expect_equal(descriptives["x","pooled.sd","--"], sd_x*10)
descriptives <- designToDescriptives(simple,
covariate.scales=c(x=sd_x*10, y=Inf))
expect_equal(descriptives["x","pooled.sd","--"], sd_x*10)
expect_warning(designToDescriptives(simple, covariate.scales=sd_x),
"name")
expect_warning(designToDescriptives(simple, covariate.scales=c(x="foo")),
"numeric")
})
test_that("descriptives for NotMissing variables",
{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=rep(c(NA, 1), c(3,2)),
x2=c(1:5),
fac=factor(c(rep(1:2,2), NA))
)
dat$'(weights)' <- 1
simple <- RItools:::makeDesigns(z ~ x1 + x2 + strata(strat) + cluster(clus), data = dat)
expect_false(any(grepl("NA", colnames(simple@Covariates))))
dsimple <- RItools:::designToDescriptives(simple)
expect_match(dimnames(dsimple)[[1]], "(x1)", all=FALSE)
expect_false(any(grepl("(x2)", dimnames(dsimple)[[1]], fixed=TRUE)))
simple2 <- RItools:::makeDesigns(z ~ x1 + x2 + fac+ strata(strat) + cluster(clus), data = dat)
expect_false(any(grepl("NA", colnames(simple2@Covariates))))
dsimple2 <- RItools:::designToDescriptives(simple2)
expect_match(dimnames(dsimple2)[[1]], "(x1)", all=FALSE)
expect_match(dimnames(dsimple2)[[1]], "(fac)", all=FALSE)
expect_false(any(grepl("(x2)", dimnames(dsimple2)[[1]], fixed=TRUE)))
}
)
test_that("Issue 36: Descriptives with NAs, appropriate weighting", {
### Descriptives with missing covariates ###
set.seed(20130801)
replicate(nreps_,{
d <- data.frame(
x = rnorm(500),
f = factor(sample(c("A", "B", "C"), size = 500, replace = T)),
c = rep(1:100, 5),
s = rep(c(1:4, NA), 100),
paired = rep(c(0,1), each = 250),
z = rep(c(0,1), 250))
d$'(weights)' = 1
d.missing <- d
d.missing$x[d.missing$x < -1] <- NA
simple.all <- RItools:::makeDesigns(z ~ x + f + strata(s) + cluster(c), data = d)
descriptives.all <- RItools:::designToDescriptives(simple.all)
expect_equal(descriptives.all["x", "Treatment", "--"], mean(d$x[d$z == 1]))
expect_equal(descriptives.all["x", "Treatment", "s"], mean(d$x[d$z == 1 & !is.na(d$s)]))
simple.missing <- RItools:::makeDesigns(z ~ x + f + strata(s) + cluster(c), data = d.missing)
descriptives.missing <- RItools:::designToDescriptives(simple.missing)
with(d.missing,
expect_equal(descriptives.missing["x", "Treatment", "--"],
mean(x[z == 1], na.rm = TRUE)))
with(d.missing,
expect_equal(descriptives.missing["x", "Treatment", "s"],
mean(x[z == 1 & !is.na(s)], na.rm = TRUE)))
expect_false(identical(descriptives.all, descriptives.missing))
# ETT weighting
design.paired <- RItools:::makeDesigns(z ~ x + f + strata(paired) + strata(s), data = d)
descriptives.paired <- RItools:::designToDescriptives(design.paired)
with(d, expect_equal(descriptives.paired["x", "Control", "paired"], mean(x[z == 0])))
with(d, expect_equal(descriptives.paired["x", "Control", "--"], mean(x[z == 0])))
with(d, expect_false(identical(descriptives.paired["x", "Control", "s"], mean(x[z == 0]))))
})
})
test_that("Aggegating designs by clusters", {
set.seed(20130801)
replicate(nreps_,{
d <- data.frame(
x = rnorm(500),
f = factor(sample(c("A", "B", "C"), size = 500, replace = T)),
c = rep(1:100, 5),
s = rep(c(1:4, NA), 100),
z = rep(c(0,1), 250))
d$'(weights)' <- 1
# grab a bunch of rows and duplicate them
d <- rbind(d, d[sample(1:dim(d)[1], size = 100), ])
# swapping around the data to make sure order doesn't mask bugs
d <- d[sample(1:dim(d)[1]), ]
design <- RItools:::makeDesigns(z ~ x + f + strata(s) + cluster(c), data = d)
aggDesign <- RItools:::aggregateDesigns(design)
# one row per cluster, with columns x, fa, fb, fc
expect_equal(dim(aggDesign@Covariates), c(100, 4))
# now spot check some cluster totals of totals
expect_equal(aggDesign@Covariates[1, ], colMeans(design@Covariates[design@Cluster == 1,]))
# Z's roll up as they should
Zs <- tapply(design@Z, design@Cluster, mean)
dim(Zs) <- NULL
Zs <- as.logical(Zs)
expect_equivalent(aggDesign@Z, Zs)
# extraneous levels in the Cluster slot are ignored
design2 <- design
levels(design2@Cluster) <- c(levels(design2@Cluster), letters)
aggDesign2 <- RItools:::aggregateDesigns(design2)
expect_equal(dim(aggDesign2@Covariates), c(100, 4))
})
})
test_that("aggregateDesigns treats NA covariates as 0's" ,{
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x=rep(c(NA, 1), c(3,2))
)
dat$'(weights)' <- 1
design <- RItools:::makeDesigns(z~x+strata(strat)+cluster(clus)-1, dat)
aggDesign <- RItools:::aggregateDesigns(design)
expect_equivalent(aggDesign@Covariates[,'x'], c(0, 0, 1, 1) )
nm.column.for.x <- aggDesign@NM.Covariates[match( 'x', colnames(aggDesign@Covariates))]
expect_equal(aggDesign@NotMissing[,nm.column.for.x], c(0,0,1,1) )
})
test_that("Aggregation of unit weights to cluster level",{
set.seed(20130801)
replicate(nreps_,{
d.short <- data.frame(
x = rnorm(500),
f = factor(sample(c("A", "B", "C"), size = 500, replace = T)),
c = rep(1:100, 5),
s = rep(c(1:4, NA), 100),
z = rep(c(0,1), 250))
## grab a bunch of rows and duplicate them
newrows <- c(1:nrow(d.short), sample(1:nrow(d.short), size=100, replace=T) )
d.tall <- d.short[newrows,]
d.tall$'(weights)' <- 1
d.short$'(weights)' <- as.vector(table(newrows))
design.tall <- RItools:::makeDesigns(z ~ x + f + strata(s) + cluster(c), data = d.tall)
aggDesign.tall <- RItools:::aggregateDesigns(design.tall)
design.short <- RItools:::makeDesigns(z ~ x + f + strata(s) + cluster(c), data = d.short)
aggDesign.short <- RItools:::aggregateDesigns(design.short)
## we should wind up in the same place.
expect_equal(aggDesign.tall, aggDesign.short)
d2 <- d.tall
d2$'(weights)' <- 2
design.d2 <- RItools:::makeDesigns(z ~ x + f + strata(s) + cluster(c), data = d2)
aggDesign.d2 <- RItools:::aggregateDesigns(design.d2)
expect_equal(2*aggDesign.tall@UnitWeights,aggDesign.d2@UnitWeights)
expect_equal(aggDesign.tall@NotMissing, aggDesign.d2@NotMissing)
expect_equal(aggDesign.tall@Covariates, aggDesign.d2@Covariates)
})
})
context("ModelMatrixPlus S4 class enriches model.matrix-value 'class' ")
test_that("R core hasn't revised conventions we may depend on",
{
ff <- log(Volume) ~ log(Height) + log(Girth)
m <- model.frame(ff, trees)
mm <- model.matrix(ff, m)
expect_equal("(Intercept)", colnames(mm)[1])
expect_equivalent(mm[,1,drop=TRUE], rep(1, nrow(mm)))
expect_false("(Intercept)" %in% colnames(model.matrix(update(ff, .~.-1), m)))
expect_equal(dim(mm[,integer(0)]), c(nrow(mm), 0))
expect_equal(complete.cases(mm[,integer(0)]), rep(TRUE, nrow(mm)))
expect_equivalent(weighted.mean(c(1:4, NA), c(rep(1,4), 0)), 2.5)
})
test_that("Model matrix material is properly formed",
{
ff <- log(Volume) ~ log(Height) + log(Girth)
trees1 <- trees
trees1$'(weights)' <- 1
DM0 <- model_matrix(ff, trees1, remove.intercept=FALSE)
expect_is(DM0, "ModelMatrixPlus")
m <- model.frame(ff, trees)
expect_equivalent(model.matrix(ff, m), as.matrix(DM0))
fff <- update(ff, .~.-1)
trees2 <- trees
trees2[1, "Volume"] <- NA # LHS variable, shouldn't affect anything
m2a <- model.frame(fff, trees2, na.action = na.pass)
m2a$'(weights)' <- 1
expect_equivalent(model.matrix(fff, m), as.matrix(model_matrix(fff, m2a)))
trees2[1, "Height"] <- NA # RHS variable, but still shouldn't cause rows to be dropped
m2b <- model.frame(fff, trees2, na.action = na.pass)
m2b$'(weights)' <- 1
expect_equal(dim(model.matrix(fff, m)),
dim(as.matrix(model_matrix(fff,m2b) ) )
)
## specified contrasts
dd <- data.frame(a = gl(3,4), b = gl(4,1,12)) # balanced 2-way
dd$'(weights)' <- 1
expect_equal(model.matrix(~ a + b, dd),
as.matrix(model_matrix(~ a + b, dd, remove.intercept=FALSE)))
expect_equal(model.matrix(~ a + b-1, dd),
as.matrix(model_matrix(~ a + b-1, dd)))
expect_equal(model.matrix(~ a + b, dd, contrasts = list(a = "contr.sum")),
as.matrix(model_matrix(~ a + b, dd,
contrasts = list(a = "contr.sum"), remove.intercept=FALSE))
)
expect_equal(model.matrix(~ a + b, dd,
contrasts = list(a = "contr.sum", b = "contr.poly")),
as.matrix(model_matrix(~ a + b, dd,
contrasts =
list(a = "contr.sum", b = "contr.poly"),
remove.intercept=FALSE))
)
}
)
context("alignDesignsByStrata")
test_that("alignDesigns, designToDescriptives output alignment", {
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=rep(c(NA, 1), c(3,2)),
x2=c(1:5),
fac=factor(c(rep(1:2,2), NA))
)
dat$'(weights)' <- 1
simple2 <- RItools:::makeDesigns(z ~ x1 + x2 + fac+ strata(strat) + cluster(clus), data = dat)
expect_equal(colnames(simple2@StrataFrame), c("strat", "--"))
simple2 <- as(simple2, "StratumWeightedDesignOptions")
simple2@Sweights <- RItools:::DesignWeights(simple2, # Have to aggregate 1st to figure stratum weights
RItools:::effectOfTreatmentOnTreated)
expect_true(setequal(names(simple2@Sweights), c("strat", "--")))
dsimple2 <- RItools:::designToDescriptives(simple2)
asimple2a <- RItools:::alignDesignsByStrata("--", simple2)
expect_equivalent(dimnames(dsimple2)[[1]], colnames(asimple2a@Covariates))
asimple2b <- RItools:::alignDesignsByStrata("strat", simple2)
expect_equivalent(dimnames(dsimple2)[[1]], colnames(asimple2b@Covariates))
})
test_that("alignDesigns centers covars by stratum", {
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=rep(c(NA, 1), c(3,2)),
x2=c(1:5),
fac=factor(c(rep(1:2,2), NA))
)
dat$'(weights)' <- 1
## first unweighted case
simple0 <- RItools:::makeDesigns(z ~ x1 + x2 + fac+ strata(strat) + cluster(clus), data = dat)
simple0 <- as(simple0, "StratumWeightedDesignOptions")
simple0@Sweights <- RItools:::DesignWeights(simple0, # Placeholder strat weights, shouldn't affect
RItools:::effectOfTreatmentOnTreated) # this test
asimple0u <- RItools:::alignDesignsByStrata("--",simple0)
expect_equivalent(colSums(asimple0u@Covariates),
rep(0,ncol(asimple0u@Covariates)))
asimple0s <- RItools:::alignDesignsByStrata("strat",simple0)
expect_equivalent(colSums(asimple0s@Covariates[simple0@StrataFrame[["strat"]]=="a",]),
rep(0,ncol(asimple0s@Covariates)))
expect_equivalent(as.matrix(t(asimple0s@StrataMatrix) %*% asimple0s@Covariates),
matrix(0,2,ncol(asimple0s@Covariates)))
## now with weights
dat1 <- dat
dat1$'(weights)' <- rpois(nrow(dat1), lambda=10)
while (any(dat1$'(weights)'==0)) dat1$'(weights)' <- rpois(nrow(dat1), lambda=10)
simple1 <- RItools:::makeDesigns(z ~ x1 + x2 + fac+ strata(strat) + cluster(clus), data = dat1)
simple1 <- as(simple1, "StratumWeightedDesignOptions")
simple1@Sweights <- RItools:::DesignWeights(simple1, # Placeholder strat weights, shouldn't affect
RItools:::effectOfTreatmentOnTreated) # this test
asimple1u <- RItools:::alignDesignsByStrata("--",simple1)
expect_equivalent(colSums(asimple1u@Covariates),
rep(0,ncol(asimple1u@Covariates)))
asimple1s <- RItools:::alignDesignsByStrata("strat",simple1)
tmp1 <- asimple1s@Covariates
expect_equivalent(colSums(tmp1[simple1@StrataFrame[["strat"]]=="a",]),
rep(0,ncol(asimple1s@Covariates)))
expect_equivalent(as.matrix(t(asimple1s@StrataMatrix) %*% tmp1),
matrix(0,2, ncol(asimple1s@Covariates)))
myrank <- function(x, weights) rank(x)
## now with weights, post alignment transform
asimple2u <- RItools:::alignDesignsByStrata("--", simple1,
post.align.transform = myrank)
expect_equivalent(colSums(asimple2u@Covariates ),
rep(0,ncol(asimple2u@Covariates)))
asimple2s <- RItools:::alignDesignsByStrata("strat", simple1,
post.align.transform = myrank)
tmp2 <- asimple2s@Covariates
expect_equivalent(colSums(tmp2[simple1@StrataFrame[["strat"]]=="a",]),
rep(0,ncol(asimple2s@Covariates)))
expect_equivalent(as.matrix(t(asimple2s@StrataMatrix) %*% tmp2),
matrix(0,2, ncol(asimple2s@Covariates)))
} )
test_that("scale() method wrapping to alignDesignsByStrata()",{
# at first pass, we're testing form but not content here.
dat <- data.frame(strat=rep(letters[1:2], c(3,2)),
clus=factor(c(1,1,2:4)),
z=c(TRUE, rep(c(TRUE, FALSE), 2)),
x1=rep(c(NA, 1), c(3,2)),
x2=c(1:5),
fac=factor(c(rep(1:2,2), NA))
)
dat$'(weights)' <- 1:5
simple2 <- RItools:::makeDesigns(z ~ x1 + x2 + fac+ strata(strat) + cluster(clus), data = dat)
scl2_scaleF <- scale(simple2, center=TRUE, scale=FALSE)
simple2b <- as(simple2, "StratumWeightedDesignOptions")
simple2b@Sweights <- RItools:::DesignWeights(simple2b, # need stratum weights to be present, even if ignored
RItools:::effectOfTreatmentOnTreated)
asimple2u <- RItools:::alignDesignsByStrata("--", simple2b,
post.align.transform = NULL)
expect_identical(scl2_scaleF, asimple2u@Covariates)
asimple2s <- RItools:::alignDesignsByStrata("strat", simple2b,
post.align.transform = NULL)
simple2c <- RItools:::makeDesigns(z ~ x1 + x2 + fac+ strata(strat) + cluster(clus) - 1, data = dat)
scl2c_scaleF <- scale(simple2c, center=TRUE, scale=FALSE)
expect_identical(scl2c_scaleF, asimple2s@Covariates)
scl2_scaleF_centerF <- scale(simple2, center=FALSE, scale=FALSE) # if it's a logical,
expect_identical(scl2_scaleF, scl2_scaleF_centerF) # `center` param is ignored
scl2_scaleT <- scale(simple2, center=TRUE, scale=TRUE)
expect_equal(length(dim(scl2_scaleT)), 2L)
expect_equivalent(is.na(scl2_scaleT),
matrix(FALSE, nrow(scl2_scaleT), ncol(scl2_scaleT)))
myrank <- function(x, weights) rank(x)
scl2_scaleF_centerrank <- scale(simple2, center = myrank, scale=FALSE)
expect_identical(dim(scl2_scaleF_centerrank), dim(scl2_scaleF))
expect_false(isTRUE(all.equal(scl2_scaleF, scl2_scaleF_centerrank, check.attributes=FALSE)),
"post alignment transform ignored")
wcent <- function(x, w) {x[order(x * w)]}
expect_silent(scl2_scaleF_wcent <-
scale(simple2, center=wcent, scale=FALSE)
)
expect_identical(dim(scl2_scaleF_wcent), dim(scl2_scaleF))
})
test_that("Issue #89: Proper strata weights", {
set.seed(20180208)
replicate(nreps_,{
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- 0.5 + 0.25 * x1 - 0.25 * x2 + rnorm(n)
idx <- 0.25 + 0.1 * x1 + 0.2 * x2 - 0.5 * x3 + rnorm(n)
y <- sample(rep(c(1,0), n/2), prob = exp(idx) / (1 + exp(idx)))
xy <- data.frame(x1, x2, x3, idx, y)
xy$m[y == 1] <- order(idx[y == 1])
xy$m[y == 0] <- order(idx[y == 0])
xy$"(weights)" <- 1
xy.wts <- xy
xy.wts$"(weights)" <- (1 + exp(idx)) / exp(idx) # inverse propensity score weights
design.nowts <- RItools:::makeDesigns(y ~ x1 + x2 + x3 + strata(m), data = xy)
design.wts <- RItools:::makeDesigns(y ~ x1 + x2 + x3 + strata(m), data = xy.wts)
## ETT weights are determined by assignment probabilities, not element counts
## or cluster masses. Accordingly presence/absence of unit weights shouldn't matter
## for their sweights. (But since they do affect h_b * m-bar_b, the corresponding
## wtratio's will be affected.)
ett.nowts <- RItools:::DesignWeights(design.nowts, stratum.weights=effectOfTreatmentOnTreated)
ett.wts <- RItools:::DesignWeights(design.wts, stratum.weights=effectOfTreatmentOnTreated)
expect_equal(ett.wts$m$sweights, ett.nowts$m$sweights)
## With a single stratum, sweights has to be 1, since it's normalized.
## wtratio is its ratio with h_b * m-bar_b, thus will generally be much
## less than 1. Check this:
h <- with(xy, 1/(1/sum(y) + 1/sum(!y)))
expect_equal(ett.nowts[['--']][,'wtratio'], 1/h)
h <- with(xy.wts, 1/(1/sum(y) + 1/sum(!y)))
expect_equal(ett.wts[['--']][,'wtratio'], 1/(h*mean(xy.wts$"(weights)")))
## split up into strata, use harmonic strata weights
## again unit weights shouldn't enter into this, although they
## would affect harmonic_times_mean_weight
dw.nowts <- RItools:::DesignWeights(design.nowts, stratum.weights=harmonic)
dw.wts <- RItools:::DesignWeights(design.wts, stratum.weights=harmonic)
expect_equal(dw.wts$m$sweights, dw.nowts$m$sweights)
## in this example by-stratum harmonic mean cluster counts are always 1 --
expect_equivalent(as.vector(dw.wts$m$sweights),
rep(1, nlevels(survival::strata(xy.wts$m)))/
nlevels(survival::strata(xy.wts$m))
)
## -- so we can check the calculation of the mean cluster mass factor as
## follows.
dw.wts2 <- RItools:::DesignWeights(design.wts)
clus_mean_weights <- tapply(xy.wts$"(weights)", xy.wts$m, mean)
expect_equivalent(dw.wts2$m$sweights, clus_mean_weights/sum(clus_mean_weights))
})
})
test_that("In inferentials, NAs imputed to stratum means",{
dat <- data.frame(z=rep(0:1, 2), x=c(0, NA, 2, 10),
x_median_imputed=c(0,2,2,10),
x_grandmean_imputed=c(0,4,2,10),
m=factor(rep(c("a", "b"), each=2)),
clus=factor(c("a", "b", "c", "b")),
w=1)
mf <- model.frame(z~x+x_median_imputed+x_grandmean_imputed+
m+clus, dat, weights=w, na.action=na.pass)
simple1 <- RItools:::makeDesigns(z ~ x+x_median_imputed+x_grandmean_imputed +
strata(m), data = mf)
simple1 <- as(simple1, "StratumWeightedDesignOptions")
simple1@Sweights <- RItools:::DesignWeights(simple1)
asimple1 <- sapply(colnames(simple1@StrataFrame),
RItools:::alignDesignsByStrata, design=simple1,
simplify=FALSE, USE.NAMES=TRUE)
expect_equivalent(asimple1[['--']]@Covariates[2,"x"],
mean(asimple1[['--']]@Covariates[c(1,3,4),"x"])
)
expect_lt(var(asimple1[['--']]@Covariates[,"x"]),
var(asimple1[['--']]@Covariates[,"x_median_imputed"])
)
expect_equivalent(var(asimple1[['--']]@Covariates[,"x"]),
var(asimple1[['--']]@Covariates[,"x_grandmean_imputed"])
)
expect_lt(var(asimple1[['m']]@Covariates[,"x"]),
var(asimple1[['m']]@Covariates[,"x_median_imputed"])
)
expect_lt(var(asimple1[['m']]@Covariates[,"x"]),
var(asimple1[['m']]@Covariates[,"x_grandmean_imputed"])
)
expect_equivalent(asimple1[['m']]@Covariates[2,"x"],
asimple1[['m']]@Covariates[1,"x"])
simple2 <- RItools:::makeDesigns(z ~ x+x_median_imputed+x_grandmean_imputed + cluster(clus), data = mf)
simple2 <- aggregateDesigns(simple2)
simple2 <- as(simple2, "StratumWeightedDesignOptions")
simple2@Sweights <- RItools:::DesignWeights(simple2)
asimple2 <- alignDesignsByStrata("--", design=simple2)
expect_equivalent(var(asimple2@Covariates[,"x"]),
var(asimple2@Covariates[,"x_grandmean_imputed"])
)
})
context("HB08*")
test_that("HB08 agreement w/ xBal()", {
set.seed(20180605)
replicate(nreps_,{
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- 0.5 + 0.25 * x1 - 0.25 * x2 + rnorm(n)
idx <- 0.25 + 0.1 * x1 + 0.2 * x2 - 0.5 * x3 + rnorm(n)
y <- sample(rep(c(1,0), n/2), prob = exp(idx) / (1 + exp(idx)))
xy <- data.frame(x1, x2, x3, idx, y)
xy$m[y == 1] <- order(idx[y == 1])
xy$m[y == 0] <- order(idx[y == 0])
## this mimics matched pairs:
expect_true(all(table(xy$y, xy$m)==1))
xy$'(weights)' <- rep(1L, n)
## first unweighted case
simple0 <- RItools:::makeDesigns(y ~ x1 + x2 + x3 + strata(m), data = xy)
simple0 <- as(simple0, "StratumWeightedDesignOptions")
simple0@Sweights <- RItools:::DesignWeights(simple0) # this test
asimple0 <- sapply(colnames(simple0@StrataFrame), RItools:::alignDesignsByStrata,
design=simple0, simplify=FALSE, USE.NAMES=TRUE)
btis0 <- lapply(asimple0, HB08)
xb0 <- xBalance(y ~ x1 + x2 + x3, data = xy,
strata = list(unmatched = NULL, matched = ~ m), report = 'all')
expect_equivalent(btis0[['--']]$adj.diff.of.totals[-4], # remove '(_non-null record_)' entry
xb0$results[,'adj.diff',"unmatched"])
expect_equivalent(btis0[['--']]$tcov[1:3,1:3], # remove '(_non-null record_)' entries
attr(xb0$overall, 'tcov')$unmatched)
expect_equivalent(btis0[['--']][c('Msq', 'DF')],
xb0[['overall']]["unmatched",c('chisquare', 'df'), drop=TRUE])
expect_equivalent(btis0[['m']]$adj.diff.of.totals[-4], # remove '(_non-null record_)' entry
xb0$results[,'adj.diff',"matched"])
expect_equivalent(btis0[['m']]$tcov[1:3,1:3], # remove '(_non-null record_)' entries
attr(xb0$overall, 'tcov')$matched)
expect_equivalent(btis0[['m']][c('Msq', 'DF')],
xb0[['overall']]["matched",c('chisquare', 'df'), drop=TRUE])
## now with weights.
## first, an example with weights that don't vary within the strata().
xy_wted1 <- xy; mwts <- 0
while (any(mwts==0)) mwts <- rpois(n/2, lambda=10)
xy_wted1$'(weights)' <- unsplit(mwts, xy$m)
simple1 <- RItools:::makeDesigns(y ~ x1 + x2 + x3 + strata(m), data = xy_wted1)
simple1 <- as(simple1, "StratumWeightedDesignOptions")
simple1@Sweights <- RItools:::DesignWeights(simple1) # this test
asimple1 <- sapply(colnames(simple1@StrataFrame),
RItools:::alignDesignsByStrata, design=simple1,
simplify=FALSE, USE.NAMES=TRUE)
btis1 <- lapply(asimple1, HB08)
wts.scaled <- xy_wted1$'(weights)' / mean(xy_wted1$'(weights)')
xy_xbwts1 <- transform(xy_wted1, x1=x1*wts.scaled,
x2=x2*wts.scaled, x3=x3*wts.scaled,
w=wts.scaled)
xb1u <- xBalance(y ~ x1 + x2 + x3 + w, data = xy_xbwts1,
strata = list(unmatched = NULL), report = 'all')
expect_equivalent(btis1[['--']]$adj.diff.of.totals,
xb1u$results[,'adj.diff',"unmatched"])
expect_equivalent(btis1[['--']]$tcov,
attr(xb1u$overall, 'tcov')$unmatched)
expect_equivalent(btis1[['--']][c('Msq', 'DF')],
xb1u[['overall']]["unmatched",c('chisquare', 'df'), drop=TRUE])
wts.scaled <- xy_wted1$'(weights)' / mean( mwts )
xy_xbwts1 <- transform(xy_wted1, x1=x1*wts.scaled,
x2=x2*wts.scaled, x3=x3*wts.scaled)
xb1m <- xBalance(y ~ x1 + x2 + x3, data = xy_xbwts1,
strata = list(matched = ~ m), report = 'all')
expect_equivalent(btis1[['m']]$adj.diff.of.totals[-4], # remove '(_non-null record_)' entry
xb1m$results[,'adj.diff',"matched"])
expect_equivalent(btis1[['m']]$tcov[1:3,1:3], # remove '(_non-null record_)' entries
attr(xb1m$overall, 'tcov')$matched)
expect_equivalent(btis1[['m']][c('Msq', 'DF')],
xb1m[['overall']]["matched",c('chisquare', 'df'), drop=TRUE])
## second and example where weights vary arbitrarily (over positive numbers)
xy_wted2 <- xy; wts <- 0
while (any(wts==0)) wts <- rpois(n, lambda=10)
xy_wted2$'(weights)' <- wts
simple2 <- RItools:::makeDesigns(y ~ x1 + x2 + x3 + strata(m), data = xy_wted2)
simple2 <- as(simple2, "StratumWeightedDesignOptions")
simple2@Sweights <- RItools:::DesignWeights(simple2) # this test
asimple2 <- sapply(colnames(simple2@StrataFrame),
RItools:::alignDesignsByStrata, design=simple2,
simplify=FALSE, USE.NAMES=TRUE)
btis1 <- lapply(asimple2, HB08)
wts.scaled <- xy_wted2$'(weights)' / mean(xy_wted2$'(weights)')
xy_xbwts2 <- transform(xy_wted2, x1=x1*wts.scaled,
x2=x2*wts.scaled, x3=x3*wts.scaled,
w=wts.scaled)
xb1u <- xBalance(y ~ x1 + x2 + x3 + w, data = xy_xbwts2,
strata = list(unmatched = NULL), report = 'all')
expect_equivalent(btis1[['--']]$adj.diff.of.totals,
xb1u$results[,'adj.diff',"unmatched"])
expect_equivalent(btis1[['--']]$tcov,
attr(xb1u$overall, 'tcov')$unmatched)
expect_equivalent(btis1[['--']][c('Msq', 'DF')],
xb1u[['overall']]["unmatched",c('chisquare', 'df'), drop=TRUE])
wt2.scaled <- xy_wted2$'(weights)' /
mean( tapply(xy_wted2$'(weights)', xy_wted2$m, mean) )
xy_xbwts2 <- transform(xy_wted2, x1=x1*wts.scaled,
x2=x2*wts.scaled, x3=x3*wts.scaled,
w=wts.scaled)
xb1m <- xBalance(y ~ x1 + x2 + x3 + w, data = xy_xbwts2,
strata = list(matched = ~ m), report = 'all')
expect_equivalent(btis1[['m']]$adj.diff.of.totals,
xb1m$results[,'adj.diff',"matched"])
expect_equivalent(btis1[['m']]$tcov,
attr(xb1m$overall, 'tcov')$matched)
expect_equivalent(btis1[['m']][c('Msq', 'DF')],
xb1m[['overall']]["matched",c('chisquare', 'df'), drop=TRUE])
})
} )
test_that("HB08_2016 agreement w/ xBal()", {
set.seed(20180605)
replicate(nreps_,{
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- 0.5 + 0.25 * x1 - 0.25 * x2 + rnorm(n)
idx <- 0.25 + 0.1 * x1 + 0.2 * x2 - 0.5 * x3 + rnorm(n)
y <- sample(rep(c(1,0), n/2), prob = exp(idx) / (1 + exp(idx)))
xy <- data.frame(x1, x2, x3, idx, y)
xy$m[y == 1] <- order(idx[y == 1])
xy$m[y == 0] <- order(idx[y == 0])
## this mimics matched pairs:
expect_true(all(table(xy$y, xy$m)==1))
xy$'(weights)' <- rep(1L, n)
## first unweighted case
simple0 <- RItools:::makeDesigns(y ~ x1 + x2 + x3 + strata(m), data = xy)
simple0 <- as(simple0, "StratumWeightedDesignOptions")
simple0@Sweights <- RItools:::DesignWeights(simple0) # this test
asimple0 <- sapply(colnames(simple0@StrataFrame),
RItools:::alignDesignsByStrata, simple0,
simplify=TRUE, USE.NAMES=TRUE)
btis0 <- lapply(asimple0, HB08_2016)
xb0 <- xBalance(y ~ x1 + x2 + x3, data = xy,
strata = list(unmatched = NULL, matched = ~ m), report = 'all')
expect_equivalent(btis0[['--']]$adj.diff.of.totals[-4], # remove '(_any Xs recorded_)' entry
xb0$results[,'adj.diff',"unmatched"])
expect_equivalent(btis0[['--']]$tcov[1:3,1:3], # remove '(_any Xs recorded_)' entries
attr(xb0$overall, 'tcov')$unmatched)
expect_equivalent(btis0[['--']][c('Msq', 'DF')],
xb0[['overall']]["unmatched",c('chisquare', 'df'), drop=TRUE])
expect_equivalent(btis0[['m']]$adj.diff.of.totals[-4], # remove '(_any Xs recorded_)' entry
xb0$results[,'adj.diff',"matched"])
expect_equivalent(btis0[['m']]$tcov[1:3,1:3], # remove '(_any Xs recorded_)' entries
attr(xb0$overall, 'tcov')$matched)
expect_equivalent(btis0[['m']][c('Msq', 'DF')],
xb0[['overall']]["matched",c('chisquare', 'df'), drop=TRUE])
## now with weights. Comparing adjusted diffs based on totals will only work
## if the weights don't vary with by stratum, at least in stratified case.
xy_wted <- xy; mwts <- 0
while (any(mwts==0)) mwts <- rpois(n/2, lambda=10)
## centering of variables is needed for unstratified mean diffs comparison.
xy_wted <- transform(xy_wted, x1=x1-weighted.mean(x1,unsplit(mwts, xy$m)),
x2=x2-weighted.mean(x2,unsplit(mwts, xy$m)),
x3=x3-weighted.mean(x3,unsplit(mwts, xy$m)))
xy_wted$'(weights)' <- unsplit(mwts, xy$m)
simple1 <- RItools:::makeDesigns(y ~ x1 + x2 + x3 + strata(m), data = xy_wted)
simple1 <- as(simple1, "StratumWeightedDesignOptions")
simple1@Sweights <- RItools:::DesignWeights(simple1) # this test
asimple1 <- sapply(colnames(simple1@StrataFrame),
RItools:::alignDesignsByStrata, simple1,
simplify=TRUE, USE.NAMES=TRUE)
btis1 <- lapply(asimple1, HB08_2016)
wts.scaled <- xy_wted$'(weights)' / mean(xy_wted$'(weights)')
xy_xbwts <- transform(xy_wted, x1=x1*wts.scaled,
x2=x2*wts.scaled, x3=x3*wts.scaled,
w=wts.scaled)
xb1u <- xBalance(y ~ x1 + x2 + x3 + w, data = xy_xbwts,
strata = list(unmatched = NULL), report = 'all')
expect_equivalent(btis1[['--']]$adj.diff.of.totals, # remove '(_any Xs recorded_)' entry
xb1u$results[,'adj.diff',"unmatched"])
expect_equivalent(btis1[['--']]$tcov, # remove '(_any Xs recorded_)' entries
attr(xb1u$overall, 'tcov')$unmatched)
expect_equivalent(btis1[['--']][c('Msq', 'DF')],
xb1u[['overall']]["unmatched",c('chisquare', 'df'), drop=TRUE])
wts.scaled <- xy_wted$'(weights)' / mean( mwts )
xy_xbwts <- transform(xy_wted, x1=x1*wts.scaled,
x2=x2*wts.scaled, x3=x3*wts.scaled,
w= wts.scaled)
xb1m <- xBalance(y ~ x1 + x2 + x3 + w, data = xy_xbwts,
strata = list(matched = ~ m), report = 'all')
expect_equivalent(btis1[['m']]$adj.diff.of.totals, # remove '(_any Xs recorded_)' entry
xb1m$results[1:3,'adj.diff',"matched"])
expect_equivalent(btis1[['m']]$tcov, # remove '(_any Xs recorded_)' entries
attr(xb1m$overall, 'tcov')$matched)
expect_equivalent(btis1[['m']][c('Msq', 'DF')],
xb1m[['overall']]["matched",c('chisquare', 'df'), drop=TRUE])
})
} )
test_that("HB08 and HB08_2016 flag degenerate statistics", {
set.seed(0303022134)
# pairs
s <- 20
n <- 2 * s
z <- rep(c(0,1), s)
b <- rep(1:s, each = 2)
## theory indicates that the statistic will be degenerate when r - (n - s) = 0
## where r is the rank of the matrix in the quadratic form of the test statistic
x <- replicate(n - s, runif(n, 0, 100))
colnames(x) <- paste0("x", 1:(n-s))
df_bad <- data.frame(z = z, x = x, b = b, '(weights)' = 1, check.names = FALSE)
expect_warning(balanceTest(z ~ . + strata(b), data = df_bad, inferentials.calculator = RItools:::HB08), "degenerate")
expect_warning(balanceTest(z ~ . + strata(b), data = df_bad, inferentials.calculator = RItools:::HB08_2016), "degenerate")
})
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.