Format wejściowy dla odpowiedzi w dwumianowym glm w R.

W Ristnieją trzy sposoby, aby sformatować dane wejściowe dla regresji logistycznej z wykorzystaniem glmfunkcji:

Dane mogą być w formacie „binarnym” dla każdej obserwacji (np. Y = 0 lub 1 dla każdej obserwacji);
Dane mogą być w formacie „Wilkinson-Rogers” (np. y = cbind(success, failure)), Przy czym każdy wiersz reprezentuje jeden zabieg; lub
Dane mogą mieć format ważony dla każdej obserwacji (np. Y = 0,3, wagi = 10).

Wszystkie trzy podejścia dają te same oszacowania współczynników, ale różnią się stopniami swobody i wynikającymi z nich wartościami dewiacji i wynikami AIC. Dwie ostatnie metody mają mniej obserwacji (a zatem stopni swobody), ponieważ wykorzystują każde leczenie dla liczby obserwacji, podczas gdy pierwsze wykorzystują każdą obserwację dla liczby obserwacji.

Moje pytanie: Czy korzystanie z jednego formatu wejściowego nad innym ma przewagę liczbową lub statystyczną? Jedyną zaletą, jaką widzę, jest to, że nie muszę ponownie formatować danych, Raby użyć z modelem.

Przejrzałem dokumentację glm , przeszukałem Internet i tę stronę i znalazłem jeden stycznie powiązany post , ale nie ma wskazówek na ten temat.

Oto symulowany przykład, który demonstruje takie zachowanie:

# Write function to help simulate data
drc4 <- function(x, b =1.0, c = 0, d = 1, e = 0){
    (d - c)/ (1 + exp(-b * (log(x)  - log(e))))
}
# simulate long form of dataset
nReps = 20
dfLong <- data.frame(dose = rep(seq(0, 10, by = 2), each = nReps))
dfLong$mortality <-rbinom(n = dim(dfLong)[1], size = 1,
                              prob = drc4(dfLong$dose, b = 2, e = 5))

# aggregate to create short form of dataset
dfShort <- aggregate(dfLong$mortality, by = list(dfLong$dose), 
                     FUN = sum)
colnames(dfShort) <- c("dose", "mortality")
dfShort$survival <- nReps - dfShort$mortality 
dfShort$nReps <- nReps
dfShort$mortalityP <- dfShort$mortality / dfShort$nReps

fitShort <- glm( cbind(mortality, survival) ~ dose, 
                 data = dfShort, 
                 family = "binomial")
summary(fitShort)

fitShortP <- glm( mortalityP ~ dose, data = dfShort, 
                  weights = nReps,     
                  family = "binomial")
summary(fitShortP)

fitLong <- glm( mortality ~ dose, data = dfLong, 
                family = "binomial")
summary(fitLong)

r logistic generalized-linear-model Richard Erickson
źródło

W twoim przykładzie różnica między odchyleniem zerowym a odchyleniem resztkowym jest równa dla wszystkich trzech modeli. Jeśli dodasz lub usuniesz parametr, zmiana w AIC będzie również taka sama dla wszystkich trzech.

Jonny Lomond

Odpowiedziałeś sobie: jeśli korzystasz z tych samych danych, z tymi samymi wynikami, to jak mogą się różnić? Co więcej, jeśli metoda dałaby inne wyniki tylko z powodu dostarczenia danych w innym formacie, coś byłoby głęboko nie tak z nią lub z jej implementacją.

Tim

Format WR jest ostatecznie ważonym prawdopodobieństwem. Problem z wagami polega na tym, że R nie może stwierdzić, czy są to wagi częstotliwości, wagi prawdopodobieństwa lub inne. W przypadku ważenia ankietowego np. Możesz mieć tylko n obserwacji, ale reprezentują one segmenty populacji / ramki próbkowania. Tak więc stopnie swobody wynoszą w rzeczywistości 100. svyglmz pakietu ankiet daje lepsze metody radzenia sobie z argumentem wagi.

AdamO,

Ale jeśli pasowałbyś do quasibinomialnego modelu przy użyciu wszystkich trzech sposobów kodowania, przyniosłyby one różne wyniki, prawda, ponieważ można mieć dodatnią naddyspersję, gdy kodowane są jako dane dwumianowe, ale nie, gdy są kodowane jako dane logistyczne / binarne. A może się mylę?

Tom Wenseleers,

Odpowiedzi:

Nie ma statystycznego powodu, aby preferować jeden od drugiego, oprócz przejrzystości pojęciowej. Chociaż zgłaszane wartości odchyleń są różne, różnice te są całkowicie spowodowane nasyconym modelem. Zatem żadne porównanie z wykorzystaniem względnego odchylenia między modelami pozostaje nienaruszone, ponieważ prawdopodobieństwo dziennika nasyconego modelu zostaje anulowane.

