統計のコツのこつ(60)

このブログは「すぐに役立つ統計のコツ」(オーム社)をもとに書いています。

 

 

今回は本書(114p)の判別分析についてご紹介します。
本書の「例題19」では腎障害患者を"慢性糸球体腎炎患者"と"尿細管障害患者"の２群に判別する方法として、ロジスティック回帰分析による判別分析をご紹介しています。
すなわち、独立変数の効果に注目した判別方法と言えます。しかし、
今回ここでは、データ解析環境「R」の関数式を用いた線形判別分析をご紹介します。例題は、

本書(例題19)のデータに「U.BGM(尿中β2-Microgloburin)」を加えた表1のデータを使って見ましょう。

 

表1　腎障害患者のExcelデータ

 

表1 の Excelデータを表２の様に縦長のフォームに加工し、コピー＆ペーストで「R」に取り込みます。

 

表２　縦長にしたデータのフォーム

 

（fDiag の「０：慢性糸球体腎炎患者、１：尿細管障害患者」としています）

 

「R」プログラムの実行

***
dat<- read.delim("clipboard", header=T)　# データの取り込み
head(dat)
s.dat<- scale(dat[, 2:4])　# データの標準化との違いをみるため
dat1<- cbind(dat, s.dat)   # 原データと標準化データの結合
head(dat1) # 結合データの確認

y1<- dat1[, 1]　# 目的変数(fDiag)
x1<- dat1[, 2]  # 説明変数(U.Pro)
x2<- dat1[, 3]  # 説明変数(S.BMG)
x3<- dat1[, 4]  # 説明変数(U.BMG)

library(MASS) # パッケージを事前にインストールしておく
z1<- lda(y1~ x1+ x2+ x3, CV=F)　# lda関数による判別分析(CV=Fは交差確認なし)

z1

 

出力結果：

Call:
lda(y1 ~ x1 + x2 + x3, CV = F)

Prior probabilities of groups:
      0       1
0.53125 0.46875

Group means:
        x1       x2        x3
0 248.5882 2.857059  1054.235
1 120.0000 8.098667 35288.667

 

Coefficients of linear discriminants:
             LD1
x1 -0.0030586273
x2 -0.0146857850
x3  0.0001160874
***

 

よって判別式は、
z1=-0.0030586273*U.Pro- 0.014685785*S.BMG+ 0.0001160874*U.BMG- c0

 

となります。ここで、c0 は次により求めることが出来ます。

 

c0<- apply(z1$means%*%z1$scaling,2,mean)
c0=1.465343

 

上記で求めた判別得点と下記の判別得点は若干異なる様です・・。

実際には、次により判別結果を求めることが出来ます。

***
predict(z1)$x　　　# 判別得点
predict(z1)$class　# 判別されたグループ
round(predict(z1)$posterior,3)　# 事後確率
lda.tab<- table(dat[, 1], predict(z1)$class)　
lda.tab 　　　　　 # 判別結果
Correct<- sum(lda.tab[row(lda.tab)==col(lda.tab)])/sum(lda.tab)*100 # 判別率(%)
Correct
100-Correct        # 誤班別率(%)
***

 

ここで、未知のデータでの求め方の一例です。

***
 # 未知データの判別

newdata<- rbind(c(320,1.36,2200),c(120,12.7,48000))
dimnames(newdata)<- list(NULL, c("U.Pro","S.BMG","U.BMG"))
newdata<- data.frame(newdata)
predict(z1, newdata=newdata)

***

 

なお、
標準化したデータを使うなら次のようにして下さい。

***
y12<- dat1[, 1] # 目的変数(fDiag)
x12<- dat1[, 5] # 標準化した説明変数(U.Pro)
x22<- dat1[, 6] # 標準化した説明変数(S.BMG)
x32<- dat1[, 7] # 標準化した説明変数(U.BMG)

z12<- lda(y12~ x12+ x22+ x32, CV=F)
z12

***

 

また、非線形判別分析は次の関数があります。

***
zq<- qda(y1~ x1+ x2+ x3, data=dat)
zq
***

 

ここで、
z1(z12)の CV を CV=T とした場合は"交差確認"の結果を得ます。
これに付いては、
「統計学入門、杉本典夫 先生」(第１章８)を参考にして下さい。

http://www.snap-tck.com/room04/c01/stat/stat01/stat0108.html

 

次回は、SPSS（正準判別分析）でやって見ましょう。

 

情報統計研究所はここから！

 

 

　

 

統計のコツのこつ(59)

昨年末の突然の"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) 

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 # 全体（用量）の平方和
***

 

用量反応直線の回帰式は次の通りです。

***
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

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
***

 

すなわち、
用量反応直線：p=0.47273 + 0.255*log10(dose)
寄与率：r2=0.9668(96.7％)

となります。

 

ここで、切片と傾きは次の様になります。

***
B<- fit$coefficients[1]  # 切片
A<- fit$coefficients[2] # 傾き
***

 

そして、
回帰の直線性は、次の「prop.trend.test」で求めることが出来ます。

 

***
output<- prop.trend.test(res1, n)  # デフォルトのスコアを使う

output    # prop.trend.testの結果
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
***

 

ここで、
「表13.4.3」のズレ(異質性)は次により求めることが出来ます。

 

***
chiLOF<- chi$statistic - output$statistic
chiLOF                      # ズレ(異質性)のカイ二乗値
1-pchisq(chiLOF, 3)  # p値

 

出力結果：

> chiLOF<- chi$statistic - output$statistic
> chiLOF
X-squared

1.79139      # ズレ(異質性)のカイ二乗値

 

> 1-pchisq(chiLOF, 3)
X-squared
0.6168103 # p値
***

 

ここでの「prop.trend.test」は外的基準は省略し、「1」から始まるスコア(整数値)となっています。
また、
例題にもとずいていることにご留意して下さい。そして、
コクラン・アーミテージの傾向分析によって求めた用量反応直線によるD50の推定は「統計学入門(第13章４)」から筆算で容易に求められます。

 

情報統計研究所はここから！