Rmd source

Analiza współzależności pomiędzy cechami

Pomiędzy zjawiskami występują związki (zależności.) Nauki formułują te związki w postaci praw. Jak takie prawo naukowe powstaje? Typowo w dwu etapach, najpierw za pomocą dedukcji stawia się hipotezę, potem konfrontuje się hipotezę z danymi (podejście hipotetyczno-dedukcyjne). Na tym drugim etapie używa się statystyki (lub matematyki jeżeli prawo ma charakter deterministyczny)

Upraszczając metoda hypodedukcji sprowadza się do dedukcyjnego sformułowania hipotezy, która następnie jest empirycznie falsyfikowana, tj. próbuje się wykazać, że jest ona nieprawdziwa. Konsekwencje: nie można dowieść prawdziwości żadnej hipotezy, można natomiast wykazać, że hipoteza jest fałszywa.

Związki między cechami mogą być: funkcyjne (nauki przyrodnicze) – wartościom jednej zmiennej odpowiada tylko jedna wartość drugiej zmiennej lub stochastyczne – wartościom jednej zmiennej odpowiadają z pewnym przybliżeniem wartości innej zmiennej.

Problem: czy istnieje związek (zależność) pomiędzy cechami? Dla uproszczenie pomiędzy dwoma cechami, np. czy istneje związek pomiędzy paleniem a chorobą nowotworową, wiekem a prawdopodobieństem zgonu z powodu COVID19 itd

Jaki jest charakter zależności? Jaka jest siła zależności?

Korelacyjny wykres rozrzutu (korelogram, wykres XY w Excelu, scatter plot)

Obie cechy muszą być mierzalne. W układzie kartezjańskim każdej obserwacji odpowiada kropka o współrzędnych XY.

O występowaniu związku świadczy układanie się kropek według jakiegoś kształtu (krzywej). O braku związku świadczy chmura punktów niepodobna do żadnej krzywej.

Punkty układające się według prostej świadczą o zależności liniowej (wyjątek: linia pozioma lub pionowa) Punkty układające się według krzywej świadczą o zależności nieliniowej.

Indianie Pima (Arizona/Meksyk)

The Pima Indians of Arizona have the highest reported prevalences of obesity and non-insulin-dependent diabetes mellitus (NIDDM). In parallel with abrupt changes in lifestyle, these prevalences in Arizona Pimas have increased to epidemic proportions during the past decades (https://care.diabetesjournals.org/content/17/9/1067)

Przykład BMI a poziom glukozy dla Indianek z plemienia Pima (USA/Meksyk)

Testy a zakażeni wg województw w PL

Wg teorii lansowanej przez pewną antyrządową telewizję (wszyscy wiemy o jaką TV chodzi:-)) Liczba zarażonych w PL jest zaniżona bo wykonuje się ciągle zbyt mało testów. Można zweryfikować (w naiwny sposób) to stwierdzenie wykreślając liczę testów oraz liczbę zarażonych i/lub zmarłych wg województw

oraz

Wniosek: twierdzenie antyrządowej TV wydaje się słabo uzasadnione zwłaszcza w odniesieniu do zależności testy vs zmarli. Dla pewności to sama zależność dla danych z wybranych krajów świata (dane pobrane ze strony ECDC)

Hmm…

Rodzaje zależności (brak, liniowa, nieliniowa)

Spożycie mięsa a GDP (Świat 2013 rok)

Obserwacje nietypowe (odstające)

Tabele korelacyjne

Łączny rozkład dwóch (dwudzielna) lub większej (wielodzielna) liczby zmiennych można przedstawić w tabeli (zwanej także korelacyjną). W języku angielskim tego typu tabele noszą nazwę two-way tables albo contingency tables.

Przykładowo plik smoker.csv zawiera dane dotyczące palenia (pali, palił-nie-pali, nigdy-nie-palił) oraz statusu społeczno-ekonomicznego (wysoki, średni, niski). Liczebności poszczególnych kategorii są następujące:

smokerData <- read.csv(file='smoker.csv',sep=';',header=T)
summary(smokerData)
##     Smoke               SES           
##  Length:356         Length:356        
##  Class :character   Class :character  
##  Mode  :character   Mode  :character
nrow(smokerData)
## [1] 356

Tabele korelacyjna dla danych z pliku smoker.csv wygląda następująco:

smoke <- table(smokerData)
smoke
##          SES
## Smoke     High Low Middle
##   current   51  43     22
##   former    92  28     21
##   never     68  22      9
## rozkłady brzegowe
margin.table(smoke,1)
## Smoke
## current  former   never 
##     116     141      99
margin.table(smoke,2)
## SES
##   High    Low Middle 
##    211     93     52

Albo wszystko razem ładnie wydrukowane:

smoke.table <- addmargins(smoke)
kable(smoke.table)
High Low Middle Sum
current 51 43 22 116
former 92 28 21 141
never 68 22 9 99
Sum 211 93 52 356

Tabela korelacyjna wyrażona w udziałach (%):

total <- sum(smoke)
smoke.p <- smoke/total *100
smoke.p
##          SES
## Smoke          High       Low    Middle
##   current 14.325843 12.078652  6.179775
##   former  25.842697  7.865169  5.898876
##   never   19.101124  6.179775  2.528090
## rozkłady brzegowe
margin.table(smoke,1)/total * 100
## Smoke
##  current   former    never 
## 32.58427 39.60674 27.80899
margin.table(smoke,2)/total * 100
## SES
##     High      Low   Middle 
## 59.26966 26.12360 14.60674

