I am trying to create a faceted survival plot with confidence intervals and p-values. The plot should look essentially equivalent to the below but with p-values.
However, I keep running into the error variable lengths differ.
I have found several posts with similar problems (Post_1, Post_2, Post_3), but none of these have yet helped me to solve my issue. I have also tried sub-setting my data in various ways so that all included variables were evenly balanced, but this didn't help either.
Any advice appreciated.
library(survival)
library(survminer)
library(ggsurvfit)
# Read in data
dat <- structure(list(Rabbit = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L,
10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L,
23L, 24L, 25L, 27L, 26L, 34L, 35L, 28L, 29L, 32L, 36L, 30L, 31L,
33L, 37L, 38L, 39L, 40L, 41L, 42L), Treatment = structure(c(3L,
3L, 4L, 4L, 2L, 2L, 1L, 1L, 5L, 3L, 3L, 4L, 4L, 2L, 2L, 1L, 1L,
5L, 3L, 3L, 4L, 4L, 2L, 2L, 1L, 1L, 5L, 3L, 3L, 4L, 4L, 2L, 2L,
1L, 1L, 5L, 4L, 4L, 5L, 4L, 4L, 5L), levels = c("Meat bait",
"Soil spray", "Carrot bait", "Oat bait", "Control"), class = "factor"),
Survival.Time.Rabbit = c(68.5, 75, 51, 51, 99, 120, 240,
219, 336, 53, 29, 77, 77, 96, 149, 91.5, 77, 336, 336, 336,
77.67, 92.67, 336, 336, 336, 336, 336, 336, 336, 53.5, 73,
336, 336, 336, 336, 336, 336, 336, 336, 336, 336, 336), Status.Rabbit = c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), Timepoint = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L), levels = c("1 Day post deployment",
"5 Days post deployment", "10 Days post deployment", "20 Days post deployment",
"40 Days post deployment", "60 Days post deployment"), class = "factor"),
Survival.Time.Bait1 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 0L, NA, NA, 10L, 10L, NA,
NA, NA, 10L, 0L, NA, NA, 20L, 20L, NA, NA, NA, NA, 0L, NA,
NA, 0L, NA, NA, 0L), Survival.Time.Bait2 = c(NA, NA, NA,
NA, NA, NA, NA, NA, 1L, NA, NA, NA, NA, NA, NA, NA, NA, 5L,
10L, 10L, NA, NA, 10L, 10L, 10L, NA, 10L, 20L, 20L, NA, NA,
20L, 20L, 20L, 20L, 20L, 40L, 40L, 40L, 60L, 60L, 60L)), row.names = c(NA,
-42L), class = "data.frame")
# Define factors
dat$Timepoint <- factor(dat$Timepoint, levels = c("1 Day post deployment",
"5 Days post deployment",
"10 Days post deployment",
"20 Days post deployment",
"40 Days post deployment",
"60 Days post deployment"))
dat$Treatment <- factor(dat$Treatment, levels = c("Meat bait",
"Soil spray",
"Carrot bait",
"Oat bait",
"Control"))
# Create survival objects
surv <- Surv(time = dat$Survival.Time.Rabbit, event = dat$Status.Rabbit)
survfit.treatment.timepoint.rabbit <- survfit2(surv ~ Treatment, conf.int = T, data = dat)
# Plot
ggsurvplot_facet(survfit.treatment.timepoint.rabbit,
facet.by = "Timepoint", size = 1.5,
data = dat, conf.int = T, font.x = c(12, face = "bold"),
font.y = c(12, face = "bold"), font.tickslab = c(10), pval = T) +
labs(x = "Time to death (hours)",
y = "Survival probability")
# Note code for plot above runs fine in absence of pval = T
