成蹊大学理工学研究報告
J. Fac. Sci. Tech., Seikei Univ.
Vol.43 No.1 (2006) pp.17-21

定量限界，検出限界を持つデータによる回帰パラメータの推定
阿部貴行*1，岩崎学*2

Estimation of Regression Parameters from Data with Detection and Quantitation Limits
Takayuki ABE *1

and Manabu IWASAKI *2

ABSTRACT：A regression model with detection limit (DL) d and quantitation limit (QL) q is examined.
Specifically, our goal is to estimate the regression parameters in a linear regression model Y = α + βx + ε,under the condition that the actual values of Y will be observed only if Y is greater than q. For observations
less than q, only the numbers of observations md and mq fallen in the intervals (0, d) or (d, q), respectively,will be reported. First, we show that a conventional method that gives the value q for the mq observationsin (d, q) and d for the md observations in (0, d) produces inadequate estimates. This type of data is, in thecontext of incomplete data analysis, called NMAR (Not Missing At Random), and hence we have to takethe missingness mechanism into account for the estimation. The maximum likelihood estimationprocedure using the EM algorithm is given. This procedure is statistically sound but numericallyintractable in practice. Therefore a simplified method is given, which is quite simple but can give adequate
estimates. A numerical example is shown to illustrate the usefulness of the procedure proposed.
Keywords：censoring, chemical constituents, EM algorithm, maximum likelihood estimates(Received February 25, 2006)

１．はじめにここにひとつ問題がある。それは，微量な化学物質の計測では，検出限界（detection limit = DL もしくは limit of

環境問題および健康への関心が高まる中，工業製品にdetection = LOD）d および定量限界（quantitation limit =含まれ人体に影響を与える懸念のある化学物質の量を把QL もしくは limit of quantitation = LOQ）q がある点であ握し，客観的な数値データとその影響の大きさを示すこる。検出限界と定量限界は次のように説明される（厚生とが，特に近年，日本を含む世界各国の規制当局側から労働省の分析法バリデーションに関するテキスト各種企業あるいは学術機関などに対して求められている。

www.nihs.go.jp/g/ich/quality/q2a/q2a_jp. html 参照）：

もちろん当局からの要請の有無にかかわらず，その種の(i) 検出限界 d：試料中に存在する分析対象物の検出可

分析は科学的にきわめて重要で，企業や各団体はそれを能な最低の量のことである。ただし，このとき必義務として捉える必要がある。しかし，製品中に含まれずしも定量できる必要はない。

る化学物質の多くはきわめて微量であり，その計測・定(ii) 定量限界 q：適切な精度と真度を伴って定量でき量のためには特殊な計測装置および環境を必要とし，結る，試料中に存在する分析対象物の最低の量のこ果として多くのコストがかかってしまう。そこで，計測とである。定量限界は，試料中に存在する低濃度のが比較的容易な別の物質の値 x から当該化学物質の値物質を定量する場合の分析能パラメータであり，特y を回帰式 y = a + bx などにより予測する試みがなされに，不純物や分解生成物の定量において評価される。ている。

この種の限界がある場合，データは完全な計測値としては得られず，回帰パラメータの推定に工夫が必要となる。

＊1

本論では，次の第 2 節で典型的な実際例および説明の＊2ための人工例を示し，問題の所在を明らかにする。第 3：万有製薬株式会社，情報処理専攻博士学生：情報処理専攻教授 (iwasaki@st.seikei.ac.jp)－17－あてはめたのが図 2 である。このときの回帰直線は節では回帰モデルを定式化し，パラメータ推定のための妥当な方法論（EM アルゴリズム）を与える。この推定y = 5.6485 + 1.2895 x法は統計的には妥当であるものの，実際の計算は反復をであった。図 2 では全 50 個中 d = 3.8 となった観測値が必要とする複雑なものであり，多くの実際家にとってパ4 個，q = 12.8 となったものが 15 個あった。

ラメータの推定値の導出はかなり難しい。そこで，第 4節でそれに代わる簡便法を示し，その妥当性を吟味する。

30.025.0

そして最後の第 5 節で簡単なまとめを行なう。

20.015.0

２．実例と典型的な数値例

