Użycie glm () jako zamiennika prostego testu chi kwadrat

Interesuje mnie zmiana hipotez zerowych za pomocą glm()R.

Na przykład:

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

sprawdza hipotezę, że $p = 0.5$ . Co jeśli chcę zmienić wartość null na $p$ = jakąś dowolną wartość, w obrębie glm()?

Wiem, że można to zrobić również za pomocą prop.test()i chisq.test(), ale chciałbym zbadać pomysł użycia glm()do testowania wszystkich hipotez dotyczących danych kategorycznych.

— Bill Ravenwood
źródło

+1.

ewidentnie odnosi się do parametru dwumianowego wyrażonego jako prawdopodobieństwo. Ponieważ linkiem naturalnym (i domyślnie używanym ) jest logit, aby uniknąć nieporozumień, ważne jest, aby odróżnić

od logit, który jest logarytmem

p

$p$ glm

p

$p$

\log (p / (1 - p))

$\log(p/(1-p))$

— whuber

Odpowiedzi:

Możesz użyć przesunięcia : glmz family="binomial"parametrami oszacowań na log-odds lub skali logit, więc odpowiada log-odds 0 lub prawdopodobieństwu 0,5. Jeśli chcesz porównać z prawdopodobieństwem , chcesz, aby wartość wyjściowa wynosiła . Model statystyczny jest teraz $\beta_0=0$ $p$ $q = \textrm{logit}(p)=\log(p/(1-p))$

\begin{aligned} Y & \sim Binom (μ) \\ μ & = 1 / (1 + \exp (- η)) \\ η & = β_{0} + q \end{aligned}

$\begin{split} Y & \sim \textrm{Binom}(\mu) \\ \mu & =1/(1+\exp(-\eta)) \\ \eta & = \beta_0 + q \end{split}$

gdzie tylko ostatnia linia zmieniła się od standardowej konfiguracji. W kodzie R:

użyj offset(q)w formule
funkcja logit / log-odds to qlogis(p)
nieco irytujące, musisz podać wartość przesunięcia dla każdego elementu w zmiennej odpowiedzi - R nie będzie automatycznie replikować stałej wartości dla ciebie. Odbywa się to poniżej poprzez ustawienie ramki danych, ale możesz po prostu użyć 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"))

— Ben Bolker
źródło

(+1) da ci to test Walda. LRT można wykonać dopasowując model zerowy glm(y ~ offset(q)-1, family=binomial, data=dd)i używając lrtestz lmtestpakietu. Test chi-kwadrat Pearsona jest testem punktowym dla modelu GLM. Wald / LRT / Score są spójnymi testami i powinny zapewniać równoważne wnioskowanie w stosunkowo dużych próbach.

— AdamO

Myślę, że możesz także użyć anova()z bazy R na glm, aby przejść test LR

— Ben Bolker

Co ciekawe, straciłem nawyk używania ANOVA. Jednak obserwuję, że anova odmawia wydrukowania wartości testu, podczas gdy lrtestrobi to.

— AdamO

może anova(.,test="Chisq")?

— Ben Bolker

Spójrz na przedział ufności dla parametrów twojego 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

Jest to przedział ufności dla logarytmicznych szans.

Dla mamy $p=0.5$ . Testując hipotezę, że $\log(odds) = \log \frac{p}{1-p} = \log 1 = 0$ $p=0.5$ jest równoważne sprawdzeniu, czy przedział ufności zawiera 0. Ta nie ma, więc hipoteza jest odrzucana.

Teraz, dla dowolnego dowolnego , możesz obliczyć log-odds i sprawdzić, czy jest on w przedziale ufności. $p$

— Łukasz Deryło
źródło

p < 0.05

$p<0.05$

confint

p < 0, 05

$p<0,05$

Używanie wartości p opartych na wartościach z / t w funkcji glm.summary jako testu hipotezy nie jest (całkowicie) poprawne / dokładne.

To jest mylący język. Podane wartości są nazywane wartościami Z. Ale w tym przypadku używają szacowanego błędu standardowego zamiast prawdziwego odchylenia. Dlatego w rzeczywistości są one bliższe wartościom t . Porównaj następujące trzy dane wyjściowe:
1) Podsumowanie. Glm
2) Test t
3) Test 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

Nie są to dokładne wartości p. Dokładne obliczenie wartości p przy użyciu rozkładu dwumianowego działałoby lepiej (przy dzisiejszej mocy obliczeniowej nie jest to problemem). Rozkład t, zakładając rozkład błędu Gaussa, nie jest dokładny (przecenia p, przekroczenie poziomu alfa występuje rzadziej w „rzeczywistości”). Zobacz następujące porównanie:

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

Czarna krzywa przedstawia równość. Czerwona krzywa znajduje się poniżej. Oznacza to, że dla danej obliczonej wartości p za pomocą funkcji podsumowania glm znajdujemy tę sytuację (lub większą różnicę) w rzeczywistości rzadziej niż wskazuje wartość p.

— Sextus Empiricus
źródło

Hmm .. Może się mylę co do uzasadnienia zastosowania rozkładu T dla GLM. Czy możesz znieść szczyt w powiązanym pytaniu, które właśnie tutaj zadałem ?

— AdamO

Ta odpowiedź jest interesująca, ale problematyczna. (1) OP w rzeczywistości nie pytał o różnicę między podejściem punktowym, chi-kwadratowym, „dokładnym” lub podejściem opartym na GLM do testowania hipotez dotyczących odpowiedzi dwumianowych ( mogliby już znać wszystkie te rzeczy), więc to nie „ t odpowiedzieć na zadane pytanie; (2) szacunki wariancji rezydualnej itp. Mają inny zestaw założeń i rozkładów próbkowania z modeli liniowych (jak w pytaniu @ AdamO), więc stosowanie testu t jest dyskusyjne; ...

— Ben Bolker,

(3) „dokładne” przedziały ufności dla odpowiedzi dwumianowych są w rzeczywistości trudne (przedziały „dokładne” [Clopper-Wilson] są konserwatywne; testy punktowe mogą działać lepiej w niektórych zakresach

— Ben Bolker,

@Ben Masz rację, że test Z jest w rzeczywistości lepszy niż test T. Wykres wyświetlany w odpowiedzi dotyczy testu Z. Wykorzystuje dane wyjściowe funkcji GLM. Najważniejsze w mojej odpowiedzi było to, że „wartość p” jest trudna. Dlatego uważam, że lepiej jest to obliczyć jawnie, np. Używając rozkładu normalnego, zamiast wyodrębnić wartość p z funkcji glm, która bardzo wygodnie została przesunięta z przesunięciem, ale ukrywa początki obliczeń dla wartości p .

— Sextus Empiricus

@BenBolker, uważam, że dokładny test jest rzeczywiście konserwatywny, ale ... tylko dlatego, że w rzeczywistości nie pobieramy próbek z idealnych rozkładów dwumianowych. Alternatywny test Z jest lepszy tylko z empirycznego punktu widzenia. Chodzi o to, że dwa „błędy” wzajemnie się znoszą 1) rozkład dwumianowy nie będący rzeczywistym rozkładem reszt w sytuacjach praktycznych, 2) rozkład z nie będący dokładnym wyrażeniem dla rozkładu dwumianowego. Wątpliwe jest, czy powinniśmy preferować niewłaściwą dystrybucję dla niewłaściwego modelu, tylko dlatego, że w praktyce okazuje się „ok”.

— Sextus Empiricus