¿Puedo asumir (log-) normalidad para esta muestra?

11

Aquí hay un gráfico QQ para mi muestra (observe el eje Y logarítmico); :n=1000

ingrese la descripción de la imagen aquí
Como señaló Whuber, esto indica que la distribución subyacente está sesgada hacia la izquierda (la cola derecha es más corta).

Utilizando shapiro.test(sobre los datos transformados logarítmicamente) en R, tengo una prueba estadística de y un valor de p de , lo que significa que formalmente rechazar la hipótesis nula al nivel de confianza del 95%.5.172 10 - 13 H 0 : la muestra está distribuida normalmenteW=0.97185.1721013H0:the sample is normal distributed

Mi pregunta es: ¿es esto lo suficientemente bueno en la práctica para un análisis posterior asumiendo (log-) normalidad? En particular, me gustaría calcular los intervalos de confianza para las medias de muestras similares utilizando el método aproximado de Cox y Land (descrito en el documento: Zou, GY, cindy Yan Huo y Taleban, J. (2009). Intervalos de confianza simples para medios lognormales y sus diferencias con las aplicaciones ambientales. Environmetrics 20, 172–180):

ci <- function (x) {
        y <- log(x)
        n <- length(y)
        s2 <- var(y)
        m <- mean(y) + s2 / 2
        z <- qnorm(1 - 0.05 / 2) # 95%
        #z <- qnorm(1 - 0.10 / 2) # 90%
        d <- z * sqrt(s2 / n + s2 * s2 / (2 * (n - 1)))

        return(c(exp(m - d), exp(m + d)))
}

Me di cuenta de que los intervalos de confianza tienden a centrarse en un punto que está ligeramente por encima de la media muestral real. Por ejemplo:

> mean(x)
[1] 82.3076
> y <- log(x)
> exp(mean(y) + var(y) / 2)
[1] 91.22831

Creo que estos dos valores deberían ser iguales bajo .H0

Vegard
fuente
1
La distribución definitivamente no encaja bien en la cola derecha.
Michael R. Chernick
1
Este gráfico QQ muestra que los datos tienen una cola derecha mucho más corta que una distribución lognormal: se deja sesgada en comparación con una lognormal. Por lo tanto, debe desconfiar del uso de procedimientos basados ​​en lognormal.
whuber
@whuber sí, tienes razón acerca de que quede sesgado en lugar de estar sesgado a la derecha. ¿Debo actualizar la pregunta?
Vegard
Claro: apreciamos las mejoras a las preguntas.
whuber
2
NB: tenga en cuenta que por "sesgado a la izquierda" me refería explícitamente a que la cola derecha es corta, no que la cola izquierda sea larga. Esto es evidente por cómo los puntos a la derecha del diagrama caen debajo de la línea de referencia. Debido a que los puntos a la izquierda de la gráfica están (relativamente) cerca de la línea de referencia, es incorrecto caracterizar esta distribución como que tiene una "cola izquierda más larga". La distinción es importante aquí, porque la cola derecha debería tener una influencia mucho mayor en la media estimada que la cola izquierda (mientras que ambas colas influyen en su intervalo de confianza).
whuber

Respuestas:

12

Estos datos tienen una cola corta en comparación con una distribución lognormal, no muy diferente de una distribución Gamma:

set.seed(17)
par(mfcol=c(1,1))
x <- rgamma(500, 1.9)
qqnorm(log(x), pch=20, cex=.8, asp=1)
abline(mean(log(x)) + .1,1.2*sd(log(x)), col="Gray", lwd=2)

QQPlot

Sin embargo, debido a que los datos están fuertemente sesgados, podemos esperar que los valores más grandes desempeñen un papel importante en la estimación de la media y su intervalo de confianza. Por lo tanto , debemos anticipar que un estimador lognormal (LN) tenderá a sobreestimar la media y los dos límites de confianza .

Verifiquemos y, para comparación, usemos los estimadores habituales: es decir, la media muestral y su intervalo de confianza de la teoría normal. Tenga en cuenta que los estimadores habituales se basan únicamente en la normalidad aproximada de la media de la muestra , no de los datos, y, con un conjunto de datos tan grande, se puede esperar que funcionen bien. Para hacer esto, necesitamos una ligera modificación de la cifunción:

ci <- function (x, alpha=.05) {
  z <- -qnorm(alpha / 2)
  y <- log(x); n <- length(y); s2 <- var(y)
  m <- mean(y) + s2 / 2
  d <- z * sqrt(s2 / n + s2 * s2 / (2 * (n - 1)))
  exp(c(mean=m, lcl=m-d, ucl=m+d))
}

Aquí hay una función paralela para las estimaciones de la teoría normal:

