Cara menggambar plot corong menggunakan ggplot2 di R?

12

Sebagai judul, saya perlu menggambar sesuatu seperti ini:

teks alternatif

Bisakah ggplot, atau paket lain jika ggplot tidak mampu, digunakan untuk menggambar sesuatu seperti ini?

lokheart
sumber
2
Saya punya beberapa ide tentang bagaimana melakukan dan mengimplementasikan ini, tetapi akan sangat menghargai memiliki beberapa data untuk dimainkan. Ada ide tentang itu?
Chase
1
Ya, ggplot dapat dengan mudah menggambar plot yang terdiri dari titik dan garis;) geom_smooth akan memberi Anda 95% jalan - jika Anda menginginkan lebih banyak saran, Anda perlu memberikan rincian lebih lanjut.
Hadley
2
Ini bukan plot corong. Alih-alih, garis-garis jelas dibangun dari perkiraan kesalahan standar berdasarkan jumlah penerimaan. Mereka tampaknya dimaksudkan untuk melampirkan proporsi data tertentu, yang akan membuatnya menjadi batas toleransi. Mereka kemungkinan dari bentuk y = baseline + konstanta / Sqrt (# penerimaan * f (baseline)). Anda dapat memodifikasi kode dalam respons yang ada untuk membuat grafik garis, tetapi Anda mungkin perlu menyediakan rumus Anda sendiri untuk menghitungnya: contoh saya telah melihat interval kepercayaan plot untuk garis yang dipasang itu sendiri . Itu sebabnya mereka terlihat sangat berbeda.
whuber
@whuber (+1) Memang itu poin yang sangat bagus. Saya harap ini mungkin memberikan titik awal yang baik (bahkan jika kode R saya tidak dioptimalkan).
chl
Ggplot masih menyediakan stat_quantile()untuk menempatkan kuantil bersyarat pada sebar sebaran. Anda kemudian dapat mengontrol bentuk fungsional regresi kuantil dengan parameter rumus. Saya menyarankan hal-hal seperti formula = y~ns(x,4)untuk mendapatkan kecocokan yang halus.
Shea Parkes

Jawaban:

12

Meskipun ada ruang untuk perbaikan, berikut adalah upaya kecil dengan data simulasi (heteroscedastic):

library(ggplot2)
set.seed(101)
x <- runif(100, min=1, max=10)
y <- rnorm(length(x), mean=5, sd=0.1*x)
df <- data.frame(x=x*70, y=y)
m <- lm(y ~ x, data=df) 
fit95 <- predict(m, interval="conf", level=.95)
fit99 <- predict(m, interval="conf", level=.999)
df <- cbind.data.frame(df, 
                       lwr95=fit95[,"lwr"],  upr95=fit95[,"upr"],     
                       lwr99=fit99[,"lwr"],  upr99=fit99[,"upr"])

p <- ggplot(df, aes(x, y)) 
p + geom_point() + 
    geom_smooth(method="lm", colour="black", lwd=1.1, se=FALSE) + 
    geom_line(aes(y = upr95), color="black", linetype=2) + 
    geom_line(aes(y = lwr95), color="black", linetype=2) +
    geom_line(aes(y = upr99), color="red", linetype=3) + 
    geom_line(aes(y = lwr99), color="red", linetype=3)  + 
    annotate("text", 100, 6.5, label="95% limit", colour="black", 
             size=3, hjust=0) +
    annotate("text", 100, 6.4, label="99.9% limit", colour="red", 
             size=3, hjust=0) +
    labs(x="No. admissions...", y="Percentage of patients...") +    
    theme_bw() 

teks alternatif

chl
sumber
20

Jika Anda mencari jenis alur corong (meta-analisis) ini , maka yang berikut ini mungkin menjadi titik awal:

library(ggplot2)

set.seed(1)
p <- runif(100)
number <- sample(1:1000, 100, replace = TRUE)
p.se <- sqrt((p*(1-p)) / (number))
df <- data.frame(p, number, p.se)

## common effect (fixed effect model)
p.fem <- weighted.mean(p, 1/p.se^2)

## lower and upper limits for 95% and 99.9% CI, based on FEM estimator
number.seq <- seq(0.001, max(number), 0.1)
number.ll95 <- p.fem - 1.96 * sqrt((p.fem*(1-p.fem)) / (number.seq)) 
number.ul95 <- p.fem + 1.96 * sqrt((p.fem*(1-p.fem)) / (number.seq)) 
number.ll999 <- p.fem - 3.29 * sqrt((p.fem*(1-p.fem)) / (number.seq)) 
number.ul999 <- p.fem + 3.29 * sqrt((p.fem*(1-p.fem)) / (number.seq)) 
dfCI <- data.frame(number.ll95, number.ul95, number.ll999, number.ul999, number.seq, p.fem)

## draw plot
fp <- ggplot(aes(x = number, y = p), data = df) +
    geom_point(shape = 1) +
    geom_line(aes(x = number.seq, y = number.ll95), data = dfCI) +
    geom_line(aes(x = number.seq, y = number.ul95), data = dfCI) +
    geom_line(aes(x = number.seq, y = number.ll999), linetype = "dashed", data = dfCI) +
    geom_line(aes(x = number.seq, y = number.ul999), linetype = "dashed", data = dfCI) +
    geom_hline(aes(yintercept = p.fem), data = dfCI) +
    scale_y_continuous(limits = c(0,1.1)) +
  xlab("number") + ylab("p") + theme_bw() 
fp

teks alternatif

Bernd Weiss
sumber
1
Kehadiran linetype=2argumen di dalam aes()tanda kurung - memplot garis 99% - memunculkan kesalahan "variabel kontinu tidak dapat dipetakan ke linetype" dengan ggplot2 saat ini (0.9.3.1). Mengubah geom_line(aes(x = number.seq, y = number.ll999, linetype = 2), data = dfCI)ke geom_line(aes(x = number.seq, y = number.ll999), linetype = 2, data = dfCI)bekerja untuk saya. Jangan ragu untuk mengubah jawaban asli dan kehilangan ini.
2

Kode Bernd Weiss sangat membantu. Saya membuat beberapa amandemen di bawah ini, untuk mengubah / menambahkan beberapa fitur:

  1. Digunakan kesalahan standar sebagai ukuran presisi, yang lebih khas dari plot corong yang saya lihat (dalam psikologi)
  2. Tukar sumbu, sehingga presisi (kesalahan standar) ada pada sumbu y, dan ukuran efek ada pada sumbu x
  3. Digunakan geom_segmentsebagai ganti geom_lineuntuk garis yang membatasi rata-rata meta-analitik, sehingga akan sama tingginya dengan garis yang membatasi wilayah kepercayaan 95% dan 99%
  4. Alih-alih memplot rata-rata meta-analitik, saya merencanakan interval kepercayaan 95%

Kode saya menggunakan rata-rata meta-analitik 0,0892 (se = 0,0035) sebagai contoh, tetapi Anda dapat mengganti nilai Anda sendiri.

estimate = 0.0892
se = 0.0035

#Store a vector of values that spans the range from 0
#to the max value of impression (standard error) in your dataset.
#Make the increment (the final value) small enough (I choose 0.001)
#to ensure your whole range of data is captured
se.seq=seq(0, max(dat$corr_zi_se), 0.001)

#Compute vectors of the lower-limit and upper limit values for
#the 95% CI region
ll95 = estimate-(1.96*se.seq)
ul95 = estimate+(1.96*se.seq)

#Do this for a 99% CI region too
ll99 = estimate-(3.29*se.seq)
ul99 = estimate+(3.29*se.seq)

#And finally, calculate the confidence interval for your meta-analytic estimate 
meanll95 = estimate-(1.96*se)
meanul95 = estimate+(1.96*se)

#Put all calculated values into one data frame
#You might get a warning about '...row names were found from a short variable...' 
#You can ignore it.
dfCI = data.frame(ll95, ul95, ll99, ul99, se.seq, estimate, meanll95, meanul95)


#Draw Plot
fp = ggplot(aes(x = se, y = Zr), data = dat) +
  geom_point(shape = 1) +
  xlab('Standard Error') + ylab('Zr')+
  geom_line(aes(x = se.seq, y = ll95), linetype = 'dotted', data = dfCI) +
  geom_line(aes(x = se.seq, y = ul95), linetype = 'dotted', data = dfCI) +
  geom_line(aes(x = se.seq, y = ll99), linetype = 'dashed', data = dfCI) +
  geom_line(aes(x = se.seq, y = ul99), linetype = 'dashed', data = dfCI) +
  geom_segment(aes(x = min(se.seq), y = meanll95, xend = max(se.seq), yend = meanll95), linetype='dotted', data=dfCI) +
  geom_segment(aes(x = min(se.seq), y = meanul95, xend = max(se.seq), yend = meanul95), linetype='dotted', data=dfCI) +
  scale_x_reverse()+
  scale_y_continuous(breaks=seq(-1.25,2,0.25))+
  coord_flip()+
  theme_bw()
fp

masukkan deskripsi gambar di sini

jsakaluk
sumber