Uživatel:Jkl~cswikiversity/Studuji biostatistiku

Úvod

V mém postgraduálnickém životě potřebuji občas něco statisticky zpracovat a protože jsem přestal používat MS-Windows, vzešla otázka co používat místo programu EpiInfo. V práci mám přístup k balíkům SPSS a MatLab, protože však chci podpořit open-source projekty (a pracovat i doma a legálně), byl jsem nucen začít samostudovat biostatistiku v Linuxu.

Aplikace

Dva nejrozšířenější balíky jsou Octave a R-project. Na fakultě mi před lety biostatistik doporučil spíš "Rko", což je odvozenina z jazyka S. Takže jsem na svém desktopu dal apt-get install octave r-base a s chutí do toho. Vzhledem k tomu, že se oba jazyky celkem liší, rozhodl jsem se začít s Rkem.

Rko I - první spuštění

Balíčkovací systém nahlásil, že je vše nainstalované. Faj, tak jak to spustit. Stáhnul jsem hafo PDF souborů (viz odkazy) a četl a četl. Všechny návody začínaly: poklepejte na ikonu na ploše ... Jenže v ubuntu žádná ikona. Nikde. Tak jsem musel chvíli googlovat abych nalezl informaci, že R se spouští tím, že v konzoli napíšete "R". Ach jo

První funkce

Každá proměnná je pole. Když chci do pole něco přiřadit, použiju

pole=c("oves","žito")

pokud chci pole vyplnit hodnotami, je tu

pole=seq(odkolika,dokolika,step)

nebo také platí, že

seq(1,7,1) = 1:7

Dělení polí: (1,2,3,4,5)/(1,2) odpovídá (1,2,3,4,5)/(1,2,1,2,1) Pojmenování sloupců: names(jmenopole)=pole_nadpisů

Výpis proměnné se prování buď jen jejím názvem - "pole" nebo - a v skriptech výlučně - funkcí "print(pole)"

Vytvoření tabulky

Buďto nadatlíme data ručně:

vzorek <- c("A", "B", "C", "D", "E")
hustota <- c(612, 671, 635, 601, 677)
tabulka = data.frame(vzorek, hustota)

... nebo máme data v nějaké tabulce. Asi většina z nás používá "české csv" - t.j. hodnoty s desetinou čárkou oddělené středníkem, anglicky comma-separated values ;-). Pak jsou následující ekvivalentní cesty načtení:

tabulka=read.table("tabulka.csv",sep=";",header = T)
tabulka=read.csv2("tabulka.csv")

read.csv2 je připraveno tak, aby četlo podle místních zvyklostí. Pokud chceme vytvořenou tabulku nebo jinou dvojrozměrnou datovou strukturu uložit, tak použijeme

write.csv2(tabulka, file="tabulka.csv")

Vypsání hodnot

Fajn, data máme načteny. Jak je vypsat:

• "tabulka" - vypíše celou tabulku
• tabulka$hustota - vypíše jen hustoty
• výpočet hodnoty (třeba BMI z jiné tabulky): "JinaTabulka$vaha/(JinaTabulka$vyska/100)^2"
• nechci-li opisovat název tabulky, lze použít attach(tabulka)/detach(tabulka)

Podmnožina

Občas (docela často) se stává, že potřebujeme vypsat jen část tabulky, například jen pacienty a ne kontroly. Nejjednodušší postup, jak dosáhnout vypsání podmnožiny je funkce subset:

subset(nazevtabulky,pacient==1,c(jmeno,vyska,cislobot))

Výpis podmnožiny je na obrazovku, není však nutné podmnožinu vypisovat na konzoli, výsledek funkce lze samozřejmě užít i jako vstupu pro jinou funkci. c(...) nám definuje, které sloupce se mají zahrnout do výstupu..

Pseudonáhodná čísla

• Pokud si chceme vsadit 5 ze 40, pak zavoláme "sample(1:40,5)"
• pomocí parametrů lze ovlivnit pravděpodobnost výberu jednotlivých hodnot - zatím to nepotřebuji, takže víc neřeknu

Rko II - Základní statistika

Testy normality

Abychom zjistili, jestli mají data normální rozložení, potřeujeme na to nějaký test. Podle normality se pak rozhodneme, který test použít dál ...

Shapirův-Wilkův test

• shapiro.test(x) -> výsledek je W a p
• Hodnota W je ve své podstatě korelační koeficient, který nám říká jak těsně naše data korelují s křivkou normálního rozdělení. p-value pak udává jaké chyby se dopouštíme, pokud zamítneme nulovou hypotézu H0: „Zkoumaná data pochází ze základního souboru s normálním rozdělením.“. Čili pro p<alfa zamítáme hypotézu o normálním rozdělení dat

Testování hypotéz - kudy dál

Takže jsme udělali Shapirův-Wilkův test a už víme, jestli jsme (teda pardon - jestli jsou data) normální. A pak už s trochou drzosti můžeme použít následující rozcestníkovou tabulku:

