¿Cómo puedo estimar la densidad de un parámetro inflado a cero en R?

Question

Tengo un conjunto de datos con muchos ceros que tiene este aspecto:

set.seed(1)
x <- c(rlnorm(100),rep(0,50))
hist(x,probability=TRUE,breaks = 25)

Me gustaría trazar una línea para su densidad, pero el density() utiliza una ventana móvil que calcula los valores negativos de x.

lines(density(x), col = 'grey')

Hay una density(... from, to) pero parece que sólo truncan el cálculo, no alteran la ventana para que la densidad en 0 sea coherente con los datos, como puede verse en el siguiente gráfico:

lines(density(x, from = 0), col = 'black')

(si se cambiara la interpolación, esperaría que la línea negra tuviera mayor densidad en 0 que la línea gris)

¿Existen alternativas a esta función que proporcionen un mejor cálculo de la densidad en cero?

Answer 1

Estoy de acuerdo con Rob Hyndman en que hay que tratar los ceros por separado. Hay algunos métodos para tratar la estimación de la densidad del núcleo de una variable con soporte limitado, incluyendo la "reflexión", la "rernormalización" y la "combinación lineal". Estos no parecen haber sido implementados en el programa R density pero están disponibles en Benn Jann's kdens para Stata .

Answer 2

Puede intentar reducir el ancho de banda (la línea azul es para adjust=0.5 ),

pero probablemente KDE no es el mejor método para tratar esos datos.

Answer 3

La densidad es infinita en cero porque incluye un pico discreto. Hay que estimar el pico utilizando la proporción de ceros, y luego estimar la parte positiva de la densidad asumiendo que es suave. La KDE causará problemas en el extremo izquierdo porque dará cierto peso a los valores negativos. Un enfoque útil es transformar a logaritmos, estimar la densidad utilizando KDE, y luego volver a transformar. Véase Wand, Marron & Ruppert (JASA 1991) para una referencia.

La siguiente función de R hará la densidad transformada:

logdensity <- function (x, bw = "SJ") 
{
    y <- log(x)
    g <- density(y, bw = bw, n = 1001)
    xgrid <- exp(g$x)
    g$y <- c(0, g$y/xgrid)
    g$x <- c(0, xgrid)
    return(g)
}

Entonces lo siguiente dará la trama que quieres:

set.seed(1)
x <- c(rlnorm(100),rep(0,50))
hist(x,probability=TRUE,breaks = 25)
fit <- logdensity(x[x>0]) # Only take density of positive part
lines(fit$x,fit$y*mean(x>0),col="red") # Scale density by proportion positive
abline(v=0,col="blue") # Add spike at zero.

Answer 4

Otra opción cuando se tienen datos con un límite inferior lógico (como 0, pero podrían ser otros valores) que se sabe que los datos no bajarán y la estimación regular de la densidad del núcleo sitúa los valores por debajo de ese límite (o si se tiene un límite superior, o ambos) es utilizar estimaciones logspline. El paquete logspline para R las implementa y las funciones tienen argumentos para especificar los límites, de modo que la estimación irá hasta el límite, pero no más allá, y seguirá escalando a 1.

También existen métodos (el oldlogspline ) que tendrá en cuenta la censura por intervalos, de modo que si esos 0 no son 0 exactos, sino que se redondean de modo que se sabe que representan valores entre 0 y algún otro número (un límite de detección, por ejemplo), entonces se puede dar esa información a la función de ajuste.

Si los 0 extra son verdaderos 0 (no redondeados) entonces la estimación de la masa de picos o puntos es el mejor enfoque, pero también se puede combinar con la estimación logspline.