分散成分分析¶

このノートブックでは、2レベルのネストかつ交差された設計の分散コンポーネント分析について説明します。

[1]:

import numpy as np
import statsmodels.api as sm
from statsmodels.regression.mixed_linear_model import VCSpec
import pandas as pd

ノートブックを再現可能にします

[2]:

np.random.seed(3123)

ネスト分析¶

以下の議論では、「グループ 2」は「グループ 1」の中にネストされています。具体的な例として、「グループ 1」は学校区で、「グループ 2」は個々の学校である場合があります。以下の関数は、そのような母集団からデータを生成します。ネストされた分析では、異なる「グループ 1」のラベルにネストされた「グループ 2」のラベルは、同じラベルであっても独立したグループとして扱われます。例えば、異なる学校区にある「学校 1」とラベル付けされた二つの学校は、同じラベルが付けられていても独立した学校として扱われます。

[3]:

def generate_nested(
    n_group1=200, n_group2=20, n_rep=10, group1_sd=2, group2_sd=3, unexplained_sd=4
):

    # Group 1 indicators
    group1 = np.kron(np.arange(n_group1), np.ones(n_group2 * n_rep))

    # Group 1 effects
    u = group1_sd * np.random.normal(size=n_group1)
    effects1 = np.kron(u, np.ones(n_group2 * n_rep))

    # Group 2 indicators
    group2 = np.kron(np.ones(n_group1), np.kron(np.arange(n_group2), np.ones(n_rep)))

    # Group 2 effects
    u = group2_sd * np.random.normal(size=n_group1 * n_group2)
    effects2 = np.kron(u, np.ones(n_rep))

    e = unexplained_sd * np.random.normal(size=n_group1 * n_group2 * n_rep)
    y = effects1 + effects2 + e

    df = pd.DataFrame({"y": y, "group1": group1, "group2": group2})

    return df

分析するデータセットを生成します。

[4]:

df = generate_nested()

generate_nestedのすべてのデフォルト引数を使用した場合、」group 1 Var」と「group 2 Var」の母集団値はそれぞれ \(2^2=4\)と\(3^2=9\) です。概要テーブルの上部に「scale」として記載されている説明されていない分散は、母集団値が \(4^2=16\) です。

[5]:

model1 = sm.MixedLM.from_formula(
    "y ~ 1",
    re_formula="1",
    vc_formula={"group2": "0 + C(group2)"},
    groups="group1",
    data=df,
)
result1 = model1.fit()
print(result1.summary())

          Mixed Linear Model Regression Results
==========================================================
Model:            MixedLM Dependent Variable: y
No. Observations: 40000   Method:             REML
No. Groups:       200     Scale:              15.8825
Min. group size:  200     Log-Likelihood:     -116022.3805
Max. group size:  200     Converged:          Yes
Mean group size:  200.0
-----------------------------------------------------------
            Coef.   Std.Err.    z     P>|z|  [0.025  0.975]
-----------------------------------------------------------
Intercept   -0.035     0.149  -0.232  0.817  -0.326   0.257
group1 Var   3.917     0.112
group2 Var   8.742     0.063
==========================================================

数式インターフェースを避けたい場合は、デザイン行列を手動で構築して同じモデルを適合させることができます。

[6]:

def f(x):
    n = x.shape[0]
    g2 = x.group2
    u = g2.unique()
    u.sort()
    uv = {v: k for k, v in enumerate(u)}
    mat = np.zeros((n, len(u)))
    for i in range(n):
        mat[i, uv[g2.iloc[i]]] = 1
    colnames = ["%d" % z for z in u]
    return mat, colnames

次に、VCSpec クラスを使用して分散コンポーネントを設定します。

[7]:

vcm = df.groupby("group1").apply(f).to_list()
mats = [x[0] for x in vcm]
colnames = [x[1] for x in vcm]
names = ["group2"]
vcs = VCSpec(names, [colnames], [mats])

最後にモデルをフィットさせます。2 つのフィッティングの結果は同一であることがわかります。

[8]:

oo = np.ones(df.shape[0])
model2 = sm.MixedLM(df.y, oo, exog_re=oo, groups=df.group1, exog_vc=vcs)
result2 = model2.fit()
print(result2.summary())

          Mixed Linear Model Regression Results
==========================================================
Model:            MixedLM Dependent Variable: y
No. Observations: 40000   Method:             REML
No. Groups:       200     Scale:              15.8825
Min. group size:  200     Log-Likelihood:     -116022.3805
Max. group size:  200     Converged:          Yes
Mean group size:  200.0
-----------------------------------------------------------
            Coef.   Std.Err.    z     P>|z|  [0.025  0.975]