Jakou použít metodu na testování hypotézy

Data	normální rozložení	není normalita
dvě skupiny, shoda středních hodnot	t- test pro dva nezávislé výběry	Mann-Whitney či Dvouvýběrový Wilcoxonův test
více dat	ANOVA	Kruskalův-Wallisův test
párové porovnání	párový t-test	test mediánový, znaménkový, párový test Wilcoxonův
rozložení proměnné u dvou nebo více skupin	testy dobré shody	?

Mimo to na porovnání 2 kategoriálních proměnných lze použít Pearsonův chí kvadrát test, Fisherův test, Spearmanův korelační koeficient (nezávislost proměnných) či korelační koeficient Pearsonův.

Jestli máme dobrou velikost vzorku nám řekne power analysis (jak se dělá ??).

Vývoj kontinuální proměnné řeší regresní analýza, diskriminační nalýza (jestli už lék zabral), analýza časových dat - jestli je nějaký bod zvratu

Faktorová analýza.

Test homogenity - F test

• var.test(x, y, ratio = 1,alternative = c("two.sided", "less", "greater"),conf.level = 0.95, ...)
• p-value udává, jaké chyby se dopustíme, když zamítneme nulovou hypotézu H0: „Zkoumané soubory nevykazují statisticky významné rozdíly mezi rozptyly.“
• opět by mělo být p<alfa

T-test (jednovýběrový parametrický)

• t.test(x, y = NULL,alternative = c("two.sided", "less", "greater"),mu = 0, paired = FALSE, var.equal = FALSE,conf.level = 0.95, ...)
• t.test(a, b, var.equal = T) - T test rovnosti rozptylů
  • je to Studentův t-test.

Jednovýběrový neparametrický test

wilcox.test(x, mu =, alternative =, conf.level =)

Dvouvýběrové testy

Dvovýběrové testy parametrické

F-test na porovnání rozptylů (normalita obou výběrů):

var.test(x, y, ratio =)
var.test(y factor, ratio =)

Dvouvýběrový T-test, shodné rozptyly (normalita obou výběrů):

t.test(x, y, mu =, var.equal = T)
t.test(y factor, mu =, var.equal = T)

Dvouvýběrový T-test (zobecněný, Welch), různé rozptyly (normalita obou výběrů)

t.test(x, y, mu =, var.equal = F)
t.test(y factor, mu =, var.equal = F)

Neparametrické

Dvouvýběrové testy

Dvouvýběrový Wilcoxonův test

Dvouvýběrový Wilcoxonův test (spojité výběry):

wilcox.test(x, y, mu =)
wilcox.test(y factor, mu =)

Mann-Whitney

• nevadí mu "tied rank" - více měření se stejnou hodnotou (třeba dva stejně vysocí pacienti)
• označuje se jako U-test

Kolmogorovův-Smirnovův test

ks.test(x, y)

Párové testy

Párový t-test

normalita rozdílu páru:

t.test(x, y, mu =, paired = T)

Párový Wilcoxonův test

spojitost a přibližná symetrie rozdílu párů:

wilcox.test(x, y, mu =, paired = T)

Jednofaktorová ANOVA

• analýza rozptylu neboli ANOVA je statistický test, který testuje nulovou hypotézu H0: Střední hodnoty jednotlivých výběrů jsou shodné
• Pokud je Pr<alfa, tak zavrhujeme nulovou hypotézu a pokračujeme metodou vícenásobného porovnání - např.Tukeyova

Příklad: máme tabulku o dvou sloupcích, v prvním je číslo bot a ve druhém dosažené vzdělání (kategoriální veličina 1-4. Ptáme se, jestli dosažené vzdělání nějak souvisí s velikostí nohy. Vytvoříme soubor skola.csv

"bota";"skola"
42;"ZŠ"
38;"ZŠ"
40;"ZŠ"
31;"SŠ"
40;"SŠ"
44;"SŠ"
42;"SŠ"
33;"VŠ"
40;"VŠ"
35;"VŠ"
45;"VŠ"
44;"VŠ"
34;"VŠ"
46;"VŠ"
44;"VŠ"

a pustíme R:

> b=read.csv2("skola.csv")
> anova(lm(bota~skola, data=b))
Analysis of Variance Table 
.
Response: bota
          Df  Sum Sq Mean Sq F value Pr(>F)
skola      2   2.108   1.054  0.0409 0.9601
Residuals 12 309.625  25.802

Takže je to dobré, žádná hvězdička, takže se příjimačky nebudou dělat podle čísla bot :-) . Pokud chcete, aby Vám něco vyšlo, tak si stáhněte tenhle soubor a dejte

> mojedata=read.csv2("data-anova.txt")
> anova(lm(hustota~lokalita, data=mojedata))

Na jednoprocentní hladině tam vychází Lokalita 3.

Tukey

• Pokud 0 nepatří do intervalu (lwr, upr), pak existuje statisticky významný rozdíl mezi středními hodnotami uvedené dvojice skupin.
• čili pokud je to statisticky významné, tak tam nula NENÍ