Myślę, że warto przejść przez obliczenie wyraźnego odchylenia.

$i$ $n_i$ $i$ $p_i$ $i$ $y_{ij}$ $j$ $i$

Długa forma

\sum_{ja} \sum_{jot} (\log (p_{ja}) y_{ja jot} + \log (1 - p_{ja}) (1 - y_{ja jot}))

$\sum_i\sum_j\left(\log(p_i)y_{ij} + \log(1 - p_i)(1 - y_{ij})\right)$

\sum_{ja} \sum_{jot} (\log (y_{ja jot}) y_{ja jot} + \log (1 - y_{ja jot}) (1 - y_{ja jot})) .

$\sum_i\sum_j \left(\log(y_{ij})y_{ij} + \log(1 - y_{ij})(1-y_{ij})\right).$

y_{i j}

$y_{ij}$

\log (0)

$\log(0)$

0 \log (0)

$0\log(0)$

lim_{x \to 0^{+}} x \log (x)

$\lim_{x \to 0^+}x\log(x)$

Skrócona forma (ważona)

Zauważ, że rozkład dwumianowy nie może faktycznie przyjmować wartości niecałkowitych, ale mimo to możemy obliczyć „logarytm prawdopodobieństwa”, wykorzystując ułamek zaobserwowanych sukcesów w każdej komórce jako odpowiedź i ważąc każdy summand w obliczeniu logarytmu prawdopodobieństwa przez liczba obserwacji w tej komórce.

\sum_{ja} n_{ja} (\log (p_{ja}) \sum_{jot} y_{ja jot} / n_{ja} + \log (1 - p_{ja}) (1 - \sum_{jot} (y_{ja jot} / n_{ja}))

$\sum_in_i \left(\log(p_i)\sum_jy_{ij}/n_i + \log(1 - p_i)(1 - \sum_j(y_{ij}/n_i)\right)$

Jest to dokładnie równe obliczonemu powyżej odchyleniu modelu, które można zobaczyć, wyciągając w miarę możliwości sumę nad w równaniu długiej postaci. $j$

Tymczasem nasycone odchylenie jest inne. Ponieważ nie mamy już odpowiedzi 0-1, nawet z jednym parametrem na obserwację nie możemy uzyskać dokładnie 0. Zamiast tego nasycone prawdopodobieństwo logarytmu modelu wynosi

\sum_{ja} n_{ja} (\log (\sum_{jot} y_{ja jot} / n_{ja}) \sum_{jot} y_{ja jot} / n_{ja} + \log (1 - \sum_{jot} y_{ja jot} / n_{ja}) (1 - \sum_{jot} y_{ja jot} / n_{ja})) .

$\sum_i n_i\left(\log(\sum_jy_{ij}/n_i)\sum_jy_{ij}/n_i + \log(1 - \sum_jy_{ij}/n_i)(1-\sum_jy_{ij}/n_i)\right).$

W twoim przykładzie możesz zweryfikować, że dwukrotność tej kwoty stanowi różnicę między zgłaszanymi wartościami zerowymi i odchyleniami resztkowymi dla obu modeli.

ni = dfShort$nReps
yavg = dfShort$mortalityP
sum.terms <-ni*(log(yavg)*yavg + log(1 - yavg)*(1 - yavg))
# Need to handle NaN when yavg is exactly 0
sum.terms[1] <- log(1 - yavg[1])*(1 - yavg[1])

2*sum(sum.terms)
fitShortP$deviance - fitLong$deviance

Jonny Lomond
źródło

Myślę, że będziesz musiał wyjaśnić wyrażenie dewiacji nasyconych modeli. Log 0 nie działa tak dobrze.

AdamO,

Dzięki, powinienem był wyjaśnić, co miałem na myśli. Dodałem edycję, aby wyjaśnić, że przez 0log (0) mam na myśli 0 w tym kontekście.

Jonny Lomond

OK, ale jestem właściwie zdezorientowany (wybacz mi, nigdy nie opisywałem dewiacji w najdrobniejszych szczegółach): jeśli masz log (y) y - log (1-y) (1-y) jako dewiację modelu nasyconego, nie każdy obserwacja tylko 0?

AdamO,

„Model nasycony” to model wyobrażony z jednym parametrem na obserwację. Tak więc przewidywane prawdopodobieństwo każdej obserwacji wynosi 1 lub 0, w zależności od faktycznie obserwowanej wartości. Zatem w tym przypadku prawdopodobieństwo dziennika modelu nasyconego wynosi 0, dane są jedynymi danymi, które mogą być wygenerowane przez model nasycony.

Jonny Lomond

Tom Wenseleers,