Dodaj równanie linii regresji i R ^ 2 na wykresie


227

Zastanawiam się, jak dodać równanie linii regresji i R ^ 2 na ggplot. Mój kod to:

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

Jakakolwiek pomoc będzie doceniona.


1
Aby zapoznać się z grafiką sieci , patrz latticeExtra::lmlineq().
Josh O'Brien

Odpowiedzi:


234

Oto jedno rozwiązanie

# 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)

EDYTOWAĆ. Zrozumiałem źródło, z którego wybrałem ten kod. Oto link do oryginalnego postu w grupach google ggplot2

Wynik


1
@ Komentarz JonasRaedle na temat uzyskiwania lepiej wyglądających tekstów annotatebył poprawny na moim komputerze.
IRTFM

2
Nie wygląda to jak publikacja na moim komputerze, gdzie etykieta jest nadpisywana tyle razy, ile wywoływane są dane, co skutkuje grubym i rozmytym tekstem etykiety. Najpierw działa przekazywanie etykiet do ramki data.frame (zobacz moją sugestię w komentarzu poniżej.
PatrickT

@PatrickT: usuń aes(i odpowiedni ). aessłuży do mapowania zmiennych ramki danych na zmienne wizualne - tutaj nie jest to potrzebne, ponieważ jest tylko jedna instancja, więc możesz umieścić to wszystko w głównym geom_textwywołaniu. Zmienię to w odpowiedzi.
naught101

Problem z tym rozwiązaniem wydaje się polegać na tym, że jeśli zbiór danych jest większy (moje to 370000 obserwacji), funkcja wydaje się nie działać. Poleciłbym rozwiązanie z @kdauria, które robi to samo, ale znacznie szybciej.
Benjamin

3
dla tych, którzy chcą wartości r i p zamiast R2 i równania: eq <- substitute (italic (r) ~ "=" ~ rvalue * "," ~ italic (p) ~ "=" ~ pvalue, list (rvalue = sprintf ("% .2f", znak (coef (m) [2]) * sqrt (podsumowanie (m) $ r.squared)), pvalue = format (podsumowanie (m) $ współczynniki [2,4], cyfry = 2 )))
Jerry T

135

stat_poly_eq()W swoim pakiecie umieściłem statystyki , ggpmiscktóre pozwalają na tę odpowiedź:

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

wprowadź opis zdjęcia tutaj

Ta statystyka działa z dowolnym wielomianem bez brakujących terminów i, mam nadzieję, ma wystarczającą elastyczność, aby być ogólnie użytecznym. Etykiety R ^ 2 lub dostosowane R ^ 2 mogą być używane z dowolną formułą modelu wyposażoną w lm (). Będąc statystyką ggplot, zachowuje się zgodnie z oczekiwaniami zarówno w przypadku grup, jak i aspektów.

Pakiet „ggpmisc” jest dostępny za pośrednictwem CRAN.

Wersja 0.2.6 została właśnie zaakceptowana do CRAN.

Adresuje komentarze @shabbychef i @ MYaseen208.

@ MYaseen208 pokazuje, jak dodać kapelusz .

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

wprowadź opis zdjęcia tutaj

@shabbychef Teraz można dopasować zmienne w równaniu do zmiennych używanych dla etykiet osi. Zastąpić X powiedzmy z Z i Y, z h można by wykorzystać:

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

wprowadź opis zdjęcia tutaj

Będąc tymi normalnymi wyrażeniami R, pisma greckie mogą być teraz używane zarówno w lewej, jak i prawej stronie równania.

[2017-03-08] @elarry Edytuj, aby dokładniej odpowiedzieć na oryginalne pytanie, pokazując, jak dodać przecinek między etykietami równania i 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

wprowadź opis zdjęcia tutaj

[2019-10-20] @ helen.h Podam poniżej przykłady użycia stat_poly_eq()z grupowaniem.

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

wprowadź opis zdjęcia tutaj

wprowadź opis zdjęcia tutaj

[2020-01-21] @Herman Na pierwszy rzut oka może to być nieco sprzeczne z intuicją, ale aby uzyskać jedno równanie podczas korzystania z grupowania, należy postępować zgodnie z gramatyką grafiki. Ogranicz mapowanie, które tworzy grupowanie do poszczególnych warstw (pokazane poniżej), lub zachowaj domyślne mapowanie i zastąp je stałą wartością w warstwie, w której nie chcesz grupowania (np colour = "black".).

Kontynuując z poprzedniego przykładu.

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

wprowadź opis zdjęcia tutaj

[2020-01-22] Dla kompletności przykład z aspektami, pokazujący, że również w tym przypadku oczekiwania gramatyki grafiki są spełnione.

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

wprowadź opis zdjęcia tutaj


1
Należy zauważyć, że xi ywe wzorze odnoszą się do xi ydanych w warstwach wykresu, a niekoniecznie do tych, które w tym czasie my.formulasą konstruowane. Zatem formuła powinna zawsze używać zmiennych xiy?
shabbychef

To jest bardzo prawdziwe, że xi yodnoszą się do dowolnych zmiennych są mapowane do tych estetyki. Jest to oczekiwanie również dla geom_smooth () i tego, jak działa gramatyka grafiki. Używanie różnych nazw w ramce danych mogło być łatwiejsze, ale zachowałem je tak, jak w pierwotnym pytaniu.
Pedro Aphalo,

Będzie to możliwe w następnej wersji ggpmisc. Dzieki za sugestie!
Pedro Aphalo,

3
Dobra uwaga @elarry! Jest to związane ze sposobem działania funkcji parse () R. Dzięki próbom i błędom odkryłem, że to aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~"))działa.
Pedro Aphalo

1
@HermanToothrot Zwykle R2 jest preferowane do regresji, więc nie ma predefiniowanej wartości referencyjnej w danych zwracanych przez stat_poly_eq(). Możesz użyć stat_fit_glance()również z pakietu „ggpmisc”, który zwraca R2 jako wartość liczbową. Zobacz przykłady na stronie pomocy i zamień stat(r.squared)na sqrt(stat(r.squared)).
Pedro Aphalo

99

Zmieniłem kilka wierszy źródła stat_smoothi powiązanych funkcji, aby utworzyć nową funkcję, która dodaje równanie dopasowania i wartość R do kwadratu. Działa to również na wykresach aspektów!

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)

