GLMM の基礎

GLMM タブで使われている統計理論の背景です。操作方法は GLMM のページを参照してください。

階層データと独立性の問題

統計モデルの多くは、観測が互いに独立であることを仮定しています。しかし現実のデータには、観測がグループに属する階層構造を持つものが多くあります。

複数の学校に通う生徒のテスト成績（生徒はそれぞれの学校に所属）
複数の病院で治療を受けた患者の回復日数（患者はそれぞれの病院に所属）
同じ被験者に対する反復測定（測定は被験者に所属）

同じグループに属する観測は、グループ固有の要因（学校の教育方針、病院の設備、被験者の体質）を共有するため、互いに似た値を取る傾向があります。この相関を無視して GLM で分析すると、モデルの誤特定（misspecification）になります。データの生成過程にグループ構造があるのに、モデルがそれを表現していない状態です。

この誤特定は、推定結果を 2 つの形で損ないます。

標準誤差の過小評価: グループ構造を無視した GLM は、観測間の共分散をゼロと仮定したまま標準誤差を計算します。実際には同じグループ内の観測に正の相関があるため、切片やグループ間で変動する説明変数の係数で標準誤差が過小になり、95% 信頼区間は実際には 95% の被覆率を持ちません。過小の度合いは説明変数の変動がグループ間とグループ内のどちらにあるかで変わります（ICC の節のデザイン効果を参照）。

固定効果の偏り: グループレベルの交絡要因があると、固定効果の推定に欠落変数バイアス（omitted variable bias）が生じます。たとえば学校の教育方針が勉強時間と成績の両方に影響する場合、学校差を無視した回帰では勉強時間の効果を正しく推定できません。ランダム切片はグループ間の平均レベルの違いをモデル化しますが、その推定はランダム切片が説明変数と無相関であることを仮定しています。未測定のグループレベル交絡が説明変数と相関している場合はこの仮定が崩れ、固定効果の推定は偏ったままです。測定された交絡要因は固定効果として投入する必要があります。

混合モデル（mixed model）は、グループ内の相関をランダム効果として明示的にモデル化することで、標準誤差の過小評価を解消します。あわせて、ばらつきをグループ間とグループ内に分けて推定するため、グループの違いに由来するばらつきの大きさを定量的に把握できます。この分解を要約する指標が後述の ICC です。

固定効果とランダム効果

混合モデルの「混合」は、固定効果とランダム効果の両方を含むことを指します。

固定効果（fixed effects） は、どのグループにも共通に働くと仮定する、説明変数と応答の体系的な関係です。係数は未知の定数としてモデル化されます。たとえば「勉強時間が 1 時間増えるとテスト成績が何点上がるか」という効果は固定効果です。

ランダム効果（random effects） は、直接観測できないグループレベルのばらつきをモデル内で表現するための装置で、変量効果とも呼ばれます。各学校の平均成績は教育方針や生徒層によって異なります。このばらつきは確率的なメカニズムで生まれたものではなく、具体的な要因に起因しますが、それらの要因をすべて説明変数として測定・投入するのは現実的ではありません。ランダム効果はこの「観測できないグループ差」を確率分布で近似的に表現します。

グループ $j$ の差を $u_j$ と書くと、 $u_j \sim N(0, \sigma_u^2)$ という仮定は「学校差が正規分布から生成されている」という主張ではありません。グループ間のばらつきの大きさ $\sigma_u^2$ を推定し、それに基づいて縮小推定（後述）を行うための、モデリング上の仮定です。

両者の実質的な違いは、パラメータの数え方に現れます。グループをダミー変数にして固定効果として扱うと、グループの数だけ自由なパラメータが増え、各グループの効果はそれ自体が推定対象になります。ランダム効果では、推定する分散成分は $\sigma_u^2$ の 1 つだけで、各グループの $u_j$ は自由なパラメータではなく、この分散を前提として予測される値です。固定効果/ランダム効果の区別は、グループがランダムに抽出されたかどうかではなく、グループ間の差をこのどちらの形で扱うかの選択です（Gelman, 2005）。

ランダム切片モデル

最も基本的な混合モデルがランダム切片モデルです。各グループに固有の切片（ベースライン）を与えます:

g(\mu_i) = x_i'\beta + u_{j[i]}, \quad u_j \sim N(0, \sigma_u^2)

