Posso assumir (log-) normalidade para esta amostra?

11

Aqui está um gráfico de QQ para minha amostra (observe o eixo Y logarítmico); :n=1000

insira a descrição da imagem aqui
Conforme apontado pelo whuber, isso indica que a distribuição subjacente é inclinada para a esquerda (a cauda direita é mais curta).

Usando shapiro.test(nos dados transformados em log) em R, recebo uma estatística de teste e um valor p de 5,172 \ cdot10 ^ {- 13} , o que significa que rejeitamos formalmente a hipótese nula H_0: \ text { a amostra é distribuída normalmente} no nível de confiança de 95%.5,172 10 - 13W=0.97185.1721013H0:the sample is normal distributed

Minha pergunta é: isso é bom o suficiente na prática para análises posteriores assumindo (log-) normalidade? Em particular, eu gostaria de calcular intervalos de confiança para as médias de amostras semelhantes usando o método aproximado de Cox e Land (descrito no artigo: Zou, GY, cindy Yan Huo e Taleban, J. (2009). meios lognormal e suas diferenças com aplicações ambientais (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)))
}

Percebi que os intervalos de confiança tendem a se centrar em torno de um ponto ligeiramente acima da média real da amostra. Por exemplo:

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

Eu acho que esses dois valores devem ser os mesmos em .H0

Vegard
fonte
1
A distribuição definitivamente não se encaixa bem na cauda direita.
Michael R. Chernick
1
Esse gráfico QQ mostra que os dados têm uma cauda direita muito mais curta do que uma distribuição lognormal: é deixada inclinada em comparação com uma lognormal. Portanto, você deve desconfiar do uso de procedimentos baseados em lognormal.
whuber
@ whuber sim, você está certo sobre isso ficar inclinado em vez de inclinado à direita. Devo atualizar a pergunta?
Vegard
Claro: agradecemos as melhorias nas perguntas.
whuber
2
NB: observe que por "esquerda inclinada" eu quis dizer explicitamente que a cauda direita é curta, não que a cauda esquerda seja longa. Isso é evidente pela forma como os pontos à direita do gráfico ficam abaixo da linha de referência. Como os pontos à esquerda do gráfico estão (relativamente) próximos à linha de referência, é incorreto caracterizar essa distribuição como tendo uma "cauda esquerda mais longa". A distinção é importante aqui, porque a cauda direita deve ter uma influência muito maior na média estimada do que a cauda esquerda (enquanto as duas caudas influenciam seu intervalo de confiança).
whuber

Respostas:

12

Esses dados têm uma cauda curta em comparação com uma distribuição lognormal, não muito diferente de uma distribuição gama:

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

No entanto, como os dados são fortemente inclinados à direita, podemos esperar que os maiores valores desempenhem um papel importante na estimativa da média e seu intervalo de confiança. Portanto , devemos antecipar que um estimador lognormal (LN) tenderá a superestimar a média e os dois limites de confiança .

Vamos verificar e, para comparação, usar os estimadores usuais: ou seja, a média da amostra e seu intervalo de confiança da teoria normal. Observe que os estimadores usuais dependem apenas da normalidade aproximada da média da amostra , não dos dados, e - com um conjunto de dados tão grande - pode funcionar bem. Para fazer isso, precisamos de uma ligeira modificação da cifunção:

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

Aqui está uma função paralela para as estimativas da teoria 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 esse conjunto de dados simulado, as saídas são

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

As estimativas da teoria normal produzidas ci.uparecem um pouco mais próximas da média real de , mas é difícil distinguir de um conjunto de dados qual procedimento tende a funcionar melhor. Para descobrir, vamos simular muitos conjuntos de dados: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 interessados ​​em comparar os resultados com a média real de . Um painel de histogramas é revelador a esse respeito: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

Agora está claro que os procedimentos lognormal tendem a superestimar a média e os limites de confiança, enquanto os procedimentos usuais fazem um bom trabalho. Podemos estimar as coberturas dos procedimentos de intervalo de confiança:

> 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 diz:

  • O limite inferior do LN não cobrirá a média verdadeira em cerca de 22,3% do tempo (em vez dos 2,5% pretendidos).

  • O limite inferior usual falhará em cobrir a média real cerca de 2,3% das vezes, próximo aos 2,5% pretendidos.

  • O limite superior do LN sempre excederá a média verdadeira (em vez de ficar abaixo dele 2,5% do tempo, conforme o planejado). Isso o torna um intervalo de confiança de 100% - (22,3% + 0%) = 77,7% em vez de um intervalo de confiança de 95%.

  • O limite superior usual falhará em cobrir a média verdadeira cerca de 100 - 96,5 = 3,5% do tempo. Isso é um pouco maior que o valor pretendido de 2,5%. Os limites usuais, portanto, compreendem um intervalo de confiança de 100% - (2,3% + 3,5%) = 94,2% nos dois lados, em vez de um intervalo de confiança de 95%.

A redução da cobertura nominal de 95% para 77,7% para o intervalo lognormal é terrível. A redução para 94,2% para o intervalo usual não é ruim e pode ser atribuída ao efeito da assimetria (dos dados brutos, não de seus logaritmos).

Temos que concluir que análises adicionais da média não devem assumir normalidade do logaritmo.

Seja cuidadoso! Alguns procedimentos (como limites de previsão) serão mais sensíveis à assimetria do que esses limites de confiança para a média, portanto, sua distribuição distorcida pode precisar ser considerada. No entanto, parece improvável que os procedimentos lognormal tenham um bom desempenho com esses dados para praticamente qualquer análise pretendida.

whuber
fonte
Uau, essa resposta me surpreende. Muito obrigado! Como você usa em abline()vez de qqline()(que produz uma linha diferente) no primeiro exemplo?
Vegard
Sua trial()função não usa seus argumentos.
Vegard
1
Bom trabalho! Para bootstrapping, modificar trial: trial <- function(y) { x <- sample(y, length(y), TRUE); cbind(ci(x), ci.u(x)) }. Em seguida, emita apenas um comando sim <- sapply(1:5000, function(i) trial(x)). Você pode explorar os histogramas das seis linhas simseguintes.
whuber
1
+1, gosto especialmente do ponto sutil de que os intervalos de previsão serão mais sensíveis à forma distributiva do que os intervalos de confiança para a média.
gung - Restabelece Monica