Agregue la ecuación de línea de regresión y R ^ 2 en el gráfico

228

Me pregunto cómo agregar la ecuación de línea de regresión y R ^ 2 en el ggplot. Mi código es:

library(ggplot2)

df <- data.frame(x = c(1:100))
df$y <- 2 + 3 * df$x + rnorm(100, sd = 40)
p <- ggplot(data = df, aes(x = x, y = y)) +
            geom_smooth(method = "lm", se=FALSE, color="black", formula = y ~ x) +
            geom_point()
p

Cualquier ayuda será muy apreciada.

MYaseen208
fuente
1
Para gráficos de celosía , vea latticeExtra::lmlineq().
Josh O'Brien el

Respuestas:

234

Aquí hay una solución

# GET EQUATION AND R-SQUARED AS STRING
# SOURCE: https://groups.google.com/forum/#!topic/ggplot2/1TgH-kG5XMA

lm_eqn <- function(df){
    m <- lm(y ~ x, df);
    eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
         list(a = format(unname(coef(m)[1]), digits = 2),
              b = format(unname(coef(m)[2]), digits = 2),
             r2 = format(summary(m)$r.squared, digits = 3)))
    as.character(as.expression(eq));
}

p1 <- p + geom_text(x = 25, y = 300, label = lm_eqn(df), parse = TRUE)

EDITAR. Descubrí la fuente de donde elegí este código. Aquí está el enlace a la publicación original en los grupos de ggplot2 google

Salida

Ramnath
fuente
1
El comentario de @JonasRaedle sobre cómo obtener textos con mejor aspecto annotateera correcto en mi máquina.
IRTFM
2
Esto no se parece en nada al resultado publicado en mi máquina, donde la etiqueta se sobrescribe tantas veces como se llaman los datos, lo que resulta en un texto de etiqueta grueso y borroso. Pasar las etiquetas a un data.frame funciona primero (vea mi sugerencia en un comentario a continuación.
PatrickT
@PatrickT: elimine el aes(y el correspondiente ). aeses para asignar variables de marco de datos a variables visuales; eso no es necesario aquí, ya que solo hay una instancia, por lo que puede ponerlo todo en la geom_textllamada principal . Editaré esto en la respuesta.
naught101
El problema con esta solución parece ser que si el conjunto de datos es más grande (el mío fue de 370000 observaciones), la función parece fallar. Recomendaría la solución de @kdauria que hace lo mismo, pero mucho más rápido.
Benjamin
3
para aquellos que desean valores r y p en lugar de R2 y ecuación: eq <- sustituto (cursiva (r) ~ "=" ~ rvalue * "," ~ cursiva (p) ~ "=" ~ pvalue, list (rvalue = sprintf ("% .2f", signo (coef (m) [2]) * sqrt (resumen (m) $ r.squared)), pvalue = formato (resumen (m) $ coeficientes [2,4], dígitos = 2 )))
Jerry T
135

Incluí una estadística stat_poly_eq()en mi paquete ggpmiscque permite esta respuesta:

library(ggplot2)
library(ggpmisc)
df <- data.frame(x = c(1:100))
df$y <- 2 + 3 * df$x + rnorm(100, sd = 40)
my.formula <- y ~ x
p <- ggplot(data = df, aes(x = x, y = y)) +
   geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula) +
   stat_poly_eq(formula = my.formula, 
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +         
   geom_point()
p

ingrese la descripción de la imagen aquí

Esta estadística funciona con cualquier polinomio sin términos faltantes, y con suerte tiene suficiente flexibilidad para ser generalmente útil. Las etiquetas R ^ 2 o R ^ 2 ajustadas se pueden usar con cualquier fórmula modelo equipada con lm (). Al ser una estadística ggplot, se comporta como se esperaba tanto con grupos como con facetas.

El paquete 'ggpmisc' está disponible a través de CRAN.

La versión 0.2.6 acaba de ser aceptada en CRAN.

Aborda los comentarios de @shabbychef y @ MYaseen208.

@ MYaseen208 esto muestra cómo agregar un sombrero .

library(ggplot2)
library(ggpmisc)
df <- data.frame(x = c(1:100))
df$y <- 2 + 3 * df$x + rnorm(100, sd = 40)
my.formula <- y ~ x
p <- ggplot(data = df, aes(x = x, y = y)) +
   geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula) +
   stat_poly_eq(formula = my.formula,
                eq.with.lhs = "italic(hat(y))~`=`~",
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE) +         
   geom_point()
