0
$\begingroup$

It's been years since I've looked at this or the R language, so I'd appreciate some help and simple explanations (I'm not a statistician).

Using bootstrap analysis, I'm trying to derive the p-value of a Cox regression for a particular coefficient. Unfortunately, the derived p-value is far from the p-value returned from the coxph survival R function. Yet the HR (exp of the coefficient of interest) is very close, so I know I'm not far off.

This is the R code:

library(boot)
library(survival)

bootCox <- function(data, indices)  {
  sampled_indices <- sample(nrow(data), 7684, replace = TRUE)  
  data <- data[sampled_indices, ]
  
  y <- "Surv(time_to_event, status==1)"
  form <- as.formula(paste(y, "~", "exposed + sex"))
  coxModel <- coxph(form, data = data)
  exposedCoef <- exp(coef(coxModel)["exposed"])
  return(exposedCoef)
}

#preparing the survival data
data <- read.table("atenolol_DF.txt_psm.txt", header=TRUE, sep="\t", stringsAsFactors=TRUE, fill=TRUE)

#run bootstrap analysis
set.seed(42)
bootResult <- boot(data, statistic = bootCox, R = 10000)

#calculate the HR once
y <- "Surv(time_to_event, status==1)"
form <- as.formula(paste(y, "~", "exposed + sex"))
coxModel <- coxph(form, data = data)
observedCoefficient <- exp(coef(coxModel)["exposed"])

p_value <- sum(bootResult$t >= observedCoefficient) / length(bootResult$t)

The HR for exposed is 0.82 using the bootstrap approach and a single call of coxph. Yet, the corresponding P-values are very different. The p-value for the single call is 0.02, whilst for the bootstrap approach, it's 0.51.

The total sample size is 7684 patients, of which around one-fifth were exposed. Inside the bootstrap function bootCox, I randomly select that number of patients with replacement. Then, I call the bootstrap function 10,000 times. I've played around with these two values, and unless I make them very small, the derived HR and p-value don't change much.

Any idea why I'm seeing such a difference? Thanks.

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ I may be misinterpreting your code, but it looks to me as though your bootstrap p-value is calculated by counting the number of bootstrap values > the observed value, which we'd expect to be close to 50%, rather than > (or <) 0. $\endgroup$
    – jbowman
    Commented Sep 1, 2023 at 23:28
  • $\begingroup$ Thanks for your reply, @jbowman. I'm trying to understand why we'd expect it to be close to 50%. $\endgroup$
    – Anthony Nash
    Commented Sep 2, 2023 at 8:49
  • $\begingroup$ 1. We'd expect the observed value to be near the center of the bootstrapped values because they are based, largely, on the same data, with randomness added in one way or another. Given that we don't expect the randomness to be entirely "one-sided"... 2. Typically a p-value is calculated with respect to a hypothesized parameter value, often (usually) zero. Try replacing observedCoefficient with 0 and see what the result is. $\endgroup$
    – jbowman
    Commented Sep 2, 2023 at 15:08

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.