Trellis graphs using the lattice package

Line plots

Data

Manipulate the data from wide format to long format. We have 43 horses that have 8 observations each:

dat.uni <- dat[rep(1:43, rep(8,43)), 1:2]
epgmat <- as.data.frame(dat[,3:10])
dat.uni <- cbind(dat.uni, time = (rep(seq(0,14,2), 43)), epg = as.numeric(t(epgmat)))

Create a grouped data object:

dat.uni <- groupedData(epg ~ time | id, outer = ~ group, data = dat.uni,
   labels = list(x = "Weeks post treatment", y = "Faecal egg count"),
   units = list(y = "(epg)"))

Make a line plot of days (x axis) versus faecal egg count (y axis), conditioned by treatment group:

background <- trellis.par.get("background")
background$col <- "white"
trellis.par.set("background", background)
superpose.line <- trellis.par.get("superpose.line")
superpose.line$col <- "black"
trellis.par.set("superpose.line", superpose.line)

plot(dat.uni, outer = ~ group,
   col = "black",
   aspect = 3,
   auto.key = FALSE,
   par.strip.text = list(cex = 1.1))

\TeX \begin{figure}[h] \centering \includegraphics[width=8.01cm]{../images/line_plot_facet_lattice.png} \caption{Line plot showing faecal egg count (expressed as the number of eggs per gram of faeces) as a function of the number of weeks post treatment, conditioned by treatment group.} \label{fig:line_plot_facet_lattice} \end{figure}

Figure 2: Line plot showing faecal egg count (expressed as the number of eggs per gram of faeces) as a function of the number of weeks post treatment, conditioned by treatment group.

Frequency histograms

Note here the use of strip.custom to manipulate the appearance of the panel strips in a lattice plot.

ename <- c(rep("Enterprise A", times = 50), rep("Enterprise B", times = 50))
ename <- factor(ename, levels = c("Enterprise A", "Enterprise B"))
distance <- c(rnorm(n = 50, mean = 20, sd = 2), rnorm(n = 50, mean = 15, sd = 1))
dat <- data.frame(ename, distance)

histogram(~ distance | ename, data = dat,
   breaks = seq(from = 0, to = 40, by = 1),
   type = "count",
   xlim = c(-0.5,40.5),
   ylim = c(0, 50),
   layout = c(1,2),
   strip = strip.custom(bg = "light green", par.strip.text = list(cex = 0.70)),
   xlab = "Distance (km)",
   ylab = "Number of poultry-positive land parcels",
   scales = list(x = list(at = seq(from = 0, to = 40, by = 5), labels = seq(from = 0, to = 40, by = 5))),
   col = "lightgray")

\TeX \begin{figure}[h] \centering \includegraphics[width=8.01cm]{../images/histogram_facet_lattice.png} \caption{Frequency histogram showing the distance separating poultry-positive enterprises, conditioned by enterprise type.} \label{fig:histogram_facet_lattice} \end{figure}

Figure 3: Frequency histogram showing the distance separating poultry-positive enterprises, conditioned by enterprise type.

Bar charts

herd <- rep(c("A","B","C"), times = 4)
test <- rep(c("MCT", "CCT"), each = 6)
tx <- rep(rep(c("Vaccinated", "Control"), each = 3), times = 2)
val <- c(34.0,34.0,63.0,19.0,15.0,35.0,3.2,0,0,1.6,1.7,1.7)
dat <- data.frame(herd, test, tx, val)

barchart(val ~ herd | test, groups = tx, data = dat,
   ylim = c(0, 100),
   cex.lab = 0.80,
   xlab = "Herd",
   ylab = "Percentage of test positives",
   layout = c(1,2),
   auto.key = list(x = 0, y = 0.80, corner = c(0,0)))

\TeX \begin{figure}[h] \centering \includegraphics[width=8.01cm]{../images/bar_chart_facet_lattice.png} \caption{Bar plot showing the percentage of test positives for control and vaccinated animals, conditioned by herd and CCT-MCT group.} \label{fig:bar_chart_facet_lattice} \end{figure}