p

ingrese la descripción de la imagen aquí

@shabbychef Ahora es posible hacer coincidir las variables de la ecuación con las utilizadas para las etiquetas de eje. Para reemplazar la x con digamos z e y con h, se usaría:

p <- ggplot(data = df, aes(x = x, y = y)) +
   geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula) +
   stat_poly_eq(formula = my.formula,
                eq.with.lhs = "italic(h)~`=`~",
                eq.x.rhs = "~italic(z)",
                aes(label = ..eq.label..), 
                parse = TRUE) + 
   labs(x = expression(italic(z)), y = expression(italic(h))) +          
   geom_point()
p

ingrese la descripción de la imagen aquí

Al ser estas expresiones analizadas en R normales, las letras griegas ahora también se pueden usar tanto en lhs como en rhs de la ecuación.

[2017-03-08] @elarry Edit para abordar con mayor precisión la pregunta original, que muestra cómo agregar una coma entre las etiquetas de ecuación y R2.

p <- ggplot(data = df, aes(x = x, y = y)) +
  geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula) +
  stat_poly_eq(formula = my.formula,
               eq.with.lhs = "italic(hat(y))~`=`~",
               aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~")), 
               parse = TRUE) +         
  geom_point()
p

ingrese la descripción de la imagen aquí

[2019-10-20] @ helen.h A continuación, doy ejemplos de uso de stat_poly_eq()agrupación.

library(ggpmisc)
df <- data.frame(x = c(1:100))
df$y <- 20 * c(0, 1) + 3 * df$x + rnorm(100, sd = 40)
df$group <- factor(rep(c("A", "B"), 50))
my.formula <- y ~ x
p <- ggplot(data = df, aes(x = x, y = y, colour = group)) +
  geom_smooth(method = "lm", se=FALSE, formula = my.formula) +
  stat_poly_eq(formula = my.formula, 
               aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               parse = TRUE) +         
  geom_point()
p

p <- ggplot(data = df, aes(x = x, y = y, linetype = group)) +
  geom_smooth(method = "lm", se=FALSE, formula = my.formula) +
  stat_poly_eq(formula = my.formula, 
               aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               parse = TRUE) +         
  geom_point()
p

ingrese la descripción de la imagen aquí

ingrese la descripción de la imagen aquí

[2020-01-21] @Herman Puede ser un poco contraintuitivo a primera vista, pero para obtener una ecuación única cuando se usa la agrupación, uno debe seguir la gramática de los gráficos. Restrinja la asignación que crea la agrupación a capas individuales (que se muestra a continuación) o mantenga la asignación predeterminada y anúlela con un valor constante en la capa donde no desea la agrupación (por ejemplo colour = "black").

Continuando del ejemplo anterior.

p <- ggplot(data = df, aes(x = x, y = y)) +
  geom_smooth(method = "lm", se=FALSE, formula = my.formula) +
  stat_poly_eq(formula = my.formula, 
               aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               parse = TRUE) +         
  geom_point(aes(colour = group))
p

ingrese la descripción de la imagen aquí

[2020-01-22] En aras de la exhaustividad, un ejemplo con facetas, que demuestra que también en este caso se cumplen las expectativas de la gramática de los gráficos.

library(ggpmisc)
df <- data.frame(x = c(1:100))
df$y <- 20 * c(0, 1) + 3 * df$x + rnorm(100, sd = 40)
df$group <- factor(rep(c("A", "B"), 50))
my.formula <- y ~ x

p <- ggplot(data = df, aes(x = x, y = y)) +
  geom_smooth(method = "lm", se=FALSE, formula = my.formula) +
  stat_poly_eq(formula = my.formula, 
               aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               parse = TRUE) +         
  geom_point() +
  facet_wrap(~group)
p

ingrese la descripción de la imagen aquí

