Nothing
### 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
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.