Albo wszystko razem ładnie wydrukowane:

smoke.table.p <- addmargins(smoke.p)
kable(smoke.table.p)
High Low Middle Sum
current 14.32584 12.078652 6.179775 32.58427
former 25.84270 7.865169 5.898876 39.60674
never 19.10112 6.179775 2.528090 27.80899
Sum 59.26966 26.123596 14.606742 100.00000

Rozkłady oczekiwane i test chi-kwadrat

Można udowodnić że jeżeli cechy \(X\) i \(Y\) sa niezależne to liczebność \(n_{ij}\) w \(ij\) tej kratce tabeli powinna się równać iloczynowi odpowiednich wartości rozkładów brzegowych podzielonemu przez \(n\), tj. \(n_{ij} = n_{i} \cdot n_{j} / n\) . Przykładowo w tabeli palenie-SES wszystkich jest 356 a palących 116, tj stanowią oni 32.58% całej populacji. Jeżeli SES nie ma wpływu na palenie to w każdej kolumnie ten udział powinien być identyczny czyli wynosić \(211 \times 116/356\), \(93 \times 116/356\) oraz \(52 \times 116/356\).

Tak obliczone liczebności nazywają się oczekiwane (lub teoretyczne)

chi_smoke <- chisq.test(smoke)
##chi_smoke
smoke.expected <- chi_smoke$expected
smoke.expected
##          SES
## Smoke         High      Low   Middle
##   current 68.75281 30.30337 16.94382
##   former  83.57022 36.83427 20.59551
##   never   58.67697 25.86236 14.46067
smoke - smoke.expected
##          SES
## Smoke            High         Low      Middle
##   current -17.7528090  12.6966292   5.0561798
##   former    8.4297753  -8.8342697   0.4044944
##   never     9.3230337  -3.8623596  -5.4606742
## rozkłady brzegowe
margin.table(smoke, margin = 1)
## Smoke
## current  former   never 
##     116     141      99
margin.table(smoke, margin = 2)
## SES
##   High    Low Middle 
##    211     93     52
# prop.table = Express Table Entries as Fraction of Marginal Table
prop.table(smoke, margin = 1)
##          SES
## Smoke           High        Low     Middle
##   current 0.43965517 0.37068966 0.18965517
##   former  0.65248227 0.19858156 0.14893617
##   never   0.68686869 0.22222222 0.09090909
prop.table(smoke, margin = 2)
##          SES
## Smoke          High       Low    Middle
##   current 0.2417062 0.4623656 0.4230769
##   former  0.4360190 0.3010753 0.4038462
##   never   0.3222749 0.2365591 0.1730769

Albo wszystko razem ładnie wydrukowane:

smoke.table.e <- addmargins(smoke.expected)
kable(smoke.table.e)
High Low Middle Sum
current 68.75281 30.30337 16.94382 116
former 83.57022 36.83427 20.59551 141
never 58.67697 25.86236 14.46067 99
Sum 211.00000 93.00000 52.00000 356

Jak widać rozkłady brzegowe dla wartości teoretycznych są identyczne z rozkładami dla wartości empirycznych. Nie może być inaczej…

Im większa jest różnica pomiędzy wartościami zaobserwowanymi a teoretycznymi tym bardziej wątpimy ze faktycznie cechy X oraz Y są nieskorelowane. Jest to istota testu chi-kwadrat (albo \(\chi\)-kwadrat jeżeli chcemy żeby wyglądało mądrzej, lub \(\chi^2\) jeżeli ma być super-mądrze).

Można bowiem obliczyć prawdopodobieństwo wystąpienia różnicy (definiowanej jako suma kwadratów różnic) równej zaobserwowanej lub większej. W naszym przykładzie prawdopodobieństwo to wynosi 0.0009808. Czyli jest bardzo małe (przykładowo 10-krotne wyrzucenie orła w 10 rzutach monetą zdarza się z prawdopodobieństwem 0.0009765625).

W tej sytuacji uznajemy że założenie o braku korelacji pomiędzy paleniem a SES za niepoprawne. Mówiąc ściśle arbitralnie przyjmujemy pewien poziom błędu (zwany poziomem istotności) zwykle na poziomie 0,01 lub 0,05 i jeżeli wartość prawdopodobieństwa zdarzenie, które przeczy postawionej hipotezie (tutaj że palenie i SES sa nieskorelowane) jest od tego prawdopodobieństwa mniejsza to znaczy że hipoteza jest nieprawdziwa. Czyli istnieje zależność korelacyjna pomiędzy SES i paleniem.

Tabele dwudzielne dla danych ilościowych

Plik rwc2019.csv zawiera m.in. dane dotyczące wagi i wzrostu zawodników którzy brali udział w Pucharze Świata w Rugby 2019 roku

rwc <- read.csv(file='rwc2019.csv',sep=';',header=T)
## https://stackoverflow.com/questions/48478054/how-to-create-a-crosstab-table-using-two-tables-in-r
## W arkuszu PIVOT TABLES vel tabele przestawne
rwc_n <- nrow(rwc)
rwc_h <- nlevels(as.factor(rwc$height))
rwc_w <- nlevels(as.factor(rwc$weight))

