Data Generation for Landmark Survival Models

Awan Afiaz

1 Data Generation

1.1 Source

1.2 Model and Notations

This paper studies the following landmark Cox model by Li et al 2007:

  • Let \(R(s)\) be the population of subjects at risk at time \(s\), i.e., subjects with lifetime \(T > s\). The landmark Cox model specifies that, for a subject in risk set \(R(s)\), the conditional distribution of the residual survival time \(T(s) = T − s\) given predictors \(X(s)\) follows a Cox model. Its survival function, evaluated at time \(s+u \ (u ≥ 0)\) and conditional on the covariates at time \(s\), is given by: \[ \begin{equation} \begin{split} \Pr[T(s) > u \vert \mathbf{X}(s)] & = -\mathbf{X}(s)^{T}\mathbf{\beta}(s) + \epsilon(s) \\ & = -\log(\theta(s) + \epsilon(s) \end{split} \tag{19.2} \end{equation} \]

  • Here \(s\) is called the landmark time, which is the time that has elapsed since the baseline.

  • The time \(u\) denotes the time elapsed since \(s\). For a given \(s\), model (1) is a Cox model with time-independent covariates in the sense that \(\mathbf{X}(s)\) does not vary with \(u\). But \(\mathbf{X}(s)\) may vary with the landmark time \(s\). All the parameters of this model, including the baseline hazard function \(h_0(u; s)\) and the log hazard ratio \(\mathbf{\beta}(s)\), may vary with the landmark time \(s\).

  • The paper derives that there exists a joint distribution of longitudinal and survival data that satisfies the landmark Cox model (19.2) without additional restrictions other than the model assumptions themselves and provides an algorithm to generate data from this joint distribution.

1.3 Joint Distribution of the Longitudinal and Time-to-Event Data for the Landmark Cox Model

  • According to model (19.2), the Cox model at landmark time \(s\) can be expressed as: \[ \begin{equation} G_0(T(s); s) = \exp[-\exp{\mathbf{X}(s)^{T}\mathbf{\beta}(s)} \int_0^u h_0(v, s)dv], (u \ge 0), \end{equation} \]

  • where \(G_0(T(s); s) = \log{H_0(t ; s)}\) and \(\epsilon(s)\) has a standard extreme value distribution with density function \(f(x) = \exp(x − e^x )\), \(x \in (-\infty, \infty)\) Cheng et al. 1995.

  • We first consider the situation where the exponentiated linear predictor \(\theta(s) = \exp{\{\mathbf{X}(s)^{T}\mathbf{\beta}(s)\}}\) among the subjects in \(R(s)\) follows a Gamma distribution with shape parameter \(\alpha(s)\) and rate parameter \(\eta(s)\).

  • Under this situation, the marginal survival function of \(T(s)\) can be calculated as \[ \begin{equation} \Pr[T(s) \gt u] = \left[ \frac{\eta(s)}{H_0(u;s) + \eta{(s)}} \right]^{\alpha(s)}, \forall u \ge 0 \tag{19.3} \end{equation} \]

  • The conditional distribution of \(\theta(s)\) given \(W(s) = G_0(T(s); s)\) is \[ \begin{equation} \begin{split} f(\theta(s) \vert W(s)) & \propto f(W(s) \vert \theta(s)) \ f(\theta(s)) \\ & \propto \theta(s)^{(\alpha(s) + 1) - 1} \exp{\{-(\eta(s) + e^{W(s)} \theta(s))\}} \end{split} \tag{19.4} \end{equation} \]

  • This is a gamma distribution with shape parameter \(\alpha(s) + 1\) and rate parameter \(\eta(s) + e^{W(s)} = \eta(s) + H_0(T(s); s)\).

  • A key observation is that the following equality (the memoryless property) holds here: \[ \begin{equation} \Pr[T(s) | u] = \Pr[T > s+u \vert T > s], \forall \ u,s \ge 0. \tag{19.5} \end{equation} \]

  • The left-hand side of this equality is the survival probability of the residual survival time among subjects in the risk set \(R(s)\). The conditional probability on the right-hand side is specified for all the subjects in the population. The probabilities on both sides of the equality are marginal probabilities involving survival time \(T\) only, i.e., not conditional on \(\mathbf{X}(.)\).

  • Both (19.3) and (19.5) imply that \[ \begin{equation} \left[ \frac{H_0(u;s) + \eta(s)}{\eta(s)} \right]^{\alpha(s)} = \left[ \frac{H_0(s+u;0) + \eta(0)}{H_0(s;0) + \eta(0)} \right]^{\alpha(0)} \end{equation} \]

  • Hence, given \(\alpha(s)\), \(\eta(s)\) and \(H_0(u; 0)\), the baseline hazard at landmark time \(s\), \(H_0(u; s)\) is uniquely determined by, \[ \begin{equation} H_0(u; s) = \eta(s) \Biggl( \left[ \frac{H_0(s+u;0) + \eta(0)}{H_0(s;0) + \eta(0)} \right]^{\frac{\alpha(0)}{\alpha(s)}} - 1 \Biggr) \tag{19.6} \end{equation} \]

  • It can be shown that \(H_0(u; s)\) is a monotone increasing function in \(u\) at any fixed \(s\) and equals \(0\) when \(u = 0\). Therefore, \(H_0(u; s)\) as given by this equation is a proper cumulative hazard function.

  • The two parameter functions \(\alpha(s)\) and \(\eta(s)\) characterize the distribution of \(\theta(s)\), the exponentiated linear predictor, in the time-varying risk set \(R(s)\). The univariate function \(H_0(t ; 0)\) is the baseline cumulative hazard function of the Cox model at baseline (i.e., \(s = 0\)). Equation \((19.6)\) shows that these three functions uniquely determine the baseline cumulative hazard function of all the subsequent Cox models (i.e., \(\forall s > 0\)).

1.4 Algorithm to Generate Data from the Joint Distribution

  1. Prespecify a time grid for the longitudinal measurement times of all subjects, denoted by \(s_1 = 0 < s_2 <, \ldots,< s_K\).

  2. Prespecify \(\alpha(s)\), \(\eta(s)\), \(\beta(s)\), \(H_0(t; 0)\), and calculate \(H_0(u; s)\) using Eq. \((19.6)\).

  3. For subjects \(i = 1, 2,\ldots, n\),

  1. Simulate \(T_i\) from its marginal distribution. Since \(T_i = T_i(0)\), the survival function of this distribution is given by \((19.3)\), with \(s = 0\). Let \(K_i\) be the number of longitudinal measurement times that fall within \([0, Ti )\). Calculate \(T_i(s) = T_i − s\) with \(s_1 = 0 < s_2 <, \ldots,< s_K\).

  2. We generate \(\{ \theta_i(s); s = s_1, s_2, \ldots,< s_K \}\) from a \(K_i\)-dimensional Gaussian copula distribution with correlation parameter \(\rho\). Based on \((19.4)\), the marginal distribution of \(\theta_i(s_j) (j = 1, 2, \ldots,K_i )\) in the copula is a gamma distribution with shape parameter \(\alpha(s_j) + 1\) and rate parameter \(\eta(s_j) + H_0(T_i(s_j ); s_j )\).

  3. When there is a single covariate in \(\mathbf{X}(s)\), calculate \(X_i(s) = \log{\theta_i}(s)/\beta(s)\) for \(s = s_1, s_2, \ldots, s_{K_i}\). When there are \(M\) covariates \((M > 1)\), write \(\log{\theta_i}(s) = \sum_{m=1}^{M} \beta(s) X_{mi}(s)\). Since \(\theta_i(s)\) is generated in Step \(3(b)\), and all the \(\beta_m(.)\) functions are pre-specified, any time-invariant or time-dependent longitudinal covariate processes \(\{X_{mi}(s); m = 1, 2, \ldots , M \}\) that satisfy the linear constraint above are sufficient.

  1. Generate a random censoring time \(C_i\) for each subject and censor both lifetime \(T_i\) and the longitudinal data process.

1.5 Simulation

We now illustrate the proposed data generating algorithm with the following simulation:

1.5.1 Prespecified values

  • The landmark times are pre-specified at s = 0, 2, 4, 6.
  • The marginal distribution of \(\theta(s)\) is a gamma distribution with \(\alpha(s) = 1 - 0.05 \sqrt{s}\) and \(\eta(s) = 1.5 + 0.05s\).
  • The cumulative hazard function at \(s = 0\) is from a Weibull distribution with \(H_0(u; 0) = (\lambda u)^{\kappa}\), \(\lambda = 0.2\) and \(\kappa = .3\).
  • The administrative end of follow-up time is \(15\).

1.5.2 Procedure

  • Given \(\alpha(s)\), \(\eta(s)\) and \(H_0(u; 0)\), following the algorithm described in Sect. \(19.2\), we first generate the event time \(T\) at \(s = 0\) by \((19.3)\).

  • Based on (19.4), \(\theta(s)\) is simulated from the gamma distribution with shape parameter \(\alpha(s) + 1\) and rate parameter \(\eta(s) + H_0(T(s); s)\)), where \(T(s) = T − s\).

  • We use an exchangeable correlation structure with \(\rho = 0.5\) in the Gaussian copula to impose serial correlation on \(\theta(s)\).

  • We set the regression coefficients for 4 covariates (with \(1\) time fixed variable and \(3\) time-varying) as a function of the landmark time \(s\),

    • \(\beta_1(s) = 0.3 + 0.1s\)
    • \(\beta_2(s) = -0.2 + 0.005 * s\)
    • \(\beta_3(s) = 0.2 + 0.01 * s\) and
    • \(\beta_4(s) = -0.3\)

  • In addition, random drop-out times are generated from a Weibull(shape = \(1.5\), scale = \(13\)) distribution to keep the censoring rate at around \(30\%\).

  • For now, a sample size of \(n = 100\) is examined in this simulation study and \(N = 1000\) data sets are simulated.

  • We fit a separate Cox model at each landmark time on the at-risk subjects, and compare the estimated \(\beta(s)\) with its true value at that landmark time.

