Lineární regresní model (VJ REGMOD-2) Základní informace V rámci této výukové jednotky si nadefinujeme lineární regresní model a seznámíme se s typy proměnných využitelných jako prediktory (vysvětlující proměnné) v takovém modelu. Ukázané principy budou užitečné v průběhu celého předmětu Regresní modelování, neboť principy lineárního regresního modelování se uplatní i v modelovacích přístupech, které budou ukázány v následujících výukových jednotkách. U studentů se předpokládá znalost základních pojmů biostatistiky, které byly rekapitulovány v předcházející výukové jednotce. K pochopení výpočetních principů je nezbytná základní znalost počítání s maticemi a vektory.
Výstupy z výuky konkrétní výukové jednotky Po prostudování učebního textu této výukové jednotky studenti • • • •
definují lineární regresní model , vysvětlí předpoklady lineárního regresního modelu, použijí různé typy proměnných (spojité i kategoriální) při konstrukci modelu, uvedou příklady běžně užívaných regresních modelů.
1. Motivace Lineární regresní modely využíváme pro modelování (vysvětlení pozorovaných hodnot) spojité výsledkové proměnné (výsledku, závisle proměnné). Výsledek se snažíme vysvětlit prostřednictvím jednoho nebo více prediktorů (nezávisle proměnných, vysvětlujících proměnných). Prediktory mohou být buď rovněž spojité, nebo kategoriální. Uvažujme nejprve nejjednodušší situaci, kdy se snažíme určit vztah mezi dvěma spojitými proměnnými Příklad: Našim úkolem je vysvětlit u pacientů sérovou koncentrací 25-hydroxyvitaminu D (pro jednoduchost dále značíme jako vitamin D) prostřednictvím indexu tělesné hmotnosti (body mass index, BMI). Údaje pacientů pocházejí z datového souboru [vitaminD]. Situace je znázorněna na obrázku 2.1. Pro každého pacienta máme změřeny obě zmíněné veličiny, které poté znázorníme na x-y grafu. Takovou sadou bodů můžeme proložit (v tomto jednoduchém případě i intuitivně) přímku, která co nejlépe reprezentuje vztah koncentrace vitaminu D a BMI. Je přirozené volit přímku tak, aby byla vzdálenost jednotlivých bodů od přímky co nejmenší. V praxi se jako výpočetní metoda pro určení koeficientů takové přímky uplatňuje takzvaná metoda nejmenších čtverců, která minimalizuje součet druhých mocnin vzdáleností mezi jednotlivými body (naměřenými pozorováními) a hledanou přímkou. Běžný statistický software nám tedy umožňuje najít „ideální“ polohu přímky, kterou lze popsat následovně: koncentrace vitaminu D = 111,1 – 2,4·BMI Regresní koeficient odhadnutý jako 111,1 se nazývá absolutní člen (posun, anglicky intercept) regresní přímky. Druhý koeficient této jednoduché rovnice se nazývá směrnice (sklon, anglicky slope). V této výukové jednotce se seznámíme se základními definicí lineárního regresního modelu a s jeho předpoklady. Ukážeme si principy odhadů koeficientů tohoto modelu a testování hypotéz založených na parametrech tohoto modelu.
100 80 60 20
40
Vitamin D
20
25
30
35
BMI Obr. 2.1 Znázornění vztahu mezi indexem tělesné hmotnosti (BMI) a sérovou koncentrací vitaminu D.
2. Jak definujeme lineární regresní model? 2.1. Lineární regresní model Předpokládejme na chvilku, že existuje pro všechna pozorování přesný vztah mezi dvěma (nenáhodnými) veličinami y (výsledek) a x (prediktor): y = β 0 + β1 x
Takto definovaný vztah mezi veličinami však na reálných datech (zejména z biologie nebo medicíny) v praxi pozorujeme stěží. Pro regresní modelování se proto využívá následujícího vztahu, který v sobě již zahrnuje náhodnou veličinu ε (reziduum) reprezentující odchylku od uvedeného ideálního vztahu. Y označuje výsledek (náhodnou veličinu), x označuje prediktor (nenáhodnou, přesně změřenou veličinu). Předpokládejme tedy, že pro jednotlivá pozorování (např. pacienty, lokality, apod.) číslované prostřednictvím indexu i od 1 do n (celkový počet pozorování) platí: Yi = β 0 + β1 xi + ε i
O reziduích budeme předpokládat, že jsou •
nesystematické – střední hodnota reziduí je rovna 0: Eε i = 0 pro i = 1,...,n
•
homogenní v rozptylu – rozptyl reziduí je pro všechna pozorování stejný: Dε i = σ 2 > 0 pro i = 1,...,n
•
jsou vzájemně nekorelované: C (ε i , ε j ) = 0 pro i ≠ j; i, j = 1,...,n
(2.1)
Pro jeden prediktor x se regresní koeficienty značí β0 a β1, jedná se o zmíněný absolutní člen a směrnici regresní přímky. Uvedený vztah lze jednoduše rozšířit na větší počet (p) prediktorů (pak máme celkem k = p + 1 parametrů včetně β0, absolutního členu). Dostáváme definici vícenásobného regresního modelu (multiple regression): Yi = β 0 + β1 xi1 + ... + β p xip + ε i
(2.2)
Rozepsáno do vztahů pro očekávané hodnoty (predikce) jednotlivých pozorování i = 1,...,n: EY1 = β 0 + β1 x11 + ... + β p x1 p EY2 = β 0 + β1 x21 + ... + β p x2 p M EYn = β 0 + β1 xn1 + ... + β p xnp
Tuto soustavu vztahů můžeme zapsat jako následující vztah využívající násobení matic: systematická složka
výsledek
náhodná složka
Y1 1 x11 L x1 p β 0 ε 1 M M + M M = M M Y 1 x L x β ε np p n n n1 matice plánu
regresní koeficienty
Vektor výsledků, matici plánu, vektor regresních koeficientů a vektor reziduí označíme po řadě Y, X, β a ε. Maticový zápis regresních rovnic nám umožní zjednodušit definice potřebných statistik. Y = Xβ + ε
(2.3)
2.2. Normálně rozdělený výsledek Definice lineárního regresního modelu doposud neobsahovala specifikaci konkrétní náhodné veličiny. Doplnění rozdělení reziduí do definice regresního modelu nám umožní předvídat rozdělení výsledků, konstruovat intervaly spolehlivosti a testovat statistické hypotézy. Uvažujme tedy, že rezidua mají normální rozdělení s nulovou střední hodnotou a rozptylem σ2. Dále budeme předpokládat, že rozdělení reziduí pro jednotlivá pozorování jsou vzájemně nezávislá. Yi = β 0 +
p
∑β
j X ij
+ ε i , i = 1,..., n
j =1
ε i ~ N (0, σ 2 )
(2.4)
2.3. Odhady parametrů regresního modelu Pro odhad parametrů lineárního regresního modelu zpravidla využíváme metodu nejmenších čtverců. Pro takový odhad musí platit, že minimalizuje druhé mocniny rozdílů mezi jednotlivými pozorováními výsledku a regresní přímkou (nebo obecně regresní nadrovinou – přímka ve dvourozměrném prostoru,
rovina ve třírozměrném prostoru, atd.). Následující vztah umožňuje odhadnout parametry regresního modelu metodou nejmenších čtverců (třeba i „ručně“, spíše však s pomocí maticového kalkulátoru). Budeme jej značit jako βˆ OLS (stříška označuje, že se jedná o odhad příslušného parametru, index OLS pochází z anglického ordinary least squares, tedy metoda (obyčejných) nejmenších čtverců). βˆ OLS = ( X ′X) −1 X ′Y
(2.5)
Důkaz tohoto tvrzení viz Statistické modelování – Lineární regresní model, věta 3.3. Lze dále dokázat, že odhad metodou nejmenších čtverců je nejlepší (ve smyslu nejmenšího rozptylu odhadu) nestranný (střední hodnota odhadu je rovna hledanému parametru) lineární odhad. Rozptyl odhadu regresních koeficientů je:
Dβˆ OLS = σ 2 ( X′X) −1
(2.6)
a tedy známe rozdělení odhadu parametrů modelu metodou nejmenších čtverců: βˆ OLS ~ N k (β, σ 2 ( X′X) −1 )
(2.7)
Veličina, kterou jsme minimalizovali prostřednictvím metody nejmenších čtverců, se nazývá reziduální součet čtverců a je definována následovně. S e = Y ′Y − β′OLS X′Y
(2.8) 2
Reziduální součet čtverců nám umožní odhadnout rozptyl regresního modelu (reziduí) σ : s2 =
Se n−k
(2.9)
2.4. Základní statistické testy v regresním modelu Klíčovou úlohou v biostatistice představuje testování statistických hypotéz. V této kapitole se seznámíme se dvěma základními třídami hypotéz, které můžeme v regresních modelech testovat. Příklad: Nejjednodušší (a zřejmě nejčastěji testovanou) nulovou hypotézou je rovnost některého z regresních koeficientů 0. Tak můžeme například testovat klinickou hypotézu, že se koncentrace vitaminu D v krevním séru mění s rostoucím indexem tělesné hmotnosti. Nulová (H0) a alternativní (H1) hypotézy tedy vypadají následovně: H0: β1 = 0 H1: β1 ≠ 0 Tento jednoduchý příklad reprezentuje první třídu hypotéz, kdy testujeme rovnost lineární kombinace regresních koeficientů a libovolné konstanty. To lze v případě modelu s jedním prediktorem (dvěma parametry včetně absolutního členu) zapsat jako vektorový součin c′ ⋅ βˆ OLS , kde c′ = (0 1)
βˆ βˆ OLS = OLS ,0 βˆ OLS ,1
c′ ⋅ βˆ OLS = βˆ OLS ,1
Obecně můžeme zkonstruovat testovou statistiku T pro regresní model s k parametry, v tomto případě mají oba vektory c, βˆ OLS velikost k ×1. Statistika vypadá následovně: T=
c′βˆ OLS − c′β s c′(X′X)−1 c
~ t (n − k )
(2.10)
a má studentovo rozdělení s n – k stupni volnosti. Pokud testujeme hypotézu H0: c′ ⋅ βˆ = 0 H1: c′ ⋅ βˆ ≠ 0,
pak nulovou hypotézu zamítáme na hladině významnosti α, pokud c′βˆ OLS s c′(X′X) −1 c
≥ t1−α / 2 (n − k )
kde t1−α / 2 (n − k ) je 1 - α/2 kvantil studentova rozdělení s n – k stupni volnosti. Uvedený test nám však nepomůže, pokud chceme otestovat rovnost celého vektoru parametrů nulovému vektoru. Takovýto test nám pomůže s rozhodnutím, zda model jako celek dokáže vysvětlit významnou část variability výsledkové veličiny, případně zda umí významnou míru variability vysvětlit kategoriální výsledková proměnná, která je v matici plánu reprezentována několika sloupci (viz dále) a tedy je popsána několika parametry ve vektoru β. Zavedeme nejprve blokové značení pro následující vektory a matice: β = ( β1 ,..., β m , β m +1 ,..., β k )'
= β1′
= β′2
βˆ β V V βˆ OLS = OLS,1 β = 1 ( X′X) −1 = 11 12 βˆ V V 22 21 β2 OLS,2 Můžeme testovat následující hypotézu H0: β2 = x H1: β2 ≠ x kde x je vektor reálných čísel. Testovou statistikou je
F=
1 (βˆ OLS,2 − βˆ 2 )' V22−1 (βˆ OLS,2 − βˆ 2 ) ~ F (k − m, n − k ) s ( k − m) 2
(2.11)
2.5 Koeficient determinace Jako veličina užitečná pro první náhled na přínosnost modelu se zavádí tzv. koeficient determinace. Tato veličina odpovídá na otázku, jakou část z celkové variability přítomné ve výsledkové proměnné
jsme dokázali prostřednictvím vytvořeného regresního modelu vysvětlit. Nejprve si tedy zavedeme veličinu udávající celkovou variabilitu výsledkové proměnné: n
ST = ∑ (Yi − Y ) 2 i =1
(2.12)
S veličinou udávající nevysvětlenou variabilitu výsledkové proměnné jsme se již setkali – jedná se o reziduální součet čtverců: n
S e = ∑ (Yi − Yˆi ) 2 i =1
(2.13)
Můžeme si všimnout, že celkovou variabilitu výsledkové proměnné lze považovat za reziduální součet čtverců nejjednoduššího modelu, kde jediným parametrem je střední hodnota proměnné Y, kterou můžeme odhadnout pomocí výběrového průměru (viz kapitola 4.1). Koeficient determinace již pak jednoduše získáme výpočtem z předchozích dvou veličin:
R2 = 1−
Se ST
(2.14)
3. Předpoklady regresních modelů Ze samotné definice lineárního regresního modelu vyplývá několik předpokladů. Tyto předpoklady, které mohou být v praxi často omezující, se v dalších výukových jednotkách tohoto kurzu naučíme překonávat. Následující přehled je tak zároveň rekapitulací klíčových předpokladů lineárních regresních modelů i motivací ke studiu dalších výukových jednotek, ve kterých budou představeny pokročilejší modelovací postupy. 1. Linearita modelu U popsaného regresního modelu předpokládáme, že očekávaná hodnota výsledku je dána lineární kombinací popsaných parametrů. V následující kapitole (VJ2 – kapitola 4.3) si nicméně ukážeme, že není bezpodmínečně nutná linearita s ohledem na hodnoty prediktorů – hodnoty prediktoru můžeme vložit jako transformované druhou nebo vyšší mocninou a dosáhnout tak polynomiální závislosti výsledku na prediktoru. 2. Aditivita účinků jednotlivých prediktorů Prozatím jsme předpokládali, že účinek nějakého prediktoru je nezávislý na hodnotách ostatních prediktorů. To však v praxi nemusí platit a tento předpoklad může být omezující. Ve výukové jednotce VJ3 – Praktické otázky vícenásobné lineární regrese si ukážeme, jak lze toto omezení překlenout prostřednictvím tzv. interakčních členů. 3. Rezidua mají normální rozdělení s nulovou střední hodnotou a konstantním rozptylem Ve třídě lineárních modelů předpokládáme normální rozdělení reziduí (a z toho vyplývající rozdělení výsledku podmíněné hodnotami prediktorů). To opět nemusí být vždy vyhovující, zejména pro výsledkové proměnné kategoriálního typu. Ve výukové jednotce VJ5 – Logistický regresní model a jiné zobecněné lineární modely se setkáme s třídou zobecněných lineárních modelů, které nám dávají mnohem větší flexibilitu s ohledem na rozdělení výsledkové proměnné.
4. Pozorování jsou vzájemně nezávislá Základní biostatistické metody včetně lineárních regresních modelů předpokládají, že rezidua jsou vzájemně nezávislé proměnné. To opět v praxi nemusí být pravda. Například při dlouhodobém sledování pacientů jsou hodnoty nějakého znaku (např. krevního tlaku, biochemického ukazatele) získané od jednoho pacienta v různých časech zřejmě vzájemně podobnější než hodnoty získané od různých pacientů. To obnáší jistou korelaci mezi různými pozorováními u stejného pacienta a tedy porušení tohoto předpokladu. Řešením je v takovém případě použít třídu tzv. smíšených modelů, které umožňují modelovat korelaci v rámci shluků podobnějších pozorování (nejen pacientů, ale například jednoho lékaře, zdravotnického zařízení apod.). S těmi se seznámíme ve výukové jednotce VJ6 – Smíšené modely.
4. Prediktory různých datových typů 4.1. Konstanta Ve většině prakticky využívaných lineárních regresních modelů je přítomen absolutní člen. Absolutnímu členu odpovídá sloupec jedniček v matici plánu X. Nejjednodušší model je pak samozřejmě takový, který v matici plánu obsahuje právě jen sloupec jedniček a jehož regresním koeficientem je právě jen absolutní člen β0. Parametr tohoto modelu můžeme odhadnout a tento odhad je roven výběrovému průměru hodnot výsledkové proměnné Y. Regresní model a očekávané hodnoty výsledku jsou dány následujícími vztahy:
EY1 = β 0
Y1 1 ε1 β M = M + 0 M Y 1 ε n n
EY2 = β 0 M EYn = β 0
30 15
20
25
BMI
35
40
45
Ukázku regresní přímky pro takový jednoduchý model naleznete na obrázku 2.2 (využitá data ve čtvrté kapitole pochází z datového souboru [heartdisease]). Všimněte si, že na ose x není uveden žádný prediktor (což by nedávalo v tomto případě smysl), ale pouze pořadové číslo příslušného pozorování.
0
100
200
300
400
Pořadové číslo Obr. 2.2 Modelování indexu tělesné hmotnosti: odhad absolutního členu.
4.2. Spojitý prediktor V kapitole 1 jsme si na příkladu již ukázali modelování jednoduché závislosti výsledku na spojité proměnné. Závislost je v takovém případě vyjádřena regresní přímkou – s každou jednotkou prediktoru narůstá (nebo klesá) očekávaná hodnota výsledkové proměnné o hodnotu regresního koeficientu – sklonu regresní přímky. V následujícím zápisu regresního modelu je tímto koeficientem β1.
EY1 = β 0 + β1 x1
Y1 1 x1 ε β 0 1 = + M M M M Y 1 x β1 ε n n n
EY2 = β 0 + β1 x2 M EYn = β 0 + β1 xn
30 15
20
25
BMI
35
40
45
Obrázek 2.3 ukazuje závislost hodnoty indexu tělesné hmotnosti (BMI) na procentuálním vyjádření podílu tukové tkáně. Závislost je v souladu s definovaným modelem vyjádřena přímkou.
10
15
20
25
30
35
40
Podíl tukové tkáně Obr. 2.3 Závislost indexu tělesné hmotnosti (BMI) na podílu tukové tkáně: lineární model.
V praxi se setkáme i se situací, kdy chceme modelovat závislost výsledku na spojité proměnné, lineární model reprezentovaný přímkou však nemusí být adekvátní. Regresní modelování umožňuje snadno modelovat regresní křivku polynomem vyššího stupně. Ukážeme si modelování kvadratické závislosti. V takovém případě zahrneme do matice plánu další sloupec s hodnotami druhé mocniny původní proměnné. Na obrázku 2.4 je znázorněna obdobná závislost jako na předchozím obrázku, nyní však již s nelineárním kvadratickým modelem. ε 1 M β1 + M x12 β 2 ε n
EY2 = β 0 + β1 x2 + β 2 x22 M EYn = β 0 + β1 xn + β 2 xn2
30 15
20
25
BMI
35
40
45
Y1 1 x1 M = M Y 1 x n n
EY1 = β 0 + β1 x1 + β 2 x12
x12 β 0
10
15
20
25
30
35
40
Podíl tukové tkáně Obr. 2.4 Závislost indexu tělesné hmotnosti (BMI) na podílu tukové tkáně: kvadratický model.
4.3. Kategoriální prediktor Neméně užitečný model zahrnuje prediktor kategoriální. Ukažme si takový model na příkladu, ve kterém se snažíme modelovat podíl tukové tkáně v procentech v závislosti na kategorii dle indexu tělesné hmotnosti (podváha, normální váha, nadváha, obezita). Příslušná data jsou znázorněna na obrázku 2.5. Do matice plánu samozřejmě není možné vložit přímo kategoriální proměnnou. Proto musíme tuto kategoriální proměnnou před použitím v regresním modelu převést na sadu indikátorových (dummy) proměnných. Pro jednotlivé kategorie původní proměnné (s výjimkou první) zavedeme indikátorové proměnné, které nabývají hodnoty 1, pokud původní proměnná nabývá příslušné hodnoty, a 0 jinak. První kategorie původní proměnné je pak reprezentována nulovou hodnotou všech indikátorových proměnných zároveň. Situace je na příkladu ukázána v tabulce 2.1.
Tabulka 2.1 Příklad převodu kategoriální proměnné na sadu nových indikátorových proměnných. V posledním sloupci je uveden vztah pro očekávanou hodnotu výsledku pro příslušné pozorování. Nové proměnné Původní proměnná kategorie BMI
Indikátor: Normální váha
Indikátor: Nadváha
Indikátor: Obezita
Podváha
0
0
0
EYi = β 0
Normální váha
1
0
0
EYi = β 0 + β1
Nadváha
0
1
0
EYi = β 0 + β 2
Obezita
0
0
1
EYi = β 0 + β 3
Příslušný řádek matice plánu pak pro jednotlivá pozorování obsahuje jedničku ve druhém, třetím, nebo čtvrtém sloupci pro pacienty s normální váhou, nadváhou a obezitou. Pacienti s podváhou mají tedy očekávanou hodnotu výsledku rovnu koeficientu β0, u pacientů s normální váhou, nadváhou nebo obezitou se přidává ještě regresní koeficient β1, β2 nebo β3. x12 M xn 2
β x13 0 ε 1 β M 1 + M β xn3 2 ε n β3
Podváha
EYi = β 0
β 0 + β1 Nadváha EYi = β 0 + β 2 Obezita EY = β + β i 0 3 Normální EYi =
40 30 20
β1
10
β3
β2
β0
0
Podíl tukové tkáně [%]
50
Y1 1 x11 M = M M Y 1 x n1 n
Podváha
Normální
Nadváha
Obezita
Obr. 2.5 Závislost podílu tukové tkáně na kategorii tělesné hmotnosti: znázornění odhadnutých koeficientů v modelu s kategoriálním prediktorem.
5. Příklady základních biostatistických modelů 5.1 T-test Prostřednictvím t-testu se snažíme testovat klinickou hypotézu o rozdílnosti střední hodnoty mezi dvěma skupinami. Představme si například randomizovanou klinickou studii nového léku, který má za cíl snížení krevního tlaku. Pacienty náhodně rozdělíme do dvou skupin, první skupině podáváme placebo (neaktivní látku), druhé skupině studovaný lék. Po nějaké době naměříme v obou skupinách pacientů hodnoty krevního tlaku. Tyto hodnoty tedy představují hodnoty výsledkové proměnné (Y) v pomyslném regresním modelu. Co jsou prediktory? Jako obvykle budou v prvním sloupci matice plánu X jedničky. Pokud by to byl jediný sloupec matice plánu, jediný regresní koeficient by odpovídal výběrovému průměru krevního tlaku všech pacientů. My přidáme do matice plánu ještě druhý sloupec, který bude vlastně indikátorovou proměnnou příslušnosti ke druhé skupině. Pokud si prvky vektoru β (regresní koeficienty) označíme jako µ (střední hodnota krevního tlaku ve skupině, které bylo podáváno placebo) a α (změna krevního tlaku po podání léku), dostáváme následující vztahy pro očekávané hodnoty krevního tlaku v obou skupinách:
1 M 1 X= 1 M 1
0 M EYi = µ 0 1 M EYi = µ + α 1
Nulovou hypotézou je pak samozřejmě nulová změna v souvislosti s podáváním léku: H0: α = 0 H1: α ≠ 0 Tuto hypotézu lze otestovat prostřednictvím testové statistiky uvedené v kapitole 2.4. 5.2 Analýza rozptylu Klinická hypotéza pro analýzu rozptylu je velmi podobná – snažíme se ukázat rozdíl mezi zkoumanými skupinami, v rámci analýzy rozptylu však zkoumáme více než dvě skupiny. Místo jediného parametru α tedy zavádíme pro m skupin m – 1 parametrů α1,..., αm-1:
1 M 1 1 M X= 1 M 1 M 1
0 L 0 M L M 0 L 0 1 L 0 M L M 1 L 0 M L M 0 L 1 M L M 0 L 1
EYi = µ EYi = µ + α1
EYi = µ + α m −1
Nulovou hypotézou je nulový rozdíl ve středních hodnotách první a kterékoliv následující skupiny:
α1 0 H0: M = M α 0 m −1 α1 0 H1: M ≠ M α 0 m −1 Tuto hypotézu lze otestovat prostřednictvím testové statistiky uvedené v kapitole kapitole 2.4.
Řešený praktický příklad: závislost koncentrace vitaminu D na BMI Vraťme se k příkladu, kterým jsme tuto výukovou jednotku začínali: modelujeme závislost sérové koncentrace vitaminu D (proměnná vitd) na indexu tělesné hmotnosti (proměnná bmi). K dispozici máme následující datovou tabulku se vzorkem 41 irských žen z datového souboru [vitaminD]. vitd
bmi
37,6
26,391
53,0
20,540
66,7
23,500
...
...
43,7
25,723
35,2
21,107
17,0
30,978
Proměnná vitd představuje výsledkovou proměnnou, proměnná bmi prediktor: Takto tedy vypadá vektor výsledků (Y) Y1 37,6 M = M Y 17,0 n
A takto matice plánu X 1 x11 1 26,391 M M M = M 1 x 1 30,978 n1
Nyní odhadneme parametry tohoto jednoduchého modelu prostřednictvím programu R. Nejprve si ukážeme zdlouhavější postup kopírující výpočty popsané v kapitole 2, poté okomentujeme syntaxi a výsledky volání funkce lm(), která se v programu R používá pro odhad parametrů lineárního regresního modelu. Následující rámečky uvádějí v levém sloupci kód v programu R a jeho slovní popis, v pravém sloupci pak symbolická reprezentace uvedeného postupu. 1080,462 41,000 X′X =& 1080,462 29146,261
XX <- t(X) %*% X
součin transponované matice X a původní matice X
1,056 − 0,039 ( X′X) −1 =& − 0,039 0,001
XX_inv <- solve(XX)
inverze maticového součinu
111,053 βˆ = ( X′X) −1 ( X′Y) =& − 2,392
Beta_hat <- XX_inv %*% ( t(X) %*% Y )
výpočet odhadu regresních koeficientů
Y_hat <- X %*% Beta_hat
ˆ = Xβˆ Y
predikované hodnoty výsledku
s2 <- t(Y_hat - Y) %*% (Y_hat - Y) / (41-2) reziduální součet čtverců
s2 =
Se (Yˆ − Y ) ′(Yˆ − Y ) = =& 320,632 n−k n−k
s = 17,906
c <- matrix(c(0,1),ncol=1)
sloupcový vektor T <- abs(t(c) %*% beta_hat) / (sqrt(s2) * sqrt( t(c) %*% XX_inv %*% c))
T=
c′βˆ OLS
s c′(X′X) −1 c
testové kritérium pro t-test
qt(0.975,41-2)
97,5% kvantil studentova rozdělení s 39 stupni volnosti
2,023
= 3,466
2*(1-pt(T,41-2))
0,0013
p-hodnota
Prakticky je samozřejmě odhad parametrů v programu R výrazně jednodušší – využijeme připravené funkce lm(). Základní výsledky získáme odesláním následujících funkcí: model <- lm(vitd ~ bmi, data = irlwomen) summary(model)
Dostáváme následující výsledek: Nejprve je zopakována formulace regresního modelu ve funkci lm(): Call: lm(formula = vitd ~ bmi, data = irlwomen)
Následuje základní popisná statistika vektoru reziduí: Residuals: Min 1Q Median -25.36 -12.96 -0.41
3Q 11.09
Max 52.83
Zde je uveden samotný odhad modelových parametrů, spolu s potřebnými testovými statistikami (postupně jsou uvedeny bodové odhady, směrodatné chyby těchto odhadů, hodnoty t-statistiky a příslušné p-hodnoty pro nulovou hypotézu rovnosti koeficientu 0): Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 111.0535 18.4014 6.035 4.63e-07 *** bmi -2.3924 0.6902 -3.466 0.0013 ** --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Reziduální součet čtverců a počet stupňů volnosti: Residual standard error: 17.91 on 39 degrees of freedom
Koeficient determinace a jeho adjustovaná varianta (její hodnota se snižuje s rostoucím počtem prediktorů, může tedy být vhodnější pro srovnávání predikční síly modelů): Multiple R-squared:
0.2355, Adjusted R-squared:
0.2159
A konečně F-statistika pro významnost všech prediktorů zároveň. Všimněte si, že v tomto případě (pouze jediný spojitý prediktor) je p-hodnota totožná s p-hodnotou t-testu významnosti prediktoru bmi. F-statistic: 12.02 on 1 and 39 DF,
p-value: 0.001299
Problémy k řešení Praktická práce s regresními modely: 1. Stáhněte z Internetu datový soubor [vitamin D] a odhadněte parametry modelu závislosti koncentrace vitaminu D na indexu tělesné hmotnosti (BMI). [β0 = 111,1; β1 = -2.4] 2. Odhadněte parametry takového modelu v případě, že index tělesné hmotnosti bude zadán v rámci kategorií (BMI < 20: podváha, BMI 20-24: normální hmotnost, BMI 25-29: nadváha, BMI ≥ 30: obezita). [β0 = 63,2; β1 = -7,5; β2 = -17,4; β3 = -25,8] Pochopení matematických vztahů definujících regresní model: 3. Ukažte, že odhad parametru β0 pro regresní model bez prediktorů (jen s absolutním členem) je roven výběrovému průměru hodnot výsledku. 4. Ukažte ekvivalenci vztahů pro reziduální součet čtverců 2.8 a 2.13. 5. Odvoďte testové statistiky pro regresní modely reprezentující klasický t-test a analýzu rozptylu. Srovnejte se vztahy, které znáte z klasické biostatistiky.
Literatura Použitá literatura [1] Forbelská, M.: Studijní materiály k předmětu Lineární statistické modely. Přírodovědecká fakulta Masarykovy univerzity, Brno (2009). [2] Andersen, P.K., Skovgaard, L.T.: Regression with Linear Predictors. Springer, New York (2010) Použité datové soubory • [heartdisease] dostupný z http://statweb.stanford.edu/~tibs/ElemStatLearn/ • [vitaminD] dostupný z http://staff.pubhealth.ku.dk/~linearpredictors/
