Регрессия: доверительный интервал? Опции

[stok1946]

Насколько я помню из теории, различают два вида доверительных областей:
а) Доверительная область для линии регрессии RD(x), (т.е. точнее, для прогнозов модели)
б) Доверительная область для самих значений зависимой переменной YD(x)

С понятием а) все понятно: если многократно извлекать из генеральной совокупности различные выборки из N пар (x,y) значений и строить по ним модели регрессии, то за пределами "доверительной трубы" окажется 100alfa% таких линий (alfa, например, равно 0.05).
Назовем эту "трубу" RD(x), поскольку она зависит от текущего значения х. Эту самую RD(x) можно лихо и точно просчитать бутстрепом для самых различных функций, включая сплайны и ядерные.

Теперь относительно ДИ под буковкой б). Это - интервал, определяющий границы, за пределами которых могут оказаться не более 100alfa% экспериментальных точек наблюдений при Х = х.
Он, вообще говоря (как пишет, в частности Гланц на стр. 243), складывается из разброса значений вокруг линии регресии и неопределенности положения самой этой линии (второе мы уже посчитали как RD(x)).
Характеристикой разброса значений y вокруг линии регрессии является только остаточное стандартное отклонение sy|x.
Я могу ширину этой части доверительного интервала оценить только по эмпирической выборке SD = t(1-alfa/2, N-2) * sqr(RSS/(N - K)), где RSS - сумма квадратов остатков. И это - постоянная величина на всем интервале определения независимой переменной.
И тут такие вопросы:
а) когда в STATISTICA и др. прикладных программах считают доверительные интервалы, то что имеют в виду RD(x), YD(x) или SD?
Например, Гайдышев в Approximations (APX) AtteStat выводит SD и называет это "Доверительные интервалы оценок модели" (но ведь "оценки" - это и есть прогнозы модели?);
б) справедливо ли выражение , YD(x) = RD(x) + SD?
в) почему на картинке к
http://ru.wikipedia.org/wiki/%D0%A0%D0%B5%...%D0%BA%D0%B0%29
95%-е доверительные области для линии регрессии RD(x) показаны криволинейными, и для значений YD(x) - прямыми красненькими линиями, хотя, согласно там же приведенных формул, YD(x) включает RD(x) .
Или мне это мерещится?

Сообщение отредактировал stok1946 - 6.03.2013 - 13:08

[nokh]

Насколько я помню из теории, различают два вида доверительных областей:
а) Доверительная область для линии регрессии RD(x), (т.е. точнее, для предикторов модели)
б) Доверительная область для самих значений зависимой переменной YD(x)
...
И тут такие вопросы:
а) когда в STATISTICA и др. прикладных программах считают доверительные интервалы, то что имеют в виду RD(x), YD(x) или SD?
Например, Гайдышев в Approximations (APX) AtteStat выводит SD и называет это "Доверительные интервалы оценок модели" (но ведь "оценки" - это как раз и есть предикторы модели?);
б) справедливо ли выражение , YD(x) = RD(x) + SD?
в) почему на картинке к
http://ru.wikipedia.org/wiki/%D0%A0%D0%B5%...%D0%BA%D0%B0%29
95%-е доверительные области для линии регрессии RD(x) показаны криволинейными, и для значений YD(x) - прямыми красненькими линиями, хотя, согласно там же приведенных формул, YD(x) включает RD(x) .
Или мне это мерещится?