$g$ はリンク関数、 $x_i'\beta$ は固定効果の線形予測子です。 $j[i]$ は観測 $i$ が属するグループを表す添字で、 $u_{j[i]}$ はそのグループのランダム切片です。 $u_j \sim N(0, \sigma_u^2)$ は、グループごとの値 $u_1, \dots, u_J$ がそれぞれ独立に同じ正規分布に従うことを表します。

直感的には、全グループに共通する回帰直線（ $x_i'\beta$ ）を引き、各グループの直線はそこから $u_j$ だけ上下にずれる、という構造です。正規分布 $N(0, \sigma_u^2)$ の仮定は、このずれの大きさを $\sigma_u^2$ という1つのパラメータで要約するための道具です。各 $u_j$ がどんな値をとるかは、データから推測します（ランダム効果の予測と縮小推定を参照）。

Gaussian 分布族の場合、これは線形混合モデル（LMM）です:

y_i = x_i'\beta + u_{j[i]} + \varepsilon_i, \quad u_j \sim N(0, \sigma_u^2), \quad \varepsilon_i \sim N(0, \sigma_e^2)

$\sigma_u^2$ はグループ間のばらつき、 $\sigma_e^2$ はグループ内（個人レベル）のばらつきです。残差分散 $\sigma_e^2$ は全グループ共通と仮定されます。ランダム切片モデルが表すのはグループ間の平均レベルの違いであり、グループ内のばらつきの大きさがグループごとに異なる構造はモデル化しません。

同じグループの観測が互いに似るという性質は、この式では共分散として表れます。同じグループの 2 つの観測は $u_j$ を共有するため、その共分散は $\sigma_u^2$ です。異なるグループの観測どうしの共分散は 0 です。各観測の分散は $\sigma_u^2 + \sigma_e^2$ なので、同一グループ内の 2 観測の相関は $\sigma_u^2 / (\sigma_u^2 + \sigma_e^2)$ となり、この量が後述の ICC です。

リンク関数が恒等でない場合、係数 $\beta$ の解釈は条件付きになります。 $\beta$ が表すのは、同じグループの中で説明変数を動かしたときのリンク尺度上の効果です。 $g^{-1}$ が非線形なため、グループを均した集団平均の応答に対する効果は、 $\sigma_u^2 > 0$ である限り一般にこれと一致しません。たとえば Binomial + logit で係数から得られるオッズ比は、グループをまたいだ平均的なオッズの比ではなく、同一グループ内での比較として読みます。

パラメータ推定

混合モデルのパラメータ推定は、固定効果 $\beta$ と分散成分（ $\sigma_u^2$ , $\sigma_e^2$ ）を同時に推定する必要がある点で、通常の回帰モデルの推定より複雑です。

REML（制限付き最尤法）

最尤法（ML）は分散成分を過小推定する傾向があります。 $\beta$ の推定に自由度を消費することを尤度関数が考慮しないためです。分母に $n$ を使った標本分散が分散を過小推定するのと同じ現象で、標本分散では分母を $n-1$ にして補正します。

REML はこの偏りを避けるため、 $\beta$ の影響を受けない量だけから分散成分を推定します。応答ベクトル $y$ を $n$ 次元空間の点とみると、固定効果 $X\beta$ が動けるのは計画行列 $X$ の列が張る $p$ 次元の部分空間です。 $y$ をその直交補空間へ射影した残差は、 $X\beta$ の成分がちょうど消えるため、分布が $\beta$ に依存せず分散成分だけで決まります。REML はこの残差の尤度を最大化します。 $\beta$ の推定に自由度を消費しないので、分散の推定は実質 $n - p$ 次元のデータに基づき、偏りが補正されます。標本分散の $n - 1$ は、 $X$ が切片のみ（ $p = 1$ ）の場合にあたります。

MIDAS は REML の最大化を、探索を 1 次元に落とした profile REML で行います。分散比 $\theta^2 = \sigma_u^2/\sigma_e^2$ を 1 つ与えるごとに、その下で最適な $\beta$ と $\sigma_e^2$ は閉形式で決まります。数値的に探索するのは $\log\theta$ だけで、黄金分割法で最大点を求めます。

Laplace 近似と PIRLS

Gaussian 分布族と identity リンクの組み合わせではランダム効果の積分が解析的に実行できるため REML で十分ですが、それ以外の組み合わせ（Binomial, Poisson, Gamma, および Gaussian + log）では、ランダム効果を積分除去する操作:

L(\beta, \sigma_u^2) = \int \prod_i f(y_i \mid \beta, u) \cdot f_u(u) \, du

に閉形式（closed form）の解がありません。 $f$ は選択した分布族における各観測の密度関数（Binomial や Poisson などの離散分布族では確率質量関数）、 $f_u$ はランダム効果の密度関数、すなわち $N(0, \sigma_u^2)$ の密度です。

Laplace 近似はこの積分を、被積分関数の最大点周辺での2次近似で置き換えます。近似の精度は、グループあたりの観測数が少ない場合や、Binomial でイベントが稀な場合に低下することがあります。

MIDAS の推定アルゴリズムは、Bates et al. (2015) が記述した入れ子の 2 ループ構造にもとづき、2 つのループを重ねて動きます。外側ループの最適化手法は MIDAS 独自の選択です。

外側ループ: 相対共分散パラメータ $\theta$ の対数を黄金分割法で最適化
内側ループ: PIRLS（Penalized IRLS）で $(\beta, u)$ を同時推定

PIRLS は GLM の IRLS にランダム効果のペナルティ項を加えたもので、次のペナルティ付き逸脱度を最小化します:

D_p = \sum_i d(y_i, \mu_i) + \sum_j \frac{u_j^2}{\theta^2}

$\theta$ は外側ループが最適化する相対共分散パラメータで、ランダム効果の標準偏差を残差の尺度に対する比で表した量です。Gaussian では $\theta = \sigma_u/\sigma_e$ 、Poisson/Binomial では $\theta = \sigma_u$ 、Gamma では $\theta = \sigma_u/\sqrt\phi$ （ $\phi$ は分散パラメータ）です。第1項はデータへの適合度、第2項は $u_j$ が正規分布から大きく外れることへのペナルティです。このペナルティにより、データの少ないグループの推定値は全体平均に引き寄せられます（後述の縮小推定を参照）。

境界解（singular fit）

分散 $\sigma_u^2$ の推定には 0 という下限があります。どちらの推定経路でも、最適化の解がこの下限に張り付いた状態を singular fit と呼びます。グループ間の変動が残差の変動に比べて無視できるほど小さいデータで起きます。

境界解は、推定結果の読み方を変えます。 $\hat\sigma_u^2 = 0$ は「グループ差が存在しない」という結論ではなく、このデータではグループ間変動を残差変動と区別できなかったという状態です。 $\hat\sigma_u^2$ に依存する ICC やランダム効果の予測値も、ほぼ 0 の値になります。また、分散成分に関する漸近的な近似は真のパラメータが境界の内部にあることを前提とするため、境界解ではこの前提が成り立ちません。この状態の検出と対処については GLMM タブのページを参照してください。

AIC/BIC の制限

REML の射影は固定効果の構成に依存するため、REML ベースの AIC/BIC は固定効果が異なるモデル間で比較できません。MIDAS が REML で推定するのは Gaussian + identity だけなので、この制約が当たるのはこの組み合わせです。それ以外の組み合わせでは AIC/BIC を最尤法ベースの量（Laplace 近似した周辺対数尤度）から導出するため、固定効果が異なるモデルどうしを比較できます。

ただし、どちらの経路でも比較できるのは同一の分布族とリンクの中だけです。分布族やリンクが変わると、AIC/BIC の差にモデルの当てはまり以外の要因が混ざります。推定経路が変われば対数尤度の基準（REML か最尤か）が揃わず、分布族が変われば対数尤度の測る対象（離散分布の確率質量か連続分布の密度か）や分散パラメータの扱いが変わり、リンクが変われば Laplace 近似の誤差の出方が変わるためです。

ランダム効果の予測と縮小推定

当てはめの結果には、固定効果の係数だけでなく、グループごとの $u_j$ の値も含まれます。 $u_j$ の値は、そのグループのデータだけから計算されるのではありません。 $u_j$ には確率分布 $N(0, \sigma_u^2)$ を仮定しているため、データとその分布という 2 つの情報源を組み合わせて値を推測します。この組み合わせが、この節で説明する縮小推定を生みます。

固定効果との扱いの違いは、用語にも現れます。固定効果 $\beta$ は未知の定数なので「推定（estimation）」と呼び、 $u_j$ はモデル内で確率変数として扱われるため「予測（prediction）」と呼び分けます。 $\beta$ ではデータの有限性だけが不確実性の原因なのに対し、 $u_j$ では仮定した分布も情報源として働く、という違いを反映した形式的な区別で、 $u_j$ が本当にランダムに生成されたという主張ではありません。

