diff --git a/apple-touch-icon-120x120.png b/apple-touch-icon-120x120.png index 07624d2d..0ec52c21 100644 Binary files a/apple-touch-icon-120x120.png and b/apple-touch-icon-120x120.png differ diff --git a/apple-touch-icon-152x152.png b/apple-touch-icon-152x152.png index 7a6cfe88..4de9048c 100644 Binary files a/apple-touch-icon-152x152.png and b/apple-touch-icon-152x152.png differ diff --git a/apple-touch-icon-180x180.png b/apple-touch-icon-180x180.png index 2da753ab..b8203993 100644 Binary files a/apple-touch-icon-180x180.png and b/apple-touch-icon-180x180.png differ diff --git a/apple-touch-icon-60x60.png b/apple-touch-icon-60x60.png index 56f0bad2..24f08533 100644 Binary files a/apple-touch-icon-60x60.png and b/apple-touch-icon-60x60.png differ diff --git a/apple-touch-icon-76x76.png b/apple-touch-icon-76x76.png index 0ecee3fc..c8993701 100644 Binary files a/apple-touch-icon-76x76.png and b/apple-touch-icon-76x76.png differ diff --git a/apple-touch-icon.png b/apple-touch-icon.png index 2da753ab..dfba1f9b 100644 Binary files a/apple-touch-icon.png and b/apple-touch-icon.png differ diff --git a/articles/basics.html b/articles/basics.html index a1254bec..bb71a1c8 100644 --- a/articles/basics.html +++ b/articles/basics.html @@ -190,10 +190,8 @@

Subset with unique event timeslibrary(pammtools) library(dplyr) -data("veteran", package = "survival") # load veteran data -
## Warning in data("veteran", package = "survival"): data set 'veteran' not found
-
-# remove ties to illustrate equivalence with Cox approach
+veteran <- survival::veteran
+# remove ties to illustrate equivalence with Cox approach
 vetu <- filter(veteran, !duplicated(time))
 ped_vetu <- vetu %>%
   as_ped(Surv(time, status)~., cut = unique(vetu$time), id = "id")
@@ -212,7 +210,7 @@ 

Subset with unique event times

Full data set

-
+
 ## Using the full data set (with ties) yields slightly different results
 # when comparing PEM to Cox-PH
 ped_vet <- veteran %>%
@@ -251,7 +249,7 @@ 

Piece-wise exponential additive m performed semi-parametrically. In our experience the differences are negligible.

Using the example above (data without ties) we get:

-
+
 ## compare to PAM
 pam_age <- gam(ped_status ~ s(tend) + age, data = ped_vetu,
     family = "poisson", offset = offset)
diff --git a/articles/competing-risks.html b/articles/competing-risks.html
index 67b066b6..4d682b3f 100644
--- a/articles/competing-risks.html
+++ b/articles/competing-risks.html
@@ -218,15 +218,15 @@ 

Cause specific haza as_ped(Surv(time, status)~., id = "id", cut = cut, combine = FALSE) str(ped, 1)

## List of 2
-##  $ cause = 1:Classes 'ped' and 'data.frame': 29807 obs. of  9 variables:
-##   ..- attr(*, "breaks")= num [1:94] 0.192 0.274 0.312 0.32 0.329 ...
+##  $ cause = 1:Classes 'ped' and 'data.frame': 30438 obs. of  9 variables:
+##   ..- attr(*, "breaks")= num [1:96] 0.0684 0.0767 0.1752 0.1834 0.1944 ...
 ##   ..- attr(*, "id_var")= chr "id"
 ##   ..- attr(*, "intvars")= chr [1:6] "id" "tstart" "tend" "interval" ...
 ##   ..- attr(*, "trafo_args")=List of 3
 ##   ..- attr(*, "time_var")= chr "time"
 ##   ..- attr(*, "status_var")= chr "status"