Мне представляется более удачным рассматривать (а) ДИ для регрессии и (б) ДИ для прогноза (так короче и часто в литературе речь и идёт о prediction). Первые указывают границы для среднего значения y при данном х, вторые - для отдельного значения у при данном х. Бегло посмотрел Гланца и Sokal & Rohlf (http://yadi.sk/d/j-B_6FDx2dp6y )- вроде всё понятно и об одном и том же.
Вопросы.
(a) Разные пакеты под одним и тем же термином могут понимать разные вещи. Я с ужасом это обнаружил когда разбирался с выбросами и влияющими наблюдениями в регрессии ( http://forum.disser.ru/index.php?showtopic=2434&st=71 , сообщение #72). В принципе, в Statistica - неплохая помощь, но тогда пришлось просто просчитывать пример руками и в пакетах, чтобы понять что где выдаётся. В данном случае про ДИ всё написано. Если в Help в Указателе набить "Confidence Interval vs. Prediction Interval" и в нижней части открывшейся страницы уйти по ссылке "Multiple Regression Results - Residuals/Assumptions/Prediction Tab" то можно увидеть формулы по которым считает пакет и ссылку на Neter, Wasserman, & Kutner, 1985. Учитывая что программисты чётко различают между Confidence Interval и Prediction Interval, есть основания полагать, что в качестве ДИ для регрессии выдаётся именно Confidence Interval. Но лучшая проверка - сверка с ручным расчётом или примером из книги.
(б) не знаю, формулы везде разные; чтобы привести их к единообразнию нужно поработать.
(в) и ДИ для регрессии, и ДИ для прогноза должны быть криволинейными, демонстрирующими нарастание ошибки от центра системы к периферии. В примере Википедии степень кривизны красных линий мала, но она есть - прикладывал линейку. Выходит - мерещится .

Сообщение отредактировал nokh - 16.02.2013 - 14:28

[stok1946]

Мне представляется более удачным рассматривать (а) ДИ для регрессии и (б) ДИ для прогноза (так короче и часто в литературе речь и идёт о prediction). Учитывая что программисты чётко различают между Confidence Interval и Prediction Interval,
(б) не знаю, формулы везде разные; чтобы привести их к единообразнию нужно поработать.

Большое спасибо за разъяснение. Я никогда толком не понимал, что такое ?предиктор?, но стеснялся спросить. И был категорически не прав :
Predictor or predictor variable - In a regression analysis, the variable that serves as the basis for prediction is called the predictor variable. It is also called the independent or explanatory variable. It is always plotted as the x variable in a scatter diagram. ( Предикторная переменная ? в регрессионном анализе независимая переменная, появляющиеся на правой стороне уравнений. Соответствует переменной x на диаграммах рассеяния).
Predicted variable - In regression analysis, the variable that is being regressed is called the dependent or predicted variable. It is always plotted as the y variable in a scatter diagram (Предсказываемая переменная - в регрессионном анализе зависимая переменная или отклик. Соответствует переменной y на диаграммах рассеяния).
Prediction interval - In a regression analysis, a confidence interval within which a future observation of the dependent variable, for a given value of the independent variable,lies with a given probability. Compare confidence band. (Определяет доверительные границы, в которых будут находиться с заданной вероятностью будущие значения прогноза зависимой переменной для данного значения независимой переменной).
Но почему тогда термин "предикторные" интервалы относят не к предикторам, а к зависимой переменной?

б) В каждой формуле всегда изначально заложена ошибка. Нет ничего практичнее хорошей теории. Для простой линейной регрессии формулы для рrediction interval и confidence band известны. Но я не нашел таковых в случае, например, полинома.
Тогда как, исходя из общих теоретических соображений оценить рrediction interval, если известны confidence band и standard error of predicted means на исходной выборке?

Сообщение отредактировал stok1946 - 6.03.2013 - 13:11

[Larina Tatjana]

Насколько я помню из теории, различают два вида доверительных областей:
а) Доверительная область для линии регрессии RD(x), (т.е. точнее, для предикторов модели)
б) Доверительная область для самих значений зависимой переменной YD(x)

И тут такие вопросы:
а) когда в STATISTICA и др. прикладных программах считают доверительные интервалы, то что имеют в виду RD(x), YD(x) или SD?
Например, Гайдышев в Approximations (APX) AtteStat выводит SD и называет это "Доверительные интервалы оценок модели" (но ведь "оценки" - это как раз и есть предикторы модели?);

"оценки" - это не только точечные оценки предикторов модели, а все совокупности оценок. Т.е. все комбинации, включая как точечные "оценки", так и доверительные интервалы.