ci.u <- function(x, alpha=.05) {
 mean(x) + sd(x) * c(mean=0, lcl=1, ucl=-1) / sqrt(length(x)) * qnorm(alpha/2)
}

Aplicado a este conjunto de datos simulado, las salidas son

> ci(x)
   mean     lcl     ucl 
2.03965 1.87712 2.21626 
> ci.u(x)
   mean     lcl     ucl 
1.94301 1.81382 2.07219 

Las estimaciones de la teoría normal producidas por ci.umirar un poco más cerca de la media real de , pero es difícil saber a partir de un conjunto de datos qué procedimiento tiende a funcionar mejor. Para averiguarlo, simulemos muchos conjuntos de datos:1.9

trial <- function(n=500, k=1.9) {
  x <- rgamma(n, k)
  cbind(ci(x), ci.u(x))
}
set.seed(17)
sim <- replicate(5000, trial())

Estamos interesados ​​en comparar los resultados con la media real de . Un panel de histogramas es revelador a ese respecto:1.9

xmin <- min(sim)
xmax <- max(sim)
h <- function(i, ...) {
  b <- seq(from=floor(xmin*10)/10, to=ceiling(xmax*10)/10, by=0.1)
  hist(sim[i,], freq=TRUE, breaks=b, col="#a0a0FF", xlab="x", xlim=c(xmin, xmax), ...)
  hist(sim[i,sim[i,] >= 1.9], add=TRUE,freq=TRUE, breaks=b, col="#FFa0a0",
                              xlab="x", xlim=c(xmin, xmax), ...)
}
par(mfcol=c(2,3))
h(1, main="LN Estimate of Mean")
h(4, main="Sample Mean")
h(2, main="LN LCL")
h(5, main="LCL")
h(3, main="LN UCL")
h(6, main="UCL")

Histogramas

Ahora está claro que los procedimientos lognormales tienden a sobreestimar la media y los límites de confianza, mientras que los procedimientos habituales hacen un buen trabajo. Podemos estimar las coberturas de los procedimientos de intervalo de confianza:

> sapply(c(LNLCL=2, LCL=5, LNUCL=3, UCL=6), function(i) sum(sim[i,] > 1.9)/dim(sim)[2])
 LNLCL    LCL  LNUCL    UCL 
0.2230 0.0234 1.0000 0.9648 

Este cálculo dice:

  • El límite inferior de LN no cubrirá la media real aproximadamente el 22,3% del tiempo (en lugar del 2,5% previsto).

  • El límite inferior habitual no cubrirá la media real aproximadamente el 2,3% del tiempo, cerca del 2,5% previsto.

  • El límite superior de LN siempre excederá la media real (en lugar de caer por debajo del 2.5% del tiempo según lo previsto). Esto lo convierte en un 100% de dos lados - (22.3% + 0%) = 77.7% de intervalo de confianza en lugar de un intervalo de confianza de 95%.

  • El límite superior habitual no cubrirá la media verdadera aproximadamente 100 - 96.5 = 3.5% del tiempo. Esto es un poco mayor que el valor previsto de 2.5%. Por lo tanto, los límites habituales comprenden un 100% de dos lados - (2.3% + 3.5%) = 94.2% de intervalo de confianza en lugar de un intervalo de confianza de 95%.

La reducción de la cobertura nominal del 95% al ​​77,7% para el intervalo lognormal es terrible. La reducción al 94.2% para el intervalo habitual no es mala en absoluto y puede atribuirse al efecto de la asimetría (de los datos en bruto, no de sus logaritmos).

Tenemos que concluir que los análisis posteriores de la media no deben suponer lognormalidad.

¡Ten cuidado! Algunos procedimientos (como los límites de predicción) serán más sensibles a la asimetría que estos límites de confianza para la media, por lo que es posible que se deba tener en cuenta su distribución sesgada. Sin embargo, parece poco probable que los procedimientos logarítmicos funcionen bien con estos datos para prácticamente cualquier análisis previsto.

whuber
fuente
Wow, esta respuesta me deja sin aliento. Muchas gracias! ¿Cómo es que usas en abline()lugar de qqline()(que produce una línea diferente) en el primer ejemplo?
Vegard
Su trial()función no usa sus argumentos.
Vegard
1
El código que usé .
Vegard
1
¡Buen trabajo! Para el arranque, modificar trial: trial <- function(y) { x <- sample(y, length(y), TRUE); cbind(ci(x), ci.u(x)) }. A continuación, emita solo un comando, sim <- sapply(1:5000, function(i) trial(x)). Es posible que desee explorar los histogramas de las seis filas de más simadelante.
whuber
1
+1, me gusta especialmente el punto sutil de que los intervalos de predicción serán más sensibles a la forma de distribución que los intervalos de confianza para la media.
gung - Restablece a Monica