Uživatel:Jkl~cswikiversity/Studuji biostatistiku/R/Příklad 1
Spadnul mi FF, takže se obsah této stránky ztratil. Klatě, Časem doplním ...
Vstupy
[editovat]- soubor vstup.csv (klasické "české" CSV, neboli hodnoty oddělené středníkem s desetinnou čárkou). Předpokládám soubor s pacienty a zdravými dobrovolníky, identifikátor Case označuje zdravé dobrovolníky (0) a nemocné (1).
- soubor head.txt
tabulka=read.csv2("vstup.csv")
case=subset(tabulka,Case==1)
ctrl=subset(tabulka,Case==0)
Předzpracování
[editovat]head vstup.csv -n1|sed -e s/";"/"\n"/g>promenne.list
cp promenne.list nezav.list
ze souboru nezav.list odstraníme závisle proměnné a ponecháme jen nezávisle proměnné
diff nezav.list promenne.list |grep ">"|sed -e s/"> "/""/g>zav.list
Nyní máme tři soubory - všechny, závisle a nezávisle proměnné. ještě budeme potřebovat soubor bez kategoriálních veličin a smetí jako je třeba identifikátor pacienta. Ten vytvoříme z promenne.list a pojmenujeme nekategorialni.list
Test normality
[editovat]cp head.txt normalita.R
sed -e s/"\(.*\)"/"shapiro.test\(case\$\1\)\nshapiro.test\(ctrl\$\1\)"/ zav.list>>normalita.R
R CMD BATCH normalita.R
egrep "p-value|data" normalita.Rout>normalita_zajimave.txt
cat normalita_zajimave.txt |./spoj_po_dvou_radcich|sed -e s/"data:"/""/g -e s/"\\$"/";"/g -e s/"W = "/";"/g -e s/", p-value = "/";"/g -e s/" "/""/g|awk -F ";" '{print $4 ";" $1 ";" $2 ";" $3}'|sort -n>normalita.csv
V souboru normalita_zajimave.txt nyní máme "zajímavé" řádky z výstupu Rka, t.j. co se porovná, W a p-hodnotu. Podle p hodnoty si vytvoříme řádky normalni.list a nenormalni.list se seznamem příslušných závisle proměnných.
Poslední řádek použije níže popsaný a použitý skript Zdroják:spoj po dvou řádcích a vytvoří CSV soubor ve formátu p-value;(pod)soubor;proměnná;W-hodnota
Test shody rozptylů
[editovat]Nyní nás budou zajímat jen proměnné s normálním rozložením. T-test se totiž spouští s parametrem varequal, který značí, jestli se má předpokládat rovnost rozptylů nebo ne. Výchozí nastavení je, že ne. takže...
cp head.txt rovnost_rozptylu.R
sed -e s/"\(.*\)"/"var.test\(case\$\1,ctrl\$\1\)"/ normalni.list>>rovnost_rozptylu.R
R CMD BATCH rovnost_rozptylu.R
A dále se v rovnost_rozptylu.Rout podíváme, které testy skončily s p>0,05 - u nich si můžeme dovolit varequal=TRUE, jinak zůstaneme u FALSE
Korelace
[editovat]I když nám to asi pravověrní statistici budou mít za zlé, už z didaktických důvodů se vydáme na ryby. Nejprve si připravíme všechny kombinace proměnných pomocí programu Uživatel:Jkl/Studuji cpp/Kombinace proměnných. Z toho padá soubor typu [a,b]. Pak trochu [sed]u a jdeme na to:
./kombinace_promennych nekategorialni.list nekategorialni.list >korelace.list
cp head.txt ryby.all.R
cp head.txt ryby.case.R
sed -e s/"\[\(.*\),\(.*\)\]"/"cor.test\(tabulka\$\1,tabulka\$\2\,use=\"complete.obs\")"/g korelace.list >>ryby.all.R
sed -e s/"\[\(.*\),\(.*\)\]"/"cor.test\(case\$\1,case\$\2\,use=\"complete.obs\")"/g korelace.list >>ryby.case.R
R CMD BATCH ryby.all.R
R CMD BATCH ryby.case.R
Na výstupu máme zase soubor s příponou .Rout, který je dlouhý a nepřehledný. Použijeme program pro přepracování korelace - přeložíme jako "korela" a spustíme v kombinaci se sedem a awk:
./korela ryby.all.Rout|awk -F ";" '{print $1 ";xxxxx" $2 "xxxxx;" $3 ";" $4 ";" $5 ";" $6 ";" $7";" $8}'|sed -e s/"xxxxx \(.*\) \(.*\) xxxxx"/"\1;\2"/g -e s/" and "/";"/g -e s/"data\:"/";"/g -e s/"t = "/""/g -e s/"df = "/""/g -e s/"p-value = "/""/g -e s/","/";"/g -e s/"\."/","/g -e s/" "/""/g -e s/"xxxxx"/""/g -e s/";;;;;"/""/g>korela.csv
Na výstupu máme .csv soubor s řadou zajímavých proměnných Je možné ho dále elegantně zpracovat na kombinatorickou tabulku typu "každý s každým" například cvičným C++ programem nextstep, který využívá dvojrozměrnou mapu. Ke zvýraznění použijeme například podmíněné formátování v OpenOffice ...
Porovnání hodnot pacienti - zdraví
[editovat]Tak jak máme předpřipravená data můžeme použít některý z dvojvýběrových testů. Pokud budeme mít někde rovnost rozptylů, můžeme si dovolit příslušné příkazy modifikovat s varequal. Taky lze přemýšlet o snižování nebo zvyšování (nikoliv double-tailed).
sed -e s/"\(.*\)"/"t.test\(case\$\1,ctrl\$\1\,paired=FALSE)"/g normalni.list >normalni.R
sed -e s/"\(.*\)"/"wilcox.test\(case\$\1,ctrl\$\1\,paired=FALSE)"/g nenormalni.list >nenormalni.R
R CMD BATCH nenormalni.R
R CMD BATCH normalni.R
Výsledky jsou opět celkem nečitelné, zase nám pomůže prográmek v C++ nebo vlastní oko ;-) Protože normálních dat jsem měl málo, ukáži na nenormálních jak elegantně na to:
egrep "p-value =|data:" nenormalni.Rout|./spoj_po_dvou_radcich|sed -e s/"data:"/""/g -e s/" and "/";"/g -e s/"W = "/";"/g -e s/", p-value = "/";"/g -e s/" "/""/g|awk -F ";" '{print $4 ";" $1 ";" $2}'|sort -n>nenormalni.csv
Vzhledem k tomu, že mne nezajímá W hodnota, tak jsem trochu zpřeházel sloupce a ve finále seřadil výsledky podle p hodnoty vzestupně. Prohlížíme pak třeba v oocalc zezhora ;-)