W pliku jest 620 zawodników, którzy mają 40 różnych wartości wzrostu oraz 61 różnych wartości wagi. Bezpośrednie utworzenie tabeli dwudzielnej skutkowałoby tabelą o liczbie 2440 komórek, z których większość zawierałaby zero. Sensowne przedstawienie takich danych wymaga wcześniejszego ich pogrupowania

rwc.t <- rwc %>% mutate(Weight_kg = cut(rwc$weight, breaks = c(60, 80,100,120,140,160)),
  Height_cm = cut(rwc$height, breaks = c(150, 160, 170, 180,190, 200, 210)))  %>%
  select (Weight_kg, Height_cm)

rwc.table <- table(rwc.t)
rwc.table.x <- addmargins(rwc.table)
kable(rwc.table.x)
(150,160] (160,170] (170,180] (180,190] (190,200] (200,210] Sum
(60,80] 1 5 27 2 0 0 35
(80,100] 0 0 79 140 19 0 238
(100,120] 0 1 30 162 94 10 297
(120,140] 0 0 4 27 14 4 49
(140,160] 0 0 0 1 0 0 1
Sum 1 6 140 332 127 14 620

Wykres mozaikowy (mosaic plot)

W przypadku danych niemierzalnych do wizualizacji stosuje się wykres mozaikowy.

Wykres mozaikowy wygląda zaś następująco (kolorem wyróżniono SES):

Wysokości słupków wizualizują proporcje liczebności rozkładów warunkowych. Szerokości słupków pokazują z kolei rozkład brzegowy. Czyli w tym przykładzie brzegowy jest SES (high/low/middle równe odpowiednio 211, 93, 52)
a warunkowe dla zmiennej Palenie. Przy założeniu że zmienne są niezależne wysokości słupków powinny być identyczne. Tymczasem widać, że nie są.

Dla utrwalenie wykres mozaikowy wagi/wzrostu zawodników na turnieju o Puchar Świata w Ruby

