Przejdź do zawartości

Uogólniony model liniowy

Z Wikipedii, wolnej encyklopedii

Uogólniony model liniowy (ang. generalized linear model, GLM) – model statystyczny będący uogólnieniem zwykłej regresji liniowej. GLM uogólnia regresję liniową poprzez powiązanie kombinacji liniowej zmiennych objaśniających ze zmienną objaśnianą za pomocą funkcji wiążącej (ang. link function) oraz poprzez umożliwienie przedstawienia wariancji każdego pomiaru jako funkcji wartości przewidywanej.

Uogólnione modele liniowe zostały sformułowane przez Johna Neldera i Roberta Wedderburna jako sposób na ujednolicenie różnych innych modeli statystycznych, w tym regresji liniowej, regresji logistycznej i regresji Poissona[1]. Zaproponowali oni iteracyjną ważoną metodę najmniejszych kwadratów do estymacji parametrów modelu metodą największej wiarygodności (MLE). MLE jest popularną metodą, która stsoowana jest jako domyślna w wielu pakietach statystycznych. Opracowano również inne podejścia, takie jak regresja bayesowska i dopasowanie metodą najmniejszych kwadratów z wykorzystaniem stabilizacji wariancji.

Idea[edytuj | edytuj kod]

Zwykła regresja liniowa przewiduje wartość oczekiwaną danej nieznanej wielkości (zmiennej objaśnianej, nazywanej też zmienną odpowiedzi) na podstawie kombinacji liniowej obserwowanych wartości (zmiennych objaśniających). Oznacza to, że określona zmiana predyktora prowadzi średnio do proporcjonalnej zmiany zmiennej objaśnianej (model odpowiedzi liniowej). Jest to odpowiednie, gdy zmienna objaśniana zmienia się (z dobrym przybliżeniem) w nieskończoność w dowolnym kierunku lub zmienia się jedynie o stosunkowo niewielką wartość w porównaniu ze zmianami zmiennych predykcyjnych, np. w przypadku wzrostu człowieka.

Jednak założenia te mogą być niewłaściwe w przypadku niektórych typów zmiennych objaśnianych. Na przykład gdy oczekuje się, że zmienna objaśniająca będzie się zmieniać w szerokim zakresie, ale zawsze będzie dodatnia, zmiany predyktorów mogą prowadzić do geometrycznie (tzn. wykładniczo) – a nie liniowo – zmieniających się wartości na wyjściu. Załóżmy na przykład, że liniowy model predykcyjny dowiaduje się na podstawie cząstkowych danych (być może pochodzących głównie z dużych plaż), że spadek temperatury o 10 stopni spowoduje, że plażę będzie odwiedzać o 1000 osób mniej. Jest mało prawdopodobne, aby ten model można było uogólnić na plaże o różnej wielkości. Mówiąc dokładniej, problem polega na tym, że jeśli użyje się modelu do przewidzenia nowej frekwencji przy spadku temperatury o 10 dla plaży, na której regularnie pojawia się 50 plażowiczów, przewidywana będzie niemożliwa liczba plażowiczów wynosząca -950. Logicznie rzecz biorąc, bardziej realistyczny model przewidywałby zamiast tego stałą stopę wzrostu frekwencji na plażach (np. wzrost o 10 stopni prowadziłby do podwojenia frekwencji na plaży, a spadek o 10 stopni prowadziłby do zmniejszenia frekwencji o połowę). Taki model nazywany jest modelem odpowiedzi logarytmicznej lub modelem logarytmiczno-liniowym, ponieważ przewiduje, że logarytm zmiennej objaśnianej będzie się zmieniał się liniowo.

Analogicznie, model liniowy przewidujący prawdopodobieństwo dokonania wyboru tak/nie (zmienna zero-jedynkowa) byłby jeszcze mniej odpowiedni, ponieważ prawdopodobieństwa są ograniczone z obu stron (muszą mieścić się w przedziale od 0 do 1). Wyobraźmy sobie na przykład model, który przewiduje prawdopodobieństwo wyjścia danej osoby na plażę w zależności od temperatury. Rozsądny model mógłby na przykład przewidywać, że zmiana o 10 stopni powoduje, że dana osoba będzie dwa razy bardziej lub mniej skłonna do wyjścia na plażę. Co jednak oznacza „dwa razy większa skłonność” w tym kontekście? Nie może dosłownie oznaczać podwojenia wartości prawdopodobieństwa (np. 50% staje się 100%, 75% staje się 150% itd.). Raczej podwajają się szanse: od szans 2:1, przez 4:1, aż do 8:1 itd. Taki model jest nazywany jest modelem logarytmu szansy lub modelem logistycznym.

