A question about a hazard model (discrete, fixed intervals, fixed effect time)

Question

We need to calculate the sample size for our proposal. The appropriate model to apply is the hazard model (discret), with a fixed effect of time, fixed and same time intervals for all subjects.

Our study focuses on patients with neurological diseases, specifically investigating factors related to time-to-dementia. Due to the costs associated with neurological testing, we will assess all patients in the prospective cohort on a monthly basis, up to 18 months, with assessments every three months (T0, T3, T6, T9, T12, T18). Dementia will be assessed at each visit. The study involves a single event with the same interval for all subjects. In this study there are two covariate: age (x1) and assessment result which is a binary time dependent covariate (x2) (with the same interval as event). In each visit we check the status of assesment result.

Now my questions are:

Question 1. I found the following simulation code to simulate a discrete Hazard model (link=logit). How can I determine the correct functions for the baseline hazard and the effect of variables (steps 6 and 7)? How can I add the random variability you would see in practice? How can model the time varying covariate?

Simulation -Discrete Time Hazard model.

Question 2. It seems that there are similarities between a discrete time proportional hazard model (with fixed intervals and fixed time effect) and grouped proportional hazard data. Do you think the sample size formula proposed in this publication can be used in our study?

https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.7847

3. How should we present the data for a discrete time model with a time dependent covariate? For example, if we want to present the record for the second subject (please find the picture) who experienced the event at T6 which presentation is the right ones? Please show me the well structure format for presenting the records of ID=2 (as an example).

time of X1=1 (level of interest in covariate) (filled circle)

Id=2
         time   event    covariate
         T0      0         0  
         T3      0         0
         T6      1         1

or

ID=2 
       time   event    covariate
         T0      0         0  
         T3      0         0
         T6      1         0

Here is the Simulation code (logit function) for the question 1.

###  1. Build some useful functions
# create the logit function, and its inverse (over values of x # from 0–1), the logistic function
 
logistic <- function(x) { return(1/(1+exp(-1*x))) }
logit    <- function(x) { return(log(x/(1-x))) }

### 2. Specify N,the total number of individuals (units) at-risk at T=0
N <- 1000 
   
### 3. Specify T,the maximum value of t in your simulated study.
T <- 20 
  
### 4. Create a data frame called Data starting with the ID variable
ID <- 1:N
Data <- data.frame(ID)

# Create conditioning variables (e.g., two dichotomous vars here, but
# these could be whatever number and kind), and bind these to your 
# data set
var1 <- c(rep(0,(N/2)),rep(1,N/2))
var2 <- c(rep(0,(N/4)),rep(1,(N/2)),rep(0,(N/4)))
Data <- cbind(Data,var1,var2)


### 5. Format Data as a person-period structure

expand <- function(x,t) {
  data <- x[rep(1:length(x[,1]), each = t), ]
  rownames(data) <- NULL
  return(data)
  }

Data <- expand(Data, T)
 Data$period <- rep(1:T,N)
     for (t in 1:T) {
       varname <- paste("t",t,sep="")
       Data[,varname] <- as.integer(0 + Data$period == t)
   }

### 6. Prepare effect of time on discrete-time hazard function
# Based on your above model assumptions about how ht relates to time:

## Constant baseline hazard:
conslogit   <- logit(0.1)   # logit(.1)   = -2.1972246    # for a nominal hazard of 0.1

## Baseline hazard as a linear function of time
linearlogit   <- logit(logistic(0) + 0.05)  # for a linear effect of time on hazard of 0.05 per 1-unit increase in t 
   
### 7. Prepare effects of conditioning variables on ht
var1logit   <- logit(logistic(0) + 0.06)   # for a nominal effect of var1 of 0.06 increase in hazard:
var2logit   <- logit(logistic(0) + 0.04)   # for a nominal hazard effect of var2 of 0.04:
   
8. Simulate your discrete-time survival data

# Example for baseline fully discrete effect of time:
for (t in 1:T) {
  Data$hlogit[Data$period==t] <- logistic(discretelogit[t])
  }
   
# Example for conditional logit hazard model with constant effect of time:
  # (The [Data$period==t] on the conditioning variables are probably only
  # necessary if your conditioning variables are *time varying*: if the values of var1 or var2 were to change over time *within individuals*

Here are the results of a similar study that we want to consider as input for calculating the sample size.

-- assesment result: HR=4
-- Age : HR = 1
--  Event  vs X2   
                                    X2=1                   X2=0
Event=1 : Dementia                  n=7                    n=10
Event=0 : Censor                    n=4                    n=6
Note: X2=1 is the level of interest 

Answer 1

I very strongly recommend that you work with an experienced local statistician to develop your proposal. I'm pretty good at survival analysis; nevertheless, my group still relies on our colleagues in biostatistics to refine designs and to write up concise and compelling power analyses for project proposals. As you need to understand the underlying issues well enough to have fruitful discussions with the statistician, here's an outline of what you need to consider.

First, you have to base the power analysis on considerations that are specific to your own situation. An earlier iteration of this question indicates that the other study you are using as a guide is J Neurooncol. 2019; 144(3): 511–518. There is no reason to suppose that study, of whether a decline in cognitive function is associated with cancer disease progression, will have anything to inform your study of how some time-varying covariate is associated with time-to-dementia.*

Second, most power calculations for survival analysis, like those in a paper you link, are based on treatment or biomarker groups that are defined at the start of the study. The time-varying binary covariate in your study adds further complication. Based on your and your colleagues' understanding of the subject matter, you have to estimate the rate of progression to dementia, how that might be associated with a change in that covariate, and the characteristics of how that covariate itself changes over time. To start:

