Plotting MA / SMA / RMA regressions lines (smatr or lmodel2 packages) into ggplot

Question

I would be very grateful for advise on how to plot Standardised Major Axis (SMA) regressions lines into a faceted ggplot. I used the following code:

Run the SMA analysis and create a data frame with the SMA reg line coefficients (intercept and slope) I want to plot

smaReg = sma(Y ~ X * Type, data = ExampleData)
summary(smaReg)
smaSummary <- data.frame(Type = 1:6,coef(smaReg))

ggplot code using geom_abline to plot SMA regressions

ModFit <- ggplot(ExampleData, aes(y = Y, x = X, color = Level)) +
  geom_point() +
  theme_bw() +
  theme_classic() + 
  facet_wrap(~ Type, nrow = 2, ncol = 3) +
  theme(strip.background = element_blank(), strip.text = element_text(face = 'bold', size = 12)) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf, color = 'black', size = 1) +
  annotate("segment", x = -Inf, xend = -Inf, y = -Inf, yend = Inf, color = 'black', size = 1) +
  scale_x_continuous(breaks = seq(from = 0, to = 60, by = 20)) + 
  scale_y_continuous(breaks = seq(from = 0, to = 120, by = 20)) +
  geom_abline(data = smaSummary, aes(intercept = elevation, slope = slope)) +
  labs(x = expression(paste("Predicted (",mu,"mol m"^{-2},"s"^{-1},")")), y = expression(paste("Observed (",mu,"mol m"^{-2},"s"^{-1},")"))) +

ModFit

This code has two remaining issues I need to resolve but my beginner coding skills are not good enough yet to tackle them successfully:

1. I used annotate() and scale_x_continuous to plot the same axes and scales across all faceted plots, however, this solution does not plot the X axis ticks and I haven't found a way to do this without something else going wrong when I make a change.

2. When I run this plot code, I get the error message below:

Error in wrap_dims(n, params$nrow, params$ncol) : nrow * ncol >= n is not TRUE

In trying different ways of solving this error, I noticed that if I change the labs() layer to the very simplified version shown below:

labs(x = expression(X), y = expression(Y), color = "Level") +

This change produces a faceted plot but with all the SMA regressions on each plot. I have no idea why changing the labs() layer allows the plot to be produced! I have run out of ideas (and google searches) on how to plot only the corresponding SMA reg line for each plot while also adding the detailed axis labels I need without something else going wrong.

Faceted plot with simplified labels and all SMA reg lines on each plot

Many thanks in advance for any advise on how to solve these two remaining issues!

Answer 1

(UPDATED: Added greek math symbols)

Let's tackle the labels, which are pretty thorny in your example. I'm using the iris dataset here as a sample, but using your axis titles.

The key is using bquote() to get the math, the dynamic variables and everything else ready.

• A lot more detail available here: https://trinkerrstuff.wordpress.com/2018/03/15/2246/
• And on SO: Evaluate expression given as a string

library(tidyverse)
mu_val <- 5.1
bq_x <- bquote("Predicted (" ~ mu ~ "=" ~ .(mu_val) ~ " mol " ~ m^-2 ~ s^-1 ~ ")")
bq_y <- bquote("Observed (" ~ mu ~ "=" ~ .(mu_val) ~ "mol m" ~ m^-2 ~ s^-1 ~ ")")

iris %>%
  ggplot(aes(Sepal.Length, Petal.Length)) +
  geom_point() +
  labs(title = "test", x = bq_x, y = bq_y)

^{Created on 2019-11-19 by the reprex package (v0.3.0)}

@ravic_ thanks to your advise I am now able to produce a faceted plot with the correct labels using the revised code below:

My updated my ggplot code:

bq_x <- bquote("Predicted (" ~mu~"mol" ~ m^-2 ~ s^-1 ~ ")")
bq_y <- bquote("Observed (" ~mu~"mol" ~ m^-2 ~ s^-1 ~ ")")

