> library(dplyr) # %>%, mutate(), select(), filter()
> serodata <- read.csv("serology data.csv")
This function plots titers (x) as a function of time (y), together with a loess smoothing and its 95% confidence interval. Age of original data is slightly jittered on the plot to make it easier to vizualize.
> plot_titer <- function(x, y, ...) {
+ x2 <- jitter(x)
+ plot(x2, y, type = "n", ...)
+ loess_model <- loess(y ~ x)
+ x_val <- seq(min(x), max(x), le = 512)
+ predictions <- predict(loess_model, data.frame(x = x_val), se = TRUE)
+ with(predictions, {
+ yval <- sapply(c(-1.96, 1.96), `*`, se.fit) + fit
+ polygon(c(x_val, rev(x_val)), c(yval[, 1], rev(yval[, 2])),
+ col = "lightgrey", border = NA)
+ lines(x_val, fit)
+ })
+ points(x2, y)
+ }
This function uses plot_titer on a given pathogen (var), with the original titers (left column) and the logged titers (right column) and with healthy, all and sick flocks from top to bottom:
> plot_pathogen <- function(var) {
+ serodata2 <- serodata %>%
+ mutate(y = !! sym(var)) %>%
+ select(Types, Age, y) %>%
+ na.exclude()
+ opar <- par(mfrow = c(3, 2))
+ with(filter(serodata2, Types == "Healthy"), {
+ plot_titer(Age, y, xlab = "age (week)", ylab = "titre")
+ plot_titer(Age, log(y), xlab = "age (week)", ylab = "log(titre)")
+ })
+ with(serodata2, {
+ plot_titer(Age, y, xlab = "age (week)", ylab = "titre")
+ plot_titer(Age, log(y), xlab = "age (week)", ylab = "log(titre)")
+ })
+ with(filter(serodata2, Types != "Healthy"), {
+ plot_titer(Age, y, xlab = "age (week)", ylab = "titre")
+ plot_titer(Age, log(y), xlab = "age (week)", ylab = "log(titre)")
+ })
+ par(opar)
+ }
Which gives:
> plot_pathogen("CAV.Titre")
> plot_pathogen("IBD.Titre")
> plot_pathogen("ORT.Titre")
> plot_pathogen("PM.Titre")