Uogólnione modele liniowe obejmują powyższe sytuacje: umożliwiają inne niż rozkład normalny rozkłady zmiennej objaśnianej oraz zakładają, że kombinacją liniową predyktorów i składnika losowego jest argument funkcji wiążącej, a nie sama zmienna objaśniana. Odwołując się do powyższych przykładów: przewidywana liczba plażowiczów będzie zazwyczaj modelowana za pomocą rozkładu Poissona i wiązania logarytmicznego, a prawdopodobieństwo obecności na plaży będzie zazwyczaj modelowane za pomocą rozkładu zero-jedynkowego (lub rozkładu dwumianowego, w zależności od tego jak dokładnie sformułowano problem) i logitowej funkcji wiążącej.

Przegląd[edytuj | edytuj kod]

W uogólnionym modelu liniowym (GLM) zakłada się, że każdy wynik Y zmiennej objaśnianej jest generowany z określonego rozkładu z rodziny wykładniczej – dużej klasy rozkładów prawdopodobieństwa, która obejmuje między innymi rozkład normalny, dwumianowy, Poissona i gamma. Średnia warunkowa μ rozkładu zależy od zmiennych niezależnych X zgodnie z następującym wzorem:

gdzie E(Y | X) jest wartością oczekiwaną Y pod warunkiem X; iloczyn Xβ,nazywany predyktorem liniowym, jest liniową kombinacją nieznanych parametrów β; g jest funkcją wiążącą.

W tym kontekście wariancja jest zazwyczaj funkcją V warunkowej wartości oczekiwanej:


Nieznane parametry β są zazwyczaj szacowane przy użyciu metody największej wiarygodności, największej quasi-wiarygodności lub technik bayesowskich.

Elementy modelu[edytuj | edytuj kod]

GLM składa się z trzech elementów:

1. określonego rozkładu wykorzystywanego do modelowania należącego do rodziny wykładniczej,
2. predyktora liniowego ,
3. funkcji wiążącej , takiej że .

Rozkład prawdopodobieństwa[edytuj | edytuj kod]

Rodzina rozkładów wykładniczych z nadmierną dyspersją jest uogólnieniem rodziny wykładniczej i dyspersyjnej rodziny rozkładów wykładniczych i obejmuje te rodziny rozkładów prawdopodobieństwa, sparametryzowane przez i , których funkcje gęstości f (lub funkcje masy prawdopodobieństwa, dla rozkładów dyskretnych) można wyrazić w postaci:

Parametr dyspersji, , zazwyczaj jest znany i powiązany z wariancją rozkładu. Funkcje , , , , i są również znane.

Do tej rodziny należy wiele powszechnie stosowanych rozkładów, w tym rozkład normalny, wykładniczy, gamma, Poissona, zero-jedynkowy i (dla ustalonej liczby prób) dwumianowy, wielomianowy i ujemny dwumianowy.

Gdy i są skalarami (stosujemy wtedy oznaczenia i ) funkcja gęstości sprowadza się do:

Parametr jest powiązany ze średnią rozkładu. Jeśli , mówimy, że rozkład ma postać kanoniczną (lub postać naturalną). Należy pamiętać, że dowolny rozkład można przekształcić do postaci kanonicznej poprzez przepisanie jako i zastosowanie transformacji ; zawsze można przekształcić w zgodzie z nową parametryzacja, nawet jeśli nie jest iniekcją. Jeżeli dodatkowo i jest znane, nazywa się parametrem kanonicznym (lub parametrem naturalnym) powiązanym ze średnią w następujący sposób:

Dla skalarnych wartości i powyższe sprowadza się się do

W tym scenariuszu można wykazać, że wariancja rozkładu wynosi

i dla skalarnych i

Predyktor liniowy[edytuj | edytuj kod]

Predyktor liniowy η (greckieeta”) to wielkość, która włącza do modelu informacje o zmiennych objaśniających. Predyktor liniowy jest powiązany z wartością oczekiwaną przez funkcję wiążącą.

η jest kombinacją liniową (stąd „liniowy”) nieznanych parametrów β. Rolę współczynników kombinacji liniowej pełni macierz zmiennych objaśniających X. η można zatem wyrazić jako

Funkcja wiążąca[edytuj | edytuj kod]

