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 <- round(runif(n,4,m))
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!
