Tambahkan persamaan garis regresi dan R ^ 2 pada grafik

228

Saya ingin tahu bagaimana cara menambahkan persamaan garis regresi dan R ^ 2 pada ggplot. Kode saya adalah:

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

Bantuan apa pun akan sangat dihargai.

MYaseen208
sumber
1
Untuk grafik kisi , lihat latticeExtra::lmlineq().
Josh O'Brien

Jawaban:

234

Inilah salah satu solusinya

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

EDIT. Saya menemukan sumber dari mana saya mengambil kode ini. Berikut ini tautan ke pos asli di grup google ggplot2

Keluaran

Ramnath
sumber
1
Komentar @ JonasRaedle tentang mendapatkan teks yang terlihat lebih baik annotatesudah benar di mesin saya.
IRTFM
2
Ini tidak terlihat seperti output yang diposting di mesin saya, di mana label ditimpa sebanyak data yang dipanggil, menghasilkan teks label yang tebal dan buram. Melewati label ke data.frame berfungsi pertama (lihat saran saya dalam komentar di bawah.
PatrickT
@ Patrickrick: menghapus aes(dan yang sesuai ). aesadalah untuk memetakan variabel kerangka data ke variabel visual - itu tidak diperlukan di sini, karena hanya ada satu contoh, sehingga Anda dapat menempatkan semuanya dalam geom_textpanggilan utama . Saya akan mengedit ini untuk jawabannya.
naught101
Masalah dengan solusi ini tampaknya, bahwa jika dataset lebih besar (milik saya adalah 370000 pengamatan) fungsi tersebut tampaknya gagal. Saya akan merekomendasikan solusi dari @ kdauria yang melakukan hal yang sama, tetapi jauh lebih cepat.
Benjamin
3
bagi mereka yang menginginkan nilai r dan p daripada R2 dan persamaan: eq <- subtitusi (italic (r) ~ "=" ~ rvalue * "," ~ italic (p) ~ "=" ~ pvalue, daftar (rvalue = sprintf ("% .2f", tanda (coef (m) [2]) * sqrt (ringkasan (m) $ r.squared)), pvalue = format (ringkasan (m) $ koefisien [2,4], digit = 2 )))
Jerry T
135

Saya menyertakan statistik stat_poly_eq()dalam paket saya ggpmiscyang memungkinkan jawaban ini:

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

masukkan deskripsi gambar di sini

Statistik ini berfungsi dengan semua polinomial tanpa istilah yang hilang, dan mudah-mudahan memiliki fleksibilitas yang cukup untuk berguna secara umum. Label R ^ 2 atau R ^ 2 yang disesuaikan dapat digunakan dengan formula model apa pun yang dilengkapi dengan lm (). Menjadi statistik ggplot berperilaku seperti yang diharapkan baik dengan kelompok dan segi.

Paket 'ggpmisc' tersedia melalui CRAN.

Versi 0.2.6 baru saja diterima oleh CRAN.

Ini alamat komentar oleh @shabbychef dan @ MYaseen208.

@ MYaseen208 ini menunjukkan cara menambahkan topi .

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

masukkan deskripsi gambar di sini

@shabbychef Sekarang dimungkinkan untuk mencocokkan variabel dalam persamaan dengan yang digunakan untuk label sumbu. Untuk mengganti x dengan mengatakan z dan y dengan h satu akan menggunakan:

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

masukkan deskripsi gambar di sini

Menjadi ekspresi parsed R normal ini huruf yunani sekarang juga dapat digunakan baik dalam lhs dan rhs dari persamaan.

[2017-03-08] @elarry Edit untuk lebih tepatnya menjawab pertanyaan awal, menunjukkan cara menambahkan koma antara persamaan- dan R2-label.

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

masukkan deskripsi gambar di sini

[2019-10-20] @ helen.h Saya berikan contoh penggunaan stat_poly_eq()pengelompokan di bawah ini.

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

masukkan deskripsi gambar di sini

masukkan deskripsi gambar di sini

[2020-01-21] @Herman Mungkin agak kontra-intuitif pada pandangan pertama, tetapi untuk mendapatkan persamaan tunggal saat menggunakan pengelompokan, seseorang harus mengikuti tata bahasa grafik. Batasi pemetaan yang membuat pengelompokan ke lapisan individual (diperlihatkan di bawah) atau pertahankan pemetaan default dan timpa dengan nilai konstan di layer di mana Anda tidak ingin pengelompokan (misalnya colour = "black").

Melanjutkan dari contoh sebelumnya.

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

masukkan deskripsi gambar di sini

[2020-01-22] Demi kelengkapan contoh dengan segi, menunjukkan bahwa juga dalam hal ini harapan tata bahasa grafis terpenuhi.

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

masukkan deskripsi gambar di sini

Pedro Aphalo
sumber
1
Perlu dicatat bahwa xdan ydalam rumus mengacu pada xdan ydata dalam lapisan plot, dan belum tentu untuk orang-orang di lingkup pada saat my.formuladibangun. Jadi formula harus selalu menggunakan variabel x dan y?
shabbychef
Sangat benar bahwa xdan ymerujuk pada variabel apa pun yang dipetakan ke estetika ini. Itulah harapan juga untuk geom_smooth () dan bagaimana tata bahasa grafis bekerja. Itu bisa lebih jelas untuk menggunakan nama yang berbeda dalam bingkai data tetapi saya hanya menyimpannya seperti pada pertanyaan awal.
Pedro Aphalo
Akan dimungkinkan di versi selanjutnya ggpmisc. Terima kasih untuk sarannya!
Pedro Aphalo
3
Poin bagus dari @elarry! Ini terkait dengan bagaimana fungsi p's parse () bekerja. Melalui trial and error saya menemukan bahwa aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~"))melakukan pekerjaan.
Pedro Aphalo
1
@HermanToothrot Biasanya R2 lebih disukai untuk regresi, sehingga tidak ada r.label yang telah ditentukan dalam data yang dikembalikan oleh stat_poly_eq(). Anda dapat menggunakan stat_fit_glance(), juga dari paket 'ggpmisc', yang mengembalikan R2 sebagai nilai numerik. Lihat contoh di halaman bantuan, dan ganti stat(r.squared)dengan sqrt(stat(r.squared)).
Pedro Aphalo
99

Saya mengubah beberapa baris sumber stat_smoothdan fungsi terkait untuk membuat fungsi baru yang menambahkan persamaan fit dan nilai kuadrat R. Ini juga akan bekerja pada plot facet!

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)