Funkcja wiążąca (ang. link function) zapewnia związek między predyktorem liniowym a warunkową wartością oczekiwaną zmiennej objaśnianej. Wykorzystuje się wiele funkcji wiążących – ich wybór wynika z kilku czynników. Zawsze istnieje dobrze określona kanoniczna funkcja wiążąca, która jest wyprowadzana z funkcji gęstości zmiennej objaśnianej. Jednak w niektórych przypadkach może mieć sens próba dopasowania dziedziny funkcji wiążącej do zakresu średniej rozkładu zmiennej objaśnianej albo wykorzystanie niekanonicznej funkcji wiążącej ze względów algorytmicznych, na przykład w ramach bayesowskiej regresji probitowej.

Gdy funkcja rozkładu ma parametr kanoniczny wówczas kanoniczna funkcja wiążąca to funkcja wyrażającą za pomocą : W przypadku najpopularniejszych rozkładów średnia jest jednym z parametrów funkcji gęstości rozkładu, a jest funkcją zdefiniowaną powyżej, która odwzorowuje funkcję gęstości na jej postać kanoniczną. Jeżeli funkcja wiążąca jest kanoniczna, to co pozwala być wystarczającą statystyką dla .

W poniższej tabeli przedstawiono kilka powszechnie używanych rozkładów z rodziny wykładniczej, sytuacje, w których są wykorzystywane, ich kanoniczne funkcje wiążące i odwrotności.

Typowe rozkłady, ich wykorzystanie oraz kanoniczne funkcje wiążące
Rozkład Nośnik Typowe wykorzystanie Nazwa funkcji wiążącej Funkcja wiążąca, Odwrotność funkcji wiążącej
normalny zależności liniowe tożsamość
wykładniczy model odpowiedzi wykładniczej, parametry skali ujemna odwrotność
gamma
odwrotny gaussowski odwrotność do kwadratu
Poissona nieujemne liczby całkowite: liczba wystąpień w określonym czasie/przestrzeni logarytm
zero-jedynkowy tak/nie logit
dwumianowy liczba wystąpień w N próbach
wielopunktowy liczba całkowita z przedziału wystąpienie jednej z K kategorii w pojedynczej próbie
K-wymiarowy wektor, gdzie dokładnie jeden element ma wartość 1, a pozostałe to zera
wielomianowy K-wymiarowy wektor o wartościach całkowitych z przedziału liczba wystąpień każdej z K kategorii w N próbach

W przypadku rozkładów wykładniczego i gamma należy w ramach dopasowywania modelu podjąć odpowiednie środki ostrożności, ponieważ dziedzina kanonicznej funkcji wiążącej nie jest tożsama z dopuszczalnym zakresem średniej: predyktor liniowy może być dodatni, co dałoby niemożliwą do uzyskania średnią ujemną. Alternatywą jest użycie niekanonicznej funkcji wiążącej.

W przypadku rozkładów zero-jedynkowego, dwumianowego, wielopunktowego i wielomianowego nośniki rozkładów nie są zgodne z wartościami przyjmowanymi przez zmienną objaśnianą. We wszystkich tych przypadkach dokonuje się predykcji jednego lub większej liczby prawdopodobieństw w przedziale . Taki model jest znany jako regresja logistyczna (lub wielomianowa regresja logistyczna jeśli przewidywane są więcej niż dwie wartości).

Dla rozkładu zero-jedynkowego i dwumianowego parametrem jest prawdopodobieństwo pojedyncze, wskazujące prawdopodobieństwo wystąpienia pojedynczego zdarzenia. W tym przypadku nadal spełniony jest podstawowy warunek uogólnionego modelu liniowego: chociaż pojedynczy wynik będzie zawsze wynosił 0 lub 1, to wartość oczekiwana będzie prawdopodobieństwem wystąpienia 1 („tak”). Podobnie w rozkładzie dwumianowym wartością oczekiwaną jest Np, co oznacza, że będzie się przewidywać oczekiwaną proporcję sukcesów.

W przypadku rozkładów wielopunktowych i wielomianowych przewidywany jest wektor K prawdopodobieństw, z zastrzeżeniem, że prawdopodobieństwa muszą sumować się do 1. W przypadku rozkładu wielomianowego i wektorowej postaci rozkładu wielopunktowego oczekiwane wartości elementów wektora można powiązać z przewidywanymi prawdopodobieństwami podobnie jak w przypadku rozkładu dwumianowego i rozkładu zero-jedynkowego.

