Puis-je supposer la normalité (log) pour cet échantillon?

11

Voici un tracé QQ pour mon échantillon (notez l'axe Y logarithmique); :n=1000

entrez la description de l'image ici
Comme l'a souligné whuber, cela indique que la distribution sous-jacente est asymétrique à gauche (la queue droite est plus courte).

En utilisant shapiro.test(sur les données transformées par log) dans R, j'obtiens une statistique de test de et une valeur de p de , ce qui signifie que nous rejetons formellement l'hypothèse nulle au niveau de confiance de 95%.W=0.97185.1721013H0:the sample is normal distributed

Ma question est la suivante: est-ce suffisamment bon dans la pratique pour une analyse plus approfondie en supposant la (log-) normalité? En particulier, je voudrais calculer les intervalles de confiance pour les moyennes d'échantillons similaires en utilisant la méthode approximative de Cox et Land (décrite dans l'article: Zou, GY, cindy Yan Huo et Taleban, J. (2009). Intervalles de confiance simples pour moyennes lognormales et leurs différences avec les applications environnementales. 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)))
}

J'ai remarqué que les intervalles de confiance ont tendance à être centrés autour d'un point légèrement supérieur à la moyenne réelle de l'échantillon. Par exemple:

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

Je pense que ces deux valeurs devraient être les mêmes sous .H0

Vegard
la source
1
La distribution ne correspond certainement pas bien à la queue droite.
Michael R. Chernick
1
Ce graphique QQ montre que les données ont une queue droite beaucoup plus courte qu'une distribution log-normale: elles sont asymétriques à gauche par rapport à une log-normale. Vous devez donc vous méfier de l'utilisation de procédures basées sur lognormal.
whuber
@whuber oui, vous avez raison de dire qu'il est asymétrique à gauche plutôt qu'à droite. Dois-je mettre à jour la question?
Vegard
Bien sûr: nous apprécions les améliorations apportées aux questions.
whuber
2
NB: veuillez noter que par "biais gauche", je voulais dire explicitement que la queue droite est courte, pas que la queue gauche est longue. Cela est évident par la façon dont les points à droite du tracé tombent sous la ligne de référence. Parce que les points à gauche du tracé sont (relativement) proches de la ligne de référence, il est incorrect de caractériser cette distribution comme ayant une "queue gauche plus longue". La distinction est importante ici, car la queue droite devrait avoir une influence beaucoup plus grande sur la moyenne estimée que la queue gauche (alors que les deux queues influencent son intervalle de confiance).
whuber

Réponses:

12

Ces données ont une queue courte par rapport à une distribution log-normale, un peu comme une distribution 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

Néanmoins, comme les données sont fortement asymétriques à droite, nous pouvons nous attendre à ce que les valeurs les plus élevées jouent un rôle important dans l'estimation de la moyenne et de son intervalle de confiance. Par conséquent, nous devons prévoir qu'un estimateur lognormal (LN) aura tendance à surestimer la moyenne et les deux limites de confiance .

Vérifions et, pour comparaison, utilisons les estimateurs habituels: c'est-à-dire la moyenne de l'échantillon et son intervalle de confiance en théorie normale. Il est à noter que les estimateurs habituels reposent uniquement sur la normalité approximative de la moyenne de l' échantillon , et non sur les données, et - avec un ensemble de données aussi volumineux - devraient fonctionner correctement. Pour ce faire, nous avons besoin d'une légère modification de la cifonction:

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

Voici une fonction parallèle pour les estimations de la théorie normale:

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

Appliquées à cet ensemble de données simulées, les sorties sont

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

ci.u1.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())

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

Histogrammes

Il est maintenant clair que les procédures lognormales ont tendance à surestimer la moyenne et les limites de confiance, alors que les procédures habituelles font du bon travail. Nous pouvons estimer les couvertures des procédures d'intervalle de confiance:

> 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 

Ce calcul dit:

  • La limite inférieure LN ne couvrira pas la vraie moyenne environ 22,3% du temps (au lieu des 2,5% prévus).

  • La limite inférieure habituelle ne couvrira pas la vraie moyenne environ 2,3% du temps, près des 2,5% prévus.

  • La limite supérieure LN dépassera toujours la vraie moyenne (au lieu de tomber en dessous de 2,5% du temps comme prévu). Cela en fait un intervalle de confiance bilatéral de 100% - (22,3% + 0%) = 77,7% au lieu d'un intervalle de confiance de 95%.

  • La limite supérieure habituelle ne couvrira pas la vraie moyenne environ 100 - 96,5 = 3,5% du temps. C'est un peu plus que la valeur prévue de 2,5%. Les limites habituelles comprennent donc un intervalle de confiance bilatéral de 100% - (2,3% + 3,5%) = 94,2% au lieu d'un intervalle de confiance de 95%.

La réduction de la couverture nominale de 95% à 77,7% pour l'intervalle lognormal est terrible. La réduction à 94,2% pour l'intervalle habituel n'est pas mal du tout et peut être attribuée à l'effet de l'asymétrie (des données brutes, pas de leurs logarithmes).

Nous devons conclure que de nouvelles analyses de la moyenne ne doivent pas supposer une lognormalité.

Faites attention! Certaines procédures (telles que les limites de prédiction) seront plus sensibles à l'asymétrie que ces limites de confiance pour la moyenne, il faudra donc peut-être tenir compte de leur distribution asymétrique. Cependant, il semble peu probable que les procédures lognormales fonctionnent bien avec ces données pour pratiquement n'importe quelle analyse prévue.

whuber
la source
Wow, cette réponse me souffle. Merci beaucoup! Comment se fait-il que vous utilisiez abline()au lieu de qqline()(ce qui produit une ligne différente) dans le premier exemple?
Vegard
Votre trial()fonction n'utilise pas ses arguments.
Vegard
1
Bon travail! Pour bootstrapping, modifier trial: trial <- function(y) { x <- sample(y, length(y), TRUE); cbind(ci(x), ci.u(x)) }. Ensuite , émettre juste une commande, sim <- sapply(1:5000, function(i) trial(x)). Vous voudrez peut-être explorer les histogrammes des six rangées d’ simaprès.
whuber
1
+1, j'aime particulièrement le fait subtil que les intervalles de prédiction seront plus sensibles à la forme distributionnelle que les intervalles de confiance pour la moyenne.
gung - Réintégrer Monica