最小2乗平均の利用 | 実験計画法

JMP および R を使って、最小2乗平均(least square mean, lsmean, LSM)を説明します。

最小2乗平均(調整平均、調整済平均)は、最小2乗法によって線形モデルをデータにあてはめ、その推定式から得られた効果の線形結合として求められます。したがって、あてはめる線形モデルが異なると、最小2乗平均も異なります。

バランスがとれたデータの場合、通常の平均と最小2乗平均は同値になります。
アンバランスなデータの場合、平均と最小2乗平均は一致しません。最小2乗平均は平均の過小評価あるいは過大評価を調整した値なので、この場合は最小2乗平均を用います。

また、共分散分析などでも最小2乗平均が用いられます。

以下、具体的にデータを使って説明します。

データ

2因子実験(質的因子 A ×質的因子 B)のデータを使って、最小2乗平均を説明します。因子 A が3水準(A1、A2、A3)、因子 B が2水準(B1、B2)、繰り返し数が3の完全無作為化法で実験を行い、表1の観測値を得ました。元々はバランスがとれた実験ですが、A2 と B1 の組合せと、A3 と B2 の組合せ(オレンジのセル)で1つずつ欠測が生じ、アンバランスなデータになりました。
図1は、2因子の水準組合せごとの平均をプロットしたグラフです。

説明用データ

表2に、水準ごとの平均と総平均を示しました。表3に、水準の組合せごとの「平均」とその平均(以下、「平均の平均」)および総平均を示しました。欠測がない水準 A1 の「平均」と「平均の平均」は 106.83 で一致します(ブルーのセル)。しかし、他の水準では欠測の影響で、「平均」と「平均の平均」が一致しません。たとえば、A2 の「平均」は 108.00 ですが、「平均の平均」は 108.58 です(オレンジのセル)。総平均も一致しません。この「平均の平均」の値は、後で使います。

平均と平均の平均

モデル式

この2因子実験のデータに、次ぎのモデル式をあてはめます。y は、A の i 番目の水準と B の j 番目の水準の組合せで実験し、k 番目に得られた観測値です。m、a、b、ab は、それぞれ全体平均、A の主効果、B の主効果、A と B の交互作用を示すパラメータ、e はランダムな誤差です。また、制約条件として、a、b、ab を、i または j について加えたものを 0 とします。すなわち、効果の合計が 0 ということです。

$$y_{ijk}=m+a_i+b_j+(ab)_{ij}+e_{ijk}$$

$$\sum_{i}a_i=0,\ \sum_{j}b_j=0,\ \sum_{i} (ab)_{ij} = 0,\ \sum_{j} (ab)_{ij}=0 $$

なお、繰り返し数が水準によって異なる場合、一般には、次の制約条件を設定します。これが大きな相違点です。ここで、添え字 ij の付いた n は、A の i 番目の水準と B の j 番目の水準の組合せにおける繰り返し数です。ドットの付いた n はその周辺度数です。詳細は、高橋ら(1989)を参照してください。

$$\sum_{i}n_{i \cdot} a_i=0,\ \sum_{j}n_{\cdot j} b_j=0,\ \sum_{i}n_{ij} (ab)_{ij} = 0,\ \sum_{j}n_{ij} (ab)_{ij}=0 $$

JMPの解析結果と最小2乗平均

このモデル式を使って JMP で解析した結果のうち、最小2乗平均に関する部分として、モデル式のパラメータの推定値(全水準の推定値)と、各水準の最小2乗平均表を示します。

JMP全水準の推定

JMP最小2乗平均表

ある水準の最小2乗平均は、他の因子を平均に固定したときの期待値です。
例として、下の推定式により、水準 A2 の最小2乗平均を計算します。因子 B と交互作用を平均に固定します。推定式にある y の添え字のドットは、平均であることを示します。モデルの制約条件から、b の項の合計は 0、ab の項の合計も 0 なので、それぞれの平均も 0 になります。したがって、JMP 出力の「全水準の推定値」で表示されている切片 110.75 と A[A2] の値 -2.17(オレンジ枠)を代入して、108.58 が得られます。最小2乗平均表の A2 の値(オレンジ枠)と一致します。

