Bonnes méthodes pour les diagrammes de densité de variables non négatives dans R?

36
plot(density(rexp(100))

De toute évidence, toute densité à gauche de zéro représente un biais.

Je cherche à résumer certaines données relatives aux non-statisticiens et à éviter de se demander pourquoi les données non négatives ont une densité inférieure à zéro. Les parcelles sont destinées à la vérification de la randomisation; Je veux montrer la distribution des variables par traitement et par groupe de contrôle. Les distributions sont souvent exponentielles. Les histogrammes sont délicats pour diverses raisons.

Une recherche rapide sur Google me donne le travail de statisticiens sur les noyaux non négatifs, par exemple: this .

Mais est-ce que tout cela a été implémenté dans R? Parmi les méthodes mises en œuvre, certaines d'entre elles sont-elles "meilleures" d'une manière ou d'une autre pour les statistiques descriptives?

EDIT: même si la fromcommande peut résoudre mon problème actuel, il serait bon de savoir si quelqu'un a implémenté des noyaux basés sur la littérature sur l'estimation de densité non négative

utilisateur_générique
la source
3
Ce n’est pas ce que vous demandez, mais je n’appliquerais pas l’estimation de la densité du noyau à quelque chose qui devrait être exponentiel, en particulier pour une présentation à un public non statistique. J'utiliserais un graphique quantile-quantile et expliquerais que le graphique devrait être rectiligne si la distribution est exponentielle.
Nick Cox
6
plot(density(rexp(100), from=0))?
Stéphane Laurent
4
Une chose que j’ai parfois fait assez bien est d’obtenir un kde sur les journaux, puis de transformer l’estimation de la densité (sans oublier le jacobien). Une autre possibilité consisterait à utiliser une estimation de la densité log-spline configurée afin de connaître les limites.
Glen_b -Reinstate Monica
1
J'ai discuté de la méthode de transformation mentionnée par @Glen_b dans stata-journal.com/sjpdf.html?articlenum=gr0003 (voir p. 76-78). Les zéros peuvent être adaptés en utilisant log (x + 1) plutôt que de se connecter et en modifiant le jacobien.
Nick Cox

Réponses:

21

Une solution, empruntée aux méthodes de pondération des statistiques spatiales, consiste à tronquer la densité de gauche à zéro mais à pondérer les données les plus proches de zéro. L'idée est que chaque valeur est "étalée" dans un noyau d'aire totale d'unité centrée sur x ; toute partie du noyau qui déborderait en territoire négatif est supprimée et le noyau est renormalisé en unité de surface.XX

Par exemple, avec un noyau gaussien, , le poids de renormalisation estKh(y,X)=exp(-12((y-X)/h)2)/2π

w(X)=1/0K(y,X)y=11-ΦX,h(0)

est la fonction de distribution cumulative d'une variable normale de la moyenne x et de l'écart type h . Des formules comparables sont disponibles pour les autres noyaux.ΦXh

C’est plus simple - et beaucoup plus rapide dans les calculs - que d’essayer de réduire la largeur de bande proche de . Il est difficile de prescrire exactement comment les largeurs de bande devraient être modifiées près de 0 , de toute façon. Néanmoins, cette méthode est également ad hoc : il restera toujours un biais près de 0 . Il semble fonctionner mieux que l'estimation de densité par défaut. Voici une comparaison utilisant un jeu de données assez long:000

Figure

Le bleu indique la densité par défaut tandis que le rouge indique la densité ajustée pour le bord à . La vraie distribution sous-jacente est tracée en trait pointillé pour référence.0


Code R

La densityfonction de Rse plaindra que la somme des poids n'est pas unitaire, car elle veut que l'intégrale de tous les nombres réels soit unitaire, alors que cette approche rend l'intégrale sur des nombres positifs égale à l'unité. À titre de contrôle, cette dernière intégrale est estimée sous forme de somme de Riemann.

set.seed(17)
x <- rexp(1000)
#
# Compute a bandwidth.
#
h <- density(x, kernel="gaussian")$bw # $
#
# Compute edge weights.
#
w <- 1 / pnorm(0, mean=x, sd=h, lower.tail=FALSE)
#
# The truncated weighted density is what we want.
#
d <- density(x, bw=h, kernel="gaussian", weights=w / length(x))
d$y[d$x < 0] <- 0
#
# Check: the integral ought to be close to 1:
#
sum(d$y * diff(d$x)[1])
#
# Plot the two density estimates.
#
par(mfrow=c(1,1))
plot(d, type="n", main="Default and truncated densities", xlim=c(-1, 5))
polygon(density(x, kernel="gaussian", bw=h), col="#6060ff80", border=NA)
polygon(d, col="#ff606080", border=NA)
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
whuber
la source
21

Une alternative est l'approche de Kooperberg et ses collègues, basée sur l'estimation de la densité en utilisant des splines pour approximer la densité logarithmique des données. Je vais donner un exemple en utilisant les données de la réponse de @ whuber, ce qui permettra une comparaison des approches.

set.seed(17)
x <- rexp(1000)

Pour cela, vous aurez besoin du paquet logspline . installez-le s'il ne l'est pas:

install.packages("logspline")

Chargez le package et estimez la densité à l'aide de la logspline()fonction:

require("logspline")
m <- logspline(x)

Dans ce qui suit, je suppose que l'objet dde la réponse de @ whuber est présent dans l'espace de travail.

plot(d, type="n", main="Default, truncated, and logspline densities", 
     xlim=c(-1, 5), ylim = c(0, 1))
polygon(density(x, kernel="gaussian", bw=h), col="#6060ff80", border=NA)
polygon(d, col="#ff606080", border=NA)
plot(m, add = TRUE, col = "red", lwd = 3, xlim = c(-0.001, max(x)))
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
rug(x, side = 3)

Le tracé obtenu est présenté ci-dessous, avec la densité de la log en boucle indiquée par la ligne rouge.

Densités par défaut, tronquées et log-logs

De plus, la prise en charge de la densité peut être spécifiée via des arguments lboundet ubound. Si nous voulons supposer que la densité est 0 à gauche de 0 et qu'il existe une discontinuité à 0, nous pourrions utiliser lbound = 0dans l'appel à logspline(), par exemple

m2 <- logspline(x, lbound = 0)

Production de l’estimation de densité suivante (illustrée ici avec l’ majustement initial de la ligne de journal car la figure précédente était déjà occupée).

plot.new()
plot.window(xlim = c(-1, max(x)), ylim = c(0, 1.2))
title(main = "Logspline densities with & without a lower bound",
      ylab = "Density", xlab = "x")
plot(m,  col = "red",  xlim = c(0, max(x)), lwd = 3, add = TRUE)
plot(m2, col = "blue", xlim = c(0, max(x)), lwd = 2, add = TRUE)
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
rug(x, side = 3)
axis(1)
axis(2)
box()

Le graphique résultant est présenté ci-dessous

Comparaison des estimations de densité de logspline avec et sans limite inférieure sur le support

xX=0x

Rétablir Monica - G. Simpson
la source
1
01
@ Whuber Bonne question. Je ne suis tombé sur cette approche que récemment. Je suppose qu’une bonne question à poser ici est la suivante: étant donné que les méthodes tronquée et logpline ne sont que des estimations de la densité réelle, les différences d’ajustement sont-elles significatives sur le plan statistique? Je ne sais pas exactement pourquoi il fait si bien à zéro, cependant. J'apprécierais savoir pourquoi aussi.
Réintégrer Monica - G. Simpson le
@GavinSimpson, Merci pour cette belle réponse. Pouvez-vous reproduire la dernière parcelle avec la dernière version de logspline? Pour moi, la densité des deux versions, la version bornée et la version non bornée va à zéro à x = 0.
cel
4

Pour comparer les distributions par groupes (dont vous dites que l'objectif est l'un de vos commentaires), pourquoi ne pas quelque chose de plus simple? Les boîtes à moustaches parallèles fonctionnent bien si N est grand; Les tracés en bande parallèles fonctionnent si N est petit (et les deux montrent bien les valeurs aberrantes, ce que vous dites être un problème dans vos données).

Peter Flom - Rétablir Monica
la source
1
Oui, merci, ça marche. Mais j'aime les parcelles de densité. Ils montrent plus sur les données que les boîtes à moustaches. Je suppose que je suis un peu surpris que rien ne semble avoir déjà été mis en œuvre. Peut-être que je mettrai moi-même en œuvre l'une de ces choses un jour. Les gens le trouveraient probablement utile.
generic_user
1
J'aime aussi les parcelles de densité; mais vous devez considérer votre public.
Peter Flom - Rétablir Monica
1
Je suis d'accord avec @PeterFlom sur celui-ci. Ne soyez pas trop compliqué si votre public n’est pas statistiquement informé. Vous pouvez également créer des boîtes à moustaches comparatives / parallèles avec une superposition de parcelles à papillons. Ainsi, le résumé de la boîte à moustaches est visible, ainsi que toutes les données.
Doug.numbers
La suggestion que différentes personnes comprennent différemment les parcelles agrégées est certainement correcte. Bien que je comprenne ce qu'est un diagramme de densité (et que ce ne soit pas une probabilité), je ne comprends pas ce que pourrait être un "diagramme à barres parallèle". Cela suggère un tracé de coordonnées parallèles mais je soupçonne que ce n’est pas correct.
DWin
2

Comme Stéphane commente, vous pouvez utiliser from = 0et, en plus, vous pouvez représenter vos valeurs sous la courbe de densité avecrug (x)

Aghila
la source
4
Corrigez-moi si je me trompe, mais from=0semble supprimer le tracé pour les valeurs inférieures à 0; cela ne corrige pas le calcul du fait qu'une partie de la distribution a été passée en dessous de 0.
Nick Cox
1
C'est correct. L'utilisation de la fromcommande génère un tracé qui semble avoir le pic juste à zéro. Mais si vous examinez les histogrammes avec des bacs toujours plus petits, de nombreuses données indiquent le pic AT zéro. La fromest juste une astuce graphique.
generic_user
@ NickCox Je ne suis pas sûr, mais je ne pense pas qu'il from=0supprime quoi que ce soit. Il commence juste la "grille" à zéro.
Stéphane Laurent
La différence est de savoir si la densité estimée est non nulle pour les valeurs négatives, et non si elle est tracée ou non. Les chercheurs peuvent décider de ne pas s’inquiéter à ce sujet si tout ce qu’ils veulent, c’est une visualisation.
Nick Cox
@NickCox La commande density(rexp(100), from=0)n'a rien à voir avec le graphique
Stéphane Laurent