doc update: vignettes: add example for constant ir with scalar baseline

vignettes/survival-data-with-a-constant-baseline-hazard.Rmd

@@ -22,7 +22,7 @@ We simulate 100 participants, all entering at age 50 and administratively censor
 
 ```r
 set.seed(0)
-n <- 100
+n <- 1000
 entry_age <- 50
 censor_age <- 70
 beta <- c(0.5, 0.5)
@@ -35,39 +35,99 @@ sim <- HDccAnalysis::sim_cox_age_data(
   censor_age = censor_age,
   beta = beta,
   covariates = covariates,
-  baseline_hazard = 0.005
+  baseline_hazard = 0.01
 )
 
 head(sim)
 # # A tibble: 6 × 6
 #      id entry_age censor_age event   Exp1    Exp2
 #   <int>     <dbl>      <dbl> <int>  <dbl>   <dbl>
-# 1     1        50         70     0  1.26   0.782 
-# 2     2        50         70     0 -0.326 -0.777 
-# 3     3        50         70     0  1.33  -0.616 
-# 4     4        50         70     0  1.27   0.0466
-# 5     5        50         70     0  0.415 -1.13  
-# 6     6        50         70     0 -1.54   0.577 
+# 1     1        50       70       0  1.26  -0.287 
+# 2     2        50       70       0 -0.326  1.84  
+# 3     3        50       50.5     1  1.33  -0.157 
+# 4     4        50       70       0  1.27  -1.39  
+# 5     5        50       70       0  0.415 -1.47  
+# 6     6        50       70       0 -1.54  -0.0695
 ```
 
 ## Quick checks
 
 ```r
 fit <- survival::coxph(survival::Surv(entry_age, censor_age,event) ~ Exp1 + Exp2, data = sim)
 print(fit)
+# Call:
 # survival::coxph(formula = survival::Surv(entry_age, censor_age, 
 #     event) ~ Exp1 + Exp2, data = sim)
 
-#        coef exp(coef) se(coef)     z       p
-# Exp1 0.5133    1.6708   0.2475 2.074 0.03806
-# Exp2 0.6465    1.9088   0.2323 2.783 0.00539
+#         coef exp(coef) se(coef)     z        p
+# Exp1 0.50790   1.66180  0.06965 7.293 3.04e-13
+# Exp2 0.61586   1.85124  0.06960 8.849  < 2e-16
 
-# Likelihood ratio test=13.06  on 2 df, p=0.001456
-# n= 100, number of events= 20 
+# Likelihood ratio test=128  on 2 df, p=< 2.2e-16
+# n= 1000, number of events= 207 
 ```
 
-## Next steps
+The achieved incidence rate is stored as an attribute (events / person-time):
 
-- Replace the constant hazard with `baseline_hazard_by_age` for piecewise rates.
-- Tune `target_avg_baseline_hazard` to hit a desired marginal incidence.
+```r
+attr(sim, "achieved_incidence_rate")
+# [1] 0.01163777
+```
+
+## adjust the baselinehazard to aquire a target incidence rate
+
+We can rerun the previous simulation with the `target_avg_baseline_hazard` to obtain a required incidence rate. This feature is mainly used when you wish to vary beta but adjust teh baseline hazard to maintain the incidence rate:
+
+```r
+set.seed(0)
+sim <- HDccAnalysis::sim_cox_age_data(
+  n = n,
+  entry_age = entry_age,
+  censor_age = censor_age,
+  beta = beta,
+  covariates = covariates,
+  baseline_hazard = 0.01,
+  target_avg_baseline_hazard = 0.02
+)
+
+fit <- survival::coxph(survival::Surv(entry_age, censor_age,event) ~ Exp1 + Exp2, data = sim)
+print(fit)
+# Call:
+# survival::coxph(formula = survival::Surv(entry_age, censor_age, 
+#     event) ~ Exp1 + Exp2, data = sim)
+
+#         coef exp(coef) se(coef)      z      p
+# Exp1 0.53005   1.69902  0.05563  9.529 <2e-16
+# Exp2 0.57055   1.76924  0.05447 10.475 <2e-16
+
+# Likelihood ratio test=189.6  on 2 df, p=< 2.2e-16
+# n= 1000, number of events= 336 
+print(attr(sim, "achieved_incidence_rate"))
+# [1] 0.02058137
+
+beta_2 <- c(1.5, 1.5)
+sim_2 <- HDccAnalysis::sim_cox_age_data(
+  n = n,
+  entry_age = entry_age,
+  censor_age = censor_age,
+  beta = beta_2,
+  covariates = covariates,
+  baseline_hazard = 0.01,
+  target_avg_baseline_hazard = 0.02
+)
+
+fit_2 <- survival::coxph(survival::Surv(entry_age, censor_age,event) ~ Exp1 + Exp2, data = sim_2)
+print(fit_2)
+# Call:
+# survival::coxph(formula = survival::Surv(entry_age, censor_age, 
+#     event) ~ Exp1 + Exp2, data = sim_2)
+
+#         coef exp(coef) se(coef)     z      p
+# Exp1 1.47632   4.37679  0.07301 20.22 <2e-16
+# Exp2 1.46515   4.32819  0.07185 20.39 <2e-16
 
+# Likelihood ratio test=739.5  on 2 df, p=< 2.2e-16
+# n= 1000, number of events= 323 
+print(attr(sim_2, "achieved_incidence_rate"))
+# [1] 0.02068137
+```