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

15

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
7
+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ć p od logit, który jest logarytmem logarytmicznym ( p / ( 1 - p ) ) . pglmplog(p/(1p))
whuber

Odpowiedzi:

19

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 p , chcesz, aby wartość wyjściowa wynosiła q = logit ( p ) = log ( p / ( 1 - p ) ) . Model statystyczny jest terazβ0=0pq=logit(p)=log(p/(1p))

YBinom(μ)μ=1/(1+exp(η))η=β0+q

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
2
(+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
1
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
2
może anova(.,test="Chisq")?
Ben Bolker
6

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 log ( o d d s ) = log pp=0.5. Testując hipotezę, żep=0,5log(odds)=logp1p=log1=0p=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
1
p<0.05
2
confintp<0,05
2

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.

  1. 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
  2. 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")

    CDF odrzucenia przez alfa

    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
2
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,
2
(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
1
@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