Gaussian（LMM）では、ランダム効果の予測に BLUP（Best Linear Unbiased Predictor、最良線形不偏予測量）を使います。BLUP は、 $u_j$ に対する不偏性 $E[\hat u_j] = E[u_j] = 0$ が未知の固定効果 $\beta$ のどの値でも成り立つという制約のもとで、平均二乗予測誤差 $E[(\hat u_j - u_j)^2]$ を最小化する線形予測量で、Henderson (1975) の混合モデル方程式から導かれます。この予測量は縮小推定（shrinkage estimator）の形をとります:

\hat{u}_j = \frac{n_j}{n_j + \sigma_e^2 / \sigma_u^2} \times \bar{r}_j

$n_j$ はグループ $j$ のサイズ、 $\bar{r}_j$ はグループ $j$ の平均残差（固定効果で説明できない部分）です。

係数 $n_j / (n_j + \sigma_e^2/\sigma_u^2)$ は0から1の値をとり、グループサイズが大きいほど1に近づきます。つまり:

大きなグループ: データが十分にあるので、そのグループ固有の推定値をほぼそのまま使う
小さなグループ: データが少ないので、全体平均（ゼロ）に向かって引き寄せる

これは「情報の借用（borrowing strength）」とも呼ばれます。この言い回しは、Tukey が 1960 年代の選挙速報予測の仕事で縮小推定の考え方に与えた呼び名に由来します（Brillinger, 2002）。データの少ないグループは他のグループの情報を借りて推定を安定させます。グループ固有の推定値をそのまま使うと分散が大きくなりますが、全体平均だけを使うとグループの特性を無視してしまいます。BLUP はこのバイアスとバリアンスのトレードオフを最適にバランスさせます。

予測値には不確かさの評価が付きます。縮小の係数を $c_j = n_j / (n_j + \sigma_e^2/\sigma_u^2)$ と書くと、 $\hat u_j$ の標準誤差は、観測データで条件付けた $u_j$ の分散 $\sigma_u^2 (1 - c_j)$ の平方根です。推定した分散成分は真の値として扱います。グループのデータが少ないほど $c_j$ は 0 に近づくため、予測値は 0 へ縮み、標準誤差は分布の標準偏差 $\sigma_u$ に近づきます。データが増えるほど標準誤差は 0 へ向かいます。

Non-Gaussian の場合、ランダム効果は条件付きモード（conditional mode）として推定されます。観測された応答で条件付けたときの $u_j$ の分布において、密度が最大になる値のことです。仮定した分布 $N(0, \sigma_u^2)$ を事前分布とみなせば、事後分布の最頻値にあたります。条件付きモードは PIRLS の反復解として得られ、BLUP と同様に全体平均への縮小を示しますが、線形予測量ではないため厳密には BLUP ではありません。予測値の不確かさは、条件付きモードの周りの正規近似から評価します。

ICC（級内相関係数）

ICC（Intraclass Correlation Coefficient）は、固定効果で説明される分を除いた残差側の分散のうち、グループ間の違いが占める割合です:

\text{ICC} = \frac{\sigma_u^2}{\sigma_u^2 + \sigma_e^2}

$\sigma_u^2$ と $\sigma_e^2$ は当てはめたモデルの分散成分です。固定効果 $X\beta$ を含むモデルでは、どちらも固定効果で説明された分散を除いた残差側の量になります。

ICC は、同一グループ内の 2 観測の相関でもあります。ランダム切片モデルの節で導いたとおり、同じグループの 2 観測の相関は $\sigma_u^2 / (\sigma_u^2 + \sigma_e^2)$ で、この定義式と同じ量です。級内相関係数という名前はこの性質を指しています。

この形の ICC を条件付き ICC（conditional ICC）と呼び、MIDAS の GLMM が計算する ICC もこれです。対になる無条件 ICC は、説明変数を含めない切片のみのモデルから求めます。条件付き ICC は固定効果の指定に依存するため、説明変数を足し引きすると値が変わり、無条件 ICC とは一般に一致しません。

Non-Gaussian の場合