-##  $ cause = 2:Classes 'ped' and 'data.frame': 29807 obs. of  9 variables:
-##   ..- attr(*, "breaks")= num [1:94] 0.192 0.274 0.312 0.32 0.329 ...
+##  $ cause = 2:Classes 'ped' and 'data.frame': 30438 obs. of  9 variables:
+##   ..- attr(*, "breaks")= num [1:96] 0.0684 0.0767 0.1752 0.1834 0.1944 ...
 ##   ..- attr(*, "id_var")= chr "id"
 ##   ..- attr(*, "intvars")= chr [1:6] "id" "tstart" "tend" "interval" ...
 ##   ..- attr(*, "trafo_args")=List of 3
@@ -238,13 +238,13 @@ 

Cause specific haza
 # data set for event type 1 (death from cardiovascular events)
 head(ped[[1]])
-
##     id    tstart      tend                    interval    offset ped_status
-## 1 5002 0.0000000 0.1916496            (0,0.1916495551] -1.652087          0
-## 2 5002 0.1916496 0.2737851 (0.1916495551,0.2737850787] -2.499385          0
-## 3 5002 0.2737851 0.3121150 (0.2737850787,0.3121149897] -3.261525          0
-## 4 5002 0.3121150 0.3203285 (0.3121149897,0.3203285421] -4.801970          0
-## 5 5002 0.3203285 0.3285421 (0.3203285421,0.3285420945] -4.801970          0
-## 6 5002 0.3285421 0.3969884 (0.3285420945,0.3969883641] -2.681706          0
+
##     id     tstart       tend                    interval    offset ped_status
+## 1 5002 0.00000000 0.06844627            (0,0.0684462697] -2.681706          0
+## 2 5002 0.06844627 0.07665982  (0.0684462697,0.076659822] -4.801970          0
+## 3 5002 0.07665982 0.17522245  (0.076659822,0.1752224504] -2.317063          0
+## 4 5002 0.17522245 0.18343600 (0.1752224504,0.1834360027] -4.801970          0
+## 5 5002 0.18343600 0.19438741 (0.1834360027,0.1943874059] -4.514288          0
+## 6 5002 0.19438741 0.26009582 (0.1943874059,0.2600958248] -2.722528          0
 ##    sex age cause
 ## 1 Male  60     1
 ## 2 Male  60     1
@@ -255,13 +255,13 @@ 

Cause specific haza
 # data set for event type 2 (death from other causes)
 head(ped[[2]])
-
##     id    tstart      tend                    interval    offset ped_status
-## 1 5002 0.0000000 0.1916496            (0,0.1916495551] -1.652087          0
-## 2 5002 0.1916496 0.2737851 (0.1916495551,0.2737850787] -2.499385          0
-## 3 5002 0.2737851 0.3121150 (0.2737850787,0.3121149897] -3.261525          0
-## 4 5002 0.3121150 0.3203285 (0.3121149897,0.3203285421] -4.801970          0
-## 5 5002 0.3203285 0.3285421 (0.3203285421,0.3285420945] -4.801970          0
-## 6 5002 0.3285421 0.3969884 (0.3285420945,0.3969883641] -2.681706          0
+
##     id     tstart       tend                    interval    offset ped_status
+## 1 5002 0.00000000 0.06844627            (0,0.0684462697] -2.681706          0
+## 2 5002 0.06844627 0.07665982  (0.0684462697,0.076659822] -4.801970          0
+## 3 5002 0.07665982 0.17522245  (0.076659822,0.1752224504] -2.317063          0
+## 4 5002 0.17522245 0.18343600 (0.1752224504,0.1834360027] -4.801970          0
+## 5 5002 0.18343600 0.19438741 (0.1834360027,0.1943874059] -4.514288          0
+## 6 5002 0.19438741 0.26009582 (0.1943874059,0.2600958248] -2.722528          0
 ##    sex age cause
 ## 1 Male  60     2
 ## 2 Male  60     2
@@ -285,21 +285,21 @@ 

