# tests/testthat/test-wkm.R In mbrueckner/AHR: Estimation and Testing of Average Hazard Ratios

```library(testthat)

context("wkm")

## TODO: test extreme cases (no censoring, all obs. censored, ...)

## TODO: test with data from Kalbfleisch/Prentice paper

test_that("single stratum wkm estimate reduces to Kaplan-Meier estimator", {
test.data <- function(n) {
Z <- rep(0, n)
T <- rexp(n, 0.1)
C <- runif(n, 0, 10)
X <- pmin(T, C)
D <- T <= C
V <- runif(n, 0, X/4)

data.frame(Y=X, D=D, W=Z, V=V)
}

data <- test.data(100)
times <- seq(0, 5, length.out=10) ##c(1.5, 3)

data <- data[order(data\$Y),]

S <- exp(-times/10)

fs <- survfit(Surv(V, Y, D) ~ 1, data=data)
fit <- wkm(times, formula=Surv(V, Y, D) ~ strata(W), data=data)

f <- approxfun(fs\$time, fs\$surv, method="constant", yleft=1, rule=2, f=0)
g <- approxfun(fs\$time, 100 * (fs\$surv * fs\$std.err)^2, method="constant", yleft=0, rule=2, f=0)

expect_true(all.equal(f(times), fit\$S))
expect_true(all.equal(g(times), fit\$V))
})

test_that("covariance matrix is symmetric", {
X <- sample(c(0,1), 100, replace=TRUE)
T <- numeric(100)
T[X == 0] <- rexp(sum(X == 0), 0.1)
T[X == 1] <- rexp(sum(X == 1), 0.5)

C <- runif(100, 0, 10)

Y <- pmin(T, C)
D <- T <= C
V <- runif(100, 0, Y/2)

data <- data.frame(Y=Y, D=D, W=X, V=V)

expect_true(isSymmetric(wkm(sort(Y), formula=Surv(V, Y, D) ~ strata(W), data=data, param=list(cov=TRUE, var=TRUE, alpha=1, left.limit=FALSE))\$COV))
})

test_that("variance and covariance calculations match", {
T <- rexp(200)
C <- rexp(200)
D <- T <= C
Y <- pmin(T, C)

Z <- rbinom(200, 1, 0.5)
X <- rbinom(200, 1, 0.5)

data <- data.frame(Y=Y, D=D, Trt=Z, W=factor(X))

times <- sort(data\$Y[data\$D == 1])

fit <- wkm(times, formula=Surv(Y, D) ~ strata(W), data=data, param=list(cov=TRUE, var=FALSE, alpha=1, left.limit=FALSE))
fit2 <- wkm(times, formula=Surv(Y, D) ~ strata(W), data=data, param=list(cov=FALSE, var=TRUE, alpha=1, left.limit=FALSE))

expect_true(all.equal(diag(fit\$COV), fit2\$V))
})
```
mbrueckner/AHR documentation built on May 22, 2017, 12:16 a.m.