2 Simulation

2.1 Setting up values and functions

## Sample sizes
n = c(50, 100, 200, 500)
## Landmark times; used as "s.time" in the function of data.fn. 
s.time = seq(from = 0, to = 6, by = 2)
## Andministrative end of follow up time
K = 15
# the correlation parameter used in AR(1) structure to be used in Gaussian Copula
rho1 = 0.5 

## User defined functions 

# 1. We assume the log of linear predictor theta(s) among the subjects in the Risk set R(s) has a 
# Gamma (alpha_s, eta_s) distribution

alpha.fn <- function(s) {1 - 0.05 * sqrt(s)}
eta.fn <- function(s) {1.5 + 0.05 * s}

# ## Pre-specify the beta values 
beta1.fn <- function(s) {0.3 + 0.1 * sqrt(s)}
beta2.fn <- function(s) {-0.2 + 0.005 * s}
beta3.fn <- function(s) {0.2 + 0.01 * s}
beta4.fn <- function(s) {-0.3}

## The Cumulative Hazard function at s = 0 is Weibull (lambda = 0.2, kappa = 3)
H0fn.T <- function(t) {
  lambda = .2
  kappa = 3
  return( (lambda*t)^kappa )
}

2.2 Data Frame Generation

data.fn <- function(n, K, s.time, 
                    cens = FALSE, cens.scale = NULL, 
                    X2.tfix = FALSE, X3.tvary = FALSE, X4.tvary = FALSE, 
                    beta1.fn = NULL, beta2.fn = NULL, beta3.fn = NULL, beta4.fn = NULL, 
                    alpha.fn, eta.fn, rho1, H0fn.T) {
  
## Part 01: Calculate cumulative baseline Hazard at time T for a given landmark s and theta_{i}(s)
## Using equation 19.6
  
  H0fn.Ts <- function(u,s){
        H0.Ts.value <- eta.fn(s) * (((H0fn.T(u+s) + eta.fn(0))/(H0fn.T(s) + eta.fn(0)))^(alpha.fn(0)/alpha.fn(s)) - 1)
        return(H0.Ts.value)
  }
  
#### Part 02: Generate failure times T_{i} and landmark specific times T_{j}(s) #############
## We use 4 covariates: X1 - X4, with X2 being time fixed binary factor and others are time varying
  
  simData <- data.frame(matrix(0, nrow = n*length(s.time), ncol = 12))
  colnames(simData) <- c("id", "stime", "alpha", "eta", "X1", "X2", "X3", "X4", "beta1", "beta2", "beta3", "beta4")
  
  simData$id <- sort(rep(seq(from = 1, to = n, by = 1), length(s.time)))
  simData$stime <- rep(s.time, n)
  simData$alpha <- rep(alpha.fn(s.time), n)
  simData$eta   <- rep(eta.fn(s.time), n)
  simData$beta1 <- rep(beta1.fn(s.time), n)
  X2.all <-  rbinom(n, size = 1, prob = 0.5)
  C <- rep(K,n)
  if(cens){ # if random censoring is true;
      C <- rweibull(n, shape = 1.5, scale = cens.scale) 
      C[ C > K] <- K 
  }
  S.T <- runif(n, 0, 1)
  for (i in 1:n){
      surv.T <- function(t){
             (eta.fn(0)/(eta.fn(0)+H0fn.T(t)))^alpha.fn(0)-S.T[i] # marginal survival function of T; (Equation 19.3)
       }
      simData$T.real[simData$id == i] <- uniroot(f = surv.T, lower = 1E-6, upper = 10000)$root
      
      simData$C[simData$id == i] <- C[i]
      
      if(X2.tfix){
          simData$X2[simData$id == i] <- X2.all[i]
          simData$beta2[simData$id == i] <- beta2.fn(s.time)
      }
      if(X3.tvary){
          simData$X3[simData$id==i] <- rnorm(length(s.time), mean = -1+2*simData$X2[simData$id == i][1], sd = 1)    
          simData$beta3[simData$id==i] <- beta3.fn(s.time)
      }   
      if(X4.tvary){
          cp.X4 <- copula::normalCopula(param = 0.5, dim = length(s.time), dispstr = "ex")  
          F.X4 <- copula::rCopula(n = 1, copula = cp.X4)
          simData$X4[simData$id == i] <- qnorm(p = F.X4, mean = 0, sd = 1)
          simData$beta4[simData$id == i] <- beta4.fn(s.time)
          rm(cp.X4,F.X4)
      }   
  }

  simData$Ts.real <- simData$T.real - simData$stime  # (uncensored) residual lifetime;
  simData$estate <- 1*(simData$T.real <= simData$C)  ## event status (Failed or Censored)
  
  simData$T <- simData$T.real*simData$estate + simData$C*( 1 - simData$estate)  
  
  simData$Ts <- simData$T - simData$stime   #(censored) residual lifetime after landmark time s;
  simData <- subset(simData, Ts > 0)
  simData <- simData[order(simData$id, simData$stime), ]
  Ts <- simData$Ts[simData$estate == 1]
  ties = length(Ts) - length(unique(Ts))
  
  ## Warning in case of ties (very unlikely, but kept as a check)
  if(ties > 0){
       cat("***** Warning *****: Number of ties among Ts in the sample: ", ties, "\n")
  }
  rm(S.T, C, X2.all, Ts, ties)
  
########### Part 03: Generate theta(s) and X_{i}(s) ###########

  for (i in sort(unique(simData$id))){
      nrow.sub <- nrow(subset(simData,id == i))
      if (nrow.sub>1){
        cp <- copula::normalCopula(param = rho1, dim = nrow.sub, dispstr = "ex")
        u <- copula::rCopula(n = 1, copula = cp)
      } else{
        u <- runif(n = 1)
      }
     
      alpha.s <- simData$alpha[simData$id == i] + 1
      eta.s <- simData$eta[simData$id == i] + H0fn.Ts(simData$Ts.real[simData$id == i], simData$stime[simData$id == i])
      simData$theta[simData$id == i] <- qgamma(p = u, shape = alpha.s, rate = eta.s)
                                            
      rm(u, alpha.s, eta.s)
  }

  simData$X1 <- (log(simData$theta) - simData$beta2*simData$X2 - simData$beta3*simData$X3 - simData$beta4*simData$X4)/simData$beta1
  
  simData <- subset(simData, select = c(id, T, stime, Ts, estate, X1, X2, X3, X4))
  simData <- simData[order(simData$id,simData$stime), ]
  
  # ### Introducing Treatment based on observed X's (to ensure confounding We again use the X's to generate the treatment probabilities)
  # logit.fn <- function(X1, X2, X3, X4) {
  # xb = 0.5*X1 + 0.4*X2 + 1*X3 + 1*X4
  # prob = exp(xb)/(1 + exp(xb))
  # return(round(prob, 2))
  # }
  # simData <- simData %>% 
  #   dplyr::mutate(p = logit.fn(X1, X2, X3, X4), Treat = ifelse(p > 0.5, 1, 0))
  
  return(simData)
}