Non-Gaussian 分布族では $\sigma_e^2$ をデータから直接推定できません。ICC が分散分解として成立するには、 $\sigma_u^2$ と $\sigma_e^2$ が同じ尺度上の同じ量のばらつきを分け合う成分である必要があります。この条件を満たす $\sigma_e^2$ は、リンク関数と分布の仮定から理論的に導出できる場合にしか得られません。

Binomial + logit/probit には閾値モデル（Goldstein et al. 2002）があります。観測されない連続的な潜在変数 $y^* = X\beta + u + e$ が $y = \mathbf{1}(y^* > 0)$ として二値応答を生成します。 $e$ の分布はリンク関数の選択から演繹的に決まります。logit リンクでは $e$ がロジスティック分布に従い分散は $\pi^2/3$ 、probit リンクでは $e$ が $N(0, 1)$ に従い分散は $1$ です。この枠組みでは $\sigma_u^2 + \sigma_e^2$ が潜在変数の全分散を表し、ICC は分散分解として成立します。

分布族 + リンク	残差分散	根拠
Binomial + Logit	$\pi^2/3 \approx 3.29$	ロジスティック分布（閾値モデル）
Binomial + Probit	$1$	標準正規分布（閾値モデル）

Binomial の ICC は潜在変数 $y^*$ の尺度上の値です。観測された二値応答どうしのグループ内相関はこれとは別の量で、二値化で情報が落ちる分だけ潜在尺度の値より小さく、全体の発生確率にも依存します。

ICC が分散分解として成立するのは Gaussian + identity（REML で $\sigma_e^2$ を直接推定）と上記 2 つの Binomial の組み合わせのみです。それ以外の family+link では ICC が定義されません:

Poisson（全リンク）: Poisson 分布では $\operatorname{Var}(Y \mid \mu) = \mu$ であり、分散は平均から完全に規定されます。独立した残差分散パラメータが存在しません。ICC の分母に $1$ を入れる規約もありますが、 $\sigma_u^2 / (\sigma_u^2 + 1)$ は $\sigma_u^2$ の単調変換に過ぎず分散分解の解釈を持ちません。
Gamma（全リンク）: profiled な分散パラメータ $\phi$ は条件付き分散のパラメータであり、潜在尺度の残差分散ではありません。
Gaussian + log: $\sigma_u^2$ はリンク尺度上の量で、 $\sigma_e^2$ は応答尺度で推定された量です。両者を同じ尺度に揃える理論残差分散が存在しません。

デザイン効果と GLM との使い分け

グループ構造を無視した GLM で標準誤差がどれだけ過小になるかは、ICC とグループサイズから見積もれます。同じグループ内の観測は互いに似ていて独立な観測ほどの情報を持たないため、係数の推定量の分散は、同じサイズの独立なデータから推定した場合より大きくなります。この倍率がデザイン効果（DEFF）で、 $\text{DEFF} = 1 + (\bar n - 1) \times \text{ICC}$ （ $\bar n$ はグループの平均サイズ）と表せます。この式はグループサイズが均等な場合に正確で、不均等な場合は近似値です。独立を仮定した標準誤差の計算はこの増加を織り込まないため、真の標準誤差は報告値のおよそ $\sqrt{\text{DEFF}}$ 倍です。

分散がどれだけ増えるかは説明変数によって異なります。グループ平均やグループレベルで変動する説明変数の係数ではこの目安がそのまま当てはまり、グループ内で変動する説明変数の係数では増加は小さく、グループ間の変動をほとんど持たない説明変数ではほとんど増えません。

この見積もりは GLM と混合モデルの使い分けの材料になります。ICC が小さくグループサイズも大きくなければ DEFF は 1 に近く、グループ構造を無視して GLM で分析してもほぼ同じ結果が得られます。DEFF が 1 を大きく超えるなら、混合モデルが必要です。

固定効果の推測と正規近似

