Tady jsou příklady ze Základů biostatistiky prof. Herbena, které jsem (vy)řešil v R. Kromě posledního příkladu jsem si tím celkem jist (pro něj bych uvítal nějakou dobrou radu:-). Původně jsem měl tendenci napsat i obecný návod (základy práce v R), ale nějak jsem se k tomu nedokopal... Krom toho už se to snad bere na přednášce. Rád bych jen připomenul linuxový program Rkward (grafická nadstavba nad R) a doplněk do OpenOffice.org R4CalcAž budu mít čas a energii, rád bych v R udělal vzorová řešení i příkladů z biostatistiky II (Canoco) a z Fenetiky a kladistiky I.
Příklad č. 1. Otázka: Liší se přežívání semenáčků zkoumané trávy na louce sekané a pasené?
Celkem bylo sledováno 54 semenáčků o stejné výchozí velikosti na sekané louce a 68 semenáčků na pasené louce. Za měsíc zbylo 12 semenáčků na sekané louce a 8 semenáčků na pasené louce. > tab<-data.frame() > fix(tab) > tab var1 var2 1 54 68 2 12 8 > chisq.test(tab) Pearsons Chi-squared test with Yates continuity correction data: tab X-squared = 1.1367, df = 1, p-value = 0.2864
Příklad č. 2. Otázka: Liší se vzájemně mezi sebou rychlosti růstu semenáčků zkoumané trávy na louce sekané, pasené a neobhospodařované?
Byly sledovány semenáčky o stejné výchozí velikosti v těchto typech obhospodařování. Velikost semenáčků na konci pokusu: louka sekaná:18,17,18,16,21, 18,16,18,19,17 louka pasená: 26,25,26,29,22, 30,27,23,24,29 louka neobhospodařovaná:13,9,8,10,11,9,8,10,14 >vyska<-c(po řadě všechny hodnoty) >faktor<-factor(rep(c("A","B","C"),c(10,10,9))) >boxplot(vyska~faktor) >hist(vyska) >model<-aov(vyska~faktor) >plot(resid(model)) >plot(model,2) >summary(model) >anova(model)
Příklad č. 3. Otázka: Liší se mezi sebou významně dvě populace zkoumané trávy v rychlosti růstu?
Z každé populace bylo odebráno pět genotypů, a každý byl vegetativně namnožen, odnože o standardní velikosti byly pěstovány v čtyřech květináčích. Po měsíci pokusu byly rostliny odebrány a zváženy. Biomasa odnoží (relativní přírůstek) Populace 1 Klon 1: 1.4,1.5,1.2,1.5 Populace 1 Klon 2: 1.8,1.5,1.3,1.7 Populace 1 Klon 3: 1.8,1.9,1.8,1.7 Populace 1 Klon 4: 1.1,1.3,1.4,1.6 Populace 1 Klon 5: 1.9,1.8,1.7,1.9 Populace 2 Klon 1: 2.0,2.2,1.9,1.8 Populace 2 Klon 2: 1.8,1.9,1.8,1.7 Populace 2 Klon 3: 2.0,1.4,1.9,1.8 Populace 2 Klon 4: 2.2,2.3,2.0,1.8 Populace 2 Klon 5: 2.3,2.0,1.8,1.9 > biomasa<-c(1.4,1.5,1.2,1.5,1.8,1.5,1.3,1.7,1.8,1.9,1.8,1.7,1.1,1.3,1.4,1.6,1.9,1.8,1.7,1.9, 2.0,2.2,1.9,1.8,1.8,1.9,1.8,1.7,2.0,1.4,1.9,1.8,2.2,2.3,2.0,1.8,2.3,2.0,1.8,1.9) > hist(biomasa) > klon<-factor(rep(rep(c(1,2,3,4,5),c(4,4,4,4,4)),2)) > pop<-factor(rep(c(1,2),c(20,20))) > boxplot(biomasa~klon) > boxplot(biomasa~pop) > boxplot(biomasa~pop*klon) > model<-aov(biomasa ~ pop + pop/klon) > plot(resid(model)) > plot(model,2) > summary(model) > anova(model)
Příklad č. 4. Otázka: liší se průměrná velikost křídla mezi dvěma druhy octomilek?
Data jsou z přirozených populací. Od každého druhu byly sledovány čtyři populace, v každém z nich měřeno šest jedinců. Hodnoty jsou v relativních jednotkách. Druh 1 Populace 1: 14,15,12,15,18,13 Druh 1 Populace 2: 18,15,13,17,19,14 Druh 1 Populace 3: 18,19,18,17,17,12 Druh 1 Populace 4: 11,13,14,16,19,16 Druh 2 Populace 1: 20,22,19,18,23,20 Druh 2 Populace 2: 18,19,18,17,19,21 Druh 2 Populace 3: 20,14,19,18,18,19 Druh 2 Populace 4: 22,23,20,18,20,19 > kridlo<-c(14,15,12,15,18,13,18,15,13,17,19,14,18,19,18,17,17,12,11,13,14,16,19,16,20,22, 19,18,23,20,18,19,18,17,19,21,20,14,19,18,18,19,22,23,20,18,20,19) > kridlo > length(kridlo) > druh<-factor(rep(c(1,2),c(24,24))) > druh > length(druh) > pop<-factor(rep(rep(c(1,2,3,4),c(6,6,6,6)),2)) > pop > hist(kridlo) > boxplot(kridlo~druh) > boxplot(kridlo~druh) > boxplot(kridlo~druh*pop) > octom<-aov(kridlo~druh+druh/pop) > octom > plot(octom,2) > plot(resid(octom)) > summary(octom) > anova(octom)
Příklad č. 5. Otázka: Liší se velikost dospělých octomilek v závislosti na výživě a genotypu?
Samičky od obou genotypů byly pěstovány na třech druzích výživy. Jedinci následující generace byli změřeni. Hodnoty jsou v relativních jednotkách. (Doplňující otázka: navrhněte, jak správně založit takový pokus.) genotyp1, výživa 1: 18,17,18,16,21,20,18,19 genotyp1, výživa 2: 26,25,26,29,22,22,20,23 genotyp1, výživa 3: 28,21,29,31,26,25,23,28 genotyp2, výživa 1: 13,19,18,15,16,18,14,17 genotyp2, výživa 2: 18,16,18,19,17,19,20,18 genotyp2, výživa 3: 22,26,17,19,27,21,24,23 > sam<-c(18,17,18,16,21,20,18,19,26,25,26,29,22,22,20,23,28,21,29,31,26,25,23,28, 13,19,18,15,16,18,14,17,18,16,18,19,17,19,20,18,22,26,17,19,27,21,24,23) > length(sam) > hist(sam) > gen<-factor(rep(c(1,2),c(24,24))) > gen > length(gen) > vyziva<-factor(rep(rep(c(1,2,3),c(8,8,8)),2)) > vyziva > length(vyziva) > boxplot(sam~vyziva) > boxplot(sam~gen) > boxplot(sam~gen*vyziva) > speky<-aov(sam~gen*vyziva) nebo > speky<-aov(sam~gen+vyziva+gen:vyziva) (je to totožné) > speky > plot(resid(speky)) > plot(speky,2) > summary(speky) > anova(speky)
6. Otázka: Liší se podíl včelích dělnic infikovaných roztočem mezi třemi roji?
Roj 1: odebráno 100 včel, z toho je 11 infikovaných, Roj 2: odebráno 100 včel, z toho je 19 infikovaných. Roj 2: odebráno 100 včel, z toho je 28 infikovaných. > roje<-data.frame() > fix(roje) > roje roj1 roj2 roj3 1 11 19 28 2 89 81 72 > chisq.test(roje) Pearsons Chi-squared test data: roje X-squared = 9.2761, df = 2, p-value = 0.009676
Příklad č. 7. Otázka: Liší se poměr pohlaví u studovaného druhu mravence (mezi samečky a plodnými samičkami) mezi dvěma koloniemi?
Kolonie 1: 58 samců, 32 samic Kolonie 2: 38 samců, 27 samic > mravenci<-data.frame() > fix(mravenci) > mravenci samci samice 1 58 32 2 38 27 > chisq.test(mravenci) Pearsons Chi-squared test with Yates continuity correction data: mravenci X-squared = 0.3474, df = 1, p-value = 0.5556
Příklad č. 8. Otázka: Jak a zda závisí výška rostliny fazole po dvou týdnech kultivace na velikosti semene?
Semena byla zvážena, rostliny vysety po jednom semeni do k větináčů do standardních podmínek. výška po řadě: 58, 66, 64, 45, 36, 21, 32, 40, 19, 25 biomasa semene po řadě: 13, 12, 13, 10, 10, 9, 11, 11, 9, 10 > vyska<-c(58, 66, 64, 45, 36, 21, 32, 40, 19, 25) > biomasa<-c(13, 12, 13, 10, 10, 9, 11, 11, 9, 10) > length(vyska) [1] 10 > length(biomasa) [1] 10 > hist(vyska) > hist(biomasa) > plot(vyska,biomasa) > abline(lm(vyska~biomasa)) > fazole<-lm(vyska~biomasa) > fazole Call: lm(formula = vyska ~ biomasa) Coefficients: (Intercept) biomasa -72.47 10.47 > plot(fazole,2) > summary(fazole) Call: lm(formula = vyska ~ biomasa) Residuals: Min 1Q Median 3Q Max -10.694 -4.913 -1.724 2.923 12.837 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -72.469 20.509 -3.534 0.007692 ** biomasa 10.469 1.883 5.559 0.000535 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 8.337 on 8 degrees of freedom Multiple R-Squared: 0.7944, Adjusted R-squared: 0.7687 F-statistic: 30.91 on 1 and 8 DF, p-value: 0.0005351
Příklad č. 9. Otázka: Jak a zda závisí počet jedinců studovaného druhu mšice na velikosti rostliny a vzdálenosti ke konspecifické rostlině?
V terénu bylo náhodně vybráno 15 rostlin. Na nich byl spočten počet jedinců mšic a změřen parametr velikosti (průměr stonku) a vzdálenost k sousedovi. počet jedinců mšic po řadě: 55, 84, 20, 60, 47, 42, 35, 65, 32, 90, 75, 60, 40, 55, 87 velikost rostliny po řadě: 8, 10, 4, 6, 6, 5, 5, 8, 4, 10, 6, 7, 4, 6, 11 vzdálenost ke konspecifickému sousedu po řadě: 14, 9, 18, 12, 13, 11, 12, 15, 12, 8, 9, 14, 18, 14, 15 > pocet<-c(55, 84, 20, 60, 47, 42, 35, 65, 32, 90, 75, 60, 40, 55, 87) > length(pocet) [1] 15 > hist(pocet) > velikost<-c(8, 10, 4, 6, 6, 5, 5, 8, 4, 10, 6, 7, 4, 6, 11) > length(velikost) [1] 15 > hist(velikost) > vzdalenost<-c(14, 9, 18, 12, 13, 11, 12, 15, 12, 8, 9, 14, 18, 14, 15) > length(vzdalenost) [1] 15 > hist(vzdalenost) > cor.test(velikost,vzdalenost) Pearsons product-moment correlation data: velikost and vzdalenost t = -1.3395, df = 13, p-value = 0.2034 (nesig. >> proměnné nezávislé) alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: -0.7302401 0.1996270 sample estimates: cor -0.3482482 > msice<-lm(pocet~velikost+vzdalenost) > msice Call: lm(formula = pocet ~ velikost + vzdalenost) Coefficients: (Intercept) velikost vzdalenost 29.266 7.462 -1.743 > plot(msice,2) > summary(msice) Call: lm(formula = pocet ~ velikost + vzdalenost) Residuals: Min 1Q Median 3Q Max -10.65684 -5.79742 0.05929 4.13654 16.65114 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 29.2657 14.8961 1.965 0.0730 . velikost 7.4621 1.0724 6.958 1.52e-05 *** vzdalenost -1.7433 0.8216 -2.122 0.0554 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 8.609 on 12 degrees of freedom Multiple R-Squared: 0.857, Adjusted R-squared: 0.8332 F-statistic: 35.96 on 2 and 12 DF, p-value: 8.547e-06
Příklad č. 10. Otázka: Jak a zda závisí biomasa obilek na průměru stébla nad povrchem půdy?
Bylo vybráno náhodně dvacet rostlin. U nich byly změřeny oba parametry. Průměr stébla: 6, 7, 8, 6, 4, 6, 5, 7, 8, 9, 9, 8, 6, 5, 6, 6, 4, 8, 7, 5 Biomasa obilek: 8,17,22, 10,2,9, 4, 15,24, 33,40,28,8, 5, 10,8, 1, 23, 14,6 > steblo<-c(6, 7, 8, 6, 4, 6, 5, 7, 8, 9, 9, 8, 6, 5, 6, 6, 4, 8, 7, 5) > length(steblo) [1] 20 > hist(steblo) > obilka<-c(8,17,22, 10,2,9, 4, 15,24, 33,40,28,8, 5, 10,8, 1, 23, 14,6) > length(obilka) [1] 20 > hist(obilka) > plot(obilka,steblo) > biomasa<-lm(obilka~steblo) > plot(resid(biomasa)) > plot(biomasa,2) > biomasa3<-lm(sqrt(obilka)~steblo) > plot(biomasa3,2) > plot(resid(biomasa3)) > biomasa3 Call: lm(formula = sqrt(obilka) ~ steblo) Coefficients: (Intercept) steblo -2.616 0.944 > summary(biomasa3) Call: lm(formula = sqrt(obilka) ~ steblo) Residuals: Min 1Q Median 3Q Max -0.24996 -0.17474 -0.07575 0.13169 0.44513 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -2.61573 0.23080 -11.33 1.26e-09 *** steblo 0.94391 0.03464 27.25 4.36e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.2271 on 18 degrees of freedom Multiple R-Squared: 0.9763, Adjusted R-squared: 0.975 F-statistic: 742.6 on 1 and 18 DF, p-value: 4.364e-16
Příklad č. 11. Liší se obsah celkového dusíku v keřovitém stepním porostu v závislosti na tom, zda (1) je půda ve svrchním nebo spodním horizontu, (2) zda je v blízkosti keře nebo v otevřeném porostu?
(Doplňující otázka: navrhněte, jak korektně sebrat tato data.) Půda v blízkosti keře: (a) svrchní horizont: 56, 48, 60, 49, 52, 67, 62; (b) spodní horizont: 39, 37, 46, 51, 49, 41, 44 Půda v otevřeném porostu: (a) svrchní horizont: 51, 48, 41, 38, 60, 45, 52; (b) spodní horizont: 28, 27, 38, 16, 31, 22, 18 > dusik<-c(56, 48, 60, 49, 52, 67, 62, 39, 37, 46, 51, 49, 41, 44, 1, 48, 41, 38, 60, 45, 52, 28, 27, 38, 16, 31, 22, 18) > length(dusik) [1] 28 > hist(dusik) > stanoviste<-factor(rep(c("K","P"),c(14,14))) > stanoviste [1] K K K K K K K K K K K K K K P P P P P P P P P P P P P P Levels: K P > length(stanoviste) [1] 28 > horizont<-factor(rep(rep(c("H","D"),c(7,7)),2)) > horizont [1] H H H H H H H D D D D D D D H H H H H H H D D D D D D D Levels: D H > length(horizont) [1] 28 > boxplot(dusik~stanoviste) > boxplot(dusik~horizont) > boxplot(dusik~horizont*stanoviste) > interaction.plot(dusik,horizont,stanoviste) > kere<-aov(dusik~horizont*stanoviste) > plot(resid(kere)) > plot(kere,2) > kere Call: aov(formula = dusik ~ stanoviste * horizont) Terms: stanoviste horizont horizont:stanoviste Residuals Sum of Squares 1316.5714 1989.1429 11.5714 2971.1429 Deg. of Freedom 1 1 1 24 Residual standard error: 11.12644 Estimated effects may be unbalanced > summary(kere) Df Sum Sq Mean Sq F value Pr(>F) stanoviste 1 1989.14 1989.14 16.0677 0.0005158 *** horizont 1 1316.57 1316.57 10.6349 0.0033112 ** horizont:stanoviste 1 11.57 11.57 0.0935 0.7624468 Residuals 24 2971.14 123.80 --- Signif. Codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (vyšlo mi to signifikantně, ale číselně úplně jinak než Herbenovi, jiný model mě nenapadá...) Tady je výstup S-plus se stejnými daty: > summary(model) Df Sum of Sq Mean Sq F Value Pr(F) ker 1 1235.571 1235.571 25.95349 0.00003273 horizont 1 2091.571 2091.571 43.93398 0.00000074 ker:horizont 1 165.143 165.143 3.46887 0.07482072 Residuals 24 1142.571 47.607 Vše samozřejmě bez záruky a volně šířitelné pod open licencí GPL! :-))
