前回は生理学系？のデータセット。今回は病気ネタです。ブラジルでcarvao da cana-de-acucar という恐ろしい病気が流行っていた（いる）そうな。罹るのはサトウキビっす。新芽が黒い炭のようになってしまって収量激減。日本では「黒穂病」と呼ばれているみたい。抵抗力のある品種を選抜するためのデータらしいです。
※「データのお砂場」投稿順Indexはこちら
R言語向けのパッケージ boot には、サンプルデータセット多数が含まれております。それを ABC順に「舐めて」ます。今回のサンプルデータセットはcaneとな。
cane データセット
sugar-cane(サトウキビ) に関する農学的なサンプルデータセットのようです。データはブラジル産。ブラジルで猛威を振るっていた（いる？）サトウキビ黒穂病（carvao da cana-de-acucar）についてのデータです。検索すればヒットしてくるので多分今もサトウキビ生産者様には頭の痛い問題なのだと思います。ただ、ブラジルの文献はポルトガル語っす。Google様のお力にすがってもイマイチ腑に落ちませぬ。
日本で調べてみると、沖縄、伊是名村役場様の以下文書が目につきました。
どうもサトウキビの生産農家宛てに注意を呼び掛けるチラシみたいです。さらに検索すると沖縄県でも数年に１度くらい流行があるようで警戒情報が掲載されてました。日本でもあるみたいっす。
さてデータセットは実験農園？的なところに45種類（バラエティ？）のサトウキビを植えて、出てきた新芽の数、そのうち病気のやつなど数えているのです。その45のバラエティの中で、抵抗力のあるやつ、かかりやすいやつ、中間と評価するためみたいっす。４５種をランダムに区画ごとに植えたものを、ブロックA、B、C、Dの４ブロックに植えているみたいなので合計180点のデータです。
-
- n 各区画の「芽」の総数
- r 病気の「芽」の数
- x 植えた「茎（苗？）」の数
- var 各区画の品種（特性？）を示す数字
- block ブロックの種別（A,B,C,D）
正直、分からんのがxっす。各区間とも50本の茎（苗？）を病原体に晒しているのですが、そのくせ植えているのは50本中 x 本だけです。これもバラつかせるための技？素人にはよう分からんです。
先ずは生データ
いつものようにbootパッケージとggplot2パッケージをライブラリとしてとりこんであります。その前提でまずは生データをロード。
これをどうしろ、って感じがしないでもないデータっす。
病気に罹った比率を求めてみる
仕方がないので、素人が勝手な思い込みで「芽」の中での罹患率みたいなものを求めてみます。処理はこんな感じ。
cane.work <- cane
cane.work$ratio <- cane$r / cane$n
cane.var <- aggregate(cane.work$ratio, by=list(var=cane.work$var), mean)
colnames(cane.var) <- c("var", "ratio")
plot(cane.var)
プロットしてみると以下のようです。ｘ軸が品種（var)で、縦軸が病気に罹っている芽の比率っす。ただし、比率はＡＢＣＤの４ブロックの平均値でしかないので、ブロックによる傾向などあっても無視してます。また、上で述べたｘ（植えた本数）も無視してます。
上のグラフみると、「よわよわ」な２種は一目瞭然だけれども、「強い」と「中間」の間はどこに線引きしますか？と言う感じだわな。
なおチャンピオンを決めるために以下のようにしてみました。
varRank <- cane.var[order(cane.var$ratio),]
頭の方だけダンプします。この計算法だと、チャンプは２５番さんね。知らんけど。
BOXプロットで足元確認
上ではテキトーに処理してみましたが、ブロックＡＢＣＤで何か傾向があるのか否かも気になります。とりあえず以下のようにしてＢＯＸプロットを描いてみます。
p0 <- ggplot(data=cane.work, mapping=aes(block, y=ratio, fill=block)) + geom_boxplot()
p1 <- p0 + ylab("ratio") + ggtitle("Diseased ratio by block")
p1
描いたものが以下に。
本来だったら、ここで「統計的に有意な差があるのかないのか」など吟味せにゃならんのでしょうが、素人は、まあいいんじゃね、と多少の差は踏みつぶします。
つづいて肝心の品種（var)の平均値ではなくバラツキを同じくBOXプロットで眺めてみます。
p1 <- ggplot(data=cane.work, mapping=aes(var, y=ratio, fill=var)) + geom_boxplot()
p2 <- p1 + ylab("ratio") + ggtitle("Diseased ratio by var")
p2
描いたものが以下に。
これでみると
-
- 強い奴ら、下の方にギュッと詰まっている
- 中間派、下から真ん中くらいまでバラツク
- 弱弱、全滅（1.00)に近いところにいる
ってな感じがするよの。そんなんで良いのか？