Lineární regrese v R
mgr. Andrea Kraus m.sc.,pii.d. & mgr. Vojtěch Šindlář
Tento text vznikl jako doprovodný ke cvičením kurzu M5120 Lineární statistické modely I
v rámci projektu FRMU1211/2019.
Obsah
1 Seznámení s R 2
2 Seznámení s daty 7
2.1 Data swim ..................................... 7
2.2 Data f ev...................................... 9
3 Deskriptívni statistiky 13
3.1 Jednotlivé proměnné............. ................... 14
3.1.1 Kvantitativní proměnná........................... 14
3.1.2 Kategoriální proměnná........................... 16
3.2 Vztahy mezi proměnnými.............................. 20
3.2.1 Vztah mezi dvěmi kvantitativními proměnnými.............. 23
3.2.2 Vztah mezi kvantitativní a kategoriální proměnnou............ 24
3.2.3 Vztah mezi dvěmi kategoriálními proměnnými.............. 27
4 Lineární model 35
5 Lineární model v R 47
5.1 Zadání modelu do R................................. 47
5.2 Celkové charakteristiky modelu .......................... 54
5.3 Inference pro jednotlivé koeficienty ........................ 57
5.4 Inference pro lineární kombinace koeficientů ................... 60
5.5 Inference pro více (lineárních kombinací) koeficientů............... 62
5.6 Odhady střední hodnoty a predikce na základě modelu .............. 63
6 Výběr modelu 74
6.1 Srovnání vnořených modelů pomocí testování................... 74
6.2 Srovnání modelů pomocí kritérií.......................... 76
6.3 Diagnostika modelu................................. 76
6.4 Multikolinearita................................... 88
7 Ilustrační analýza časů přeplavání jezera 91 Literatura 118
1
Kapitola 1 Seznámení s R
V této kapitole si připomeneme základní postupy při práci s €8, jelikož se v následujících kapitolách základní znalost Cl předpokládá. Čtenář, který je na práci s Ct zvyklý, může tuto kapitolu bez újmy vynechat. Čtenář, pro kterého je následující přehled naopak příliš stručný, si jistě vybere z velkého množství mnohdy i volně dostupných knih, výukových materiálů, návodů a diskusí týkajících se
Při práci s ® je rozumné ukládat si příkazy do skriptů, např. skript. R, které pak můžeme upravovat a opakovaně používat. Také objekty, se kterými v <S$ pracujeme, je rozumné si ukládat.
V d ukládáme pomocí příkazu <- .
> a <- 1    # store number 1 in a
> b <- 2    # store number 2 in b
Znak # označuje v 0$ komentář, tj. to, co následuje po #, slouží jako poznámka pro autora nebo čtenáře kódu a Q ji nevnímá jako příkaz.
Do a a b jsme si uložili čísla. Teď můžeme <S požádat, aby nám řeklo, o jaká čísla se jedná, nebo aby s nimi pracovalo.
> a    # What is in a? [1] 1
>a + b    # a + b = ? [1] 3
Očekáváme-li, že výsledek operace a + b budeme později potřebovat, raději si jej také uložíme.
3
> d <- a + b    # save the result of a + b in d
Jelikož se ve statistice hodně využívá lineární algebra a O je statistický software, je speciálně připraveno pro práci s vektory a maticemi. Vektor lze v <& vytvořit několika způsoby:
> e <- c(l,  2,   3,   4)     # a vector with elements 1,  2,   3, 4
> e
[1]   12 3 4
> f <- c(l:4)     # another way of defining the same vector
> f
[1]   12  3 4
> g <- rep(l,   4)     # a vector with element 1 repeated 4 times
> g
[i] i i i i
V kódu výše jsme poprvé použili funkce, c () a rep () . Uvnitř závorek jsme zadali požadované hodnoty jejich argumentů. Detaily o dané funkci najdeme v nápovědě, o kterou můžeme v O požádat příkazem ? . Například ? c otevře nápovědu k funkci c () .
K jednotlivým složkám vektoru přistupujeme pomocí hranatých závorek:
> f[3]     # show the third element of f [1] 3
Operace s vektory vykonává €j| po složkách (nepožádáme-li explicitně o jiný přístup):
> # componentwise operát ions
> f + g
[1]   2  3  4 5
> f + 1
[1]   2  3  4 5
> f * g
[1]   12  3 4
> 1/f
[1]   1.0000000  0.5000000  0.3333333 0.2500000
4
Všimněme si, že v kódu výše jsme pro sčítání
		
1	1	2
1	+	3
v)		w
vlastně nepotřebovali definovat vektor g. Stačilo k vektoru f přičíst obyčejnou 1 a ^ si z ní samo vytvořilo vektor vhodné délky. Všimněme si také, že jsme v kódu výše zadávali násobení a dělení vektorem, f g a 1/f, což matematicky nedává smysl. Tohoto značení se v O využívá pro násobení a dělení po složkách, jak vidíme i z výstupů výše. Toto značení má za úlohu zjednodušit uživateli zadávání operací, které statistici často potřebují.
Kdybychom, naopak, chtěli „opravdové" násobení mezi vektory, fTg nebo f gT, použili bychom v ® maticové násobení % * % (a transpozici t () ):
> t(f]
;*% g    # transpose and matrix multiplication
[, 1] [1,] 10
> f %*% t(g)
[,1] [,2] [,3] [,4]
[1, ]         1 1 1 1
[2,]         2 2 2 2
[3,]         3 3 3 3
[4,]         4 4 4 4
Všimněme si, že <§ ve skutečnosti vnímá vektory jako sloupcové a do řádku je vypisuje jenom pro úsporu místa.
Matici v £j$ vytvoříme pomocí funkce matrix () .
> A <- matrix(c(1:4),  nrow = 2,  ncol = 2)
> # 2x2 matrix with elements 1:4 filled in by columns
> A
[,1] [,2] [1,] 1 3
[2,] 2 4
> A <- matrix(c(1:4),  nrow = 2,  ncol = 2,  byrow = T)
> # 2x2 matrix with elements 1:4 filled in by rows
> A
[,1] [,2] [1,] 1 2
[2,] 3 4
5
K jejím složkám přistoupíme opět pomocí hranatých závorek:
> A[l,   1]     # element  of A at position   (1, 1) [1] 1
Regulární matici můžeme invertovat pomocí příkazu s o 1 ve () :
> solve(A)     # inversion of A
[,1] [,2] [1,] -2.0 1.0 [2,]     1.5 -0.5
Při psaní skriptů se mnohdy hodí také možnosti jakožto programovacího jazyka. Můžeme využít třeba for cyklů,
> # using a for cycle
> for   (i in 1:10) { + print(i)
+ }
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
zkonstruovat větvení programu pomocí if/else
> # using if/else
> if   (i > 5) {
+ print ("i is greater than 5")
+ } else {
+ print("i is less or equal to 5")
+ }
[1]   "i  is greater than 5"
6
anebo si vytvořit vlastní funkci pomocí příkazu f unct i on () .
> # defining a funct ion
> mysům <- function(a,  b) { + return(a + b)
+ }
> mysům(1,  2)
[1] 3
Do naší funkce mysům () vstupují dva argumenty, a a b, a funkce vrací jejich součet, a + b. Sumační funkce, pochopitelně, v Cl již existuje (jmenuje se sum() ) a nebylo ji potřeba vyrábět. Funkci mysům proto smažeme:
> rm(mysům)     # remove unneeded objects
Samotné CS obsahuje mnoho funkcí a další jsou dostupné ve specializovaných knihovnách (balíčcích). V následujících kapitolách se seznámíme s postupy, funkcemi a balíčky, které se hodí pro analýzu dat pomocí lineárních modelů.
Kapitola 2 Seznámení s daty
V této kapitole si představíme data, na kterých budeme v následujících kapitolách ilustrovat použití jednotlivých metod.
2.1   Data swim
Data swim obsahují údaje o plavcích, kteří překonali jezero Ontario mezi kanadským Torontem a ústím řeky Niagara na hranici mezi Kanadou a Spojenými státy, od prvního úspěšného pokusu v roce 1954 do konce roku 2018. Původní data jsou dostupná pod odkazem [2].
Máme-li data uložena v souboru swim.txt v podadresáři Data adresáře, ve kterém <í$ momentálně pracuje, načteme je do příkazem
> # read data from a  .txt file into R object swim
> swim <- read.table("Data/swim.txt",  header = TRUE)
Na adresár, ve kterém <5 momentálně pracuje, se můžeme © zeptat příkazem getwd () , případně jej můžeme nastavit příkazem setwd() , ve kterém do závorky zadáme cestu k požadovanému adresáři.
Příkaz read.table () načítá data z textových souborů (.txt). Pro data ve fromátu . csv lze použít příkaz read. cs v . Pro načítání dat z jiných formátů existují v ©jiné příkazy nebo balíčky. Výše jsme volbou parametru header = TRUE řekli, že první řádek souboru obsahuje názvy jednotlivých proměnných, hodnoty pak začínají od druhého řádku. Cj| počítá s formátem, kde jsou hodnoty jednotlivých proměnných ve sloupcích a hodnoty pro jednotlivé jednotky/jedince/pozorování v řádcích. Všechny volby (options) příkazu read.txt si můžeme prohlédnout pomocí 7read.txt .
Prvních několik řádků načtených dat prohlédneme příkazem
> head(swim)     # view the top of the data
Name Sex
Age Start.Day Month Year Time.min.
2.1. DATA SWIM
1 Marilyn Bell
2 John Jaremey
3 Brenda Fisher
4 Bill Sadlo
5 Jim Woods
6 Jim Woods Direction
1 SN
2 SN
3 SN
4 SN
5 SN
6 SN
Všimneme si, že jde opravdu o formát, kde každý sloupec odpovídá jedné proměnné a každý řádek jednomu jedinci (v tomto případě plavci). Pomocí příkazu
> str(swim)     # overview of the data
data.frame'		65 obs.  of    8 variables:			
$	Name	Factor w/ 54 levels  "Angela Kondrak",..	: 32	24 7 5 22 22	11
$	Sex	Factor w/ 2 levels "f","M":   12  12 2 2	1 1	11...	
$	Age	num    16 36 28 57.3 41 ...			
$	Start.Day	int    8 23 12 23 26 2  17 30  16 22 ...			
$	Month	Factor w/ 3 levels "Aug","July","Sep":	3 2 1	113 111	1
$	Year	int    1954 1956 1956 1957  1957  1961 1974	1974	1975 1976 .	
$	Time.min.	num    1255 1273 1131  1501  1115   . . .			
$	Direction	Factor w/ 3 levels "NS","NSN","SN":   3 3	3 3	3 3 3 1 3 3	
zjistíme, jak jednotlivé proměnné vnímá <3§. Z výstupu vidíme, že data vidí j ako data . frame se 65 pozorováními o 8 proměnných. Jedná se v podstatě o matici, u data. frame je ale povoleno, aby jednotlivé sloupce byly různých typů. Z výstupu str (swim) vidíme, že proměnné Name, Sex, Month a Direction vidíjako diskrétní (načetlo jejako/afctory s jistými úrovněmi), zatím co zbylé proměnné vidí jako spojité (načetlo je jako čísla). Základní přehled o datech získáme příkazem
> summary(swim)     # descriptive statistics
Name		Sex	Age			Start	Day
Vicki Keith	4	F:39	Min.	14	. 05	Min.	1
Colleen Shields	3	M: 26	1st Qu.	21	.00	1st Qu.	8
Kim Middleton	3		Median	28	.00	Median	13
Jim Woods	2		Mean	31	. 03	Mean	14
John Scott	2		3rd Qu.	38	.00	3rd Qu.	19
F 16.00000 M 36.00000 F 28.00000 M 57.31781 M 41.00000 M 45.00000
8 Sep 1954
23 July 1956
12 Aug 1956
23 Aug 1957
26 Aug 1957
2 Sep 1961
1255 1273 1131 1501 1115 1027
2.2. datafev
9
Kim Lumsdon :   2 Max.       :66.57      Max. :31
(Other) :49
Month Year Time.min. Direction
Aug	50	Min.	1954	Min.	829	NS	5
July	6	1st Qu.	1979	1st Qu.	1027	NSN	1
Sep	9	Median	1993	Median	1199	SN	59
		Mean	1992	Mean	1279		
		3rd Qu.	2007	3rd Qu.	1413		
		Max.	2018	Max.	3370		
Všimněme si, že €jt pro každou proměnnou vybere deskriptívni statistiky podle toho, zda ji vnímá jako diskrétní, nebo spojitou.
Na první pohled není vidět žádné problémy s načítáním datové sady. Ty by nastaly zejména v případě, načetlo-li by    numerickou proměnnou jako znak nebo faktor.
2.2   Data f e v
Datový soubor f ev obsahuje výběr z údajů o dětech a mladých lidech, kteří se účastnili studie Childhood Respiratory Desease Study v roce 1980 ve Spojených Státech, která zkoumala vývoj plicní funkce u dětí. Data i jejich popis lze najít pod odkazem [1]. Klíčová proměnná Forced Expiratory Volume (FEV) měří objem vzduchu vydechnutého za první sekundu intenzivního výdechu. Kromě této proměnné data obsahují údaje o věku, výšce, pohlaví a o tom, zda se jedná o kuřáka nebo ne.
Data načteme příkazem
> fev <- read.table("Data/fev.txt",  header = TRUE) Začátek databáze prohlédneme příkazem
> head(fev)
	ID	Age		FEV	Height	Sex	Smoker
1	301	9	1	.708	57 . 0	Female	Non
2	451	8	1	. 724	67 . 5	Female	Non
3	501	7	1	. 720	54 . 5	Female	Non
4	642	9	1	. 558	53.0	Male	Non
5	901	9	1	.895	57 . 0	Male	Non
6	1701	8	2	.336	61 . 0	Female	Non
Základní přehled o databázi získáme příkazem
2.2. DATA FEV
10
> summary(fev)
ID Age
Min.	201	Min.	3.	000
Ist Qu.	15811	1st Qu.	8.	000
Median	36071	Median	10.	000
Mean	37170	Mean	9.	931
3rd Qu.	53639	3rd Qu.	12 .	000
Max.	90001	Max.	19.	000
Sex		Smoker		
Female:318 Current: 65 Male    :336      Non :589
FEV Height
Min.	0.	791	Min.	46	.00
1st Qu.	1 .	981	1st Qu.	57	.00
Median	2 .	547	Median	61	.50
Mean	2 .	637	Mean	61	. 14
3rd Qu.	3.	119	3rd Qu.	65	.50
Max.	5 .	793	Max.	74	.00
Ani na těchto datech nejsou zřejmé žádné problémy s načítáním. Nyní si data prohlédneme podrobněji.
Začneme s rozsahem dat. Funkce dim () vrátí rozměr dat.
> dim(fev)     # dimensions of the data frame [1]   654 6
Jedná se teda o rozsáhlejší data (6 proměnných o 654 jedincích).
Někdy je praktické přistoupit jen k části dat. Konkrétní řádky nebo sloupce vyžádáme následujícím způsobem.
> fev[l:10,   ]     # first ten rows and all columns of the data
	ID	Age		FEV	Height		Sex	Smoker
1	301	9	1	.708	57	. 0	Female	Non
2	451	8	1	. 724	67	. 5	Female	Non
3	501	7	1	. 720	54	. 5	Female	Non
4	642	9	1	. 558	53	. 0	Male	Non
5	901	9	1	.895	57	. 0	Male	Non
6	1701	8	2	.336	61	. 0	Female	Non
7	1752	6	1	. 919	58	. 0	Female	Non
8	1753	6	1	.415	56	. 0	Female	Non
9	1901	8	1	. 987	58	. 5	Female	Non
10	1951	9	1	. 942	60	. 0	Female	Non
> fev[c(l,   3,   5),   ]     # rows 1,   3,   5 and all columns
2.2. DATAFEV
11
	ID Age			FEV	Height Sex	Smoker
1	301	9	1.	708	57.0 Female	Non
3	501	7	1.	720	54.5 Female	Non
5	901	9	1.	895	57.0 Male	Non
>	fev [1	:10,		1:2]	# first ten	rows and
	ID	Age				
