Packages:
> library(dplyr)
> library(purrr)
> library(magrittr)
> library(tidyr)Reading the data:
> (data <- readr::read_csv("Figure1.csv") %>%
+ select(-CODE) %>%
+ group_by(FlockID, WEEKNO) %>%
+ summarise_all(sum) %>%
+ ungroup())
Parsed with column specification:
cols(
FlockID = col_character(),
WEEKNO = col_double(),
Fowl.Pox.vaccinated = col_double(),
HPAI.vaccinated = col_double(),
NDV.vaccinated = col_double(),
IBV.vaccinated = col_double(),
IBD.vaccinated = col_double(),
PM.vaccinated = col_double(),
CODE = col_character()
)
# A tibble: 4,476 x 8
FlockID WEEKNO Fowl.Pox.vaccinated HPAI.vaccinated NDV.vaccinated IBV.vaccinated IBD.vaccinated PM.vaccinated
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 75-0011 1 0 0 1 0 1 0
2 75-0011 2 1 0 0 0 0 0
3 75-0011 3 0 1 1 0 0 0
4 75-0011 4 0 0 0 0 1 0
5 75-0011 5 0 0 0 0 0 1
6 75-0011 6 0 1 0 0 0 0
7 75-0011 7 0 0 0 0 0 0
8 75-0011 8 0 0 0 0 0 0
9 75-0011 9 0 0 0 0 0 0
10 75-0011 10 0 0 0 0 0 0
# … with 4,466 more rowsChecking that there are no missing weeks:
> sel <- data %>%
+ arrange(WEEKNO) %>%
+ group_by(FlockID) %>%
+ summarise(x = length(unique(diff(WEEKNO)))) %>%
+ filter(x > 1) %>%
+ purrr::pluck("FlockID")
>
> data %>%
+ filter(FlockID %in% sel)
# A tibble: 17 x 8
FlockID WEEKNO Fowl.Pox.vaccinated HPAI.vaccinated NDV.vaccinated IBV.vaccinated IBD.vaccinated PM.vaccinated
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 75-0046 1 0 0 1 0 0 0
2 75-0046 2 1 0 0 0 1 0
3 75-0046 3 0 1 0 0 0 0
4 75-0046 4 0 0 1 0 1 0
5 75-0046 5 0 0 0 0 0 0
6 75-0046 6 0 0 0 0 0 1
7 75-0046 7 0 1 0 0 0 0
8 75-0046 8 0 0 0 0 0 0
9 75-0046 9 0 0 0 0 0 0
10 75-0046 10 0 0 0 0 0 0
11 75-0046 11 0 0 0 0 0 0
12 75-0046 12 0 0 0 0 0 0
13 75-0046 13 0 0 0 0 0 0
14 75-0046 14 0 0 0 0 0 0
15 75-0046 15 0 0 0 0 0 0
16 75-0046 16 0 0 0 0 0 0
17 75-0046 18 0 0 0 0 0 0Prepare the data:
> tmp <- data %>%
+ group_by(FlockID) %>%
+ tally() %>%
+ arrange(desc(n))The code below draws the figure, directely to a A4 PDF file:
> # opening an A4 PDF file:
> pdf("figure.pdf", width = 8, height = 11.7)
>
> # tuning some graphical parameters:
> opar <- par(plt = c(.1, .95, .05, 1))
>
> # plotting the weeks:
> a <- table(tmp$n)
> b <- rev(cumsum(rev(a)))
> x <- as.integer(names(b)) + .5
> plot(x, b, type = "h", xlim = c(0, max(x)), ann = FALSE, axes = FALSE)
> axis(1, line = -1)
> mtext("farm", 2, -1)
> mtext("week", 1, .5)
> x_rev <- rev(x)
> x2 <- rep(x_rev, c(2, diff(rev(b))))
> y <- 0:nrow(tmp)
> segments(.5, y, x2, y)
> x3 <- .5:(min(x) - 1)
> segments(x3, 0, x3, max(b))
>
> # adding vaccination information:
> coords <- seq(.15, .85, le = 6)
> vaccines <- grep("vaccinated", names(data), value = TRUE)
> cols <- setNames(c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3", "#ff7f00", "#ffff33"), vaccines)
> data %>%
+ filter(select(., ends_with("vaccinated")) %>% rowSums() %>% as.logical) %>%
+ mutate(Fowl.Pox.vaccinated = Fowl.Pox.vaccinated * coords[1],
+ HPAI.vaccinated = HPAI.vaccinated * coords[2],
+ NDV.vaccinated = NDV.vaccinated * coords[3],
+ IBV.vaccinated = IBV.vaccinated * coords[4],
+ IBD.vaccinated = IBD.vaccinated * coords[5],
+ PM.vaccinated = PM.vaccinated * coords[6],
+ y_val = setNames(seq_len(nrow(tmp)) - .5, tmp$FlockID)[FlockID]) %>%
+ mutate_at(vars(ends_with("vaccinated")), funs(`*`(. > 0, `+`(WEEKNO, .)))) %>%
+ pivot_longer(cols = ends_with("vaccinated"), names_to = "vaccine", values_to = "x_val") %>%
+ mutate(col = cols[vaccine]) %>%
+ filter(x_val > 0) %$%
+ points(x_val - .5, y_val, pch = 15, col = col, cex = .5)
>
> # adding the legend:
> legend("right", legend = sub("\\.", " ", sub(".vaccinated", "", vaccines)),
+ fill = cols, bty = "n", title = "Vaccines:")
>
> # closing the PDf file:
> dev.off()
>
> # setting graphical parameters back to their initial values
> par(opar)