2.3 Generate one data set

set.seed(420)
simData <- data.fn(n = 100, K = 15, s.time, cens = T, cens.scale = 13, 
                   X2.tfix = T, X3.tvary = T, X4.tvary = T, 
                   beta1.fn = beta1.fn, beta2.fn = beta2.fn, beta3.fn = beta3.fn, 
                   beta4.fn = beta4.fn, alpha.fn, eta.fn, rho1 = 0.5, H0fn.T)

2.4 Checking censoring proportion to see if they are at \(30\%\)

simData %>% group_by(stime) %>% 
  count(estate) %>% 
  mutate(prop = n/sum(n),
         stime = if_else(estate == 1, true = NA_real_, stime)) %>%
  kable(digits = 2, format = "html", booktabs = T, 
        caption = "Censoring proportions at each landmark time (Event = 0 implies censored)",
        col.names = c("Landmark time", "Event status", "Observations", "Proportion"), 
        align = "c", escape = F) %>%
  kable_styling(font_size = 14)
Censoring proportions at each landmark time (Event = 0 implies censored)
Landmark time Event status Observations Proportion
0 0 31 0.31
1 69 0.69
2 0 27 0.29
1 67 0.71
4 0 22 0.33
1 45 0.67
6 0 11 0.30
1 26 0.70

