分散成分分析¶
このノートブックでは、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ラベルは、たとえ同じラベルを持っていても、独立したグループとして扱われます。たとえば、2つの異なる学区にある「学校1」というラベルの2つの学校は、同じラベルを持っていても独立した学校として扱われます。
[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 のすべてのデフォルト引数を使用すると、「グループ1の分散」と「グループ2の分散」の母集団値はそれぞれ 2^2=4 と 3^2=9 になります。要約表の先頭に「スケール」としてリストされている説明されない分散の母集団値は 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])
/tmp/ipykernel_4815/1119967950.py:1: DeprecationWarning: DataFrameGroupBy.apply operated on the grouping columns. This behavior is deprecated, and in a future version of pandas the grouping columns will be excluded from the operation. Either pass `include_groups=False` to exclude the groupings or explicitly select the grouping columns after groupby to silence this warning.
vcm = df.groupby("group1").apply(f).to_list()
最後に、モデルを適合させます。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つだけにし、すべてのランダム効果を分散成分として指定することで、交差モデルを適合させることができます。多くの交差モデルは、この方法で適合させることができますが、すべてではありません。以下の関数は、2レベルのランダム構造を持つ交差データセットを生成します。
[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
==========================================================
最終更新日: 2024年10月03日