Stéphane Laurent
Stéphane Laurent

Reputation: 84599

Coefficients per facet with output.type="numeric" in ggpmisc::stat_poly_eq

ggpmisc::stat_poly_eq has an option output.type = "numeric" allowing to get the estimates of the parameters of the fitted model. Below is my attempt to use it with facet_wrap. I get a different per facet but the coefficients are the same in the two facets. Do I do something wrong, or is it a bug?

library(ggpmisc)

set.seed(4321)
x <- 1:100
y <- (x + x^2 + x^3) + rnorm(length(x), mean = 0, sd = mean(x^3) / 4)
my.data <- data.frame(x = x, 
                      y = y,
                      group = c("A", "B"))
my.data[my.data$group=="A",]$y <- my.data[my.data$group=="A",]$y + 200000

formula <- y ~ poly(x, 1, raw = TRUE)

myformat <- "Intercept: %s\nSlope: %s\nR²: %s"
ggplot(my.data, aes(x, y)) + 
  facet_wrap(~ group) + 
  geom_point() +
  geom_smooth(method = "lm", formula = formula) +
  stat_poly_eq(formula = formula, output.type = "numeric",
               mapping = aes(label = 
                               sprintf(myformat,
                                       formatC(stat(coef.ls)[[1]][[1, "Estimate"]]),
                                       formatC(stat(coef.ls)[[1]][[2, "Estimate"]]),
                                       formatC(stat(r.squared))))) 

enter image description here


Edit

We have to catch the panel number. It is strange that formatC(stat(as.integer(PANEL))) returns the panel number per facet:

enter image description here

but however formatC(stat(coef.ls)[[stat(as.integer(PANEL))]][[1, "Estimate"]]) does not work, because here PANEL = c(1,2).

Upvotes: 1

Views: 989

Answers (2)

Pedro J. Aphalo
Pedro J. Aphalo

Reputation: 6528

Version 0.3.2 of 'ggpmisc' is now in CRAN. Submitted earlier this week. In the documentation I now give some examples of the use of geom_debug() from my package 'gginnards' to have a look at the data frame returned by stats (usable with any ggplot stat or by itself). For your example, it would work like this:

library(ggpmisc)
library(gginnards)

set.seed(4321)
x <- 1:100
y <- (x + x^2 + x^3) + rnorm(length(x), mean = 0, sd = mean(x^3) / 4)
my.data <- data.frame(x = x, 
                      y = y,
                      group = c("A", "B"))
my.data[my.data$group=="A",]$y <- my.data[my.data$group=="A",]$y + 200000

formula <- y ~ poly(x, 1, raw = TRUE)

myformat <- "Intercept: %s\nSlope: %s\nR²: %s"
ggplot(my.data, aes(x, y)) + 
  facet_wrap(~ group) + 
  geom_point() +
  geom_smooth(method = "lm", formula = formula) +
  stat_poly_eq(formula = formula, output.type = "numeric",
               aes(label = ""),
               geom = "debug") 

Which prints to the console, two tibbles, one for each panel: enter image description here

Example below added to address comment:

ggplot(my.data, aes(x, y)) + 
  facet_wrap(~ group) + 
  geom_point() +
  geom_smooth(method = "lm", formula = formula) +
  stat_poly_eq(formula = formula, output.type = "numeric",
               aes(label = ""),
               summary.fun = function(x) {x[["coef.ls"]][[1]]})

prints just the coefs.ls.

enter image description here

I added the "numeric" option recently in response to a suggestion and with this example I noticed a bug: aes(label = "") should not have been needed, but is needed because the default mapping for the label aesthetic is wrong. I will fix this for the next release.

Upvotes: 0

St&#233;phane Laurent
St&#233;phane Laurent

Reputation: 84599

Ok, I figured it out.

ggplot(my.data, aes(x, y)) + 
  facet_wrap(~ group) + 
  geom_point() +
  geom_smooth(method = "lm", formula = formula) +
  stat_poly_eq(
    formula = formula, output.type = "numeric",
    mapping = aes(label = 
                    sprintf(myformat,
                            c(formatC(stat(coef.ls)[[1]][[1, "Estimate"]]), 
                              formatC(stat(coef.ls)[[2]][[1, "Estimate"]])),
                            c(formatC(stat(coef.ls)[[1]][[2, "Estimate"]]), 
                              formatC(stat(coef.ls)[[2]][[2, "Estimate"]])),
                            formatC(stat(r.squared))))) 

enter image description here

Upvotes: 2

Related Questions