anova.rq() в пакете quantreg в R

Мне интересно сравнить оценки из разных квантилей (один и тот же результат, одни и те же ковариаты) с использованием функции anova.rqlist, вызываемой anova в среде пакета quantreg в R. Однако математика в функции выходит за рамки моего элементарного опыта. Допустим, я подобрал 3 модели в разных квантилях;

library(quantreg)
data(Mammals) # data in quantreg to be used as a useful example
fit1 <- rq(weight ~ speed + hoppers + specials, tau = .25, data = Mammals)
fit2 <- rq(weight ~ speed + hoppers + specials, tau = .5, data = Mammals)
fit3 <- rq(weight ~ speed + hoppers + specials, tau = .75, data = Mammals)

Затем я сравниваю их, используя;

anova(fit1, fit2, fit3, test="Wald", joint=FALSE)

Мой вопрос: какая из этих моделей используется в качестве основы для сравнения?

Мое понимание теста Вальда (статья в вики)

введите здесь описание изображения

где θ^ — оценка интересующего(их) параметра(ов) θ, которая сравнивается с предлагаемым значением θ0.

Итак, мой вопрос: что функция anova в quantreg выбирает в качестве θ0?

Основываясь на pvalue, возвращенном из anova, я могу предположить, что он выбирает самый низкий указанный квантиль (т.е. tau=0.25). Есть ли способ указать медиану (tau = 0.5) или, что еще лучше, среднюю оценку, полученную с использованием lm(y ~ x1 + x2 + x3, data)?

anova(fit1, fit2, fit3, joint=FALSE)

на самом деле производит

Quantile Regression Analysis of Deviance Table

Model: weight ~ speed + hoppers + specials
Tests of Equality of Distinct Slopes: tau in {  0.25 0.5 0.75  }

             Df Resid Df F value  Pr(>F)  
speed         2      319  1.0379 0.35539  
hoppersTRUE   2      319  4.4161 0.01283 *
specialsTRUE  2      319  1.7290 0.17911  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

пока

anova(fit3, fit1, fit2, joint=FALSE)

дает точно такой же результат

Quantile Regression Analysis of Deviance Table

Model: weight ~ speed + hoppers + specials
Tests of Equality of Distinct Slopes: tau in {  0.5 0.25 0.75  }

             Df Resid Df F value  Pr(>F)  
speed         2      319  1.0379 0.35539  
hoppersTRUE   2      319  4.4161 0.01283 *
specialsTRUE  2      319  1.7290 0.17911  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Порядок моделей явно меняется в анова, но как получается, что значение F и Pr(>F) одинаковы в обоих тестах?


person JustGettinStarted    schedule 23.09.2015    source источник
comment
Обратите внимание, что этот вопрос частично мотивирован тем, что я прочитал книгу (пусть) «Квантильная регрессия» Линсинь Хао и Дэниела К. Наймана. На страницах 49, 50 и 60 они подробно описывают, как бутстрапированные квантильные оценки можно сравнивать по квантилям с использованием критерия Вальда. Они четко описывают, как интересующий квантиль можно сравнить с медианой, наименьшим или ближайшим квантилем, что заставило меня задуматься о том, что делает функция анова в пакете quantreg для оценок без начальной загрузки и какой квантиль выбран в качестве эталона для таких сравнений. .   -  person JustGettinStarted    schedule 23.09.2015


Ответы (1)


Используются все введенные вами квантили, и в качестве эталона не используется ни одна модель.

Я предлагаю вам прочитать этот пост и связанный с ним ответ, чтобы понять, какова ваша "тэта-регрессия". .0" есть.

Я считаю, что вы пытаетесь проверить, параллельны ли линии регрессии. Другими словами, является ли влияние переменных-предикторов (здесь только доход) одинаковым для всех квантилей.

Вы можете использовать anova() из пакета quantreg, чтобы ответить на этот вопрос. Вы действительно должны использовать несколько подгонок для каждого квантиля.

Когда вы используете joint=FALSE, как вы это делали, вы получаете сравнения с коэффициентами. Но у вас есть только один коэффициент, поэтому есть только одна линия! И ваши результаты говорят о том, что влияние дохода неодинаково по квантилям в вашем примере. Используйте несколько переменных-предикторов, и вы получите несколько p-значений.

Вы можете выполнить общий тест на равенство всех наборов коэффициентов, если не используете joint=FALSE, и это даст вам «Совместный тест на равенство наклонов» и, следовательно, только одно p-значение.

РЕДАКТИРОВАТЬ:

Я думаю, что theta.0 - это средний наклон для всех значений «тау» или фактическая оценка из «lm()», а не конкретный наклон любой из моделей. Я полагаю, что 'anova.rq()' не требует какого-либо конкретного низкого значения «тау» или даже медианы «тау».

Есть несколько способов проверить это. Либо делайте расчеты вручную, приравнивая тета.0 к среднему значению, либо сравнивайте множество комбинаций, потому что тогда вы можете столкнуться с ситуацией, когда некоторые из ваших моделей близки к модели с низкими значениями «тау», но не к «lm». ()' ценность. Таким образом, если theta.0 - это наклон первой модели с самым низким значением тау, то ваш Pr(>F) будет высоким, тогда как в другом случае он будет низким.

Возможно, этот вопрос следовало задать на перекрестной проверке.

person Julian Wittische    schedule 23.09.2015
comment
1) Сообщение, на которое вы ссылаетесь, было информативным: в нем говорится, что наш интересующий параметр обозначается θ0, и обычно это 0, поскольку мы хотим проверить, отличается ли коэффициент от 0 или нет. Меня не интересуют коэффициенты, отличные от 0 (значение p в summary.rq(fit1) скажет мне об этом). Я пытаюсь проверить, значительно ли влияние предикторов при tau = 0,1 и tau = 0,9 отличается от tau = 0,5. В базе anova первая модель, указанная в аргументе, используется как θ0, почему это не так с anova.rq и как сделать так? - person JustGettinStarted; 24.09.2015
comment
2) скорректированный пример для включения нескольких ковариатов, чтобы можно было применить joint=FALSE. Для моего приложения меня интересует, как эффекты одного конкретного предиктора варьируются в разных квантилях. Даже включение модели с несколькими ковариатами, изменение порядка моделей в аргументах для anova.rq не меняет результат. - person JustGettinStarted; 24.09.2015
comment
Вы уже можете сказать, что эффекты предиктора различны для разных значений тау. Если вам нужна информация о конкретных значениях «тау», вы всегда можете использовать оригинальную функцию «anova()», я думаю. - person Julian Wittische; 25.09.2015
comment
В ответ на ваше редактирование. На самом деле я подозреваю, что anova.rq использует все возможные комбинации θ.hat и θ0, чтобы проверить, является ли θhat - θ0 = 0 и НЕ, если θhat = θ0. Таким образом, эталона (θ0) как такового нет, потому что каждый тау-белок тестируется против любого другого тау во всех комбинациях. Я предполагаю, что это из книги Koenker Quantile Regression 2005 pg 75-76, а также из кодирования в anova.rq, показанного getAnywhere(anova.rq). Это объясняет, почему порядок моделей ничего не меняет, а lm не вызывается anova.rq. - person JustGettinStarted; 25.09.2015