## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
library(Durga)
library(datasets)
## -----------------------------------------------------------------------------
# Inspect the insulin data set, paired data in long format
head(insulin[insulin$id < 4, ])
# Calculate differences in group means
d <- DurgaDiff(insulin, data.col = "sugar", group.col = "treatment", id.col = "id")
d
## -----------------------------------------------------------------------------
# Use a formula to define the model
DurgaDiff(sugar ~ treatment, insulin, id.col = "id")
## -----------------------------------------------------------------------------
# Inspect the wide insulin data set
head(insulin.wide)
# Calculate differences in group means
(d <- DurgaDiff(insulin.wide, groups = c("sugar.before", "sugar.after")))
## ----fig.height=4, fig.width=4------------------------------------------------
DurgaPlot(d)
## ----fig.height=4, fig.width=8------------------------------------------------
# Show two plots in the same figure, and increase margins so they don't look crowded
par(mfrow = c(1, 2), mar = c(5, 5, 4, 5) + 0.1)
DurgaPlot(d,
# Display box plots
box = TRUE,
# Don't display data points
points = FALSE,
# Don't display paired-value lines
paired = FALSE,
main = "Box plot")
par(mar = c(5, 5, 4, 1) + 0.1)
DurgaPlot(d,
# Display bar charts
bar = TRUE,
# Position effect size below
ef.size.position = "below",
# No paired-value lines
paired = FALSE,
main = "Bar chart")
## ----fig.height=4, fig.width=5------------------------------------------------
# No title on this plot, so reduce the top margin size
par(mar = c(5, 4, 1, 1) + 0.1)
# Immature males are yellow, mature are red. Display juveniles first, and label
# with colour rather than developmental stage
d <- DurgaDiff(mass ~ maturity, damselfly, na.rm = TRUE,
groups = c("Yellow" = "juvenile", "Red" = "adult"))
DurgaPlot(d,
# Custom y-axis label
left.ylab = "Body mass (mg)",
# Specify the colours to (approximately) match the data
group.colour = c("orange", "red"),
box = TRUE, box.fill = FALSE)
## ----fig.height=4, fig.width=5------------------------------------------------
# No title on this plot, so reduce the top margin size
par(mar = c(5, 4, 1, 1) + 0.1)
# Immature males are yellow, mature are red. Display juveniles first, and label
# with colour rather than developmental stage
d <- DurgaDiff(mass ~ maturity, damselfly, na.rm = TRUE,
groups = c("Juveniles" = "juvenile", "Adults" = "adult"))
d$group.names <- c(expression(italic("Juveniles")), expression(bold("Adults")))
DurgaPlot(d,
# Custom y-axis label
left.ylab = "Body mass (mg)",
# Specify the colours to (approximately) match the data
group.colour = c("orange", "red"),
box = TRUE, box.fill = FALSE)
## ----fig.height=5, fig.width=8.5----------------------------------------------
par(mfrow = c(1, 2), mar = c(11, 4, 3, 6) + 0.1)
# Only show 2 groups even though there are 3 in the data
d <- DurgaDiff(iris, data.col = "Sepal.Length", group.col = "Species",
groups = c("versicolor", "setosa"))
DurgaPlot(d,
# Custom colours for points - the points are partially transparent
points = c(DurgaTransparent("coral", .6), DurgaTransparent("darkturquoise", .8)),
# Different colurs for violin (border) and violin fill
violin = c("red", "blue"), violin.fill = c("#f8a8c840", "#a8c8f840"),
main = "Data subset; Two groups")
par(mar = c(5, 4, 3, 1) + 0.1)
d <- DurgaDiff(iris, data.col = "Sepal.Length", group.col = "Species")
# Define colours for the 3 species
bgCols <- c(DurgaTransparent("red", 0.5),
DurgaTransparent("blue", 0.5),
DurgaTransparent("green", 0.5))
DurgaPlot(d,
# Make all boxes the same colour
box = "grey50", box.fill = "grey90",
# Black for border/outline colour of points
points = "black",
# Additional data point parameters: small squares filled with colour for species
points.params = list(pch = 22, cex = 0.8, bg = bgCols[iris$Species]),
main = "Coloured point fill")
## ----fig.height=6, fig.width=9------------------------------------------------
par(mfrow = c(1, 2))
d <- DurgaDiff(iris, data.col = "Sepal.Length", group.col = "Species")
# Only show error bars (95% CI) and central tendency
DurgaPlot(d, error.bars.type = "CI",
# Custom colour for central tendency
central.tendency = rainbow(3),
points = FALSE, violin = FALSE, main = "Mean and 95% CI")
par(mar = c(5, 2.4, 4, 2) + 0.1)
DurgaPlot(d, bar = TRUE, error.bars.type = "SD", points = FALSE,
left.ylab = "", ef.size.label = "",
main = "Bar chart with std. deviation")
## ----fig.height=4, fig.width=8------------------------------------------------
par(mfrow = c(1, 2))
di <- DurgaDiff(iris, "Sepal.Length", "Species")
DurgaPlot(di, contrasts = "*", main = "Problematic effect sizes")
mtext("a)", cex = 1.5, col = "brown", font = 2, adj = -0.3, line = 1.4)
DurgaPlot(di, main = "Default effect sizes (subset)")
mtext("b)", cex = 1.5, col = "brown", font = 2, adj = -0.3, line = 1.4)
## ----fig.height=4, fig.width=7------------------------------------------------
d <- DurgaDiff(petunia, 1, 2)
# Use the top margin for brackets
op <- par(mar = c(2, 4, 4, 1) + 0.1)
# Custom group colours
cols = rainbow(3)
# Save return value from DurgaPlot
p <- DurgaPlot(d,
# Don't draw effect size because we will draw brackets
ef.size = FALSE,
# Don't draw frame because brackets will appear in the upper margin
frame.plot = FALSE,
# Custom group colours
group.colour = cols)
# Customise colours by drawing brackets with the colour of the taller group.
# Get the taller group from each difference
tallerIdx <- sapply(d$group.differences, function(diff) {
# If the group difference > 0, the first group is taller
ifelse(diff$t0 > 0, diff$groupIndices[1], diff$groupIndices[2])
})
# Draw the brackets
DurgaBrackets(p, labels = "level CI", br.col = cols[tallerIdx])
## ----fig.height=4, fig.width=8------------------------------------------------
par(mfrow = c(1, 2), mar = c(3, 4, 3, 1) + 0.1)
# Construct some data with multiple groups
set.seed(1) # Ensure we always get the exact same data set
multiGroupData <- data.frame(Measurement = unlist(lapply(c(10, 12, 8, 8.5, 9, 9.1),
function(m) rnorm(40, m))),
Group = rep(paste0("g", 1:6), each = 40),
Id = rep(1:6, 40))
d <- DurgaDiff(multiGroupData, 1, 2)
DurgaPlot(d,
# Visually group each pair
group.dx = c(0.18, -0.18),
# Show box plots
box = TRUE,
# Make the boxes narrow
box.params = list(boxwex = 0.4),
# Don't display effect size or data points
ef.size = FALSE, points = FALSE,
xlim = c(0.8, 6.2), # Reduce white space on plot left and right
main = "Grouped into pairs")
DurgaPlot(d,
# Shift violins left
violin.dx = -0.12,
# Shift points right
points.dx = 0.1,
# Make points small and jittered
points.params = list(cex = 0.5), points.method = "jitter",
# No central tendency indicator
central.tendency = FALSE,
ef.size = FALSE, main = "Component shifts")
## ----eval=FALSE---------------------------------------------------------------
# # Example snippet, make *.dx arguments relative to group.dx
# group.dx <- c(0.75, 0) # Shift group 1 right, leave group 2 unchanged
# DurgaPlot(..., group.dx = group.dx,
# central.tendency.dx = group.dx + 0.1, # Shift mean & error bars right
# violin.dx = group.dx + c(-0.4, -0.5), # Shift violins left
# ef.size.dx = -0.5)
## ----fig.height=6, fig.width=8------------------------------------------------
par(mar = c(5, 9, 3, 1) + 0.1)
# Define effect size labels and their positions on the y-axis
effectSizes <- c("No effect" = 0, "Large positive effect" = 0.8, "Huge positive effect" = 2)
# Give the groups friendly names
groups <- c("Self" = "self_fertilised", "Westerham" = "westerham_cross", "Inter" = "inter_cross")
# Calculate bias-corrected Cohen's d, named Hedges' g, rather than unstandardised difference in means
d <- DurgaDiff(petunia, 1, 2, effect.type = "hedges g", groups = groups)
DurgaPlot(d, violin = FALSE, main = "Hedges' g, labelled effect size",
points.method = "jitter",
# Use our ticks and labels instead of default
ef.size.ticks = effectSizes,
# Horizontal tick labels
ef.size.params = list(las = 1),
# Don't label the effect size y-axis because it won't fit with the long tick labels
ef.size.label = "")
## ----fig.height=4.5, fig.width=6----------------------------------------------
# Define a function that calculates the statistic for two vectors
log2Ratio <- function(x1, x2) log2(mean(x2) / mean(x1))
# Generate some random data
set.seed(1)
data <- data.frame(species = rep(c("Sp1", "Sp2"), each = 10),
volume = c(rnorm(10, mean = 5, sd = 2), rnorm(10, mean = 10, sd = 2)))
# Apply log2Ratio as the effect type
d <- DurgaDiff(data, "volume", "species", effect.type = log2Ratio)
# Apply custom labels to describe the effect type, since it is not Sp2 - Sp1
d$group.differences[[1]]$label.plot <- expression("log"[2] ~ italic("Sp 2") * ":" * italic("Sp 1"))
# Plot, and set the y-axis label for effect size
DurgaPlot(d, ef.size.label = expression("log"[2]~"ratio"))
# It is also possible to set the label to be used for printing to the console
d$group.differences[[1]]$label.print <- "log2 Sp 2:Sp 1"
d
## ----eval=FALSE---------------------------------------------------------------
# # Create a new function that calls effectsize::cohens_d and returns just the calculated statistic
# myCohensAv <- function(x1, x2) {
# # Note that we swap the arguments because Durga expects the calculation to be x2 - x1
# cd <- effectsize::cohens_d(x2, x1, pooled_sd = FALSE)
# # Extract and return just the Cohen's d value
# cd$Cohens_d
# }
# # Pass it in to DurgaDiff as the effect.type
# d <- DurgaDiff(petunia, 1, 2, effect.type = myCohensAv)
#
# # Specify the effect type label
# DurgaPlot(d, ef.size.label = expression("Cohen's d"))
## ----fig.height=5, fig.width=7------------------------------------------------
par(mar = c(5, 4, 3, 1) + 0.1)
# Define display order and display pretty group names
groups = c("Control 1" = "g1", "Treatment 1" = "g2",
"Control 2" = "g3", "Treatment 2" = "g4",
"Control 3" = "g5", "Treatment 3" = "g6")
d <- DurgaDiff(multiGroupData, 1, 2, 3, groups = groups,
contrasts = "g2 - g1, g4 - g3, g6 - g5")
# Reduce text size
par(cex = 0.7)
DurgaPlot(d, box = TRUE,
# Group into pairs
group.dx = c(0.15, -0.15),
# Don't display points
points = FALSE,
# Adjust plot width to reduce left/right wasted space
xlim = c(1, 6),
# Don't display pair lines
paired = FALSE,
# Hide the boxplot median line (medcol = NA) and make the boxes narrow (boxwex = 0.4)
box.params = list(medcol = NA, boxwex = 0.4),
# Display mean and error bar
central.tendency = "#8060b0a0", # Display and colour central tendency
central.tendency.symbol = "segment",
central.tendency.width = 0.07,
error.bars = "#8060b0a0", # Same colour for error bars
# Don't display the effect size violin
ef.size.violin = FALSE,
main = "Group offsets"
)
## ----fig.height=5, fig.width=7------------------------------------------------
# Add in some fake sex data to the insulin data set
data <- cbind(insulin, Sex = sample(c("Male", "Female"), nrow(insulin), replace = TRUE))
# Thin the data so that individual symbols are visible.
# Obviously these data manipulations are for plot demonstration purposes only
data <- data[data$id %in% 1:15, ]
d <- DurgaDiff(data, "sugar", "treatment", "id",
groups = c("Before insulin" = "before", "After insulin" = "after"), na.rm = TRUE)
par(mar = c(2, 4, 4, 1))
# Use different colours to show sex
cols <- RColorBrewer::brewer.pal(3, "Set1")
DurgaPlot(d,
# Don't draw containing box
frame.plot = FALSE,
# Relabel the contrasts
contrasts = c("Blood sugar change" = "after - before"),
left.ylab = "Blood sugar level",
# Customise the violins
violin.shape = c("left", "right"), violin.dx = c(-0.055, 0.055), violin.width = 0.3,
# Customise the points to display sex
points = "black",
points.params = list(bg = ifelse(data$Sex == "Female", cols[1], cols[2]),
pch = ifelse(data$Sex == "Female", 21, 24)),
# Customise the effect size
ef.size.pch = 19,
ef.size.violin = "blue",
ef.size.violin.shape = "full",
# Make mean lines match group colour
ef.size.line.col = RColorBrewer::brewer.pal(3, "Set2")[2:1],
ef.size.line.lty = 3,
central.tendency = FALSE,
main = "Effects of insulin on blood sugar")
# Add a legend
legend("topright", c("Female", "Male"), pch = c(21, 24), pt.bg = cols)
## ----fig.height=5, fig.width=6------------------------------------------------
d <- DurgaDiff(insulin, "sugar", "treatment", "id",
groups = c("Before insulin" = "before", "After insulin" = "after"), na.rm = TRUE)
# Change effect size tick label
d$group.differences[[1]]$label.plot <- "Change in sugar"
# Horizontal axis tick labels, move y-axis labels outwards to make space for the
# horizontal ticks, and ensure the margins are large enough to show them
par(las = 1, mgp = c(4, 1, 0), mar = c(3, 5.5, 4, 5.5) + 0.1)
cols <- RColorBrewer::brewer.pal(3, "Set2")
colsTr <- DurgaTransparent(cols, 0.7)
DurgaPlot(d, paired = FALSE,
# Move 2nd group to be at the same position as 1st
group.dx = c(0.2, -0.8),
# X-axis ticks are meaningless since both groups are at the same location
x.axis = FALSE,
# Adjust plot width and spacing
xlim = c(0.7, 2.1),
violin = TRUE, violin.shape = "right", violin.fill = colsTr,
# No need to specify violin position because violin.dx defaults to group.dx
# Bumpy violins
violin.adj = 1,
# Box plots are black outlines
box = "black", box.dx = c(0.06, -0.895), box.fill = colsTr,
# Don't show outliers
box.outline = FALSE,
# Thin boxes, solid thick lines
box.params = list(boxwex = 0.06, lty = 1, lwd = 2,
# normal line weight for median line
medlwd = 2,
# no end of whisker line
staplecol = NA),
# Points as solid colours, slightly smaller than default
points = cols, points.params = list(pch = 16, cex = 0.8),
points.method = "jitter", points.dx = c(-0.14, -1.14), points.spread = 0.09,
# Move effect size left to near where 2nd group would have been
ef.size.dx = -1.2,
left.ylab = "Blood sugar",
main = "Rain cloud with effect size"
)
legend("topright", c("Before insulin", "After insulin"), col = cols, fill = colsTr, bty = "n")