Risk Comparison Over Time
Comparison of risk differences or risk ratios over all timepoints.
## S3 method for class 'ate' anova( object, allContrast = NULL, type = "diff", estimator = object$estimator[1], test = "CvM", transform = NULL, alternative = "two.sided", n.sim = 10000, print = TRUE, ... )
object |
A |
allContrast |
[matrix] contrast for which the risks should be compared. Matrix with two rows, the first being the sequence of reference treatments and the second the sequence of alternative treatments. |
type |
[character vector] the functionnal used to compare the risks: |
estimator |
[character] The type of estimator relative to which the comparison should be performed. |
test |
[character] The type of statistic used to compare the risks over times:
|
transform |
[character] Should a transformation be used, e.g. the test is performed after log-transformation of the estimate, standard error, and influence function. |
alternative |
[character] a character string specifying the alternative hypothesis, must be one of |
n.sim |
[integer, >0] the number of simulations used to compute the p-values. |
print |
[logical] should the results be displayed? |
... |
Not used. |
Experimental!!!
library(survival)
library(data.table)
## Not run:
## simulate data
set.seed(12)
n <- 200
dtS <- sampleData(n,outcome="survival")
dtS$X12 <- LETTERS[as.numeric(as.factor(paste0(dtS$X1,dtS$X2)))]
dtS <- dtS[dtS$X12!="D"]
## model fit
fit <- cph(formula = Surv(time,event)~ X1+X6,data=dtS,y=TRUE,x=TRUE)
seqTime <- 1:10
ateFit <- ate(fit, data = dtS, treatment = "X1", contrasts = NULL,
times = seqTime, B = 0, iid = TRUE, se = TRUE, verbose = TRUE, band = TRUE)
## display
autoplot(ateFit)
## inference (two sided)
statistic <- ateFit$diffRisk$estimate/ateFit$diffRisk$se
confint(ateFit, p.value = TRUE, method.band = "bonferroni")$diffRisk
confint(ateFit, p.value = TRUE, method.band = "maxT-simulation")$diffRisk
anova(ateFit, test = "KS")
anova(ateFit, test = "CvM")
anova(ateFit, test = "sum")
## manual calculation (one sided)
n.sim <- 1e4
statistic <- ateFit$diffRisk[, estimate/se]
iid.norm <- scale(ateFit$iid$GFORMULA[["1"]]-ateFit$iid$GFORMULA[["0"]],
scale = ateFit$diffRisk$se)
ls.out <- lapply(1:n.sim, function(iSim){
iG <- rnorm(NROW(iid.norm))
iCurve <- t(iid.norm) %*% iG
data.table(max = max(iCurve), L2 = sum(iCurve^2), sum = sum(iCurve),
maxC = max(iCurve) - max(statistic),
L2C = sum(iCurve^2) - sum(statistic^2),
sumC = sum(iCurve) - sum(statistic),
sim = iSim)
})
dt.out <- do.call(rbind,ls.out)
dt.out[,.(max = mean(.SD$maxC>=0),
L2 = mean(.SD$L2C>=0),
sum = mean(.SD$sumC>=0))]
## permutation
n.sim <- 250
stats.perm <- vector(mode = "list", length = n.sim)
pb <- txtProgressBar(max = n.sim, style=3)
treatVar <- ateFit$variables["treatment"]
for(iSim in 1:n.sim){ ## iSim <- 1
iData <- copy(dtS)
iIndex <- sample.int(NROW(iData), replace = FALSE)
iData[, c(treatVar) := .SD[[treatVar]][iIndex]]
iFit <- update(fit, data = iData)
iAteSim <- ate(iFit, data = iData, treatment = treatVar,
times = seqTime, verbose = FALSE)
iStatistic <- iAteSim$diffRisk[,estimate/se]
stats.perm[[iSim]] <- cbind(iAteSim$diffRisk[,.(max = max(iStatistic),
L2 = sum(iStatistic^2),
sum = sum(iStatistic))],
sim = iSim)
stats.perm[[iSim]]$maxC <- stats.perm[[iSim]]$max - max(statistic)
stats.perm[[iSim]]$L2C <- stats.perm[[iSim]]$L2 - sum(statistic^2)
stats.perm[[iSim]]$sumC <- stats.perm[[iSim]]$sum - sum(statistic)
setTxtProgressBar(pb, iSim)
}
dtstats.perm <- do.call(rbind,stats.perm)
dtstats.perm[,.(max = mean(.SD$maxC>=0),
L2 = mean(.SD$L2C>=0),
sum = mean(.SD$sumC>=0))]
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.