Comment dessiner un entonnoir en utilisant ggplot2 dans R?

12

En tant que titre, je dois dessiner quelque chose comme ceci:

texte alternatif

Est-ce que ggplot, ou d'autres packages si ggplot n'est pas capable, peut être utilisé pour dessiner quelque chose comme ça?

lokheart
la source
2
J'ai quelques idées sur la façon de le faire et de l'implémenter, mais j'aimerais avoir quelques données avec lesquelles jouer. Des idées à ce sujet?
Chase
1
Oui, ggplot peut facilement dessiner un tracé composé de points et de lignes;) geom_smooth vous donnera 95% du chemin - si vous voulez plus de conseils, vous devrez fournir plus de détails.
hadley
2
Ce n'est pas un complot en entonnoir. Au lieu de cela, les lignes sont évidemment construites à partir d'estimations d'erreurs standard basées sur le nombre d'admissions. Ils semblent destinés à contenir une proportion spécifiée de données, ce qui en ferait des limites de tolérance. Ils sont probablement de la forme y = ligne de base + constante / Sqrt (# admissions * f (ligne de base)). Vous pouvez modifier le code dans les réponses existantes pour représenter graphiquement les lignes, mais vous devrez probablement fournir votre propre formule pour les calculer: les exemples que j'ai vus tracer des intervalles de confiance pour la ligne ajustée elle-même . Voilà pourquoi ils sont si différents.
whuber
@whuber (+1) C'est vraiment un très bon point. J'espère que cela pourrait fournir un bon point de départ de toute façon (même si mon code R n'est pas optimisé).
chl
Ggplot prévoit toujours stat_quantile()de mettre des quantiles conditionnels sur un nuage de points. Vous pouvez ensuite contrôler la forme fonctionnelle de la régression quantile avec le paramètre de formule. Je suggérerais des choses comme la formule = y~ns(x,4)pour obtenir un ajustement cannelé lisse.
Shea Parkes

Réponses:

12

Bien qu'il y ait place à amélioration, voici une petite tentative avec des données simulées (hétéroscédastiques):

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

texte alternatif

chl
la source
20

Si vous recherchez ce type d' entonnoir (méta-analyse) , alors ce qui suit peut être un point de départ:

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

texte alternatif

Bernd Weiss
la source
1
La présence de l' linetype=2argument à l'intérieur des aes()parenthèses - traçant les lignes à 99% - donne lieu à une erreur «la variable continue ne peut pas être mappée au type de ligne» avec ggplot2 actuel (0.9.3.1). Modifier geom_line(aes(x = number.seq, y = number.ll999, linetype = 2), data = dfCI)pour geom_line(aes(x = number.seq, y = number.ll999), linetype = 2, data = dfCI)travailler pour moi. N'hésitez pas à modifier la réponse originale et à la perdre.
2

Le code de Bernd Weiss est très utile. J'ai apporté quelques modifications ci-dessous, pour changer / ajouter quelques fonctionnalités:

  1. Erreur standard utilisée comme mesure de la précision, ce qui est plus typique des tracés en entonnoir que je vois (en psychologie)
  2. Inversé les axes, de sorte que la précision (erreur standard) est sur l'axe des y et la taille de l'effet sur l'axe des x
  3. Utilisé geom_segmentau lieu de geom_linepour la ligne délimitant la moyenne méta-analytique, de sorte qu'elle serait de la même hauteur que les lignes délimitant les régions de confiance à 95% et 99%
  4. Au lieu de tracer la moyenne méta-analytique, j'ai tracé son intervalle de confiance à 95%

Mon code utilise une moyenne méta-analytique de 0,0892 (se = 0,0035) à titre d'exemple, mais vous pouvez remplacer vos propres valeurs.

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

entrez la description de l'image ici

jsakaluk
la source