2.5 Run 1000 simulations to check estimated Beta’s

2.5.1 Sample size n = 50

## Placeholders for Beta estimates
B = 1000
s0 = matrix(0, nrow = B, ncol = 4)
s2 = matrix(0, nrow = B, ncol = 4)
s4 = matrix(0, nrow = B, ncol = 4)
s6 = matrix(0, nrow = B, ncol = 4)

## For loop
for(i in 1:B) {
  DF = data.fn(n = 50, K = 15, s.time, cens = T, cens.scale = 13,
               X2.tfix = T, X3.tvary = T, X4.tvary = T, 
               beta1.fn = beta1.fn,beta2.fn = beta2.fn,
               beta3.fn = beta3.fn, beta4.fn = beta4.fn, alpha.fn,
               eta.fn, rho1 = 0.5, H0fn.T)
  ## For s = 0 (at baseline)             
  s0[i, ] = survival::coxph(Surv(time = Ts, event = estate) ~ X1 + X2 + X3 + X4,
                            data = DF %>% filter(stime == 0))$coefficients
  
  ## For s = 2 (First follow up)
  s2[i, ] = survival::coxph(Surv(time = Ts, event = estate) ~ X1 + X2 + X3 + X4,
                            data = DF %>% filter(stime == 2))$coefficients
  
  ## For s = 4 (Second follow up)
  s4[i, ] = survival::coxph(Surv(time = Ts, event = estate) ~ X1 + X2 + X3 + X4,
                            data = DF %>% filter(stime == 4))$coefficients
  
  ## For s = 6 (Third follow up)
  s6[i, ] = survival::coxph(Surv(time = Ts, event = estate) ~ X1 + X2 + X3 + X4,
                            data = DF %>% filter(stime == 6))$coefficients
}
tab <- tibble(betas = c("$\\beta_1$", "$\\beta_2$", "$\\beta_3$", "$\\beta_4$"),
              s0_T = c(beta1.fn(0), beta2.fn(0), beta3.fn(0), beta4.fn(0)), 
              s0 = colMeans(s0, na.rm = T), 
              s2_T = c(beta1.fn(2), beta2.fn(2), beta3.fn(2), beta4.fn(2)), 
              s2 = colMeans(s2, na.rm = T),
              s4_T = c(beta1.fn(4), beta2.fn(4), beta3.fn(4), beta4.fn(4)), 
              s4 = colMeans(s4, na.rm = T),
              s6_T = c(beta1.fn(6), beta2.fn(6), beta3.fn(6), beta4.fn(6)), 
              s6 = colMeans(s6, na.rm = T))

