統計ソフトRの初心者です。

質問

Ｒを使って解くんですけどわからない問題があるので教えてください。

次のデータは害虫蛾を駆除するときのピレスロイド系殺虫剤に対する語毒性の抵抗力を調べたものである。雄雌格２０匹を集め、様々な投薬量の箱に２、３日入れて死亡した蛾の数を調べたところ次のようになった。

　　　　　　　　　薬物量
性別　１　２　４　８　１６　３２
オス　１　４　９　１３　１８　２０
メス　０　２　６　１０　１２　１６

このデータを使ってロジスティックモデルを当てはめたい。ただし、投薬量は２のべき乗なので、log2（投薬量）を使う。データ入力は次のようにする。

#オスの場合
> x=as.numeric(gl(6,20,label=c(1:6))) #投薬量
> numdead=c(1,4,9,13,18,20) #死亡数
> y=vector(length=0) #２値のデータの作成
> for(i in 1:6)
+ {
+ y=append(y,c(rep(1,numdead[i]),rep(0,20-numdead[i])))
+ }

(a)オス、メスの死亡率のロジスティックモデルを当てはめよ
(b)推定した曲線と、データを合わせてプロットし、特徴を述べよ
(c)オス、メスのオッズ比の９５％信頼区間を構成し、その値を解釈せよ

お願いします。

回答No.1

backs · ベストアンサー

まず　#オスの場合　として示されているコードを実行させても何も起きないですよ？xには投薬量とコメントがつけられていますが、これはどう見ても投薬量のデータベクトルではないようですが。。。死亡数というのは分かりますが、2値データの作成というコメントとyに代入されているのは何なのか？最後の繰り返し処理も理解できません（私の実力不足で）。

よく分かりませんが、致死量の推定と同じような状況ですかね。Rに最初から実装されているnls()を使うのなら次のように。

dose <- c(1.691, 1.724, 1.755, 1.784, 1.811, 1.837, 1.861, 1.884) #投与量
ind.num <- c(59, 60, 62, 56, 63, 59, 62, 60) #個体数
dead.num <- c(4, 10, 19, 31, 52, 53, 60, 60) #死亡数
dead.rate <- dead.num / ind.num #死亡率
#データフレームの作成
myData <- data.frame(DOSE = dose, INDIVIDUAL = ind.num, DEAD = dead.num, RATE = dead.rate)
attach(myData) #データフレームにパスを通す
plot(DOSE, RATE, pch=16, col="green") #死亡率をグラフ化
result <- nls( #nlsを使ってモデル解析
RATE ~ gamma / (1 + exp(alpha + beta*DOSE)), #ロジスティックモデル式
start = c(gamma = 1, alpha = 64.77219, beta = -36.49183) #初期値の設定
)
x <- seq(1.6, 1.9, 0.01) #ロジスティック曲線を描くためのx軸の値
points(x, 1.027 / (1 + exp(61.100 - 34.383*x)), type="l") #ロジスティック曲線を描く

LD50を求めるのならhttp://aoki2.si.gunma-u.ac.jp/R/ed50.htmlなど参考になるでしょう。シンプレックス法を用いてやるならhttp://aoki2.si.gunma-u.ac.jp/R/simplex.htmlが参考になります。いずれも群馬大学の青木先生のwebサイトにあるものです（Rを使うのならここはぜひ見ておくべきです）。

お礼

親切、丁寧に教えていただきありがとうございました。参考にします。