wprowadź opis zdjęcia tutaj

Użyłem kodu w odpowiedzi @ Ramnatha, aby sformatować równanie. Ta stat_smooth_funcfunkcja nie jest bardzo niezawodna, ale nie powinno być trudno się nią bawić.

https://gist.github.com/kdauria/524eade46135f6348140 . Spróbuj zaktualizować, ggplot2jeśli pojawi się błąd.


2
Wielkie dzięki. Ten działa nie tylko dla aspektów, ale nawet dla grup. Uważam, że jest to bardzo przydatne w przypadku regresji częściowej, np. stat_smooth_func(mapping=aes(group=cut(x.val,c(-70,-20,0,20,50,130))),geom="text",method="lm",hjust=0,parse=TRUE)W połączeniu z EvaluateSmooths z stackoverflow.com/questions/19735149/…
Julian

1
@aelwan, zmień te linie: gist.github.com/kdauria/... jak chcesz. Następnie sourcecały plik w skrypcie.
kdauria

1
@kdauria Co jeśli mam kilka równań w każdym z aspektów fac_wraps i mam różne wartości y_ w każdym z aspektów fac_wrap. Jakieś sugestie, jak naprawić pozycje równań? Wypróbowałem kilka opcji hjust, vjust i angle, korzystając z tego przykładu dropbox.com/s/9lk9lug2nwgno2l/R2_facet_wrap.docx?dl=0, ale nie mogłem doprowadzić wszystkich równań do tego samego poziomu w każdym z aspektów fac_wrap
błyszczący

3
@aelwan, pozycja równania jest określona przez następujące linie: gist.github.com/kdauria/… . Zrobiłem xposi yposargumenty funkcji w Gist. Więc jeśli chcesz, aby wszystkie równania zachodziły na siebie, po prostu ustaw xposi ypos. W przeciwnym razie xposi ypossą obliczane na podstawie danych. Jeśli chcesz czegoś bardziej wymyślnego, dodanie logiki wewnątrz funkcji nie powinno być zbyt trudne. Na przykład, być może mógłbyś napisać funkcję określającą, która część wykresu ma najwięcej pustej przestrzeni i umieścić tam funkcję.
kdauria

6
Wystąpił błąd związany z source_gist: Error in r_files [[which]]: niepoprawny typ indeksu „zamknięcie”. Zobacz ten post dotyczący rozwiązania: stackoverflow.com/questions/38345894/r-source-gist-not-working
Matifou

73

Zmodyfikowałem post Ramnatha na a) uczynić bardziej ogólnym, więc akceptuje model liniowy jako parametr, a nie ramkę danych, i b) wyświetla odpowiednio negatywy.

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));                 
}

Użycie zmieni się na:

p1 = p + geom_text(aes(x = 25, y = 300, label = lm_eqn(lm(y ~ x, df))), parse = TRUE)

17
To wygląda świetnie! Ale rysuję punkty geom na wielu aspektach, gdzie df różni się w zależności od zmiennej facet. W jaki sposób mogę to zrobić?
bshor,

24
Rozwiązanie Jaydena działa całkiem dobrze, ale krój pisma wygląda bardzo brzydko. Polecam zmianę sposobu użycia na: p1 = p + annotate("text", x = 25, y = 300, label = lm_eqn(lm(y ~ x, df)), colour="black", size = 5, parse=TRUE)edytuj: rozwiązuje to również wszelkie problemy z literami pojawiającymi się w legendzie.
Jonas Raedle,

1
@ Jonas, z jakiegoś powodu dostaję "cannot coerce class "lm" to a data.frame". Ta alternatywa działa: df.labs <- data.frame(x = 25, y = 300, label = lm_eqn(df))i p <- p + geom_text(data = df.labs, aes(x = x, y = y, label = label), parse = TRUE)
PatrickT

1
@PatrickT - taki komunikat o błędzie pojawia się, gdy zadzwonisz lm_eqn(lm(...))z rozwiązaniem Ramnatha. Prawdopodobnie próbowałeś tego po wypróbowaniu, ale zapomniałeś upewnić się, że przedefiniowałeślm_eqn
Hamy

@PatrickT: czy możesz uczynić swoją odpowiedź osobną odpowiedzią? Z przyjemnością zagłosuję!
JelenaČuklina

11

naprawdę kocham rozwiązanie @Ramnath. Aby umożliwić użycie w celu dostosowania formuły regresji (zamiast ustalonej jako y i x jako dosłowne nazwy zmiennych) i dodania wartości p również do wydruku (jak skomentował @Jerry T), oto 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')

wprowadź opis zdjęcia tutaj Niestety nie działa to z facet_wrap ani facet_grid.


Bardzo schludnie, wspomniałem tutaj . Wyjaśnienie - czy brakuje twojego kodu ggplot(mtcars, aes(x = wt, y = mpg, group=cyl))+przed geom_point ()? Semi-pokrewne pytanie - jeśli odwołujemy się do KM i wag w aes()dla ggplot, możemy chwycić je używać w wywołaniu lm_eqn, tak to mamy tylko do kodu w jednym miejscu? Wiem, że moglibyśmy skonfigurować xvar = "hp"przed wywołaniem ggplot () i użyć xvara w obu lokalizacjach, aby zastąpić hp , ale wydaje się, że to powinno być niepotrzebne.
Mark Neal

9

Za pomocą 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)

wprowadź opis zdjęcia tutaj

# 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

wprowadź opis zdjęcia tutaj


Czy widziałeś schludny programowy sposób na określenie liczby label.y?
Mark Neal

@MarkNeal może uzyskać maks. Y, a następnie pomnożyć przez 0,8. label.y = max(df$y) * 0.8
zx8754

1
@MarkNeal dobre punkty, może zgłosić problem jako żądanie funkcji w GitHub ggpubr.
zx8754

1
Kwestia lokalizacji auto składać tutaj
Mark Neal

1
@ zx8754, na twojej działce pokazane jest rho a nie R², jakikolwiek łatwy sposób pokazać R²?
matmar

5

Oto najprostszy kod dla każdego

Uwaga: Pokazuje Rho Pearsona, a nie 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

Jeden taki przykład z moim własnym zestawem danych


Ten sam problem jak powyżej, na twojej działce pokazane jest rho, a nie R²!
matmar

3

Zainspirowany stylem równania podanym w tej odpowiedzi , bardziej ogólne podejście (więcej niż jeden predyktor + wyjście lateksu jako opcja) może być:

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)
}

modelArgumentem oczekuje lmcel, latexargument jest logiczna, aby poprosić o charakterze prostego równania lub lateksu sformatowany, a ...argumentem przekazać swoje wartości do formatfunkcji.

Dodałem również opcję wyświetlania go jako lateksu, abyś mógł użyć tej funkcji w rmarkdown jak poniżej:


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

Teraz go używam:

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)

Ten kod daje: y = 11.3382963933174 + 2.5893419 * x + 0.1002227 * z

A jeśli poprosimy o równanie lateksowe, zaokrąglając parametry do 3 cyfr:

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

Daje to: równanie lateksowe


0

Mam wątpliwości, jak umieścić znamienne statystyki t.test dla bhety w równaniu, używając ggpmisc::stat_poly_eq()?

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

Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.