統計ブログはじめました！

各専門分野の統計技術、方法、テクニックなどを気ままに分かり易く例題をもとに解説します。

統計のコツのこつ(48)

2017-06-28 18:58:19 | 日記・エッセイ・コラム
このブログは「すぐに役立つ統計のコツ」(オーム社)をもとに書いています。
杉本典夫 先生の
「統計学入門　第5章（5.5 各種手法の相互関係）」(下記URL)も参考にしています。
http://www.snap-tck.com/room04/c01/stat/stat05/stat0505.html
 
ところで、
上記URLに、「主軸回帰(MA回帰：Major axis regression」についての解説がありましたので、ひとつ下記の例題でやって見ましょう。例題は、
血清CPK活性値についてです。ちなみに、
CPK活性値は心筋梗塞の発症初期に異常に高い値を示しますので、心筋梗塞の診断によく用いられます。しかし、
採血し遠心分離した血清を窓際などの日射下に不用意に放置しておくとCPK値が減少しますので、採血後は速やかに冷蔵庫に保管し、なるべく、速やかに測定すべきでしょう。
 
さて、
例題は次の条件で測定した血清CPK活性値です。
・A1h：冷蔵庫に1時間放置
・B1h：室温28℃の日射下に1時間放置
・A2h：冷蔵庫に2時間放置
・B2h：室温28℃の日射下に1時間放置
 
図1　血清CPK活性値の比較
 
 
 
上記の２群間の相関関係は、単純に、相関回帰分析によって見ることが出来ます(すぐに役立つ統計のコツ、第６章 参照)。
ここでは、下記について既存の統計ソフトによる方法をご紹介します。
 
・単純な線形回帰の方法（OLS：Ordinary Least Squares）：
　　ⅹに誤差がなくＹに誤差がある。
・Deming の方法(重み付き線形回帰分析)：
　　ⅩとＹの分散比で補正。
・主軸回帰の方法(主軸回帰(MA回帰：Major axis regression)：
　　Ⅹに対するＹ、Ｙに対するⅩの傾きの幾何平均による回帰分析。
・Passing＆Bablokの方法(Passing–Bablok regression)：
　　すべての傾きの中央値による回帰分析(ノンパラメトリック法)。
 
ご紹介するソフトは下記URLからダウンロード出来ます。
 
ダウンロード先のご紹介（2017.6.24 現在確認）：
「ValidationｰSupport処理プログラム」
 
http://www.jscc-jp.gr.jp/?page_id=1145
 
図２　日本臨床化学会
 
 
図３　ダウンロード
 

 
図３から、
「Validation-Support-V318-改良ソフト3 (zipファイル)」をダウンロードし展開(解凍)して本ソフトを立ち上げて下さい。
そして、
Sheet名「相関分析」を開き→「解説」(１行Ｃ列)にマウスを重ねて説明に目を通して下さい。
次に、
図4の様に、例題(CPK)のデータ(A1h-B1h、A2h-B2h)のみを入力しますと、即座に結果が表示されます。
 
図４　データの入力(コピー＆ペースト)
 

 
なお、
出力結果の編集は出来ない様ですので、項目名をそのまま使用して下さい。
・試薬１→A1h
・試薬２→B1h
・Ⅹ2 　→A2h
・Ｙ2   →B2h

●「ブートストラップ」をクリックし「Demingと標準主軸回帰の信頼限界」を求めて下さい。
● 誤差分散比は「１」を入力して下さい。
● 線形回帰(古典的回帰式)の「ｙ→X」の係数(b)と切片(a)は、次により計算されている様です。
　 勾配(b)＝Syy／Sxy＝101474.6／126861.4≒0.8
    切片(a)＝Ymean－b*Xmean＝81.65－0.8*96.05≒4.8
 
次回に、
データ解析環境「R」での方法とあわせて、分析結果を検証したいと思います。
 
情報統計研究所はここから！
 

統計のコツのこつ(47)

2017-06-17 19:13:46 | 日記・エッセイ・コラム
前回の「R]による検量線の逆推定では「R」パッケージ(investr)の"invest()"による方法をご紹介しました。
そして、「直線回帰からの逆推定(95%信頼限界)」(前回の図１)での信頼限界(inversion interval)は、「192.88～220.67」となっています。
これは、FCA(Filler's Confidence Interval,1954)の式を用いていない様です。Fillerに付いては杉本典夫先生の下記webページに詳しく紹介されていますので、ご参照下さい。
*****
第13章　用量反応解析
http://www.snap-tck.com/room04/c01/stat/stat13/stat1301.html#note02
*****
 
ここでは、
Fiellerの式による"FCI"の結果と一致する「R」での１方法をご紹介しておきます(ただし、確証はありませんのでご注意ください)。
 
「R」プログラム
*********************************************
X<- c(50, 100, 200, 300, 400)
Y<- c(0.09, 0.15, 0.29, 0.42, 0.52)
dat<- data.frame(x=X, y=Y)
dat
library(investr)
mod <- lm(y ~ x , data = dat)
(res <- calibrate(mod, y0 = 0.29　, interval = "inversion", level = 0.95))
plotFit(mod, interval = "prediction", level = 0.95, shade = TRUE, col.pred = "red")
abline(h = 0.29, v = c(res$lower, res$estimate, res$upper), lty = 2)

出力結果
> res$lower
[1] 172.7452　　
> res$estimate
[1] 206.8093
> res$upper
[1] 240.7984
*********************************************
 
それでは、今回の本題です。
前回・前前回と検量線(Calibration)についてご紹介してきました(すぐに役立つ統計のコツ, 83ページ参照）。
前回の例題の場合ですが、
検量線として、2次曲線からの逆推定はあまり良くありませんでした。
そこで、
図１の様な線形補完(linear interpolation)を適用して見ようと思います。
 
ここで言う線形補間は"折れ線グラフの線形補間"であり、線形多項式(一次式)による回帰分析を検量線として利用したもので、1つの方法としてご理解下さい。
 
図１　線形補間の検量線
 

　図２　MS-Excel による区間ごとの計算
 
 
図３　未知濃度(Ⅹ0)の計算
 

図４　エクセル関数式(図２、図3の計算式)
 

各区間を２点の直線で結び、その直線上にある吸光度(Y0)から濃度(X0)を推測する訳です。
検量線は直線であれば逆推定も容易ですが、RIA や EIA など曲線の場合も多々あります。
統計的方法を理解して、最も適切な方法を選ぶことが肝要かと思います。
自動分析機器の時代でもcalibration の方法を確認しておきたいものです。

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

　
 
 
 
 

統計のコツのこつ(46)

2017-06-15 10:38:21 | 日記・エッセイ・コラム
前回はMS-Excelでの検量線についてご紹介しました。
今回は、
同じ例題を使ってデータ解析環境「R」でやって見ましょう。
初めに、
package"investr" をインストールしておいて下さい。では、
前回の図１のデータ(直線回帰)を使って見ましょう。
 
「R」プログラム
***
# 標準物質の濃度(Ⅹ)と吸光度(Ｙ)
 
XY
datdat
 
fit1summary(fit1)
library(investr)
INV.fit1round(INV.fit1$estimate,1)
plotFit(fit1, interval = "confidence")
abline(h = 0.29, v = c(INV.fit1$lower, INV.fit1$estimate, INV.fit1$upper), lty = 2, col = "red")

「R」の実行結果
> dat
    x    y
1  50 0.09
2 100 0.15
3 200 0.29
4 300 0.42
5 400 0.52
 
> summary(fit1)
Call:
lm(formula = y ~ x, data = dat)
Residuals:
......1..............2...............3.............4.............5
-0.003415 -0.006098  0.008537  0.013171 -0.012195
 
Coefficients:
..................Estimate....Std. Error...t value...Pr(>|t|)   
(Intercept)..3.073e-02..1.045e-02....2.941...0.0604 . 
x................1.254e-03..4.248e-05...29.512..8.54e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01216 on 3 degrees of freedom
Multiple R-squared:  0.9966,    Adjusted R-squared:  0.9954
F-statistic:   871 on 1 and 3 DF,  p-value: 8.544e-05

> round(INV.fit1$estimate,1)　＃ 吸光度(y0=0.29)の逆推定値
[1] 206.8
 
図１ 直線回帰からの逆推定(95%信頼限界)　
 
***

次に、
前回(図３)のデータ(2次方程式の当てはめ)を使って見ましょう。
 
「R」プログラム
***
# 標準物質の濃度(Ⅹ)と吸光度(Ｙ)
XY 
datdat
 
fit2summary(fit2)
 
library(investr)
INV.fit2

round(INV.fit2$estimate,1)　＃ 吸光度(y0)の逆推定値
 
plotFit(fit2, interval = "confidence")
abline(h = 0.835, v = c(INV.fit2$lower, INV.fit2$estimate, INV.fit2$upper), lty = 2, col = "red")
 
「R」の実行結果
> dat
     x     y
1   23 0.056
2   46 0.126
3   92 0.224
4  180 0.419
5  370 0.601
6  740 0.835
7 1500 1.177

> summary(fit2)
Formula: y ~ a * x^2 + b * x + c
 
Parameters:
......Estimate...Std. Error...t value...Pr(>|t|)  
a..-4.989e-07..1.285e-07..-3.883....0.01780 *
b...1.469e-03..1.980e-04....7.423....0.00176 **
c....8.395e-02..4.207e-02....1.995....0.11673  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.06677 on 4 degrees of freedom
Number of iterations to convergence: 1
Achieved convergence tolerance: 5.894e-09
 
> round(INV.fit2$estimate,1) ＃ 吸光度(y0=0.835)の逆推定値
[1] 658.3

図１ 2次曲線回帰からの逆推定(95%信頼限界)　
 
 
***
 
検量線として、2次方程式からの逆推定はあまり良くないようです。
この様な、曲線になる検量線では片対数変換とか両対数変換を試みることもあります。
(すぐに役立つ統計のコツ、85ページ 参照)
 
次回は、線形補間による方法をご紹介したいと思います。
 
情報統計研究所はここから！
 
 

統計のコツのこつ(45)

2017-06-10 11:58:04 | 日記・エッセイ・コラム

だいぶ間が開いたけど、「すぐに役立つ統計のコツ」(第6章)では"相関と回帰"の手法を紹介しています。
今回は、
医学・生物学などの実験や臨床検査の現場でよく用いられる検量線について、ご紹介しましょう。

本書の83ページの"４. 臨床検査で大切な検量線"を見て下さい。通常、
検量線は校正直線(曲線)と呼ばれ、標準物質濃度(x)に対応する分光光度計などの吸光度(y)を、グラフ用紙の横軸(x)と濃度を縦軸(y)にプロットし、測定物質の未知の吸光度(y0)を読み取り、未知物質の濃度(x0)を求めます。こでは昔のことです。
現在では、
ｘとｙの直線(曲線)回帰式、例えば、ｙ＝aⅹ + b(1次方程式) や ｙ＝aⅹ^2＋bⅹ＋ｃ(2次方程式) を使い、吸光度(ｙ0)を知って濃度(ｘ0)をパソコンなどで推定します。これを回帰式からの逆推定と言います。
これは、
回帰式が「ⅹを基準にｙを回帰するので、回帰誤差はｙの方にである」からです。
それでは、
MS-Excelで実際に例題を使ってやって見ましょう。
血清中の総コレステロール値を測定する場合、
標準物質(既知濃度のもの)を倍々希釈などして、図１の様な濃度(ⅹ)と吸光度(ｙ)が得られたとします(図１)。
 
図１ コレステロールの濃度(x)と吸光度(y)
 
図１の緑色セル部分をクリックし、「挿入→グラフ→散布図→散布図」で図２の散布図を得ることが出来ます。
 
図２ 濃度(x)と吸光度(y)の散布図と近似直線
 
 
そして、
図２のプロットのどれかを、「右クリック→近似曲線の追加→近似曲線のオプション：
[直線近似◎、グラフに数式を表示する◎、グラフにR-2乗値を表示する◎]→[×](終了)]
とすれば、
図２の様な近似直線が得られます。盲検(ブランク)をとれば原点近くを通る直線となります。
ここで、
直線回帰式は 吸光度(ｙ)＝0.0013×濃度(ⅹ)＋0.0307 となりますので、その解はⅹ0＝(ｙ0－0.0307)／0.0014 となります。
よって、
未知吸光度ｙ0＝0.42 の濃度はⅹ0＝299.5 と逆推定出来ます。
 
検量線は曲線である場合も多々あります。
本書の85ページの上段にある「濃度と吸光度」を使って見ましょう。
 
図３ 濃度(x)と吸光度(y)のデータ
 

 
図３の緑色セル部分をクリックし、同様の方法で図４の様な散布図を作ります。
 
図４ 濃度(x)と吸光度(y)の散布図と近似曲線
 

 
図４の回帰式は、
ｙ＝aⅩ^2＋bⅩ＋c＝(－5E－07)・Ⅹ^2＋0.0015・Ⅹ＋0.084
 
ですので、
この２次方程式の解から、未知吸光度(ｙ0)に対する濃度(ⅹ0)の逆推定は、
ⅹ0＝(－b±√(b^2－4a(c-y0))／(2a)
 
です。よって、
未知吸光度 ｙ0＝0.024 のとき、
 
Ⅹ0＝(-0.0015+√(0.0015^2-4*-0.0000005*(0.084-0.224)))/(2*-0.0000005)＝96.4
 
となります。
 
検量線としては当てはめがあまり良くありませんので、別の関数へのあてはめを検討する必要がありそうです。
 
対数変換による当てはめは「すぐに役立つ統計のコツ」(85ページ)を見て下さい。
本書の記載(校正)ミスの訂正は"情報統計研究所"(TOPページ)の「新刊書紹介(正誤表あり)」にあります。
***
訂正：
本書(85ページ)の常用対数変換濃度は次の次の通り訂正します。
1.362, 1.663, 1.964, 2.255, 2.568, 2.869, 3.176
***
 
なお、
「R」の環境があれば、もっと簡単に詳しい情報を得ることが出来ます。それは、次回にご紹介したいと思います。
 
次回は、
「R」による２次曲線からの逆推定です。
 
情報統計研究所はここから！