Pedro Aphalo
fuente
1
Cabe señalar que xy yen la fórmula se refieren a los datos xy yen las capas de la trama, y ​​no necesariamente a aquellos en el alcance en el momento en que my.formulase construye. Entonces, ¿la fórmula siempre debe usar las variables x e y?
shabbychef
Es muy cierto que xy se yrefieren a las variables que se asignan a esta estética. Esa es la expectativa también para geom_smooth () y cómo funciona la gramática de los gráficos. Podría haber sido más claro usar diferentes nombres dentro del marco de datos, pero los mantuve como en la pregunta original.
Pedro Aphalo
Será posible en la próxima versión de ggpmisc. ¡Gracias por la sugerencia!
Pedro Aphalo
3
Buen punto @elarry! Esto está relacionado con el funcionamiento de la función parse () de R. A través de prueba y error descubrí que aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~"))hace el trabajo.
Pedro Aphalo
1
@HermanToothrot Por lo general, se prefiere R2 para una regresión, por lo que no hay un r.label predefinido en los datos devueltos por stat_poly_eq(). Puede usar stat_fit_glance(), también del paquete 'ggpmisc', que devuelve R2 como un valor numérico. Vea ejemplos en la página de ayuda y reemplácelos stat(r.squared)por sqrt(stat(r.squared)).
Pedro Aphalo
99

Cambié algunas líneas de la fuente de stat_smoothfunciones relacionadas para crear una nueva función que agregue la ecuación de ajuste y el valor R al cuadrado. ¡Esto también funcionará en gráficos de facetas!

library(devtools)
source_gist("524eade46135f6348140")
df = data.frame(x = c(1:100))
df$y = 2 + 5 * df$x + rnorm(100, sd = 40)
df$class = rep(1:2,50)
ggplot(data = df, aes(x = x, y = y, label=y)) +
  stat_smooth_func(geom="text",method="lm",hjust=0,parse=TRUE) +
  geom_smooth(method="lm",se=FALSE) +
  geom_point() + facet_wrap(~class)

ingrese la descripción de la imagen aquí

Usé el código en la respuesta de @ Ramnath para formatear la ecuación. La stat_smooth_funcfunción no es muy robusta, pero no debería ser difícil jugar con ella.

https://gist.github.com/kdauria/524eade46135f6348140 . Intente actualizar ggplot2si obtiene un error.

kdauria
fuente
2
Muchas gracias. Este no solo funciona para facetas, sino incluso para grupos. Me resulta muy útil para las regresiones por partes, por ejemplo stat_smooth_func(mapping=aes(group=cut(x.val,c(-70,-20,0,20,50,130))),geom="text",method="lm",hjust=0,parse=TRUE), en combinación con EvaluateSmooths de stackoverflow.com/questions/19735149/…
Julian
1
@aelwan, cambia estas líneas: gist.github.com/kdauria/… como quieras. Luego, sourceel archivo completo en su secuencia de comandos.
kdauria
1
@kdauria ¿Qué pasa si tengo varias ecuaciones en cada una de facet_wraps y tengo diferentes valores de y_ en cada una de facet_wrap. ¿Alguna sugerencia de cómo arreglar las posiciones de las ecuaciones? Probé varias opciones de hjust, vjust y angle usando este ejemplo dropbox.com/s/9lk9lug2nwgno2l/R2_facet_wrap.docx?dl=0 pero no pude llevar todas las ecuaciones al mismo nivel en cada una de las facet_wrap
brillante
3
@aelwan, la posición de la ecuación está determinada por estas líneas: gist.github.com/kdauria/… . Hice xposy yposargumentos de la función en el Gist. Entonces, si desea que todas las ecuaciones se superpongan, simplemente configure xposy ypos. De lo contrario, xposy yposse calculan a partir de los datos. Si quieres algo más elegante, no debería ser demasiado difícil agregar algo de lógica dentro de la función. Por ejemplo, tal vez podría escribir una función para determinar qué parte del gráfico tiene el espacio más vacío y colocar la función allí.
kdauria
66
Me encontré con un error con source_gist: Error en r_files [[which]]: tipo de subíndice inválido 'cierre'. Vea esta publicación para la solución: stackoverflow.com/questions/38345894/r-source-gist-not-working
Matifou
73

Modifiqué la publicación de Ramnath para a) hacer más genérico para que acepte un modelo lineal como parámetro en lugar del marco de datos yb) muestra los negativos de manera más adecuada.

lm_eqn = function(m) {

  l <- list(a = format(coef(m)[1], digits = 2),
      b = format(abs(coef(m)[2]), digits = 2),
      r2 = format(summary(m)$r.squared, digits = 3));

  if (coef(m)[2] >= 0)  {
    eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,l)
  } else {
    eq <- substitute(italic(y) == a - b %.% italic(x)*","~~italic(r)^2~"="~r2,l)    
  }

  as.character(as.expression(eq));                 
}

