Simulating time-fixed confounder in time-varying survival model

Question

I am attempting to simulate a survival Cox PH setting with time-varying exposure in R. The goal is assess the effect of a time-fixed confounder in the relationship exposure-outcome.

What I am trying to do with the following code is to define a longitudinal dataset that can be partially explained by a group variable (grp), the confounder. With these dat I then simulate a time-to-event dataset also partially explained by the variable grp.

Problem: The grp variable does not seem to exert any confounding effect and I don't understand why. From what I know its exclusion from the model should lead to biased estimates of the exposure coefficient, but it does not seem to have that influence.

My attempt:

Simulation of longitudinal data (based on previous question)

###############################################################################
#SIMULATE SURVIVAL DATA WITH TIME_VARYING CONTINUOUS EXPOSURE
##############################################################################

set.seed(1234)


#Load packages
library(MASS)
library(nlme)
library(survival)
library(PermAlgo)
library(data.table)
library(fastDummies)

### set number of individuals
n <- 10000

### average intercept and slope
beta0 <- 10.0
beta1 <- -2.0

### true autocorrelation
ar.val <- -.8

### true error SD, intercept SD, slope SD, and intercept-slope cor
sigma <- 1.5
tau0  <- 2.5
tau1  <- 2.0
tau01 <- 0.3

### maximum number of possible observations (e.g. number of years)
m <- 10

### simulate number of observations for each individual
p <- rep(m,n)

### simulate observation moments (assume everybody has 10 obs)
obs<-rep(1:m, n)

### set up data frame with id and binary time-fixed covariate X1
dat <- data.frame(id=rep(1:n, times=p), obs=obs,
                  x1=rep(rbinom(n, 1, 0.3), each=m))

### simulate confounder as a group variable with 4 levels and corresponding values (grpnum)
grp<-c(rep(1,n*0.2),
       rep(2,n*0.3),
       rep(3,n*0.4),
       rep(4,n*0.1))
dat$grp<-as.factor(grp)

grpnum<-c(rnorm(n = n*0.2,mean = 0,sd = 3),
          rnorm(n = n*0.3,mean=3,sd = 1),
          rnorm(n = n*0.4,mean=-1,sd = 1),
          rnorm(n = n*0.1,mean=-2,sd = 4))
dat$grpnum<-grpnum

### simulate (correlated) random effects for intercepts and slopes
mu  <- c(0,0)
S   <- matrix(c(1, tau01, tau01, 1), nrow=2)
tau <- c(tau0, tau1)
S   <- diag(tau) %*% S %*% diag(tau)
U   <- mvrnorm(n, mu=mu, Sigma=S)

### simulate AR(1) errors and then the actual outcomes
dat$eij <- c(sapply(p, function(x) arima.sim(model=list(ar=ar.val), n=x) * sqrt(1-ar.val^2) * sigma))
dat$yij <- round((beta0 + rep(U[,1], times=p)) +
                (beta1 + rep(U[,2], times=p))  * log(obs) * dat$grpnum  + dat$eij ,2)

### note: use arima.sim(model=list(ar=ar.val), n=x) * sqrt(1-ar.val^2) * sigma
### construction, so that the true error SD is equal to sigma

#Define longitudinal data
dat<-as.matrix(dat[c("x1","grp","yij")])

#Define dummy variables from "grp" column, becaues Permalgo does not accept factors
dat <- dummy_cols(dat,select_columns = "grp",remove_selected_columns = T)
dat <- as.matrix(sapply(dat, as.numeric)) 

Simulate event times conditional on three covariates

# generate vectors of event and censoring times prior to calling the
# function for the algorithm
eventRandom <- round(runif(n, 1,10),0)
censorRandom <- round(runif(n, 1,10),0)

# Generate the survival data conditional on the three covariates
data <- permalgorithm(n, m, dat[,1:5], 
                      XmatNames=c("sex", "exposure","grp1","grp2","grp3"), 
                      eventRandom=eventRandom,
                      censorRandom=censorRandom,
                      betas=c(log(1.5), log(1.3), 
                              log(0.8), log(0.8),
                              log(1.5)))


#Run model

#Including grp
coxph(Surv(Start, Stop, Event) ~ sex + exposure + grp1 + grp2 + grp3 , data=data)

