Использование glm () вместо простого теста хи-квадрат

15

Я заинтересован в изменении нулевых гипотез, используя glm()в R.

Например:

x = rbinom(100, 1, .7)  
summary(glm(x ~ 1, family = "binomial"))

проверяет гипотезу, что . Что если я захочу изменить значение null на = какое-то произвольное значение внутри ? рпзнак равно0,5пglm()

Я знаю, что это можно сделать также с помощью prop.test()и chisq.test(), но я хотел бы изучить идею использования glm()для проверки всех гипотез, относящихся к категориальным данным.

Билл Равенвуд
источник
7
+1. очевидно, относится к биномиальному параметру, выраженному как вероятность. Поскольку естественная ссылка (и та, которая используется по умолчанию) - это logit, во избежание путаницы важно отличать p от его logit, который представляет собой log odds log ( p / ( 1 - p ) ) . пglmпжурнал(п/(1-п))
whuber

Ответы:

19

Вы можете использовать смещение : glmс family="binomial"оценочными параметрами по лог-коэффициентам или по шкале логитов, так что соответствует лог-коэффициентам 0 или вероятности 0,5. Если вы хотите сравнить с вероятностью p , вы хотите, чтобы базовое значение было равно q = logit ( p ) = log ( p / ( 1 - p ) ) . Статистическая модель сейчасβ0знак равно0пQзнак равнологит(п)знак равножурнал(п/(1-п))

Y~Бином(μ)μзнак равно1/(1+ехр(-η))ηзнак равноβ0+Q

где только последняя строка изменилась из стандартной настройки. В коде R:

  • использовать offset(q)в формуле
  • функция logit / log-odds qlogis(p)
  • Немного досадно, но вы должны указать значение смещения для каждого элемента в переменной ответа - R не будет автоматически копировать для вас постоянное значение. Это делается ниже путем настройки фрейма данных, но вы можете просто использовать rep(q,100).
x = rbinom(100, 1, .7)
dd <- data.frame(x, q = qlogis(0.7)) 
summary(glm(x ~ 1 + offset(q), data=dd, family = "binomial"))
Бен Болкер
источник
2
(+1) это даст вам тест Вальда. LRT может быть сделан, приспосабливая нулевую модель glm(y ~ offset(q)-1, family=binomial, data=dd)и используя lrtestиз lmtestпакета. Критерий хи-квадрат Пирсона является тестом для модели GLM. Wald / LRT / Score - это последовательные тесты, которые должны обеспечивать эквивалентный вывод при достаточно больших размерах выборки.
AdamO
1
Я думаю, что вы также можете использовать anova()базу R на glm, чтобы пройти тест LR
Бен Болкер,
Интересно, я потерял привычку использовать ANOVA. Тем не менее, я наблюдаю, как anova отказывается печатать значение для теста, тогда как lrtestделает.
AdamO
2
может быть anova(.,test="Chisq")?
Бен Болкер
6

Посмотрите на доверительный интервал для параметров вашего GLM:

> set.seed(1)
> x = rbinom(100, 1, .7)
> model<-glm(x ~ 1, family = "binomial")
> confint(model)
Waiting for profiling to be done...
    2.5 %    97.5 % 
0.3426412 1.1862042 

Это доверительный интервал для лог-шансов.

Для имеем log ( o d d s ) = log pпзнак равно0,5журнал(оdds)знак равножурналп1-пзнак равножурнал1знак равно0пзнак равно0,5

п

Лукаш Дерило
источник
1
п<0,05
2
confintп<0,05
2