10.0
y = 1.5326x + 2.9174
R2 = 0.8005

5.0Counts, et al. (2004) は，紙巻たばこ中の化学物質の予0.00.0測モデルを扱っている。たばこの煙中には 4000 種類以上の化学物質が同定されていて，それらの物質の計量には2.04.06.08.0

10.012.014.016.0

図 1. 人工的に生成した全データと回帰直線専門のラボを必要とし多額の費用もかかる。そこで，計測が容易なたばこの煙中のタールやニコチン量を用いて30.0

種々の化学物質の量を予測しようという試みがなされて25.0いる。Counts, et al. (2004) では主としてタールの量から20.0各化学物質の量を予測する単回帰モデルを扱っている。

15.0

論文中では 48 種類の紙巻たばこが選ばれ，種々の化学物10.0

質の計測値が示されている。表 1 は化学物質とその単位，5.0および DL と QL の値とその個数の一部である。化学物

0.0y = 1.2895x + 5.6485
R2 = 0.7454
0.0

2.0

4.0

6.0

8.0

10.0

12.0

14.0

16.0

質の中にはきわめて微量しか含有されていないため，定図 2. 検出限界，定量限界のあるデータ

量限界，検出限界を超える測定値が得られないものもあった（Selenium，Nickel はその例）。

図 2 は，検出限界，定量限界の値をそのままデータと表 1.
Constituent
Selenium
NickelChromium
Arsenic
Lead
Resorcinol
Crotonaldehyde

して用い，それらが限界値であることを考慮しなかった

DL と QL の例（48 種中）

unitng/cig
ng/cigng/cigng/cigng/cigµg/cigµg/cigDL
2.216.47
5.94
1.13
3.85
0.158
0.98

QL
#(DL)
7.37
48
21.6
46
19.8
37
3.75
812.8
40.52633.29
1ときの計算で，全データでの回帰式と比較すると次の 2

#(QL)
21137
191112つの特徴が読み取れる。
（ア）回帰直線の傾き b が 1.5326 から 1.2895 と小さくなっている。
（イ）決定係数 R2 が 0.8005 から 0.7454 と減少している。
これらは両方とも，限界を考慮しないパラメータ推定では不適切な結果を生み出すことを意味する。そもそもの以下では，人工的な数値例を基に具体的な議論を行な分析の目的は目的変数 y（この場合は鉛を想定）の予測う。人工例ではあるが，Counts, et al. (2004) における鉛であり，あまり大きくないとはいえ（イ）の予測精度の(lead) に類したものである。図 1 は乱数により生成した低下は好ましいものではない。問題は（ア）である。y は全データ 50 個とそれにあてはめた回帰直線である。回人体への有害物質であり，その過小評価は我々の健康に

帰直線の式は対し重大な影響を及ぼしかねない。回帰直線の傾きが小y = 2.9174 + 1.5326 x

さいと，大きな x に対する y の値は本来の予測値よりであった（散布図の作成ならびに回帰直線の導出には小さくなってしまい，悪影響が懸念される y の大きな部

MS EXCEL を用いた。表示は EXCEL のデフォルトの設分での過小評価につながる。図 3 は図 1 の全データと図定によったので必要以上に桁数が多い。以下同様）。それ2 の限界を無視した回帰分析での残差をプロットしたもに対し，検出限界を d = 3.8，定量限界を q = 12.8 と設のである。説明変数 x の大きな部分で，限界を無視した定し，観測値が d 以下のときは d と置き，d より大きときの残差 (R-CEN) の方が全データの残差 (R-ALL)

く q 以下のときは q として散布図を作成し回帰直線をよりも大きく，－18－

y i − yˆ i (CEN ) > y i − yˆ i ( ALL )⇓補題 3.1

Y ~ N(µ, σ2) のとき，yˆ i ( ALL ) > yˆ i (CEN )a* =より，限界を無視すると過小評価につながることが実証a−µ, b* =σb−µσと置くと，以下が成り立つ：

されている。

E[Y | a ≤ Y ≤ b] = µ +ϕ (a*) − ϕ (b*)σ
Φ (b*) − Φ (a*)

8

 a * ϕ (a*) − b * ϕ (b*)