в) почему на картинке к
http://ru.wikipedia.org/wiki/%D0%A0%D0%B5%...%D0%BA%D0%B0%29
95%-е доверительные области для линии регрессии RD(x) показаны криволинейными, и для значений YD(x) - прямыми красненькими линиями, хотя, согласно там же приведенных формул, YD(x) включает RD(x) .
Или мне это мерещится?
Нет, не мерещится. Так и должно быть. Криволинейные области получаются как результат одновременного изменения величина свободного члена (колебания линии регрессии вверх-вниз), так и значений коэффициентов регрессии (угол наклона линии регрессии). В итоге и получаются криволинейные области.

Рекомендую почитать хорошие книги по регрессии, которых немало.

[stok1946]

Рекомендую почитать хорошие книги по регрессии, которых немало.

Буду очень признателен, если Вы мне порекомендуете конкретный источник, где приводилась бы схема расчета предикторных интервалов для полиномиальной регрессии (и вообще, для общего нелинейного случая). Для простой линейной регрессии формулы есть на каждом шагу...
А точнее, вот код R:

CODE

# Определение векторов с активностями ферментов
x <- c(19.1, 22.4, 21.7, 20.2, 25.1, 18.6, 15.6, 17.6, 13.5, 22.4, 14.2, 13.3, 17.9,
24.6, 25.8, 16.1, 16.6, 31.2, 15.8, 11.4, 20.3, 6.2, 8.9, 9.4, 6.4, 9.21, 2.6, 8.2)
y <- c(24.5, 19.6, 22.2, 24.3, 25.5, 25.5, 25.2, 23.9, 26.1, 25.3, 21.9, 20.6, 27.1,
43.2, 25.6, 24.9, 25.8, 23.8, 23.5, 21.1, 23.8, 0.8, 15.5, 2.6, 2.8, 0.9, 5.7, 1.4)
xy <- data.frame(x,y)
model1 <- lm(y ~ x) # Линейная
model2 <- lm( y ~ x + I(x*x)) # Полином 2 степени
sr <- 1.96*summary(model2)$sigma ; x.order <- order(x)
par(mfrow = c(1,2))
plot(x, y, pch=16, main="95% CI и PI линейной модели")
lines(x, model1$fitted, lwd=2)
model1.pred <- predict(model1, level=0.95, interval="conf")
lines(x[x.order],model1.pred[,2][x.order], col="red")
lines(x[x.order],model1.pred[,3][x.order], col="red")
model1.pred2 <- predict(model1, level=0.95, interval="pred")
lines(x[x.order],model1.pred2[,2][x.order], col="blue")
lines(x[x.order],model1.pred2[,3][x.order], col="blue")
plot(x, y, pch=16, main="95% CI и PI полиномиальной модели")
lines(x[x.order], model2$fitted[x.order], lwd=2)
model2.pred <- predict(model2, level=0.95, interval="conf") # Доверительные интервалы
lines(x[x.order],model2.pred[,2][x.order], col="red")
lines(x[x.order],model2.pred[,3][x.order], col="red")
model2.pred2 <- predict(model2, level=0.95, interval="pred") # Предикторные интервалы
lines(x[x.order],model2.pred2[,2][x.order], col="blue")
lines(x[x.order],model2.pred2[,3][x.order], col="blue")
lines(x[x.order],model2$fitted[x.order]+sr, col="green") # Версия Гайдышева
lines(x[x.order],model2$fitted[x.order]-sr, col="green")

Получаем картинку с ДИ и ПИ для линейной и полиномиальной модели. Файл с картинкой приложил к сообщению


Зелененьким показаны "доверительные интервалы оценки модели" по Гайдышеву, постоянные на всем диапазоне Х-в, которые находятся как-то посередине между ПИ и ДИ.
А ведь есть еще и толерантные интервалы регрессии - см. стр. 665. - Прикладная математическая статистика. Кобзарь А.И.
А правда где, брат?

Сообщение отредактировал stok1946 - 6.03.2013 - 13:11

[nokh]