Не является (полностью) правильным / точным использование p-значений на основе z- / t-значений в функции glm.summary в качестве проверки гипотезы.

  1. Это запутанный язык. Сообщаемые значения называются z-значениями. Но в этом случае они используют оценочную стандартную ошибку вместо истинного отклонения. Поэтому в действительности они ближе к т-значениям . Сравните следующие три выхода:
    1) summary.glm
    2) t-тест
    3) z-тест

    > set.seed(1)
    > x = rbinom(100, 1, .7)
    
    > coef1 <- summary(glm(x ~ 1, offset=rep(qlogis(0.7),length(x)), family = "binomial"))$coefficients
    > coef2 <- summary(glm(x ~ 1, family = "binomial"))$coefficients
    
    > coef1[4]  # output from summary.glm
    [1] 0.6626359
    > 2*pt(-abs((qlogis(0.7)-coef2[1])/coef2[2]),99,ncp=0) # manual t-test
    [1] 0.6635858
    > 2*pnorm(-abs((qlogis(0.7)-coef2[1])/coef2[2]),0,1) # manual z-test
    [1] 0.6626359
  2. Они не являются точными значениями р. Точное вычисление p-значения с использованием биномиального распределения будет работать лучше (с вычислительной мощностью в настоящее время это не проблема). T-распределение, предполагающее гауссово распределение ошибки, не является точным (оно завышает p, превышение уровня альфа встречается реже в «реальности»). Смотрите следующее сравнение:

    # trying all 100 possible outcomes if the true value is p=0.7
    px <- dbinom(0:100,100,0.7)
    p_model = rep(0,101)
    for (i in 0:100) {
      xi = c(rep(1,i),rep(0,100-i))
      model = glm(xi ~ 1, offset=rep(qlogis(0.7),100), family="binomial")
      p_model[i+1] = 1-summary(model)$coefficients[4]
    }
    
    
    # plotting cumulative distribution of outcomes
    outcomes <- p_model[order(p_model)]
    cdf <- cumsum(px[order(p_model)])
    plot(1-outcomes,1-cdf, 
         ylab="cumulative probability", 
         xlab= "calculated glm p-value",
         xlim=c(10^-4,1),ylim=c(10^-4,1),col=2,cex=0.5,log="xy")
    lines(c(0.00001,1),c(0.00001,1))
    for (i in 1:100) {
      lines(1-c(outcomes[i],outcomes[i+1]),1-c(cdf[i+1],cdf[i+1]),col=2)
    #  lines(1-c(outcomes[i],outcomes[i]),1-c(cdf[i],cdf[i+1]),col=2)
    }
    
    title("probability for rejection as function of set alpha level")

    CDF отклонения альфа

    Черная кривая представляет равенство. Красная кривая ниже этого. Это означает, что для данного вычисленного значения p с помощью функции суммирования glm мы находим эту ситуацию (или большую разницу) реже, чем указывает значение p.

Секст Эмпирик
источник
Хм .. Я могу быть смущен обоснованием использования T-распределения для GLM. Можете ли вы взять пик по связанному вопросу, который я только что задал здесь ?
AdamO
2
Этот ответ интересный, но проблемный. (1) ОП на самом деле не спрашивал о разнице между оценкой, хи-квадратом, «точным» или основанным на GLM подходом к проверке гипотез о биномиальных ответах (они могли бы уже все это знать), так что это не t ответить на вопрос, который был задан; (2) оценки остаточной дисперсии и т. Д. Имеют другой набор допущений и распределений выборки из линейных моделей (как в вопросе @ AdamO), поэтому использование t-критерия является спорным; ...
Бен Болкер,
2
(3) «точные» доверительные интервалы для биномиальных ответов на самом деле хитры («точные» интервалы [Клоппера-Уилсона] консервативны; тесты с
оценками
@Ben Вы правы, что z-тест на самом деле лучше, чем t-тест. График, отображаемый в ответе, предназначен для z-теста. Он использует вывод функции GLM. Суть моего ответа заключалась в том, что «p-значение» - хитрая вещь. Поэтому я считаю, что лучше вычислить его явно, например, используя нормальное распределение, а не извлекать p-значение из функции glm, которая очень удобно сдвигается со смещением, но скрывает начало вычислений для p-значения. ,
Секст Эмпирик
1
@BenBolker, я считаю, что точный тест действительно консервативный, но ... только потому, что на самом деле мы не выбираем из идеальных биномиальных распределений. Альтернативный z-тест, только лучше с эмпирической точки зрения. Дело в том, что две «ошибки» взаимно компенсируют 1) биномиальное распределение, не являющееся реальным распределением невязок в практических ситуациях, 2) z-распределение не является точным выражением для биномиального распределения. Сомнительно, следует ли нам выбирать неправильное распределение для неправильной модели, просто потому, что на практике получается «хорошо».
Секст Эмпирик