Figure 4: Bar plot showing the percentage of test positives for control and vaccinated animals, conditioned by herd and CCT-MCT group.

Dotplots

farm <- c(rep("A", times = 6), rep("B", times = 6))
spore <- c(runif(n = 6, min = 0, max = 1500), runif(n = 6, min = 0, max = 500))
dat <- data.frame(farm, spore)

dotplot(farm ~ spore, data = dat,
   layout = c(1,1),
   aspect = c(0.7),
   xlab = "Spore count",
   ylab = "Farm",
   auto.key = list(space = "right"))

Annotated barcharts

mnth <- rep(c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"), times = 4)
dist <- c(rep("Co Do", times = 12), rep("Vinh Thanh", times = 12), rep("Vinh Loi", times = 12), rep("Phuoc Long", times = 12))
n <- c(47,68,50,61,50,42,45,26,20,20,18,18,26,35,20,16,12,8,10,2,4,6,6,8,8,20,13,13,13,13,13,10,8,
8,8,0,29,16,15,11,9,18,41,74,69,60,29,25)
dat <- data.frame(mnth, dist, n)

dat$dist <- factor(dat$dist, levels = c("Vinh Loi", "Phuoc Long", "Co Do", "Vinh Thanh"))
dat$mnth <- factor(dat$mnth, levels = c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"))

ref <- c(70, 50, 20, 60)

barchart(n ~ mnth | dist, data = dat,
   ylim = c(0, 100),
   horizontal = FALSE,
   cex.lab = 0.80,
   xlab = "Month",
   ylab = "Percentage of households",
   col = "lightgray",
   scales = list(x = list(cex = 0.70), rot = 90),
   panel = function(..., panel.number) {
      panel.barchart(...)
      panel.text(x = 2, y = 90, labels = paste("n = ", ref[packet.number()], sep = ""))
   })

\TeX \begin{figure}[h] \centering \includegraphics[width=8.01cm]{../images/annotated_bar_chart_facet_lattice.png} \caption{Bar chart showing the number of households visited by calendar month, conditioned by commune.} \label{fig:annotated_bar_chart_facet_lattice} \end{figure}

Figure 5: Bar chart showing the number of households visited by calendar month, conditioned by commune.

Error bar plots

meetn <- seq(from = 1, to = 365, by = 30)
venue <- rep("Track A", times = length(meetn))
nstart <- round(runif(n = length(meetn), min = 100, max = 300), digits = 0)
ninj <- rbinom(n = length(meetn), size = nstart, prob = 0.03)
tracka <- data.frame(venue, meetn, nstart, ninj)

meetn <- seq(from = 10, to = 365, by = 60)
venue <- rep("Track B", times = length(meetn))
nstart <- round(runif(n = length(meetn), min = 100, max = 100), digits = 0)
ninj <- rbinom(n = length(meetn), size = nstart, prob = 0.01)
trackb <- data.frame(venue, meetn, nstart, ninj)

meetn <- seq(from = 5, to = 365, by = 15)
venue <- rep("Track C", times = length(meetn))
nstart <- round(runif(n = length(meetn), min = 100, max = 400), digits = 0)
ninj <- rbinom(n = length(meetn), size = nstart, prob = 0.005)
trackc <- data.frame(venue, meetn, nstart, ninj)

dat <- data.frame(rbind(tracka, trackb, trackc))
dat$meetd <- as.Date("2018-01-01", format = "%Y-%m-%d") + dat$meetn
dat$venue <- factor(dat$venue, levels = c("Track C", "Track B", "Track A"))

irisk <- epi.conf(as.matrix(cbind(dat$ninj, dat$nstart)), ctype = "prevalence", method = "exact", N = 1000, design = 1, conf.level = 0.95) * 100
dat <- cbind(dat, irisk)

x.points <- seq(from = 1, to = 365, by = 30)
x.lab <- as.Date("2018-01-01", format = "%Y-%m-%d") + x.points
x.lab <- as.character(format(x.lab, format = "%b-%y"))

segplot(meetn ~ lower + upper | venue, data = dat,
   centers = est,
   horizontal = FALSE,
   layout = c(1,3),
   xlim = c(0, 400),
   ylim = c(0, 20),
   scales = list(x = list(at = x.points, labels = x.lab)),
   xlab = "Meeting date",
   ylab = "Percentage of starts resulting in injury")

\TeX \begin{figure}[h] \centering \includegraphics[width=8.01cm]{../images/error_bar_plot_facet_lattice.png} \caption{Error bar plot showing the percentage of race starts resulting in injury as a function of calendar date, conditioned by track.} \label{fig:error_bar_plot_facet_lattice} \end{figure}

Figure 6: Error bar plot showing the percentage of race starts resulting in injury as a function of calendar date, conditioned by track.

Ranked box and whisker plots

var <- c("movn", "mrs", "movc")
est <- c(0.4, 0.7, 0.9)
low <- c(0.3, 0.6, 0.8)
hig <- c(0.5, 0.8, 1.0)
dat <- data.frame(var, est, low, hig)

Sort the data in order of est and define factors. The order of factors defines the order in which the individual box and whisker plots appear on the graph:

dat <- dat[sort.list(dat$est, decreasing = TRUE),]
dat$var <- factor(dat$var, levels = c("movc","mrs","movn"),
   labels = c("Number of contacts following movement", "Date of movement restriction start", "Number of movements per day"))

var <- rep(dat$var, times = 3)
vals <- c(dat[,2], dat[,3], dat[,4])
dat <- data.frame(var, vals)

bwplot(var ~ vals, data = dat,
   aspect = 1/2,
   box.ratio = 0.6,
   xlab = "Partial rank correlation coefficient",
   ylab = "Parameter",
   cex = 0.5, col = 1,
   scales = list(x = list(tick.number = 5), y = list(tick.number = 5)),
   panel = function(...){
      panel.bwplot(...)
      panel.abline(v = 0, lty = 2, col = 1)
   })

Time series plots

edate <- seq(from = as.Date("1/7/08", format = "%d/%m/%y"), to = as.Date("30/6/09", format = "%d/%m/%y"), by = "1 day")
bw2yo <- rnorm(n = 365, mean = 400, sd = 35)
bw3yo <- rnorm(n = 365, mean = 450, sd = 15)
bw46yo <- rnorm(n = 365, mean = 525, sd = 30)

dat <- data.frame(edate = rep(edate, times = 3),
   class = c(rep("2yo", times = 365), rep("3yo", times = 365), rep("4+yo", times = 365)),
   bw = c(bw2yo, bw3yo, bw46yo))

Custom strip function:

my.strip <- function(which.given, which.panel, ...){
   strip.labels <- c("2 year olds", "3 year olds", "4 - 6 year olds")
   panel.rect(0,0,1,1, col = "#ffe5cc", border = 1)
   panel.text(x = 0.5, y = 0.5, adj = c(0.5,0.55), cex = 0.95, lab = strip.labels[which.panel[which.given]])
   }

Define x axis date range:

xlim <- range(edate)
x.points <- seq(from = as.Date("1/7/08", format = "%d/%m/%y"), to = as.Date("30/6/09", format = "%d/%m/%y"), by = "2 months")
x.lab <- as.character(format(x.points, format = "%b-%y"))

xyplot(bw ~ edate | class, data = dat,
   xlim = xlim,
   ylim = c(250, 750),
   strip = my.strip,
   outer = TRUE,
   layout = c(1,3,1),
   xlab = "Date",
   ylab = "Bodyweight (kg)",
   scales = list(x = list(at = x.points, labels = x.lab)),
   type = c("p", "smooth"), grid = TRUE, col.line = "darkorange", lwd = 4)

\TeX \begin{figure}[h] \centering \includegraphics[width=8.01cm]{../images/time_series_plot_facet_lattice.png} \caption{Time series plot showing bodyweight (kg) as a function of calendar date, conditioned by age class.} \label{fig:time_series_plot_facet_lattice} \end{figure}

Figure 7: Time series plot showing bodyweight (kg) as a function of calendar date, conditioned by age class.

Comparing risk ratios from several models

library(epiR); library(lattice); library(latticeExtra); library(nlme); library(tidyr); library(dplyr)

exp.a <- gl(n = 3, k = 3); exp.b <- gl(3,1,9); counts <- c(18,17,15,20,10,20,25,13,12)
dat <- data.frame(exp.a, exp.b, counts)
dat.glm01 <- glm(counts ~ exp.a + exp.b, family = poisson, data = dat)

exp.a <- gl(n = 3, k = 3); exp.b <- gl(3,1,9); counts <- c(9,8,7,10,5,10,12,7,6)
dat <- data.frame(exp.a, exp.b, counts)
dat.glm02 <- glm(counts ~ exp.a + exp.b, family = poisson, data = dat)

Create a data frame of risk ratios:

model <- c(rep("Model 1", times = nrow(summary(dat.glm01)$coefficients)), rep("Model 2", times = nrow(summary(dat.glm01)$coefficients)))
vname <- c(row.names(summary(dat.glm01)$coefficients), row.names(summary(dat.glm02)$coefficients))
m1.est <- c(exp(summary(dat.glm01)$coefficients[,1]), exp(summary(dat.glm02)$coefficients[,1]))
m1.low <- c(exp(confint(dat.glm01))[,1], exp(confint(dat.glm02))[,1])
m1.upp <- c(exp(confint(dat.glm01))[,2], exp(confint(dat.glm02))[,2])
dat <- data.frame(model = model, vname = vname, est = m1.est, low = m1.low, upp = m1.upp)

Add risk ratios for each of the reference categories for each explanatory variable to the data frame:

exp.a <- data.frame(model = c("Model 1", "Model 2"), vname = c("exp.a1", "exp.a1"), est = c(1,1), low = c(1,1), upp = c(1,1))
exp.b <- data.frame(model = c("Model 1", "Model 2"), vname = c("exp.b1", "exp.b1"), est = c(1,1), low = c(1,1), upp = c(1,1))
dat <- rbind(dat, exp.a, exp.b)
dat <- subset(dat, subset = vname != "(Intercept)")
dat$vname <- factor(dat$vname, levels = c("exp.a1", "exp.a2", "exp.a3", "exp.b1", "exp.b2", "exp.b3"), labels = c("Size: small", "Size: medium", "Size: large", "Type: mixed", "Type: dairy", "Type: beef"))

Make the plot:

segplot(vname ~ low + upp | model, centers = est, data = dat,
   draw.bands = FALSE,
   segments.fun = panel.arrows,
   ends = "both",
   angle = 90,
   length = 0.03,
   layout = c(1,2),
   scales = list(x = list(at = c(0.5,1,2,4), lim = c(0.25,2))),
   xlab = expression("Incidence rate ratio"),
   ylab = "Explanatory variable",
   panel = function(x, y, ...) {
      panel.abline(v = 1, col = "grey", lty = 2)
      panel.segplot(x, y, ...)})

\TeX \begin{figure}[h] \centering \includegraphics[width=8.01cm]{../images/compare_models_facet_lattice.png} \caption{Error bar plot showing the point estimate of adjusted incidence rate ratios (and their 95\% confidence intervals) for herd size and herd type for two Poisson regression models.} \label{fig:compare_models_facet_lattice} \end{figure}

Figure 1: Error bar plot showing the point estimate of adjusted incidence rate ratios (and their 95% confidence intervals) for herd size and herd type for two Poisson regression models.