Проблема с нелинейной регрессией оказалась куда сложнее, чем я думал. Пусть не сразу, но формулы и для доверительного, и для прогнозного интервалов я нашёл. Но дело осложняется тем, что в эти формулы входит гессиан используемой нелинейной функции, а почти нигде не написано, как его рассчитать. Точнее написано, но с использованием матричной формы записи + насколько я понял этот расчёт проводится итерационно. Также найденные формулы - это аппроксимации, а в ряде источников убедительно показано, что они могут быть весьма далеки от реальности. Вот ссылка на хорошую книгу по нелинейной регрессии в биологии: http://yadi.sk/d/esAZDJ2b2rXZe . Хотя здесь не показано построение ДИ для всей регрессии, но зато показано построение ДИ для параметров нелинейной функции аж тремя способами, приводящими к разным значениям.

[100$]

Проблема с нелинейной регрессией оказалась куда сложнее, чем я думал. Пусть не сразу, но формулы и для доверительного, и для прогнозного интервалов я нашёл. Но дело осложняется тем, что в эти формулы входит гессиан используемой нелинейной функции, а почти нигде не написано, как его рассчитать. Точнее написано, но с использованием матричной формы записи + насколько я понял этот расчёт проводится итерационно. Также найденные формулы - это аппроксимации, а в ряде источников убедительно показано, что они могут быть весьма далеки от реальности. Вот ссылка на хорошую книгу по нелинейной регрессии в биологии: http://yadi.sk/d/esAZDJ2b2rXZe . Хотя здесь не показано построение ДИ для всей регрессии, но зато показано построение ДИ для параметров нелинейной функции аж тремя способами, приводящими к разным значениям.

Но дело осложняется тем, что в эти формулы входит гессиан используемой нелинейной функции, а почти нигде не написано, как его рассчитать.

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

Хотя здесь не показано построение ДИ для всей регрессии, но зато показано построение ДИ для параметров нелинейной функции аж тремя способами, приводящими к разным значениям.

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

P.S. Кстати, о полиномах. Полиномы хорошо подгоняют данные (н-р, сплайны, представляющие собой, чаще всего, кубические полиномы) на основе широко известного факта, что через n точек проходит полином n-1 порядка. Но вот прогнозные войства полиномов - отвратительные.
Кроме того, полиномы являются линейной (по параметрам) моделью, в связи с чем ее оценивание вообще не представляет сложностей.
Наконец, многое нелинейные модели после логарифмирования становятся линейными и оцениваются на ура.

Сообщение отредактировал 100$ - 25.02.2013 - 18:48

[nokh]

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

Для не слишком сложной функции я смогу найти частные производные и составить систему уравнений для оценки параметров методом наименьших квадратов (МНК), но и это - самообразование. А, к сожалению, в векторной алгебре не понимаю даже на уровне троечника Буду очень признателен если Вы объясните чем заменить гессиан, если мы будем действовать в рамках именно МНК.

Наконец, многое нелинейные модели после логарифмирования становятся линейными и оцениваются на ура.

С этим есть известные проблемы смещения оценок. Всё зависит от того что в действительности должно быть распределено нормально: ошибка исходной нелинейной функции или ошибка линеаризованной преобразованием функции. Допустим, что нормально распределена ошибка именно исходной функции. Тогда если подгонять модель путём предварительного преобразование к линейному виду, то после ретрансформации ошибка будет распределена ненормально: например, в случае предварительного логарифмирования - логнормально; при этом оценки самих параметров окажутся немного смещены. Поэтому если раньше (до прихода ЭВМ) "на ура" действовали через преобразования шкал, то сейчас обычная практика - подгонка исходных нелинейных моделей по итерационным алгоритмам (хотя я думаю, что для ряда зависимостей естественней было бы работать "по старинке", исходя именно из природы данных). Если бы не было этой проблемы смещения, то через преобразование оси х по Боксу-Коксу можно было бы оптимально линеаризировать огромную массу нелинейных зависимостей самого распространённого степенного семейства.

Сообщение отредактировал nokh - 26.02.2013 - 22:31

[100$]