V [Y | a ≤ Y ≤ b] = σ 2 1 +
Φ (b*) − Φ (a*)
，2 ϕ (a*) − ϕ (b*)  − Φ (b*) − Φ ( a*)  

6420-20.0

5.0

10.0

15.0

R-CENR-ALL-4-6

ここで，ϕ(y) および Φ(y) はそれぞれ標準正規分布 N(0,-81) の確率密度関数および累積分布関数である。
図 3. 全データと限界を無視した場合の残差
この補題より次の定理が成り立つ。

３．EM アルゴリズムによるパラメータの最尤推定
定理 3.1 モデル (3.1) および (3.2) の下で，まず，回帰モデルおよび種々の仮定を述べておく。説

E[Y | x; a ≤ Y ≤ b] = α + β x +明変数の値が xi のときの目的変数を表わす確率変数をであり，V[Y | x; a ≤ Y ≤ b] は補題 1 の V[Y | a ≤ Y ≤ b] とYi とし，それらの間に単回帰モデル
Yi = α + β xi + εi (i = 1, . . . , n)，ϕ (a*) − ϕ (b*)σ
Φ (b*) − Φ (a*)(3.1)

同じとなる。

を想定する。ここで，α および β が推定対象となる未知パラメータであり，誤差項 εi には

εi ~ N(0, σ ), R[εi, εj] = 0 (i ≠ j)2

推定すべきパラメータ α，β，σ の対数尤度関数は(3.2)

nn1l = − log(2πσ 2 ) −{ yi − (α + β xi )}222σ 2 i =1∑を仮定する。目的変数は化学物質の量であるので必然的に非負であり，厳密な意味では誤差項の正規性は仮定でである。若干の計算の結果，同時十分統計量はきないが，ここでは近似的なモデルとして正規分布を仮定し，モデルでの Yi が負となる確率は無視できるほどwy =小さいとする。目的変数 Yi の実現値を yi とし，yi の計測には検出限界 d および定量限界 q が存在すると仮定

nnni =1i =1i =1∑ yi , w yy = ∑ yi 2 , wxy = ∑ xi yiする (d < q)。目的変数 y の値が定量限界 q 以下の場合であることが分かる。n 組の全データが観測されたときには，その測定値そのものは分からないものの，それが dの回帰パラメータの最尤推定値は以下であるかもしくは d < y ≤ q であるかは分かるもの

nx=とする。
ここで扱っている問題は，不完全データ解析の言葉で

1xi ,n i =1は，目的変数 y の値に基づく打ち切りであり，欠測メカ

∑

ny=

wy
1
yi =,
n i =1n∑n

ニズムは NMAR (not missing at random) で，無視できな

s xx =1( xi − x ) 2 ,n i =1s yy =

11( yi − y ) 2 = ( w yy − ny 2 )n i =1
ns xy =1
1
( xi − x )( yi − y ) = ( w xy − nx y )n i =1nい (nonignorable) 欠測となる（Little and Rubin (2002) も∑n

しくは岩崎 (2002) を参照）。したがって，欠測メカニズ
ムを考慮した推測を行なわなければならない。
ここでは，回帰パラメータ α と β の最尤推定値を∑およびEM アルゴリズム (Dempster, et al., 1977) により求める。
そのため補題を用意する（Johnson and Kotz (1970) もしくは岩崎 (2002) 参照）。

－19－n∑まず，目的変数の値が得られた m 個のデータのみかとしてb=s xy2, a = y − bx , s = s yy −s xxs xy 2s xx

ら求めた回帰直線は妥当でないことに注意する。欠測が，x には依存するが y には依存しないいわゆる MAR(3.3)(Missing At Random) の場合には欠測データを無視した回帰分析は妥当な結果を与えるが（Little and Rubin (2002)と求められる。
n 組の測定において，検出限界 d 以下の観測値が mdもしくは岩崎 (2002) を参照），ここでの欠測は NMAR個，d より大きく定量限界 q 以下の観測値が mq 個あであるので，欠測は無視可能でない。

り，q 以上で値が観測された観測値が y1, . . . , ym と m

