Uživatel:Jkl~cswikiversity/Studuji biostatistiku

ÚvodEditovat

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.

AplikaceEditovat

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íEditovat

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í funkceEditovat

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í tabulkyEditovat

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í hodnotEditovat

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žinaEditovat

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á číslaEditovat

  • 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í statistikaEditovat

Testy normalityEditovat

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 testEditovat

  • 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álEditovat

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 testEditovat

  • 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ý)Editovat

  • 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ý testEditovat

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

Dvouvýběrové testyEditovat

Dvovýběrové testy parametrickéEditovat

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éEditovat

Dvouvýběrové testyEditovat
Dvouvýběrový Wilcoxonův testEditovat

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

wilcox.test(x, y, mu =)
wilcox.test(y factor, mu =)
Mann-WhitneyEditovat
  • 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 testEditovat
ks.test(x, y)

Párové testyEditovat

Párový t-testEditovat

normalita rozdílu páru:

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

Párový Wilcoxonův testEditovat

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

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

Jednofaktorová ANOVAEditovat

  • 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 :-) [1]. 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.

TukeyEditovat

  • 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-WallisEditovat

"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 grafikouEditovat

Základní typy grafůEditovat

plot(x,y)Editovat

kreslí bubliny

boxplotEditovat

  • 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)

pieEditovat

2D koláčový graf

barplotEditovat

sloupcový graf. V podstatě no comment.

histEditovat

  • 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ýstupuEditovat

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íEditovat

  • "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í modulyEditovat

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říkladyEditovat

Další užitečné programyEditovat

Použitá a doporučená literaturaEditovat

Poznámky pod čarouEditovat

  1. Já vím, že je soubor moc malý na to, aby něco signifikantně vyšlo. Delší tabulka by ale vypadala nehezky ...