El uso cambiaría a:

p1 = p + geom_text(aes(x = 25, y = 300, label = lm_eqn(lm(y ~ x, df))), parse = TRUE)
Jayden
fuente
17
¡Esto se ve genial! Pero estoy trazando geom_points en múltiples facetas, donde el df difiere según la variable de faceta. ¿Cómo puedo hacer eso?
bshor
24
La solución de Jayden funciona bastante bien, pero el tipo de letra se ve muy feo. Recomendaría cambiar el uso a esto: p1 = p + annotate("text", x = 25, y = 300, label = lm_eqn(lm(y ~ x, df)), colour="black", size = 5, parse=TRUE)editar: esto también resuelve cualquier problema que pueda tener con las letras que aparecen en su leyenda.
Jonas Raedle
1
@ Jonas, por alguna razón me estoy poniendo "cannot coerce class "lm" to a data.frame". Esta alternativa funciona: df.labs <- data.frame(x = 25, y = 300, label = lm_eqn(df))y p <- p + geom_text(data = df.labs, aes(x = x, y = y, label = label), parse = TRUE)
PatrickT
1
@PatrickT: ese es el mensaje de error que obtendría si llamara lm_eqn(lm(...))con la solución de Ramnath. Probablemente probaste este después de probarlo, pero olvidaste asegurarte de haber redefinidolm_eqn
Hamy
@PatrickT: ¿podría hacer que su respuesta sea una respuesta separada? ¡Me encantaría votarlo!
JelenaČuklina
11

Realmente amo la solución @Ramnath. Para permitir el uso de personalizar la fórmula de regresión (en lugar de fijarla como y y x como nombres de variables literales) y agregar el valor p también en la impresión (como comentó @Jerry T), aquí está el mod:

lm_eqn <- function(df, y, x){
    formula = as.formula(sprintf('%s ~ %s', y, x))
    m <- lm(formula, data=df);
    # formating the values into a summary string to print out
    # ~ give some space, but equal size and comma need to be quoted
    eq <- substitute(italic(target) == a + b %.% italic(input)*","~~italic(r)^2~"="~r2*","~~p~"="~italic(pvalue), 
         list(target = y,
              input = x,
              a = format(as.vector(coef(m)[1]), digits = 2), 
              b = format(as.vector(coef(m)[2]), digits = 2), 
             r2 = format(summary(m)$r.squared, digits = 3),
             # getting the pvalue is painful
             pvalue = format(summary(m)$coefficients[2,'Pr(>|t|)'], digits=1)
            )
          )
    as.character(as.expression(eq));                 
}

geom_point() +
  ggrepel::geom_text_repel(label=rownames(mtcars)) +
  geom_text(x=3,y=300,label=lm_eqn(mtcars, 'hp','wt'),color='red',parse=T) +
  geom_smooth(method='lm')

ingrese la descripción de la imagen aquí Desafortunadamente, esto no funciona con facet_wrap o facet_grid.

XX
fuente
Muy ordenado, he hecho referencia aquí . Una aclaración: ¿falta su código ggplot(mtcars, aes(x = wt, y = mpg, group=cyl))+antes de geom_point ()? Una cuestión relacionada con la semi - si nos referimos a HP y peso en el aes()de ggplot, entonces podemos agarrar a utilizar en la llamada a lm_eqn, por lo que entonces sólo tiene que codificar en un solo lugar? Sé que podríamos configurar xvar = "hp"antes de la llamada a ggplot (), y usar xvar en ambas ubicaciones para reemplazar hp , pero parece que esto debería ser innecesario.
Mark Neal
9

Usando ggpubr :

library(ggpubr)

# reproducible data
set.seed(1)
df <- data.frame(x = c(1:100))
df$y <- 2 + 3 * df$x + rnorm(100, sd = 40)

# By default showing Pearson R
ggscatter(df, x = "x", y = "y", add = "reg.line") +
  stat_cor(label.y = 300) +
  stat_regline_equation(label.y = 280)

ingrese la descripción de la imagen aquí