そこで，妥当な結果を与える EM アルゴリズムと不適個あるとき（すなわち n = md + mq + m），最尤推定値の切な結果となった検出限界および定量限界の値をそのま計算のための EM アルゴリズムは以下のようになる（Mま用いた方法を今一度見直し，妥当と思われる方法を見ステップの計算は容易であるが E ステップは，多くの EM出す。検出限界 d より大きく定量限界 q 以下となったアルゴリズムがそうであるように，計算が複雑となる）。

mq 個のデータに対応した目的変数 Yi* の条件付き確率密度関数を fi*(x) とし，累積分布関数を Fi*(x) とすると
(i = 1, . . . , mq)，Yi* の条件付き期待値は・パラメータ推定のための EM アルゴリズム
(0)

(0)

(0)

(1) 初期値 a ，b ，sを適当に定める。
(t)(2) E ステップ：第 t ステップ目でのパラメータ値が a ，b(t)，s(t) のとき，説明変数の値が xi のデータに対し
d i (t ) =d − ( a ( t ) + b ( t ) xi )s(t ), q i (t ) =E[Yi * | d < Yi * ≤ q ] =q1

xf i ( x)dx
Fi (q ) − Fi (d ) ∫d***(4.1)となる。当然 d < E[Yi* | d < Yi* ≤ q] ≤ q であり，この条件q − ( a ( t ) + b ( t ) xi )付き期待値を最大値の q とした図 2 に基づく解析は妥(t )当ではない。また，EM アルゴリスムは (4.1) の条件付sとする。そして，補題 3.1 および定理 3.1 における a* とき期待値を逐次計算するものとなっている。十分統計量b* を，md 個の検出限界以下のデータについてはのひとつ wy（第 3 節参照）における mq 個の寄与分は

a* = – ∞, b* = di(t)mdqi =1

タについてはa* = di(t),mdf * ( x)

∑ E[Yi* | d < Yi* ≤ q] = ∫d x∑ F * (q)i − F * (d ) dx (4.2)とし，mq 個の検出限界より大きく底流限界以下のデー
b* = qi(t)i =1iiである。そしてこの値は md(d + q)/2 に近いことが期待2として yi および yi の期待値を求め，それらにより十される。その理由は 2 つある。第 1 に，fi*(x) が区間 (d, q)

分統計量 wy，wyy，wxy の期待値を計算する。

で左右対称に近ければその期待値は区間の中点 (d +(3) M ステップ：E ステップで求めた十分統計量を用いて，q)/2 に近い。第 2 に，区間 (d, q) で歪度が正のものと負(t+1)

完全データの公式 (3) から a，b(t+1)(t+1)，sを求める。

mdのものとが相殺し(4) 収束判定により，終了するか (2) に戻る。

f * ( x)∑ F * (q)i − F * (d ) は全体として左i =1ii

右対称に近くなる。したがって，(4.2) 式の値が md(d +第 2 節の人工例（図 2 のデータ）に上記の EM アルゴq)/2 となるよう各 yi* の値を定めればよいことになる。

リズムを適用した結果，推定された回帰直線は
y = 1.1511 + 1.7019 x

求めようとする回帰直線は左下がりであり，x が小さであった。傾きの推定値 1.7019 は限界を考慮しない場合くなるほど y も小さくなる傾向にある。そこで，そのよの値 1.2895 よりも大きく，傾きの過小評価が解消されてうに y の値を決めればよいと思われるが，それは，y のいる。

条件付き分散の過小評価につながり，得策とはいえない。
実際，EM アルゴリズムでも，Y の期待値と同時に Y2 の4．簡便法とその性能期待値の推定値を用いて条件付き分散の過小評価を回避している。これらを考え合わせ，最も簡単な方策として，

EM アルゴリズムによる最尤推定は，統計理論的には「定量限界以下のデータの値をすべて (d + q)/2 とす妥当性を持つものであるが，通常の（市販の）ソフトウる」ことが考えられる。

ェアでは対処できず，計算ははなはだ面倒である。そこ

次に，検出限界 d 以下の md 個のデータについて考でここでは，計算が容易でしかも妥当な結果を与える簡える。第 3 節のモデルの定式化の際に述べたが，正規分便法を考える。