#Excluding grp
coxph(Surv(Start, Stop, Event) ~ sex + exposure , data=data)

#The exposure coefficient does not change!

Answer 1

As per my comment to the question, you have these 2 lines:

p <- round(runif(n, 4, m)) 
p <- rep(m, n)

which seems likely to be an oversight, and not the root cause of the problem.

The first thing I note is:

#The exposure coefficient does not change!

But it does change. Albeit not by a great deal but it does change by around 4%.

> #Including grp
> coxph(Surv(Start, Stop, Event) ~ sex + exposure + grp1 + grp2 + grp3 , data=data)
Call:
coxph(formula = Surv(Start, Stop, Event) ~ sex + exposure + grp1 + 
    grp2 + grp3, data = data)

              coef exp(coef)  se(coef)       z        p
sex       0.244747  1.277298  0.028769   8.507  < 2e-16
exposure  0.070450  1.072991  0.000845  83.375  < 2e-16
grp1     -0.048445  0.952710  0.051892  -0.934    0.351
grp2     -0.660520  0.516582  0.059480 -11.105  < 2e-16
grp3      0.343343  1.409652  0.048284   7.111 1.15e-12

Likelihood ratio test=6292  on 5 df, p=< 2.2e-16
n= 39889, number of events= 5556 
> #Excluding grp
> coxph(Surv(Start, Stop, Event) ~ sex + exposure , data=data)
Call:
coxph(formula = Surv(Start, Stop, Event) ~ sex + exposure, data = data)

              coef exp(coef)  se(coef)       z        p
sex      0.2263063 1.2539597 0.0287483   7.872 3.49e-15
exposure 0.0731306 1.0758711 0.0006808 107.422  < 2e-16

Likelihood ratio test=5598  on 2 df, p=< 2.2e-16
n= 39889, number of events= 5556 

So we have a change of (0.0731306 - 0.070450) / 0.070450 = 3.8%

To investigate a little further I changed the dependency of both the outcome and exposure on grp, by changing this:

# dat$yij <- round((beta0 + rep(U[,1], times=p)) +
                   #(beta1 + rep(U[,2], times=p))  * log(obs) * dat$grpnum  + dat$eij ,2)

to this:

# Debugging
dat$yij <- round((beta0 + rep(U[,1], times=p)) +
                   (beta1 + rep(U[,2], times=p))  * log(obs) * 5 * dat$grpnum  + dat$eij, 2)
# End of Debugging

and adding this line:

# Debugging
data$exposure <- data$exposure + 2 * (data$grp1 + data$grp2 + data$grp3)
# End of Debugging

immediately after the call to permalgorithm

These changes resulted in this output:

> coxph(Surv(Start, Stop, Event) ~ sex + exposure + grp1 + grp2 + grp3 , data=data)
Call:
coxph(formula = Surv(Start, Stop, Event) ~ sex + exposure + grp1 + 
    grp2 + grp3, data = data)

               coef  exp(coef)   se(coef)       z        p
sex       0.1227734  1.1306281  0.0290448   4.227 2.37e-05
exposure  0.0136026  1.0136955  0.0001807  75.275  < 2e-16
grp1      0.1108150  1.1171882  0.0517749   2.140   0.0323
grp2     -0.9389425  0.3910411  0.0613747 -15.299  < 2e-16
grp3      0.2414734  1.2731236  0.0495794   4.870 1.11e-06

Likelihood ratio test=6724  on 5 df, p=< 2.2e-16
n= 39889, number of events= 5556 
> #Excluding grp
> coxph(Surv(Start, Stop, Event) ~ sex + exposure , data=data)
Call:
coxph(formula = Surv(Start, Stop, Event) ~ sex + exposure, data = data)

              coef exp(coef)  se(coef)       z        p
sex      0.1227088 1.1305552 0.0290455   4.225 2.39e-05
exposure 0.0149415 1.0150537 0.0001402 106.601  < 2e-16

Likelihood ratio test=5831  on 2 df, p=< 2.2e-16
n= 39889, number of events= 5556 

And now the change in the estimate for the exposure is:

(0.0149415 - 0.0136026) / 0.0149415 = 9%

So the change is larger, as we would expect.

So to sum up, I can't find anything (obviously) wrong with your code apart from oversight about p