Pythonで2×2の二要因分散分析をやってみた【心理統計】

こんにちは、ヒロムです。

前回は「Pythonで一要因分散分析をやってみた」ので、今回は2要因のケースで試してみます。

一要因分散分析は1つの要因が3水準以上ある場合に行いました。多くの場合、グラフで表すと棒グラフとなります。

今回扱う二要因分散分析では、2つの要因が2水準以上ある場合に行います。多くの場合、グラフで表すと折れ線グラフとなります。

では、早速やっていきましょう。

Pythonで二要因分散分析をやってみる

例題

「ただしイケメンに限る」

このセリフはネット界ではよく見かけるワードだ。最近聞いた例では「頭ナデナデ」である。女性は男性に頭をナデナデされるとキュンとするというが、まさにあれは「イケメンに限る」行為である。ブサメンにやられたら顔が引きつること間違いなし (偏見)。もう1つ例を挙げよう。次は「医者」である。というのは、医者という職業は多くの人が敬う職業で、それなりに給料も良い。医者の男性と付き合いたいという女性は多数派であろう。

ここで、医者の容姿について考えてみる。たとえ同じ医者でも、イケメンとブサメンではその魅力度は異なる可能性がある。つまり、「お医者さんってカッコイイ〜」というが「ただしイケメンに限る」話かもしれない。これについて、容姿(イケメン・ブサメン)、職業(医者・非医者)を要因とし、調査した。調査協力者は20名であり、協力者には「イケメン・医者 / イケメン・非医者 / ブサメン・医者 / ブサメン・非医者」の4通りにおいて、それぞれの魅力度を1点を最低点、5点を最高点とし評価させた。

 

情報整理

・要因:容姿(イケメン・ブサメンの2水準)、職業(医者・非医者の2水準)

・計画:被験者内計画

・被験者数:20名 (女性20名)

 

必要最低限の情報は以上の通りですかね。

なお、サンプルデータは独自で作成したものです。実際に調査したデータではありませんので、結果は鵜呑みにしないように。

 

二要因分散分析を実施する

まずはデータを読み込みます。

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

data = pd.read_csv("2wayAnova.csv", encoding="shift-jis")
data.head(10)

表の並べ方は心理統計で見るようなものとは異なりますが、Pythonで用いる場合はこのスタイルが適切です。

次に、分布を視覚的に把握するために折れ線グラフを描いてみましょう。

from statsmodels.graphics.factorplots import interaction_plot

fig = interaction_plot(data.job, data.appearance, data.score, colors = ['red','blue'], markers=['D','^'], ms=10)

グラフを見る限り、職業の主効果、容姿の主効果、職業×容姿の交互作用がありそうですな。

ということで、被験者内計画による二要因分散分析をやってみましょう。

Pythonで分散分析 (対応なし・二元配置)」の記事を参考にしました。

from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm

formula = 'score ~ (job) + (appearance) + (job):(appearance)'
model = ols(formula, data).fit()
aov_table = anova_lm(model, typ=2)
print(aov_table)

結果、職業の主効果、容姿の主効果、職業×容姿の交互作用が全て有意でした。

職業×容姿の交互作用が有意であったので、多重比較をやっていきます (交互作用が有意な場合、主効果について多重比較する必要はない)。

多重比較

まずは分解。

coolD = data[(data.job == "doctor") & (data.appearance == "cool")] 
coolN = data[(data.job == "non_doctor")&(data.appearance == "cool")]
uglyD = data[(data.job == "doctor") & (data.appearance == "ugly")] 
uglyN = data[(data.job == "non_doctor")&(data.appearance == "ugly")]
coolD.head(10) # 中身の確認
A = coolD["score"]
B = coolN["score"]
C = uglyD["score"]
D = uglyN["score"]

# 関数(グループ数の羅列,グループ名)
def tukey_hsd( ind, *args ):
    from statsmodels.stats.multicomp import pairwise_tukeyhsd
    data_arr = np.hstack( args )
    ind_arr = np.array([])
    for x in range(len(args)):
        ind_arr = np.append(ind_arr, np.repeat(ind[x], len(args[x])))
    print(pairwise_tukeyhsd(data_arr,ind_arr))

tukey_hsd(list('ABCD'), A, B, C, D)

(因みに、各条件のscoreをわざわざ変数化させましたが、こうしなかったら最後のtukey_hsd(list(‘ABCD’), A, B, C, D)のところでなぜかエラーが吐かれました。)

結果、CーD間以外でそれぞれ有意差があったようですね。

結果の解釈

さて、グラフを見れば一目瞭然なのですが、以下のことがわかりました。

・職業に関わらず、ブサメンよりもイケメンの方が魅力度が高い

・イケメン×医者の組み合わせが最強。次点でイケメン×非医者

結論!やっぱり「ただしイケメンに限る」でしたね。

以上。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です