布は負の値も取り得るが y は非負であり，厳密には矛盾－20－が生じている。そこで，定量限界以下のデータのアナロの直線回帰式および (3.2) の仮定がそのまま成立しない

ジーとして，便宜的に「検出限界以下のデータの値をす恐れがある。特に yi の分散は，xi の値が 0 に近いときべて d /2 とする」とする（この部分の理論的な正当化は小さく，xi が大きくなるに連れて大きくなることが多い。
このときの対処法として y に適当な変数変換（たとえは困難である）。
上記の簡便法「検出限界 d 以下の y の値は d/2 とし，ば対数変換 log y）を施し，モデル (3.1) および仮定 (3.2)

定量限界 q 以下で検出限界より大きな y の値は (d +の両方が成立するようにする。そして第 4 節で述べた簡q)/2 とそれぞれ置き直して回帰分析にかける」を適用す便法を適用する。この方法も，定量限界値および検出限る。第 2 節の人工例に簡便法を適用した結果を図 4 に示

界値そのものを用いるのに比べより妥当な結果を与えるす。と推察されるが，本論ではこれに言及する余裕がない。
場を改めて議論したい。

30.025.0
20.0

本論の第 3 節で用いた回帰パラメータの最尤推定値導15.0

出のための EM アルゴリズムは SAS を用いて書いた。興
味をもたれた読者にはそれを提供する用意があるので，

10.0
y = 1.6587x + 1.2925
2R = 0.7894

5.0

著者まで連絡されたい。

0.00.0

2.0

4.0

6.0

8.0

10.012.0

14.0

16.0

謝辞

図 4. 簡便法による回帰分析
情報処理専攻修士学生の中戸川智彦君には本論文での求められた回帰直線は計算をお手伝いいただいた。また，本研究には科学研究

y = 1.2925 + 1.6587 x

費基盤研究 (A) No. 16200022 の援助を受けた。ここに感であり，第 3 節の EM アルゴリズムによる推定値に近く，謝する。

限界を無視した場合の傾きの過小評価も解消されている。

参考文献

この簡便法は，まさに簡便でありしかも上記の例ではbest possible と思われる EM アルゴリズムによる最尤推Counts, M. E., Hsu, F. S., Laffoon, S. W., Dwyer, R. W. and定値に近い推定値を与えている。

Cox, R. H. (2004) Mainstream smoke constituent yields

５．おわりにand predicting relationships from a worldwide market
sample of cigarette brands: ISO smoking conditions.

不完全データに対する EM アルゴリズムによる最尤推Regulatory Toxicology and Pharmac- ology, 39, 111-134.

定や多重代入法による対処などは，近年の統計理論の整

Dempster, A. P., Laird, N. M. and Rubin, D. B. (1977)

備ならびにソフトウェアの開発により実用に供されつつMaximum Likelihood Estimation from Incomplete Data

ある（Little and Rubin (2002)，
岩崎 (2002)，Rubin (1987) な

via the EM Algorithm. Journal of the Royal Statistical

どを参照）。しかし，本研究で取り上げたような NMAR

Society, Series B, 39, 1-38 (with discussion).
岩崎 学 (2002) 不完全データの統計解析，エコノミストとなる欠測（打ち切り）では，通常のソフトウェアでは対処できず，そのためのプログラムを個別に書かなけれ社．

ばならない。それは，科学的な研究には不可欠ではある

Little, R. J. A. and Rubin, D. B. (2002) Statistical Analysis

が，現実問題では実行が困難である。そこで，それらはwith Missing Data, Second Edition. John Wiley & Sons,

理想的な対処法であることは知った上で，なるべくそれ

New York.に近い結果を与えるより簡便な方法があればそれを採用

Johnson, N. L. and Kotz, S. (1970) Distributions in Statistics.

するのが現実的な選択となるであろう。本研究の視点のContinuous Univariate Distributions – 1. John Wiley &ひとつはそこにあった。

Sons, New York.
Rubin, D. B. (1987) Multiple Imputation for Nonresponse in

実は厄介な問題がある。検出限界や定量限界が存在す

Surveys. John Wiley & Sons, New York.

る問題では，y と x の値のいずれかもしくは両方がきわめて 0 に近いことが予想され，生の観測値に対して (3.1)－21－