變異數成分分析¶
此筆記本說明了用於雙層巢狀和交叉設計的變異數成分分析。
[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 的所有預設引數，「群組 1 變異數」和「群組 2 變異數」的母體值分別為 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])
/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()
最後我們擬合模型。 可以看到，這兩個擬合的結果是相同的。
[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 中的群組始終是巢狀的，但可以通過只有一個群組並將所有隨機效應指定為變異數成分來擬合交叉模型。 許多但並非所有交叉模型都可以用這種方式擬合。 下面的函數會產生一個具有兩個隨機結構層級的交叉數據集。
[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 日