Nothing
tests/testthat/test-BuyseTest-standardization.R
In BuyseTest: Generalized Pairwise Comparisons
### test-BuyseTest-standardization.R ---
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: apr 11 2025 (15:42)
## Version:
## Last-Updated: May 8 2025 (10:27)
## By: Brice Ozenne
## Update #: 105
##----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
##----------------------------------------------------------------------
##
### Code:
if(FALSE){
library(testthat)
library(BuyseTest)
library(data.table)
library(survival)
library(prodlim)
}
context("Check standardization of the net benefit across strata \n")
## * Setting
BuyseTest.options(order.Hprojection = 1,
method.inference = "u-statistic",
trace = 0)
## * Binary outcome and binary covariate
## ** simulate data
seed <- 1
n <- 100
tau <- 0
set.seed(seed)
A <- rbinom(n, size = 1, prob = 0.75)
X <- sample(0:1, size = n, replace = TRUE)
df <- data.frame(id = 1:n,
Y = 0.5 * A + X + rnorm(n),
arm = A, X = X)
## ** standardization (first order)
test_that("BuyseTest - standarization (first order, binary outcome)",{
## 1 - marginal
## BT.unstrat <- BuyseTest(data = df, arm ~ cont(Y, threshold = tau), trace = F)
## BT.unstrat2 <- BuyseTest(data = df, arm ~ cont(Y, threshold = tau), trace = TRUE, method.inference = "bootstrap")
## rbind(confint(BT.unstrat, transformation = FALSE),
## confint(BT.unstrat2, method.ci = "gaussian"))
## estimate se lower.ci upper.ci null p.value
## Y 0.1811263 0.1188404 -0.05179662 0.4140493 0 0.1274803
## Y1 0.1811263 0.1220947 -0.06875553 0.4096230 0 0.1543344
## 2- standardized
BT.std <- BuyseTest(data = df, arm ~ cont(Y, threshold = tau) + X, pool.strata = "standardization", trace = FALSE)
expect_equal(mean(getIid(BT.std)),0, tol = 1e-6)
expect_equivalent(confint(BT.std, transformation = FALSE),
data.frame("estimate" = c(0.23088881),"se" = c(0.10551968),"lower.ci" = c(0.02407403),"upper.ci" = c(0.43770359),"null" = c(0),"p.value" = c(0.02866149)),,
tol = 1e-6)
expect_equivalent(coef(BT.std, strata = TRUE), c("0" = 0.03658537, "1.0" = -0.30406504, "0.1" = 0.75, "1" = 0.49583333))
expect_error(confint(BT.std, strata = TRUE)) ## cannot compute strata specific uncertainty since the same individual appear in sevearl strata
BT.std2 <- BuyseTest(data = df, arm ~ cont(Y, threshold = tau) + X, pool.strata = "standardization", trace = FALSE,
method.inference = "bootstrap", strata.resampling = c("arm","X"), seed = 10)
expect_equivalent(confint(BT.std2, method.ci = "gaussian", transformation = FALSE),
data.frame("estimate" = c(0.23088881),"se" = c(0.10915709),"lower.ci" = c(0.01694484),"upper.ci" = c(0.44483278),"null" = c(0),"p.value" = c(0.03441312)),
tol = 1e-6)
## 3- stratified to stratified (by hand)
df00 <- df[df$X==0,]
df10 <- rbind(df[df$X==0 & df$arm==1,],df[df$X==1 & df$arm==0,])
df01 <- rbind(df[df$X==0 & df$arm==0,],df[df$X==1 & df$arm==1,])
df11 <- df[df$X==1,]
BT.strat00 <- BuyseTest(arm ~ cont(Y, threshold = tau), method.inference = "U-statistic", trace = F, data = df00)
BT.strat10 <- BuyseTest(arm ~ cont(Y, threshold = tau), method.inference = "U-statistic", trace = F, data = df10)
BT.strat01 <- BuyseTest(arm ~ cont(Y, threshold = tau), method.inference = "U-statistic", trace = F, data = df01)
BT.strat11 <- BuyseTest(arm ~ cont(Y, threshold = tau), method.inference = "U-statistic", trace = F, data = df11)
grid.strata <- expand.grid(X = sort(unique(df$X)), Xstar = sort(unique(df$X)))
grid.strata$X.X <- paste(grid.strata$X,grid.strata$Xstar,sep=".")
df.extended <- do.call(rbind,apply(grid.strata, MARGIN = 1, function(iRow){
cbind(rbind(df[df$arm == 0 & df$X == iRow[1],],
df[df$arm == 1 & df$X == iRow[2],]),
X.X = paste(iRow[1],iRow[2],sep="."))
}))
grid.strata$nC.X <- tapply(df.extended$arm==0,df.extended$X.X,sum)[grid.strata$X.X]
grid.strata$nE.Xstar <- tapply(df.extended$arm==1,df.extended$X.X,sum)[grid.strata$X.X]
grid.strata$n.X <- as.numeric(table(df$X)[as.character(grid.strata$X)])
grid.strata$n.Xstar <- as.numeric(table(df$X)[as.character(grid.strata$Xstar)])
grid.strata$w.unadjusted <- grid.strata$nC.X * grid.strata$nE.Xstar / prod(table(df$arm))
grid.strata$w.standardize <- grid.strata$n.X * grid.strata$n.Xstar / NROW(df$arm)^2
## direct
Myiid <- matrix(0, nrow = NROW(df), ncol = 4, dimnames = list(NULL,unique(df.extended$X.X)))
Myiid[df00$id,"0.0"] <- getIid(BT.strat00)
Myiid[df10$id,"1.0"] <- getIid(BT.strat10)
Myiid[df01$id,"0.1"] <- getIid(BT.strat01)
Myiid[df11$id,"1.1"] <- getIid(BT.strat11)
MyWiid <- sweep(Myiid, FUN = "*", MARGIN = 2, STATS = grid.strata$w.standardize)
## mean(rowSums(MyWiid))
expect_equivalent(BT.std@covariance[,"netBenefit"], crossprod(rowSums(MyWiid)), tol = 1e-6)
## from the partial win/loss sums
Myiid2 <- matrix(0, nrow = NROW(df), ncol = 4, dimnames = list(NULL,unique(df.extended$X.X)))
Myiid2[df00$id,"0.0"] <- getIid(BT.strat00, scale = FALSE, center = FALSE) * table(df00$arm)[as.character(1-df00$arm)]
Myiid2[df10$id,"1.0"] <- getIid(BT.strat10, scale = FALSE, center = FALSE) * table(df10$arm)[as.character(1-df10$arm)]
Myiid2[df01$id,"0.1"] <- getIid(BT.strat01, scale = FALSE, center = FALSE) * table(df01$arm)[as.character(1-df01$arm)]
Myiid2[df11$id,"1.1"] <- getIid(BT.strat11, scale = FALSE, center = FALSE) * table(df11$arm)[as.character(1-df11$arm)]
## MyWiid2 <- sweep(Myiid2, FUN = "*", MARGIN = 2, STATS = grid.strata$w.standardize)
## MyWiid2.norm <- matrix(0, nrow = NROW(df), ncol = 4, dimnames = list(NULL,unique(df.extended$X.X)))
## MyWiid2.norm[df00$id,"0.0"] <- ((MyWiid2[df00$id,"0.0"] / table(df00$arm)[as.character(1-df00$arm)]) - coef(BT.strat00)*BT.std@weightStrata[1])/table(df00$arm)[as.character(df00$arm)]
## MyWiid2.norm[df10$id,"1.0"] <- ((MyWiid2[df10$id,"1.0"] / table(df10$arm)[as.character(1-df10$arm)]) - coef(BT.strat10)*BT.std@weightStrata[2])/table(df10$arm)[as.character(df10$arm)]
## MyWiid2.norm[df01$id,"0.1"] <- ((MyWiid2[df01$id,"0.1"] / table(df01$arm)[as.character(1-df01$arm)]) - coef(BT.strat01)*BT.std@weightStrata[3])/table(df01$arm)[as.character(df01$arm)]
## MyWiid2.norm[df11$id,"1.1"] <- ((MyWiid2[df11$id,"1.1"] / table(df11$arm)[as.character(1-df11$arm)]) - coef(BT.strat11)*BT.std@weightStrata[4])/table(df11$arm)[as.character(df11$arm)]
## mean(rowSums(MyWiid2.norm))
## crossprod(rowSums(MyWiid2.norm))
## this is what is implemented in the C++ function
MyWiid3 <- rowSums(sweep(Myiid2, FUN = "*", MARGIN = 2, STATS = grid.strata$w.standardize/(grid.strata$w.unadjusted*prod(table(df$arm)))))
## grid.strata$w.unadjusted*prod(table(df$arm)) - c(prod(table(df00$arm)),prod(table(df10$arm)),prod(table(df01$arm)),prod(table(df11$arm)))
MyWiid3[df00$id] <- MyWiid3[df00$id] - coef(BT.strat00)*BT.std@weightStrata[1]/table(df00$arm)[as.character(df00$arm)]
MyWiid3[df10$id] <- MyWiid3[df10$id] - coef(BT.strat10)*BT.std@weightStrata[2]/table(df10$arm)[as.character(df10$arm)]
MyWiid3[df01$id] <- MyWiid3[df01$id] - coef(BT.strat01)*BT.std@weightStrata[3]/table(df01$arm)[as.character(df01$arm)]
MyWiid3[df11$id] <- MyWiid3[df11$id] - coef(BT.strat11)*BT.std@weightStrata[4]/table(df11$arm)[as.character(df11$arm)]
expect_equal(mean(MyWiid3), 0, tol = 1e-6)
expect_equivalent(BT.std@covariance[,"netBenefit"], crossprod(MyWiid3), tol = 1e-6)
})
## ** standardization (second order)
BuyseTest.options(order.Hprojection = 2)
test_that("BuyseTest - standarization (second order, binary outcome, balanced)",{
## with balanced dataset standardization equals marginal analysis
dfB <- do.call(rbind,by(df, interaction(df$arm,df$X), function(iDF){iDF[1:min(table(df$arm,df$X)),,drop=FALSE]}))
dfB$id <- 1:NROW(dfB)
rownames(dfB) <- NULL
BTB.marginal <- BuyseTest(data = dfB, arm ~ cont(Y, threshold = tau), trace = FALSE)
BTB.marginal2 <- BuyseTest(data = dfB, arm ~ cont(Y, threshold = tau), trace = F, keep.pairScore = TRUE)
BTB.std <- BuyseTest(data = dfB, arm ~ cont(Y, threshold = tau) + X, pool.strata = "standardization", trace = F, keep.pairScore = FALSE)
BTB.std2 <- BuyseTest(data = dfB, arm ~ cont(Y, threshold = tau) + X, pool.strata = "standardization", trace = F, keep.pairScore = TRUE)
expect_equal(coef(BTB.marginal), coef(BTB.std), tol = 1e-5)
expect_equal(coef(BTB.std), mean(coef(BTB.std, strata = TRUE)), tol = 1e-5)
expect_true(confint(BTB.marginal)$se > confint(BTB.marginal, order.Hprojection = 1)$se) ## confint(BTB.marginal)$se^2 - confint(BTB.marginal, order.Hprojection = 1)$se^2
expect_true(confint(BTB.std)$se > confint(BTB.std, order.Hprojection = 1)$se) ## confint(BTB.std)$se^2 - confint(BTB.std, order.Hprojection = 1)$se^2
expect_equivalent(confint(BTB.std), confint(BTB.std2), tol = 1e-6)
## by hand (marginal)
sigmaMB_h1fav <- tapply(getIid(BTB.marginal, statistic = "favorable"), dfB$arm, crossprod)
termM1 <- (coef(BTB.marginal, statistic = "favorable") - coef(BTB.marginal, statistic = "favorable")^2)/prod(table(dfB$arm)) - sum(sigmaMB_h1fav/rev(table(dfB$arm)))
sigmaMB_h1unfav <- tapply(getIid(BTB.marginal, statistic = "unfavorable"), dfB$arm, crossprod)
termM2 <- (coef(BTB.marginal, statistic = "unfavorable") - coef(BTB.marginal, statistic = "unfavorable")^2)/prod(table(dfB$arm)) - sum(sigmaMB_h1unfav/rev(table(dfB$arm)))
sigmaMB_h1cross <- tapply(getIid(BTB.marginal, statistic = "favorable")*getIid(BTB.marginal, statistic = "unfavorable"), dfB$arm, sum)
termM3 <- 0 - coef(BTB.marginal, statistic = "favorable")*coef(BTB.marginal, statistic = "unfavorable")/prod(table(dfB$arm)) - sum(sigmaMB_h1cross/rev(table(dfB$arm)))
expect_equal(termM1+termM2-2*termM3, confint(BTB.marginal, order.Hprojection = 2)$se^2 - confint(BTB.marginal, order.Hprojection = 1)$se^2, tol = 1e-6)
mytableM <- getPairScore(BTB.marginal2)
HprojM_fav <- getIid(BTB.marginal2, statistic = "favorable", scale = FALSE)
estM_fav <- coef(BTB.marginal2, statistic = "favorable")
mytableM[, hproj.0 := HprojM_fav[index.0]]
mytableM[, hproj.1 := HprojM_fav[index.1]]
expect_equal(termM1, mytableM[, sum((favorable - hproj.0 - hproj.1 - estM_fav)^2)/.N^2], tol = 1e-6)
## by hand (standardized)
sigmaB_h1fav <- tapply(getIid(BTB.std, statistic = "favorable"), dfB$arm, crossprod)
term1 <- (coef(BTB.std, statistic = "favorable") - coef(BTB.std, statistic = "favorable")^2)/prod(table(dfB$arm)) - sum(sigmaB_h1fav/rev(table(dfB$arm)))
sigmaB_h1unfav <- tapply(getIid(BTB.std, statistic = "unfavorable"), dfB$arm, crossprod)
term2 <- (coef(BTB.std, statistic = "unfavorable") - coef(BTB.std, statistic = "unfavorable")^2)/prod(table(dfB$arm)) - sum(sigmaB_h1unfav/rev(table(dfB$arm)))
sigmaB_h1cross <- tapply(getIid(BTB.std, statistic = "favorable")*getIid(BTB.std, statistic = "unfavorable"), dfB$arm, sum)
term3 <- 0 - coef(BTB.std, statistic = "favorable")*coef(BTB.std, statistic = "unfavorable")/prod(table(dfB$arm)) - sum(sigmaB_h1cross/rev(table(dfB$arm)))
expect_equal(term1+term2-2*term3, confint(BTB.std, order.Hprojection = 2)$se^2 - confint(BTB.std, order.Hprojection = 1)$se^2, tol = 1e-6)
mytable <- getPairScore(BTB.std2)
mytable[, weightStrata := stats::setNames(BTB.std2@weightStrata,BTB.std2@level.strata)[strata]]
mytable[, m.X := length(unique(index.0)), by = "strata"]
mytable[, n.X := length(unique(index.1)), by = "strata"]
mytable[, m := length(attr(BTB.std2@level.treatment,"indexC"))]
mytable[, n := length(attr(BTB.std2@level.treatment,"indexT"))]
mytable[, Wfavorable := favorable * weightStrata/((m.X*n.X)/(n*m))]
## mean(mytable$Wfavorable)
## coef(BTB.std2, statistic = "favorable")
Hproj_fav <- getIid(BTB.std2, statistic = "favorable", scale = FALSE)
mytable[, hproj.0 := Hproj_fav[index.0]]
mytable[, hproj.1 := Hproj_fav[index.1]]
mytable[, Tfavorable := coef(BTB.std2, statistic = "favorable")]
expect_equal(term1, mytable[, sum((Wfavorable - hproj.0 - hproj.1 - Tfavorable)^2)/.N^2], tol = 1e-6)
})
test_that("BuyseTest - standarization (second order, binary outcome, unbalanced)",{
## 2- standardized
BTUB.std <- BuyseTest(data = df, arm ~ cont(Y, threshold = tau) + X, pool.strata = "standardization", trace = FALSE)
BTUB.std2 <- BuyseTest(data = df, arm ~ cont(Y, threshold = tau) + X, pool.strata = "standardization", keep.pairScore = TRUE, trace = FALSE)
expect_equivalent(confint(BTUB.std), confint(BTUB.std2), tol = 1e-6)
expect_equivalent(confint(BTUB.std), data.frame("estimate" = c(0.23088881),"se" = c(0.10675249),"lower.ci" = c(0.01411401),"upper.ci" = c(0.42693406),"null" = c(0),"p.value" = c(0.03705691)), tol = 1e-6)
## by hand (standardized)
mytableUB <- getPairScore(BTUB.std)
mytableUB[, weightStrata := stats::setNames(BTUB.std@weightStrata,BTUB.std@level.strata)[strata]]
mytableUB[, m.X := length(unique(index.0)), by = "strata"]
mytableUB[, n.X := length(unique(index.1)), by = "strata"]
mytableUB[, m := length(attr(BTUB.std@level.treatment,"indexC"))]
mytableUB[, n := length(attr(BTUB.std@level.treatment,"indexT"))]
mytableUB[, Wfavorable := favorable * weightStrata/((m.X*n.X)/(n*m))]
## mean(mytableUB$Wfavorable)
## coef(BTUB.std, statistic = "favorable")
Hproj_fav <- getIid(BTUB.std2, statistic = "favorable", scale = FALSE)
mytableUB[, hproj.0 := Hproj_fav[index.0]]
mytableUB[, hproj.1 := Hproj_fav[index.1]]
mytableUB[, Tfavorable := coef(BTUB.std2, statistic = "favorable")]
term1 <- confint(BTUB.std, statistic = "favorable", order.Hprojection = 2)$se^2 - confint(BTUB.std, statistic = "favorable", order.Hprojection = 1)$se^2
expect_equal(term1, mytableUB[, sum((Wfavorable - hproj.0 - hproj.1 - Tfavorable)^2)/.N^2], tol = 1e-6)
})
BuyseTest.options(order.Hprojection = 1)
## * check Peron scoring rule in a perfectly balanced dataset
n <- 50
set.seed(10)
dt.strata <- simBuyseTest(n,
n.strata = 2,
names.strata = "stratum")
dtB.strata <- do.call(rbind,by(dt.strata, interaction(dt.strata$treatment,dt.strata$stratum), function(iDT){iDT[1:min(table(dt.strata$treatment,dt.strata$stratum))]}))
rownames(dtB.strata) <- NULL
test_that("BuyseTest - standarization (tte outcome)",{
BT.raw <- BuyseTest(treatment ~ tte(eventtime, status), trace = FALSE, data = dt.strata)
confint(BT.raw)
## estimate se lower.ci upper.ci null p.value
## eventtime -0.008226969 0.1286523 -0.2546679 0.2392174 0 0.9490145
BT.std <- BuyseTest(treatment ~ tte(eventtime, status) + stratum, trace = FALSE, data = dt.strata,
pool.strata = "standardization")
## estimate se lower.ci upper.ci null p.value
## eventtime -0.009499187 0.127931 -0.2545409 0.2366887 0 0.9408131
expect_equivalent(confint(BT.std), data.frame("estimate" = c(-0.00949919),"se" = c(0.12793103),"lower.ci" = c(-0.25454086),"upper.ci" = c(0.23668869),"null" = c(0),"p.value" = c(0.9408131)), tol = 1e-5)
BTB.raw <- BuyseTest(treatment ~ tte(eventtime, status), trace = FALSE, data = dtB.strata,
model.tte = prodlim(Hist(eventtime, status) ~ treatment, data = dtB.strata))
confint(BTB.raw)
## estimate se lower.ci upper.ci null p.value
## eventtime 0.02881503 0.1417496 -0.2441965 0.2975942 0 0.8390032
BTB.std <- BuyseTest(treatment ~ tte(eventtime, status) + stratum, trace = FALSE, data = dtB.strata,
pool.strata = "standardization",
model.tte = prodlim(Hist(eventtime, status) ~ treatment, data = dtB.strata))
## estimate se lower.ci upper.ci null p.value
## eventtime 0.02881503 0.1325542 -0.2271614 0.2810671 0 0.8280037
expect_equivalent(confint(BTB.std), data.frame("estimate" = c(0.02881503),"se" = c(0.13255421),"lower.ci" = c(-0.22716139),"upper.ci" = c(0.28106715),"null" = c(0),"p.value" = c(0.82800367)), tol = 1e-5)
})
##----------------------------------------------------------------------
### test-BuyseTest-standardization.R ends here
Try the BuyseTest package in your browser
Any scripts or data that you put into this service are public.
BuyseTest documentation built on June 8, 2025, 11:19 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.
### test-BuyseTest-standardization.R ---
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: apr 11 2025 (15:42)
## Version:
## Last-Updated: May 8 2025 (10:27)
## By: Brice Ozenne
## Update #: 105
##----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
##----------------------------------------------------------------------
##
### Code:
if(FALSE){
library(testthat)
library(BuyseTest)
library(data.table)
library(survival)
library(prodlim)
}
context("Check standardization of the net benefit across strata \n")
## * Setting
BuyseTest.options(order.Hprojection = 1,
method.inference = "u-statistic",
trace = 0)
## * Binary outcome and binary covariate
## ** simulate data
seed <- 1
n <- 100
tau <- 0
set.seed(seed)
A <- rbinom(n, size = 1, prob = 0.75)
X <- sample(0:1, size = n, replace = TRUE)
df <- data.frame(id = 1:n,
Y = 0.5 * A + X + rnorm(n),
arm = A, X = X)
## ** standardization (first order)
test_that("BuyseTest - standarization (first order, binary outcome)",{
## 1 - marginal
## BT.unstrat <- BuyseTest(data = df, arm ~ cont(Y, threshold = tau), trace = F)
## BT.unstrat2 <- BuyseTest(data = df, arm ~ cont(Y, threshold = tau), trace = TRUE, method.inference = "bootstrap")
## rbind(confint(BT.unstrat, transformation = FALSE),
## confint(BT.unstrat2, method.ci = "gaussian"))
## estimate se lower.ci upper.ci null p.value
## Y 0.1811263 0.1188404 -0.05179662 0.4140493 0 0.1274803
## Y1 0.1811263 0.1220947 -0.06875553 0.4096230 0 0.1543344
## 2- standardized
BT.std <- BuyseTest(data = df, arm ~ cont(Y, threshold = tau) + X, pool.strata = "standardization", trace = FALSE)
expect_equal(mean(getIid(BT.std)),0, tol = 1e-6)
expect_equivalent(confint(BT.std, transformation = FALSE),
data.frame("estimate" = c(0.23088881),"se" = c(0.10551968),"lower.ci" = c(0.02407403),"upper.ci" = c(0.43770359),"null" = c(0),"p.value" = c(0.02866149)),,
tol = 1e-6)
expect_equivalent(coef(BT.std, strata = TRUE), c("0" = 0.03658537, "1.0" = -0.30406504, "0.1" = 0.75, "1" = 0.49583333))
expect_error(confint(BT.std, strata = TRUE)) ## cannot compute strata specific uncertainty since the same individual appear in sevearl strata
BT.std2 <- BuyseTest(data = df, arm ~ cont(Y, threshold = tau) + X, pool.strata = "standardization", trace = FALSE,
method.inference = "bootstrap", strata.resampling = c("arm","X"), seed = 10)
expect_equivalent(confint(BT.std2, method.ci = "gaussian", transformation = FALSE),
data.frame("estimate" = c(0.23088881),"se" = c(0.10915709),"lower.ci" = c(0.01694484),"upper.ci" = c(0.44483278),"null" = c(0),"p.value" = c(0.03441312)),
tol = 1e-6)
## 3- stratified to stratified (by hand)
df00 <- df[df$X==0,]
df10 <- rbind(df[df$X==0 & df$arm==1,],df[df$X==1 & df$arm==0,])
df01 <- rbind(df[df$X==0 & df$arm==0,],df[df$X==1 & df$arm==1,])
df11 <- df[df$X==1,]
BT.strat00 <- BuyseTest(arm ~ cont(Y, threshold = tau), method.inference = "U-statistic", trace = F, data = df00)
BT.strat10 <- BuyseTest(arm ~ cont(Y, threshold = tau), method.inference = "U-statistic", trace = F, data = df10)
BT.strat01 <- BuyseTest(arm ~ cont(Y, threshold = tau), method.inference = "U-statistic", trace = F, data = df01)
BT.strat11 <- BuyseTest(arm ~ cont(Y, threshold = tau), method.inference = "U-statistic", trace = F, data = df11)
grid.strata <- expand.grid(X = sort(unique(df$X)), Xstar = sort(unique(df$X)))
grid.strata$X.X <- paste(grid.strata$X,grid.strata$Xstar,sep=".")
df.extended <- do.call(rbind,apply(grid.strata, MARGIN = 1, function(iRow){
cbind(rbind(df[df$arm == 0 & df$X == iRow[1],],
df[df$arm == 1 & df$X == iRow[2],]),
X.X = paste(iRow[1],iRow[2],sep="."))
}))
grid.strata$nC.X <- tapply(df.extended$arm==0,df.extended$X.X,sum)[grid.strata$X.X]
grid.strata$nE.Xstar <- tapply(df.extended$arm==1,df.extended$X.X,sum)[grid.strata$X.X]
grid.strata$n.X <- as.numeric(table(df$X)[as.character(grid.strata$X)])
grid.strata$n.Xstar <- as.numeric(table(df$X)[as.character(grid.strata$Xstar)])
grid.strata$w.unadjusted <- grid.strata$nC.X * grid.strata$nE.Xstar / prod(table(df$arm))
grid.strata$w.standardize <- grid.strata$n.X * grid.strata$n.Xstar / NROW(df$arm)^2
## direct
Myiid <- matrix(0, nrow = NROW(df), ncol = 4, dimnames = list(NULL,unique(df.extended$X.X)))
Myiid[df00$id,"0.0"] <- getIid(BT.strat00)
Myiid[df10$id,"1.0"] <- getIid(BT.strat10)
Myiid[df01$id,"0.1"] <- getIid(BT.strat01)
Myiid[df11$id,"1.1"] <- getIid(BT.strat11)
MyWiid <- sweep(Myiid, FUN = "*", MARGIN = 2, STATS = grid.strata$w.standardize)
## mean(rowSums(MyWiid))
expect_equivalent(BT.std@covariance[,"netBenefit"], crossprod(rowSums(MyWiid)), tol = 1e-6)
## from the partial win/loss sums
Myiid2 <- matrix(0, nrow = NROW(df), ncol = 4, dimnames = list(NULL,unique(df.extended$X.X)))
Myiid2[df00$id,"0.0"] <- getIid(BT.strat00, scale = FALSE, center = FALSE) * table(df00$arm)[as.character(1-df00$arm)]
Myiid2[df10$id,"1.0"] <- getIid(BT.strat10, scale = FALSE, center = FALSE) * table(df10$arm)[as.character(1-df10$arm)]
Myiid2[df01$id,"0.1"] <- getIid(BT.strat01, scale = FALSE, center = FALSE) * table(df01$arm)[as.character(1-df01$arm)]
Myiid2[df11$id,"1.1"] <- getIid(BT.strat11, scale = FALSE, center = FALSE) * table(df11$arm)[as.character(1-df11$arm)]
## MyWiid2 <- sweep(Myiid2, FUN = "*", MARGIN = 2, STATS = grid.strata$w.standardize)
## MyWiid2.norm <- matrix(0, nrow = NROW(df), ncol = 4, dimnames = list(NULL,unique(df.extended$X.X)))
## MyWiid2.norm[df00$id,"0.0"] <- ((MyWiid2[df00$id,"0.0"] / table(df00$arm)[as.character(1-df00$arm)]) - coef(BT.strat00)*BT.std@weightStrata[1])/table(df00$arm)[as.character(df00$arm)]
## MyWiid2.norm[df10$id,"1.0"] <- ((MyWiid2[df10$id,"1.0"] / table(df10$arm)[as.character(1-df10$arm)]) - coef(BT.strat10)*BT.std@weightStrata[2])/table(df10$arm)[as.character(df10$arm)]
## MyWiid2.norm[df01$id,"0.1"] <- ((MyWiid2[df01$id,"0.1"] / table(df01$arm)[as.character(1-df01$arm)]) - coef(BT.strat01)*BT.std@weightStrata[3])/table(df01$arm)[as.character(df01$arm)]
## MyWiid2.norm[df11$id,"1.1"] <- ((MyWiid2[df11$id,"1.1"] / table(df11$arm)[as.character(1-df11$arm)]) - coef(BT.strat11)*BT.std@weightStrata[4])/table(df11$arm)[as.character(df11$arm)]
## mean(rowSums(MyWiid2.norm))
## crossprod(rowSums(MyWiid2.norm))
## this is what is implemented in the C++ function
MyWiid3 <- rowSums(sweep(Myiid2, FUN = "*", MARGIN = 2, STATS = grid.strata$w.standardize/(grid.strata$w.unadjusted*prod(table(df$arm)))))
## grid.strata$w.unadjusted*prod(table(df$arm)) - c(prod(table(df00$arm)),prod(table(df10$arm)),prod(table(df01$arm)),prod(table(df11$arm)))
MyWiid3[df00$id] <- MyWiid3[df00$id] - coef(BT.strat00)*BT.std@weightStrata[1]/table(df00$arm)[as.character(df00$arm)]
MyWiid3[df10$id] <- MyWiid3[df10$id] - coef(BT.strat10)*BT.std@weightStrata[2]/table(df10$arm)[as.character(df10$arm)]
MyWiid3[df01$id] <- MyWiid3[df01$id] - coef(BT.strat01)*BT.std@weightStrata[3]/table(df01$arm)[as.character(df01$arm)]
MyWiid3[df11$id] <- MyWiid3[df11$id] - coef(BT.strat11)*BT.std@weightStrata[4]/table(df11$arm)[as.character(df11$arm)]
expect_equal(mean(MyWiid3), 0, tol = 1e-6)
expect_equivalent(BT.std@covariance[,"netBenefit"], crossprod(MyWiid3), tol = 1e-6)
})
## ** standardization (second order)
BuyseTest.options(order.Hprojection = 2)
test_that("BuyseTest - standarization (second order, binary outcome, balanced)",{
## with balanced dataset standardization equals marginal analysis
dfB <- do.call(rbind,by(df, interaction(df$arm,df$X), function(iDF){iDF[1:min(table(df$arm,df$X)),,drop=FALSE]}))
dfB$id <- 1:NROW(dfB)
rownames(dfB) <- NULL
BTB.marginal <- BuyseTest(data = dfB, arm ~ cont(Y, threshold = tau), trace = FALSE)
BTB.marginal2 <- BuyseTest(data = dfB, arm ~ cont(Y, threshold = tau), trace = F, keep.pairScore = TRUE)
BTB.std <- BuyseTest(data = dfB, arm ~ cont(Y, threshold = tau) + X, pool.strata = "standardization", trace = F, keep.pairScore = FALSE)
BTB.std2 <- BuyseTest(data = dfB, arm ~ cont(Y, threshold = tau) + X, pool.strata = "standardization", trace = F, keep.pairScore = TRUE)
expect_equal(coef(BTB.marginal), coef(BTB.std), tol = 1e-5)
expect_equal(coef(BTB.std), mean(coef(BTB.std, strata = TRUE)), tol = 1e-5)
expect_true(confint(BTB.marginal)$se > confint(BTB.marginal, order.Hprojection = 1)$se) ## confint(BTB.marginal)$se^2 - confint(BTB.marginal, order.Hprojection = 1)$se^2
expect_true(confint(BTB.std)$se > confint(BTB.std, order.Hprojection = 1)$se) ## confint(BTB.std)$se^2 - confint(BTB.std, order.Hprojection = 1)$se^2
expect_equivalent(confint(BTB.std), confint(BTB.std2), tol = 1e-6)
## by hand (marginal)
sigmaMB_h1fav <- tapply(getIid(BTB.marginal, statistic = "favorable"), dfB$arm, crossprod)
termM1 <- (coef(BTB.marginal, statistic = "favorable") - coef(BTB.marginal, statistic = "favorable")^2)/prod(table(dfB$arm)) - sum(sigmaMB_h1fav/rev(table(dfB$arm)))
sigmaMB_h1unfav <- tapply(getIid(BTB.marginal, statistic = "unfavorable"), dfB$arm, crossprod)
termM2 <- (coef(BTB.marginal, statistic = "unfavorable") - coef(BTB.marginal, statistic = "unfavorable")^2)/prod(table(dfB$arm)) - sum(sigmaMB_h1unfav/rev(table(dfB$arm)))
sigmaMB_h1cross <- tapply(getIid(BTB.marginal, statistic = "favorable")*getIid(BTB.marginal, statistic = "unfavorable"), dfB$arm, sum)
termM3 <- 0 - coef(BTB.marginal, statistic = "favorable")*coef(BTB.marginal, statistic = "unfavorable")/prod(table(dfB$arm)) - sum(sigmaMB_h1cross/rev(table(dfB$arm)))
expect_equal(termM1+termM2-2*termM3, confint(BTB.marginal, order.Hprojection = 2)$se^2 - confint(BTB.marginal, order.Hprojection = 1)$se^2, tol = 1e-6)
mytableM <- getPairScore(BTB.marginal2)
HprojM_fav <- getIid(BTB.marginal2, statistic = "favorable", scale = FALSE)
estM_fav <- coef(BTB.marginal2, statistic = "favorable")
mytableM[, hproj.0 := HprojM_fav[index.0]]
mytableM[, hproj.1 := HprojM_fav[index.1]]
expect_equal(termM1, mytableM[, sum((favorable - hproj.0 - hproj.1 - estM_fav)^2)/.N^2], tol = 1e-6)
## by hand (standardized)
sigmaB_h1fav <- tapply(getIid(BTB.std, statistic = "favorable"), dfB$arm, crossprod)
term1 <- (coef(BTB.std, statistic = "favorable") - coef(BTB.std, statistic = "favorable")^2)/prod(table(dfB$arm)) - sum(sigmaB_h1fav/rev(table(dfB$arm)))
sigmaB_h1unfav <- tapply(getIid(BTB.std, statistic = "unfavorable"), dfB$arm, crossprod)
term2 <- (coef(BTB.std, statistic = "unfavorable") - coef(BTB.std, statistic = "unfavorable")^2)/prod(table(dfB$arm)) - sum(sigmaB_h1unfav/rev(table(dfB$arm)))
sigmaB_h1cross <- tapply(getIid(BTB.std, statistic = "favorable")*getIid(BTB.std, statistic = "unfavorable"), dfB$arm, sum)
term3 <- 0 - coef(BTB.std, statistic = "favorable")*coef(BTB.std, statistic = "unfavorable")/prod(table(dfB$arm)) - sum(sigmaB_h1cross/rev(table(dfB$arm)))
expect_equal(term1+term2-2*term3, confint(BTB.std, order.Hprojection = 2)$se^2 - confint(BTB.std, order.Hprojection = 1)$se^2, tol = 1e-6)
mytable <- getPairScore(BTB.std2)
mytable[, weightStrata := stats::setNames(BTB.std2@weightStrata,BTB.std2@level.strata)[strata]]
mytable[, m.X := length(unique(index.0)), by = "strata"]
mytable[, n.X := length(unique(index.1)), by = "strata"]
mytable[, m := length(attr(BTB.std2@level.treatment,"indexC"))]
mytable[, n := length(attr(BTB.std2@level.treatment,"indexT"))]
mytable[, Wfavorable := favorable * weightStrata/((m.X*n.X)/(n*m))]
## mean(mytable$Wfavorable)
## coef(BTB.std2, statistic = "favorable")
Hproj_fav <- getIid(BTB.std2, statistic = "favorable", scale = FALSE)
mytable[, hproj.0 := Hproj_fav[index.0]]
mytable[, hproj.1 := Hproj_fav[index.1]]
mytable[, Tfavorable := coef(BTB.std2, statistic = "favorable")]
expect_equal(term1, mytable[, sum((Wfavorable - hproj.0 - hproj.1 - Tfavorable)^2)/.N^2], tol = 1e-6)
})
test_that("BuyseTest - standarization (second order, binary outcome, unbalanced)",{
## 2- standardized
BTUB.std <- BuyseTest(data = df, arm ~ cont(Y, threshold = tau) + X, pool.strata = "standardization", trace = FALSE)
BTUB.std2 <- BuyseTest(data = df, arm ~ cont(Y, threshold = tau) + X, pool.strata = "standardization", keep.pairScore = TRUE, trace = FALSE)
expect_equivalent(confint(BTUB.std), confint(BTUB.std2), tol = 1e-6)
expect_equivalent(confint(BTUB.std), data.frame("estimate" = c(0.23088881),"se" = c(0.10675249),"lower.ci" = c(0.01411401),"upper.ci" = c(0.42693406),"null" = c(0),"p.value" = c(0.03705691)), tol = 1e-6)
## by hand (standardized)
mytableUB <- getPairScore(BTUB.std)
mytableUB[, weightStrata := stats::setNames(BTUB.std@weightStrata,BTUB.std@level.strata)[strata]]
mytableUB[, m.X := length(unique(index.0)), by = "strata"]
mytableUB[, n.X := length(unique(index.1)), by = "strata"]
mytableUB[, m := length(attr(BTUB.std@level.treatment,"indexC"))]
mytableUB[, n := length(attr(BTUB.std@level.treatment,"indexT"))]
mytableUB[, Wfavorable := favorable * weightStrata/((m.X*n.X)/(n*m))]
## mean(mytableUB$Wfavorable)
## coef(BTUB.std, statistic = "favorable")
Hproj_fav <- getIid(BTUB.std2, statistic = "favorable", scale = FALSE)
mytableUB[, hproj.0 := Hproj_fav[index.0]]
mytableUB[, hproj.1 := Hproj_fav[index.1]]
mytableUB[, Tfavorable := coef(BTUB.std2, statistic = "favorable")]
term1 <- confint(BTUB.std, statistic = "favorable", order.Hprojection = 2)$se^2 - confint(BTUB.std, statistic = "favorable", order.Hprojection = 1)$se^2
expect_equal(term1, mytableUB[, sum((Wfavorable - hproj.0 - hproj.1 - Tfavorable)^2)/.N^2], tol = 1e-6)
})
BuyseTest.options(order.Hprojection = 1)
## * check Peron scoring rule in a perfectly balanced dataset
n <- 50
set.seed(10)
dt.strata <- simBuyseTest(n,
n.strata = 2,
names.strata = "stratum")
dtB.strata <- do.call(rbind,by(dt.strata, interaction(dt.strata$treatment,dt.strata$stratum), function(iDT){iDT[1:min(table(dt.strata$treatment,dt.strata$stratum))]}))
rownames(dtB.strata) <- NULL
test_that("BuyseTest - standarization (tte outcome)",{
BT.raw <- BuyseTest(treatment ~ tte(eventtime, status), trace = FALSE, data = dt.strata)
confint(BT.raw)
## estimate se lower.ci upper.ci null p.value
## eventtime -0.008226969 0.1286523 -0.2546679 0.2392174 0 0.9490145
BT.std <- BuyseTest(treatment ~ tte(eventtime, status) + stratum, trace = FALSE, data = dt.strata,
pool.strata = "standardization")
## estimate se lower.ci upper.ci null p.value
## eventtime -0.009499187 0.127931 -0.2545409 0.2366887 0 0.9408131
expect_equivalent(confint(BT.std), data.frame("estimate" = c(-0.00949919),"se" = c(0.12793103),"lower.ci" = c(-0.25454086),"upper.ci" = c(0.23668869),"null" = c(0),"p.value" = c(0.9408131)), tol = 1e-5)
BTB.raw <- BuyseTest(treatment ~ tte(eventtime, status), trace = FALSE, data = dtB.strata,
model.tte = prodlim(Hist(eventtime, status) ~ treatment, data = dtB.strata))
confint(BTB.raw)
## estimate se lower.ci upper.ci null p.value
## eventtime 0.02881503 0.1417496 -0.2441965 0.2975942 0 0.8390032
BTB.std <- BuyseTest(treatment ~ tte(eventtime, status) + stratum, trace = FALSE, data = dtB.strata,
pool.strata = "standardization",
model.tte = prodlim(Hist(eventtime, status) ~ treatment, data = dtB.strata))
## estimate se lower.ci upper.ci null p.value
## eventtime 0.02881503 0.1325542 -0.2271614 0.2810671 0 0.8280037
expect_equivalent(confint(BTB.std), data.frame("estimate" = c(0.02881503),"se" = c(0.13255421),"lower.ci" = c(-0.22716139),"upper.ci" = c(0.28106715),"null" = c(0),"p.value" = c(0.82800367)), tol = 1e-5)
})
##----------------------------------------------------------------------
### test-BuyseTest-standardization.R ends here
Try the BuyseTest 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.