### test-phreg.R ---
#----------------------------------------------------------------------
## author: Brice Ozenne
## created: feb 28 2017 (09:52)
## Version:
## last-updated: Sep 17 2022 (07:02)
## By: Thomas Alexander Gerds
## Update #: 31
#----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
#----------------------------------------------------------------------
##
### Code:
## * Settings
library(testthat)
library(riskRegression)
library(rms)
library(survival)
library(prodlim)
library(mets)
context("Compatibility with the phreg function from the mets package")
## * Survival
## ** Data
set.seed(10)
d <- sampleData(5e1, outcome = "survival")[,.(eventtime,event,X1,X2,X3,X4,X6)]
d[ , Xcat2 := as.factor(paste0(X1,X2))]
d[, entry:= 0]
d[, id := 1:.N]
d[, group := factor(rbinom(.N,2,0.5))]
## ** Tests
test_that("survival - no strata", {
m.phreg <- phreg(Surv(entry,eventtime,event)~X1+X2,data=d)
m.coxph <- coxph(Surv(eventtime,event)~X1+X2,data=d, x = TRUE, y = TRUE) ## do not center anymore factor or binary variables
m.cph <- cph(Surv(eventtime,event)~X1+X2,data=d, x = TRUE, y = TRUE)
## slight difference in estimated coefficients: coef(m.phreg) - coef(m.cph)
expect_equal(ignore_attr=TRUE,predictCox(m.phreg, center = FALSE),predictCox(m.coxph, center = FALSE), tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,predictCox(m.phreg, center = TRUE),predictCox(m.cph, center = TRUE), tolerance = 1e-4)
pred.phreg <- predictCox(m.phreg, newdata = d, time = 1:5, se = TRUE)
pred.coxph <- predictCox(m.coxph, newdata = d, time = 1:5, se = TRUE)
expect_equal(ignore_attr=TRUE,pred.phreg$survival, pred.coxph$survival, tolerance = 1e-8)
})
test_that("survival - one strata variable", {
mS.phreg <- phreg(Surv(entry,eventtime,event)~X1+X2+strata(group)+cluster(id),data=d)
mS.coxph <- coxph(Surv(eventtime,event)~X1+X2+strata(group),data=d, x = TRUE, y = TRUE) ## do not center anymore factor or binary variables
mS.cph <- cph(Surv(eventtime,event)~X1+X2+strat(group),data=d, x = TRUE, y = TRUE)
## slight difference in estimated coefficients: coef(mS.phreg) - coef(mS.cph)
expect_equal(ignore_attr=TRUE,predictCox(mS.phreg, center = FALSE),predictCox(mS.coxph, center = FALSE), tolerance = 1e-8)
expect_equal(ignore_attr=TRUE,predictCox(mS.phreg, center = TRUE)[c("cumhazard","survival")],predictCox(mS.cph, center = TRUE)[c("cumhazard","survival")], tolerance = 1e-4)
predS.phreg <- predictCox(mS.phreg, newdata = d, time = 1:5, se = TRUE)
predS.coxph <- predictCox(mS.coxph, newdata = d, time = 1:5, se = TRUE)
expect_equal(ignore_attr=TRUE,predS.phreg$survival, predS.coxph$survival, tolerance = 1e-8)
## bug in previous version, could not create design matrix when group had only one value
predS.phreg <- predictCox(mS.phreg, newdata = d[group==0], time = 1:5, se = TRUE)
predS.coxph <- predictCox(mS.coxph, newdata = d[group==0], time = 1:5, se = TRUE)
expect_equal(ignore_attr=TRUE,predS.phreg$survival, predS.coxph$survival, tolerance = 1e-8)
})
## several strata variables
test_that("survival - several strata variables", {
mS.phreg <- phreg(Surv(entry,eventtime,event)~X1+X2+strata(group,X3)+cluster(id),data=d)
mS.coxph <- coxph(Surv(eventtime,event)~X1+X2+strata(group)+strata(X3),data=d, x = TRUE, y = TRUE) ## do not center anymore factor or binary variables
mS.cph <- cph(Surv(eventtime,event)~X1+X2+strat(group)+strat(X3),data=d, x = TRUE, y = TRUE)
## slight difference in estimated coefficients: coef(mS.phreg) - coef(mS.cph)
expect_equal(ignore_attr=TRUE,predictCox(mS.phreg, center = FALSE)[c("cumhazard","survival")],
predictCox(mS.coxph, center = FALSE)[c("cumhazard","survival")],
tolerance = 1e-8)
## FAIL BECAUSE DIFFERENT ORDERING OF THE STRATA
## expect_equal(ignore_attr=TRUE,predictCox(mS.phreg, center = TRUE)[c("cumhazard","survival")],
## predictCox(mS.cph, center = TRUE)[c("cumhazard","survival")],
## tolerance = 1e-4)
predS.phreg <- predictCox(mS.phreg, newdata = d, time = 1:5, se = TRUE)
predS.coxph <- predictCox(mS.coxph, newdata = d, time = 1:5, se = TRUE)
expect_equal(ignore_attr=TRUE,predS.phreg$survival, predS.coxph$survival, tolerance = 1e-8)
})
## * Competing risks
## ** Data
set.seed(10)
d <- sampleData(1e2, outcome = "competing.risks")[,.(time,event,X1,X2,X3,X4,X6)]
d[ , Xcat2 := as.factor(paste0(X1,X2))]
d[, id := 1:.N]
d[, group := factor(rbinom(.N,2,0.5))]
## ** Tests
test_that("competing risk - no strata", {
m.phreg <- CSC(Hist(time, event)~X1+X2,data=d,
fitter = "phreg")
m.coxph <- CSC(Hist(time,event)~X1+X2,data=d,
fitter = "coxph")
expect_equal(ignore_attr=TRUE,class(m.phreg$models[[1]]),"phreg")
pred.phreg <- predict(m.phreg, newdata = d, times = 1:5, cause = 1, se = TRUE)
pred.coxph <- predict(m.coxph, newdata = d, times = 1:5, cause = 1, se = TRUE)
expect_equal(ignore_attr=TRUE,pred.phreg,pred.coxph, tolerance = 1e-8)
})
test_that("competing risk - one strata variable", {
mS.phreg <- CSC(Hist(time, event)~X1+X2+strata(group)+cluster(id),data=d,
fitter = "phreg")
mS.coxph <- CSC(Hist(time,event)~X1+X2+strata(group),data=d,
fitter = "coxph")
expect_equal(ignore_attr=TRUE,class(mS.phreg$models[[1]]),"phreg")
predS.phreg <- predict(mS.phreg, newdata = d, times = 1:5, cause = 1, se = TRUE)
predS.coxph <- predict(mS.coxph, newdata = d, times = 1:5, cause = 1, se = TRUE)
expect_equal(ignore_attr=TRUE,predS.phreg,predS.coxph, tolerance = 1e-3)
})
test_that("competing risk - several strata variables", {
mS.phreg <- CSC(Hist(time, event)~X1+X2+strata(group,X3)+cluster(id),data=d,
fitter = "phreg")
mS.coxph <- CSC(Hist(time,event)~X1+X2+strata(group)+strata(X3),data=d,
fitter = "coxph")
expect_equal(ignore_attr=TRUE,class(mS.phreg$models[[1]]),"phreg")
predS.phreg <- predict(mS.phreg, newdata = d, times = 1:5, cause = 1, se = TRUE)
predS.coxph <- predict(mS.coxph, newdata = d, times = 1:5, cause = 1, se = TRUE)
expect_equal(ignore_attr=TRUE,predS.phreg,predS.coxph, tolerance = 1e-3)
})
##----------------------------------------------------------------------
### test-phreg.R ends here
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.