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

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

統計のコツのこつ(47)

2017-06-17 19:13:46 | 日記・エッセイ・コラム
前回の「R]による検量線の逆推定では「R」パッケージ(investr)の"invest()"による方法をご紹介しました。
そして、「直線回帰からの逆推定(95%信頼限界)」(前回の図1)での信頼限界(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」での1方法をご紹介しておきます(ただし、確証はありませんのでご注意ください)。
 
「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次曲線からの逆推定はあまり良くありませんでした。
そこで、
図1の様な線形補完(linear interpolation)を適用して見ようと思います。
 
ここで言う線形補間は"折れ線グラフの線形補間"であり、線形多項式(一次式)による回帰分析を検量線として利用したもので、1つの方法としてご理解下さい。
 
図1 線形補間の検量線
 

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

図4 エクセル関数式(図2、図3の計算式)
 

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

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

 
 
 
 
 

この記事をはてなブックマークに追加

統計のコツのこつ(46)

2017-06-15 10:38:21 | 日記・エッセイ・コラム
前回はMS-Excelでの検量線についてご紹介しました。
今回は、
同じ例題を使ってデータ解析環境「R」でやって見ましょう。
初めに、
package"investr" をインストールしておいて下さい。では、
前回の図1のデータ(直線回帰)を使って見ましょう。
 
「R」プログラム
***
# 標準物質の濃度(Ⅹ)と吸光度(Y)
 
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
 
図1 直線回帰からの逆推定(95%信頼限界) 
 
***

次に、
前回(図3)のデータ(2次方程式の当てはめ)を使って見ましょう。
 
「R」プログラム
***
# 標準物質の濃度(Ⅹ)と吸光度(Y)
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

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

この記事をはてなブックマークに追加

統計のコツのこつ(45)

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

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

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

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

 
図4の回帰式は、
y=aⅩ^2+bⅩ+c=(-5E-07)・Ⅹ^2+0.0015・Ⅹ+0.084
 
ですので、
この2次方程式の解から、未知吸光度(y0)に対する濃度(ⅹ0)の逆推定は、
ⅹ0=(-b±√(b^2-4a(c-y0))/(2a)
 
です。よって、
未知吸光度 y0=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」による2次曲線からの逆推定です。
 
情報統計研究所はここから!
 
 
 
 
 
 

この記事をはてなブックマークに追加

統計のコツのこつ(44)

2017-05-19 10:59:08 | 日記・エッセイ・コラム
杉本典夫 先生から大変参考になるコメントを頂きましたのでご紹介させて頂きます。
 
コメント要旨:
2×2分割表には色々な手法を適用できるため、統計学解説書等ではそれらの手法が体系的に整理されていない傾向があると思います。例えば2×2分割表の場合は次のようになります。
2種類の2分類化された計量尺度項目の相関性を検討したい時。
 
今回の話題であるピポットテーブルのように、
2種類の計量尺度項目を境界値で2分類にした時、2つの項目の間にどの程度の相関性があるかを検討したい場合です。
 
・研究デザイン:横断的研究
・評価指標:四分点相関係数=φ(ファイ)係数
・検定手法:四分点相関係数の検定=Mantel-Haenszelの検定
・推定手法:Fisherのz変換を利用した相関係数の区間推定
 
χ2乗検定は本来は連関係数の検定ですが、これに連続修正(Yatesの修正)を施すと、出現率の正規近似検定と同じものになります。そのため検定対象の評価指標の種類によって、連続修正を施すか施さないかを使い分けると良いと思います。
 
SASやSPSSやJMPやRでは、通常は出現率の区間推定に連続修正を施しません。
 
しかし
連続修正は超幾何分布を正規分布でより正確に近似するためのものですから、区間推定でも施す必要があります。
出現率が低い時は、2つの項目の関係が直線ではなく指数関係に近くなります。
そのため出現率の差の代わりに、リスク比とその検定と推定を用いる時があります。
リスク比を対数変換すると、対数変換した出現率の差つまり対数リスク差になります。
そのためリスク比の検定と推定は、対数変換した出現率の差を正規近似して検定と推定を行っていることに相当します。
 
 
 
 

この記事をはてなブックマークに追加

統計のコツのこつ(43)

2017-05-10 13:29:08 | 日記・エッセイ・コラム
前回の例題「Excel_Sample(1).xlsx」について、ピボットテーブルを「R」で作って見ましょう。
インターラクティブ(ブラウザ上で操作)でピボットテーブルを作成できる便利なツールです(筆者おススメ)。
それでは、前回と同様に「情報統計研究所」のホームページ(Top)から、「Excel_Sample(1).xlsx」をダウンロードして使用しましょう。
 
「R」の実行
***
dat<- read.delim("clipboard", header=T)
 # 前回と同様にして取り込みます
head(dat) # データの確認
 
 # 事前に "rpivotTable" をインストールしておいて下さい。
library(rpivotTable) 
 
  # Age(年齢)とTG(中性脂肪)で試してみましょう。
rpivotTable(dat, row = "Gender", col = "TG") # たったこの一行でOKです。
***
 
実行結果:
 
図1:Age と TG のピボットテーブル
 
 図1の[Count]を変えてみましょう。
 ・[Count as Fraction of Total]      →全体の%
 ・[Count as Fraction of Rows ]     →行  の%
 ・[Count as Fraction of Columuns]→列  の%
 
また、
図1の[Table]を変えてみましょう。
 ・[Horizontal Bar Chart]      :ヨコ棒
   ・[Horizontal Stached Bar Cart]    :積み重ねヨコ棒
 ・[Bar Chart]            :タテ棒
 ・[Stacked Bar Chart]       :積み重ねタテ棒
 
その他にも、
 [Line Chart]や[Area Chart]などがあります。
 
また、
変数も自由に変えられます。
「Age」を「Gender」に変えてみましょう。
 
灰色エリアの「Gender」をクリック保持し、青色エリアに移動させ「Age」を灰色エリアに戻せば
図2のテーブルが得られます。
 
図2:Gender と TG のピボットテーブル
 
如何でしょう・・・、Excelのピボットテーブルより使いやすいかも知れませんね・・・!?
 
情報統計研究所はここから!
 
 
 
 
 
 

この記事をはてなブックマークに追加