Cause specific haza ## ped_status ~ s(tend) + sex + age ## ## Parametric coefficients: -## Estimate Std. Error z value Pr(>|z|) -## (Intercept) -3.17557 0.58527 -5.426 5.77e-08 *** -## sexMale -0.13712 0.13226 -1.037 0.2999 -## age 0.02013 0.00847 2.376 0.0175 * +## Estimate Std. Error z value Pr(>|z|) +## (Intercept) -3.178995 0.585412 -5.430 5.62e-08 *** +## sexMale -0.136964 0.132265 -1.036 0.3004 +## age 0.020164 0.008471 2.380 0.0173 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Approximate significance of smooth terms: ## edf Ref.df Chi.sq p-value -## s(tend) 1.003 1.007 7.327 0.00687 ** +## s(tend) 1.004 1.007 7.867 0.0051 ** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## -## R-sq.(adj) = -0.00892 Deviance explained = 0.545% -## UBRE = -0.91143 Scale est. = 1 n = 29807 +## R-sq.(adj) = -0.00732 Deviance explained = 0.594% +## UBRE = -0.91774 Scale est. = 1 n = 30438 ## ## $`cause = 2` ## @@ -311,20 +311,20 @@

Cause specific haza ## ## Parametric coefficients: ## Estimate Std. Error z value Pr(>|z|) -## (Intercept) -6.95012 0.88239 -7.877 3.37e-15 *** -## sexMale 0.16840 0.18029 0.934 0.35 -## age 0.06415 0.01244 5.158 2.49e-07 *** +## (Intercept) -6.95344 0.88250 -7.879 3.29e-15 *** +## sexMale 0.16884 0.18029 0.936 0.349 +## age 0.06417 0.01244 5.160 2.47e-07 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Approximate significance of smooth terms: -## edf Ref.df Chi.sq p-value -## s(tend) 1.02 1.039 10.13 0.00155 ** +## edf Ref.df Chi.sq p-value +## s(tend) 1.017 1.033 10.51 0.00126 ** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## -## R-sq.(adj) = -0.00383 Deviance explained = 2.3% -## UBRE = -0.94888 Scale est. = 1 n = 29807

+## R-sq.(adj) = -0.00461 Deviance explained = 2.32% +## UBRE = -0.95002 Scale est. = 1 n = 30438

Cause specific hazards with shared effects @@ -339,13 +339,13 @@

Cause specific hazards with as_ped(Surv(time, status)~., id = "id", cut = cut) %>% mutate(cause = as.factor(cause)) head(ped_stacked)

-
##     id    tstart      tend                    interval    offset ped_status
-## 1 5002 0.0000000 0.1916496            (0,0.1916495551] -1.652087          0
-## 2 5002 0.1916496 0.2737851 (0.1916495551,0.2737850787] -2.499385          0
-## 3 5002 0.2737851 0.3121150 (0.2737850787,0.3121149897] -3.261525          0
-## 4 5002 0.3121150 0.3203285 (0.3121149897,0.3203285421] -4.801970          0
-## 5 5002 0.3203285 0.3285421 (0.3203285421,0.3285420945] -4.801970          0
-## 6 5002 0.3285421 0.3969884 (0.3285420945,0.3969883641] -2.681706          0
+
##     id     tstart       tend                    interval    offset ped_status
+## 1 5002 0.00000000 0.06844627            (0,0.0684462697] -2.681706          0
+## 2 5002 0.06844627 0.07665982  (0.0684462697,0.076659822] -4.801970          0
+## 3 5002 0.07665982 0.17522245  (0.076659822,0.1752224504] -2.317063          0
+## 4 5002 0.17522245 0.18343600 (0.1752224504,0.1834360027] -4.801970          0
+## 5 5002 0.18343600 0.19438741 (0.1834360027,0.1943874059] -4.514288          0
+## 6 5002 0.19438741 0.26009582 (0.1943874059,0.2600958248] -2.722528          0
 ##    sex age cause
 ## 1 Male  60     1
 ## 2 Male  60     1
@@ -368,25 +368,25 @@ 