# Use R2 instead of R
ggscatter(df, x = "x", y = "y", add = "reg.line") +
  stat_cor(label.y = 300, 
           aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~"))) +
  stat_regline_equation(label.y = 280)

## compare R2 with accepted answer
# m <- lm(y ~ x, df)
# round(summary(m)$r.squared, 2)
# [1] 0.85

ingrese la descripción de la imagen aquí

zx8754
fuente
¿Has visto una forma programática ordenada para especificar un número label.y?
Mark Neal
@MarkNeal tal vez obtenga el máximo de y luego multiplique por 0.8. label.y = max(df$y) * 0.8
zx8754
1
@MarkNeal buenos puntos, tal vez envíe el problema como solicitud de función en GitHub ggpubr.
zx8754
1
Problema sobre la ubicación automática enviado aquí
Mark Neal
1
@ zx8754, en su diagrama se muestra rho y no R², ¿alguna forma fácil de mostrar R²?
Matmar
6

Aquí está el código más simple para todos

Nota: Mostrar Rho de Pearson y no R ^ 2.

library(ggplot2)
library(ggpubr)

df <- data.frame(x = c(1:100)
df$y <- 2 + 3 * df$x + rnorm(100, sd = 40)
p <- ggplot(data = df, aes(x = x, y = y)) +
        geom_smooth(method = "lm", se=FALSE, color="black", formula = y ~ x) +
        geom_point()+
        stat_cor(label.y = 35)+ #this means at 35th unit in the y axis, the r squared and p value will be shown
        stat_regline_equation(label.y = 30) #this means at 30th unit regresion line equation will be shown

p

Uno de esos ejemplos con mi propio conjunto de datos

Sork-kal
fuente
¡El mismo problema que el anterior, en su diagrama se muestra rho y no R²!
Matmar
3

Inspirado por el estilo de ecuación proporcionado en esta respuesta , un enfoque más genérico (más de un predictor + salida de látex como opción) puede ser:

print_equation= function(model, latex= FALSE, ...){
    dots <- list(...)
    cc= model$coefficients
    var_sign= as.character(sign(cc[-1]))%>%gsub("1","",.)%>%gsub("-"," - ",.)
    var_sign[var_sign==""]= ' + '

    f_args_abs= f_args= dots
    f_args$x= cc
    f_args_abs$x= abs(cc)
    cc_= do.call(format, args= f_args)
    cc_abs= do.call(format, args= f_args_abs)
    pred_vars=
        cc_abs%>%
        paste(., x_vars, sep= star)%>%
        paste(var_sign,.)%>%paste(., collapse= "")

    if(latex){
        star= " \\cdot "
        y_var= strsplit(as.character(model$call$formula), "~")[[2]]%>%
            paste0("\\hat{",.,"_{i}}")
        x_vars= names(cc_)[-1]%>%paste0(.,"_{i}")
    }else{
        star= " * "
        y_var= strsplit(as.character(model$call$formula), "~")[[2]]        
        x_vars= names(cc_)[-1]
    }

    equ= paste(y_var,"=",cc_[1],pred_vars)
    if(latex){
        equ= paste0(equ," + \\hat{\\varepsilon_{i}} \\quad where \\quad \\varepsilon \\sim \\mathcal{N}(0,",
                    summary(MetamodelKdifEryth)$sigma,")")%>%paste0("$",.,"$")
    }
    cat(equ)
}

El modelargumento espera un lmobjeto, el latexargumento es un booleano para pedir un carácter simple o una ecuación con formato de látex, y el ...argumento pasa sus valores aformat función.

También agregué una opción para generarlo como latex para que pueda usar esta función en un rmarkdown como este:


```{r echo=FALSE, results='asis'}
print_equation(model = lm_mod, latex = TRUE)
```

Ahora usándolo:

df <- data.frame(x = c(1:100))
df$y <- 2 + 3 * df$x + rnorm(100, sd = 40)
df$z <- 8 + 3 * df$x + rnorm(100, sd = 40)
lm_mod= lm(y~x+z, data = df)

print_equation(model = lm_mod, latex = FALSE)

Este código produce: y = 11.3382963933174 + 2.5893419 * x + 0.1002227 * z

Y si pedimos una ecuación de látex, redondeando los parámetros a 3 dígitos:

print_equation(model = lm_mod, latex = TRUE, digits= 3)

Esto produce: ecuación de látex

rvezy
fuente
0

Tengo una duda, ¿cómo poner en una ecuación estadísticas significativas de t.test para bheta, usando ggpmisc::stat_poly_eq() ?

ex: expression(hat(Y)== 0000*"**"+0000*"x"*"*"-0000*"x"^2*"**"~~~~"R"^2*":"~~0.000)

Jean Karlos
fuente