I want to study the pvalue of the variables included in my Cox model.
I compare this two results :
library(survival)
library(rms)
n <- 1000
set.seed(731)
age <- 50 + 12*rnorm(n)
label(age) <- "Age"
sex <- factor(sample(c('Male','Female'), n,
rep=TRUE, prob=c(.6, .4)))
var1 = rnorm(n)*2
var2 = rnorm(n)
var3 = rnorm(n)*5 + 1*(sex=="Female") - 0.5
var4 = rbinom(n,1,0.3)
var5 = rbinom(n,3,0.6)
var6 = rbinom(n,3,0.2)
dta = data.frame(age=(age-50),sex=(sex=='Female'),var1,var2,var3,var4,var5,var6)
coef = sample(size = ncol(dta),replace = T,x = c(-1,-0.5,0,0.5,1))
cens <- 5*runif(n)
h <- 0.5*exp(rowSums(sapply(1:length(coef),function(i) coef[i]*dta[,i] )))
dt <- -log(runif(n))/h
dt = dt/365
dt[dt>10] = 10
hist(dt,breaks = 20)
label(dt) <- 'Follow-up Time'
e <- ifelse(dt <= cens,1,0)
dt <- pmin(dt, cens)
units(dt) <- "Year"
dta = data.frame(dt,e,age=(age-50),sex=(sex=='Female'),var1,var2,var3,var4,var5,var6)
dd <- datadist(dta)
options(datadist='dd')
my_formula <- as.formula(paste("Surv(dt, e) ~ ", paste(colnames(dta[,3:ncol(dta)]), collapse = " + ")))
f2 <- coxph(my_formula, data = dta)
f <- cph(my_formula , data = dta, x=TRUE, y=TRUE)
summary(f2)
anova(f)
When I look at the output of the "rms" package (comparing with the results given by the "survival" package), I get two different overall Wald test values. In order to study my variables, I use this function:
anova(f)
How can I do my test based on my Wald value from the survival package?