Cause specific hazards with ## ped_status ~ s(tend, by = cause) + sex + sex:cause + age + age:cause ## ## Parametric coefficients: -## Estimate Std. Error z value Pr(>|z|) -## (Intercept) -3.17558 0.58527 -5.426 5.77e-08 *** -## sexMale -0.13712 0.13226 -1.037 0.299858 -## age 0.02013 0.00847 2.376 0.017487 * -## sexFemale:cause2 -3.77440 1.05884 -3.565 0.000364 *** -## sexMale:cause2 -3.46889 1.01238 -3.426 0.000611 *** -## age:cause2 0.04402 0.01505 2.926 0.003439 ** +## Estimate Std. Error z value Pr(>|z|) +## (Intercept) -3.179001 0.585413 -5.430 5.62e-08 *** +## sexMale -0.136962 0.132265 -1.036 0.300428 +## age 0.020164 0.008471 2.380 0.017299 * +## sexFemale:cause2 -3.774307 1.059010 -3.564 0.000365 *** +## sexMale:cause2 -3.468520 1.012503 -3.426 0.000613 *** +## age:cause2 0.044004 0.015047 2.924 0.003451 ** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Approximate significance of smooth terms: ## edf Ref.df Chi.sq p-value -## s(tend):cause1 1.001 1.002 7.333 0.00680 ** -## s(tend):cause2 1.026 1.051 10.095 0.00161 ** +## s(tend):cause1 1.001 1.003 7.873 0.00504 ** +## s(tend):cause2 1.026 1.051 10.462 0.00134 ** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## -## R-sq.(adj) = -0.00658 Deviance explained = 2.01% -## UBRE = -0.93015 Scale est. = 1 n = 59614

+## R-sq.(adj) = -0.00582 Deviance explained = 2.1% +## UBRE = -0.93388 Scale est. = 1 n = 60876

The results indicate that cause specific terms (interactions) are necessary in this case and the two models largely agree. For example, the age effect for the two causes are very similar for both models:

@@ -416,14 +416,14 @@

Cumulative Incidence Function (CIF) slice(1:3)

## # A tibble: 6 × 5
 ## # Groups:   cause [2]
-##    tend cause     cif cif_lower cif_upper
-##   <dbl> <fct>   <dbl>     <dbl>     <dbl>
-## 1 0.192 1     0.0202    0.0161     0.0254
-## 2 0.274 1     0.0289    0.0230     0.0361
-## 3 0.312 1     0.0329    0.0263     0.0412
-## 4 0.192 2     0.00973   0.00681    0.0139
-## 5 0.274 2     0.0139    0.00977    0.0199
-## 6 0.312 2     0.0159    0.0112     0.0226
+## tend cause cif cif_lower cif_upper +## <dbl> <fct> <dbl> <dbl> <dbl> +## 1 0.0684 1 0.00723 0.00565 0.00924 +## 2 0.0767 1 0.00810 0.00633 0.0103 +## 3 0.175 1 0.0185 0.0145 0.0235 +## 4 0.0684 2 0.00339 0.00235 0.00485 +## 5 0.0767 2 0.00380 0.00263 0.00543 +## 6 0.175 2 0.00871 0.00607 0.0123

 ggplot(ndf, aes(x = tend, y = cif)) +
   geom_line(aes(col = cause)) +