1	301		9			
2	451		8			
3	501		7			
4	642		9			
5	901		9			
6	1701		8			
7	1752		6			
8	1753		6			
9	1901		8			
10	1951		9			
O konkrétní sloupec lze požádat také jeho jménem a následně s ním pracovat jako s vektorem.
> range(fev$FEV)     # minimum and maximum [1]   0.791 5.793
> mean(fev$FEV)     # mean [1] 2.63678
> sd(fev$FEV)     # standard deviation [1] 0.8670591
Rovněž je možné požádat jenom o část dat zadanou jejími vlastnostmi.
> # descriptive statistics for males
> summary(fev[fev$Sex == "Male", ])
ID
Age
FEV
Height
Min.
1st Qu Median Mean
201 15317 34171 36233
Min.
1st Qu Median Mean
3.00 8.00
10.00
10.01
Min.
1st Qu Median Mean
0.796 2.009 2. 606 2 .812
Min.
1st Qu Median Mean
47 .00 57 .00 62 .00 62 . 03
2.2. DATAFEV
12
3rd Qu.:52286 3rd Qu.:12.00 3rd Qu.:3.535 3rd Qu.:67.50 Max.       :90001      Max.       :19.00      Max.       :5.793      Max. :74.00
Sex Smoker Female:     0      Current:  2 6 Male    :336      Non :310
> # descriptive statistics for height in females
> summary(fev$Height[fev$Sex == "Female"])
Min.   1st Qu.    Median        Mean 3rd Qu. Max. 46.00      57.50      61.00      60.21      63.50 71.00
Proměnná ID obsahuje identifikační čísla jedinců v rámci studie, ze které data pocházejí. V našich datech, která jsou podmnožinou původních dat ze studie, máme od každého jedince jen jedno pozorování, jak snadno ověříme třeba sečtením výstupů z funkce duplicated () . Tato funkce přiřadí hodnotu TRUE, tj. 1, těm prvkům vektoru, které se již vyskytly na předchozích pozicích, a FALSE, tj. 0, všem ostatním.
> sum(duplicated(fev$ID))
[1] 0
Identifikační čísla jedinců dále nebudeme potřebovat, proto se nám bude hodit tuto proměnnou z dat vymazat. Na výstupech výše jsme si možná také uvědomili, že výška je zadaná v palcích. Abychom se v ní lépe orientovali, převedeme ji na centimetry.
> # remove the first column of the data
> fev <- fev[, -1]
> # replace height in inches by height in centimeters
> fev$Height <- fev$Height * 2.54
Výsledná data uložíme, abychom je příště nemuseli opět upravovat.
> # save the modified data for later use
> save(fev,   file = "fev.RData")
Upravená data posléze do eft načteme příkazem load("fev.RData") . Pomocí funkce save () je možné do souboru typu . RDat a uložit více objektů současně. Ty se pak příkazem load ()  načtou přímo do prostředí <H, kde budou figurovat s takovými názvy a vlastnostmi, jaké měly při ukládání. Do souboru typu . t xt ukládáme data pomocí příkazu write.table ()
Kapitola 3
Deskriptívni statistiky
Když začínáme pracovat s nějakou datovou sadou, je rozumné si data nejdříve prohlédnout: zjistit, jaké obsahují proměnné, jakých hodnot tyto proměnné nabývají a jak jsou v datech vzájemně propojeny. Takový předběžný pohled nám umožní najít a odstranit případné chyby v datech, utvořit si představu o vhodném postupu samotné analýzy a předjímat případné problémy. Provést jej můžeme s využitím různých metod deskriptívni statistiky. V této fázi tedy nebudeme využívat žádných modelů, jenom vhodně zvolených charakteristik dat. Ty mohou být numerické nebo grafické a jejich výběr závisí od typu proměnné nebo proměnných, na které se chceme zaměřit.
V matematické statistice jsme se zabývali zejmnéna diskrétními a spojitými náhodnými veličinami, které se mezi sebou liší především počtem hodnot, kterých mohou nabývat (spočetně, nebo nespočetně mnoho). Jelikož i realizace spojitých veličin obvykle měříme na diskrétní stupnici (dané například rozlišením měřící technologie), v aplikované statistice se proměnné spíše dělí na kvalitativní (kategoriální) a kvantitativní (číselné). Kategoriální proměnné pak mohou být nominální (kategorie nemají uspořádání), anebo ordinální (kategorie lze uspořádat). V <í8 kódujeme kategoriální proměnné třídou factor, případně ordered factor a kvantitativní proměnné třídou numeric, případně integer. Jsou-li data na vstupu dobře připravena, O často při jejich načítání samo přiřadí proměnným vhodné třídy. Je-li kategoriální proměnná v datech kódována číselně, lze jí v ď změnit na f act or pomocí funkce f act or () .
Třídy jednotlivých proměnných v konkrétním data . frame zjistíme například z prvního sloupce výstupu příkazu
> # types of variables in the first column of the output below
> str(fev)
'data.frame':   654 obs.  of    5 variables: $ Age      :  int 9879986689... $ FEV      :  num    1.71  1.72  1.72  1.56 1.9 ... $ Height:  num    145 171 138  135 145 ...
$ Sex : Factor w/ 2 levels "Female","Male": 1112211111 ... $ Smoker:  Factor w/ 2 levels  "Current","Non":   2222222222 ...
Proměnná Age je tedy typu integer, proměnné FEV a Height jsou třídy numeric aproměnné Sex a Smoker jsou třídy factor, jak se nám hodí. Na třídu konkrétní proměnné se můžeme
3.1. JEDNOTLIVÉ PROMĚNNÉ
14
zeptat také příkazem class () .
> # type of variable Sex in the fev data
> class(fev$Sex)
[1] "factor"
Jsou-li proměnné v ® načteny se správnými třídami, vrátí nám příkaz summary () pro každou z nich vhodné číselné deskriptívni statistiky:
> # basic descriptive statistics chosen by the type of variable
> summary(fev)
Age FEV Height Sex
Min.	3.	000	Min.	0.	791	Min.	116 .	8	Female:318
1st Qu.	8.	000	1st Qu.	1.	981	1st Qu.	144 .	8	Male :336
Median	10.	000	Median	2 .	547	Median	156 .	2	
Mean	9.	931	Mean	2 .	637	Mean	155 .	3	
3rd Qu.	12 .	000	3rd Qu.	3.	119	3rd Qu.	166 .	4	
Max.	19.	000	Max.	5.	793	Max.	188 .	0	
Smoker Current: 65 Non :589
Nyní se na vhodné číselné i grafické charakteristiky podle typu proměnné podíváme podrobněji.
3.1 Jednotlivé proměnné 3.1.1   Kvantitativní proměnná
Rozložení kvantitativní proměnné můžeme číselně charakterizovat pomocí měr polohy nebo měr variability. Mezi oblíbené míry polohy patří průměr, a různé kvantily, zejména medián, maximum, minimum, dolní a horní kvartil. To jsou také charakteristiky na výstupu z funkce summary () aplikované na kvantitativní proměnnou. Mezi populární míry charakteristiky patří (mezi)kvartilové rozpětí, rozptyl a směrodatná odchylka. U všech uvedených charakteristik jde, pochopitelně, o výběrovou verzi.
3.1. JEDNOTLIVÉ PROMĚNNÉ
15
> IQR(fev$FEV)     # interquartile range [1] 1.1375
> var(fev$FEV)     # variance [1] 0.7517915
> sd(fev$FEV)     # standard deviation [1] 0.8670591
Graficky lze rozložení kvantitativní proměnné znázornit pomocí histogramu nebo pomocí krabicového diagramu.
> # histogram
> hist(fev$FEV,   freq = FALSE,
+ main =  "FEV",   xlab = "FEV [1]",
+ col = "lightgreen",  border = "palegreen3")
FEV
o
CO
o
CD
Q
CvJ
o
o
o o
1
2
3
4
5
6
FEV [I]
3.1. JEDNOTLIVÉ PROMĚNNÉ
16
> # boxplot
> boxplot(fev$FEV,  horizontal = TRUE,
+ main =  "FEV",     xlab =  "FEV [1]",
+ col = "lightgreen",  border = "palegreen3",  pen = 19)
FEV
-1-1-1-1-1-
1 2 3 4 5
FEV [I]
Rozdělení proměnné FEV je symetrické, až na několik vzdálenějších (možná odlehlých) větších hodnot.
3.1.2   Kategoriální proměnná
Rozložení kategoriální proměnné lze plně popsat i číselnými charakteristikami, konkrétně procentuálním rozložením dat mezi jednotlivými kategoriemi. Užitečný může být i počet pozorování v jednotlivých kategoriích. Ten také najdeme na výstupu z funkce summary () aplikované na kategoriální proměnnou. Pro ordinální kategoriální proměnné je informativní také kumulativní procentuální zastoupení v uspořádaných kategoriích.
> # nominal variable: percentages per category
> prop.table(table(fev$Sex) )
Female Male 0.4862385 0.5137615
3.1. JEDNOTLIVÉ PROMĚNNÉ
17
> # ordinal variable
> summary(swim$Month)
Aug July Sep 50        6 9
> # first  the right order
> swim$Month <- factor(swim$Month,
+ levels = levels(swim$Month) [c (2,   1, 3)],
+ ordered = TRUE)
> summary(swim$Month)
July    Aug Sep 6      50 9
> prop.table(table(swim$Month))   # percentages per category
July Aug Sep
0.09230769 0.76923077 0.13846154
> cumsum(prop.table(table(swim$Month)))   # cumulative percentages
July Aug Sep
0.09230769 0.86153846 1.00000000
Pro grafické znázornění rozložení kategoriální proměnné lze využít sloupcový nebo koláčový graf.
> # bar plot for Sex
> barplot(
+      100*prop.table(table(fev$Sex) ) , +      col = c("lightpink",   "cornflowerblue"), +      border = c("lightpink2",   "dodgerblue2") + )
3.1. JEDNOTLIVÉ PROMĚNNÉ
18
> # bar plot for Sex
> barplot(
+      100*matrix(prop.table(table(fev$Sex)),  nrow = 2,  ncol = 1),
+      ylab = "(Cumulative) percentage",
+      col = c("lightpink",   "cornflowerblue"),
+      border = c("lightpink2",   "dodgerblue2")
+ )
> text(x =0.7,  y = 4,   labels = "Female",  cex = 1.1,  pos = 3)
> text(x = 0.7,  y = 55,   labels = "Male",   cex = 1.1,  pos = 3)
3.1. jednotlivé proměnné
19
> # pie chart for Sex
> pie (
+ summary(fev$Sex),
+ col = c("lightpink",   "cornflowerblue"),
+ border = c("lightpink2",   "dodgerblue2")
+ )
3.2. VZTAHY MEZI PROMĚNNÝMI
20
Data obsahují zhruba stejný počet děvčat a chlapců.
3.2   Vztahy mezi proměnnými
Číselnné charakteristiky pro všechny proměnné v datech získáme pomocí příkazu summary () . Příkaz pairs() vrátí grafické znázornění vztahů mezi proměnnými.
> # relationships between variables
> pairs(fev,   col = "lightgreen",  pen = 19)
3.2. VZTAHY MEZI PROMĚNNÝMI
21
Age
1 2 3 4 5 J_I_'_ '_L_
1.0       1.4 1.8 i_I_I_I_I_1
FEV
-1
Height
n      i r
5       10 15
n—i—i—i—i—i—r-
120   140   160 180
Sex
Smoker
i—i—i—i—i—r ^
1.0       1.4 1.8
Tyto grafy jsou smysluplné především pro kvantitativní proměnné.
> pairs(fev[,  c(2:4)],   col = "lightgreen",  pch = 19)
3.2. VZTAHY MEZI PROMĚNNÝMI 22
o
cm
FEV
120
140 160 _l_I_1_
180
Height
Sex
- lo
- -st
- co
- cm
-1-1-1-1-
1.0   1.2   1.4   1.6   1.8 2.0
Funkce ggpairs () z knihovny GGally vybírá vhodné grafy pro kvantitativní i kate-goriální proměnné.
> # relationships between variables   (quantitative and categorical)
> library("GGally")
> ggpairs(fev)
3.2. VZTAHY MEZI PROMĚNNÝMI
23
Age FEV Height Sex Smoker
Nyní se zaměříme na vhodné grafické i numerické popisy vztahů mezi proměnnými podle jejich typů.
3.2.1   Vztah mezi dvěmi kvantitativními proměnnými
Graficky můžema vztah mezi dvěmi kvantitativními proměnnými zobrazit pomocí klasického grafu (funkce plot () ). Funkce lowess () neparametrický (bez parametrického modelu) odhadne závislost mezi proměnnými. Tu pak do grafu přidáme pomocí funkce 1 ines () .
> plot(
+      fev$Height~fev$Age,
+      main = "Height by age",
+      ylab = "Height   [cm]",  xlab = "Age [y]", +      pch = 19,   col = "lightgreen"
+ )
>
> lineš(lowess(fev$Height~fev$Age),   lwd = 3,   col = "palegreen3")
3.2. VZTAHY MEZI PROMĚNNÝMI
24
Height by age
i-1-r
5 10 15
Age [y]
Číselným vyjádřením (lineární) závislosti mezi proměnnými je korelace.
> cor(fev$Height, fev$Age) [1] 0.7919436
Těsnou závislost mezi věkem a výškou jsme očekávali. Také jsme očekávali postupné spo-malování růstu výšky s věkem až po jeho pozvolné zastavení.
3.2.2   Vztah mezi kvantitativní a kategoriální proměnnou
Rozložení kvantitativní proměnné v závislosti na proměnné kategoriální můžeme studovat porovnáním podmíněného rozložení kvantitativní proměnné při daných úrovních kategoriální proměnné. To můžeme popsat pomocí charakteristik z části 3.1.1 jednotlivě vyhodnocených na datech rozdělených podle úrovní kategoriální proměnné.
> # descriptive statistics for Height
> # separately for girls and boys
> summary(fev$Height[fev$Sex == "Female"])
Min. 1st Qu. Median Mean 3rd Qu. Max. 116.8      146.1      154.9      152.9      161.3 180.3
3.2. VZTAHY MEZI PROMĚNNÝMI
25
> summary(fev$Height[fev$Sex == "Male"])
Min.   1st Qu.    Median        Mean 3rd Qu. Max. 119.4      144.8      157.5      157.5      171.4 188.0
> # boxplots for Height
> # separately for girls and boys
> boxplot(
+      fev$Height~fev$Sex,
+      main = "Height by Sex",     ylab = "Height [cm]", +      col = c("lightpink",   "cornflowerblue"), +      border = c("lightpink2",   "dodgerblue2"), +      pch =19 + )
Height by Sex
E
o
o
00
o
CD
O)
'<u o
o
C\J
Female
Male
fev$Sex
> par(mfrow = c(l,  2))   # divide the plot
> # into left and right panels
> # draw in the left panel
> hist(
+      fev$Height[fev$Sex ==  "Female"],   freq = FALSE,
+      xlim = c(min(fev$Height)   - 1,  max(fev$Height)   + 1),
3.2. VZTAHY MEZI PROMĚNNÝMI
26
+      ylim = c(0, 0.04),
+      main = "Height:  girls",   xlab = "Height [cm]", +      col =  "lightpink",  border = "lightpink2" + )
> # draw in the right panel
> hist(
+      fev$Height[fev$Sex ==  "Male"],   freq = FALSE,
+      xlim = c(min(fev$Height)   - 1,  max(fev$Height)   + 1),
+      ylim = c(0, 0.04),
+      main = "Height: boys",  xlab = "Height [cm]", +      col = "cornflowerblue",  border = "dodgerblue2" + )
Height: girls
Height: boys
CD
Q
o
CO
o
CvJ
o
o o
o o
á
^ r
I-1-1-1
120   140   160 180 Height [cm]
	
	o
	o
	CO
	o
	o
	
Hsu	OJ o
CD	o
Q	
	T—
	o
	o
	o
	o
	o
I-1-1-1
120   140   160 180 Height [cm]
Zatímco krabicové diagramy ® vykreslilo do společného grafu, histogramy jsme museli nakreslit do dvou sousedících grafů. Aby byly porovnatelné, nastavili jsme jim stejné rozsahy os (možnosti xlim () a ylim ()). Čtenáři seznámeni se syntaxí knihovny ggplot2 mohou histogram podle skupin vykreslit s její pomocí.
> # boxplots for Height
> # separately for girls and boys using ggplot
> library (ggplot2)
> ggplot(data = fev,  aes(x = Height,   color = Sex,
3.2. VZTAHY MEZI PROMĚNNÝMI
27
+ stat(density))) +
+ geom_histogram(position = "identity",
+ aes(fill = Sex,  color = Sex),
+ binwidth=5,   alpha =0.5) +
+ ylab("Density" )
0.03-
"S 0.02-
CD
Q
0.01
0.00-
Sex
Female Male
125
150
Height
175
Chlapci j sou podle očekávání o něco vyšší než dívky. Jejich výšky j sou v porovnání s děvčaty o něco více rozptýlené.
3.2.3   Vztah mezi dvěmi kategoriálními proměnnými
Vzájemný vztah mezi dvěma kategoriálními proměnnými také můžeme popsat pomocí podmíněných rozdělení. Za tímto účelem využijeme funkce table () , která vrátí počty jedinců pro jednotlivé kombinace kategorií. Na výsledek můžeme dále aplikovat funkci pr op . table () , která počty přepočítá na proporce. V základní formě se procenta v celé tabulce sčítají na 1 (popisují sdružené rozdělení). Pomocí voleb margin = 1 a margin = 2 se na 1 sčítají procenta v řádcích nebo sloupcích (popisují marginální rozdělení).
> # joint distribution of Smoking and Sex
> table(fev$Smoker, fev$Sex)
3.2. VZTAHY MEZI PROMĚNNÝMI
28
Non
Current
Female Male 39 26 279 310
> # conditional  distribution of Sex for given Smoking status
> prop.table(table(fev$Smoker,   fev$Sex),  margin = 1)
> # conditional  distribution of Smoking for given Sex
> prop.table(table(fev$Smoker,   fev$Sex) ,  margin = 2)
Sdružené rozdělení můžeme vizualizovat pomocí funkce mo s a i cp 1 o t () .
> # mosaic plot for Sex and Smoking
> mosaicplot(
+      table (fev$Sex,   fev$Smoker),
+      main = "Smoking and Sex",
+      col = c("cornflowerblue", "lightpink")
+ )
Female Male Current 0.6000000 0.4000000 Non 0.4736842 0.5263158
Female Male Current 0.12264151 0.07738095 Non 0.87735849 0.92261905
3.2. VZTAHY MEZI PROMĚNNÝMI
29
Smoking and Sex
Female
o
Male
Podmíněné rozdělení můžeme také vizualizovat pomocí sloupcového grafu rozděleného podle jednotlivých úrovní kategoriální proměnné.
> # bar plot for Smoking by Sex
> barplot(
+      height = prop.table(table(fev$Smoker,   fev$Sex),  margin = 2),
+      beside = TRUE,
+      main = "Smoking by Sex",
+      ylab = "Percentage",  xlab = "Sex",
+      col = c("cornflowerblue" ,   "lightpink"),
+      border = c("dodgerblue2",   "lightpink2")
+ )
> legend(
+      x = 3.3,  y = 0.8,
+      legend = c("Current", "Non"),
+      col = c("cornflowerblue", "lightpink"),
+      pch = 15
+ )
3.2. VZTAHY MEZI PROMĚNNÝMI
30
> # bar plot for Smoking by Sex
> barplot(
+      height = prop.table(table(fev$Smoker,   fev$Sex),  margin = 2),
+      beside = FALSE,
+      main = "Smoking by Sex",
+      ylab = "Percentage",  xlab = "Sex",
+      col = c("cornflowerblue" ,   "lightpink"),
+      border = c("dodgerblue2",   "lightpink2")
+ )
> text(
+      x = 1.2 * c(0:2)   + 0.7,  y = c(-.01, -.03),
+      labels = "Current",
+      cex = c(l.l,   1),  pos = 3
+ )
> text(
+      x = 1.2  * c(0:2)   +  0.7,   y = 0.5,
+      labels = "Non",
+      cex = 1.1,  pos = 3
+ )
3.2. VZTAHY MEZI PROMĚNNÝMI
31
Smoking by Sex
o
Female Male
Sex
Všimněme si, že představu o kouření mezi chlapci a děvčaty jsme si mohli udělat i na základě mozaikového grafu výše. Mezi děvčaty najdeme kuřačky o něco častěji než kuřáky mezi chlapci.
Rozdělení pohlaví mezi kuřáky a nekuřáky získáme výměnou pořadí proměnných v kódech uvedených výše.
> # mosaic plot for Sex and Smoking
> mosaicplot(
+      table(fev$Smoker,   fev$Sex),
+      main = "Sex and Smoking",
+      col = c("lightpink",   "cornflowerblue")
+ )
3.2. VZTAHY MEZI PROMĚNNÝMI 32
Sex and Smoking
Current
Non
> # bar plot for Sex by Smoking
> barplot(
+      height = prop.table(table(fev$Sex,   fev$Smoker),  margin = 2),
+      beside = TRUE,
+      main = "Sex by Smoking",
+      ylab = "Percentage",  xlab = "Smoking",
+      col = c("lightpink",   "cornflowerblue"),
+      border = c("lightpink2",   "dodgerblue2")
+ )
> legend(
+ x = 3,  y = 0 . 6,
+ legend = c("Female", "Male"),
+ col = c("lightpink",   "cornflowerblue"),
+ pch = 15
+ )
3.2. VZTAHY MEZI PROMĚNNÝMI
33
Sex by Smoking
O)
CÖ -t—'
CD
ü s_
CD
Q_
CD Ö
Ö
C\J Ö
O Ö
i Female		
i Male		
		