Estymacja[edytuj | edytuj kod]

Metoda maksymalnej wiarygodności[edytuj | edytuj kod]

Oszacowania największej wiarygodności można znaleźć za pomocą iteracyjnie przeważonego algorytmu najmniejszych kwadratów/metody Newtona z aktualizacjami postaci:

lub za pomocą metody ocen Fishera (Fisher scoring):

gdzie jest funkcją wynikową (gradientem log-wiarygodności); jest obserwowaną macierzą informacyjną (ujemną macierzą Hessego), a jest macierzą informacyjną Fishera. Warto zauważyć, że jeśli korzystamy z kanonicznej funkcji wiążącej, to te macierze takie same[2].

Metody bayesowskie[edytuj | edytuj kod]

Zwykle rozkładu a posteriori nie można wyrazić w sposób jawny, dlatego aproksymuje się go za pomocą przybliżeń Laplace'a lub z wykorzystaniem jakiejś formy próbkowania Monte Carlo łańcuchami Markowa,np. próbkowania Gibbsa.

Przykłady[edytuj | edytuj kod]

Ogólne modele liniowe[edytuj | edytuj kod]

Podobieństwo nazwy między uogólnionymi modelami liniowymi i ogólnymi modelami liniowymi, dwoma szerokimi grupami modeli statystycznych, może prowadzić do nieporozumień. Współtwórca John Nelder wyraził ubolewanie z powodu zastosowania tego nazewnictwa[3].

Ogólny model liniowy można postrzegać jako szczególny przypadek uogólnionego modelu liniowego z tożsamościową funkcją wiążącą i warunkowym rozkładem normalnym zmiennej objaśnianej. Ogólne modele liniowe powstały wcześniej pod względem historycznym. Umożliwiają one uzyskanie bardziej precyzyjnych wyników, podczas gdy wyniki dla uogólnionego modelu liniowego z funkcją wiążącą inną niż tożsamość są asymptotyczne (zwykle dobrze sprawdzają się w przypadku dużych próbek).

Regresja liniowa[edytuj | edytuj kod]

Prostym, bardzo ważnym przykładem uogólnionego modelu liniowego (będącym jednocześnie przykładem ogólnego modelu liniowego) jest regresja liniowa. W regresji liniowej zastosowanie estymatora najmniejszych kwadratów uzasadnione jest twierdzeniem Gaussa-Markowa, które nie zakłada, że rozkład jest normalny.

Z perspektywy uogólnionych modeli liniowych warto jednak założyć, że rozkład jest normalny o stałej wariancji, a funkcja wiążąca to tożsamość, która jest funkcją kanoniczną, jeżeli wariancja jest znana. Przy tych założeniach estymator metodą najmniejszych kwadratów jest estymatorem największej wiarygodności.

W przypadku rozkładu normalnego uogólniony model liniowy ma wzór jawny dla szacunków największej wiarygodności, co jest wygodne. Większość innych modeli nie ma wzoru jawnego.

Dane binarne[edytuj | edytuj kod]

Gdy zmienna objaśniana Y jest binarna (przyjmuje tylko wartości 0 i 1), na ogół jako funkcję rozkładu wybiera się rozkład zero-jedynkowy, a μ jest wówczas prawdopodobieństwem, że Y przyjmie wartość jeden.

Istnieje kilka popularnych funkcji wiążących dla rozkładu dwumianowego.

Logitowa funkcja wiążąca[edytuj | edytuj kod]

Najbardziej typową funkcją wiążącą jest kanoniczne wiązanie logitowe:

Modele z taką konfiguracją nazywane są modelami regresji logistycznej (lub modelami logitowymi).

Probitowa funkcja wiążąca[edytuj | edytuj kod]

Alternatywnie jako funkcję wiążącą można zastosować odwrotność dowolnej dystrybuanty rozkładu ciągłego, ponieważ dystrybuanta przyjmuje wartości z przedziału . Jeżeli wybierzemy dystrybuantę standaryzowanego rozkładu normalnego , to otrzymamy model probitowy. Jego funkcja wiążąca to

Modele probitowe zwracają wyniki podobne do modeli logitowych, lecz ich zaletą jest to, że są w niektórych sytuacjach bardziej praktyczne. W paradygmacie bayesowskim, w którym parametry mają normalne rozkłady a priori, relacja pomiędzy rozkładami normalnymi a probitową funkcją wiążącą sprawia, że model probitowy można dopasować przy użyciu próbkowania Gibbsa, podczas nie jest to na ogół możliwe dla modelu logitowego.