diff --git a/articles/competing-risks_files/figure-html/unnamed-chunk-10-1.png b/articles/competing-risks_files/figure-html/unnamed-chunk-10-1.png
index 79f6c44f..b4671912 100644
Binary files a/articles/competing-risks_files/figure-html/unnamed-chunk-10-1.png and b/articles/competing-risks_files/figure-html/unnamed-chunk-10-1.png differ
diff --git a/articles/competing-risks_files/figure-html/unnamed-chunk-11-1.png b/articles/competing-risks_files/figure-html/unnamed-chunk-11-1.png
index 4fd29e7b..3fbaf2ed 100644
Binary files a/articles/competing-risks_files/figure-html/unnamed-chunk-11-1.png and b/articles/competing-risks_files/figure-html/unnamed-chunk-11-1.png differ
diff --git a/articles/competing-risks_files/figure-html/unnamed-chunk-8-1.png b/articles/competing-risks_files/figure-html/unnamed-chunk-8-1.png
index d9f2d58b..152ba677 100644
Binary files a/articles/competing-risks_files/figure-html/unnamed-chunk-8-1.png and b/articles/competing-risks_files/figure-html/unnamed-chunk-8-1.png differ
diff --git a/articles/model-evaluation.html b/articles/model-evaluation.html
index 44cbac57..6a0c5f1e 100644
--- a/articles/model-evaluation.html
+++ b/articles/model-evaluation.html
@@ -226,10 +226,10 @@ 

Prediction e ## Integrated Brier score (crps): ## ## IBS[0.01;time=203) IBS[0.01;time=533) IBS[0.01;time=1120) -## Reference 0.065 0.119 0.169 -## pam1 0.063 0.115 0.166 -## pam2 0.058 0.108 0.159 -## pam3 0.055 0.108 0.162

+## Reference 0.062 0.116 0.168 +## pam1 0.061 0.115 0.166 +## pam2 0.057 0.110 0.163 +## pam3 0.055 0.108 0.163

C-Index @@ -256,8 +256,8 @@

C-Index## ## Pattern: ## Freq -## event 192 -## right.censored 184 +## event 182 +## right.censored 194 ## ## Censoring model for IPCW: marginal model (Kaplan-Meier for censoring distribution) ## @@ -267,9 +267,9 @@

C-Index## ## $AppCindex ## time=203 time=533 time=1120 -## pam1 67.3 60.6 59.4 -## pam2 76.7 66.5 62.9 -## pam3 76.4 64.6 56.0 +## pam1 60.9 59.9 59.3 +## pam2 70.4 64.3 60.6 +## pam3 74.3 61.8 55.3
## Warning in summary.Cindex(x, print = TRUE, ...): The C-index is not proper for t-year predictions. Blanche et al. (2018), Biostatistics, 20(2): 347--357.
 ## 
 ## Consider using time-dependent AUC instead: riskRegression::Score
diff --git a/articles/model-evaluation_files/figure-html/unnamed-chunk-3-1.png b/articles/model-evaluation_files/figure-html/unnamed-chunk-3-1.png index e83eed47..44e73462 100644 Binary files a/articles/model-evaluation_files/figure-html/unnamed-chunk-3-1.png and b/articles/model-evaluation_files/figure-html/unnamed-chunk-3-1.png differ diff --git a/articles/splines.html b/articles/splines.html index eed8cf25..c1e4d8da 100644 --- a/articles/splines.html +++ b/articles/splines.html @@ -166,10 +166,8 @@

Veteran’s dataWe first illustrate the estimation on non-linear smooth effects using the veteran data from the survival package:

-data("veteran", package="survival") # load veteran package
-
## Warning in data("veteran", package = "survival"): data set 'veteran' not found
-
-veteran <- veteran %>%
+veteran <- survival::veteran # load veteran package
+veteran <- veteran %>%
   mutate(
     trt   = 1+(trt == 2),
     prior = 1*(prior==10)) %>%
@@ -192,7 +190,7 @@ 

Veteran’s data\[ \lambda(t|x) = \exp(f_0(t_j) + f(x_{karno})) \] respectively:

-
+
 ## Cox-PH
 cph <- coxph(Surv(time, status) ~ pspline(karno, df=0), data=veteran)
 ## PAM
@@ -202,7 +200,7 @@ 

Veteran’s dataThe Figure below visualizes the smooth estimates from both models:

Expand here to see R-Code -
+
 # crate data set with varying karno values (from min to max in 100 steps)
 # and add term contribution of karno from PAM and Cox-PH models
 karno_df <- ped %>% make_newdata(karno = seq_range(karno, n = 100)) %>%
@@ -235,10 +233,10 @@ 