-----------------------------------------------------------
const       -0.035     0.149  -0.232  0.817  -0.326   0.257
x_re1 Var    3.917     0.112
group2 Var   8.742     0.063
==========================================================

交差分析¶

交差分析では、一つのグループのレベルが別のグループのレベルと任意の組合せで発生する可能性があります。Statsmodels MixedLM のグループは常にネストされていますが、1つのグループのみを持ち、すべての変量効果を分散コンポーネントとして指定することで、交差モデルを適合させることができます。すべてではありませんが、多くの交差モデルをこの方法で適合させることができます。次の関数は、二つのレベルのランダム構造を持つ交差データ・セットを生成します。

[9]:

def generate_crossed(
    n_group1=100, n_group2=100, n_rep=4, group1_sd=2, group2_sd=3, unexplained_sd=4
):

    # Group 1 indicators
    group1 = np.kron(
        np.arange(n_group1, dtype=int), np.ones(n_group2 * n_rep, dtype=int)
    )
    group1 = group1[np.random.permutation(len(group1))]

    # Group 1 effects
    u = group1_sd * np.random.normal(size=n_group1)
    effects1 = u[group1]

    # Group 2 indicators
    group2 = np.kron(
        np.arange(n_group2, dtype=int), np.ones(n_group2 * n_rep, dtype=int)
    )
    group2 = group2[np.random.permutation(len(group2))]

    # Group 2 effects
    u = group2_sd * np.random.normal(size=n_group2)
    effects2 = u[group2]

    e = unexplained_sd * np.random.normal(size=n_group1 * n_group2 * n_rep)
    y = effects1 + effects2 + e

    df = pd.DataFrame({"y": y, "group1": group1, "group2": group2})

    return df

分析するデータセットを生成します。

[10]:

df = generate_crossed()

次に、モデルをフィットさせます。 groups ベクトルは一定であることに注意してください。 generate_crossed のデフォルトパラメータを使用すると、レベル 1 の分散は \(2^2=4\)、レベル 2 の分散は \(3^2=9\) 、原因不明の分散は \(4^2=16\) である必要があります。

[11]:

vc = {"g1": "0 + C(group1)", "g2": "0 + C(group2)"}
oo = np.ones(df.shape[0])
model3 = sm.MixedLM.from_formula("y ~ 1", groups=oo, vc_formula=vc, data=df)
result3 = model3.fit()
print(result3.summary())

          Mixed Linear Model Regression Results
==========================================================
Model:            MixedLM Dependent Variable: y
No. Observations: 40000   Method:             REML
No. Groups:       1       Scale:              15.9824
Min. group size:  40000   Log-Likelihood:     -112684.9688
Max. group size:  40000   Converged:          Yes
Mean group size:  40000.0
-----------------------------------------------------------
            Coef.   Std.Err.    z     P>|z|  [0.025  0.975]
-----------------------------------------------------------
Intercept   -0.251     0.353  -0.710  0.478  -0.943   0.442
g1 Var       4.282     0.154
g2 Var       8.150     0.291
==========================================================

数式インターフェースを避けたい場合は、設計マトリックスを手動で構築して同じモデルを適合させることができます。

[12]:

def f(g):
    n = len(g)
    u = g.unique()
    u.sort()
    uv = {v: k for k, v in enumerate(u)}
    mat = np.zeros((n, len(u)))
    for i in range(n):
        mat[i, uv[g[i]]] = 1
    colnames = ["%d" % z for z in u]
    return [mat], [colnames]


vcm = [f(df.group1), f(df.group2)]
mats = [x[0] for x in vcm]
colnames = [x[1] for x in vcm]
names = ["group1", "group2"]
vcs = VCSpec(names, colnames, mats)

ここでは数式を使用せずにモデルをフィットさせ、モデル 3 と 4 の結果が同一であることを簡単に確認できます。

[13]:

oo = np.ones(df.shape[0])
model4 = sm.MixedLM(df.y, oo[:, None], exog_re=None, groups=oo, exog_vc=vcs)
result4 = model4.fit()
print(result4.summary())

          Mixed Linear Model Regression Results
==========================================================
Model:            MixedLM Dependent Variable: y
No. Observations: 40000   Method:             REML
No. Groups:       1       Scale:              15.9824
Min. group size:  40000   Log-Likelihood:     -112684.9688
Max. group size:  40000   Converged:          Yes
Mean group size:  40000.0
-----------------------------------------------------------
            Coef.   Std.Err.    z     P>|z|  [0.025  0.975]
-----------------------------------------------------------
const       -0.251     0.353  -0.710  0.478  -0.943   0.442
group1 Var   4.282     0.154
group2 Var   8.150     0.291
==========================================================

最終更新日: 2025年01月28日