tab %>%
  kable(digits = 3, format = "html", booktabs = T, 
        caption = "Average estimated values against true betas (n = 50)",
        col.names = c("Coeffs", rep(x = c("True", "Estimated"), times = 4)),
        align = "c", escape = F) %>% 
  add_header_above(c(" ", "s = 0" = 2, "s = 2" = 2, "s = 4" = 2, "s = 6" = 2))
Average estimated values against true betas (n = 50)
s = 0
s = 2
s = 4
s = 6
Coeffs True Estimated True Estimated True Estimated True Estimated
\(\beta_1\) 0.3 0.331 0.441 0.487 0.50 0.583 0.545 3.21
\(\beta_2\) -0.2 -0.214 -0.190 -0.200 -0.18 -0.186 -0.170 -4.06
\(\beta_3\) 0.2 0.206 0.220 0.238 0.24 0.257 0.260 2.38
\(\beta_4\) -0.3 -0.349 -0.300 -0.356 -0.30 -0.356 -0.300 -1.43

2.5.2 Sample size n = 100

for(i in 1:B) {
  DF = data.fn(n = 100, K = 15, s.time, cens = T, cens.scale = 13,
               X2.tfix = T, X3.tvary = T, X4.tvary = T, 
               beta1.fn = beta1.fn,beta2.fn = beta2.fn,
               beta3.fn = beta3.fn, beta4.fn = beta4.fn, alpha.fn,
               eta.fn, rho1 = 0.5, H0fn.T)
  ## For s = 0 (at baseline)             
  s0[i, ] = survival::coxph(Surv(time = Ts, event = estate) ~ X1 + X2 + X3 + X4,
                            data = DF %>% filter(stime == 0))$coefficients
  
  ## For s = 2 (First follow up)
  s2[i, ] = survival::coxph(Surv(time = Ts, event = estate) ~ X1 + X2 + X3 + X4,
                            data = DF %>% filter(stime == 2))$coefficients
  
  ## For s = 4 (Second follow up)
  s4[i, ] = survival::coxph(Surv(time = Ts, event = estate) ~ X1 + X2 + X3 + X4,
                            data = DF %>% filter(stime == 4))$coefficients
  
  ## For s = 6 (Third follow up)
  s6[i, ] = survival::coxph(Surv(time = Ts, event = estate) ~ X1 + X2 + X3 + X4,
                            data = DF %>% filter(stime == 6))$coefficients
}
Average estimated values against true betas (n = 100)
s = 0
s = 2
s = 4
s = 6
Coeffs True Estimated True Estimated True Estimated True Estimated
\(\beta_1\) 0.3 0.313 0.441 0.461 0.50 0.540 0.545 0.634
\(\beta_2\) -0.2 -0.218 -0.190 -0.192 -0.18 -0.188 -0.170 -0.199
\(\beta_3\) 0.2 0.207 0.220 0.226 0.24 0.253 0.260 0.300
\(\beta_4\) -0.3 -0.316 -0.300 -0.312 -0.30 -0.327 -0.300 -0.330