Cox PH workflowThe example presented in the vignette goes as follows, using the standard base R workflow and termplot for visualization:

-
+
 data("mgus", package = "survival")
## Warning in data("mgus", package = "survival"): data set 'mgus' not found
-
+
 mfit <- coxph(Surv(futime, death) ~ sex + pspline(age, df = 4), data = mgus)
 mfit
## Call:
@@ -255,7 +253,7 @@ 

Cox PH workflow## Degrees of freedom for terms= 1.0 4.1 ## Likelihood ratio test=108 on 5.04 df, p=<2e-16 ## n= 241, number of events= 225

-
+
 termplot(mfit, term = 2, se = TRUE, col.term = 1, col.se = 1)

@@ -265,7 +263,7 @@

PAM workflowmgcv::gam to fit the model:

-
+
 mgus.ped <- mgus %>% as_ped(Surv(futime, death)~sex + age, id = "id")
 head(mgus.ped)
##   id tstart tend interval   offset ped_status    sex age
@@ -275,7 +273,7 @@ 

PAM workflow## 4 1 31 32 (31,32] 0.000000 0 female 78 ## 5 1 32 39 (32,39] 1.945910 0 female 78 ## 6 1 39 60 (39,60] 3.044522 0 female 78

-
+
 pamfit <- gam(ped_status ~ s(tend) + sex + s(age, bs = "ps"), data = mgus.ped,
   method = "REML", family = poisson(), offset = offset)
 summary(pamfit)
@@ -304,7 +302,7 @@

PAM workflow## -REML = 1350.2 Scale est. = 1 n = 28326

For visualization of the smooth effects we can use the default mgcv::plot.gam function:

-
+
 layout(matrix(1:2, nrow = 1))
 termplot(mfit, term = 2, se = TRUE, col.term = 1, col.se = 1)
 plot(pamfit, select = 2, se = TRUE, ylim = c(-4, 3.5))
@@ -319,12 +317,12 @@

MGUS data: Hemoglobin example
+
 fit <- coxph(Surv(futime, death) ~ age + pspline(hgb, 4), mgus2)
 mgus2.ped <- mgus2 %>% as_ped(Surv(futime, death)~age + hgb, id = "id")
 pamfit2 <- gam(ped_status~s(tend) + age + s(hgb), data = mgus2.ped,
     family = poisson(), offset = offset)
-
+
 layout(matrix(1:2, nrow = 1))
 termplot(fit, term = 2, se = TRUE, col.term = 1, col.se = 1)
 plot(pamfit2, select = 2, se = TRUE, ylim = c(-0.5, 2))
@@ -342,11 +340,11 @@

Monotonicity constraintsbs = "mpd". Note that the fit using constraints takes considerably longer compared to the unconstrained fit.

-
+
 library(scam)
 mpam <- scam(ped_status ~ s(tend) + age + s(hgb, bs = "mpd"), data = mgus2.ped,
   family = poisson(), offset = offset)
-
+
 layout(matrix(1:2, nrow = 1))
 plot(pamfit2, select = 2, se = TRUE, ylim = c(-0.75, 2), main="Unconstrained fit")
 # 1.72 = intercept difference between pamfit2 and mpam
diff --git a/articles/strata.html b/articles/strata.html
index 805dbfde..8709939e 100644
--- a/articles/strata.html
+++ b/articles/strata.html
@@ -173,10 +173,8 @@ 

Stratified baselines
-data("veteran") # load veteran data

-
## Warning in data("veteran"): data set 'veteran' not found
-
-veteran <- veteran %>%
+veteran <- survival::veteran # load veteran data
+veteran <- veteran %>%
   mutate(
     trt   = 1*(trt == 2),
     prior = 1*(prior==10)) %>%
@@ -191,11 +189,11 @@ 

Stratified Cox model
+
 cph  <- coxph(Surv(time, status) ~ strata(celltype), data=veteran)
 base <- basehaz(cph)

The different baselines are visualized below:

