-
Notifications
You must be signed in to change notification settings - Fork 0
/
simulation_recurent.Rmd
145 lines (124 loc) · 4.07 KB
/
simulation_recurent.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
---
title: "Simulating trials with multiple events"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Simulating trials with multiple events}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
bibliography: references.bib
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
```{r setup, message=FALSE, warning=FALSE}
library(survobj)
library(survival)
library(dplyr)
```
## Introduction
This document presents a way to simulate a trial with multiple events and compare the empirical power of the analysis of the first or only episode using Cox regression and the analysis of multiple episodes under an Andersen and Gill model. Data is generated using a renewal homogeneous Poisson process as described @leemis1987 .
## Trial simulation
A total of 1000 trials are simulated using a survival object of class Weibull
with shape of 0.5 and failure rate at time 1 of 40%. Each group will include
250 participants, and the hazard ratio (HR) for the intervention group will be
0.7 and the follow-up will be censored at time 1. Empirical power is defined as the
proportion of trials with a robust p-value below 0.05
```{r, eval= FALSE}
nsim = 1000
s_obj = s_weibull(fail = 0.4, t=1, shape = 0.5)
n = 250
subjid = seq(1, 2*n)
group = c(rep(0,n), rep(1,n))
hr = c(rep(1,n), rep(0.7,n))
tmax = 1
set.seed = 12345
sim <- lapply(
1:nsim,
function(x){
# simulate survival times for one trial
tsim <- matrix(rsurvhr(s_obj, hr), ncol = 1)
i = 1
while(min(tsim[,i]) < tmax) {
i = i+1
tsim<- cbind(tsim,renewhr(s_obj, hr, tsim[,i-1]))
}
# Analysis data.frame
df <- data.frame(
subjid = rep(subjid,i),
group = rep(group, i),
time = as.vector(tsim)
) |>
arrange(subjid, time) |>
group_by(subjid) |>
mutate(ncase = row_number()) |>
mutate(start = lag(time, default = 0)) |>
filter(start < tmax) |>
mutate(event = censor_event(tmax, time, 1)) |>
mutate(end = censor_time(tmax, time))
# Analysis multiple episodes
mult <- summary(
coxph(Surv(start,end,event)~group,
method = "breslow",
id = subjid,
robust = T,
data = df,
control = coxph.control(timefix = FALSE)))
# Analysis first or only episode
sing <- summary(
coxph(Surv(start,end,event)~group,
method = "breslow",
id = subjid,
robust = T,
data = filter(df, ncase == 1),
control = coxph.control(timefix = FALSE)))
# Export results for analysis
return(
data.frame(
simid = c(x,x),
res = c("recurrent","single"),
events = c(mult$nevent, sing$nevent),
hr = c(mult$coefficients[1,"exp(coef)"],sing$coefficients[1,"exp(coef)"]),
pvalue = c(mult$coefficients[1,"Pr(>|z|)"],sing$coefficients[1,"Pr(>|z|)"])
)
)
}
)
# Join all the simulations in a single data frame
sim_rec <- do.call(rbind, sim)
```
```{r, include=FALSE }
# The simulation takes to much time to be included in CRAN
# Load a previous simulation
load("sim_rec.rda")
```
## Analysis of the simulation
### For recurrent events
```{r}
rec_empirical_power =
binom.test(
sum(sim_rec$pvalue[sim_rec$res == "recurrent"] <= 0.05),
length(sim_rec$pval[sim_rec$res == "recurrent"] ))
rec_empirical_power$estimate
rec_empirical_power$conf.int
# Distribution of the simulated VEs}
summary(sim_rec$hr[sim_rec$res == "recurrent"])
# Distribution of the simulated number of events
summary(sim_rec$events[sim_rec$res == "recurrent"])
```
### For first or only event
```{r}
sing_empirical_power =
binom.test(
sum(sim_rec$pvalue[sim_rec$res == "single"] <= 0.05),
length(sim_rec$pval[sim_rec$res == "single"] ))
sing_empirical_power$estimate
sing_empirical_power$conf.int
# Distribution of the simulated VEs}
summary(sim_rec$hr[sim_rec$res == "single"])
# Distribution of the simulated number of events
summary(sim_rec$events[sim_rec$res == "single"])
```
## References