2.5.3 Sample size n = 300

for(i in 1:B) {
  DF = data.fn(n = 300, K = 15, s.time, cens = T, cens.scale = 13,
               X2.tfix = T, X3.tvary = T, X4.tvary = T, 
               beta1.fn = beta1.fn,beta2.fn = beta2.fn,
               beta3.fn = beta3.fn, beta4.fn = beta4.fn, alpha.fn,
               eta.fn, rho1 = 0.5, H0fn.T)
  ## For s = 0 (at baseline)             
  s0[i, ] = survival::coxph(Surv(time = Ts, event = estate) ~ X1 + X2 + X3 + X4,
                            data = DF %>% filter(stime == 0))$coefficients
  ## For s = 2 (First follow up)
  s2[i, ] = survival::coxph(Surv(time = Ts, event = estate) ~ X1 + X2 + X3 + X4,
                            data = DF %>% filter(stime == 2))$coefficients
  ## For s = 4 (Second follow up)
  s4[i, ] = survival::coxph(Surv(time = Ts, event = estate) ~ X1 + X2 + X3 + X4,
                            data = DF %>% filter(stime == 4))$coefficients
  ## For s = 6 (Third follow up)
  s6[i, ] = survival::coxph(Surv(time = Ts, event = estate) ~ X1 + X2 + X3 + X4,
                            data = DF %>% filter(stime == 6))$coefficients
}
Average estimated values against true betas (n = 300)
s = 0
s = 2
s = 4
s = 6
Coeffs True Estimated True Estimated True Estimated True Estimated
\(\beta_1\) 0.3 0.304 0.441 0.449 0.50 0.511 0.545 0.567
\(\beta_2\) -0.2 -0.209 -0.190 -0.195 -0.18 -0.190 -0.170 -0.174
\(\beta_3\) 0.2 0.204 0.220 0.224 0.24 0.247 0.260 0.269
\(\beta_4\) -0.3 -0.303 -0.300 -0.301 -0.30 -0.304 -0.300 -0.314
  • The simulated values are very close to the true values at each landmark time point with increasing sample sizes.