-
+
 baseline_gg <- ggplot(base, aes(x=time)) +
   geom_step(aes(y=hazard, group=strata)) +
     ylab(expression(hat(Lambda)(t))) + xlab("t")
@@ -212,13 +210,13 @@ 

Stratified PAMNote that in contrast to the model specification in coxph we need to include the celltype factor variable in the model specification:

-
+
 ped <- veteran %>% as_ped(Surv(time, status)~., id="id")
 pam <- gam(ped_status ~ celltype + s(tend, by=celltype), data=ped,
     family=poisson(), offset=offset)

Comparing the estimated baselines we get very similar results between the stratified Cox and stratified PAM models.

-
+
 pinf <- ped %>%
   group_by(celltype) %>%
   ped_info() %>%
@@ -250,7 +248,7 @@ 

Stratified Proportional Hazards M \] where \(\lambda_{0k}(t)\) is the baseline for observation with cell type \(k\), \(k\in \{squamos,\ smallcell,\ adeno,\ large\}\):

-
+
 cph2 <- update(cph, .~. + trt + age + karno)
 pam2 <- update(pam, .~. + trt + age + karno)
 cbind(pam=coef(pam2)[5:7], cox=coef(cph2))
diff --git a/articles/tveffects.html b/articles/tveffects.html index 40a13c99..aeca0e7e 100644 --- a/articles/tveffects.html +++ b/articles/tveffects.html @@ -210,10 +210,8 @@

Veteran Data example
 # for some reason the prior variable is coded 0/10 instead of 0/1
-data("veteran", package = "survival")

-
## Warning in data("veteran", package = "survival"): data set 'veteran' not found
-
-veteran <- veteran %>%
+veteran <- survival::veteran
+veteran <- veteran %>%
   mutate(
     trt   = 1L * (trt == 2),
     prior = 1L * (prior == 10)) %>%
@@ -237,7 +235,7 @@ 

Extended \] This is an instance of the “known time-variation function” case above with \(g(t) = \log(t+20).\)