Wykres nieczytelny `w okolicach’ skrajnych kolumn i wierszy, które są mało liczne. Możemy go spróbować poprawić inaczej grupując dane:

<180 180-85 185-90 190-95 >195 Sum
<80 33 2 0 0 0 35
80-100 79 83 57 15 4 238
100-120 31 71 91 49 55 297
>120 4 14 14 8 10 50
Sum 147 170 162 72 69 620

Że istnieje zależność waga wzg. wzrostu nie ulega wątpliwości?

Linus Pauling i witamina C

Linus Pauling (laureat nagrody Nobla/witamina C). Podawanie witaminy C a przeziębienie/brak przeziębienia:

A study of the effect of 1 g of ascorbic acid per day was reported […] The study was carried out in a ski resort with 279 skiers during two periods of 5-7 days.

Jak ja rozumiem ten opis to eksperyment polegał na tym, że przez 5–7 dni podawał witaminę C oraz placebo i obserwował zachorowania na przeziębienie.

nocold cold razem
C 122 17 139
P 109 31 140
Sum 231 48 279

Dla przypomnienia (wykres mozaikowy):

Dwa rozkłady brzegowe: Brał placebo/Brał witaminę C (140 vs 139) Zachorował/nie zachorował (48/231)

Cztery rozkłady warunkowe: Brał placebo oraz zachorował/nie zachorował 31 vs 109 lub \(31/(31+109)=0,2214\) oraz \(109/(109+31.0) = 0.7785\). Interpretacja: 22,14% tych którzy brali placebo zachorowało, a 77,85% nie zachorowało.

Brał witaminę C oraz zachorował/nie zachorował 17 vs 122 lub \(17/(17+122.0) = 0.1223\) oraz \(122/(17+122.0) = 0.8776\) Interpretacja: 12,23% tych którzy brali witaminę C zachorowało, a 87,76% nie zachorowało.

Zachorował oraz brał placebo/witaminę C 31 vs 17 Nie zachorował oraz brał placebo/witaminę C 109 vs 122

Statystyka medyczna: ryzyko względne i iloraz szans

Ryzyko to udział (iloraz) liczby sukcesów do liczby prób (zdarzeń pozytywnych/wyróżnionych do wszystkich). Zwykle podawany w procentach. Warto zauważyć że jest to empiryczny odpowiednik prawdopodobieństwa. Przykładowo jeżeli w grupie 139 narciarzy którym podano witaminę C zachorowało 17, to ryzyko zachorowania wyniosło 17/139 = 12,2%. Wiemy też że w grupie tych co nie brali witaminy C ryzyko to wyniosło 22,14%.

Prostymi miarami oceny siły zależności mogą być: różnica ryzyk (risk difference) ryzyko względne (relative risk) oraz iloraz szans (odds ratio).

Jeżeli \(r_e\) oznacza ryzyko w grupie eksperymentalnej (test group; grupa narażona/exposed group) a \(r_k\) w grupie kontrolnej (control group; grupa nienarażona/unexposed), to różnica ryzyk to po prostu \(r_e - r_k\). W przykładzie będzie to -9,94% Ta miara aczolwiek prosta jest rzadko stosowana.

Znacznie częściej używa sę ryzyka względnego definiowanego jako \(rr = r_e/r_k\). W przykładzie będzie to 12,2/22,14 = 0,55. Podanie witaminy C zmniejsza ryzyko o prawie połowę. Oczywiste jest że \(rr < 1\) oznacza zmniejszenie ryzyka; \(rr > 1\) zwiększenia a \(rr = 1\) oznacza brak zależności.

Zamiast ryzyka (czyli ilorazu sukcesów do prób) można używać pojęcia szansa/szansy (odds) definiowanych jako iloraz sukcesów do porażek.

Przykład #1: jeżeli w dwukrotnym rzucie monetą otrzymano orła i reszkę to ryzyko otrzymania orła wynosi 1/2 = 0,5 a szansa otrzymania orła wynosi 1. Przykład #2: ryzyko zachorowania w grupie narciarzy którym podano witaminę C wynosi 12,2 (jak wiemy); szansa że narciarz w tej grupie zachoruje wynosi 17/122 = 13,9%. (A w tej drugiej grupie 28,44%)

Jak widać dla dużych prawdopodobieństw (rzut monetą) szansa różni się znacznie od prawdopodobieństwa, ale dla małych ryzyk obie miary mają zbliżoną wartość.

Jeżeli \(o_e\) oznacza szanse w grupie eksperymentalnej a \(o_k\) w grupie kontrolnej, to iloraz szans (odds ratio), jest definiowany jako stosunek \(or = o_e/o_t\). Zatem iloraz szans dla narciarzy wyniesie 13,9/28,44 = 0,48. Podanie witaminy C zmniejsza szansę na zachorowanie o ponad połowę. Albo 1/0,48 = 2,04, narciarz który nie brał witaminy C ma ponad dwukrotnie większą szansę na zachorowanie.

Właściwości \(or\): jeżeli równe 1 to sukces/porażka równie prawdopodobne; jeżeli większe od 1 to sukces bardziej prawdopodobny; jeżeli jest mniejsze od 1 to porażka jest bardziej prawdopodobna.

Dane w badaniach wykorzystujących ryzyko/szanse mają często postać tabeli dwudzielnej o wymiarach 2x2, którą można przestawić następująco (a, b, c i d to liczebności):

                      | sukces |  porażka
----------------------+--------+---------
grupa kontrolna       | a      |  b
grupa eksperymentalna | c      |  d
----------------------+--------+---------

Dla danych w tej postaci: \(rr = c(a+b)/a(c+d)\) oraz \(or= (ad)/ (bc)\)

Przedziały ufności dla błędu względnego i ilorazu szans

Dla dużej próby można założyć, że ln(RR) ma rozkład normalny o odchyleniu standardowym równym: \(se_{rr} (ln(rr)) = \sqrt{1/a + 1/c - 1/(a+b) - 1/(c+d)}\), stąd 95% przedział określony jest przez: \(rr - 1,96 se_{rr} < ln(RR) < rr + 1,96 se_{rr}\), albo po przekształceniu: \(\exp(ln(rr) \pm 1,96 se_{rr})\) (dla przypomnienia \(\exp(x)\) to \(e^x\))

Dla dużej próby można założyć, że ln(OR) ma rozkład normalny o odchyleniu standardowym równym: \(se_{or} (ln(or)) = \sqrt{1/a + 1/b + 1/c + 1/d}\), stąd 95% przedział określony jest przez: \(or - 1,96 se_{or} < ln(OR) < or + 1,96 se_{or}\), albo po przekształceniu: \(\exp(ln(or) \pm 1,96 se_{or})\).

Jeżeli przedział ufności zawiera 1 to świadczy to o braku zależności (na 95% poziomie ufności).

Przykład: Przedział ufności dla ilorazu szans ln(or) = ln(0,48) = -0.7339691 \(se{or} (ln(or)) = 0.3293214\) stąd \(-0,7339691 \pm 1,96 \cdot 0.3293214\). Końce przedziałów: [-1.379439044; -0.0884991560]; ostatecznie [0,2517; 0,9153]

Przedział nie zawiera 1; zatem branie witaminy C zmniejsza szanse na zachorowanie; albo zwiększa na niezachorowanie od \(1/25 = 4\) do \(1/0,9 = 1,1\). Żeby to zabrzmiało ładnie i po polsku. Zwiększa na niezachorowanie od 300% do 10%.

Test chi-kwadrat (Pearsona)

Dla przypomnienia: idea jest prosta – porównujemy liczebności tabeli korelacyjnej z liczebnościami teoretycznymi, to znaczy takimi, które wynikają z przyjęcia założenia że cechy \(X\) i \(Y\) są niezależne. Innymi słowy oceną wielkości związku jest różnica pomiędzy liczebnościami teoretycznymi a liczebnościami empirycznym.

Narciarze leczeni witaminą C

nocold cold razem
C 122 17 139
P 109 31 140
Sum 231 48 279

Gdyby \(XY\) były niezależne to tabela powinna wyglądać następująco

N Y Sum
C 115.086 23.91398 139
P 115.914 24.08602 140
Sum 231.000 48.00000 279

Policzmy różnice (pierwszy wiersz) 122 - 115.0860215 + 17 - 23.9139785. W sumie zero i dla drugiego wiersza oczywiście też jest zero. Ale różnice podniesione do kwadratu już nie będą się zerowały: ((122 - 115.0860215)² + (17 - 23.9139785)²) + ((109 - 115.9139785)² + (31 - 24.0860215)²) = 191.2123945

Im większa różnica, tym nasze przekonanie, że istnieje związek pomiędzy cechami jest większe (co oczywiste). Pytanie kiedy jest dużo? I na to pytanie obiektywnie i zawsze tak samo odpowiada test \(\chi^2\). I tym się różni metoda statystystyczna od medycyny ludowej lub prognozy przysłowiowego górala:

chi_skiers$statistic
## X-squared 
##  4.140679
chi_skiers$p.value
## [1] 0.04186438

Z czego by wynikało że na poziomie \(\alpha=0,05\) hipotezę o niezależności cech należy odrzucić ale na poziomie \(\alpha = 0,01\) już nie. Dlaczego?

Pomiar siły zależności: współczynnik korelacji liniowej Pearsona

Kowariancja to suma iloczynów odchyleń wartości zmiennych \(XY\) od ich wartości średnich:

\[cov (xy) = \frac{1}{n} \sum_{i=1}^N (x - \bar x) (y - \bar y)\]

Kowariacja zależy od rozproszenia (im większe tym większa), ma też dziwną jednostkę (jednostkaX · jednostkaY) oraz zależy od wybranych skal (tony vs gramy na przykład)

https://tinystats.github.io/teacups-giraffes-and-statistics/05_correlation.html

Dlatego do pomiaru związku pomiędzy cechami nie używa się kowariancji ale współczynnika korelacji (liniowej, Pearson linear correlation coefficient):

\[r_xy = cov(xy) / (S_x \cdot S_y)\]

Współczynnik jest miarą niemianowaną, o wartościach ze zbioru \([-1;1]\); Skrajne wartości \(\pm 1\) świadczą o związku funkcyjnym (wszystkie punkty układają się na linii prostej); wartość zero świadczy o braku związku (linia pozioma/pionowa)

Interpretacja opisowa: wartości powyżej 0,9 świadczą o silnej zależności

Pomiar siły zależności: współczynnik korelacji rang

Współczynnik korelacji rang (Spearmana vel Spearman’s Rank-Order Correlation) może być stosowany w przypadku gdy cechy są mierzone w skali porządkowej (lub lepszej)

Obliczenie współczynnika Spearmana dla \(N\) obserwacji na zmiennych \(XY\) polega na zamianie wartości \(XY\) na rangi (numery porządkowe od \(1...N\)). Następnie stosowana jest formuła Pearsona, tj. (\(\tau_x\) oraz \(\tau_y\) oznaczają rangi):

\[\rho_{xy} = cov(\tau_x, \tau_y) / ( S_{\tau_x} \cdot S_{\tau_y} )\]

Współczynnik \(\rho_{xy}\) to także miara niemianowana, o wartościach ze zbioru [-1;1];

W podręcznikach można też spotkać formułę alternatywną:

\[r_s = 1 - \frac{6 \cdot \sum_\limits{i=1}^N d_i^2}{ N(N^2 -1)}\]

gdzie \(d_i\) oznacza różnicę między rangami dla cechy \(X\) oraz \(Y\), tj. \(d_i = \tau_{x_i} - \tau_{y_i}\). Wymagany jest aby \(d_i\) były różne od siebie. Można też spotkać (w podręcznikach) przykłady użycia formuły alternatywnej dla rang identycznych, wówczas (dla \(k\) identycznych rang, zamienia się je na wartości średnie) \(d_i = \frac{1}{k}\sum_\limits{i=1}^k d_i^k\). Raczej nie zalecane…

Przykład: testy z zakażeni wg województw w PL

Współczynnik Pearsona i Spearmana dla zależności między liczbą testów a liczbą zarażonych wirusem COVID19 (PL/województwa):

## [1] 0.1962708
## [1] 0.2852941

Współczynnik Pearsona i Spearmana dla zależności między liczbą testów a liczbą zmarłych z powodu COVID19:

## [1] 0.01310447
## [1] -0.005882353

Przykład: spożycie mięsa

Współczynnik Pearsona i Spearmana dla zależności między spożyciem mięsa w 1980 a spożyciem mięsa w 2013 roku (zmienna objaśniana):

## [1] 0.6823158
## [1] 0.6845429

Regresja liniowa

Cechy muszą być mierzone za pomocą skali co najmniej interwałowej.

Regresja liniowa zakłada że istnieje związek przyczyna-skutek i ten związek można opisać linią prostą (stąd liniowa). Skutek jest jeden i nazywa się go zmienną zależną a przyczyn może być wiele i noszą nazwę zmiennych niezależnych. W przypadku gdy związek dotyczy dwóch zmiennych mówi się o regresji prostej. Regresja wieloraka to zależność pomiędzy jedną zmienną objaśnianą a wieloma zmiennymi objaśniającymi

Liniowa funkcja regresji ma zatem postać: \(Y = a \cdot X + b\)

Współczynnik \(a\) ma jasną interpretację: jeżeli wartość zmiennej \(X\) rośnie o jednostkę to wartość zmiennej \(Y\) zmienia się przeciętnie o \(a\). Wyraz wolny zwykle nie ma sensu (formalnie jest to wartość zmiennej \(Y\) dla \(X=0\))

Oznaczmy przez \(y_i\) wartości obserwowane (empiryczne) a przez \(\hat y_i\) wartości teoretyczne (leżące na prostej linii regresji). Wartości \(a\) oraz \(b\) wyznacza się minimalizując: \(\sum_{i=1}^N (\hat y_i - y_i)^2\) (suma kwadratów odchyleń wartości empirycznych od wartości teoretycznych ma być minimalna)

Ocena dopasowania

Średni błąd szacunku (MSE mean squared error). Oznaczając resztę jako: \(e_i = y_i - \hat y_i\), definiujemy średni błąd szacunku jako: \(S_e = \sqrt {\sum (e_i^2 / n -k)}\). (Przy okazji \(S_e^2\) nazywamy wariancją resztową)

Współczynnik zbieżności i determinacji

Suma kwadratów reszt (albo odchyleń wartości teoretycznych od wartości empirycznych, albo suma kwadratów błędów vel resztowa suma kwadratów): \(RSK = \sum (y_i - \hat y_i)^2\). Suma kwadratów odchyleń wartości empirycznych od średniej (ogólna suma kwadratów): \(OSK = \sum (y_i - \bar y)^2\) Suma kwadratów odchyleń wartości teoretycznych od średniej (wyjaśniona suma kwadratów): \(WSK = \sum (\hat y_i - \bar y)^2\)

Można wykazać, że \(OSK = WSK + RSK\). Współczynnik zbieżności to \(R^2 = WSK/OSK\). Współczynnik determinacji to \(\Phi^2 = RSK/OSK\).

Współczynniki przyjmują wartość z przedziału \([0,1]\) lub \([0, 100]\)% Interpretacja współczynników: procent zmienność wyjaśniane/nie wyjaśniane przez linię regresji.

Zerowa wartość parametru \(a\)

Jeżeli: \(Y= 0 \cdot X + b\), to \(Y = b\) czyli nie ma zależności pomiędzy \(XY\). Wartości \(a\) bliskie zero wskazują na słabą zależność pomiędzy cechami. Przyjmijmy

\[\hat Y = a \cdot \hat X + b + e\]

w tym zapisie obserwowana linia regresji jest równa linii teoretycznej \(Y = a \cdot X + b\) plus błąd losowy. (przy czym losowość błędu przejawia się m.in w tym, że nie ma on charakteru systematycznego, tj. jest określona jako zero plus/minus cośtam)

Można teraz zadać pytanie jeżeli \(a=0\) to jakie jest prawdopodobieństwo, że współczynnik \(\hat a\) oszacowany na podstawie \(n\) obserwacji mierzonych z błędem \(e\) będzie (co do wartości bezwzględnej) większy niż \(a_e\). Albo inaczej: otrzymaliśmy \(a_e\), jakie jest prawdopodobieństwo otrzymania takiej wartości (lub mniejszej co do wartości bezwzględnej) przy założeniu że istotnie \(a=0\).

Jeżeli takie prawdopodobieństwo jest duże to zakładamy że być może \(a=0\) a jeżeli małe to że \(a \not0\). Duże/małe przyjmujemy arbitralnie, zwykle jest to \(0,1\), \(0,05\) lub \(0,01\).

Pakiety statystyczne szacują ww. prawdopodobieństwo, zwykle oznaczając je prob. Reasumując jeżeli przyjeliśmy \(0,05\) jako wartość na tyle małą że niemożliwą w realizacji oraz otrzymaliśmy \(prob=0,006\) to mówimy że współczynnik \(a\) jest istotnie różny od zera.

Testowanie istotności współczynnika regresji jest ważnym kryterium oceny jakości dopasowania. Regresja z nieistotnym współczynnikiem nie może być podstawą do interpretowania zależności pomiędzy XY.

BMI a glukoza dla Indian Pima

Szacujemy liniową funkcję regresji dla zmiennych glucose (objaśniana) oraz bmi (objaśniająca):

## (Intercept)         bmi 
##   92.212896    0.896509
## [1] 0.04887242
##              Estimate Std. Error   t value
## (Intercept) 92.212896  4.7082932 19.585207
## bmi          0.896509  0.1428985  6.273746
##                                                                                  Pr(>|t|)
## (Intercept) 0.000000000000000000000000000000000000000000000000000000000000000000001478974
## bmi         0.000000000589141252708448832510577241613298690725564199510699836537241935730

Testy a zmarli wg województw (PL)

##      (Intercept)            testy 
## 265.747361489396   0.000003530217
## [1] 0.0001717273
##                     Estimate     Std. Error    t value       Pr(>|t|)
## (Intercept) 265.747361489396 32.57264599659 8.15860528 0.000001088874
## testy         0.000003530217  0.00007199137 0.04903666 0.961582788709

Testy a zarażeni (dane ECDC)

Im więcej testów tym więcej zakażonych:

##   (Intercept)    tests_done 
## 311.421050064   0.003171251
## [1] 0.02352672
##                  Estimate   Std. Error   t value    Pr(>|t|)
## (Intercept) 311.421050064 96.543664682 3.2257016 0.003281042
## tests_done    0.003171251  0.003931861 0.8065521 0.426972229

GDP a spożycie mięsa

Im większe GDP tym więcej ludzie jedzą mięsa:

##  (Intercept)      gdp2013 
## 34.084768118  0.001107544
## [1] 0.4087518
##                 Estimate    Std. Error  t value
## (Intercept) 34.084768118 2.23247986990 15.26767
## gdp2013      0.001107544 0.00009956108 11.12427
##                                               Pr(>|t|)
## (Intercept) 0.0000000000000000000000000000000003049703
## gdp2013     0.0000000000000000000003460772431846800413

Prognozowanie: ile wyniesie liczba zmarłych gdyby wykonywano \(p\) testów. Ile wyniesie poziom glukozy dla BMI równego \(x\).

Regresja nieliniowa

Jeżeli zależność ma charakter wykładniczy można założyć, że

\[Y = a \cdot X^b\] co można sprowadzić do zależności liniowej przez obstronne zlogarytmowanie:

\[ln(Y) = ln(a) + b \cdot ln(Y)\] Interpretacja: jeżeli wartość wzrośnie o 1% to wartość \(Y\) wzrośnie/spadnie o \(b\)%. Czyli stały nie jest wzrost, ale względne tempo wzrostu.

Regresja logistyczna

Jeżeli \(Y\) jest zmienną dwuwartościową, to jest taką, która przymuje tylko dwie wartości (chory/zdrowy) to metoda regresji nie może być zastosowana. Przykładowo jeżeli zakodujemy te wartości jako 0 i 1 odpowiednio, to zastosowanie regresji doprowadzi do obliczenia (teoretycznych) wartości \(Y\) różnych od \(0\) i \(1\). Taki wynik nie ma sensownej interpretacji…

Ale zamiast szacować regresję \(Y\) względem (\(X\)/\(X\)-ów) można szacować regresję względem ryzyka dla \(Y\) (czyli prawdopodobieństwa że \(Y\) przyjmnie wartość 1). Konkretnie jeżeli to ryzyko wynosi \(r\), to model można zapisać jako:

\[ln(r/(1-r)) = a_0 + a_1 \cdot x_1 + \ldots + a_k \cdot x_k\]$

Jak łatwo zauważyć \(r/(1-r)\) to nic innego jak szansa (odds). Czemu tak? Ano głównie temu, że transformacja \(ln(o)\) ma prostą interpretację: jeżeli zmienna \(Y\) jest dwuwartościowa to \(\exp(a)\) przy zmiennej \(X\) jest ilorazem szans. Jeżeli \(X\) wzrośnie \(o\) jednostkę, to logarytm ilorazu szans wzrośnie o \(log(O_R)\) lub iloraz szans wzrośnie (lub spadnie) i \(exp^{log(O_R)} \cdot 100 - 100\) procent. Jeżeli \(X\) też jest dwuwartościowe to interpretacja jest jeszcze bardziej intuicyjna: jest to iloraz szans dla przypadku gdy \(X=1\).

Przykład #1: narciarze i witamina C

## 
## Call:
## glm(formula = factor(cold) ~ factor(treatment), family = binomial, 
##     data = skiers)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7075  -0.7075  -0.5108  -0.5108   2.0500  
## 
## Coefficients:
##                    Estimate Std. Error z value           Pr(>|z|)    
## (Intercept)         -1.9708     0.2589  -7.613 0.0000000000000268 ***
## factor(treatment)P   0.7134     0.3293   2.166             0.0303 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 256.18  on 278  degrees of freedom
## Residual deviance: 251.31  on 277  degrees of freedom
## AIC: 255.31
## 
## Number of Fisher Scoring iterations: 4
## <NA> 
##   NA

Przykład #2: kobiety z plemienia Pima

##   pregnancies        glucose        bloodpress          skin      
##  Min.   : 0.000   Min.   :  0.0   Min.   :  0.00   Min.   : 0.00  
##  1st Qu.: 1.000   1st Qu.: 99.0   1st Qu.: 62.00   1st Qu.: 0.00  
##  Median : 3.000   Median :117.0   Median : 72.00   Median :23.00  
##  Mean   : 3.845   Mean   :120.9   Mean   : 69.11   Mean   :20.54  
##  3rd Qu.: 6.000   3rd Qu.:140.2   3rd Qu.: 80.00   3rd Qu.:32.00  
##  Max.   :17.000   Max.   :199.0   Max.   :122.00   Max.   :99.00  
##     insulin           bmi        DiabetesPedigreeFunction      Age       
##  Min.   :  0.0   Min.   : 0.00   Min.   :0.0780           Min.   :21.00  
##  1st Qu.:  0.0   1st Qu.:27.30   1st Qu.:0.2437           1st Qu.:24.00  
##  Median : 30.5   Median :32.00   Median :0.3725           Median :29.00  
##  Mean   : 79.8   Mean   :31.99   Mean   :0.4719           Mean   :33.24  
##  3rd Qu.:127.2   3rd Qu.:36.60   3rd Qu.:0.6262           3rd Qu.:41.00  
##  Max.   :846.0   Max.   :67.10   Max.   :2.4200           Max.   :81.00  
##      Class      
##  Min.   :0.000  
##  1st Qu.:0.000  
##  Median :0.000  
##  Mean   :0.349  
##  3rd Qu.:1.000  
##  Max.   :1.000

Traktujemy wartości zerowe zmiennych glucose oraz bloodpress jako błędy i usuwamy je z próby:

p0 <- pima %>% filter(glucose > 0 & bloodpress > 0)
summary(pima)
##   pregnancies        glucose        bloodpress          skin      
##  Min.   : 0.000   Min.   :  0.0   Min.   :  0.00   Min.   : 0.00  
##  1st Qu.: 1.000   1st Qu.: 99.0   1st Qu.: 62.00   1st Qu.: 0.00  
##  Median : 3.000   Median :117.0   Median : 72.00   Median :23.00  
##  Mean   : 3.845   Mean   :120.9   Mean   : 69.11   Mean   :20.54  
##  3rd Qu.: 6.000   3rd Qu.:140.2   3rd Qu.: 80.00   3rd Qu.:32.00  
##  Max.   :17.000   Max.   :199.0   Max.   :122.00   Max.   :99.00  
##     insulin           bmi        DiabetesPedigreeFunction      Age       
##  Min.   :  0.0   Min.   : 0.00   Min.   :0.0780           Min.   :21.00  
##  1st Qu.:  0.0   1st Qu.:27.30   1st Qu.:0.2437           1st Qu.:24.00  
##  Median : 30.5   Median :32.00   Median :0.3725           Median :29.00  
##  Mean   : 79.8   Mean   :31.99   Mean   :0.4719           Mean   :33.24  
##  3rd Qu.:127.2   3rd Qu.:36.60   3rd Qu.:0.6262           3rd Qu.:41.00  
##  Max.   :846.0   Max.   :67.10   Max.   :2.4200           Max.   :81.00  
##      Class      
##  Min.   :0.000  
##  1st Qu.:0.000  
##  Median :0.000  
##  Mean   :0.349  
##  3rd Qu.:1.000  
##  Max.   :1.000

Przekodowujemy zmnienną BMI. Wg definicji WHO otyłości, osoba otyła ma BMI większe od 30:

p0 <- p0 %>% mutate (bmi = cut (p0$bmi, 
        breaks=c(0, 30, Inf), labels=c("no", "o")))
summary(p0)
##   pregnancies        glucose        bloodpress          skin      
##  Min.   : 0.000   Min.   : 44.0   Min.   : 24.00   Min.   : 0.00  
##  1st Qu.: 1.000   1st Qu.:100.0   1st Qu.: 64.00   1st Qu.: 0.00  
##  Median : 3.000   Median :117.0   Median : 72.00   Median :24.00  
##  Mean   : 3.863   Mean   :121.9   Mean   : 72.44   Mean   :21.39  
##  3rd Qu.: 6.000   3rd Qu.:141.2   3rd Qu.: 80.00   3rd Qu.:33.00  
##  Max.   :17.000   Max.   :199.0   Max.   :122.00   Max.   :99.00  
##     insulin         bmi      DiabetesPedigreeFunction      Age       
##  Min.   :  0.00   no  :275   Min.   :0.0780           Min.   :21.00  
##  1st Qu.:  0.00   o   :449   1st Qu.:0.2450           1st Qu.:24.00  
##  Median : 48.00   NA's:  4   Median :0.3800           Median :29.00  
##  Mean   : 84.15              Mean   :0.4765           Mean   :33.39  
##  3rd Qu.:130.00              3rd Qu.:0.6295           3rd Qu.:41.00  
##  Max.   :846.00              Max.   :2.4200           Max.   :81.00  
##      Class       
##  Min.   :0.0000  
##  1st Qu.:0.0000  
##  Median :0.0000  
##  Mean   :0.3434  
##  3rd Qu.:1.0000  
##  Max.   :1.0000

Ostatecznie model zawiera trzy zmienne bmi (zmienna dwuwartościowa), pregnancies oraz bloodpress:

## 
## Call:
## glm(formula = Class ~ factor(bmi) + pregnancies + bloodpress, 
##     family = binomial, data = p0)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7005  -0.9304  -0.5317   1.1209   2.1250  
## 
## Coefficients:
##               Estimate Std. Error z value         Pr(>|z|)    
## (Intercept)  -3.146685   0.530392  -5.933 0.00000000297890 ***
## factor(bmi)o  1.361376   0.193970   7.018 0.00000000000224 ***
## pregnancies   0.131255   0.025170   5.215 0.00000018397293 ***
## bloodpress    0.014276   0.007157   1.995           0.0461 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 931.94  on 723  degrees of freedom
## Residual deviance: 829.08  on 720  degrees of freedom
##   (4 observations deleted due to missingness)
## AIC: 837.08
## 
## Number of Fisher Scoring iterations: 4
##  (Intercept) factor(bmi)o  pregnancies   bloodpress 
##   0.04299441   3.90155884   1.14025872   1.01437873

Interpretacja: osoby otyłe mają prawie 4 razy większe szanse na cukrzycę. Każda kolejna ciąża zwiększa szansę na cukrzycę o 14% a wzrost ciśnienia krwi o 1mm Hg zwiększa te szanse o 1%.

Materiały dodatkowe (dla pasjonatów)

Two-Way Tables/RR/OR https://bookdown.org/mpfoley1973/data-sci/two-way-tables.html

Wisconsin Breast Cancer (Diagnostic) DataSet Analysis https://rstudio-pubs-static.s3.amazonaws.com/344010_1f4d6691092d4544bfbddb092e7223d2.html

High-Risk Populations: The Pimas of Arizona and Mexico https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4418458/

Machine Learning: Pima Indians Diabetes https://www.andreagrandi.it/2018/04/14/machine-learning-pima-indians-diabetes/

R Tutorial http://rstudio-pubs-static.s3.amazonaws.com/244535_7c4c81af10434a679f4a8a0d79cc4bb5.html

Repozytoria danych

https://figshare.com/ (https://en.wikipedia.org/wiki/Figshare)

https://datadryad.org/ ( https://en.wikipedia.org/wiki/Dryad ; a nymph or nature spirit who lives in trees and takes the form of a beautiful young woman.)

Data sharing in PLOS ONE: An analysis of Data Availability Statements https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0194768

https://journals.plos.org/plosone/s/data-availability

https://www.bmj.com/content/340/bmj.c181.long