Wie zeichne ich ein Trichterdiagramm mit ggplot2 in R?

12

Als Titel muss ich so etwas zeichnen:

Alt-Text

Kann ggplot oder andere Pakete, wenn ggplot nicht in der Lage ist, verwendet werden, um so etwas zu zeichnen?

lokheart
quelle
2
Ich habe ein paar Ideen, wie dies zu tun und zu implementieren ist, würde mich aber über einige Daten freuen, mit denen ich spielen kann. Irgendwelche Ideen dazu?
Chase
1
Ja, ggplot kann problemlos ein Diagramm zeichnen, das aus Punkten und Linien besteht. Mit geom_smooth erhalten Sie 95% des Weges. Wenn Sie weitere Ratschläge wünschen, müssen Sie weitere Details angeben.
Hadley
2
Dies ist kein Trichterplot. Stattdessen werden die Linien offensichtlich aus Schätzungen von Standardfehlern basierend auf der Anzahl der Zulassungen konstruiert. Sie scheinen beabsichtigt zu sein, einen bestimmten Anteil von Daten einzuschließen , wodurch sie Toleranzgrenzen erreichen würden. Sie haben wahrscheinlich die Form y = Grundlinie + Konstante / Sqrt (# Zulassungen * f (Grundlinie)). Sie könnten den Code in den vorhandenen Antworten ändern, um die Linien grafisch darzustellen, aber Sie müssten wahrscheinlich Ihre eigene Formel angeben, um sie zu berechnen: Die Beispiele, die ich gesehen habe, zeichnen Konfidenzintervalle für die angepasste Linie selbst . Deshalb sehen sie so anders aus.
whuber
@whuber (+1) Das ist in der Tat ein sehr guter Punkt. Ich hoffe, dass dies ohnehin einen guten Ausgangspunkt bietet (auch wenn mein R-Code nicht so optimiert ist).
Chl
Ggplot bietet weiterhin die Möglichkeit stat_quantile(), bedingte Quantile auf ein Streudiagramm zu setzen. Mit dem Formelparameter können Sie dann die Funktionsform der Quantilregression steuern. Ich würde Dinge wie Formel = vorschlagen y~ns(x,4), um eine glatte Passform zu erhalten.
Shea Parkes

Antworten:

12

Obwohl es Raum für Verbesserungen gibt, ist hier ein kleiner Versuch mit simulierten (heteroskedastischen) Daten:

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

Alt-Text

chl
quelle
20

Wenn Sie nach dieser Art von Trichterdiagramm (Metaanalyse) suchen , ist möglicherweise Folgendes der Ausgangspunkt:

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

Alt-Text

Bernd Weiss
quelle
1
Das Vorhandensein des linetype=2Arguments in den aes()Klammern - das Zeichnen der 99% -Zeilen - führt zu einem Fehler "Kontinuierliche Variable kann nicht dem Linientyp zugeordnet werden" mit dem aktuellen ggplot2 (0.9.3.1). Änderung geom_line(aes(x = number.seq, y = number.ll999, linetype = 2), data = dfCI)zu geom_line(aes(x = number.seq, y = number.ll999), linetype = 2, data = dfCI)Werken für mich. Fühlen Sie sich frei, die ursprüngliche Antwort zu ändern und diese zu verlieren.
2

Der Code von Bernd Weiss ist sehr hilfreich. Ich habe unten einige Änderungen vorgenommen, um einige Funktionen zu ändern / hinzuzufügen:

  1. Verwendete den Standardfehler als Maß für die Genauigkeit, was typischer für die Trichterdiagramme ist, die ich sehe (in der Psychologie).
  2. Vertauschte die Achsen, sodass die Genauigkeit (Standardfehler) auf der y-Achse und die Effektgröße auf der x-Achse liegt
  3. Wird geom_segmentanstelle geom_lineder Linie verwendet, die das metaanalytische Mittel abgrenzt, sodass es dieselbe Höhe wie die Linien hat, die die 95% - und 99% -Konfidenzbereiche abgrenzen
  4. Anstatt den metaanalytischen Mittelwert aufzuzeichnen, habe ich das 95% -Konfidenzintervall aufgezeichnet

Mein Code verwendet als Beispiel einen metaanalytischen Mittelwert von 0,0892 (se = 0,0035), aber Sie können Ihre eigenen Werte ersetzen.

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

Geben Sie hier die Bildbeschreibung ein

jsakaluk
quelle