masukkan deskripsi gambar di sini

Saya menggunakan kode dalam jawaban @ Ramnath untuk memformat persamaan. The stat_smooth_funcFungsi ini tidak sangat kuat, tetapi seharusnya tidak sulit untuk bermain-main dengan hal itu.

https://gist.github.com/kdauria/524eade46135f6348140 . Coba perbarui ggplot2jika Anda mendapatkan kesalahan.

kdauria
sumber
2
Terimakasih banyak. Yang ini tidak hanya berfungsi untuk segi, tetapi bahkan untuk kelompok. Saya merasa sangat berguna untuk regresi sedikit demi sedikit, misalnya stat_smooth_func(mapping=aes(group=cut(x.val,c(-70,-20,0,20,50,130))),geom="text",method="lm",hjust=0,parse=TRUE), dalam kombinasi dengan EvaluateSmooths dari stackoverflow.com/questions/19735149/…
Julian
1
@aelwan, mengubah baris ini: gist.github.com/kdauria/... yang Anda inginkan. Kemudian sourceseluruh file dalam skrip Anda.
kdauria
1
@kdauria Bagaimana jika saya memiliki beberapa persamaan di setiap facet_wraps dan saya memiliki nilai y_ yang berbeda di setiap facet_wrap. Ada saran bagaimana cara memperbaiki posisi persamaan? Saya mencoba beberapa opsi hjust, vjust dan angle menggunakan contoh ini dropbox.com/s/9lk9lug2nwgno2l/R2_facet_wrap.docx?dl=0 tapi saya tidak bisa membawa semua persamaan pada level yang sama di setiap facet_wrap
shiny
3
@aelwan, posisi persamaan ditentukan oleh baris-baris ini: gist.github.com/kdauria/… . Saya membuat xposdan yposargumen fungsi di Intisari. Jadi jika Anda ingin semua persamaan tumpang tindih, tetapkan saja xposdan ypos. Kalau tidak, xposdan yposdihitung dari data. Jika Anda menginginkan sesuatu yang lebih menarik, seharusnya tidak terlalu sulit untuk menambahkan beberapa logika di dalam fungsi. Sebagai contoh, mungkin Anda bisa menulis fungsi untuk menentukan bagian mana dari grafik yang memiliki ruang paling kosong dan meletakkan fungsi di sana.
kdauria
6
Saya mengalami kesalahan dengan source_gist: Kesalahan pada r_files [[yang]]: tipe subscript 'closure' tidak valid. Lihat posting ini untuk solusinya: stackoverflow.com/questions/38345894/r-source-gist-not-working
Matifou
73

Saya telah memodifikasi posting Ramnath ke a) membuat lebih umum sehingga menerima model linier sebagai parameter daripada bingkai data dan b) menampilkan negatif lebih tepat.

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

Penggunaan akan berubah menjadi:

p1 = p + geom_text(aes(x = 25, y = 300, label = lm_eqn(lm(y ~ x, df))), parse = TRUE)
Jayden
sumber
17
Ini terlihat hebat! Tapi aku merencanakan geom_points pada banyak sisi, di mana df berbeda berdasarkan pada variabel facet. Bagaimana aku melakukan itu?
bshor
24
Solusi Jayden bekerja dengan sangat baik, tetapi jenis huruf terlihat sangat jelek. Saya akan merekomendasikan mengubah penggunaannya menjadi ini: p1 = p + annotate("text", x = 25, y = 300, label = lm_eqn(lm(y ~ x, df)), colour="black", size = 5, parse=TRUE)sunting: ini juga menyelesaikan masalah yang mungkin Anda miliki dengan surat-surat yang muncul dalam legenda Anda.
Jonas Raedle
1
@ Jonas, entah kenapa aku mengerti "cannot coerce class "lm" to a data.frame". Alternatif ini berfungsi: df.labs <- data.frame(x = 25, y = 300, label = lm_eqn(df))dan p <- p + geom_text(data = df.labs, aes(x = x, y = y, label = label), parse = TRUE)
PatrickT
1
@ Patrickrick - Itulah pesan kesalahan yang akan Anda dapatkan jika Anda menelepon lm_eqn(lm(...))dengan solusi Ramnath. Anda mungkin mencoba yang ini setelah mencoba yang itu tetapi lupa untuk memastikan bahwa Anda telah mendefinisikan ulanglm_eqn
Hamy
@ Patrickrick: bisakah Anda membuat jawaban Anda sebagai jawaban yang terpisah? Saya akan dengan senang hati memilihnya!
JelenaČuklina
11

benar-benar menyukai solusi @Ramnath. Untuk mengizinkan penggunaan untuk mengkustomisasi rumus regresi (alih-alih tetap sebagai y dan x sebagai nama variabel literal), dan menambahkan nilai p ke dalam cetakan juga (seperti yang dikomentari @Jerry T), berikut adalah 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')

masukkan deskripsi gambar di sini Sayangnya, ini tidak berfungsi dengan facet_wrap atau facet_grid.

XX
sumber
Sangat rapi, saya referensi di sini . Klarifikasi - apakah kode Anda hilang ggplot(mtcars, aes(x = wt, y = mpg, group=cyl))+sebelum geom_point ()? Pertanyaan semi-terkait - jika kita merujuk ke hp dan wt di aes()for ggplot, dapatkah kita mengambilnya untuk digunakan dalam panggilan ke lm_eqn, jadi kita hanya perlu kode di satu tempat? Saya tahu kita bisa mengatur xvar = "hp"sebelum panggilan ggplot (), dan menggunakan xvar di kedua lokasi untuk mengganti hp , tetapi ini terasa tidak perlu.
Mark Neal
9

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

masukkan deskripsi gambar di sini

# 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

masukkan deskripsi gambar di sini

zx8754
sumber
Pernahkah Anda melihat cara terprogram yang rapi untuk menentukan angka label.y?
Mark Neal
@MarkNeal mungkin mendapatkan maksimum y lalu kalikan dengan 0.8. label.y = max(df$y) * 0.8
zx8754
1
@MarkNeal poin bagus, mungkin mengirimkan masalah sebagai permintaan fitur di GitHub ggpubr.
zx8754
1
Masalah tentang lokasi otomatis yang diajukan di sini
Mark Neal
1
@ zx8754, di plot Anda ditunjukkan rho dan bukan R², ada cara mudah untuk menunjukkan R²?
matmar
6

Inilah kode yang paling sederhana untuk semua orang

Catatan: Menampilkan Pearson's Rho dan bukan 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

Salah satu contohnya dengan dataset saya sendiri

Sork-kal
sumber
Masalah yang sama seperti di atas, dalam plot Anda ditampilkan rho dan bukan R²!
matmar
3

Terinspirasi oleh gaya persamaan yang disediakan dalam jawaban ini , pendekatan yang lebih umum (lebih dari satu prediktor + keluaran lateks sebagai opsi) dapat berupa:

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

The modelArgumen mengharapkan lmobjek, latexargumen adalah boolean untuk meminta karakter sederhana atau persamaan lateks-formated, dan ...argumen lulus nilai kepada formatfungsi.

Saya juga menambahkan opsi untuk menampilkannya sebagai lateks sehingga Anda dapat menggunakan fungsi ini dalam rmarkdown seperti ini:


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

Sekarang menggunakannya:

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)

Kode ini menghasilkan: y = 11.3382963933174 + 2.5893419 * x + 0.1002227 * z

Dan jika kita meminta persamaan lateks, membulatkan parameter menjadi 3 digit:

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

Ini menghasilkan: persamaan lateks

rvezy
sumber
0

Saya ragu, bagaimana menggunakan statistik t.test yang signifikan untuk bheta dalam persamaan, menggunakan ggpmisc::stat_poly_eq()?

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

Jean Karlos
sumber