Kruskal-Wallis

"Kruskal-Wallis rank sum test"

kruskal.test(x, g, ...)

• x - jeden či více numerických vektorů
• g - vektor či faktorový objekt (???)

Tímto testem se testuje hypotéza, že poloha parametrů rozložení x jsou stejnná ve všech skupinách. Ukázka (doufám, že správně chápu vztah copyrightu k citování ??)

##Hollander & Wolfe (1973), 116.
## Mucociliary efficiency from the rate of removal of dust in normal
## subjects, subjects with obstructive airway disease, and subjects
## with asbestosis.
x <- c(2.9, 3.0, 2.5, 2.6, 3.2) # normal subjects
y <- c(3.8, 2.7, 4.0, 2.4) # with obstructive airway disease
z <- c(2.8, 3.4, 3.7, 2.2, 2.0) # with asbestosis
kruskal.test(list(x, y, z))
## Equivalently,
x <- c(x, y, z)
g <- factor(rep(1:3, c(5, 4, 5)),
labels = c("Normal subjects",
"Subjects with obstructive airway disease",
"Subjects with asbestosis"))
kruskal.test(x, g)
## Formula interface.
require(graphics)
boxplot(Ozone ~ Month, data = airquality)
kruskal.test(Ozone ~ Month, data = airquality)

Rko III - Práce s grafikou

Základní typy grafů

plot(x,y)

kreslí bubliny

boxplot

• kreslí krabice ;-)
• Střední čárka v krabici představuje medián. Hranice krabice pak představují 1. a 3. kvartil.

boxplot(x, y, main = "Krabice", ylab = "rozměr [mm]", names = NULL, col = "lightgray")

• parametry - boxplot(x, ..., range = 1.5, width = NULL, varwidth = FALSE, notch = FALSE, outline = TRUE, names, plot = TRUE,border = par("fg"), col = NULL, log = "",pars = list(boxwex = 0.8, staplewex = 0.5, outwex = 0.5),horizontal = FALSE, add = FALSE, at = NULL)

pie

2D koláčový graf

barplot

sloupcový graf. V podstatě no comment.

hist

• kreslí histogram
• syntax: hist(x, breaks = "Sturges", freq = NULL, probability = !freq, include.lowest = TRUE, right = TRUE, density = NULL, angle = 45, col = NULL, border = NULL, main = paste("Histogram of" , xname), xlim = range(breaks), ylim = NULL, xlab = xname, ylab, axes = TRUE, plot = TRUE, labels = FALSE, nclass = NULL, ...)
• příklad na histogram - histogram s proloženou gaussovkou:

hist(y, freq = F, ylim = c(0, 0.5))a <- seq(-3.5, 3.5, 0.1),lines(a, dnorm(a))
a <- seq(-3.5, 3.5, 0.1)
lines(a, dnorm(a))

Ukládání výstupu

Souhrnně se

jpeg(filename = "Rplot%03d.jpg", width = 480, height = 480,pointsize = 12, quality = 75, bg = "white", res = NA)
plot(1,1)
...
dev.off()

Stejně funguje i funkce bmp,png, pdf, postscript. SVG bohužel - zdá se - moje verze Rka neumí. Na druhou stranu, když by to bylo potřeba nebude asi problém příslušnou rutinu napsat.

Rko IV - Skriptování

• "rm(list=ls())" vymaže všechny proměnné z paměti
• spuštění skriptu:
  • v kódu: source("/home/jim/psych/adoldrug/partyuse1.R")
  • z řádky: R CMD BATCH /home/jim/psych/adoldrug/partyuse1.R
• cd = setwd(dir) (vs. getwd() )
• sink jako příkaz pro otevření výstupního souboru

Rko V - externí moduly

Neljslavnějším externím modulem je Rcmdr, který se mi ale nedaří automaticky instalovat. Proto si postup ukážeme na modulu car:

install.packages("car", dependencies = TRUE)

Program se zeptá na mirror, a balíky nainstaluje (v UNIXu) do /usr/local/lib/R/site-library (je dobré, aby instalující uživatel měl právo zápisu do tohoto adresáře ;-) ).

Pokud chceme balík použít, napíšeme

library(car)

Rko VI - modelové příklady

• 1. praktický příklad

Další užitečné programy

• GNU plot - kreslení grafů

Použitá a doporučená literatura

• Karel Zvára: Programové prostředí R pro biology - perfektně se čte, nepředpokládá hluboké znalosti statistiky !
• Vladislav Bína, Arnošt Komárek, Lenka Komárková: Jak na jazyk R, instalace a základní příkazy
• Michal Kulich: R pro samouky
• odkaz na mnoho dalších užitečných zdrojů
• Höschl C, Libiger J, Švestka J (ed.): Psychiatrie, Tigis Praha 2002, 245-6
• Arnaud Delorme:STATISTICAL METHODS