$$E[y_{2 \cdot \cdot}] =  m+a_2 +\overline{b_j}+\overline{(ab)_{2j}}= 110.75+2.17 = 108.58$$

A2 と B1 の組合せの最小2乗平均は、次のようにして求められます。

$$E[y_{21 \cdot}]=m+a_2+b_1+(ab)_{21}=110.75-2.17-2.64+0.28=111.50$$

最小2乗平均の特徴

最小2乗平均は、欠測による平均の過大評価または過少評価を調整した値になっています。

この事例で、因子 B の観測値は、B1 が B2 よりも高い傾向にあります(図1)。したがって、B1 に欠測が生じると、この組合せの平均は過少評価になります。また、B2 に欠測が生じると、この組合せの平均は過大評価になります。

たとえば、水準 A2 のデータ数は欠測により5個、このうち B1 との組合せが2個、B2 との組合せが3個(表1)、これらを単純に平均した値 108.00 は過小評価になっています(表2)。この最小2乗平均は 108.58 でした。過小評価を調整して、平均よりも高い値になっています。

一方、水準 A3 の欠測は B2 との組合せにあるため、平均 117.40 は過大評価になっています。最小2乗平均 116.83 は、この過大評価を調整して平均よりも低い値になっています。

JMP の最小2乗平均表に表示されている平均は、表2の水準ごとの「平均」に一致しています。一方、最小2乗平均は、表3の水準の組合せごとの「平均の平均」に一致しています。このことから、最小2乗平均は「平均の平均」であると説明されることがあります。しかし、これは便宜的な説明です。
仮に、交互作用の項を除いたモデル式を用いて、上記と同様の解析を行うと、最小2乗平均は「平均の平均」と一致しません。

Rの解析結果と最小2乗平均

上記の「JMP の解析結果と最小2乗平均」の内容を基にして、R を使って最小2乗平均を説明します。

R スクリプトとその表示結果の一部を示します。

R最小2乗平均

表1のデータをデータフレーム df に付値します。変数 A と変数 B は因子型、変数 y は数値型です。

lm 関数でモデル式をあてはめて、結果をオブジェクト lm_out に付値します。引数 contrasts で因子 A と B を “contr.sum” に指定することにより、効果の和が 0 になるというモデルの制約条件を設定します。

lm_out から summary 関数によって、モデル式のパラメータの推定値を表示させると、JMP の「全水準の推定値」に相当する結果が得られます。

ただし、A3 のパラメータ推定値が省略されています。A1 のパラメータ推定値が  -3.9167、A2 の推定値が -2.1677、効果の和が 0 という制約条件により、A3 の推定値は -(-3.9167)-(-2.1667)=6.0833 になります。これらのパラメータ推定値から、JMP の出力で説明したように、効果の線形結合で最小2乗平均が求められます。

また、emmeans パッケージの emmeans 関数に lm_out を渡すと、最小2乗平均が得られます。さらに、emmeans 関数により、水準間の多重比較を行うことができます。

詳細は、R スクリプトと Excel ファイルをダウンロードして、RStudio でコンパイルして確認してください(ダウンロード)。

このサイトでは、芳賀敏郎著「医薬品開発のための統計解析」(サイエンティスト社)の理解を支援する PDF ファイルを提供しています。この本(グリーン本)は、第1部~第3部で構成されています。
第2部「実験計画法」の以下の節に、「最小2乗平均」の説明を加えました。

第3章「乱塊法実験」3節「欠測値のある場合」
第4章「共分散分析」2節「解析手順」
第5章「2因子実験」2節「2因子実験(質的因子×質的因子)」
第5章「2因子実験」3節「2因子実験(質的因子×量的因子)」

2因子実験のモデル式については、以下の節を参照してください。

第2章「量的因子の1因子実験」3節「ダミー変数による質的因子の効果の推定」
第5章「2因子実験」1節「2因子実験の基礎」

RとRStudio および R スクリプトのコンパイルについては、以下を参照してください。

「RとRStudioの使い方」

参考文献
高橋行雄・大橋靖雄・芳賀敏郎(1989)SAS による実験データの解析、東京大学出版会
(2022年2月21日)