Для не слишком сложной функции я смогу найти частные производные и составить систему уравнений для оценки параметров методом наименьших квадратов (МНК), но и это - самообразование. А, к сожалению, в векторной алгебре не понимаю даже на уровне троечника Буду очень признателен если Вы объясните чем заменить гессиан, если мы будем действовать в рамках именно МНК.

Если удалось составить систему нормальных уравнений МНК (т.е. найдя первые производные по параметрам, приравнять их к 0), то она (теоретически) решается с помощью обратной матрицы. Однако, в общем случае проблема нелинейного МНК не решена, и оный МНК сводится к методу макс. правдоподобия. (что забавно: в линейном случае-метод МП сводится к МНК, в нелинейном - МНК ->МП).

А аппроксимации (обратного) гессиана для каждого метода оптимизации свои: у Левенберга-Марквардта(L-M)-одна, у Дэвидона-Флетчера-Пауэлла (DFP)-другая, у Бройдена-Флетчера-Гольдфарба-Шенно-третья(BFGS), у Берндта-Холла-Холла-Хаусмана(BHHH)- четвертая. Кто во что горазд.

Сообщение отредактировал 100$ - 27.02.2013 - 11:04

[nokh]

Если удалось составить систему нормальных уравнений МНК (т.е. найдя первые производные по параметрам, приравнять их к 0), то она (теоретически) решается с помощью обратной матрицы. Однако, в общем случае проблема нелинейного МНК не решена, и оный МНК сводится к методу макс. правдоподобия. (что забавно: в линейном случае-метод МП сводится к МНК, в нелинейном - МНК ->МП).

А аппроксимации (обратного) гессиана для каждого метода оптимизации свои: у Левенберга-Марквардта(L-M)-одна, у Дэвидона-Флетчера-Пауэлла (DFP)-другая, у Бройдена-Флетчера-Гольдфарба-Шенно-третья(BFGS), у Берндта-Холла-Холла-Хаусмана(BHHH)- четвертая. Кто во что горазд.

Как всё непросто . "Добить" бы тему ветки хоть для случая параболы - я и для параболы формул ДИ не встречал!

Сообщение отредактировал nokh - 11.03.2013 - 22:06

[100$]

Как всё непросто . "Добить" бы тему ветки хоть для случая параболы - я и для параболы формул ДИ не встречал!

Я тоже. Возможно, по причине того, что полиномы по оцениваемым параметрам линейны, для них сгодятся те же ф-лы, что и для парной линейной регрессии.

[stok1946]

Мне все же кажется, что время приближенных (т.е. основанных на предположениях о распределениях) формул катастрофически уходит.
Нет никаких проблем посчитать доверительные интервалы для полиномиальной регрессии, сплайнов и проч. бутстреп-методом. Нужно только понять, что хочется считать.
Если многократно извлекать из генеральной совокупности наблюдений различные выборки из n пар (x, y) значений и строить по ним модели регрессии, то за пределы "доверительной полосы" может выйти только alfa% таких линий (alfa = 5%).
Из этих соображений легко выполнить расчет CI для регрессии:
а) выделяем m = 100 или больше опорных значений, равномерно распределенных по шкале x, относительно которых будет рассчитываться величины доверительных интервалов;
б) делаем случайную выборку с возвращениями из порядковых номеров строк исходной таблицы и по этой перевыборке рассчитываем нужную нам модель регрессии любой сложности;
в) по модели п. б) выполняем расчет m опорных значений зависимой переменной;
г) пункты б-в) повторяем В = 1000 или больше раз, после чего для каждой 100 опорных точек x вычисляем по 1000 расчетных значений ŷ , т.е. воспроизводим распределение прогноза отклика в этих точках;
д) для каждой из опорных точек x находим значения квантилей при p = 1-alfa/2 и p = alfa/2 и вычисляем доверительные интервалы по любой из 3-4 формул (процентили, Халла-Луннеборга, стьюдентизированного типа, ВСа).
Несложные скриптики приведены в "Рандомизация и бутстреп: статистический анализ в биологии и экологии с использованием R". - разделы 3.4 и 7.1 и с картинками.

