Могу ли я принять (лог) нормальность для этого образца?

11

Вот график QQ для моего образца (обратите внимание на логарифмическую ось Y); :Nзнак равно1000

введите описание изображения здесь
Как указывает whuber, это указывает на то, что базовое распределение перекошено влево (правый хвост короче).

Используя shapiro.test(на лог-преобразованных данных) в R, я получаю тестовую статистику и p-значение , что означает, что мы формально отвергаем нулевую гипотезу на уровне достоверности 95%.Wзнак равно0,97185,17210-13ЧАС0:образец нормально распределен

Мой вопрос: достаточно ли это на практике для дальнейшего анализа, предполагающего (лог) нормальность? В частности, я хотел бы рассчитать доверительные интервалы для средних значений схожих выборок, используя приближенный метод Кокса и Лэнда (описанный в статье: Zou, GY, Cindy Yan Huo and Taleban, J. (2009). Простые доверительные интервалы для Логнормальные средства и их отличия от экологических приложений. 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)))
}

Я заметил, что доверительные интервалы, как правило, сосредоточены вокруг точки, которая немного выше среднего значения фактической выборки. Например:

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

Я думаю, что эти два значения должны быть одинаковыми при .ЧАС0

Vegard
источник
1
Распределение определенно не вписывается в правильный хвост.
Майкл Р. Черник
1
Этот график QQ показывает, что данные имеют гораздо более короткий правый хвост, чем логнормальное распределение: они отклонены влево по сравнению с логнормальным. Поэтому вам следует опасаться использования процедур, основанных на логарифмических нормах.
whuber
@whuber: да, вы правы, когда вы наклонены, а не наклонены вправо. Должен ли я обновить вопрос?
Вегард
Конечно: мы ценим улучшения вопросов.
whuber
2
NB: пожалуйста, обратите внимание, что под «перекосом влево» я имел в виду, что правый хвост короткий, а не длинный. Это видно по тому, как точки справа от графика падают ниже контрольной линии. Поскольку точки в левой части графика (относительно) близки к контрольной линии, неправильно характеризовать это распределение как имеющее «более длинный левый хвост». Это различие важно здесь, потому что правый хвост должен оказывать гораздо большее влияние на оценочное среднее значение, чем левый хвост (тогда как оба хвоста влияют на его доверительный интервал).
whuber

Ответы:

12

Эти данные имеют короткий хвост по сравнению с логнормальным распределением, мало чем отличающимся от гамма-распределения:

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

Тем не менее, поскольку данные являются сильно правой перекос, мы можем ожидать , что наибольшие значения играют важную роль в оценке среднего значения и его доверительный интервал. Поэтому следует ожидать, что логнормальная (LN) оценка будет склонна переоценивать среднее значение и два доверительных интервала .

Давайте проверим и, для сравнения, используем обычные оценки: то есть среднее значение выборки и его доверительный интервал в нормальной теории. Обратите внимание, что обычные оценщики полагаются только на приблизительную нормальность среднего значения выборки , а не данных, и - при таком большом наборе данных - можно ожидать, что они будут работать хорошо. Для этого нам понадобится небольшая модификация ciфункции:

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

Вот параллельная функция для нормальных оценок:

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

Применительно к этому смоделированному набору данных, выходы

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

Гистограммы

Теперь ясно, что логнормальные процедуры имеют тенденцию переоценивать среднее значение и пределы достоверности, тогда как обычные процедуры делают хорошую работу. Мы можем оценить покрытия процедур доверительного интервала:

> 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 

Этот расчет говорит:

  • Нижний предел LN не сможет охватить истинное среднее значение примерно в 22,3% времени (вместо запланированных 2,5%).

  • Обычный нижний предел не сможет охватить истинное среднее значение около 2,3% времени, что близко к предполагаемому 2,5%.

  • Верхний предел LN всегда будет превышать истинное среднее значение (вместо того, чтобы опускаться ниже 2,5% времени, как предполагалось). Это делает двусторонний 100% - (22,3% + 0%) = 77,7% доверительный интервал вместо 95% доверительного интервала.

  • Обычный верхний предел не сможет охватить истинное среднее значение в 100 - 96,5 = 3,5% времени. Это немного больше, чем предполагаемое значение 2,5%. Таким образом, обычные пределы включают двусторонний доверительный интервал 100% - (2,3% + 3,5%) = 94,2% вместо 95% доверительного интервала.

Сокращение номинального покрытия с 95% до 77,7% для логнормального интервала ужасно. Снижение до 94,2% для обычного интервала совсем не плохо и может быть объяснено влиянием асимметрии (необработанных данных, а не их логарифмов).

Мы должны сделать вывод, что дальнейший анализ среднего значения не должен предполагать логнормальность.

Быть осторожен! Некоторые процедуры (такие как пределы прогнозирования) будут более чувствительными к асимметрии, чем эти доверительные пределы для среднего значения, поэтому, возможно, придется учитывать их искаженное распределение. Однако маловероятно, что логнормальные процедуры будут хорошо работать с этими данными практически для любого предполагаемого анализа.

Whuber
источник
Ух ты, этот ответ поражает меня. Огромное спасибо! Почему вы используете abline()вместо qqline()(который производит другую строку) в первом примере?
Вегард
Ваша trial()функция не использует свои аргументы.
Вегард
1
Хорошая работа! Для самозагрузки, изменить trial: trial <- function(y) { x <- sample(y, length(y), TRUE); cbind(ci(x), ci.u(x)) }. Затем выполните только одну команду sim <- sapply(1:5000, function(i) trial(x)). Возможно, вы захотите изучить гистограммы шести рядов simвпоследствии.
whuber
1
+1, мне особенно нравится тонкий момент, что интервалы прогнозирования будут более чувствительными к форме распределения, чем доверительные интервалы для среднего значения.
gung - Восстановить Монику