Funkcja cloglog[edytuj | edytuj kod]

Jako funkcji wiążącej można również użyć funkcji cloglog:

Funkcja cloglog jest asymetryczna i jej zastosowanie może sprawić, że model zwróci wyniki odmienne niż te otrzymywane z wykorzystaniem funkcji logit lub probit[4]. Model cloglog pasuje do zastosowań, w których obserwujemy albo zero zdarzeń (np. defektów), albo jedno lub więcej, zaś liczba zdarzeń jest zgodna z rozkładem Poissona[5]. W rozkładzie Poissona

gdzie μ jest liczbą dodatnią oznaczającą oczekiwaną liczbę zdarzeń. Jeśli p oznacza proporcję obserwacji z co najmniej jednym zdarzeniem, to jego dopełnienie (ang. complement) wynosi:

i stąd

Model liniowy wymaga, aby zmienna odpowiedzi przyjmowała wartości z całej linii rzeczywistej. Ponieważ μ musi być dodatnie, możemy to wymusić stosując logarytm. W ten sposób powstaje przekształcenie „cloglog” (ang. complementary log-log).

Wiązanie tożsamościowe[edytuj | edytuj kod]

Wiązanie tożsamościowe g(p) = p bywa również czasami używane dla danych o charakterze binarnym. W takiej sytuacji powstaje liniowy model prawdopodobieństwa. Jednak zastosowanie tożsamościowej funkcji wiążącej może prowadzić do przewidywania nonsensownych „prawdopodobieństw” mniejszych od zera lub większych niż jeden. Można tego uniknąć, stosując transformację cloglog, probit lub logit (lub dowolną inną funkcję kwantylową). Podstawową zaletą wiązania tożsamościowego jest to, że model sprowadzony jest matematycznie do prostego modelu liniowego, zaś przebieg popularnych funkcji wiążących jest w pobliżu p = 0,5 w przybliżeniu liniowy i zbliżony do funkcji tożsamościowej.

Funkcja wariancji[edytuj | edytuj kod]

Funkcja wariancji dla modelu „quasi-dwumianowego” to:

gdzie parametr dyspersji τ wynosi dokładnie 1 dla rozkładu dwumianowego. Rzeczywiście, standardowe prawdopodobieństwo dwumianu pomija τ. Jeśli jest obecny, model nazywa się „quasi-dwumianowym”, a zmodyfikowaną funkcję wiarygodności nazywa się quasi-wiarygodnością, ponieważ generalnie nie jest to prawdopodobieństwo odpowiadające żadnej rzeczywistej rodzinie rozkładów prawdopodobieństwa. Jeśli τ przekracza 1, mówimy, że model wykazuje nadmierną dyspersję .

Regresja wielomianowa[edytuj | edytuj kod]

Przypadek dwumianowy można łatwo rozszerzyć, aby uwzględnić w odpowiedzi rozkład wielomianowy (także uogólniony model liniowy dla zliczeń z ograniczoną sumą). Zwykle można to zrobić na dwa sposoby:

Porządkowa zmienna objaśniana[edytuj | edytuj kod]

Jeśli zmienna odpowiedzi jest porządkowa, można dopasować funkcję modelu w postaci:

dla m > 2. Różne powiązania g prowadzą do modeli regresji porządkowej, takich jak modele proporcjonalnych szans lub modele uporządkowanego probitu .

Nominalna zmienna objaśniana[edytuj | edytuj kod]

Jeżeli zmienna objaśniana jest nominalna lub dane nie spełniają założeń modelu porządkowego, można dopasować model o postaci:

dla m > 2. Różne wiązania g prowadzą do wielomianowych modeli logitowych lub wielomianowych modeli probitowych. Są one bardziej ogólne niż modele porządkowe; w przypadku takich modeli szacuje się więcej parametrów.

Dane pokazujące zliczenia[edytuj | edytuj kod]

Innym przykładem uogólnionych modeli liniowych jest regresja Poissona, która modeluje dane z wykorzystaniem rozkładu Poissona. Funkcją wiążącą jest zazwyczaj logarytm, wiązanie kanoniczne.

Funkcja wariancji jest proporcjonalna do średniej:

gdzie parametr dyspersji τ jest zwykle ustalony na dokładnie 1. Jeśli tak nie jest, powstały model quasi-wiarygodności jest często opisywany jako Poissona z nadmierną dyspersją lub quasi-Poissona .

