昨年末の突然の"Aufnehmen(admission)" から "Entlassen" して間がないのですが、新年早々から依頼の統計分析に追われていました。今年も健康に留意しながら投稿したいと思っていますので、よろしくご指導ご鞭撻のほどお願いします。
さて、前回(58)を思い出しながら、「コックラン・アミテージ（Cochran Amitage）」による用量反応直線について、「Ｒ」での一例をご紹介しておきます。
どうぞ・・、、
以下の「統計学入門(第13章４)」(杉本典夫先生)を見て下さい。
http://www.snap-tck.com/room04/c01/stat/stat13/stat1304.html
「表13.4.3 コクラン・アーミテージの傾向分析による用量反応解析の分散分析表」を「Ｒ」でやってみましょう。
上記HPの筆算過程を参考にやってみて下さい。
データは前回と同じで「Ｒ」では次のようになります。
***
res1<- c(0,2,16,15,19)
res0<- c(20,18,14,5,1)
dose<- c(0.01,0.1,1,10,100)
res1<- c(0,2,16,15,19)
res0<- c(20,18,14,5,1)
dose<- c(0.01,0.1,1,10,100)
tbl<- matrix(c(0,2,16,15,19,20,18,14,5,1), nc=2)
tbl
chi<- chisq.test(tbl) # 全体の平方和
chi # カイ二乗検定の結果
chi[1] # 全体(用量)のカイ二乗値
出力結果：
Pearson's Chi-squared test
data: tbl
X-squared = 53.967, df = 4, p-value = 5.348e-11 # 全体（用量）の平方和
***
tbl
chi<- chisq.test(tbl) # 全体の平方和
chi # カイ二乗検定の結果
chi[1] # 全体(用量)のカイ二乗値
出力結果：
Pearson's Chi-squared test
data: tbl
X-squared = 53.967, df = 4, p-value = 5.348e-11 # 全体（用量）の平方和
***
用量反応直線の回帰式は次の通りです。
***
n<- res1 + res0
r<- res1/n
p<- sum(res1)/sum(n)
w<- n/p/(1 - p)
fit<- lm(r~ log10(dose), weights = w, data=dat)
summary(fit)
n<- res1 + res0
r<- res1/n
p<- sum(res1)/sum(n)
w<- n/p/(1 - p)
fit<- lm(r~ log10(dose), weights = w, data=dat)
summary(fit)
出力結果：
Coefficients:
.............Estimate...Std. Error...t value...Pr(>|t|)
(Intercept)...0.47273....0.03678......12.851...0.00102 **
log10(dose)...0.25500....0.02728.......9.348...0.00259 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
.............Estimate...Std. Error...t value...Pr(>|t|)
(Intercept)...0.47273....0.03678......12.851...0.00102 **
log10(dose)...0.25500....0.02728.......9.348...0.00259 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.7727 on 3 degrees of freedom
Multiple R-squared: 0.9668, Adjusted R-squared: 0.9557
F-statistic: 87.38 on 1 and 3 DF, p-value: 0.002593
***
Multiple R-squared: 0.9668, Adjusted R-squared: 0.9557
F-statistic: 87.38 on 1 and 3 DF, p-value: 0.002593
***
すなわち、
用量反応直線：p=0.47273 + 0.255*log10(dose)
寄与率：r2=0.9668(96.7％)
用量反応直線：p=0.47273 + 0.255*log10(dose)
寄与率：r2=0.9668(96.7％)
となります。
ここで、切片と傾きは次の様になります。
***
B<- fit$coefficients[1] # 切片
A<- fit$coefficients[2] # 傾き
***
B<- fit$coefficients[1] # 切片
A<- fit$coefficients[2] # 傾き
***
そして、
回帰の直線性は、次の「prop.trend.test」で求めることが出来ます。
回帰の直線性は、次の「prop.trend.test」で求めることが出来ます。
***
output<- prop.trend.test(res1, n) # デフォルトのスコアを使う
output<- prop.trend.test(res1, n) # デフォルトのスコアを使う
output # prop.trend.testの結果
output[1] # カイ二乗値
output$statistic # カイ二乗値の取出し
output[1] # カイ二乗値
output$statistic # カイ二乗値の取出し
出力結果：
Chi-squared Test for Trend in Proportions
data: res1 out of n ,
using scores: 1 2 3 4 5
X-squared = 52.175, df = 1, p-value = 5.076e-13
***
using scores: 1 2 3 4 5
X-squared = 52.175, df = 1, p-value = 5.076e-13
***
ここで、
「表13.4.3」のズレ(異質性)は次により求めることが出来ます。
「表13.4.3」のズレ(異質性)は次により求めることが出来ます。
***
chiLOF<- chi$statistic - output$statistic
chiLOF # ズレ(異質性)のカイ二乗値
1-pchisq(chiLOF, 3) # p値
chiLOF<- chi$statistic - output$statistic
chiLOF # ズレ(異質性)のカイ二乗値
1-pchisq(chiLOF, 3) # p値
出力結果：
> chiLOF<- chi$statistic - output$statistic
> chiLOF
X-squared
> chiLOF
X-squared
1.79139 # ズレ(異質性)のカイ二乗値
> 1-pchisq(chiLOF, 3)
X-squared
0.6168103 # p値
***
X-squared
0.6168103 # p値
***
ここでの「prop.trend.test」は外的基準は省略し、「1」から始まるスコア(整数値)となっています。
また、
例題にもとずいていることにご留意して下さい。そして、
コクラン・アーミテージの傾向分析によって求めた用量反応直線によるD50の推定は「統計学入門(第13章４)」から筆算で容易に求められます。
また、
例題にもとずいていることにご留意して下さい。そして、
コクラン・アーミテージの傾向分析によって求めた用量反応直線によるD50の推定は「統計学入門(第13章４)」から筆算で容易に求められます。
情報統計研究所はここから！
