論文を書いていたらAnovaRMはp値とF値は出力するものの、エフェクトサイズは出力しなくて困ってしまった。どうやらエフェクトサイズは書くものらしい。Rで計算しなおすのも一旦はいいや...ということで、AnovaRMを拡張してanova_tableにエフェクトサイズを出力させる。
[pystatsmodels] AnovaRM sum of squares to get effect size?
肝心のエフェクトサイズを出力するようにAnovaRMをカスタマイズするコードに関してはこちらのブログのものを流用した。このままではColabで使いにくいので、Colabでも実行できるようにする。
AnovaRMのサブクラスとしてCustomAnovaRMというクラスを作成、そこでfit()をオーバーライドする。
statsmodels/statsmodels/stats/anova.py at main · statsmodels/statsmodels · GitHub
これが元コードで、ここから色々importする。
! pip install patsy
import patsy
from scipy import stats
from statsmodels.stats.anova import _not_slice
from statsmodels.stats.anova import AnovaRM
from statsmodels.stats.anova import AnovaResults
from statsmodels.regression.linear_model import OLS
class CustomAnovaRM(AnovaRM):
def fit(self):
"""estimate the model and compute the Anova table
Returns
-------
AnovaResults instance
"""
y = self.data[self.depvar].values
# Construct OLS endog and exog from string using patsy
within = ['C(%s, Sum)' % i for i in self.within]
subject = 'C(%s, Sum)' % self.subject
factors = within + [subject]
x = patsy.dmatrix('*'.join(factors), data=self.data)
term_slices = x.design_info.term_name_slices
for key in term_slices:
ind = np.array([False] * x.shape[1])
ind[term_slices[key]] = True
term_slices[key] = np.array(ind)
term_exclude = [':'.join(factors)]
ind = _not_slice(term_slices, term_exclude, x.shape[1])
x = x[:, ind]
# Fit OLS
model = OLS(y, x)
results = model.fit()
if model.rank < x.shape[1]:
raise ValueError('Independent variables are collinear.')
for i in term_exclude:
term_slices.pop(i)
for key in term_slices:
term_slices[key] = term_slices[key][ind]
params = results.params
df_resid = results.df_resid
anova_table = pd.DataFrame(
for key in term_slices:
print(key)
if self.subject not in key and key != 'Intercept':
# Independent variables are orthogonal
ssr1, df_resid1 = _ssr_reduced_model(
y, x, term_slices, params, [key])
df1 = df_resid1 - df_resid
msm = (ssr1 - ssr) / df1
if (key == ':'.join(factors[:-1]) or
(key + ':' + subject not in term_slices)):
# Interaction effect ### Edit 04/10/2018
mse = ssr / df_resid
df2 = df_resid
type_3_SS_error = ssr # Edit 04/10/2018
else:
# Simple main effect
ssr1, df_resid1 = _ssr_reduced_model(
y, x, term_slices, params,
[key + ':' + subject])
df2 = df_resid1 - df_resid
mse = (ssr1 - ssr) / df2
type_3_SS_error = ssr1 - ssr # Edit 04/10/2018
F = msm / mse
p = stats.f.sf(F, df1, df2)
term = key.replace('C(', '').replace(', Sum)', '')
anova_table.loc[term, 'Num DF'] = df1
anova_table.loc[term, 'Den DF'] = df2
anova_table.loc[term, 'Pr > F'] = p
# ---------------- Edit 04/10/2018 ----------------- #
anova_table.loc[term, 'T3 SS'] = msm
anova_table.loc[term, 'T3 SS (error)'] = type_3_SS_error
# calculate effect sizes.
ss_total = np.var(y) * (len(y) - 1)
omega_squared = (msm - (df1 * mse)) / (ss_total + mse) # omega squared
eta_squared = msm / ss_total
partial_eta_squared = msm / (msm + type_3_SS_error)
anova_table.loc[term, 'omega^2'] = omega_squared
# ---------------- Edit 04/10/2018 ---------------- #
return AnovaResults(anova_table.iloc[:, [4, 5, 1, 2, 0, 3, 6, 7, 8]])
これをセル上で実行すると、CustomAnovaRM()という関数が使えるようになる。
あとはAnovaRM()の代わりにCustomAnovaRM()を実行するだけ。
# モデルを作成する
data = df[['SubjectNumber', 'r_1_Score', 'r_2_Score', 'r_3_Score', 'r_4_Score']]
# データフレームを長い形式に変換する
data_melted = data.melt(id_vars=['SubjectNumber'], value_vars=['r_1_Score', 'r_2_Score', 'r_3_Score', 'r_4_Score'], var_name='Group', value_name='Score')
# Repeated Measures ANOVAを実行する
aovrm = CustomAnovaRM(data_melted, 'Score', 'SubjectNumber', within=['Group'])
res = aovrm.fit()
print(res.summary())
print(res.anova_table)
あとはAnovaRM()の代わりにCustomAnovaRM()を実行するだけ。
Before
計算方法自体の正確性はチェックしていないのでご注意ください。
