work on survival vignette

vignettes/survival.qmd

@@ -16,16 +16,19 @@ vignette: >
   %\VignettePackage{gamlss2}
 ---
 
+Here is the text for survival model introduction. Important is, how to estimate survival models with the GAMLSS model class.
+
 ```{r packages}
 #| message: false
 #| results: hide
 #| echo: 3:8
-pkg <- c("spBayesSurv", "gamlss.cens")
+pkg <- c("spBayesSurv", "gamlss.cens", "sf")
 for(p in pkg) {
   if(!(p %in% installed.packages())) install.packages(p)
 }
 library("gamlss2")
 library("gamlss.cens")
+library("sf")
 data("LeukSurv", package = "spBayesSurv")
 ```
 
@@ -48,9 +51,11 @@ s1 <- survfit(Surv(time, cens) ~ 1, data = LeukSurv)
 Predict survival curves.
 
 ```{r}
-## gamlss2 survival
-p1 <- predict(b1, newdata = LeukSurv)
-p1 <- 1 - family(b1)$p(Surv(LeukSurv$time, rep(1, nrow(LeukSurv))), p1)
+## gamlss2 survival, predict parameters
+par <- predict(b1, newdata = LeukSurv)
+
+## compute survival probabilities
+p1 <- 1 - family(b1)$p(Surv(LeukSurv$time, rep(1, nrow(LeukSurv))), par)
 
 ## surfit() survival
 f <- stepfun(s1$time, c(1, s1$surv))
@@ -82,12 +87,80 @@ b2 <- gamlss2(f, family = fam, data = LeukSurv)
 Plot estimated effects.
 
 ```{r}
+#| fig-height: 8
+#| fig-width: 6
+#| fig-align: center
 plot(b2)
 ```
 
 Residuals diagnostic plots.
 
 ```{r}
+#| fig-align: center
 plot(b2, which = "resid")
 ```
 
+Visualize the spatial effect.
+
+```{r}
+#| fig-align: center
+## read the map of new west england.
+file <- system.file("otherdata/nwengland.bnd", package = "spBayesSurv")
+d <- readLines(file)
+
+## transform the polygons to a list().
+id <- grep('\"', d, fixed = TRUE)
+polys <- list()
+for(i in 1:length(id)) {
+  j <- strsplit(d[id[i]], ",")[[1]][2]
+  if(i < length(id))
+    polys[[j]] <- d[(id[i] + 1):(id[i + 1] - 1)]
+  else
+    polys[[j]] <- d[(id[i] + 1):length(d)]
+}
+
+polys <- lapply(polys, function(x) {
+  tf <- tempfile()
+  writeLines(x, tf)
+  pol <- as.matrix(read.csv(tf, header = FALSE))
+  unlink(tf)
+  return(st_polygon(list(pol)))
+})
+
+## transform to sf object
+polys_sfc <- st_sfc(polys)
+map <- st_sf(geometry = polys_sfc)
+map$id <- names(polys)
+map$district <- 1:nrow(map)
+
+## plot the map
+par(mar = rep(0, 4))
+plot(st_geometry(map))
+
+## sample coordinates for plotting
+co <- st_sample(map, size = 1000, type = "regular")
+
+## create new data for prediction
+nd <- as.data.frame(st_coordinates(co))
+names(nd) <- c("xcoord", "ycoord")
+nd$sex <- 1
+nd$wbc <- mean(LeukSurv$wbc)
+nd$tpi <- mean(LeukSurv$tpi)
+nd$age <- mean(LeukSurv$age)
+
+## predict parameters
+par <- predict(b2, newdata = nd)
+
+## compute survival probabilities
+p2 <- 1 - family(b2)$p(Surv(rep(365, nrow(nd)), rep(1, nrow(nd))), par)
+
+## plot on map
+sp <- st_sf(geometry = co)
+sp$survprob <- p2
+
+par(mar = rep(0, 4))
+plot(st_geometry(map))
+plot(sp, pch = 15, cex = 1.3, add = TRUE)
+plot(st_geometry(map), add = TRUE)
+```
+