-
+
 vfit <- coxph(
   formula = Surv(time, status) ~ trt + prior + karno + tt(karno),
   data    = veteran,
@@ -247,7 +245,7 @@ 

Extended ## 0.07914694 0.12051224 -0.15466404 0.02930082

Thus the time-varying component of the effect becomes \(\beta_{\text{karno}}+\beta_{\text{karno},t}\cdot\log(t+20) = -0.155 + 0.029\cdot\log(t+20)\):

-
+
 t <- seq(0, 400, by = 10)
 plot(x = t, y = coef(vfit)["karno"] + coef(vfit)["tt(karno)"] * log(t + 20),
   type = "l", ylab = "Beta(t) for karno", las = 1, ylim = c(-.1, .05),
@@ -259,7 +257,7 @@ 

PAM (with known shape o

To fit a PAM with equivalent model specification (except for the baseline hazard) we can use

-
+
 # data transformation
 ped <- veteran %>% as_ped(Surv(time, status)~., id = "id") %>%
     mutate(logt20 = log(tstart + (tstart - tend) / 2 + 20))
@@ -271,7 +269,7 @@ 

PAM (with known shape o ## 4 (3,4] 0 0 60 69 0 3.113515 ## 5 (4,7] 0 0 60 69 0 3.113515 ## 6 (7,8] 0 0 60 69 0 3.277145

-
+
 # fit model
 pam <- gam(ped_status ~ s(tend) + trt + prior + karno + karno:logt20,
     data = ped, offset = offset, family = poisson())
@@ -283,7 +281,7 @@ 

PAM (with known shape o ## prior 0.11689786 0.12051224 ## karno -0.15921724 -0.15466404 ## karno:logt20 0.03049036 0.02930082

-
+
 # compare fits
 plot(x = t, y = coef(vfit)["karno"] + coef(vfit)["tt(karno)"] * log(t + 20),
   type = "l", ylab = "Beta(t) for karno", ylim = c(-.1, .05), las = 1,
@@ -304,14 +302,14 @@ 

PAM (penalized esti

In case we don’t want to pre-specify which shape the time-dependency should have, we can specify the effect of karno as \(f(x_{\text{karno}},t) = f(t)\cdot x_{\text{karno}}\), where \(f(t)\) is estimated from the data:

-
+
 # no need to specify main effect for karno here
 pam2 <- gam(ped_status ~ s(tend) + trt + prior + s(tend, by = karno),
     data = ped, offset = offset, family = poisson())
The figure below visualizes the estimate of the time-varying effect of karno for all three models
Expand here to see R-Code -
+
 term.df <- ped %>% ped_info() %>% add_term(pam2, term = "karno") %>%
     mutate_at(c("fit", "ci_lower", "ci_upper"), funs(. / .data$karno)) %>%
     mutate(
@@ -327,7 +325,7 @@ 

PAM (penalized esti ## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE)) ## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was ## generated.

-
+
 gg_tv_karno <- ggplot(term.df, aes(x = tend, y = fit)) +
     geom_step(aes(col = "PAM with penalized spline")) +
     geom_stepribbon(aes(ymin = ci_lower, ymax = ci_upper), alpha = 0.2) +
@@ -361,7 +359,7 @@ 

PAM includes the shape of the baseline hazard as well. The level of the baseline log hazard is given by the intercept of the model.

-
+
 # Non-linear, non-linearly time-varying effects
 pam3 <- gam(
   formula = ped_status ~ trt + prior + s(age) + te(tend, karno),
@@ -372,7 +370,7 @@ 

PAM function \(\hat{f}(x_{\text{karno}}, t)\) is highly non-linear (\(edf \approx 8.9\)):

-
+
 summary(pam3)
## 
 ## Family: poisson 
@@ -403,7 +401,7 @@ 

PAM displays the contribution of the effect to the log-hazard for each combination of \(x_{\text{karno}}\) and \(t\).1

-
+
 plot(pam3, select = 3, scheme = 1, theta = 120, ticktype = "detailed")

Such 3D plots are sometimes difficult to interpret, thus we also provide a heat-/contourplot (left panel) with respective slices for @@ -430,7 +428,7 @@

PAM over time as the information collected at the beginning of the follow-up becomes outdated (but see uncertainty).

Expand here to see R-Code -
+
 # heat map/contour plot
 te_gg <- gg_tensor(pam3) +
   geom_vline(xintercept = c(1, 51, 200), lty = 3) +
@@ -443,7 +441,7 @@ 

PAM theme(legend.position = "bottom")

## Scale for fill is already present.
 ## Adding another scale for fill, which will replace the existing scale.
-
+
 # plot f(karno, t) for specific slices
 karno_df <- ped %>%
   make_newdata(tend = unique(tend), karno = c(40, 75, 95)) %>%
@@ -487,7 +485,7 @@ 

PAM not include the estimated average time-constant log-hazard (coefficients(pam3)["(Intercept)"]=-4.966) and its uncertainty.

-
+
 gg_tensor(pam3, ci = TRUE) +
   xlab("t") + ylab(expression(x[plain(karno)]))

diff --git a/favicon-16x16.png b/favicon-16x16.png index a30da37a..de76f2ee 100644 Binary files a/favicon-16x16.png and b/favicon-16x16.png differ diff --git a/favicon-32x32.png b/favicon-32x32.png index 75cfc9e5..c24e264d 100644 Binary files a/favicon-32x32.png and b/favicon-32x32.png differ diff --git a/pkgdown.yml b/pkgdown.yml index 1b3cae73..aec9c2e7 100644 --- a/pkgdown.yml +++ b/pkgdown.yml @@ -17,7 +17,7 @@ articles: strata: strata.html tdcovar: tdcovar.html tveffects: tveffects.html -last_built: 2024-02-24T15:31Z +last_built: 2024-02-24T19:57Z urls: reference: https://adibender.github.io/pammtools/reference article: https://adibender.github.io/pammtools/articles