Но вот что нужно ловить в в море бутстрепа, чтобы посчитать Prediction Intervals? Казалось бы нужно что-то добавлять/отнимать к ŷ (или к рассчитанным CI), но что именно - ничего на ум не приходит.
Это потому, что я никак не пойму смысла PI, исходя из обычной статитической прелюдии: "Если многократно извлекать из генеральной совокупности наблюдений различные выборки из n пар (x, y) значений и строить по ним модели регрессии, то..."

[p2004r]

Но вот что нужно ловить в в море бутстрепа, чтобы посчитать Prediction Intervals? Казалось бы нужно что-то добавлять/отнимать к ŷ (или к рассчитанным CI), но что именно - ничего на ум не приходит.
Это потому, что я никак не пойму смысла PI, исходя из обычной статитической прелюдии: "Если многократно извлекать из генеральной совокупности наблюдений различные выборки из n пар (x, y) значений и строить по ним модели регрессии, то..."

Полностью согласен с высказыванием о методах ручного докомпьютерного счета. Все эти подходы в настоящее время имеют значение лишь для развития теории. На практике проще и надежнее проводить прямые измерения.

Что касается основного вопроса:

In statistical inference, specifically predictive inference, a prediction interval is an estimate of an interval in which future observations will fall, with a certain probability, given what has already been observed. Prediction intervals are often used in regression analysis.

Prediction intervals are used in both frequentist statistics and Bayesian statistics: a prediction interval bears the same relationship to a future observation that a frequentist confidence interval or Bayesian credible interval bears to an unobservable population parameter: prediction intervals predict the distribution of individual future points, whereas confidence intervals and credible intervals of parameters predict the distribution of estimates of the true population mean or other quantity of interest that cannot be observed. Prediction intervals are also present in forecasts. Some experts have shown that it is difficult to estimate the prediction intervals of forecasts that have contrary series

как я понимаю ключевое отличие это prediction intervals predict the distribution of individual future points, whereas confidence intervals and credible intervals of parameters predict the distribution of estimates of the true population mean or other quantity of interest that cannot be observed

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

Все это для регрессии сводится к определению остатков и изучение их распределения. Сводящееся в свою очередь к построению доверительного интервала который плюсминусуется от прогнозного значения модели на данном наборе значений независимых переменных (естественно это все когда модель хорошая и остатки независимы от параметров регрессии)

Поскольку остатки могут обладать весьма специфическим распределением, то для бутстрепа по моему открывается масса возможностей в построении интервалов. Особенно когда остатки "плохие". В этом случае надо построить интервал который зависит от состояния независимых переменных регрессионной модели.

Конечно надо подумать насколько такая процедура геометрически эквивалентна прямому построению доверительного интервала бутстрепом. Очень все какую то инкридибел машину напоминает.

[Енот]

Задам вопрос в этой теме, чтобы не открывать новую.
Насколько я знаю из руководств по стат анализу, в многофакторных регрессионных моделях рассчитывается OR для каждого из предикторов. Сейчас работаю с иностранной статьей, в которой указано значение OR и доверительный интервал только одно для всей модели (то есть для 3х предикторов сразу). Как такое может быть? Для анализа иностранныи авторами использовалась программа SPSS IBM 20. Мне нужно провести анализ аналогичным образом, опираясь на алгоритм этой статьи.

Сообщение отредактировал Енот - 1.07.2013 - 15:46

[TheThing]

Задам вопрос в этой теме, чтобы не открывать новую.
Насколько я знаю из руководств по стат анализу, в многофакторных регрессионных моделях рассчитывается OR для каждого из предикторов. Сейчас работаю с иностранной статьей, в которой указано значение OR и доверительный интервал только одно для всей модели (то есть для 3х предикторов сразу). Как такое может быть? Для анализа иностранныи авторами использовалась программа SPSS IBM 20. Мне нужно провести анализ аналогичным образом, опираясь на алгоритм этой статьи.

А можно ссылку на статью?