Current
Non
Smoking
> # bar plot for Sex by Smoking
> barplot( height = prop.table(table(fev$Sex, beside = FALSE, main = "Sex by Smoking", ylab = "Percentage",  xlab = "Smoking", col = c("lightpink",   "cornflowerblue") , border = c("lightpink2",   "dodgerblue2"'
fev$Smoker),  margin = 2),
text (
x = 1.2 * c(0:2)   + 0.7, labels = "Female", cex = 1.1,  pos = 3 )
text (
x = 1.2 * c(0:2)   + 0.7, labels = "Male", cex = 1.1,  pos = 3
0.2,
0 . 75,
3.2. VZTAHY MEZI PROMĚNNÝMI
34
Smoking
Mezi kuřáky j sou o něco více zastoupena děvčata než chlapci, mezi nekuřáky je tomu obráceně.
Kapitola 4 Lineární model
Lineární model popisuje závislost mezi vysvětlovanou proměnnou Y (říkáme jí také závisle proměnná nebo odezva, anglicky outcome, response anebo také dependent variable) a vysvětlujícími proměnnými xi,... ,Xk (říkáme jim také nezávisle proměnně nebo prediktory, anglicky co-variates, predictors anebo také idependent variables nebo explanatory variables). V lineárním modelu je jejich závislost ve tvaru
Yi = fa + + ... + fikxiik + e,t,   i = l,...,n, (4.1)
kde (30,... ,/3k jsou regresní koeficienty (anglicky regression coefficients) a Si jsou náhodné chyby (anglicky random errors). Jak již napovídá terminologie, náhodné chyby považujeme za náhodné veličiny. Naproti tomu prediktory považujeme za nenáhodné, respektive jsou-li prediktory náhodné, studujeme závislost odezvy Y na prediktorech X\,..., Xk podmíněně při jejich pozorovaných hodnotách xi,..., xk, a tuto závislost pak popisujeme rovnicí (4.1). Odezva Y je v obou přístupech náhodná veličina, zatímco koeficienty [30,..., (3k jsou nenáhodné.
Má-li se jednat o lineární model, pak vše, co je v závislosti odezvy na prediktorech systematické, musí být popsáno lineárním prediktorem i]i = (30 + flix^i + ... + fikxi,k, tj- El^ = rji, resp. Ee,i = 0, i = 1,... ,n. Navíc předpokládáme, že odezvy pro jednotlivá pozorování jsou nekorelované, C or (Yi, Yj) = 0 pro i ^j,i,j = 1,... ,n, a mají stejný rozptyl Var Yi = a2, i = 1,..., n. Jelikož Var Yi = Var e i, lze rovnost rozptylů chápat také tak, že odezvy pro jednotlivá pozorování kolem svých středních hodnot kolísají se stejnou přesností. Všechny předpoklady lineánŕho modelu můžeme formulovat v řeči náhodných chyb e i, od kterých požadujeme stejné rozdělení, nulovou střední hodnotu, nekorelováno st a stejný rozptyl. Někdy navíc požadujeme normální rozdělení a jelikož v normálním rozdělení nekorelovanost implikuje nezávislost, předpoklady v takovém případě můžeme zapsat ve tvaru Si ~ N(0, a2), i G {1,..., n}, kde skratka iid pochází z anglického independent, identically distributed: nezávislé, stejně rozdělené.
Lineární prediktor rji = (30 + fti^i + ... + je funkce lineární jak v prediktorech
tak v koeficientech. Linearita v prediktorech a linearita v koeficientech ale v lineárním modelu nehrají stejnou roli. Uvažujme například model
Yi = (30 + fcxi^ + (32Xíj2 + Ei,   i = 1,..., n.
36
Jde o speciální případ modelu (4.1), kde k můžeme uvažovat i model
2, a prediktory jsou x\ a x2. Se stejnými daty
nebo
anebo
Yi = fio + Pixli + (32xi:2 + Si, i Yi = fa + fixXi^ + fcxh + ÄÄ,2 + £i
Yi = (30 + /3i— +       i = l,...,n.
^,2
(4.2)
(4.3)
(4.4)
Ve všech případech se jedná o lineární model. Model (4.2) je speciálním případem modelu (4.1) s k = 2 a prediktory x\ ax2, Model (4.3) je speciálním případem modelu (4.1) s k = 3 a prediktory x1,x"lax2,a model (4.4) je speciálním případem modelu (4.1) s A; = 1 a prediktorem xi/x2. V dalším textu nebudeme modely vždy nutně přepisovat do tvaru (4.1). Dojde-li k transformaci prediktoru nebo prediktorů, například tím, že prediktor x\ nahradíme nebo doplníme prediktorem x\, anebo prediktory x\ a x2 nahradíme prediktorem xi/x2, ponecháme modely ve tvarech (4.2), (4.3), anebo (4.4). Také nutně nepřestaneme původní prediktory x\ a x2 nazývat prediktory, ani když se v modelu vyskytovat nebudou (vždy si můžeme představit, že se v modelu vyskytují, ale příslušný koeficient (3 je nulový). Linearita prediktoru v koeficientech naopak zachována být musí. Například model
Yi = /30 + /32x^\ + e,,, i
l,...,n,
už není lineárním modelem, protože jej nelze přepsat do tvaru (4.1).
Zápis (4.1) definuje lineami model po složkách. Mnohdy se hodí využít maticového zápisu
kde
				Y = X/3	
[Y1\		(l		x1)2   ...	Xl,k\
Y2		1	x2,i	x2)2   ■ ■ ■	X2,k
		\1	Xn,l	xn,2    ■ ■ ■	Xn,k /
(4.5)
/AA
£2
\£nj
Matici X v zápisu (4.5) říkáme regresní matice, matice modelu nebo také matice plánu, anglicky design matrix nebo model matrix.
Nyní se podíváme na několik lineárních modelů pro data f e v ze 2. kapitoly. Nejprve si ale připomeneme, jak data vypadají.
> head(fev)     # first rows of the data
37
	Age		FEV	Height	Sex	Smoker
1	9	1	.708	144.78	Female	Non
2	8	1	. 724	171.45	Female	Non
3	7	1	. 720	138.43	Female	Non
4	9	1	. 558	134.62	Male	Non
5	9	1	.895	144 .78	Male	Non
6	8	2	.336	154.94	Female	Non
> tail(fevl
# last rows of the data
	Age		FEV	Height	Sex	Smoker
649	16	4	. 872	182.88	Male	Current
650	16	4	.270	170 .18	Male	Current
651	15	3	. 727	172.72	Male	Current
652	18	2	. 853	152.40	Female	Non
653	16	2	.795	160 . 02	Female	Current
654	15	3	.211	168.91	Female	Non
,n
(1.) FEVj = & + Ä x Height, + eh   i = 1, ,
Odezvou je v tomto modelu indikátor objemu plic FEV a jediným prediktorem je výška Height. Máme tedy k = 1. Pomocí vektorů a matic bychom model zapsali následovně.
/1.708\ 1.724
/l 1
2.795
\3.2iiy
Pro lineárni prediktor platí rj,
144.78\ 171.45
^2
fii x Height^, střední hodnota FEV
1    160.02       v-/ £653 V 1    168.91/ \e654J
e(FEV4) = A,
podle modelu tedy lineárně roste s výškou. Koeficient fíi vyjadřuje, o kolik litrů je střední hodnota FEV vyšší u dítěte, které je o 1 cm vyšší. Obecně netvrdíme, že nárůst střední hodnoty FEV je způsobený rozdílem ve výšce: to jenom na základě modelu tvrdit nemůžeme. Proto raději mluvíme o asociaci: u vyšších dětí pozorujeme větší plíce. Koeficient [30 v modelu reprezentuje střední hodnotu FEV pro dítě s nulovou výškou: v tomto modelu tedy rozumnou interpretaci nemá. Vztah mezi střední hodnotou FEV a výškou popsaný modelem by mohl vypadat jako na obrázku níže.
# Coefficient beta_l represents the increase in expected value
# of FEV associated with a 1 cm increase in Height.
# We talk about associations rather than causal dependence:
# the increase in expected value of FEV is  *associated with*
# an increase in Height,  not  *caused by* it.
# beta_0 represents the expected value of FEV for a child
# with Height of 0 cm,   which makes little sense.
38
(2.) FEVí = /30 + Ä x Height^ + /32 x Height? + e,t,   i = l,...,n
I v tomto modelu máme odezvu FEV, ale nyní je k = 2. Technicky jsou prediktory výška Height a kvadrát výšky Height2, ve skutečnosti ale jde o funkce jediného prediktoru: výšky. Maticový zápis je
/1.708\		(l	144.78 144.782\		
1.724		1	171.45 171.452	/A>\	
	=			x   Ä +	
2.795		1	160.02 160.022	W	^653
\3.2iy			168.91 168.917		\£654/
Pro lineární prediktor platí r\i = e(FEVj) = /30+/3i x Heightj+/32 x Height2. Střední hodnota FEV podle modelu tedy kvadraticky roste s výškou. Interpretaci jednotlivých koeficientů je v takovém modelu vhodnější nahradit obrázkem anebo výpočtem či porovnáním střední hodnoty FEV pro několik zajímavých hodnot výšky Height.
# In a model like this,   it is more practical to replace
# the interpretation of individual coefficients by a picture
# or by computing and comparing expected values of FEV
# for several interesting values of Age.
39
(3.) FEVj = (30 + Ä x Sexj + eť,   i = 1,..., n
V tomto modelu máme odezvu FEV a k = 1. Prediktor Sex je tentokrát kategoriální proměnná o dvou úrovních: chlapec a dívka. Tyto úrovně budeme číselně kódovat jako 0 a 1. Nezvolíme-li si jinak, «5? přiřadí menší číslo úrovni, která je první v abecedním pořadí.
V našem případě půjde o dívku (female). Maticový zápis modelu je
/1.708\      / 1     0 \ / e1 \
1.724 1     0        fQ N e2
2.795
£653
\3-211/      \l     OJ \e654J
Pro lineární prediktor platí rji = e(FEVj) = [30 + fii x Sexj. Koeficient /30 reprezentuje střední hodnotu FEV pro dívky, zatímco koeficient fii vyjadřuje, o kolik litrů vyšší je střední hodnota FEV pro chlapce oproti dívkám. Střední hodnotu FEV pro chlapce můžeme spočítat jako A) + ft.
# Coefficient beta_0 represents the expected value of FEV
# for females   (coded as 0);  coefficient beta_l represents
# the difference between expected value of FEV for males (coded
# as 1)   and females   (coded as 0).
# The default order of coding is alphabetical   (F comes before M) .
# The expected value of FEV for males is beta_0 + beta_l.
40
(4.) FEVí = (30 + Ä x Height^ + (32 x Sex^ + eh   i = l,...,n
I v tomto modelu máme odezvu FEV, prediktory máme dva (k = 2): výška Height a pohlaví Sex. Stejně jako výše je pohlaví kódováno jako 0 a 1: 0 kóduje dívky a 1 kóduje chlapce. Maticový zápis je
/1.708\		(l	144.78	o\		
1.724		1	171.45	0	/A>\	
	=				x   Ä +	
2.795		1	160.02	1	väz	^653
\3.21j7			168.91	oy1		\£654/
Pro lineární prediktor platí tjí = e(FEVj) = (30 + fii x H eighty + /32 x Sexj. Střední hodnota FEV tedy podle modelu roste s výškou lineárně a přímky pro dívky a chlapce jsou rovnoběžné, obě se směrnicí fii. Koeficient fii vyjadřuje, o kolik litrů je vyšší střední hodnota FEV u dítěte, které je o 1 cm vyšší, porovnáváme-li dvě děti stejného pohlaví. Koeficient /32 vyjadřuje, o kolik litrů vyšší je střední hodnota FEV pro chlapce oproti stejně vysoké dívce.
# Model without interaction:
# Expected value of FEV increases linearly with Height,
# the lines for boys and girls are assumed to be parallel.
# Coefficient beta_l represents the increase in expected value
# of FEV associated with a 1 cm increase in Height when comparing
41
# two children of the same Sex.
# Coefficient beta_2 represents how much higher the expected value
# of FEV is for a boy compared to a girl of the same Height.
(5.) FEVj = (30 + /3i x Height^ + (32 x Sex^ + (33 x (Sex x Height). + eh   i = 1,..., n
Tento model se od předchozího liší přidáním nového prediktoru, který je funkcí prediktorů z předchozího modelu: výšky Height a pohlaví Sex. V kontextu regresních modelů mu říkáme interakce. Jde o funkci Sex x Height, pro všechny dívky má tedy nulovou hodnotu, zatímco pro chlapce obsahuje jejich výšky. V přítomnosti interakce mezi prediktory se na původní prediktory někdy odkazujeme jako na hlavní efekty. Máme tedy k = 3 a maticový zápis modelu je
/1.708 1.724
2.795 \3.211
Pro lineární prediktor platí rji = e(FEVj) = /?o + /?i x Height^ + /32 x S ex j + f33 x (Height x Sex)j. Střední hodnota FEV tedy podle modelu roste s výškou lineárně, tentokrát ale přímky pro dívky a pro chlapce nemusí být rovnoběžné. Směrnice pro dívky je fii,
( 1 144.78 1 171.45
1 160.02 1 168.91
0 0
1 0
0 0
160.02 0
/AA
w
£2
^653 \^654/
42
zatímco směrnice pro chlapce je {3i + (3%. Koeficient {3i vyjadřuje, o kolik litrů se u dívek zvýší střední hodnota FEV s 1 cm výšky navíc. U chlapců dostaneme obdobnou interpretaci pro koeficient fíi + /33. Rozdíl mezi střední hodnotou FEV pro chlapce oproti stejně vysoké dívce v tomto modelu závisí na výšce uvažovaných dětí.
# Model with interaction:
# Expected value of FEV increases linearly with Height,
# the lines for boys and girls are not necessarily parallel.
# Coefficient beta_l represents the increase in expected value
# of FEV associated with a 1 cm increase in Height for girls.
# beta_l  + beta_3 represents the increase in expected value
# of FEV associated with a 1 cm increase in Height for boys.
# The difference in expected values of FEV for boys and girls
# of the same Height here depends on the value of Height.
(Expected value of) FEV by Height and Sex
\       i       i       i       i       i       i r
120      130      140      150      160      170      180 190
Height [cm]
(6.) FEVj = f30 + f31xHe\ghtí + f32xI{9 < Age, < 12} + /33 xI{Age% > 12}+ehi = 1,    ... ,n
V tomto modelu je odezvou opět FEV. Prediktory jsou prakticky výška Height a věk Age, ale věk je kategorizovaný - rozdělený do tří kategorií - čímž se z něj stává kategoriální proměnná o třech úrovních. V lineárním modelu budeme jednu z těchto úrovní (v našem případě věk do 9 let včetně) brát jako referenční a zbylé úrovně (věk mezi 9 a 12 lety, věk nad 12 let) s ní budeme srovnávat. Za tímto účelem vytvoříme dva nové prediktory - jeden pro každou úroveň, kterou chceme srovnávat s referenční úrovní - a naplníme je jedničkami
43
pro pozorování, pro která věk dosahuje dané úrovně, a nulami pro zbylá pozorování. Dostaneme tak k = 3 a maticový zápis modelu bude
/1.708 1.724
2.795 \3.211
Pro lineární prediktor platí ^ = E(FEVj) = (30 + & x Height^ + (32 x I{9 < Age^ < 12} + /?3 x I{Agej > 12}. Střední hodnota FEV tedy podle modelu roste s výškou lineárně a přímky pro jednotlivé věkové kategorie jsou rovnoběžné, všechny se směrnicí fíi. Koeficient fii vyjadřuje, o kolik litrů je vyšší střední hodnota FEV u dítěte, které je o 1 cm vyšší, porovnáváme-li dvě děti ze stejné věkové kategorie. Koeficient /32 vyjadřuje, o kolik litrů vyšší je střední hodnota FEV pro dítě mezi 9 a 12 lety oproti stejně vysokému dítěti do 9 let. Koeficient /33 vyjadřuje, o kolik litrů vyšší je střední hodnota FEV pro dítě starší 12 let oproti stejně vysokému dítěti do 9 let.
/ 1 144.78 1 171.45
1 160.02 1 168.91
0 0
0 0
0\ 0
x
/AA
w
£2
^653 \^654/
# Expected value of FEV increases linearly with Height,
# the lines for the three age categories are parallel.
# Coefficient beta_l represents the increase in expected value
# of FEV associated with a 1 cm increase in Height when comparing
# two children from the same age category.
# Coefficient beta_2 represents how much higher the expected value
# of FEV is for a child aged between 9 and 12 compared to a child
# of the same Height aged 9 or less.
# Coefficient beta_3 represents how much higher the expected value
# of FEV is for a child aged 12 or more compared to a child
# of the same Height aged 9 or less.
44
(Expected value of) FEV by Height and Age
° o o
n-1-1-1-1-1-1-r
120      130      140      150      160      170      180 190
Height [cm]
(7.) FEV j = Po + ft x Height^ + /32 x Age^ + £i,i = l, ...,n
V tomto modelu je odezvou opět FEV a prediktory jsou výška Height a věk Age, tentokrát ale věk nekategorizujeme. Máme tedy k = 2 a maticový zápis modelu
/1.708\		(l	144.78   9 \		(e1\
1.724		1	171.45 8	/AA	
	=			x   Ä +	
2.795		1	160.02 16	w	^653
\3.2iy			168.91   15 J		\^654/
Pro lineárni prediktor platí tjí = e(FEVj) = (30 + Pi x Height^ + /32 x Ageí5 střední hodnota FEV podle modelu tedy lineárně roste s věkem i s výškou. Jde vlastně o rovinu, jak je vidět na obrázku niže.
45
(Expectation of) FEV by Age and Height
>
LU
Age
Pro každou pevnou hodnotu výšky roste střední hodnota s věkem lineárně a pro různé hodnoty výšky jsou tyto přímky rovnoběžné, jak je vidět na následujícím obrázku.
(Expected value of) FEV by Age and Height
>
LU co
	Height = 145 cm Height = 156 cm Height = 166 cm				0 0 0 8  8  1       ° o	o
					y    Q O	
	--~Q					
	I 5	I 10			i 15	
Age [years]
46
Analogická tvrzení i obrázek je možné vykreslit pro závislost střední hodnoty na výšce pro různé úrovně věku.
(Expected value of) FEV by Height and Age
Age = 8 years Age = 10 years - Age = 12 years		°   o o
		
o		o
I         I I 120      130 140	I 150	I        I        I I 160      170      180 190
m —
co —
CvJ —
Height [cm]
Při interpretaci koeficientu {3i uvažujeme dvě děti stejného věku. {3i vyjadřuje, o kolik litrů je střední hodnota FEV vyšší u dítěte, které je o 1 cm vyšší. Při interpretaci koeficientu f32 zase uvažujeme dvě děti stejné výšky. /32 vyjadřuje, o kolik litrů je střední hodnota FEV vyšší u dítěte, které je o 1 rok starší. Koeficient [30 v modelu reprezentuje střední hodnotu FEV pro dítě s nulovou výškou i věkem: v tomto modelu tedy nemá rozumnou interpretaci.
# Coefficient beta_l represents the increase in E(FEV) associated
# with a 1  cm increase in Height when comparing children of the sa-
# me Age.   Coefficient beta_2 represents the increase in E(FEV)
# associated with a 1 year increase in Age when comparing children
# of equal Height.  beta_0 represents E(FEV)   for a child aged 0 and
# with Height of 0 cm,   which makes little sense.
Kapitola 5 Lineární model v R
V této kapitole si ukážeme, jak daný lineárni model zadáme do a jak ® využijeme ke konstrukci odhadů parametrů modelu, posouzení kvality modelu i vyvození závěrů na základě modelu.
5.1   Zadání modelu do R
Připomeňme si první model z předchozí kapitoly.
(1.) FEVí = Ä) + Ä x Height^ + eh   i = l,...,n. Chceme-li tento model odhadnout pomocí ÍH. požádáme o to následujícím příkazem.
> lm(FEV ~  Height,   data = fev)
Příkazem lm () žádáme o odhad lineárního modelu. Prvním argumentem příkazu je takzvaná formula, ve které znakem ~ oddělujeme odezvu od prediktoru: na levou stranu píšeme odezvu (FEV) a na pravou stranu prediktor (Height). Argument data nám umožňuje odezvu a prediktory psát tak, aniž bychom u každého specifikovali název data. f rame, který je obsahuje. Bez využití tohoto argumentu bychom příkaz zadali následovně.
> lm(fev$FEV ~ fev$Height)
Na oba příkazy €S odpoví následujícím výstupem.
> lm(FEV ~  Height,   data = fev)
Call:
lm(formula = FEV ~  Height,   data = fev)
5.1. ZADÁNÍ MODELU DO R
48
Coefficients: (Intercept) Height -5.43268 0.05196
Koeficienty [30 a fii se pomocí ® odhadly na —5.43 a 0.05. Odhadujeme tedy, že střední hodnota FEV je u dítěte, které je o 1 cm vyšší, vyšší o 0.05 litrů. Pro snazší orientaci v odhadech koeficientů je @ vrací jako vektor, jehož složky mají jména. Odhad koeficientu odpovídajícího pre-diktoru pojmenuje názvem prediktoru, zatímco odhad koeficientu [30 pojmenuje Intercept, jelikož se v grafu závislosti střední hodnoty odezvy na prediktoru jedná o průsečík s osou y (y—intercept).
# Estimates of beta_0 and beta_l in the model are -5.43 and 0.05,
# respectively.   We estimate that the expected value of FEV
# increases by 0.05 1 with a 1  cm increase in a child's height.
Model s více prediktory zadáme tak, že je vyjmenujeme na pravé straně od znaku ~ argumentu formula a oddělíme od sebe znakem + (neuvažujeme-li mezi nimi interakci). Například model
(4.) FEVí = /30 + Ä x Height^ + /32 x Sex^ + eh   i = l,...,n zadáme následujícím příkazem.
> lm(FEV ~  Height + Sex,   data = fev) Call:
lm(formula = FEV ~  Height + Sex,   data = fev) Coefficients:
(Intercept) Height SexMale
-5.39026 0.05127 0.12512
Koeficienty [30, fii a /32 se pomocí C! odhadly na —5.39, 0.05 a 0.13. Odhadujeme tedy, že oproti dítěti stejného pohlaví je střední hodnota FEV u dítěte, které je o 1 cm vyšší, vyšší o 0.05 litrů. U chlapce odhadujeme střední hodnotu FEV o 0.13 litrů vyšší než u stejně vysoké dívky. Název SexMale u koeficientu /32 se skládá ze dvou částí (které ale nejsou odděleny): název prediktoru (Sex) a úroveň prediktoru, která je v matici plánu kódována jedničkou (Male). Díky tomuto značení víme, že rozdíl ve střední hodnotě FEV o 0.13 je ve prospěch chlapců, nikoliv dívek. Bez tohoto značení by bylo bezpečnější ověřit si kódování v matici plánu pomocí příkazu model. matrix () , který vrátí matici plánu modelu v argumentu příkazu.
5.1. ZADÁNÍ MODELU DO R
49
> model.matrix(lm(FEV ~ Height + Sex,  data = fev))[1:5, ] (Intercept)   Height SexMale
1			1 144.78	0
2			1 171.45	0
3			1 138.43	0
4			1 134.62	1
5			1 144.78	1
>	fev [1	:5, ;		
	Age	FEV	Height Sex	Smoker
1	9 1	.708	144.78 Female	Non
2	8 1	. 724	171.45 Female	Non
3	7 1	. 720	138.43 Female	Non
4	9 1	. 558	134.62 Male	Non
5	9 1	.895	144.78 Male	Non
# Estimates of beta_0,  beta_l and beta_2 in the model are -5.40,
# 0.05,   and 0.13,  respectively.   We estimate that the expected
# value of FEV increases by 0.05 1 with a 1  cm increase
# in height for children of the same sex.   The difference
# between the expected values of FEV for children of the same
# height is 0.13 in favour of a boy  (compared to a girl).
Chceme-li do modelu zahrnout interakci, oddělíme od sebe příslušné prediktory znakem * místo +. Pro model
(5.) FEVí = Po + ft x Height^ + (32 x Sex; + (33 x (Height x Sex); +
bychom použili následující příkaz.
> lm(FEV ~  Height  *  Sex,   data = fev) Call:
lm(formula = FEV " Height  * Sex,   data = fev) Coefficients:
(Intercept) Height SexMale    Height:SexMale
-4.31822 0.04426 -1.54563 0.01081
5.1. ZADÁNÍ MODELU DO R
50
Koeficienty (30, (3ľ, (32 a (33 se pomocí Q odhadly na -4.32, 0.04, -1.55 a 0.01. Odhadujeme tedy, že s každým centimetrem výšky je spojený nárůst střední hodnoty FEV u dívek o 0.04 litrů a u chlapců o 0.06 litrů. Název interakce Height: SexMale obsahuje názvy prediktorů, mezi kterými interakci uvažujeme, a také úroveň kategoriálního prediktorů, která je v matici plánu kódována jedničkou. Oddělovací znak : se v Cl používá pro interakci.
> model.matrix(lm(FEV ~ Height * Sex,  data = fev))[1:5, ]
(Intercept) Height
1 1 144.78
2 1 171.45
3 1 138.43
4 1 134.62
5 1 144.78
exMale Height:SexMale
0 0.00
0 0.00
0 0.00
1 134.62 1 144.78
Stejný model můžeme zadat také explicitním vyjmenováním všech prediktorů včetně interakce oddělených znakem +.
> lm(FEV ~  Height + Sex + Height:Sex,   data = fev)
Call:
lm (formula = FEV ~  Height + Sex + Height:Sex,   data = fev) Coefficients:
(Intercept) Height SexMale    Height:SexMale
-4.31822 0.04426 -1.54563 0.01081
Pro usnadnění zadávání modelů s interakcemi se v argumentu formula používá speciální syntaxe, ve které n—tá mocnina značí zahrnutí hlavních efektů a všech interakcí až do n—tého řádu. Pomocí této syntaxe můžeme výše uvedený model zadat následovně.
> lm(FEV ~   (Height + Sex)"2,   data = fev)
Call:
lm(formula = FEV ~   (Height + Sex)"2,   data = fev) Coefficients:
(Intercept) Height SexMale    Height:SexMale
-4.31822 0.04426 -1.54563 0.01081
5.1. ZADÁNÍ MODELU DO R
51
# Estimates of beta_0,  beta_l,  beta_2,   and beta_3 in the model
# are -4.30,   0.04,   -1.55,   and 0.01,  respectively.  We estimate
# that a 1 cm increase in height is associated with an increase
# in the expected value of FEV by 0.04 1 for girls and by 0.06 1
# for boys.
Pro interakce vyššího než druhého řádu bychom v modelu museli mít více proměnných. Pro tři proměnné můžeme uvažovat například následující modely s interakcemi.
> lm(FEV ~   (Height + Sex)"2 + Age,   data = fev)
Call:
lm(formula = FEV ~   (Height + Sex)"2  + Age,   data = fev) Coefficients:
(Intercept) Height SexMale Age
-3.09862 0.03199 -1.80829 0.06685
Height:SexMale 0.01276
> lm(FEV ~   (Height + Sex + Age)"2,   data = fev)
Call:
lm(formula
FEV
(Height + Sex + Age)"2,  data = fev)
Coefficients: (Intercept) -0.763023 Height:SexMale 0.005764
Height 0 .017919 Height:Age 0.002256
SexMale -0.876319 SexMale:Age 0.010684
Age ■0.303109
> lm(FEV ~   (Height + Sex + Age)"2  - Height:Age,   data = fev)
Call:
lm(formula = FEV ~   (Height + Sex + Age)"2 - Height:Age,   data = fev)
Coefficients: (Intercept) -3.329648 Height:SexMale 0.007958
Height 0 . 034311 SexMale:Age 0.028716
SexMale ■1 . 347481
Age 0 . 054187
5.1. ZADÁNÍ MODELU DO R
52
> lm(FEV ~   (Height + Age + Sex)"3,   data = fev)
Call:
lm(formula
FEV
(Height + Age + Sex)"3,  data = fev)
Coefficients:
(Intercept) -3.380e+00 SexMale 3.460e+00 Age:SexMale -5.661e-01
Height 3.464e-02 Height:Age -4.464e-05 Height:Age:SexMale 3.577e-03
Age 6 .126e-02 Height:SexMale -2.147e-02
Chceme-li jako prediktory využít funkce stávajících prediktorů, můžeme hodnoty těchto nových prediktorů nejprve spočítat a přidat jako další sloupec do příslušného data. f rame. Ten pak zadáme do argumentu formula příkazu lm () . Tak bychom postupovali například u modelu
(6.) FEV, = /30+/3iXHeight^xI{9 < Age, < 12}+/33xI{Aget > 12}+eť
,n.
> quantile(fev$Age,  probs = c(0:3)/3)
0% 33.33333% 66.66667% 100% 3 9 11 19
> fev$Age.cat <- cut(fev$Age, breaks = c(3, 9, 12, 19), + include.lowest = TRUE)
> summary(fev$Age.cat)
[3,9]     (9,12] (12,19]
309 228 117
> lm(FEV ~  Height + Age.cat,   data = fev)
Call:
lm(formula = FEV ~  Height + Age.cat,   data = fev) Coefficients:
(Intercept) Height      Age.cat(9,12]     Age.cat(12,19]
-4.50925 0.04525 0.12235 0.42566
5.1. ZADÁNÍ MODELU DO R
53
U jednodušších funkcí prediktoru je výhodnější zadat požadovanou funkci přímo uvnitř agru-mentu formula. V takovém případě je potřeba nový prediktor zadat uvnitř funkce I () . Například model
(2.) FEVí = (30 + ft x Height^ + (32 x Height? + eh   i = l,...,n zadáme následujícím příkazem.
> lm(FEV ~  Height + I(Height"2),   data = fev) Call:
lm(formula = FEV ~  Height + I(Height"2),   data = fev) Coefficients:
(Intercept) Height I(Height~2)
6.0268779      -0.0984444 0.0004891
Nakonec si ještě ukážeme zjednodušenou syntaxi pro dva speciální modely. Model bez prediktoru
FEVj = (30 + Si,   i = l,...,n
zadáme následujícím příkazem.
> lm(FEV ~  1,   data = fev)
Call:
lm (formula = FEV ~  1,   data = fev)
Coefficients: (Intercept) 2 . 637
Naopak model, který jako prediktory obsahuje všechny sloupce data . frame, ze kterého pochází odezva, zadáme následovně.
> lm(FEV ~   .,   data = fev) Call:
lm (formula = FEV ~   .,   data = fev) Coefficients:
5.2. CELKOVÉ CHARAKTERISTIKY MODELU
54
Intercept) -4 .46536 SmokerNon 0.13668
Age 0.01858 Age.cat (9,12] 0 .11516
Height 0.04253 Age.cat (12,19] 0.40769
SexMale 0 .14928
5.2   Celkové charakteristiky modelu
Vraťme se nyní k modelu
(2.) FEVí = (30 + 0! x Height^ + (32 x Sex^ + eh   i = l,...,n.
Již jsme viděli, že příkaz lm () vrátí odhady koeficientů modelu. Zároveň s odhady spočítá <Ř mnoho charakteristik užitečných pro inferenci. Abychom k nim mohli lépe přistupovat, uložíme si vše, co se po zadání příkazu lm () napočítalo, například do proměnné mod.
> mod <- lm(FEV ~  Height + Sex,   data = fev)
Funkci model. matrix () , popsanou výše, nyní můžeme uplatnit přímo na proměnnou mod. Funkce coefficients () uplatněna na proměnnou mod vrátí odhady koeficientů modelu.
> model.matrix(mod) [1:5, ]
(Intercept)   Height SexMale
1 1  144.78 0
2 1  171.45 0
3 1  138.43 0
4 1  134.62 1
5 1  144.78 1
> coefficients(mod)
(Intercept) Height SexMale
-5.39026319    0.05127185 0.12512339
Přehled nejdůležitějších charakteristik, které lze získat z proměnné mod, poskytuje příkaz summary () .
> summary(mod)
Call:
lm(formula = FEV ~  Height + Sex,   data = fev)
5.2. CELKOVÉ CHARAKTERISTIKY MODELU
55
Residuals:
Min 1Q    Median 3Q Max
-1.6763 -0.2505    0.0001    0.2347 2.0722
Coefficients
(Intercept) -Height SexMale
Estimate -5.390263 0.051272 0.125123
Std, 0 . 0 .
0 ,
Error 180082 001167 033801
t value -29.932 43. 933 3 .702
Pr
<
<
:>it i)
2e-16 2e-16
0.000232
Signif.  codes:     0   '***'   0.001   '**' 0.01
0 . 05
0 .1
Residual standard error:   0.4265 on 651 degrees of freedom Multiple R-squared:     0.7587,Adjusted R-squared: 0.758 F-statistic:     1024 on 2  and 651 DF,    p-value:   < 2.2e-16
První řádek výstupu obsahuje definici modelu. Následují empirické kvantily reziduí e = Y — X/3. Také poslední tři řádky se týkají celého modelu. Obsahují odhad a odmocniny z rozptylu a2 náhodné chyby e, koeficienty determinace a test hypotézy o nulovosti všech koeficientů modelu současně.
Odhad a je odmocninou z nevychýleného odhadu a2 = eTe/(n — p). V normálním lineárním modelu má náhodná veličina [n — p)a2/a2 rozdělení \2 s n — P stupni volnosti. Ve výstupu z fuknce summary () jsou tyto stupně volnosti uvedeny za odhadem a. Jedná se o počet pozorování snížený o rozměr vektoru /3, jak snadno ověříme.
> nrow(fev)   - length(coefficients(mod))     # n-p
[1] 651
Koeficient determinace
F?2   -   1   -   ž^ = l(^ Y) -    _ y
reprezentuje proporci variability dat vysvětlenou modelem. Náš model vysvětluje 76 % variability proměnné FE V. Koeficient determinace vzroste pokaždé, když do modelu přidáme prediktor, proto je pro srovnání vnořených modelů vhodnější adjustovaný (upravený) koeficient determinace
,2 _n Eti(Y*-Y*)2/(ri-p) Ľľ=i(*-ň7(n-l) který bere do úvahy i počet prediktorů v modelu. V našem případě je téměř stejný jako koeficient determinace, i?^dj = 0.76. Při velkém počtu dat a malém počtu koeficientů v modelu jsou si totiž hodnoty dělitelů n — p a n — 1 velmi blízké.
-^adj — 1 -r>\o // i\> (5.1)
5.2. CELKOVE CHARAKTERISTIKY MODELU 56
Na posledním řádku máme pozorovanou hodnotu testové statistiky F—testu hypotézy H0 : (Pi, faV = 0 proti alternativě Hi : ((3i, (32)T Ý 0. Stupně volnosti příslušného F—rozdělení jsou rozměr testovaného vektoru a stupně volnosti %2—rozdělení odhadu rozptylu a2 náhodné chyby e, tedy p — 1 a n — p. Nezamítnutí nulové hypotézy sice neznamená její potvrzení, ale vyvolávalo by pochybnosti o vhodnosti prediktorů pro popis střední hodnoty odezvy. V našem případě však nulovou hypotézu jednoznačně zamítáme (p—hodnota < 1CT15), tudíž tento výsledek pochybnosti o vhodnosti prediktorů nevyvolává.
Pro přístup k jednotlivým výsledkům zobrazeným ve výstupu ze summary (mod) potřebujeme znát jejich interní názvy. Ty můžeme najít ve výstupu z funkce s t r () uplatněné na objekt summary(mod) .
> str(summary(mod))
List  Of 11
$ call :   language Im(formula = FEV ~ Height + Sex,   data = fev)
$ terms :Classes   'terms',   'formula'     language FEV ~ Height + Sex
.- attr(*,   "variables")= language list(FEV,  Height, Sex) - attr(*,   "factors")= int   [1:3,   1:2]   0  10  0 0 1 - attr(*,   "dimnames")=List of 2 ..$   :   chr   [1:3]   "FEV"  "Height" "Sex" chr   [1:2]   "Height" "Sex"
"term.labels")= chr   [1:2]   "Height" "Sex" "order")= int   [1:2]   1 1 "intercept")= int 1 "response")= int 1
".Environment")=<environment: R_GlobalEnv> "predvars")= language list(FEV,   Height, Sex)
"dataClasses")= Named chr  [1:3]   "numeric"  "numeric" "factor" attr(*,   "names")= chr   [1:3]   "FEV"  "Height" "Sex" residuals        :  Named num   [1:654]   -0.3249 -1.6763 0.0127 -0.0791 -0.263 ... .- attr(*,   "names")= chr   [1:654]   "1"  "2"  "3"  "4" ...
coefficients   :  num   [1:3,   1:4]   -5.39026 0.05127  0.12512  0.18008  0.00117 ... attr(*,   "dimnames")=List of 2
..$  :  chr   [1:3]   "(Intercept)"  "Height" "SexMale" ..$   :   chr   [1:4]   "Estimate"  "Std.  Error"  "t value" "Pr(>|t|)" $ aliased :  Named logi   [1:3]  FALSE FALSE FALSE
..- attr(*,   "names")= chr   [1:3]   "(Intercept)"  "Height" "SexMale"
. .   . .$
- attr(*
- attr(*
- attr(*
- attr(*
- attr(*
- attr(*
- attr(*
$ sigma
$ df
$ r.squared
$ adj.r.squared
$ fstatistic
num 0.427 int [1:3] num 0.75 9 num 0.758 Named num ..- attr(*,   "names")= chr $ cov.unsealed  :  num   [1:3,   1:3]   1.78e-01 -1.14e-03 2 ..- attr(*,   "dimnames")=List of 2
..   ..$   :   chr   [1:3]   "(Intercept)"  "Height" "SexMale" ..   ..$   :   chr   [1:3]   "(Intercept)"  "Height" "SexMale" - attr(*,   "class")= chr "summary.lm"
3  651 3
[1:3]   1024 2 651 [1:3]   "value"  "numdf" "dendf"
13e-03 -1.14e-03 7.49e-06
Z výstupu výše vidíme, že například a, R2 a R2^ získáme pomocí následujících příkazů
5.3. INFERENCE PRO JEDNOTLIVÉ KOEFICIENTY
57
> summary(mod)$sigma    # estimate of sigma [1] 0.4265405
> summary(mod)$r.squared    # coefficient  of determination (R"2) [1] 0.7587368
> summary(mod)$adj.r.squared    # adjusted R~2 [1] 0.7579956
5.3   Inference pro jednotlivé koeficienty
Tabulka uprostřed výstupu ze summary (mod) obsahuje informace užitečné pro posouzení statistické významnosti a výpočet konfidenčních intervalů pro jednotlivé koeficienty modelu. Uložíme si ji do samostatné proměnné.
> tab <- summary(mod)$coefficients
> tab
Estimate    Std.  Error t value Pr(>|t|)
(Intercept)   -5.39026319 0.180082479 -29.932191 1.821608e-124
Height              0.05127185 0.001167051 43.932823 6.861183e-197
SexMale            0.12512339 0.033800905 3.701776 2.322264e-04
Jedná se o matici s pojmenovanými řádky i sloupci.
První sloupec (Estimate) obsahuje odhady koeficientů.
> tab[, 1]
(Intercept) Height SexMale
-5.39026319    0.05127185 0.12512339
> coefficients(mod)
(Intercept) Height SexMale
-5.39026319    0.05127185 0.12512339
Druhý sloupec (Std. Error) obsahuje odmocniny z odhadnutých rozptylů odhadů jednotlivých koeficientů. Z teorie víme, že
Var/3 = (72 (XTX
5.3. INFERENCE PRO JEDNOTLIVÉ KOEFICIENTY
58
> tab[, 2]
(Intercept) Height SexMale
0.180082479 0.001167051 0.033800905
> X <- model.matrix(mod)
> summary(mod)$sigma * sqrt(diag(solve(t(X)   %*% X)))
(Intercept) Height SexMale
0.180082479 0.001167051 0.033800905
Třetí sloupec matice (t value) obsahuje pozorované hodnoty testových štatistik pro testy hypotéz H0 : fy = 0 pro jednotlivé složky vektoru koeficientů proti alternativám Hi : (3i ^ 0. Testové statistiky jsou
71= Ä
<P (xTx),-f
jedná se tedy o podíl prvního a druhého sloupce matice tab.
> tab[, 3]
(Intercept) Height SexMale
-29.932191      43.932823 3.701776
> tab [,   1] /tab [, 2]
(Intercept) Height SexMale
-29.932191      43.932823 3.701776
V normálním lineárním modelu má každá z těchto testových statistik za platnosti příslušné nulové hypotézy ŕ—rozdělení s n — p stupni volnosti, Ti ~ tn-p. Příslušné p—hodnoty najdeme v posledním sloupci matice (P r (> 11 | )). Pro jejich výpočet můžeme využít toho, že p—hodnota představuje pravděpodobnost, že za platnosti nulové hypotézy má testová statistika hodnotu, kterou jsme pozorovali, nebo hodnoty, které by ještě více svědčily ve prospěch alternativy. V našem případě ve prospěch alternativy svědčí hodnoty testové statistiky vzdálenější od nuly, a to na kladnou i zápornou stranu.
> tab[, 4]
(Intercept) Height SexMale
1.821608e-124  6.861183e-197 2.322264e-04
5.3. INFERENCE PRO JEDNOTLIVÉ KOEFICIENTY
59
> 2 * pt(-abs(tab[,   3]),  df = summary(mod)$df[2])
(Intercept) Height SexMale
1.821608e-124  6.861183e-197 2.322264e-04
Invertováním těchto ŕ—testů získáme konfidenční intervaly pro jednotlivé koeficienty
(fit - ŕ„-P(l - a/2) y/v*(XrX)-i ; p. + tn_p{l _ a/2) ^2(XTX)r^ .
V © je vrací funkce confint() .
> confint(mod)
2.5 % 97.5 %
(Intercept) -5.74387579 -5.03665059 Height 0.04898022 0.05356349
SexMale 0.05875144 0.19149534
> alpha <- 0.05
> df <- summary(mod)$df[2]
> lwr <- tab[,   1]   - qt(l - alpha/2,  df = df)   * tab[, 2]
> upr <- tab[,   1]   + qt(l - alpha/2,  df = df)   * tab[, 2]
> data.frame(lwr, upr)
lwr upr (Intercept) -5.74387579 -5.03665059 Height 0.04898022 0.05356349
SexMale 0.05875144 0.19149534
Připomeňme si, že hustota ŕ—rozdělení je podobná hustotě standardního normálního rozdělení a s rostoucími stupni volnosti se jí blíží, avšak ŕ—rozdělení má těžší chvosty. Kvantily ŕ—rozdělení se proto s rostoucími stupni volnosti blíží kvantilům standardního normálního rozdělení, v absolutní hodnotě jsou ale o něco větší. Nemáme-li příliš malý rozsah dat, pro rychlou představu o konfidenčních intervalech díky tomu stačí vzít první sloupec matice a odečíst (resp. přičíst) k němu dvojnásobek («0.975 ~ 1.96) druhého sloupce.
Každý z výše uvedených ŕ—testů a konfidenčních intervalů se týká jediného koeficientu. Jsou-li splněny předpoklady normálního lineárního modelu, každý z těchto testů má hladinu významnosti a a každý z těchto konfidenčních intervalů má pravděpodobnost pokrytí 1 — a. Pravděpodobnost, že všechny konfidenční intervaly pokryjí odhadované parametry, může být nižší a společná hladina všech testů může být vyšší. Chceme-li kontrolovat hladinu významnosti testu o více koeficientech současně, použijeme místo jednotlivých ŕ—testů společný F—test, kterému se budeme podrobněji věnovat v části 5.5. Zajímá-li nás společné pokrytí konfidenčních intervalů, nahradíme je konfidenčními regiony nebo konfidenčními pásy, kterým se budeme věnovat v části 5.6.
5.4. INFERENCE PRO LINEÁRNÍ KOMBINACE KOEFICIENTŮ
60
5.4   Inference pro lineární kombinace koeficientů
Někdy se můžeme vedle jednotlivých koeficientů zajímat také o jejich lineární kombinace. Například v modelu
(5.) FEVj = (30 + ft x Height^ + (32 x Sex^ + (33 x (Sex x Height). + e,h   i = 1,..., n
vyjadřuje koeficient fii nárůst střední hodnoty FEV spojený s nárůstem výšky u dívek o 1 cm. Pro chlapce má stejnou interpretaci součet fii + [33.
> mod <- lm(FEV ~  Height  *  Sex,   data = fev)
> summary (mod)
Call:
lm(formula = FEV ~  Height  *  Sex,   data = fev) Residuals:
Min IQ      Median 3Q Max
-1.54654  -0.25282    0.00649    0.25666 2.00491
Coefficients:
Height:SexMale    0.010810      0.002409      4.487  8.54e-06 ***
Signif.  codes:     0   •***•   0.001   '**'   0.01   '*'   0.05   '.'   0.1   '   ' 1
Residual standard error:   0.4204 on 650 degrees of freedom Multiple R-squared:     0.766,Adjusted R-squared: 0.7649 F-statistic:  709.2 on 3 and 650 DF,    p-value:  < 2.2e-16
Tento efekt odhadneme jako J3i + J33 = 0.06. Odhad rozptylu Var(f31 + (33) = Var (5i + Var (33 + 2 x Cov(/3i, (33), který bychom potřebovali ke konstrukci testu nebo konfidenčního intervalu pro fii + [33, ale ve výstupu z funkce summary ()  nenajdeme. Za tímto účelem použijeme
funkci vcov () , která vrací odhad a2 (XTX)_1 kovarianční matice odhadů koeficientů v modelu, který je jejím argumentem.
Estimate Std.  Error t value Pr(>|t|)
(Intercept) Height SexMale
-4.318219      0.297637 -14.508    < 2e-16 *** 0.044262      0.001940    22.815    < 2e-16 *** -1.545629      0.373843    -4.134 4.02e-05 ***
> est.var <- vcov(mod)
> est.var
5.4. INFERENCE PRO LINEÁRNÍ KOMBINACE KOEFICIENTŮ
61
(Intercept) Height SexMale
Height:SexMale
(Intercept) Height SexMale
Height:SexMale
(Intercept) 0.088588026 --0.000575606 -0.088588026
0.000575606 -Height:SexMale 5.756060e-04 -3.763647e-06 -8.970667e-04 5.804094e-06
Height 5.756060e-04 3.763647e-06 5.756060e-04 3.763647e-06
SexMale -0 .0885880258 0 .0005756060 0 .1397583375 -0 .0008970667
Nyní již můžeme zkonstruovat testovou statistiku pro test hypotézy H0 : {3i + (53 = 0 proti alternativě Hi : f3i + [33 ý 0 i konfidenční interval pro fii + [33.
> est <- s vim (coefficients (mod) [c(2,4)])
> est.var.sqrt <- sqrt(sum(diag(est.var)[c(2, 4)]) + + 2  * est.var [2,   4])
> c(est,   est.var.sqrt)   #Estimate and Std.Error for beta_l+beta_3
[1]   0.055072110 0.001428442
> test.stat <- est/est.var.sqrt
> test.p <- 2  * pt(-abs(test.stat),  df = summary(mod)$df[2])
> c (test, stat,  test.p)   # t value and P(>/tl)   for beta_l  + beta_3
[1]     3.855396e+01 4.200169e-170
> ci <- est + c(-l,   1)   * qnorm(l-alpha/2)   * est.var.sqrt
> ci  # predidence interval for beta_l  + beta_3
[1]   0.05227242 0.05787181
Obecněji, pro test o lineární kombinaci aT/3 = a0[30 + a\fi\ + ... + s H0 : aT/3 = 0 a Hi : aT/3 ^ 0, použijeme testovou statistiku
T=^     ^-, (5.2)
a2 aT(XTX)-!a
která má za platnosti normálního lineárního modelu a nulové hypotézy ŕ—rozdělení sn—p stupni volnosti. Invertováním testu dostaneme konfidenční interval pro aT/3. V kódu výše bychom změnili jenom výpočet odhadu a odmocniny z jeho odhadnutého rozptylu.
5.5. INFERENCE PRO VÍCE (LINEÁRNÍCH KOMBINACÍ) KOEFICIENTŮ
62
> est <- t(a)   %*% coefficients(mod)
> est.var.sqrt <- sqrt(t(a)   %*% est.var %*% a)
5.5   Inference pro více (lineárních kombinací) koeficientů
Souvislost mezi prediktorem výška Height a střední hodnotou odezvy FEV v modelu (5.) prakticky popisují směrnice [3i a [3i + [33. Kdyby výsledky jednorozměrných testů o jejich nulovosti popsané v částech 5.3 a 5.4 nebyly tak jednoznačné, bylo by pro vyčerpávající odpověd na otázku o souvislosti výšky Height se střední hodnotou FEV potřeba testovat nulovost obou směrnic zároveň. Za tímto účelem můžeme využít F—test pro testování hypotéz o několika nezávislých lineárních kombinacích koeficientů modelu zároveň. Tvoří-li koeficienty m nezávislých lineárních kombinací řádky matice a, testová statistika pro test hypotézy H0 : A/3 = 0 proti alternativě Hi : A/3 ý 0 má tvar
F = ^^(a3)t(A(XtX)-1At)-1a3 (5.3) m a2
a za platnosti normálního lineárního modelu a nulové hypotézy má F—rozdělení s m a n — p stupni volnosti. Následujícím způsobem pak můžeme najednou provést test o nulovosti obou směrnic odpovídajících výšce Height.
> A <- rbind(c(0,   1,   0,   0),  c(0,   1,   0, 1))
> beta.hat <- coefficients(mod)
> var.beta.hat <- vcov(mod)
> est <- A %*% beta.hat
> est  # estimate of the vector  (beta_l,  beta_l  + beta_3)
[,1]
[1,] 0.04426220 [2,] 0.05507211
> F <- t(est)   %*% solve(A %*% var.beta.hat %*% t(A))   %*% est/2
> p <- l-pf(F,  dfl = 2,  df2 = summary(mod)$df[2])
> # test statistic and p-value for testing the hypothesis
> # that both beta_l = 0 and beta_l  + beta_3 = 0
> c(F, p)
[1]   1003.476 0.000
Vzhledem k jednoznačným výsledkům jednotlivých ŕ—testů není překvapivé, že také výsledek F—testu je jednoznačný. Obecně ale nemusí F—test vést ke stejnému výsledku jako jednotlivé ŕ—testy. Připomeňme ale, že použijeme-li F—test (5.3) k testování hypotézy o jediné lineární
5.6. ODHADY STŘEDNÍ HODNOTY A PREDIKCE NA ZÁKLADĚ MODELU
63
kombinaci koeficientů modelu, jeho testová statistika je kvadrátem testové statistiky ŕ—testu z části 5.4. Jelikož kvadrát náhodné veličiny s ŕ—rozdělením s n stupni volnosti má F—rozdělení s 1 a n stupni volnosti, výsledky testů budou v takovém případě stejné.
Nulovou hypotézu H0 : fii = 0 a fii + [33 = 0 můžeme ekvivalentně zapsat jako H0 : fii = 0 a /33 = 0. Za její platnosti se tedy model (5.) redukuje na model
(3.) FEVí = (30 + Ä x Sexí + eh   i = l,...,n.
Test hypotézy H0 : fíi = 0 a fíi + (53 = 0 proti alternativě Hi : fíi ^ 0 nebo fii + f33 ^ 0 proto můžeme vnímat i jako test o možnosti přechodu od modelu (5.) k modelu (3.). Takovým testům se budeme věnovat v části 6.1.
5.6   Odhady střední hodnoty a predikce na základě modelu
Důležitým speciálním případem lineárních kombinací koeficientů diskutovaných v části 5.4 jsou kombinace odpovídající konkrétním hodnotám prediktorů. Uvažujme opět model (5.) a počítejme odhad střední hodnoty FEV pro dívku vysokou 110 cm. Podle modelu (5.) platí, že
FEV = (30 + /3i x Height + (32 x Sex + (33 x (Sex x Height) + e
a E e = 0. Pro střední hodnotu FEV pak platí vztah
E(FEV) = f30 + Ä x Height + (52 x Sex + f33 x (Sex x Height).
Připomeňme, že dívky jsou v proměnné Sex kódovány nulou. Pro střední hodnotu FEV dívky vysoké 110 cm tak dostáváme vztah
E(FEV I Sex = 0, Height = 110) = #, + Ä x 110 + (52 x 0 + /33 x 0 = #, + HO Ä.
Pro lepší viditelnost hodnot prediktorů, pro které počítáme střední hodnotu odezvy, jsme v rovnici výše použili značení, které je obvyklejší v situacích, kdy jsou prediktory považovány za náhodné veličiny, na jejíchž konkrétních hodnotách podmiňujeme.
Střední hodnota odezvy za konkrétních hodnot prediktorů, v našem případě E(FEV | Sex = 0, Height = 110), je tedy lineární kombinací koeficientů modelu. V našem případě tvoří její koeficienty vektor a= (1,110,0,0)T.V tomto kontextu se pro koeficienty lineární kombinace používá spíše značení x = (1,110, 0, 0)T připomínající (generický) řádek regresní matice.
Jak jsme si již rozmysleli v části 5.4, lineární kombinaci xT/3 můžeme odhadnout pomocí xT/3. V našem případě tedy
É(FEV I Sex = 0, Height = 110) = /30 + 110 x & = -4.32 + 110 x 0.04 « 0.55.
Pro testování hypotézy o konkrétní hodnotě E(FEV | Sex = 0, Height = 110) můžeme použít testovou statistiku (5.2). Jejím invertováním pak dostaneme konfidenční interval pro E (FEV | Sex = 0, Height = 110), jak jsme již také viděli v části 5.4.
5.6. ODHADY STŘEDNÍ HODNOTY A PREDIKCE NA ZÁKLADĚ MODELU
64
> # Expected value of FEV for a 110 cm tall girl: beta_0 + 110 *
> # beta_2
> x <- c(l,   110,   0, 0)
> est <- t(x)   %*% beta.hat
> est.var.sqrt <- sqrt(t(x)   %*% var.beta.hat %*% x)
> # Estimate and Std.Error for beta_0 + 110  * beta_2
> c(est, est.var.sqrt)
[1]   0.55062368 0.08657273
> test.stat <- est/est.var.sqrt
> test.p <- 2  * pt(-abs(test.stat),  df = summary(mod)$df[2])
> # t value and P(>\t\)   for beta_0 + 110  * beta_2
> c(test.stat, test.p)
[1]   6.360244e+00 3.798466e-10
> ci <- est + c(-l,   1)   * qnorm(l  - alpha/2)   * est.var.sqrt
> # confidence interval for beta_0 +110  * beta_2
> ci
[1]   0.3809442 0.7203031
Související úlohou je predikce hodnoty nového pozorování při konkrétních hodnotách pre-diktorů. Nové pozorování je náhodná veličina, proto jej nemůžeme odhadovat, ale můžeme jej predikovat. Pro dívku vysokou 110 cm je podle modelu (5.)
Odhadneme-li ve vztahu výše koeficienty /30 a {3i a predikujeme-li e jeho střední hodnotou, dostaneme predikci FEV pro dívku vysokou 110 cm ve tvaru
predikce tedy bude stejná jako odhad příslušné střední hodnoty. Rozdíl ale nastane při tvorbě predikčního intervalu. Tam musíme vedle variability plynoucí z nahrazení koeficientů v lineární kombinaci jejími odhady vzít do úvahy také variabilitu plynoucí z nahrazení náhodné chyby její střední hodnotou. Celkový rozptyl tedy nebude o"2xT(XTX) x, nýbrž o"2xT(XTX) x + a2. Predikční interval pak bude
FEV = & + Ä x 110 + (52 x 0 + f33 x 0 + e = f30 + 110 Ä + e.
(FEV I Height = 110, Sex = 0) = $0 + 110 (3-_
5.6. ODHADY STŘEDNÍ HODNOTY A PREDIKCE NA ZÁKLADĚ MODELU
65
> # Prediction interval  for FEV for a 110 cm tall girl
> est.var.sqrt <- sqrt(t(x)   %*% var.beta.hat %*% x + + (summary(mod)$sigma)"2)
> pi <- est + c(-l,   1)   * qnorm(l  - alpha/2)   * est.var.sqrt
> pi
[1]   -0.2906489 1.3918962
Jelikož záporné hodnoty FEV nedávají smysl, můžeme predikční interval v této sitaci omezit na kladná čísla.
K výpočtu odhadu střední hodnoty i predikce odezvy pro dané hodnoty prediktorů a odpovídajících intervalů můžeme v <H využít funkce predict () . Jejími argumenty jsou model, na základě kterého chceme odhadovat nebo predikovat, a hodnoty prediktorů, pro které chceme odhady nebo predikce odezvy počítat (ty vkládáme jako data . frame, jehož sloupce mají stejné názvy jako prediktory v modelu). Volbou interval = "confidence" dále požádáme o konfidenční intervaly, zatímco volbou interval = "prediction" požádáme
0 predikční intervaly. Nepožádáme-li o interval, funkce vrátí jenom odhad střední hodnoty/predikci.
> # estimate and confidence interval for the expected value
> # of FEV for a 110 cm tall girl
> predict(mod,
+ newdata = data.frame(Height = 110,   Sex = "Female"),
+ interval = "confidence")
fit lwr upr
1 0.5506237 0.3806277 0.7206196
> # prediction and prediction interval  for FEV for a 110 cm
> # tall girl
> predict(mod,
+ newdata = data.frame(Height = 110,   Sex = "Female"),
+ interval = "prediction")
fit lwr upr
1  0.5506237 -0.2922183 1.393466
> # estimated expected value/prediction of FEV for a 110 cm
> # tall girl
> predict(mod,
+ newdata = data.frame(Height = 110,   Sex = "Female"))
1
0.5506237
5.6. ODHADY STŘEDNÍ HODNOTY A PREDIKCE NA ZÁKLADĚ MODELU
66
V argumentu newdata funkce predict () můžeme zadat data. frame s více řádky, což se hodí u vykreslování odhadů nebo predikcí na základě modelu jako funkcí prediktorů.
> xx <- seq(from = min(fev$Height),  to = max(fev$Height), + length = 501)
> # a grid o heights for which to compute estimates/predictions
>
> yyl  <- predict(mod,
+ newdata = data.frame(Height = xx,
+ Sex = "Female"))
> # estimates/predictions for girls on the grid of heights
> yy2  <- predict(mod,
+ newdata = data.frame(Height = xx,   Sex = "Male"))
> # estimates/predictions for boys on the grid of heights
>
> plot(fev$FEV ~ fev$Height,
+ main =  "(Expected value of)   FEV by Height and Sex",
+ xlab =  "Height   [cm]",   ylab = "FEV   [1]",   type = "n")
> #   (empty)  plot  of FEV against Height
> with(fev,
+ points(FEV[Sex ==  "Female"]   ~  Height[Sex == "Female"],
+ col = "lightpink"))
> # pink points for girls
> with(fev,
+ points(FEV[Sex ==  "Male"]   ~  Height[Sex == "Male"],
+ col = "cornflowerblue"))
> # blue points for boys
> lines(xx,  yyl,   col = "lightpink2",   lwd = 2)
> # estimated expected/predicted FEV by Height for girls
> lines(xx,  yy2,   col = "dodgerblue2",   lwd = 2)
5.6. ODHADY STŘEDNÍ HODNOTY A PREDIKCE NA ZÁKLADĚ MODELU
67
> # estimated expected/predicted FEV by Height for boys Přidání konfidenčních intervalů je s funkcí predict () přímočaré.
> # confidence intervals for the fitted lines
> yyl.conf <- predict(mod,
+ newdata = data.frame(Height = xx,
+ Sex = "Female"),
+ interval = "confidence")
> yy2.conf <- predict(mod,
+ newdata = data.frame (Height = xx,
+ Sex = "Male"),
+ interval = "confidence")
>
> plot(fev$FEV ~ fev$Height,
+ main =  "(Expected value of)   FEV by Height and Sex",
+ xlab =  "Height   [cm]",   ylab = "FEV   [1]",   type = "n")
> with(fev,
+ points(FEV[Sex ==  "Female"]   ~  Height[Sex == "Female"],
+ col = "lightpink"))
> with(fev,
+ points(FEV[Sex ==  "Male"]   ~  Height[Sex == "Male"],
5.6. ODHADY STŘEDNÍ HODNOTY A PREDIKCE NA ZÁKLADĚ MODELU
68
+ col = "cornflowerblue"))
> lines(xx,   yyl.conf[,   1],   col =  "lightpink2",   lwd = 2)
> # estimated expected value of FEV by Height for girls
> lines(xx,  yy2.conf[,   1],  col = "dodgerblue2",   lwd = 2)
> # estimated expected value of FEV by Height for boys
> lines(xx,  yyl.conf[,   2],  col = "lightpink2",   lwd = 2,   lty = 2)
> # lower limits of the confidence intervals for girls
> lines(xx,  yy2.conf[,   2],  col = "dodgerblue2",   lwd = 2,   lty = 2)
> # lower limits of the confidence intervals for boys
> lines(xx,  yyl.conf[,   3],  col = "lightpink2",   lwd = 2,   lty = 2)
> # upper limits of the confidence intervals for girls
> lines(xx,  yy2.conf[,   3],  col = "dodgerblue2",   lwd = 2,   lty = 2)
(Expected value of) FEV by Height and Sex
120      130      140      150      160      170      180 190
Height [cm]
> # upper limits of the confidence intervals for boys Stejně tak přidání predikčních intervalů.
> # prediction intervals for the fitted lines
> yyl.pred <- predict(mod,
+ newdata = data.frame(Height = xx,
+ Sex = "Female"),
5.6. ODHADY STŘEDNÍ HODNOTY A PREDIKCE NA ZÁKLADĚ MODELU
69
+ interval = "prediction")
> yy2.pred <- predict(mod,
+ newdata = data.frame(Height = xx,
+ Sex = "Male"),
+ interval = "prediction")
>
> plot(fev$FEV ~ fev$Height,
+ main =  "(Prediction of)   FEV by Height and Sex",
+ xlab =  "Height   [cm]",   ylab = "FEV   [1]",   type = "n")
> with(fev,
+ points(FEV[Sex ==  "Female"]   ~  Height[Sex == "Female"],
+ col = "lightpink"))
> with(fev,
+ points(FEV[Sex ==  "Male"]   ~  Height[Sex == "Male"],
+ col = "cornflowerblue"))
> lines(xx,  yyl.pred[,   1],  col = "lightpink2",   lwd = 2)
> # prediction of FEV by Height for girls
> lines(xx,  yy2.pred[,   1],  col = "dodgerblue2",   lwd = 2)
> # prediction of FEV by Height for boys
> lines(xx,  yyl.pred[,   2],  col = "lightpink2",   lwd = 2,   lty = 2)
> # lower limits of the predict ion intervals for girls
> lines(xx,  yy2.pred[,   2],  col = "dodgerblue2",   lwd = 2,   lty = 2)
> # lower limits of the predict ion intervals for boys
> lines(xx,  yyl.pred[,   3],  col = "lightpink2",   lwd = 2,   lty = 2)
> # upper limits of the predict ion intervals for girls
> lines(xx,  yy2.pred[,   3],  col = "dodgerblue2",   lwd = 2,   lty = 2)
5.6. ODHADY STŘEDNÍ HODNOTY A PREDIKCE NA ZÁKLADĚ MODELU
70
> # upper limits of the prediction intervals for boys
Konfidenční intervaly jsou poměrně úzké díky velkému rozsahu výběru (n = 654). To neplatí pro predikční intervaly, v jejichž šířce hraje významnou roli rozptyl a2 náhodné chyby e, který se s rozsahem dat nemění.
Vykreslené konfidenční (i predikční) intervaly jsou bodové, jejich (1 — a) 100% pokrytí je tedy zaručeno pro každou hodnotu výšky Height a pohlaví Sex zvláší. Společné pokrytí by zaručovaly konfidenční pásy. Konzervativní konfidenční pásy můžeme odvodit ze Scheffého věty a poté spočítat a vykreslit v Q. Podle Scheffého věty je pravděpodobnost, že
{bT(A3 - A/3)}2 <mFi_a(m,n - p) a2 bTA (XTX)_1ATb
rovna l — a pro všechna b G W11 současně, je-li A matice typu mxp plné hodnosti. V modelu (5.) je vektor koeficientů (fi0, fii, fii, fis)T • Přímka popisující závislost střední hodnoty FE V na výšce Height má pro dívky předpis y = fi0 + fiix a pro chlapce předpis y = (fi0 + ^2) + (fii + /%) x. Pro konfidenční pásy pro přímku pro dívky můžeme v Scheffého větě použít
A = (o   1  0  0) ' b=(x)'xGR' zatímco pro konfidenční pásy pro přímku pro chlapce můžeme použít
5.6. ODHADY STŘEDNÍ HODNOTY A PREDIKCE NA ZÁKLADĚ MODELU
71
Pokrytí je pak alespoň 100 x (1 — a) % pro celou přímku současně (pro každou přímku zvlášť). Konzervativní konfidenční pásy pro obě přímky mají tvar
x'/3- y2F1_a(2,n-p) a2 xT(XTX)-1x
x'/3+ \J2F1-a{2,n-p) a2 xT(XTX)-xx
kde x = (l,x,0,0)T, x E IR, pro dívky a x = (l,x, l,x)T, x E IR, pro chlapce. Konfidenční pásy pro přímku tedy vypadají podobně jako konfidenční intervaly: místo násobení kvantilem
ti-a/2(n — p) se použije násobení a/2 Fi_a(2, n — p).
> # confidence bands for the fitted lines
> df2 <- summary(mod)$df[2]
> yyl.conf.b <- yyl.conf[,   1]   +   (yyl.conf - yyl.conf[,   1]) / +      qt(l  - alpha / 2,   df = df2) *
+      sqrt(2  * qf(l  - alpha,   dfl = 2,   df2 = df2))
> yy2.conf.b <- yy2.conf[,   1]   +   (yy2.conf - yy2.conf[,   1]) / +      qt(l  - alpha / 2,   df = df2) *
+      sqrt(2  * qf(l  - alpha,   dfl = 2,   df2 = df2))
>
> plot(fev$FEV ~ fev$Height,
+ main =  "(Expected value of)   FEV by Height and Sex",
+ xlab =  "Height   [cm]",   ylab = "FEV   [1]",   type = "n")
> with(fev,
+ points(FEV[Sex ==  "Female"]   ~  Height[Sex == "Female"],
+ col = "lightpink"))
> with(fev,
+ points(FEV[Sex ==  "Male"]   ~  Height[Sex == "Male"],
+ col = "cornflowerblue"))
> lines(xx,   yyl.conf[,   1],   col =  "lightpink2",   lwd = 2)
> # estimated expected value of FEV by Height for girls
> lines(xx,  yy2.conf[,   1],  col = "dodgerblue2",   lwd = 2)
> # estimated expected value of FEV by Height for boys
> lines(xx,  yyl.conf.b[,   2],   col = "lightpink2", + lwd = 2,   lty = 2)
> # lower limits of the confidence bands for girls
> lines(xx,  yy2.conf.b[,   2],   col = "dodgerblue2", + lwd = 2,   lty = 2)
> # lower limits of the confidence bands for boys
> lines(xx,  yyl.conf.b[,   3],   col = "lightpink2", + lwd = 2,   lty = 2)
> # upper limits of the confidence bands for girls
5.6. ODHADY STŘEDNÍ HODNOTY A PREDIKCE NA ZÁKLADĚ MODELU
72
> lines(xx,  yy2.conf.b[,   3],   col = "dodgerblue2", + lwd = 2,   lty = 2)
(Expected value of) FEV by Height and Sex
>
LU
120      130      140      150      160      170      180 190
Height [cm]
> # upper limits of the confidence bands for boys Jelikož rozdíl mezi ti_a/2{n — p) a a/2 Fi-a(2, n — p) není velký,
> qt(l  - alpha/2,   df = df2) [1] 1.96362
> sqrt(2  * qf(l  - alpha,   dfl = 2,   df2 = df2)) [1] 2.453398
nejsou konfidenční pásy výrazně širší než konfidenční intervaly. K vykreslení obrázků můžeme použít také knihovnu ggplot 2.
> newdat <- expand.grid(Height = xx,
+ Sex = factor(c("Female", "Male")))
> newdat <- cbind(newdat,   FEV = NA,   Lw = NA,   Up = NA)
5.6. ODHADY STŘEDNÍ HODNOTY A PREDIKCE NA ZÁKLADĚ MODELU
73
> newdat[,  c(3:5)]   <- rbind(yyl.conf, yy2.conf)
>
> ggplot(data = fev,
+ mapping = aes(x = Height,  y = FEV,  color = Sex)) +
+ geom__point () +
+ geom_ribbon(data = newdat,
+ aes(ymin = Lw,  ymax = Up,
+ fill = Sex,   color = NULL),
+ alpha =  .15) +
+ geom_line(data = newdat) +
+ xlab("Height   [cm]")   + ylab("FEV  [1]") +
+ ggtitle("(Expected value of)   FEV by Height and Sex") +
+ theme(plot.title = element_text(hjust = 0.5))
(Expected value of) FEV by Height and Sex
6-
130 150 170 190
Height [cm]
Kapitola 6 Výběr modelu
V předchozí kapitole jsme se zabývali inferencí v daném lineárním modelu. Nemáme-li model předem daný, zpravidla je se součástí analýzy dat i výběr vhodného modelu, což je v lineárních modelech ekvivalentní výběru vhodné regresní matice. Výběr modelu je komplexní úloha, na kterou se můžeme dívat z několika úhlů. Některé aspekty jsou přitom důležitější než jiné v závislosti na tom, k čemu má výsledný model sloužit. Jediný správný model obvykle neexistuje, proto nebývá užitečné jej hledat. Užitečnější bývá pohlížet na modely jako na aproximace skutečnosti a hledat takovou, která nám umožní zodpovědět otázku, kvůli které data analyzujeme. Na zodpovězení různých otázek je možné použít i různé modely. Závěry, které na jejich základě můžeme udělat, by ale neměly být kvalitativně odlišné. Dospějeme-li k závěru, že data podporují modely s kvalitativně odlišnými závěry, je vhodné zvážit, zda je na základě takových dat možné zodpovědět na položené otázky.
V následujících částech si popíšeme několik faktorů, které mohou do výběru modelu vstoupit.
6.1   Srovnání vnořených modelů pomocí testování
V části 5.5 jsme se věnovali F—testu hypotézy o několika lineárních kombinacích složek vektoru koeficientů (3. Nyní se na tento problém podíváme jako na problém testování možnosti přechodu od většího modelu k menšímu. Připomeňme si, že o modelech Y = Xb/3b + e a Y = XS/3S + e mluvíme jako o vnořených, je-li lineární obal sloupců regresní matice menšího modelu podpro-storem lineárního obalu sloupců regresní matice většího modelu. O menším modelu pak mluvíme jako o podmodelu většího modelu. Tento vztah mezi lineárními obaly sloupců regresních matic nastane, když jsou sloupce regresní matice podmodelu lineárními kombinacemi sloupců regresní matice modelu. Platnost podmodelu pak implikuje konkrétní hodnoty pro konkrétní lineární kombinace příslušných koeficientů v modelu. Na testování možnosti přechodu od modelu k podmodelu se proto můžeme dívat jako na testování hypotézy o těchto lineárních kombinacích koeficientů v modelu.
V kontextu testování možnosti přechodu od modelu k podmodelu se testová statistika (5.3)
6.1. SROVNÁNÍ VNOŘENÝCH MODELŮ POMOCÍ TESTOVÁNÍ
75
používá spíše ve tvaru
(lle	•s 1	12	—	|eb|	\2)/m
	|eb|		|2/(ri-		-P)
(6.1)
kde eb = Y — Xb/3b je vektor reziduí v modelu a es = Y — XS/3S je vektor reziduí v podmodelu. Za platnosti podmodelu normálního lineárního modelu má testová statistika F—rozdělení s m a n — p stupni volnosti, kde m je rozdíl mezi dimenzí lineárního obalu sloupců regresní matice modelu a dimenzí lineárního obalu sloupců regresní matice podmodelu. Mají-li obě regresní matice plnou hodnost, jedná se o rozdíl mezi počtem koeficientů v modelu a v podmodelu. Vraťme se nyní ke srovnání modelů
(5.) FEVj = (30 + Ä x Height^ + (32 x Sex^ + (33 x (Sex x Height). + eh   i = 1,... ,n
a
(3.) FEVí = (30 + Ä x Sexí + eh   i = l,...,n,
které jsme již v části 5.5 provedli pomocí testové statistiky (5.3). Test o možnosti přechodu od modelu k podmodelu provedeme v ^ pomocí funkce anova () .
> mod.b <- Im(FEV ~  Height  *  Sex,   data = fev)
> mod.s <- Im(FEV ~  Sex,   data = fev)
>
> anova(mod.s,  mod.b)
Analysis of Variance Table
Model  1:  FEV ~ Sex
Model 2:  FEV ~  Height  * Sex
Res.Df        RSS Df Sum of Sq F Pr(>F)
1 652 469.60
2 650  114.88    2        354.71  1003.5 < 2.2e-16 *** Signif.  codes:     0   '***'   0.001   '**'   0.01   '*'   0.05   '.' 0.1
Hodnota testové statistiky (sloupec F ve výstupu výše) a p—hodnoty (sloupec P r (>F)) jsou, podle očekávání, stejné (až na zaokrouhlení) jako v části 5.5. Zbytek výstupu obsahuje jednotlivé stavební kameny výpočtu testové statistiky podle vzorce (6.1): normy reziduí jednotlivých modelů eb a es (sloupec RSS), jejich rozdíl (sloupec Sum of Sq.), rozdíl m v počtu koeficientů modelu a podmodelu (sloupec D f) a stupně volnosti %2—rozdělení odhadů rozptylů náhodných chyb v jednotlivých modelech (sloupec Re s . D f).
Pro korektní provedení testu je nutné, aby oba modely, model. b i model. s, byly v odhadnuty na stejných datech. Někdy se stane, že v datech pro některá pozorování chybí informace o některých prediktorech. Taková pozorování O nepoužije k odhadu modelu. Stane-li se u některých pozorování, že chybějící hodnoty se objeví v prediktorech modelu, ale už ne
6.2. SROVNÁNÍ MODELŮ POMOCÍ KRITÉRIÍ
76
v prediktorech podmodelu, budou tato pozorování použita k odhadu podmodelu ale ne k odhadu modelu. V takovém případě Cg na dotaz ano va (mod. s, mod. b) odpoví hlášením chyby. Té se můžeme vyhnout tak, že z dat, ze kterých se má odhadnout podmodel, vynecháme pozorování, u kterých chybí informace o prediktorech modelu, například pomocí funkce na . omit () .
6.2 Srovnání modelů pomocí kritérií
Kromě testování lze modely srovnávat i pomocí kritérií, jako jsou adjustovaný koeficient determinace, kterému jsme se věnovali v části 5.2, anebo Akaikeho či bayesovské informační kritérium. Na rozdíl od testování lze pomocí kritérií srovnávat i modely, které nejsou vnořené. Stále ale platí, že srovnání má smysl jenom tehdy, byly-li modely odhadnuty na stejných datech. Na tuto podmínku si v případě kritérií musíme dávat větší pozor než u testování, jelikož <S o kritéria můžeme žádat pro každý model zvláší. CJt pak „nevf, které modely chceme srovnávat, a nemůže nás upozornit na případný problém s nestejnými daty.
Z části 5.2 víme, že adjustovaný koeficient determinace (5.1) pro model model získáme pomocí příkazu summary (model) $ad j . r . squared a že větší hodnota kritéria znamená lepší model (z hlediska tohoto kritéria).
Ve vztahu (5.1) pro adjustovaný koeficient determinace je kvalita modelu reprezentována vzdáleností predikce na základě modelu od pozorované odezvy ||e||2 = XT=i(^ ~ Yi)2■ Tato kvantita se nutně zvýší, kdykoliv do modelu přidáme prediktor, který není lineární kombinací prediktorů, které již v modelu jsou. Proto je ve vztahu pro adjustovaný koeficient determinace tato vzdálenost penalizována počtem p koeficientů modelu (dělí se n — p).
Na zohlednění odhadnuté přesnosti predikce a velikosti modelu je postaveno také Akaikeho a bayesovské informační kritérium, (AIC) a BIC. Odhad přesnosti predikce je v obou případech reprezentován maximalizovanou věrohodností a penalizace má tvar násobku počtu všech parametrů modelu. V případě AIC jde o dvojnásobek počtu parametrů, zatímco v případě BIC násobíme počet parametrů logaritmem počtu pozorování: penalizace je tedy vyšší pro vyšší rozsahy dat. Znaménka jsou u obou kritérií nastavena tak, že menší hodnota kritéria znamená lepší model (z hlediska daného kritéria). AIC mívá tendenci vybírat spíše větší modely a doporučuje se k němu přihlížet spíše v situacích, kdy nás zajímá predikce na základě modelu. BIC, naopak, častěji preferuje menší modely a lze jej využít v situaci, kdy nás spíše zajímá, které prediktory opravdu souvisejí se střední hodnotou odezvy. Obě kritéria lze použít i pro jiné než lineární modely. O Akaikeho a bayesovské informační kriterium pro model model zažádáme v <© pomocí příkazů AIC (model)  a BIC (model) .
6.3 Diagnostika modelu
Inference na základě modelu je spolehlivá jenom tehdy, když model platí, tj. když jsou splněny jeho předpoklady. Připomeňme si, že u lineárního modelu předpokládáme, že Y = X/3 + e, náhodná chyba e má nulovou střední hodnotu, Es = 0, a rozptyl Var e = a2I. V normálním lineárním modelu navíc požadujeme normalitu náhodné chyby.
6.3. DIAGNOSTIKA MODELU
77
Splnění předpokladů lineárního modelu nelze ověřovat přímo, protože (3 neznáme a e nepozorujeme. Proto při ověřování předpokladů místo rovnice Y = X/3 + e pracujeme s rovnicí Y = X/3 + e, která obsahuje známé složky. Jelikož lze předpoklady lineárního modelu formulovat na náhodnou chybu e, jejich splnění ověřujeme především na vektoru reziduí e.
U ověřování předpokladů je výhodné pracovat s diagnostickými grafy. Na nich můžeme vidět, zda se nějaký předpoklad modelu zdá být porušen, a také jakým způsobem, což nás může navést na vhodné řešení problému. Nejdůležitější diagnostické grafy vrací v <3 funkce plot () uplatněna na lineární model, u kterého chceme splnění předpokladů ověřovat.
> plot(mod)
Fitted values lm(FEV ~ Height * Sex)
6.3. diagnostika modelu
78
Normal Q-Q
Theoretical Quantiles lm(FEV ~ Height * Sex)
Scale-Location
6240
Fitted values lm(FEV ~ Height * Sex)
6.3. DIAGNOSTIKA MODELU
79
Leverage lm(FEV ~ Height * Sex)
První z diagnostických grafů vykresluje rezidua e proti odhadnutým středním hodnotám Y. Červená křivka je neparametrický odhad závislosti e na Y. Z definice jsou tyto vektory na sebe kolmé a za platnosti normálního lineárního modelu jde o realizace nezávislých náhodných vektorů. Na grafu by proto neměly být patrné žádné trendy. Různé trendy, které na tomto grafu může být vidět, mohou naznačovat porušení různých předpokladů modelu.
Nej důležitějším (pro platnost inference) předpokladem lineárního modeluje nulovost střední hodnoty náhodné chyby e nebo ekvivalentně správný tvar regresní matice X. Za platnosti tohoto předpokladu mají i rezidua e nulovou střední hodnotu. Je-li tento předpoklad, naopak, porušen, rezidua nulovou střední hodnotu mít nemusí a na grafu reziduí proti odhadnutým středním hodnotám se mohou objevit nějaké trendy. Červená křivka v prvním grafu se v takovém případě může výrazněji odchylovat od nuly.
Máme-li podezření, že jsme v modelu nezachytili závislost střední hodnoty odezvy na nějakém (potenciálním) prediktoru anebo že jsme nezvolili správnou formu závislosti na nějakém prediktoru (například je potřeba transformace prediktoru), můžeme si rezidua vykreslit také proti tomuto prediktoru. Tvar závislosti, který na grafu případně uvidíme, můžeme poté zkusit zahrnout do modelu. Po každém takovém rozšíření modelu je potřeba znovu zkontrolovat diagnostiku modelu.
Dalším problémem, kterého si případně můžeme všimnout již na prvním grafu, je závislost rozptýlení reziduí na odhadnutých středních hodnotách Y. To naznačuje problém s předpokladem o stejných rozptylech náhodných chyb pro jednotlivá pozorování (homoskedasticitě). Na ověřování tohoto předpokladu je zaměřený třetí graf. Ten místo reziduí e používá standardizovaná
6.3. DIAGNOSTIKA MODELU
80
rezidua
T' = - —
\J^{\-hu)
kde h,hi jsou diagonální prvky projekční matice H = X(XTX)~1XT. Standardizovaná rezidua r i mají, na rozdíl od reziduí ei5 za platnosti modelu stejné (jednotkové) rozptyly. Standardizovaná rezidua se na třetím grafu vykreslí proti odhadnutým středním hodnotám Y. Pro zvýraznění případné závislosti se na y—novou osu vykresluje odmocnina z absolutní hodnoty reziduí a graf se prokládá neparametrickým odhadem jejich závislosti na Y (červená křivka). Výraznější odchylky od (kladné) horizontální přímky naznačují problémy s homoskedasticitou. Podobně jako u střední hodnoty můžeme i u rozptylu vyhodnotit případnou závislost na prediktorech pomocí grafů standardizovaných reziduí proti prediktorům.
Vedle homoskedasticity předpokládá lineární model i nekorelováno st náhodných chyb poslouchajících jednotlivým pozorováním. Porušení tohoto předpokladu často vyžaduje využití pokročilejšího modelu, například časové řady, plyne-li závislost z časové následnosti mezi pozorováními, anebo modelu s náhodnými efekty, je-li závislost dána opakovanými pozorováními na stejných objektech. Detailní diagnostika závislosti se pak provede nástroji obvyklými v kontextu těchto pokročilejších modelů.
Druhý graf slouží k ověření normality náhodných chyb. I tady se používají spíše standardizovaná rezidua, která mají za platnosti modelu přibližně stejné normální rozdělení. Jedná se o klasický Q-Q graf, a proto by se za platnosti předpokladů body neměly výrazně odchylovat od proložené přímky, vyskytovat se od ní jenom na jednu stranu (problémy se šikmostí rozdělení) anebo ve tvaru písmene S (problémy se špičatostí rozdělení). Mírná odchýlení od proložené přímky jsou ale vpořádku. Názorné ukázky Q-Q grafů pro data různých rozsahů pocházející z normálního i z jiných rozdělení je možné nalézt například v příloze bakalářské práce [4].
Poslední graf slouží k detekci vlivných pozorování. Leverages hi:i vykresleny na x—ové ose kvantifikují potenciál pozorování ovlivnit odhad modelu. Cookova vzdálenost, která je do grafu vnesena pomocí vrstevnic, kvantifikuje, jestli k ovlivnění odhadu opravdu dochází. U pozorování s vysokými hodnotami těchto ukazatelů (u leverages se, v závislosti od zdroje, uvádí p/n, 2p/n, anebo 3p/n, kde p je počet koeficientů v modelu a n je počet pozorování; u Cookovy vzdálenosti 0.5 anebo 1) je potřeba ověřit, zda jsou hodnoty odezvy i prediktorů správně zadané, a rozmyslet si, zda jsou tato pozorování reprezentativní pro populaci, na kterou chceme výsledky modelu zobecňovat. Pro podrobnější diagnostiku vlivných pozorování je v možné si v <S vyžádat rozšířenou nabídku diagnostických grafů.
> plot(mod, which
= c(l:6) )
6.3. diagnostika modelu
81
Residuals vs Fitted
Fitted values lm(FEV ~ Height * Sex)
Normal Q-Q
-3-2-10 1 2 3
Theoretical Quantiles lm(FEV ~ Height * Sex)
6.3. diagnostika modelu
82
Scale-Location
624Ô
Fitted values lm(FEV ~ Height * Sex)
Cook's distance
473
íLjáAlňÁÁíiáiíÁili
624
0
100
200
300
400
500
600
Obs. number lm(FEV ~ Height * Sex)
6.3. DIAGNOSTIKA MODELU
83
Residuals vs Leverage
i-r
0.000 0.005
0.010
i-r
0.015 0.020
0.025 0.030
Leverage lm(FEV ~ Height * Sex)
cd o
Ctí -t—*
w T3 </)
o o
O
o
CO
o
CvJ
o
o o
o o
Cook's dist vs Leverage hN/(1 -hN)
5o624 4'
/ 02
3
• ' 473o
Or'
i / i /
O
o,-
'/ .<&£Z$~o « o o
§o °o
iHi9QPCr-8 — #---u©--
.a
T
T
T
T
0
0.005 0.01
0.015 0.02
—I— 0.025
Leverage hN lm(FEV ~ Height * Sex)
Na diagnostických grafech pro model
0.03
(5.) FEVj = (30 + Ä x Height^ + (32 x Sex^ + (33 x (Height x Sex)j + eh   i = 1,..., n
6.3. DIAGNOSTIKA MODELU
84
vidíme problémy již s volbou regresní matice. Tvar závislosti na prvním diagnostickém grafu naznačuje nepodchycení nějaké kvadratické závislosti. Podobný trend najdeme na grafu reziuí e proti prediktoru výška Height.
> e <- residuals(mod)
> plot(e ~  fev$Height,  main = "Residuals vs Height", + xlab = "Height",  ylab = "Residuals")
> abline(h = 0,   col = "grey",   lty = 3)
> lines(lowess(e ~  fev$Height),   col = "red")
Můžeme proto do modelu zkusit přidat kvadratickou závislost střední hodnoty FEV na výšce Height
FEVí = (30 + Pi x Sexí + p2 x Height^ + p3 x (Height x Sex)í+
+ p4 x Height? + p5 x (Height2 x Sex)^ + eh   i = l,...,n.
> mod <- lm(FEV ~  Sex *   (Height + I(Height"2)),   data = fev)
> plot(mod,  which = c(l:6))
6.3. DIAGNOSTIKA MODELU
85
Residuals vs Fitted
cP 080
-ňOft
CD O
624o 648o
0 °
ZOO O @8§8  QOO o
Oooí;0roqoO ^§Qpvn 89 n 0 0
1 80S u°o08o@ g 0 o°
ä 08   0 0
@o o°80o0   o 0
02
T
T
Fitted values lm(FEV ~ Sex * (Height + l(HeightA2)))
Theoretical Quantiles lm(FEV ~ Sex * (Height + l(HeightA2)))
6.3. diagnostika modelu
86
Scale-Location
624o
1 2 3 4 5
Fitted values lm(FEV ~ Sex * (Height + l(HeightA2)))
Cook's distance
473
H 2
i.iL. u—.
.■■■.jjLU.l..j|.Lilil
636
JtuJ. Ul.LL.L-lL
0
T
T
100
200
300
400
500 600
Obs. number lm(FEV ~ Sex * (Height + l(HeightA2)))
6.3. DIAGNOSTIKA MODELU
87
^ -
c\j —
o —
C\J
I
I
Residuals vs Leverage
0.00
0.02
0.04
0.06
0.08
0.10
Leverage
lm(FEV ~ Sex * (Height + l(HeightA2)))
Cook's dist vs Leverage hN/(1 -hN)
cd o
Ctí -t—*
w T3
</)
o o
O
C\J
o
00
o
o
o o
	5/ 4 / / '/ -	3 4730	
♦ * * ♦	'• ' '• / * / /Ô2 ©636 /           *           ^"		
9 /	a-'o	O	
		hr---' ~~0~ i °	o
h-— 0	i i 0.02 0.04	1 1 0.06 0.08	1 0.1
Leverage hN lm(FEV ~ Sex * (Height + l(HeightA2)))
Dianostické grafy rozšířeného modelu již nenaznačují problémy se střední hodnotou náhodných chyb. Naopak zůstal problém s homoskedasticitou a s normalitou dat. Při takto velkém rozsahu
6.4. MUĽTIKOLINEARITA
88
výběru (n = 654) a absenci pozorování s výrazným potenciálem ovlivnit odhad modelu není normalita klíčová, jelikož v takových případech lze platnost inference doložit pomocí centrální limitní věty. Homoskedasticita je ale pro platnost testů a konfidenčních intervalů potřeba. Její nesplnění můžeme řešit například nahrazením lineárního modelu některým z jeho zobecnění. Ty jsou již nad rámec tohoto textu.
Multikolinearita, tedy vztah blízký lineární závislosti mezi prediktory, může působit problémy se stabilitou i statistickou kvalitou odhadů koeficientů i v jinak platných modelech. Závislost mezi dvojicemi prediktorů můžeme odhalit pomocí grafů, na kterých vykreslujeme jeden pre-diktor proti druhému, například pomocí funkce pairs () , kterou jsme používali ve 2. kapitole. Na závislosti nás také může upozornit vysoká výběrová korelace mezi dvojicemi (kvantitativních) prediktorů.
> cor (fev[,  c(l,   3) ] )
Age 1.0000000 0.7919436 Height 0.7919436 1.0000000
Složitější závislosti zahrnující více prediktorů snadněji odhalíme pomocí variance inflation factor VIFj = 1/(1 — R2), kde R2 je koeficient determinace v lineárním modelu, kde je j—tý prediktor z původního modelu odezvou a zbylé prediktory z původního modelu jsou prediktory i v tomto modelu. Hodnoty větší než pět se považují za problematické. Odpovídá-li jednomu prediktorů v původním modelu více koeficientů, použije se místo variance inflation factor VIF jeho zobecnění generalized variance inflation factor gVIF pocházející z lineární regrese s mnohorozměrnou odezvou. VClo tyto charakteristiky požádáme funkcí vif () z knihovny car.
> library(car)
> vif(mod)
6.4 Multikolinearita
Age
Height
Sex
11654.353 Sex:I(Height"2) 15254.133
Height 1137.165
I(Height"2) 1213.311
Sex:Height 52072.429
U proměnných, které jsou z definice propojeny (například interakce a polynomy) lze, pochopitelně, závislost očekávat.
6.4. MUĽTIKOLINEARITA
89
> vif(lm(FEV ~  Sex + Height,  data = fev))
Sex Height 1.025946 1.025946
> vif(lm(FEV ~  Sex + Age,   data = fev))
Sex Age 1.00085 1.00085
U polynomů lze problém řešit pomocí funkce poly () , která místo klasické parametrizace polynomů využívá ortogonální polynomy daného řádu. Odhady střední hodnoty a tedy i celkové charakteristiky modelu budou stejné, jako kdybychom použili klasický polynom stejného řádu. Sloupce regresní matice odpovídající prediktoru, na kterém má střední hodnota odezvy záviset polynomiálně, jsou ale na sebe kolmé, což zaručuje stabilitu odhadu.
> summary(lm(FEV ~ Height + I(Height"2),  data = fev))
Call:
lm(formula = FEV ~  Height + I(Height~2),   data = fev) Residuais:
Min 1Q      Median 3Q Max
-1.80103 -0.22928 -0.00255    0.21924 1.99699
Coefficients:
Estimate Std.  Error t value Pr(>|t|)
(Intercept)     6.027e+00 1.503e+00 4.009 6.79e-05 ***
Height            -9.844e-02 1.963e-02 -5.015 6.83e-07 ***
I(Height"2)     4.891e-04 6.372e-05 7.675 6.07e-14 ***
Signif.  codes:     0   •***•   0.001   '**'   0.01   '*'   0.05   '.'   0.1   '   ' 1
Residual standard error:   0.4127 on 651 degrees of freedom Multiple R-squared:     0.7741,Adjusted R-squared: 0.7734 F-statistic:     1115 on 2  and 651 DF,    p-value:   < 2.2e-16
> summary(Im(FEV ~ poly(Height,   2),  data = fev))
Call:
lm (formula = FEV ~ poly (Height,   2),   data = fev)
6.4. MULTIKOLINEARITA
90
Residuals:
Min IQ Median
-1.80103 -0.22928 -0.00255
3Q
0 .21924
Max 1.99699
Coefficients:
(Intercept)
Estimate Std. Error 2.63678 0.01614
t value Pr (>11I) 163.376    < 2e-16 ***
poly(Height, 2)1 19.23502 0.41274 46.604 < 2e-16 *** poly(Height,   2)2    3.16778        0.41274      7.675  6.07e-14 ***
Signif.  codes:     0   •***•   0.001   '**'   0.01   '*'   0.05   '.'   0.1   '   ' 1
Residual standard error:   0.4127 on 651 degrees of freedom Multiple R-squared:     0.7741,Adjusted R-squared: 0.7734 F-statistic:     1115 on 2  and 651 DF,    p-value:   < 2.2e-16
Funkce vif () chápe ortogonální polynom jako jediný prediktor, kterému odpovídá více koeficientů.
> vif(lm(FEV ~ poly(Height,   2)   + Sex,   data = fev))
GVIF Df GVIF"(1/(2*Df)) poly(Height, 2) 1.092802 2 1.022434 Sex 1.092802    1 1.045372
Kapitola 7
Ilustrační analýza časů přeplavání jezera
V této kapitole budeme pomocí lineárních modelů zkoumat, zda čas, za který plavci překonali jezero Ontario, souvisí s jejich věkem a pohlavím. Nejprve si připomeňme data swim z časti 2.1.
> str(swim)
1 data.frame' $ Name $ Sex $ Age
$ Start.Day $ Month $ Year $ Time.min. $ Direction
65 obs.  of    8 variables:
Factor w/ 54 levels  "Angela Kondrak",..:   32 24 7 Factor w/ 2 levels »f","M": 1212221111 num    16 36 28 57.3 41 ... int    8 23 12 23 26 2  17 30  16 22 ... Ord.factor w/ 3 levels  "July"<"Aug"<"Sep":   3 12 int    1954 1956 1956 1957  1957  1961  1974  1974  1975 1976 num    1255 1273 1131  1501  1115   . . .
Factor w/ 3 levels "NS","NSN","SN": 3333333133
5 22 22 11  16 15 1
2 2 3 2 2 2 2
Naší závislou proměnnou bude Time.min. udávající dobu proplavání jezera v minutách. Pro snadnější práci si ji přejmenujeme na Time a podíváme se na ni.
> names(swim)[7]   <- "Time"
> with(swim,
+ boxplot(Time,  xlab="Time [min]",
+ col="lightgreen",   horizontal = TRUE)
+ )
92
-1-1-1-1-1-
1000        1500        2000        2500 3000
Time [min]
Většina plavců jezero překonala do 1 500 minut, zhruba čtvrtina potřebovala více času, dva plavci plavali dokonce 2 500 a i více než 3 000 minut. Tento nejdelší čas ale patří Vicki Keith, která v roce 1987 jezero proplavala obousměrně. Její čas proto nebudeme uvažovat.
> swim <- swim[swim$Time < 3000, ]
> swim$Direction <- factor(swim$Direction)
> with(swim,
+ boxplot(Time,  xlab = "Time [min]",
+ col =  "lightgreen",  horizontal = TRUE)
+ )
93
o
i-1-1 r
1000 1500 2000 2500
Time [min]
Nyní se na data podíváme podrobněji.
> head(swim)
	Name	Sex		Age	Start.Day	Month	Year	Time	Direction
1	Marilyn Bell	F	16.	00000	8	Sep	1954	1255	SN
2	John Jaremey	M	36.	00000	23	July	1956	1273	SN
3	Brenda Fisher	F	28 .	00000	12	Aug	1956	1131	SN
4	Bill Sadlo	M	57 .	31781	23	Aug	1957	1501	SN
5	Jim Woods	M	41.	00000	26	Aug	1957	1115	SN
6	Jim Woods	M	45 .	00000	2	Sep	1961	1027	SN
> summary(swim)
Name		Sex	Age		Start	.Day
Colleen Shields	3	F:38	Min. :14	. 05	Min.	: 1.00
Kim Middleton	3	M: 26	1st Qu.:20	.75	Ist Qu.	: 8.00
Vicki Keith	3		Median :28	.00	Median	:13.50
Jim Woods	2		Mean :31	.11	Mean	:14.14
John Scott	2		3rd Qu.:38	.00	3rd Qu.	: 19 . 75
Kim Lumsdon	2		Max. :66	. 57	Max.	:31.00
(Other)	49					
Month Year Time Direction
94
July	6	Min.	1954	Min.	829
Aug	49	1st Qu.	1979	1st Qu.	1018
Sep	9	Median	1993	Median	1181
		Mean	1992	Mean	1246
		3rd Qu.	2007	3rd Qu.	1407
		Max.	2018	Max.	2461
Můžeme si všimnout, že většina plavců se rozhodla plavat z jihu na sever (Direction SN). > with(swim,
+ boxplot(Time  ~  Direction,   ylab =  "Time [min]",
+ col = c("lightpink",   "cornflowerblue"),
+ border = c("lightpink2",   "dodgerblue2"),
+ pen = 19)
+ )
Direction
> ggplot(data = swim, aes(x = +      geom_histogram(position =
+ aes(fill =
+ bins = 10,
+      xlab("Time   [min]")
Time,   color = Direction)) + "identity",
Direction, color = Direction), alpha =0.5) +
95
15-
10-
o o
Direction
I NS
SN
1000
1500 2000
Time [min]
2500
Tento směr se opravdu jeví jako rychlejší. V obráceném směru plavalo jenom pět plavců, odhad vlivu směru na délku plavání by tak nemusel být spolehlivý. Nás primárně zajímá závislost délky plavání na věku a pohlaví, proto se omezíme jenom na data o plavcích, kteří plavali z jihu na sever.
> swim <- swim[swim$Direction
> summary(swim)
"SN", -8]
Name		Sex	Age			Start	Day
Colleen Shields	3	F:33	Min.	14	. 05	Min.	1.00
Jim Woods	2	M: 26	1st Qu.	19	.50	Ist Qu.	8 .00
John Scott	2		Median	28	.00	Median	12 .00
Kim Lumsdon	2		Mean	31	.28	Mean	13.73
Vicki Keith	2		3rd Qu.	39	.00	3rd Qu.	18.50
Angela Kondrak	1		Max.	66	. 57	Max.	31.00
(Other)	47						
Month
July
Aug
Sep
5 46
Min. 1st Qu. Median Mean 3rd Qu. Max.
Year
1954 1978 1993 1993 2008 2018
Time
Min. Ist Qu. Median Mean 3rd Qu. Max.
82 98 116 122 138 246
96
> with(swim,
+ boxplot(Time,  xlab = "Time   [min]",
+ col =  "lightgreen",  horizontal = TRUE)
+ )
O
-1-1-1-r
1000 1500 2000 2500
Time [min]
Proměnné Month a Start. Day společně obsahují informaci o datu, kdy plavec zahájil úspěšný pokus o přeplavání jezera. Datum můžeme také uvažovat jako kvantitativní prediktor, například tak, že budeme uvažovat, kolik dnů uběhlo od začátku července v den, kdy plavec zahájil svůj pokus.
> swim$Date <- NA
> swim$Date[swim$Month ==  "July"]   <- 0
> swim$Date[swim$Month ==  "Aug"]   <- 31
> swim$Date[swim$Month ==  "Sep"]   <- 62
> swim$Date <- swim$Date + swim$Start.Day
Nyní se můžeme podívat na souvislost mezi časem Time a datem plavby Date.
> with(swim,  plot(Time  ~  Date,   ylab =  "Time [min]", + col="lightgreen"))
> with(swim,  lines(lowess(Time ~ Date),   col = "palegreen3"))
abline(v = 31.5, col = "grey", lty = 3) abline(v = 62.5,  col = "grey",   lty = 3)
o
C\J
o o o c\j
o o
o o o
July
Aug Month
Sep
ggplot(data = swim,  aes(x = geom_histogram(position =
aes(fill = bins = 10, xlab("Time [min]")
Time, color = Month)) + "identity",
Month,   color = Month), alpha =0.5) +
99
10-
o o
Month I July I Aug
Sep
1000
1500 2000 Time [min]
2500
Na grafech výše vidíme, že nejvíce úspěšných pokusů o překonání jezera se uskutečnilo v srpnu, což může souviset s vhodnějším počasím. Časy přeplavání jezera v červenci a v září, kdy bylo úspěšných pokusů o poznání méně, vypadají rozptýlenější než v srpnové časy. Červencové časy navíc vypadají delší. Z podobných důvodů jako u směru plavání může i u data plavání být rozumné omezit se u analýzy vlivu věku a pohlaví na čas přeplavání jezera na srpnové pokusy, kdy časy plavání nejsou natolik ovlivněny nesouvisejícími faktory, na odhadnutí jejichž efektů nemáme dostatek dat.
> swim <- swim[swim$Month
> summary(swim)
"Aug",   -c(4,   5, 8)]
Name		Sex
Colleen Shields	3	F:26
John Scott	2	M: 20
Kim Lumsdon	2	
Angela Kondrak	1	
Annaleise Carr	1	
Ashleigh Beacham	1	
(Other)	36	
Time Min.       : 829 1st Qu.: 976 Median :1142
Age Year
Min.	14 .	05	Min.	1956
1st Qu.	19.	50	1st Qu.	1980
Median	30 .	00	Median	1995
Mean	31.	66	Mean	1994
3rd Qu.	39.	50	3rd Qu.	2008
Max.	66.	57	Max.	2017
100
Mean :1187 3rd Qu.:1356 Max. :1626
Nyní se podíváme na vztahy mezi zbylými proměnnými. > ggpairs(swim[, -1])
Sex
Age	
1 -u_	
Year						
						
						
						
						
						
						
						
Time
Corr: 0 .187
1500 1250 1000
•    • •
«••->.   '  * * •   #•   • •
0   1   2  30   1   2  3    20 30 40 50 60
C/3 CD
Corr:		Corr:	>
			
0 .154		0.208	CD
			
O)
1960 1980 2000 20800 1000120014001600
Grafy nenaznačují žádné zřejmé problémy s daty ani žádné výrazné závislosti. Nejspíše se tedy nemusíme obávat multikolinearity, ale také nemůžeme očekávat vysvětlení podstatnější části variability časů přeplavání jezera pomocí věku a pohlaví plavce anebo roku, ve kterém svůj pokus uskutečnil.
Podíváme se ještě podrobněji na závislost mezi plánovanou odezvou a prediktory.
> with(swim,
+ boxplot(Time ~  Sex,  ylab = "Time   [min]",  xlab =
+ col = c("lightpink",   "cornflowerblue")
+ border = c("lightpink2",   "dodgerblue2"
+ pen = 19)
+ )
101
> +
>
with(swim, with(swim,
plot(Time  ~ Age,   ylab = "Time
col="lightgreen")) lines(lowess(Time ~ Age), col
[min]",
= "palegreen3"))
102
> with(swim,  plot(Time  ~  Year,   ylab =  "Time [min]", + col="lightgreen"))
> with(swim,  lines(lowess(Time ~ Year),   col = "palegreen3"))
103
cd
E
o o
o o o
o o
00
		0 o
o	0             o o 8 °	o o o o o
	O      ^-----—_____ o       o < --____—-"-^ o	3 O ° o o ° °
	o o °   0 oo	cP o °
1 1960	I         I         I I 1970      1980      1990 2000	I 2010
Year
Grafy nenaznačují žádné složité závislosti, proto se podíváme na jednoduchý model
Timej = /30 + Ä x Sex^ + /32 x Age^ + /33 x Yearj + eh i
1,.. .,n.
> model.1  <- lm(Time  ~  Sex + Age + Year,   data = swim)
> summary(model.1)
Call:
lm(formula = Time  ~  Sex + Age + Year,   data = swim)
Residuals:
Min 1Q Median 3Q Max
-355.2  -156.6    -48.1    172.4 491.5
Coefficients:
Estimate (Intercept) -1298.371 SexM -123.696 Age 4.084 Year 1.208
Std.  Error t value Pr(>|t|) 4063.641    -0.320 0.7509 68.047    -1.818 0.0762 2.447      1.669 0.1026 2.042      0.592 0.5571
104
Signif.  codes:     0   •***•   0.001   '**'   0.01   '*'   0.05   '.'   0.1   '   ' 1
Residual standard error:  216.6 on 42 degrees of freedom Multiple R-squared:     0.1356,Adjusted R-squared: 0.07388 F-statistic:   2.197 on 3 and 42 DF,    p-value: 0.1026
Model vysvětluje jen malé procento variability dat, což jsme ale očekávali. Test nulovosti všech koeficientů zároveň vyšel nevýznamný, prediktory tedy buďto nemají velký vliv na střední hodnotu odezvy, anebo data nemají dostatečnou sílu na prokázání tak malých efektů. Zkusíme ještě vynechat nevýznamný prediktor Year, abychom se soustředili na prediktory, které nás zajímají.
> model.2  <- lm(Time  ~  Sex + Age,   data = swim)
> summary(model.2)
Call:
lm(formula = Time  ~  Sex + Age,   data = swim)
Residuals:
Min 1Q    Median 3Q Max
-351.86 -167.51    -66.03    173.03 496.01
Coefficients
(Intercept' SexM Age
Estimate Std.  Error t value Pr(>|t|)
1106.163 81.241    13.616      <2e-16 *■
-133.967 65.298    -2.052      0.0463 *
4.383 2.376      1.845      0.0720 .
Signif.  codes:     0   '***'   0.001   '**' 0.01
0 . 05
0 .1
Residual standard error:  215 on 43 degrees of freedom Multiple R-squared:     0.1284,Adjusted R-squared: 0.08788 F-statistic:  3.168 on 2 and 43 DF,    p-value: 0.05208
I tento model vysvětluje jen malé procento variability dat, adjustovaný koeficient determinace se ale mírně zlepšil a test o nulovosti obou prediktorů současně je nyní na hranici signifikance. Také testy o nulovosti jednotlivých prediktorů dávají přesvědčivější výsledky, což může být docíleno soustředěním veškeré síly na dva prediktory. Abychom se na výsledky testů mohli spolehnout, musí model splňovat předpoklady. Podíváme se proto na diagnostiku modelu.
plot(model.2,  which = c(l:6))
Residuals vs Fitted
	038 57Q 340		
o o o o o	°o?       ° ° ° 0 O o o	o	o
© o o °o	cP   o©      ^~ o o CP    ° o	o o	o
T
T
T
1050    1100    1150    1200    1250    1300    1350 1400
Fitted values lm(Time ~ Sex + Age)
Theoretical Quantiles lm(Time ~ Sex + Age)
Scale-Location
038
57Q
1050    1100    1150    1200    1250    1300    1350 1400
Fitted values lm(Time ~ Sex + Age)
Cook's distance
Obs. number lm(Time ~ Sex + Age)
107
Residuals vs Leverage
cm —\
038
057
0.5
■— OH
C\J
I
0.00
0°0°° o o co
Q> o
Cook's distance
—I— 0.05
480
0.10
0.15
0.20
Leverage lm(Time ~ Sex + Age)
cd o
CC
w T3 </)
o
o O
C\J
o
00
o
C/5 o
o
o o
Cook's dist vs Leverage hN/(1 -hN)
2.5 2 / /'	48Qr' ✓ ✓	1.5		
/o38 osrŕ . i            * *	✓ *			
> / i    / ' i/' '	''/O			
rt * /     t  v, t i      '    ť 'n '    /   ť    ť o				
			o	
. _ At - - - - _£H§^°0 O. _° _ _	\J			
				
1 1 0 0.05	1 0.1	I 0.15		I 0.2
0.5 0
Leverage hN lm(Time ~ Sex + Age)
Na diagnostických grafech nevidíme žádné zásadní problémy, menší problémy se vyskytují jenom u nejrychlejších časů.
108
Můžeme ještě vyzkoušet, zda si nepolepšíme zahrnutím kvadratické závislosti na věku, kterou by mohl naznačovat již graf z deskriptívni fáze analýzy.
> model.3 <- lm(Time ~  Sex + poly(Age,   2),  data = swim)
> summary(model.3)
Call:
lm(formula = Time  ~  Sex + poly(Age,   2),   data = swim)
Residuals:
Min 1Q    Median 3Q Max
-350.03 -168.75    -61.16    166.29 505.89
Coefficients:
Signif.  codes:     0   •***•   0.001   '**'   0.01   '*'   0.05   '.'   0.1   '   ' 1
Residual standard error:  216.8 on 42 degrees of freedom Multiple R-squared:     0.1337,Adjusted R-squared: 0.07184 F-statistic:  2.161 on 3 and 42 DF,    p-value: 0.1069
> plot(model.3,  which = c(l:6))
(Intercept) SexM
poly(Age, 2)1 poly(Age,   2)2
Estimate Std.  Error t value Pr(>|t|) 1238.12 44.98    27.529      <2e-16 ***
-118.30 72.75    -1.626 0.1114
394.32 222.47      1.772 0.0836
121.47 239.50      0.507 0.6147
109
o o
CD
O O C\J
Residuals vs Fitted
038	m			
O o	o 8 °	o° o	o	
o ° O —.	o	o		o
6> 8S	@° o	o o	o o	o
1 1100	I 1200		I 1300	I 1400
o —
o o
I
Fitted values lm(Time ~ Sex + poly(Age, 2))
Normal Q-Q
	...••38o
	.■• 057 ..-034
	.■-"o° ..<aP
o 0ooo.°PP	
o	
1 1 -2 -1	I              I I 0              1 2
C\J —
o —
Theoretical Quantiles lm(Time ~ Sex + poly(Age, 2))
Scale-Location
	038	57o 34o	o	
	O o	S    o o o J>~ ~v J? <p	o o	
	o cP			o
	1 1100	I 1200	I 1300	I 1400
Fitted values lm(Time ~ Sex + poly(Age, 2))
Cook's distance
38
i-1-1-1-r
0 10 20 30 40
Obs. number lm(Time ~ Sex + poly(Age, 2))
57
48
111
_co
CC
■g
C/í CD
T3 cd n
T3
CC T3
CC -I—*
C/)
cd o
cc -t—* c/5
T3
c/í
o o
O
CM —\
= OH
C\J
I
00
o
o
o o
Residuals vs Leverage
	038 570 o ° o <S>     oo^--		o
	o$y° o o Cook's distance	o	
1 1 0.0 0.1		I I 0.2 0.3	I 0.4
lm(Time ~		Leverage Sex + poly(Age, 2))	
Cook's dist vs Leverage hN/(1 -hN)			
2.5572) 1;5
'' oo/'48o
p38/
1
0.5
0.5
0
0.1
0.2
0.3
0.4
Leverage hN lm(Time ~ Sex + poly(Age, 2))
Kvadratický efekt není signifikantní a ani diagnostické grafy se nezlepšily. Zůstaneme proto u lineární závislosti.
Nakonec ještě můžeme prozkoumat potřebu interakce mezi věkem a pohlavím.
> model.4  <- lm(Time  ~  Sex * Age,   data = swim)
> summary(model.4)
Call:
lm(formula = Time  ~  Sex * Age,   data = swim)
Residuals:
Min 1Q    Median 3Q Max
-304.27 -139.56    -76.11    158.36 489.78
Coefficients:
Estimate Std.  Error t value Pr(>|t|)
(Intercept)   1162.319 90.578 12.832 4e-16 ***
SexM                -373.270 188.598 -1.979 0.0544 .
Age                        2.462 2.750 0.895 0.3757
SexM:Age              7.182 5.317 1.351 0.1840
Signif.  codes:     0   '***'   0.001   '**'   0.01   '*'   0.05   '.' 0.1
Residual standard error:  212.9 on 42 degrees of freedom Multiple R-squared:     0.1647,Adjusted R-squared: 0.105 F-statistic:     2.76 on 3 and 42 DF,    p-value: 0.05391
> plot(model.4, which
= c(l:6) )
Residuals vs Fitted
038
1000 1100 1200 1300
Fitted values lm(Time ~ Sex * Age)
Normal Q-Q
Theoretical Quantiles lm(Time ~ Sex * Age)
114
Scale-Location
o
1000 1100 1200 1300
Fitted values lm(Time ~ Sex * Age)
CvJ
o
00
o
o
o o
Cook's distance
Obs. number lm(Time ~ Sex * Age)
115
Residuals vs Leverage
Cook's distance °48
S-1-1-1-1-T~
0.00       0.05       0.10       0.15       0.20 0.25
Leverage lm(Time ~ Sex * Age)
CD
o
CC -t—* CO
T3
C/í
O
o O
00
o
o
o o
Cook's dist vs Leverage hN/(1 -hN)
0
0.05 0.1
0.15
0.2
0.25
Leverage hN lm(Time ~ Sex * Age)
Ani toto rozšíření se neukázalo jako nevyhnutelné. Diagnostické grafy se o něco málo zlepšily, ale ani u menšŕho modelu nebyly výrazně problematické. Interakce přitom není ani hraničně
116
statisticky významná.
Na základě dat tedy pro čas přeplavání jezerem Ontario ve směru z jihu na sever v měsíci srpnu vybereme model
Timej = (30 + (3± x Sex^ + (32 x Age^ + Si,   i = 1,..., n.
Připomeňme si jeho výsledky.
> summary(model.2)
Call:
lm(formula = Time  ~  Sex + Age,   data = swim)
Residuals:
Min 1Q    Median 3Q Max
-351.86 -167.51    -66.03    173.03 496.01
Coefficients:
Estimate Std.  Error t value Pr(>|t|)
(Intercept)   1106.163 81.241 13.616 <2e-16 *•
SexM               -133.967 65.298 -2.052 0.0463 *
Age                        4.383 2.376 1.845 0.0720 .
Signif. codes
0
0 . 001
0 .01
0 . 05
0 .1
Residual standard error:  215 on 43 degrees of freedom Multiple R-squared:     0.1284,Adjusted R-squared: 0.08788 F-statistic:  3.168 on 2 and 43 DF,    p-value: 0.05208
Na základě modelu můžeme říci, že muži jezero přeplavou v průměru o 134 minut rychleji než ženy stejného věku, a plavcům stejného pohlaví trvá přeplavání jezera s každým rokem věku v průměru o 4 minuty více. Věk a pohlaví ale vysvětlují pouhých 13 % variability mezi časy přeplavání jezera. Podstatnější část variability mezi výkony jednotlivých plavců tedy nesouvisí s jejich věkem ani pohlavím.
Vliv pohlaví jsme prokázali pětiprocentní hladině, zatímco vlik věku jenom na desetiprocentní. I když diagnostické grafy nedávají důvod na pochybnosti o platnosti testů, těmto poměrně hraničním výsledkům nemůžeme připisovat plnou váhu, jelikož jsme si na deskriptivních statistikách pro jména plavců mojli všimnout, že dva plavci přeplavali jezero dvakrát a jeden třikrát. Jejich časy nemůžeme považovat za nezávislé, jak to vyžadují předpoklady lineárního modelu. I když se nejedná o velký počet závislostí, u celkového počtu 46 dat by bylo korektnější použít místo standardních testů pro lineární modely robustnější testy, které platí i v případě závislosti mezi pozorováními. Ty by p—hodnoty pravděpodobně o něco málo posunuly.
Na závěr ještě můžeme výsledky modelu vykreslit.
117
> xx <- seq(min(swim$Age), max(swim$Age), length=501)
> newdat <- expand.grid(Age = xx,
+ Sex = factor (c ("F", "M")))
> newdat <- cbind(newdat,   Time = NA,   Lw = NA,   Up = NA)
> yy.conf <- predict(model.2,  newdata = newdat, + interval = "confidence")
> newdat[,  c(3:5)]   <- yy.conf
>
> ggplot(data = swim,
+ mapping = aes(x = Age,  y = Time,   color = Sex) ) +
+ geom__point () +
+ geom_ribbon(data = newdat,
+ aes(ymin = Lw,  ymax = Up,
+ fill = Sex,   color = NULL),
+ alpha =  .15) +
+ geom_line(data = newdat) +
+ xlab("Age")   + ylab("Time   [min]") +
+ ggtitle("(Expected)   time to swim Lake Ontario",
+ subtitle = "in the SN direction in August") +
+ theme(plot.title = element_text(hjust = 0.5),
+ plot.subtitle = element_text(hjust = 0.5))
Literatura
[1] Ilustrační data: Childhood respiratory disease.   Dostupná z http ://www. statsci . org/data/general/f ev. html. cit. 29. 1. 2021.
[2] Ilustrační data: Successful Lake Ontario swims (Toronto).   Dostupná z http://www. soloswims . com/swims . htm. cit. 29. 1. 2021.
[3] Anděl, J. (2007). Základy matematické statistiky. Matfyzpress, Druhé vydání.
[4] Bubeníček, J. (2019). Shapirův-Wilkův test. Bakalářská práce, Masarykova Univerzita.
[5] Faraway, J. J. (2014). Linear Models with R. Chapman & Hall/CRC, Druhé vydání.
[6] R Core Team (2020). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing.
[7] Wood, S. (2006).   Generalized Additive Models: An Introduction with R.   Chapman & Hall/CRC.
[8] Zvára, K. (2008). Regrese. Matfyzpress.
118