Rozszerzenia[edytuj | edytuj kod]

Dane skorelowane lub skupione[edytuj | edytuj kod]

Standardowy GLM zakłada, że obserwacje są nieskorelowane. Opracowano rozszerzenia umożliwiające korelację między obserwacjami, co ma miejsce na przykład w badaniach podłużnych lub gdy dane są skupione w klastry:

  • Uogólnione równania estymujące (ang. generalized estimating equations, GEE) pozwalają na korelację między obserwacjami bez użycia wyraźnego modelu prawdopodobieństwa pochodzenia korelacji, więc nie ma wyraźnego prawdopodobieństwa. Są one odpowiednie, gdy efekty losowe i ich wariancje nie są przedmiotem zainteresowania, ponieważ pozwalają na korelację bez wyjaśniania jej pochodzenia. Nacisk położony jest na oszacowanie średniej reakcji w populacji („efekty „uśrednione przez populację”), a nie na parametrach regresji, które umożliwiłyby przewidywanie wpływu zmiany jednego lub większej liczby składników X na daną osobę. GEE są zwykle używane w połączeniu ze standardowymi błędami Hubera-White'a[6][7].
  • Uogólnione liniowe modele mieszane (ang. generalized linear mixed models, GLMM) stanowią rozszerzenie GLM, które uwzględnia efekty losowe w predyktorze liniowym, w efekcie zwracając model wyjaśniający pochodzenie korelacji. Uzyskane w ten sposób oszacowania parametrów „specyficznych dla podmiotu” są odpowiednie, gdy nacisk położony jest na oszacowanie wpływu zmiany jednego lub większej liczby składników X na daną jednostkę. GLMM nazywane są także modelami wielopoziomowymi i modelami mieszanymi. Ogólnie rzecz biorąc, dopasowanie GLMM jest bardziej złożone obliczeniowo i intensywniejsze niż dopasowanie GEE.

Uogólnione modele addytywne[edytuj | edytuj kod]

Uogólnione modele addytywne (ang. generalized additive models, GAM) to kolejne rozszerzenie GLM, w którym predyktor liniowy η nie jest ograniczony do liniowości we współzmiennych X, ale jest sumą funkcji wygładzających nałożonych na xi:

Funkcje wygładzające fi są estymowane na podstawie danych. Wymaga to dużej ilości danych i mocy obliczeniowej[8][9].

Przypisy[edytuj | edytuj kod]

  1. J.A. Nelder, R.W.M. Wedderburn, Generalized Linear Models, „Journal of the Royal Statistical Society. Series A (General)”, 135 (3), 1972, s. 370–384, DOI10.2307/2344614, ISSN 0035-9238, JSTOR2344614.
  2. Peter McCullagh, John Ashworth Nelder, Generalized linear models, wyd. 2nd ed, Monographs on statistics and applied probability, London [etc.]: Chapman & Hall, 1991, rozdział 2, ISBN 978-0-412-31760-6 [dostęp 2024-06-04].
  3. Stephen Senn. A conversation with John Nelder. „Statistical Science”. 18 (1), s. 118–131, 2003. DOI: 10.1214/ss/1056397489. Cytat: I suspect we should have found some more fancy name for it that would have stuck and not been confused with the general linear model, although general and generalized are not quite the same. I can see why it might have been better to have thought of something else.. .
  4. Complementary Log-log Model.
  5. Which Link Function — Logit, Probit, or Cloglog?. Bayesium Analytics, 2015-08-14. [dostęp 2019-03-17]. (ang.).
  6. Scott L. Zeger. Models for Longitudinal Data: A Generalized Estimating Equation Approach. „Biometrics”. 44 (4), s. 1049–1060, 1988. International Biometric Society. DOI: 10.2307/2531734. PMID: 3233245. JSTOR: 2531734. 
  7. James Hardin: Generalized Estimating Equations. London, England: Chapman and Hall/CRC, 2003. ISBN 1-58488-307-3. LCCN Joseph Hilbe.
  8. Trevor Hastie, Robert Tibshirani, Generalized additive models, Boca Raton, Fla: Chapman & Hall/CRC, 1999, ISBN 978-0-412-34390-2 [dostęp 2024-06-04].
  9. Simon N. Wood, Generalized additive models: an introduction with R, Texts in statistical science, Boca Raton: Chapman & Hall/CRC, 2006, ISBN 978-1-58488-474-3 [dostęp 2024-06-04].