Fitted-Model-Based Annotations
-‘ggpmisc’ 0.5.5.9003
+‘ggpmisc’ 0.5.6.9001
Pedro J. Aphalo
-2024-05-05
+2024-05-13
@@ -1849,8 +1849,9 @@
+ stat_poly_eq(aes(label = paste("atop(", after_stat(AIC.label), ",",
+ after_stat(BIC.label), ")", sep = "")),
+ formula = formula)
stat_poly_eq() and stat_poly_line()
ggplot(my.data, aes(x, y)) + geom_point() + stat_poly_line(formula = formula) + - stat_poly_eq(aes(label = paste("atop(", after_stat(AIC.label), ",", after_stat(BIC.label), ")", sep = "")), - formula = formula)
Next, one example of how to remove the left hand side (lhs).
@@ -2218,7 +2219,7 @@stat_quant_eq(), stat_quant_line() and stat_quant_band()
ggplot(my.data, aes(x, y)) +
geom_point() +
stat_quant_band(formula = formula, color = "black", fill = "grey60") +
- stat_quant_eq(aes(label = paste(after_stat(grp.label), "*\": \"*",
+ stat_quant_eq(aes(label = paste(after_stat(qtl.label), "*\": \"*",
after_stat(eq.label), sep = "")),
formula = formula) +
theme_classic()stat_quant_eq(), stat_quant_line() and stat_quant_band()
ggplot(my.data, aes(x, y, color = group)) +
geom_point() +
stat_quant_line(formula = formula) +
- stat_quant_eq(aes(label = paste(after_stat(grp.label), "*\": \"*",
+ stat_quant_eq(aes(label = paste(after_stat(qtl.label), "*\": \"*",
after_stat(eq.label), sep = "")),
formula = formula)
stat_quant_eq(), stat_quant_line() and stat_quant_band()
shape = group, grp.label = group)) + geom_point() + stat_quant_line(formula = formula, quantiles = c(0.05, 0.95), color = "black") + - stat_quant_eq(aes(label = paste(after_stat(grp.label), "*\": \"*", - after_stat(eq.label), sep = "")), - formula = formula, quantiles = c(0.05, 0.95)) + - theme_classic() + stat_quant_eq(aes(label = paste(after_stat(grp.label), "*\", \"*", + after_stat(qtl.label), "*\": \"*", + after_stat(eq.label), sep = "")), + formula = formula, quantiles = c(0.05, 0.95)) + + theme_classic()
In some cases double quantile regression is an informative method to assess reciprocal constraints. For this a fit of x on @@ -2285,112 +2287,139 @@
stat_multcomp
adj.method.tag = 3, size = 2.75) + expand_limits(y = 0) -
## Warning: Computation failed in `stat_multcomp()`.
+## Caused by error in `compute_panel()`:
+## ! object 'tstat.char' not found
Here we use a negative value to use an abbreviation of the word “adjusted” to three characters.
-# position of contrasts' bars (manual)
-ggplot(mpg, aes(factor(cyl), hwy)) +
- geom_boxplot(width = 0.33) +
- stat_multcomp(p.adjust.method = "bonferroni",
- adj.method.tag = -3,
- size = 2.75) +
- expand_limits(y = 0)
A character string passed as argument is used as is, here to set the -tag in Spanish.
# position of contrasts' bars (manual)
ggplot(mpg, aes(factor(cyl), hwy)) +
geom_boxplot(width = 0.33) +
- stat_multcomp(adj.method.tag = "ajustada",
- size = 2.75) +
- expand_limits(y = 0)
## Warning: Computation failed in `stat_multcomp()`.
+## Caused by error in `compute_panel()`:
+## ! object 'tstat.char' not found
A character string passed as argument is used as is, here to set the +tag in Spanish.
+# position of contrasts' bars (manual)
+ggplot(mpg, aes(factor(cyl), hwy)) +
+ geom_boxplot(width = 0.33) +
+ stat_multcomp(adj.method.tag = "ajustada",
+ size = 2.75) +
+ expand_limits(y = 0)## Warning: Computation failed in `stat_multcomp()`.
+## Caused by error in `compute_panel()`:
+## ! object 'tstat.char' not found
A numeric vector passed to label.y can be used to
manually set the location of each label along the y axis.
# position of contrasts' bars (manual)
-ggplot(mpg, aes(factor(cyl), hwy)) +
- geom_boxplot(width = 0.33) +
- stat_multcomp(label.y = c(7, 4, 1),
- contrasts = "Dunnet",
- size = 2.75) +
- expand_limits(y = 0)
ggplot(mpg, aes(factor(cyl), hwy)) +
- geom_boxplot(width = 0.33) +
- stat_multcomp(label.y =
- seq(from = 15,
- by = -3,
- length.out = 6),
- size = 2.5) +
- expand_limits(y = 0)
# position of contrasts' bars (manual)
+ggplot(mpg, aes(factor(cyl), hwy)) +
+ geom_boxplot(width = 0.33) +
+ stat_multcomp(label.y = c(7, 4, 1),
+ contrasts = "Dunnet",
+ size = 2.75) +
+ expand_limits(y = 0)## Warning: Computation failed in `stat_multcomp()`.
+## Caused by error in `compute_panel()`:
+## ! object 'tstat.char' not found
ggplot(mpg, aes(factor(cyl), hwy)) +
+ geom_boxplot(width = 0.33) +
+ stat_multcomp(label.y =
+ seq(from = 15,
+ by = -3,
+ length.out = 6),
+ size = 2.5) +
+ expand_limits(y = 0)## Warning: Computation failed in `stat_multcomp()`.
+## Caused by error in `compute_panel()`:
+## ! object 'tstat.char' not found
We can pre-compute the numeric vector to achieve special positioning, in this case, next to each observation.
-means <-
- aggregate(mpg$hwy,
- by = list(mpg$cyl),
- FUN = mean,
- na.rm = TRUE)[["x"]]
-
-ggplot(mpg, aes(factor(cyl), hwy)) +
- stat_summary(fun.data = mean_se) +
- stat_multcomp(label.type = "letters",
- label.y = c(18, means), # 18 is for critical P label
- position = position_nudge(x = 0.1))
means <-
+ aggregate(mpg$hwy,
+ by = list(mpg$cyl),
+ FUN = mean,
+ na.rm = TRUE)[["x"]]
+
+ggplot(mpg, aes(factor(cyl), hwy)) +
+ stat_summary(fun.data = mean_se) +
+ stat_multcomp(label.type = "letters",
+ label.y = c(18, means), # 18 is for critical P label
+ position = position_nudge(x = 0.1))
We can override the default use of geom_text() and also
remove the P critical label.
# Using other geometries
-ggplot(mpg, aes(factor(cyl), hwy)) +
- geom_boxplot(width = 0.33) +
- stat_multcomp(label.type = "letters",
- adj.method.tag = FALSE,
- geom = "label")# Using other geometries
+ggplot(mpg, aes(factor(cyl), hwy)) +
+ geom_boxplot(width = 0.33) +
+ stat_multcomp(label.type = "letters",
+ adj.method.tag = FALSE,
+ geom = "label")
With Dunnet contrasts, the use of bars can clutter a plot, and we can alternatively show the outcome for each level that has been compared to a control, the first level.
-ggplot(mpg, aes(factor(cyl), hwy)) +
- geom_boxplot(width = 0.33) +
- stat_multcomp(aes(x = stage(start = factor(cyl),
- after_stat = x.right.tip)),
- geom = "text",
- label.y = "bottom",
- vstep = 0,
- contrasts = "Dunnet")
ggplot(mpg, aes(factor(cyl), hwy)) +
- geom_boxplot(width = 0.33) +
- stat_multcomp(aes(x = stage(start = factor(cyl),
- after_stat = x.right.tip),
- label = after_stat(stars.label)),
- geom = "text",
- label.y = "bottom",
- vstep = 0,
- contrasts = "Dunnet")
ggplot(mpg, aes(factor(cyl), hwy)) +
+ geom_boxplot(width = 0.33) +
+ stat_multcomp(aes(x = stage(start = factor(cyl),
+ after_stat = x.right.tip)),
+ geom = "text",
+ label.y = "bottom",
+ vstep = 0,
+ contrasts = "Dunnet")## Warning: Computation failed in `stat_multcomp()`.
+## Caused by error in `compute_panel()`:
+## ! object 'tstat.char' not found
ggplot(mpg, aes(factor(cyl), hwy)) +
+ geom_boxplot(width = 0.33) +
+ stat_multcomp(aes(x = stage(start = factor(cyl),
+ after_stat = x.right.tip),
+ label = after_stat(stars.label)),
+ geom = "text",
+ label.y = "bottom",
+ vstep = 0,
+ contrasts = "Dunnet")## Warning: Computation failed in `stat_multcomp()`.
+## Caused by error in `compute_panel()`:
+## ! object 'tstat.char' not found
The returned value includes numeric values in addition to character
strings mapped to the label aesthetic. The numeric values
can be used to encode the outcomes using additional or different
aesthetics than the default.
# use colour to show significance
-ggplot(mpg, aes(factor(cyl), hwy)) +
- geom_boxplot(width = 0.33) +
- stat_multcomp(aes(colour = after_stat(p.value) < 0.01),
- size = 2.75) +
- scale_colour_manual(values = c("grey60", "black")) +
- theme_bw()
# add arrow heads to segments and use fill to show significance
-ggplot(mpg, aes(factor(cyl), hwy)) +
- geom_boxplot(width = 0.33) +
- stat_multcomp(aes(fill = after_stat(p.value) < 0.01),
- size = 2.5,
- arrow = grid::arrow(angle = 45,
- length = unit(1, "mm"),
- ends = "both")) +
- scale_fill_manual(values = c("white", "lightblue"))
# use colour to show significance
+ggplot(mpg, aes(factor(cyl), hwy)) +
+ geom_boxplot(width = 0.33) +
+ stat_multcomp(aes(colour = after_stat(p.value) < 0.01),
+ size = 2.75) +
+ scale_colour_manual(values = c("grey60", "black")) +
+ theme_bw()## Warning: Computation failed in `stat_multcomp()`.
+## Caused by error in `compute_panel()`:
+## ! object 'tstat.char' not found
# add arrow heads to segments and use fill to show significance
+ggplot(mpg, aes(factor(cyl), hwy)) +
+ geom_boxplot(width = 0.33) +
+ stat_multcomp(aes(fill = after_stat(p.value) < 0.01),
+ size = 2.5,
+ arrow = grid::arrow(angle = 45,
+ length = unit(1, "mm"),
+ ends = "both")) +
+ scale_fill_manual(values = c("white", "lightblue"))## Warning: Computation failed in `stat_multcomp()`.
+## Caused by error in `compute_panel()`:
+## ! object 'tstat.char' not found
stat_fit_residuals
@@ -2410,28 +2439,28 @@stat_fit_residuals
by the residuals from the fit and thus also by default mapped to the aesthetic of the same name. The default geom is"point".
-formula <- y ~ poly(x, 3, raw = TRUE)
-ggplot(my.data, aes(x, y, colour = group)) +
- geom_hline(yintercept = 0, linetype = "dashed") +
- stat_fit_residuals(formula = formula)formula <- y ~ poly(x, 3, raw = TRUE)
+ggplot(my.data, aes(x, y, colour = group)) +
+ geom_hline(yintercept = 0, linetype = "dashed") +
+ stat_fit_residuals(formula = formula)
We can of course also map the weights returned in the model fit
object to ggplot2 aesthetics, here we use size.
formula <- y ~ poly(x, 3, raw = TRUE)
-ggplot(my.data, aes(x, y, colour = group)) +
- geom_hline(yintercept = 0, linetype = "dashed") +
- stat_fit_residuals(formula = formula,
- method = "rlm",
- mapping = aes(size = sqrt(after_stat(weights))),
- alpha = 2/3)formula <- y ~ poly(x, 3, raw = TRUE)
+ggplot(my.data, aes(x, y, colour = group)) +
+ geom_hline(yintercept = 0, linetype = "dashed") +
+ stat_fit_residuals(formula = formula,
+ method = "rlm",
+ mapping = aes(size = sqrt(after_stat(weights))),
+ alpha = 2/3)
Weighted residuals are available, but in this case they do not differ
as we have not mapped any variable to the weight aesthetic
in the input to the statistic.
stat_fit_deviations
@@ -2445,36 +2474,36 @@stat_fit_deviations
y, and the fitted values are mapped to aesthetics xend and yend. As the default geom is"segment", each
deviation is displayed as a segment.
-formula <- y ~ poly(x, 3, raw = TRUE)
-ggplot(my.data, aes(x, y)) +
- geom_smooth(method = "lm", formula = formula) +
- stat_fit_deviations(formula = formula, colour = "red") +
- geom_point()formula <- y ~ poly(x, 3, raw = TRUE)
+ggplot(my.data, aes(x, y)) +
+ geom_smooth(method = "lm", formula = formula) +
+ stat_fit_deviations(formula = formula, colour = "red") +
+ geom_point()
The geometry used by default is ggplot2::geom_segment()
and additional aesthetics can be mapped or set to constant values. Here
we add arrowheads.
formula <- y ~ poly(x, 3, raw = TRUE)
-ggplot(my.data, aes(x, y)) +
- geom_smooth(method = "lm", formula = formula) +
- geom_point() +
- stat_fit_deviations(formula = formula, colour = "red",
- arrow = arrow(length = unit(0.015, "npc"),
- ends = "both"))formula <- y ~ poly(x, 3, raw = TRUE)
+ggplot(my.data, aes(x, y)) +
+ geom_smooth(method = "lm", formula = formula) +
+ geom_point() +
+ stat_fit_deviations(formula = formula, colour = "red",
+ arrow = arrow(length = unit(0.015, "npc"),
+ ends = "both"))
When weights are available, either supplied by the user, or computed
as part of the fit, they are returned in data. Having
weights available allows encoding them using colour. We here use a
robust regression fit with MASS::rlm().
my.data.outlier <- my.data
-my.data.outlier[6, "y"] <- my.data.outlier[6, "y"] * 5
-ggplot(my.data.outlier, aes(x, y)) +
- stat_smooth(method = MASS::rlm, formula = formula) +
- stat_fit_deviations(formula = formula, method = "rlm",
- mapping = aes(colour = after_stat(weights)),
- show.legend = TRUE) +
- scale_color_gradient(low = "red", high = "blue", limits = c(0, 1)) +
- geom_point()my.data.outlier <- my.data
+my.data.outlier[6, "y"] <- my.data.outlier[6, "y"] * 5
+ggplot(my.data.outlier, aes(x, y)) +
+ stat_smooth(method = MASS::rlm, formula = formula) +
+ stat_fit_deviations(formula = formula, method = "rlm",
+ mapping = aes(colour = after_stat(weights)),
+ show.legend = TRUE) +
+ scale_color_gradient(low = "red", high = "blue", limits = c(0, 1)) +
+ geom_point()
@@ -2488,37 +2517,37 @@
+
+
stat_fit_glance
statistic in a plot with grouping, and assemble a label for the P-value using a string parsed into a expression. We also change the default position of the labels. -# formula <- y ~ poly(x, 3, raw = TRUE)
-# broom::augment does not handle poly() correctly!
-formula <- y ~ x + I(x^2) + I(x^3)
-ggplot(my.data, aes(x, y, colour = group)) +
- geom_point() +
- geom_smooth(method = "lm", formula = formula) +
- stat_fit_glance(method = "lm",
- method.args = list(formula = formula),
- label.x = "right",
- label.y = "bottom",
- aes(label = sprintf("italic(P)*\"-value = \"*%.3g",
- after_stat(p.value))),
- parse = TRUE)# formula <- y ~ poly(x, 3, raw = TRUE)
+# broom::augment does not handle poly() correctly!
+formula <- y ~ x + I(x^2) + I(x^3)
+ggplot(my.data, aes(x, y, colour = group)) +
+ geom_point() +
+ geom_smooth(method = "lm", formula = formula) +
+ stat_fit_glance(method = "lm",
+ method.args = list(formula = formula),
+ label.x = "right",
+ label.y = "bottom",
+ aes(label = sprintf("italic(P)*\"-value = \"*%.3g",
+ after_stat(p.value))),
+ parse = TRUE)
It is also possible to fit a non-linear model with
method = "nls", and any other model for which a
broom::glance() method exists. Do consult the documentation
for package ‘broom’. Here we fit the Michaelis-Menten equation to
reaction rate versus concentration data.
micmen.formula <- y ~ SSmicmen(x, Vm, K)
-ggplot(Puromycin, aes(conc, rate, colour = state)) +
- geom_point() +
- geom_smooth(method = "nls",
- formula = micmen.formula,
- se = FALSE) +
- stat_fit_glance(method = "nls",
- method.args = list(formula = micmen.formula),
- aes(label = paste("AIC = ", signif(after_stat(AIC), digits = 3),
- ", BIC = ", signif(after_stat(BIC), digits = 3),
- sep = "")),
- label.x = "centre", label.y = "bottom")micmen.formula <- y ~ SSmicmen(x, Vm, K)
+ggplot(Puromycin, aes(conc, rate, colour = state)) +
+ geom_point() +
+ geom_smooth(method = "nls",
+ formula = micmen.formula,
+ se = FALSE) +
+ stat_fit_glance(method = "nls",
+ method.args = list(formula = micmen.formula),
+ aes(label = paste("AIC = ", signif(after_stat(AIC), digits = 3),
+ ", BIC = ", signif(after_stat(BIC), digits = 3),
+ sep = "")),
+ label.x = "centre", label.y = "bottom")
@@ -2540,106 +2569,106 @@
+
+
+
+
+
+
+
stat_fit_tb
The default output of stat_fit_tb() is the default
output from tidy(fm) where fm is the fitted
model.
formula <- y ~ x + I(x^2) + I(x^3)
-ggplot(my.data, aes(x, y)) +
- geom_point() +
- geom_smooth(method = "lm", formula = formula) +
- stat_fit_tb(method = "lm",
- method.args = list(formula = formula),
- tb.vars = c(Parameter = "term",
- Estimate = "estimate",
- "s.e." = "std.error",
- "italic(t)" = "statistic",
- "italic(P)" = "p.value"),
- label.y = "top", label.x = "left",
- parse = TRUE)formula <- y ~ x + I(x^2) + I(x^3)
+ggplot(my.data, aes(x, y)) +
+ geom_point() +
+ geom_smooth(method = "lm", formula = formula) +
+ stat_fit_tb(method = "lm",
+ method.args = list(formula = formula),
+ tb.vars = c(Parameter = "term",
+ Estimate = "estimate",
+ "s.e." = "std.error",
+ "italic(t)" = "statistic",
+ "italic(P)" = "p.value"),
+ label.y = "top", label.x = "left",
+ parse = TRUE)
When tb.type = "fit.anova" the output returned is that
from tidy(anova(fm)) where fm is the fitted
model. Here we also show how to replace names of columns and terms, and
exclude one column, in this case, the mean squares.
formula <- y ~ x + I(x^2) + I(x^3)
-ggplot(my.data, aes(x, y)) +
- geom_point() +
- geom_smooth(method = "lm", formula = formula) +
- stat_fit_tb(method = "lm",
- method.args = list(formula = formula),
- tb.type = "fit.anova",
- tb.vars = c(Effect = "term",
- df = "df",
- "italic(F)" = "statistic",
- "italic(P)" = "p.value"),
- tb.params = c(x = 1, "x^2" = 2, "x^3" = 3, Resid = 4),
- label.y = "top", label.x = "left",
- parse = TRUE)formula <- y ~ x + I(x^2) + I(x^3)
+ggplot(my.data, aes(x, y)) +
+ geom_point() +
+ geom_smooth(method = "lm", formula = formula) +
+ stat_fit_tb(method = "lm",
+ method.args = list(formula = formula),
+ tb.type = "fit.anova",
+ tb.vars = c(Effect = "term",
+ df = "df",
+ "italic(F)" = "statistic",
+ "italic(P)" = "p.value"),
+ tb.params = c(x = 1, "x^2" = 2, "x^3" = 3, Resid = 4),
+ label.y = "top", label.x = "left",
+ parse = TRUE)
When tb.type = "fit.coefs" the output returned is that
of tidy(fm) after selecting the term and
estimate columns.
formula <- y ~ x + I(x^2) + I(x^3)
-ggplot(my.data, aes(x, y)) +
- geom_point() +
- geom_smooth(method = "lm", formula = formula) +
- stat_fit_tb(method = "lm",
- method.args = list(formula = formula),
- tb.type = "fit.coefs", parse = TRUE,
- label.y = "center", label.x = "left")formula <- y ~ x + I(x^2) + I(x^3)
+ggplot(my.data, aes(x, y)) +
+ geom_point() +
+ geom_smooth(method = "lm", formula = formula) +
+ stat_fit_tb(method = "lm",
+ method.args = list(formula = formula),
+ tb.type = "fit.coefs", parse = TRUE,
+ label.y = "center", label.x = "left")
Faceting works as expected, but grouping is ignored as mentioned
above. In this case, the colour aesthetic is not applied to the text of
the tables. Furthermore, if label.x.npc or
label.y.npc are passed numeric vectors of length > 1,
the corresponding values are obeyed by the different panels.
micmen.formula <- y ~ SSmicmen(x, Vm, K)
-ggplot(Puromycin, aes(conc, rate, colour = state)) +
- facet_wrap(~state) +
- geom_point() +
- geom_smooth(method = "nls",
- formula = micmen.formula,
- se = FALSE) +
- stat_fit_tb(method = "nls",
- method.args = list(formula = micmen.formula),
- tb.type = "fit.coefs",
- label.x = 0.9,
- label.y = c(0.75, 0.2)) +
- theme(legend.position = "none") +
- labs(x = "C", y = "V")micmen.formula <- y ~ SSmicmen(x, Vm, K)
+ggplot(Puromycin, aes(conc, rate, colour = state)) +
+ facet_wrap(~state) +
+ geom_point() +
+ geom_smooth(method = "nls",
+ formula = micmen.formula,
+ se = FALSE) +
+ stat_fit_tb(method = "nls",
+ method.args = list(formula = micmen.formula),
+ tb.type = "fit.coefs",
+ label.x = 0.9,
+ label.y = c(0.75, 0.2)) +
+ theme(legend.position = "none") +
+ labs(x = "C", y = "V")
The data in the example below are split by ggplot into six groups
based on the levels of the feed factor. However, as
stat_fit_tb() ignores groupings, we can still fit a linear
model to all the data in the panel.
ggplot(chickwts, aes(factor(feed), weight)) +
- stat_summary(fun.data = "mean_se") +
- stat_fit_tb(tb.type = "fit.anova",
- label.x = "center",
- label.y = "bottom") +
- expand_limits(y = 0)ggplot(chickwts, aes(factor(feed), weight)) +
+ stat_summary(fun.data = "mean_se") +
+ stat_fit_tb(tb.type = "fit.anova",
+ label.x = "center",
+ label.y = "bottom") +
+ expand_limits(y = 0)
We can flip the system of coordinates, if desired.
-ggplot(chickwts, aes(factor(feed), weight)) +
- stat_summary(fun.data = "mean_se") +
- stat_fit_tb(tb.type = "fit.anova", label.x = "left", size = 3) +
- scale_x_discrete(expand = expansion(mult = c(0.2, 0.5))) +
- coord_flip()ggplot(chickwts, aes(factor(feed), weight)) +
+ stat_summary(fun.data = "mean_se") +
+ stat_fit_tb(tb.type = "fit.anova", label.x = "left", size = 3) +
+ scale_x_discrete(expand = expansion(mult = c(0.2, 0.5))) +
+ coord_flip()
It is also possible to rotate the table using angle.
Here we also show how to replace the column headers with strings to be
parsed into R expressions.
ggplot(chickwts, aes(factor(feed), weight)) +
- stat_summary(fun.data = "mean_se") +
- stat_fit_tb(tb.type = "fit.anova",
- angle = 90, size = 3,
- label.x = "right", label.y = "center",
- hjust = 0.5, vjust = 0,
- tb.vars = c(Effect = "term",
- "df",
- "M.S." = "meansq",
- "italic(F)" = "statistic",
- "italic(P)" = "p.value"),
- parse = TRUE) +
- scale_x_discrete(expand = expansion(mult = c(0.1, 0.35))) +
- expand_limits(y = 0)ggplot(chickwts, aes(factor(feed), weight)) +
+ stat_summary(fun.data = "mean_se") +
+ stat_fit_tb(tb.type = "fit.anova",
+ angle = 90, size = 3,
+ label.x = "right", label.y = "center",
+ hjust = 0.5, vjust = 0,
+ tb.vars = c(Effect = "term",
+ "df",
+ "M.S." = "meansq",
+ "italic(F)" = "statistic",
+ "italic(P)" = "p.value"),
+ parse = TRUE) +
+ scale_x_discrete(expand = expansion(mult = c(0.1, 0.35))) +
+ expand_limits(y = 0)
@@ -2657,116 +2686,116 @@
+
+
+
+
+
+
+
stat_fit_tidy
equation of reaction kinetics, usingstats::nls(). We use
the self-starting function stats::SSmicmen() available in
R.
-micmen.formula <- y ~ SSmicmen(x, Vm, K)
-ggplot(Puromycin, aes(conc, rate, colour = state)) +
- geom_point() +
- geom_smooth(method = "nls",
- formula = micmen.formula,
- se = FALSE) +
- stat_fit_tidy(method = "nls",
- method.args = list(formula = micmen.formula),
- label.x = "right",
- label.y = "bottom",
- aes(label = paste("V[m]~`=`~", signif(after_stat(Vm_estimate), digits = 3),
- "%+-%", signif(after_stat(Vm_se), digits = 2),
- "~~~~K~`=`~", signif(after_stat(K_estimate), digits = 3),
- "%+-%", signif(after_stat(K_se), digits = 2),
- sep = "")),
- parse = TRUE)micmen.formula <- y ~ SSmicmen(x, Vm, K)
+ggplot(Puromycin, aes(conc, rate, colour = state)) +
+ geom_point() +
+ geom_smooth(method = "nls",
+ formula = micmen.formula,
+ se = FALSE) +
+ stat_fit_tidy(method = "nls",
+ method.args = list(formula = micmen.formula),
+ label.x = "right",
+ label.y = "bottom",
+ aes(label = paste("V[m]~`=`~", signif(after_stat(Vm_estimate), digits = 3),
+ "%+-%", signif(after_stat(Vm_se), digits = 2),
+ "~~~~K~`=`~", signif(after_stat(K_estimate), digits = 3),
+ "%+-%", signif(after_stat(K_se), digits = 2),
+ sep = "")),
+ parse = TRUE)
Using paste we can build a string that can be parsed into an R expression, in this case for a non-linear equation.
-micmen.formula <- y ~ SSmicmen(x, Vm, K)
-ggplot(Puromycin, aes(conc, rate, colour = state)) +
- geom_point() +
- geom_smooth(method = "nls",
- formula = micmen.formula,
- se = FALSE) +
- stat_fit_tidy(method = "nls",
- method.args = list(formula = micmen.formula),
- size = 3,
- label.x = "center",
- label.y = "bottom",
- vstep = 0.12,
- aes(label = paste("V~`=`~frac(", signif(after_stat(Vm_estimate), digits = 2), "~C,",
- signif(after_stat(K_estimate), digits = 2), "+C)",
- sep = "")),
- parse = TRUE) +
- labs(x = "C", y = "V")micmen.formula <- y ~ SSmicmen(x, Vm, K)
+ggplot(Puromycin, aes(conc, rate, colour = state)) +
+ geom_point() +
+ geom_smooth(method = "nls",
+ formula = micmen.formula,
+ se = FALSE) +
+ stat_fit_tidy(method = "nls",
+ method.args = list(formula = micmen.formula),
+ size = 3,
+ label.x = "center",
+ label.y = "bottom",
+ vstep = 0.12,
+ aes(label = paste("V~`=`~frac(", signif(after_stat(Vm_estimate), digits = 2), "~C,",
+ signif(after_stat(K_estimate), digits = 2), "+C)",
+ sep = "")),
+ parse = TRUE) +
+ labs(x = "C", y = "V")
What if we would need a more specific statistic, similar to
stat_poly_eq()? We can use stat_fit_tidy() as
the basis for its definition.
stat_micmen_eq <- function(vstep = 0.12,
- size = 3,
- ...) {
- stat_fit_tidy(method = "nls",
- method.args = list(formula = micmen.formula),
- aes(label = paste("V~`=`~frac(", signif(after_stat(Vm_estimate), digits = 2), "~C,",
- signif(after_stat(K_estimate), digits = 2), "+C)",
- sep = "")),
- parse = TRUE,
- vstep = vstep,
- size = size,
- ...)
-}stat_micmen_eq <- function(vstep = 0.12,
+ size = 3,
+ ...) {
+ stat_fit_tidy(method = "nls",
+ method.args = list(formula = micmen.formula),
+ aes(label = paste("V~`=`~frac(", signif(after_stat(Vm_estimate), digits = 2), "~C,",
+ signif(after_stat(K_estimate), digits = 2), "+C)",
+ sep = "")),
+ parse = TRUE,
+ vstep = vstep,
+ size = size,
+ ...)
+}The code for the figure is now simpler, and still produces the same figure (not shown).
-micmen.formula <- y ~ SSmicmen(x, Vm, K)
-ggplot(Puromycin, aes(conc, rate, colour = state)) +
- geom_point() +
- geom_smooth(method = "nls",
- formula = micmen.formula,
- se = FALSE) +
- stat_micmen_eq(label.x = "center",
- label.y = "bottom") +
- labs(x = "C", y = "V")micmen.formula <- y ~ SSmicmen(x, Vm, K)
+ggplot(Puromycin, aes(conc, rate, colour = state)) +
+ geom_point() +
+ geom_smooth(method = "nls",
+ formula = micmen.formula,
+ se = FALSE) +
+ stat_micmen_eq(label.x = "center",
+ label.y = "bottom") +
+ labs(x = "C", y = "V")
- As a second example we show a quantile regression fit using function
quantreg::rq()from package ‘quantreg’.
my_formula <- y ~ x
-
-ggplot(mpg, aes(displ, 1 / hwy)) +
- geom_point() +
- geom_quantile(quantiles = 0.5, formula = my_formula) +
- stat_fit_tidy(method = "rq",
- method.args = list(formula = y ~ x, tau = 0.5),
- tidy.args = list(se.type = "nid"),
- mapping = aes(label = sprintf('y~"="~%.3g+%.3g~x*", with "*italic(P)~"="~%.3f',
- after_stat(Intercept_estimate),
- after_stat(x_estimate),
- after_stat(x_p.value))),
- parse = TRUE)my_formula <- y ~ x
+
+ggplot(mpg, aes(displ, 1 / hwy)) +
+ geom_point() +
+ geom_quantile(quantiles = 0.5, formula = my_formula) +
+ stat_fit_tidy(method = "rq",
+ method.args = list(formula = y ~ x, tau = 0.5),
+ tidy.args = list(se.type = "nid"),
+ mapping = aes(label = sprintf('y~"="~%.3g+%.3g~x*", with "*italic(P)~"="~%.3f',
+ after_stat(Intercept_estimate),
+ after_stat(x_estimate),
+ after_stat(x_p.value))),
+ parse = TRUE)
We can define a stat_rq_eq() if we need to add similar
equations to several plots. In this example we retain the ability of the
user to override most of the default default arguments.
stat_rq_eqn <-
- function(formula = y ~ x,
- tau = 0.5,
- method = "br",
- mapping = aes(label = sprintf('y~"="~%.3g+%.3g~x*", with "*italic(P)~"="~%.3f',
- after_stat(Intercept_estimate),
- after_stat(x_estimate),
- after_stat(x_p.value))),
- parse = TRUE,
- ...) {
- method.args <- list(formula = formula, tau = tau, method = method)
- stat_fit_tidy(method = "rq",
- method.args = method.args,
- tidy.args = list(se.type = "nid"),
- mapping = mapping,
- parse = parse,
- ...)
- }stat_rq_eqn <-
+ function(formula = y ~ x,
+ tau = 0.5,
+ method = "br",
+ mapping = aes(label = sprintf('y~"="~%.3g+%.3g~x*", with "*italic(P)~"="~%.3f',
+ after_stat(Intercept_estimate),
+ after_stat(x_estimate),
+ after_stat(x_p.value))),
+ parse = TRUE,
+ ...) {
+ method.args <- list(formula = formula, tau = tau, method = method)
+ stat_fit_tidy(method = "rq",
+ method.args = method.args,
+ tidy.args = list(se.type = "nid"),
+ mapping = mapping,
+ parse = parse,
+ ...)
+ }And the code of the figure now as simple as. Figure not shown, as is identical to the one above.
-ggplot(mpg, aes(displ, 1 / hwy)) +
- geom_point() +
- geom_quantile(quantiles = 0.5, formula = my_formula) +
- stat_rq_eqn(tau = 0.5, formula = my_formula)ggplot(mpg, aes(displ, 1 / hwy)) +
+ geom_point() +
+ geom_quantile(quantiles = 0.5, formula = my_formula) +
+ stat_rq_eqn(tau = 0.5, formula = my_formula)
@@ -2781,66 +2810,66 @@
+
+
+
+
+
+
+
stat_fit_augment
stat_fit_augment() can handle any fitted model that is
supported by package ‘broom’ or its extensions. All these statistics
support grouping and facets.
-# formula <- y ~ poly(x, 3, raw = TRUE)
-# broom::augment does not handle poly correctly!
-formula <- y ~ x + I(x^2) + I(x^3)
-ggplot(my.data, aes(x, y)) +
- geom_point() +
- stat_fit_augment(method = "lm",
- method.args = list(formula = formula))# formula <- y ~ poly(x, 3, raw = TRUE)
+# broom::augment does not handle poly correctly!
+formula <- y ~ x + I(x^2) + I(x^3)
+ggplot(my.data, aes(x, y)) +
+ geom_point() +
+ stat_fit_augment(method = "lm",
+ method.args = list(formula = formula))
formula <- y ~ x + I(x^2) + I(x^3)
-ggplot(my.data, aes(x, y, colour = group)) +
- geom_point() +
- stat_fit_augment(method = "lm",
- method.args = list(formula = formula))formula <- y ~ x + I(x^2) + I(x^3)
+ggplot(my.data, aes(x, y, colour = group)) +
+ geom_point() +
+ stat_fit_augment(method = "lm",
+ method.args = list(formula = formula))
We can override the variable returned as y to be any of
the variables in the data frame returned by
broom::augment() while still preserving the original
y values.
formula <- y ~ x + I(x^2) + I(x^3)
-ggplot(my.data, aes(x, y)) +
- stat_fit_augment(method = "lm",
- method.args = list(formula = formula),
- geom = "point",
- y.out = ".resid")formula <- y ~ x + I(x^2) + I(x^3)
+ggplot(my.data, aes(x, y)) +
+ stat_fit_augment(method = "lm",
+ method.args = list(formula = formula),
+ geom = "point",
+ y.out = ".resid")
formula <- y ~ x + I(x^2) + I(x^3)
-ggplot(my.data, aes(x, y, colour = group)) +
- stat_fit_augment(method = "lm",
- method.args = list(formula = formula),
- geom = "point",
- y.out = ".std.resid")formula <- y ~ x + I(x^2) + I(x^3)
+ggplot(my.data, aes(x, y, colour = group)) +
+ stat_fit_augment(method = "lm",
+ method.args = list(formula = formula),
+ geom = "point",
+ y.out = ".std.resid")
We can use any model fitting method for which
broom::augment() is implemented.
args <- list(formula = y ~ k * e ^ x,
- start = list(k = 1, e = 2))
-ggplot(mtcars, aes(wt, mpg)) +
- geom_point() +
- stat_fit_augment(method = "nls",
- method.args = args)args <- list(formula = y ~ k * e ^ x,
+ start = list(k = 1, e = 2))
+ggplot(mtcars, aes(wt, mpg)) +
+ geom_point() +
+ stat_fit_augment(method = "nls",
+ method.args = args)
args <- list(formula = y ~ k * e ^ x,
- start = list(k = 1, e = 2))
-ggplot(mtcars, aes(wt, mpg)) +
- stat_fit_augment(method = "nls",
- method.args = args,
- geom = "point",
- y.out = ".resid")args <- list(formula = y ~ k * e ^ x,
+ start = list(k = 1, e = 2))
+ggplot(mtcars, aes(wt, mpg)) +
+ stat_fit_augment(method = "nls",
+ method.args = args,
+ geom = "point",
+ y.out = ".resid")
Note: The tidiers for mixed models have moved to package ‘broom.mixed’!
-args <- list(model = y ~ SSlogis(x, Asym, xmid, scal),
- fixed = Asym + xmid + scal ~1,
- random = Asym ~1 | group,
- start = c(Asym = 200, xmid = 725, scal = 350))
-ggplot(Orange, aes(age, circumference, colour = Tree)) +
- geom_point() +
- stat_fit_augment(method = "nlme",
- method.args = args,
- augment.args = list(data = quote(data)))args <- list(model = y ~ SSlogis(x, Asym, xmid, scal),
+ fixed = Asym + xmid + scal ~1,
+ random = Asym ~1 | group,
+ start = c(Asym = 200, xmid = 725, scal = 350))
+ggplot(Orange, aes(age, circumference, colour = Tree)) +
+ geom_point() +
+ stat_fit_augment(method = "nlme",
+ method.args = args,
+ augment.args = list(data = quote(data)))
@@ -2860,20 +2889,20 @@
+
How to find out what variables are computed
statistics or directly as a layer function. However, unrecognised aesthetics can trigger warnings during plot rendering. These warnings can be ignored in most cases. -# formula <- y ~ poly(x, 3, raw = TRUE)
-# broom::augment does not handle poly() correctly!
-formula <- y ~ x + I(x^2) + I(x^3)
-ggplot(my.data, aes(x, y, colour = group)) +
- geom_point() +
- geom_smooth(method = "lm", formula = formula) +
- stat_fit_glance(geom = "debug",
- method = "lm",
- method.args = list(formula = formula),
- label.x = "right",
- label.y = "bottom",
- aes(label = sprintf("italic(P)*\"-value = \"*%.3g",
- after_stat(p.value))),
- parse = TRUE)# formula <- y ~ poly(x, 3, raw = TRUE)
+# broom::augment does not handle poly() correctly!
+formula <- y ~ x + I(x^2) + I(x^3)
+ggplot(my.data, aes(x, y, colour = group)) +
+ geom_point() +
+ geom_smooth(method = "lm", formula = formula) +
+ stat_fit_glance(geom = "debug",
+ method = "lm",
+ method.args = list(formula = formula),
+ label.x = "right",
+ label.y = "bottom",
+ aes(label = sprintf("italic(P)*\"-value = \"*%.3g",
+ after_stat(p.value))),
+ parse = TRUE)## Warning in stat_fit_glance(geom = "debug", method = "lm", method.args =
## list(formula = formula), : Ignoring unknown parameters: `parse`## [1] "PANEL 1; group(s) 1, 2; 'draw_function()' input 'data' (head):"
@@ -2895,21 +2924,21 @@ How to find out what variables are computed
case, function str() is more informative than
head(), so we can pass it as argument overriding the
default.
-formula <- y ~ x + I(x^2) + I(x^3)
-ggplot(my.data, aes(x, y)) +
- geom_point() +
- geom_smooth(method = "lm", formula = formula) +
- stat_fit_tb(geom = "debug",
- summary.fun = str,
- method = "lm",
- method.args = list(formula = formula),
- tb.vars = c(Parameter = "term",
- Estimate = "estimate",
- "s.e." = "std.error",
- "italic(t)" = "statistic",
- "italic(P)" = "p.value"),
- label.y = "top", label.x = "left",
- parse = TRUE)
+formula <- y ~ x + I(x^2) + I(x^3)
+ggplot(my.data, aes(x, y)) +
+ geom_point() +
+ geom_smooth(method = "lm", formula = formula) +
+ stat_fit_tb(geom = "debug",
+ summary.fun = str,
+ method = "lm",
+ method.args = list(formula = formula),
+ tb.vars = c(Parameter = "term",
+ Estimate = "estimate",
+ "s.e." = "std.error",
+ "italic(t)" = "statistic",
+ "italic(P)" = "p.value"),
+ label.y = "top", label.x = "left",
+ parse = TRUE)
## Warning in stat_fit_tb(geom = "debug", summary.fun = str, method = "lm", : Ignoring unknown parameters: `table.theme`, `table.rownames`, `table.colnames`,
## `table.hjust`, `parse`, and `summary.fun`
## [1] "PANEL 1; group(s) NULL; 'draw_function()' input 'data' (head):"
@@ -2970,7 +2999,7 @@ Volcano-plot examples
Outcomes encoded as -1, 0 or 1, as seen in the tibble below need to
be converted into factors with suitable labels for levels. This can be
easily achieved with function outcome2factor().
-
+
## tag gene outcome logFC PValue genotype
## 1 AT1G01040 ASU1 0 -0.15284466 0.35266997 Ler
## 2 AT1G01290 ASG4 0 -0.30057068 0.05471732 Ler
@@ -2996,24 +3025,24 @@ Volcano-plot examples
arguments to parameters name and labels of
these scales. These x and y scales by default squish off-limits
(out-of-bounds) observations towards the limits.
-ggplot(volcano_example.df,
- aes(logFC, PValue, colour = outcome2factor(outcome))) +
- geom_point() +
- scale_x_logFC(name = "Transcript abundance%unit") +
- scale_y_Pvalue() +
- scale_colour_outcome() +
- stat_quadrant_counts(data = function(x) {subset(x, outcome != 0)})
+ggplot(volcano_example.df,
+ aes(logFC, PValue, colour = outcome2factor(outcome))) +
+ geom_point() +
+ scale_x_logFC(name = "Transcript abundance%unit") +
+ scale_y_Pvalue() +
+ scale_colour_outcome() +
+ stat_quadrant_counts(data = function(x) {subset(x, outcome != 0)})
By default outcome2factor() creates a factor with three
levels as in the example above, but this default can be overridden as
shown below.
-ggplot(volcano_example.df,
- aes(logFC, PValue, colour = outcome2factor(outcome, n.levels = 2))) +
- geom_point() +
- scale_x_logFC(name = "Transcript abundance%unit") +
- scale_y_Pvalue() +
- scale_colour_outcome() +
- stat_quadrant_counts(data = function(x) {subset(x, outcome != 0)})
+ggplot(volcano_example.df,
+ aes(logFC, PValue, colour = outcome2factor(outcome, n.levels = 2))) +
+ geom_point() +
+ scale_x_logFC(name = "Transcript abundance%unit") +
+ scale_y_Pvalue() +
+ scale_colour_outcome() +
+ stat_quadrant_counts(data = function(x) {subset(x, outcome != 0)})
@@ -3032,7 +3061,7 @@
+
-
+
+
+
Quadrant-plot examples
circle and map one of the outcomes to colour and the other to fill, using the two matched scalesscale_colour_outcome() and
scale_fill_outcome().
-
+
## tag gene outcome.x outcome.y logFC.x logFC.y genotype
## 1 AT5G11060 TIC56 0 0 -0.17685517 -0.32956762 Ler
## 2 AT3G01280 ATWRKY48 0 0 -0.06471884 0.07771315 Ler
@@ -3042,66 +3071,66 @@ Quadrant-plot examples
## 6 AT2G16070 UBQ11 0 0 -0.22328946 -0.23210780 LerIn this plot we do not include those genes whose change in transcript abundance is uncertain under both x and y conditions.
- ggplot(subset(quadrant_example.df,
- xy_outcomes2factor(outcome.x, outcome.y) != "none"),
- aes(logFC.x, logFC.y,
- colour = outcome2factor(outcome.x),
- fill = outcome2factor(outcome.y))) +
- geom_quadrant_lines(linetype = "dotted") +
- stat_quadrant_counts(size = 3, colour = "white") +
- geom_point(shape = "circle filled") +
- scale_x_logFC(name = "Transcript abundance for x%unit") +
- scale_y_logFC(name = "Transcript abundance for y%unit") +
- scale_colour_outcome() +
- scale_fill_outcome() +
- theme_dark() ggplot(subset(quadrant_example.df,
+ xy_outcomes2factor(outcome.x, outcome.y) != "none"),
+ aes(logFC.x, logFC.y,
+ colour = outcome2factor(outcome.x),
+ fill = outcome2factor(outcome.y))) +
+ geom_quadrant_lines(linetype = "dotted") +
+ stat_quadrant_counts(size = 3, colour = "white") +
+ geom_point(shape = "circle filled") +
+ scale_x_logFC(name = "Transcript abundance for x%unit") +
+ scale_y_logFC(name = "Transcript abundance for y%unit") +
+ scale_colour_outcome() +
+ scale_fill_outcome() +
+ theme_dark()
To plot in separate panels those observations that are significant along both x and y axes, x axis, y axis, or none, with quadrants merged takes more effort. We first define two helper functions to add counts and quadrant lines to each of the four panels.
-all_quadrant_counts <- function(...) {
- list(
- stat_quadrant_counts(data = . %>% filter(outcome.xy.fct == "xy"), ...),
- stat_quadrant_counts(data = . %>% filter(outcome.xy.fct == "x"), pool.along = "y", ...),
- stat_quadrant_counts(data = . %>% filter(outcome.xy.fct == "y"), pool.along = "x", ...),
- stat_quadrant_counts(data = . %>% filter(outcome.xy.fct == "none"), quadrants = 0L, ...)
- )
-}all_quadrant_lines <- function(...) {
- list(
- geom_hline(data = data.frame(outcome.xy.fct = factor(c("xy", "x", "y", "none"),
- levels = c("xy", "x", "y", "none")),
- yintercept = c(0, NA, 0, NA)),
- aes(yintercept = yintercept),
- na.rm = TRUE,
- ...),
- geom_vline(data = data.frame(outcome.xy.fct = factor(c("xy", "x", "y", "none"),
- levels = c("xy", "x", "y", "none")),
- xintercept = c(0, 0, NA, NA)),
- aes(xintercept = xintercept),
- na.rm = TRUE,
- ...)
- )
-}all_quadrant_counts <- function(...) {
+ list(
+ stat_quadrant_counts(data = . %>% filter(outcome.xy.fct == "xy"), ...),
+ stat_quadrant_counts(data = . %>% filter(outcome.xy.fct == "x"), pool.along = "y", ...),
+ stat_quadrant_counts(data = . %>% filter(outcome.xy.fct == "y"), pool.along = "x", ...),
+ stat_quadrant_counts(data = . %>% filter(outcome.xy.fct == "none"), quadrants = 0L, ...)
+ )
+}all_quadrant_lines <- function(...) {
+ list(
+ geom_hline(data = data.frame(outcome.xy.fct = factor(c("xy", "x", "y", "none"),
+ levels = c("xy", "x", "y", "none")),
+ yintercept = c(0, NA, 0, NA)),
+ aes(yintercept = yintercept),
+ na.rm = TRUE,
+ ...),
+ geom_vline(data = data.frame(outcome.xy.fct = factor(c("xy", "x", "y", "none"),
+ levels = c("xy", "x", "y", "none")),
+ xintercept = c(0, 0, NA, NA)),
+ aes(xintercept = xintercept),
+ na.rm = TRUE,
+ ...)
+ )
+}And use these functions to build the final plot, in this case including all genes.
-quadrant_example.df %>%
- mutate(.,
- outcome.x.fct = outcome2factor(outcome.x),
- outcome.y.fct = outcome2factor(outcome.y),
- outcome.xy.fct = xy_outcomes2factor(outcome.x, outcome.y)) %>%
- ggplot(., aes(logFC.x, logFC.y, colour = outcome.x.fct, fill = outcome.y.fct)) +
- geom_point(shape = 21) +
- all_quadrant_lines(linetype = "dotted") +
- all_quadrant_counts(size = 3, colour = "white") +
- scale_x_logFC(name = "Transcript abundance for x%unit") +
- scale_y_logFC(name = "Transcript abundance for y%unit") +
- scale_colour_outcome() +
- scale_fill_outcome() +
- facet_wrap(~outcome.xy.fct) +
- theme_dark()quadrant_example.df %>%
+ mutate(.,
+ outcome.x.fct = outcome2factor(outcome.x),
+ outcome.y.fct = outcome2factor(outcome.y),
+ outcome.xy.fct = xy_outcomes2factor(outcome.x, outcome.y)) %>%
+ ggplot(., aes(logFC.x, logFC.y, colour = outcome.x.fct, fill = outcome.y.fct)) +
+ geom_point(shape = 21) +
+ all_quadrant_lines(linetype = "dotted") +
+ all_quadrant_counts(size = 3, colour = "white") +
+ scale_x_logFC(name = "Transcript abundance for x%unit") +
+ scale_y_logFC(name = "Transcript abundance for y%unit") +
+ scale_colour_outcome() +
+ scale_fill_outcome() +
+ facet_wrap(~outcome.xy.fct) +
+ theme_dark()