Řešený praktický příklad: závislost koncentrace vitamínu na BMI

Vraťme se k příkladu, kterým jsme tuto výukovou jednotku začínali: modelujeme závislost sérové koncentrace vitaminu D (proměnná vitd) na indexu tělesné hmotnosti (proměnná bmi). K dispozici máme následující datovou tabulku se vzorkem 41 irských žen z datového souboru [vitaminD].

Proměnná vitd představuje výsledkovou proměnnou, proměnná bmi prediktor:

Takto tedy vypadá vektor výsledků (Y)

A takto matice plánu X

Nyní odhadneme parametry tohoto jednoduchého modelu prostřednictvím programu R. Nejprve si ukážeme zdlouhavější postup kopírující výpočty popsané v kapitole Jak definujeme lineární regresní model, poté okomentujeme syntaxi a výsledky volání funkce lm(), která se v programu R používá pro odhad parametrů lineárního regresního modelu. Následující rámečky uvádějí v levém sloupci kód v programu R a jeho slovní popis, v pravém sloupci pak symbolická reprezentace uvedeného postupu.

XX <- t(X) %*% X součin transponované matice X a původní matice X

Beta_hat <- XX_inv %*% ( t(X) %*% Y ) výpočet odhadu regresních koeficientů

Y_hat <- X %*% Beta_hat predikované hodnoty výsledku

c <- matrix(c(0,1),ncol=1) sloupcový vektor T <- abs(t(c) %*% beta_hat) / (sqrt(s2) * sqrt( t(c) %*% XX_inv %*% c)) testové kritérium pro t-test

qt(0.975,41-2) 97,5% kvantil studentova rozdělení s 39 stupni volnosti

2*(1-pt(T,41-2)) p-hodnota

Prakticky je samozřejmě odhad parametrů v programu R výrazně jednodušší – využijeme připravené funkce lm(). Základní výsledky získáme odesláním následujících funkcí:

model <- lm(vitd ~ bmi, data = irlwomen)

summary(model)

Dostáváme následující výsledek:

Nejprve je zopakována formulace regresního modelu ve funkci lm():

Call:

lm(formula = vitd ~ bmi, data = irlwomen)

Následuje základní popisná statistika vektoru reziduí:

Residuals:

   Min     1Q Median     3Q    Max

-25.36 -12.96  -0.41  11.09  52.83

Zde je uveden samotný odhad modelových parametrů, spolu s potřebnými testovými statistikami (postupně jsou uvedeny bodové odhady, směrodatné chyby těchto odhadů, hodnoty t-statistiky a příslušné p-hodnoty pro nulovou hypotézu rovnosti koeficientu 0):

Coefficients:

            Estimate Std. Error t value Pr(>|t|)   

(Intercept) 111.0535    18.4014   6.035 4.63e-07 ***

bmi          -2.3924     0.6902  -3.466   0.0013 **

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Reziduální součet čtverců a počet stupňů volnosti:

Residual standard error: 17.91 on 39 degrees of freedom

Koeficient determinace a jeho adjustovaná varianta (její hodnota se snižuje s rostoucím počtem prediktorů, může tedy být vhodnější pro srovnávání predikční síly modelů):

Multiple R-squared:  0.2355, Adjusted R-squared:  0.2159

A konečně F-statistika pro významnost všech prediktorů zároveň. Všimněte si, že v tomto případě (pouze jediný spojitý prediktor) je p-hodnota totožná s p-hodnotou t-testu významnosti prediktoru bmi.

F-statistic: 12.02 on 1 and 39 DF,  p-value: 0.001299