固定効果 $\hat\beta_j$ の信頼区間は Wald 法に基づき、 $\hat\beta_j \pm z_{1-\alpha/2} \cdot \text{SE}(\hat\beta_j)$ として構成します。 $\text{SE}(\hat\beta_j)$ は $(X'\hat V^{-1}X)^{-1}$ の第 $j$ 対角成分の平方根で、 $\hat V$ には推定した分散成分を代入します。Gaussian + identity では応答の共分散行列 $\hat V = \hat\sigma_u^2 ZZ' + \hat\sigma_e^2 I$ （ $Z$ は各観測のグループ所属を表す指示行列）、それ以外の組み合わせでは PIRLS 収束時の重み $\hat W$ から作る作業応答の近似共分散 $\hat V = \hat W^{-1} + \hat\theta^2 ZZ'$ を使い、分散パラメータ $\phi$ を推定する分布族ではさらに $\hat\phi$ を乗じます。この構成は、 $(\hat\beta_j - \beta_j) / \text{SE}(\hat\beta_j)$ が近似的に標準正規分布に従うことを根拠とします。

通常の回帰（OLS）では誤差分散 $\sigma^2$ を推定するため $(\hat\beta_j - \beta_j) / \text{SE}(\hat\beta_j)$ は t 分布に従い、自由度は $n - p$ です。混合モデルでは分散成分が複数あり ( $\sigma_u^2$ , $\sigma_e^2$ )、対応する自由度が一意に定まりません。そのため分散成分の推定値を真の値として扱い、標準正規分布で近似します。

LMM では $u_j$ と $\varepsilon_i$ に正規分布を仮定しているため、分散成分が既知なら $(\hat\beta_j - \beta_j) / \text{SE}(\hat\beta_j)$ は正規分布に正確に従います。Non-Gaussian GLMM では正規分布は漸近的な近似です。いずれの場合も、分散成分を推定値で置き換えると追加の不確実性が生じ、 $(\hat\beta_j - \beta_j) / \text{SE}(\hat\beta_j)$ の分布は標準正規より裾が重くなります。

グループ数が多ければ分散成分の推定が安定し、この近似は十分に正確です。グループ数が少ない場合は信頼区間が実際より狭くなる傾向があります。

ランダム傾きと交差ランダム効果

MIDAS の GLMM が対応するのはランダム切片モデルだけです。ランダム切片が表すのはグループごとのベースラインの違いまでなので、この節では対応範囲の外にある 2 つの拡張を、ランダム切片モデルで足りるかどうかを判断する材料として説明します。

説明変数の効果もグループごとに異なる場合は、ランダム傾き（random slope）モデルが必要です:

g(\mu_i) = x_i'\beta + u_{0j[i]} + u_{1j[i]} x_{1i}

$u_{0j}$ はランダム切片、 $u_{1j}$ はランダム傾きで、両者は多変量正規分布に従います。

傾きに群間変動があるデータをランダム切片モデルで済ませると、その変動は誤差項に入ります。このとき群内の共分散は実際には説明変数の値に依存しますが、ランダム切片モデルはそれを定数 $\sigma_u^2$ で置き換えます。共分散構造の誤特定が残るため、当の説明変数の係数の標準誤差は信頼できなくなります。

観測が入れ子になっていない複数のグルーピング変数に属する場合は、交差ランダム効果（crossed random effects）を使います。たとえば生徒が学校と地域の両方に属するなら、グルーピング変数ごとのランダム効果を線形予測子に並べて加えます:

g(\mu_i) = x_i'\beta + u_{j[i]} + v_{k[i]}, \quad u_j \sim N(0, \sigma_u^2), \quad v_k \sim N(0, \sigma_v^2)

$j[i]$ と $k[i]$ は観測 $i$ が属する学校と地域を表す添字で、2 つのランダム効果は互いに独立と仮定します。

ランダム切片モデルはグルーピング変数を 1 つしか持たないため、交差構造のデータでは、もう一方のグルーピング変数に由来する変動が残差に残ります。残った変動はその変数のグループ内に相関を生み、階層データと独立性の問題の節で述べたのと同じ形で、無視すれば標準誤差が過小になります。

参考文献

Gelman, A. (2005). Analysis of variance—why it is more important than ever. The Annals of Statistics, 33(1), 1-53. https://www.jstor.org/stable/3448650
Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1-48. https://www.jstatsoft.org/v67/i01/
Henderson, C. R. (1975). Best linear unbiased estimation and prediction under a selection model. Biometrics, 31(2), 423-447. https://www.jstor.org/stable/2529430
Brillinger, D. R. (2002). John W. Tukey: His life and professional contributions. The Annals of Statistics, 30(6), 1535-1575. https://doi.org/10.1214/aos/1043351246
Goldstein, H., Browne, W., & Rasbash, J. (2002). Partitioning variation in multilevel models. Understanding Statistics, 1(4), 223-231. https://doi.org/10.1207/S15328031US0104_02

このページの Markdown 版もあります。