• What fraction of individuals in your study might develop dementia within the first year if they have X2 = 0 throughout?
• What fraction would develop dementia if they have X2 = 1 throughout?

If you are willing to assume a constant baseline hazard over time, the corresponding exponential model for survival over time is: $S(t)=\exp (-\lambda t)$, with $\lambda$ the hazard rate per unit time. If you use 1-year survival and rates per year, then for each of the above scenarios you can estimate $\lambda = -\log S(t=1)$, with the ratio of the two $\lambda$ values being your hazard-ratio estimate.

That still leaves the question of the time course of the change in your binary covariate. Based on the J. Neurooncol. paper, I take it that once the binary covariate X2 moves from a value of 0 to 1 it stays at the value of 1. That simplifies things somewhat, as you can think about this as a multi-state model instead of treating X2 as a covariate for the simulation. The 3 states would be X2 = 0, X2 = 1, and dementia.

Possible transitions are from either X2 state to dementia, and a transition from X2 = 0 to X2 = 1. Under the constant baseline hazard assumption, you already have estimates for the rates of developing dementia starting from either of the X2 = 0 or the X2 = 1 states. You just need to estimate how the transition from X2 = 0 to X2 = 1 occurs over time. For example, you could treat the X2 = 0 to X2 = 1 transition as another exponential survival model with a rate constant that you estimate based on your understanding of the subject matter.

Third, once you have decided on all the rates, use them to simulate a large amount of data in continuous time to represent your population of interest. With exponential models that could be easy: use the rexp() function with the 3 corresponding rates in parallel for a very large number of individuals (a few tens of thousands).** Then evaluate the simulated event times for each individual:

If the X2 = 0 to X2 = 1 transition occurs after the X2 = 0 to dementia transition, ignore the simulated X2 = 0 to X2 = 1 transition time.

If the X2 = 0 to X2 = 1 transition occurs before the X2 = 0 to dementia transition, note the time of the X2 = 0 to X2 = 1 transition. Then add to that time the simulated time of that individual's X2 = 1 to dementia transition, for the ultimate simulated time of the transition to dementia. You can get away with that here, as an exponential model has no memory of prior experience.

Fourth, sample repeatedly (1000 or more times) from your simulated population exactly as you expect to draw from the real population of your study. That includes a defined sample size, its accrual over time, dropouts, binning your continuous-time survival times and X2 = 0 to X2 = 1 transition times according to your planned observation intervals, censoring at the end of the study, etc. For each sample, perform your discrete-time analysis with X2 as a time-varying covariate. If you want X2 to be a predictor of transition to dementia, make sure that for each time interval you use its value at the end of the preceding time interval. The fraction of samples that provide a "significant" result is the power for that particular study design.

Fifth, as the initial power estimate will probably not be what you need, repeat step Four until you find a useful sample size.

In response to comments

With respect to the form of baseline hazard, the simplest thing to do is to simulate data as above based on whatever baseline hazard you have in mind. I can't predict what will change if you do so. Again, for project planning it's important to use a baseline hazard that fairly represents the situation with your own study.

You are certainly correct that "discrete-time modeling of events during a time interval works with the covariate values in place during that interval." This raises an important distinction made by Tutz and Schmid in Modeling Discrete Time-to-Event Data, page 5:

discrete time-to-event data occur as

• intrinsically discrete measurements, where the measurements represent natural numbers, or
• grouped data, which represent events in underlying time intervals, and the response refers to an interval.

Although methods are generally the same either way, that's not the case with time-varying covariates. For intrinsically discrete measurements, the time-varying covariate values to use are those in place at observation number $t$, but (page 59, Section 3.5)

if “discrete time” refers to intervals, [use the values] at the beginning of the interval $[a_{t-1},a_t)$. (Emphasis added.]

where $a_t$ represents the actual time of the $t^{th}$ observation. In practice, the covariate value at the beginning of the time interval $[a_{t-1},a_t)$ is what you saw at observation t-1, not the value at observation t.

Think about it this way. Say that there is an event, with X2 = 1, at observation t. Say that at observation t-1 you had X2 = 0 and no event. Both changed during the interval. From the perspective of prediction (let alone causal inference), you can't tell from the grouped data whether the change in X2 caused the event, or whether the event itself led to the change in X2. For the event at observation t, you need to use the X2 value in place at the beginning of the interval, which is what you found at observation t-1.

For the data shown in the picture (with ID 2 presumably the second from the bottom), that means the correct data structure is the second that you show:

ID=2 
       time   event    covariate
         T0      0         0  
         T3      0         0
         T6      1         0

Although a covariate change was seen at T6, you don't know whether that occurred before or after the disease-progression event that was also observed at T6. In this case, it's quite possible that disease progression could lead to cognitive decline. For the covariate to be a "predictor" of the event, you have to use the covariate value in place at the start of the T3-T6 time interval, which is the value observed at T3. The authors of that study evidently didn't do that coding properly, leading to the unrealistically high hazard ratio.*

---

*That study also has a fundamental flaw: if both cognitive decline and cancer progression were first found at the same visit, the way the study coded the discrete-time data nevertheless treated that cognitive decline as a predictor of the progression. Logically, a predictor needs to be in place before the event it's predicting. That logical flaw is probably why the hazard ratio estimate in that study is so unusually high. I'd be reluctant to base any clinical study design on a hazard ratio greater than 2 unless there is extremely compelling data supporting a higher value. If your plans depend on having a very high hazard ratio, you are likely to end up very disappointed when you don't get a clean result from your study because the true hazard ratio isn't quite as high as you had hoped.

**This assumes that the 3 possible transition times within an individual are independent. If not, for example if the X2 = 0 to X2 = 1 transition time has some distribution around the ultimate transition time to dementia, you need to modify the simulation accordingly.