# 1a. Plots by PFTs
ModFit <- ggplot(ExampleData, aes(y = Y, x = X, color = Level)) +
  geom_point() +
  theme_bw() +
  theme_classic() + 
  facet_wrap(~ Type, nrow = 2, ncol = 3) +
  # customise facet labels
  theme(strip.background = element_blank(), strip.text = element_text(face = 'bold', size = 12)) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf, color = 'black', size = 1) +
  annotate("segment", x = -Inf, xend = -Inf, y = -Inf, yend = Inf, color = 'black', size = 1) +
  scale_x_continuous(breaks = seq(from = 0, to = 60, by = 20)) + 
  scale_y_continuous(breaks = seq(from = 0, to = 120, by = 20)) +
  labs(x = bq_x, y = bq_y) +
  geom_abline(data = smaSummary, aes(intercept = elevation, slope = slope))

The two remaining problems now are: 1. Add Y and X axis lines with ticks across all faceted plots 2. Add the corresponding single SMA reg line to each plot

Many thanks again for all your help @ravic_!

Answer 2

This answer is close, but not exactly an answer to the axis label part of the question. It is, however, a full answer to adding a major axis regression to a facetted ggplot. Trellis plots, from which facets derive, are meant to be a matrix of similar plots and thus share the axis labels. The only way of having different axis labels is to make the plots individually and to compose the final figure from them using, for example, package 'patchwork'. Anyway, in general, annotations derived from the data are usually included in the plotting area or highlighted along the axes. The examples below achieve the highlighting, labelling and adding the MA regression using statistics from package 'ggpmisc' and geometries from package 'ggpp', re-exported by 'ggpmisc'.

library(tidyverse)
library(ggpmisc)
#> Loading required package: ggpp
#> 
#> Attaching package: 'ggpp'
#> The following object is masked from 'package:ggplot2':
#> 
#>     annotate

ggplot(iris, aes(Sepal.Length, Petal.Length)) +
  geom_point() +
  stat_ma_line() +
  stat_centroid(color = "red") +
  stat_centroid(geom = "rug", color = "red") +
  labs(title = "test", 
       x = expression("Observed "*(mu*mol~m^{-2} ~ s^{-1})), 
       y = expression("Predicted "*(mu*mol~m^{-2} ~ s^{-1}))) +
  facet_wrap(~Species)

ggplot(iris, aes(Sepal.Length, Petal.Length)) +
  geom_point() +
  stat_ma_line() +
  stat_centroid(color = "red") +
  stat_centroid(color = "red", geom = "rug") +
  stat_centroid(aes(label = sprintf("bar(x)~`=`~%.1f*\";  \"*bar(y)~`=`~%.1f", 
                                    after_stat(x), after_stat(y))), 
                geom = "text_npc", 
                color = "red", npcy = 0.95, npcx = 0.05,
                hjust = "inward", vjust = "inward",
                parse = TRUE) +
  labs(title = "test", 
       x = expression("Observed "*(mu*mol~m^{-2} ~ s^{-1})), 
       y = expression("Predicted "*(mu*mol~m^{-2} ~ s^{-1}))) +
  facet_wrap(~Species)

ggplot(iris, aes(Sepal.Length, Petal.Length)) +
  geom_point() +
  stat_ma_line() +
  stat_centroid(color = "red") +
  stat_centroid(color = "red", geom = "rug") +
  stat_centroid(aes(label = sprintf("bar(x)~`=`~%.1f*\";  \"*bar(y)~`=`~%.1f", 
                                    after_stat(x), after_stat(y))), 
                geom = "text_npc", 
                color = "red", npcx = 0.95, npcy = 0.05,
                hjust = "inward", vjust = "inward",
                parse = TRUE) +
  labs(title = "test", 
       x = expression("Observed "*(mu*mol~m^{-2} ~ s^{-1})), 
       y = expression("Predicted "*(mu*mol~m^{-2} ~ s^{-1}))) +
  facet_wrap(~Species, scales = "free")

^{Created on 2022-08-15 by the reprex package (v2.0.1)}

stat_ma_line() and its companion stat_ma_eq() are wrappers on lmodel2 and support most of the underlying options. Please, have a look at the documentation for ['ggpmisc'][1] and for ['ggpp'][2] for the details. The examples above, for simplicity use the default "lmodel2:MA", but stat_ma_line() also supports SMA and RMA as methods.

The geoms from 'ggpp' that use native plot coordinates (npc) allow the same code to be used with "free" and "fixed" panel axis scales.