特別講義DS Ch16 カテゴリーデータ処理

import pandas as pd
import matplotlib.pyplot as plt
import matplotlib_fontja
import seaborn as sns
import scipy.stats as st
import numpy as np

#グラフの保存場所
dir_fig = './result/fig/'

#グラフの設定
plt.style.use('ggplot') #グラフスタイル
plt.rcParams['figure.figsize'] = [12, 9] #グラフサイズ
plt.rcParams['font.size'] = 14 #フォントサイズ

#------------------------------------------------------------------
#データの読み込み
#------------------------------------------------------------------

df = pd.read_csv('./data/questionnaire_sample_original.csv')
print(df.columns)

"""
Index(['あなたの「スマートフォン」に対する印象について、どれくらいあてはまるかを教えてください。',
       'あなたの「スマートフォン」に対する操作や知識について、どれくらい自信があるかを教えてください。',
       'あなたの「パソコン」に対する印象について、どれくらいあてはまるかを教えてください。',
       'あなたの「パソコン」に対する操作や知識について、どれくらい自信があるかを教えてください。',
       'あなたの「SNS（FacebookやInstagramなど）」に対する印象について、どれくらいあてはまるかを教えてください。',
       'あなたの「SNS（FacebookやInstagramなど）」に対する操作や知識について、どれくらい自信があるかを教えてください。',
       'あなたの「ネットを使ったコミュニケーション（LINEやTwitterなど）」に対する印象について、どれくらいあてはまるかを教えてください。',
       'あなたの「ネットを使ったコミュニケーション（LINEやTwitterなど）」に対する操作や知識について、どれくらい自信があるかを教えてください。',
       '次の「IT用語」について「知っている」という言葉にチェックをしてください（知っているものすべてにチェック）。',
       'Windowsの操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。',
       '電子メールの操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。',
       'WWW(ウェブ)の操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。',
       'Wordの操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。',
       'Excelの操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。',
       'PowerPointの操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。',
       'Officeの応用的な操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。',
       '情報倫理について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。'],
      dtype='object')
"""

new_columns = dict([(t,'Q'+str(n + 1))for (n,t) in  enumerate(df.columns)])

"""
一つ一つ手書きすると以下のようになる.
new_columns = {'あなたの「スマートフォン」に対する印象について、どれくらいあてはまるかを教えてください。':'Q1'
              ,'あなたの「スマートフォン」に対する操作や知識について、どれくらい自信があるかを教えてください。':'Q2'
              ,'あなたの「パソコン」に対する印象について、どれくらいあてはまるかを教えてください。':'Q3'
              ,'あなたの「パソコン」に対する操作や知識について、どれくらい自信があるかを教えてください。':'Q4'
              ,'あなたの「SNS（FacebookやInstagramなど）」に対する印象について、どれくらいあてはまるかを教えてください。':'Q5'
              ,'あなたの「SNS（FacebookやInstagramなど）」に対する操作や知識について、どれくらい自信があるかを教えてください。':'Q6'
              ,'あなたの「ネットを使ったコミュニケーション（LINEやTwitterなど）」に対する印象について、どれくらいあてはまるかを教えてください。':'Q7'
              ,'あなたの「ネットを使ったコミュニケーション（LINEやTwitterなど）」に対する操作や知識について、どれくらい自信があるかを教えてください。':'Q8'
              ,'次の「IT用語」について「知っている」という言葉にチェックをしてください（知っているものすべてにチェック）。':'Q9'
              ,'Windowsの操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。':'Q10'
              ,'電子メールの操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。':'Q11'
              ,'WWW(ウェブ)の操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。':'Q12'
              ,'Wordの操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。':'Q13'
              ,'Excelの操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。':'Q14'
              ,'PowerPointの操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。':'Q15'
              ,'Officeの応用的な操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。':'Q16'
              ,'情報倫理について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。':'Q17'}
"""
print(new_columns)

"""
{'あなたの「スマートフォン」に対する印象について、どれくらいあてはまるかを教えてください。': 'Q1'
, 'あなたの「スマートフォン」に対する操作や知識について、どれくらい自信があるかを教えてください。': 'Q2'
, 'あなたの「パソコン」に対する印象について、どれくらいあてはまるかを教えてください。': 'Q3'
, 'あなたの「パソコン」に対する操作や知識について、どれくらい自信があるかを教えてください。': 'Q4'
, 'あなたの「SNS（FacebookやInstagramなど）」に対する印象について、どれくらいあてはまるかを教えてください。': 'Q5'
, 'あなたの「SNS（FacebookやInstagramなど）」に対する操作や知識について、どれくらい自信があるかを教えてください。': 'Q6'
, 'あなたの「ネットを使ったコミュニケーション（LINEやTwitterなど）」に対する印象について、どれくらいあてはまるかを教えてください。': 'Q7'
, 'あなたの「ネットを使ったコミュニケーション（LINEやTwitterなど）」に対する操作や知識について、どれくらい自信があるかを教えてください。': 'Q8'
, '次の「IT用語」について「知っている」という言葉にチェックをしてください（知っているものすべてにチェック）。': 'Q9'
, 'Windowsの操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。': 'Q10'
, '電子メールの操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。': 'Q11'
, 'WWW(ウェブ)の操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。': 'Q12'
, 'Wordの操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。': 'Q13'
, 'Excelの操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。': 'Q14'
, 'PowerPointの操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。': 'Q15'
, 'Officeの応用的な操作や事柄について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。': 'Q16'
, '情報倫理について「身に付いている」ものにチェックをしてください（身に付いているものすべてにチェック）。': 'Q17'}
"""

df.rename(new_columns, inplace=True,axis='columns')
print(df)
"""
             Q1         Q2         Q3  ...                                                Q15                         Q16                                                Q17
0     まあまあ興味がある    どちらでもない    どちらでもない  ...           スライドマスターを活用することができる;SmartArtを挿入することができる;  学内PCで作成したファイルを自宅のPCに持ち運べる;                                       どれも身についていない;
1      すごく興味がある  まあまあ自信がある   すごく興味がある  ...                                新しいスライドを追加することができる;    学術的文章の規則を踏まえてレポートを作成できる;  学内PCでログアウトせずに帰る危険性を理解している;マルウェア（不正ソフト）について理解して...
2      すごく興味がある   あまり自信がない  まあまあ興味がある  ...                                       どれも身についていない;                どれも身に付いていない;                                       どれも身についていない;
3       どちらでもない   あまり自信がない    全く興味がない  ...                                新しいスライドを追加することができる;                どれも身に付いていない;                         学内PCでログアウトせずに帰る危険性を理解している;
4     まあまあ興味がある  まあまあ自信がある   すごく興味がある  ...                                新しいスライドを追加することができる;                どれも身に付いていない;                                       どれも身についていない;
...         ...        ...        ...  ...                                                ...                         ...                                                ...
1347  まあまあ興味がある   あまり自信がない  まあまあ興味がある  ...  新しいスライドを追加することができる;スライドにページ番号を挿入することができる;スライドマ...                どれも身に付いていない;  学内PCでログアウトせずに帰る危険性を理解している;フィッシング詐欺について理解している;ウ...
1348  まあまあ興味がある   あまり自信がない  まあまあ興味がある  ...                                新しいスライドを追加することができる;                どれも身に付いていない;                         学内PCでログアウトせずに帰る危険性を理解している;
1349  まあまあ興味がある    全く自信がない    全く興味がない  ...                                新しいスライドを追加することができる;                どれも身に付いていない;                                       どれも身についていない;
1350  まあまあ興味がある  まあまあ自信がある  まあまあ興味がある  ...             アニメーション機能を使うことができる;新しいスライドを追加することができる;                どれも身に付いていない;  学内PCでログアウトせずに帰る危険性を理解している;ウェブ上の画像に関する著作権について理解...
1351  まあまあ興味がある   あまり自信がない  まあまあ興味がある  ...          新しいスライドを追加することができる;スライドにページ番号を挿入することができる;                どれも身に付いていない;  フィッシング詐欺について理解している;ウェブ上の画像に関する著作権について理解している;学内...
"""

kyomi = {'すごく興味がある':5
        ,'まあまあ興味がある':4
        ,'どちらでもない':3
        ,'あまり興味がない':2
        ,'全く興味がない':1}

zishin = {'すごく自信がある':5
         ,'まあまあ自信がある':4
         ,'どちらでもない':3
         ,'あまり自信がない':2
         ,'全く自信がない':1}

for i in [1,3,5,7]:
    df['Q' + str(i) + '_num'] = df['Q' + str(i)].apply(lambda x : kyomi[x])

for i in [2,4,6,8]:
    df['Q' + str(i) + '_num'] = df['Q' + str(i)].apply(lambda x : zishin[x])

#｢身についているもの｣の数の列を作る
#------------------------------------------------------------------

def to_list(x):
    # ';'で区切ったリストにする
    # 'どれも' で始まる回答 ('どれも身についていない' など) がある場合は空のリストを返す
    xs = x.split(';')[0:-1]
    if any(t.startswith('どれも') for t in xs):
        return []
    else:
        return xs

#テスト
print(df.at[100,'Q13'])
print(to_list(df.at[100,'Q13']))
# >>> フォントの種類やサイズを変更することができる;中央揃えや右揃えを設定することができる;
# >>> ['フォントの種類やサイズを変更することができる', '中央揃えや右揃えを設定することができる']

for i in range(9,18,1):
    df['Q' + str(i) + '_num'] = df['Q' + str(i)].apply(lambda x : len(to_list(x)))

print(df['Q17_num'])
"""
0       0
1       3
2       0
3       1
4       0
       ..
1347    3
1348    1
1349    0
1350    2
1351    3
Name: Q17_num, Length: 1352, dtype: int64
"""
num_columns = ['Q' + str(x) + '_num' for x in range(1,18,1)]
print(num_columns)
"""
['Q1_num', 'Q2_num', 'Q3_num', 'Q4_num', 'Q5_num', 'Q6_num', 'Q7_num', 'Q8_num', 'Q9_num', 'Q10_num', 'Q11_num', 'Q12_num', 'Q13_num', 'Q14_num', 'Q15_num', 'Q16_num', 'Q17_num']
"""

#ヒストグラムを作る
#------------------------------------------------------------------
## Q1~Q6
fig, ax = plt.subplots(nrows = 6
                      ,ncols = 1
                      ,tight_layout=True)
for i in range(6):
    col = 'Q' + str(i + 1) + '_num'
    xs = df[col].value_counts().sort_index()
    ax[i].bar(xs.index,xs.values, width = 1)
    ax[i].set_title(col)

plt.savefig(dir_fig + 'Q1_6_histgram.png')
plt.close(fig)

##Q7~Q12
fig, ax = plt.subplots(nrows = 6
                      ,ncols = 1
                      ,tight_layout=True)
for i in range(6):
    col = 'Q' + str(i + 7) + '_num'
    xs = df[col].value_counts().sort_index()
    ax[i].bar(xs.index,xs.values, width = 1)
    ax[i].set_title(col)

plt.savefig(dir_fig + 'Q7_12_histgram.png')
plt.close(fig)

##Q13~Q17
fig, ax = plt.subplots(nrows = 5
                      ,ncols = 1
                      ,tight_layout=True)
for i in range(5):
    col = 'Q' + str(i + 13) + '_num'
    xs = df[col].value_counts().sort_index()
    ax[i].bar(xs.index,xs.values, width = 1)
    ax[i].set_title(col)

plt.savefig(dir_fig + 'Q13_17_histgram.png')
plt.close(fig)

#順位相関係数を求める
#------------------------------------------------------------------
sns.heatmap(df[num_columns].corr(method='spearman')
           ,cmap=plt.get_cmap('bwr')
           ,vmax=1.0
           ,center =0
           ,vmin=-1.0
           ,xticklabels=1
           ,yticklabels=1
           ,linewidths=.5
           ,annot=True)

plt.title('spearman')
plt.xticks(rotation=45)
plt.savefig(dir_fig + 'spearman.png')
plt.close()

#クロス表
#------------------------------------------------------------------
#質問項目の区分間にどのような関係があるのかを検定する
#今回はスマホの自信(Q2)と,PCの自信(Q4)だけ(自分の仮説に応じて色々ためす)
# 作業手順
## クロス表を作成する
## 独立性の検定(お互いに影響があるか)を検定する
## 検定が有意であれば,どの項目が具体的にどの項目に影響しているかを残差分析で明らかにする.

#クロス表の作成
## 列と行に指定する質問項目
target_row = 'Q4'
target_col = 'Q2'
## 選択肢の順序(指定しないとデータの出てきた順になる)
order = {target_row:['すごく自信がある'
                    ,'まあまあ自信がある'
                    ,'どちらでもない'
                    ,'あまり自信がない'
                    ,'全く自信がない']
        ,target_col:['すごく自信がある'
                    ,'まあまあ自信がある'
                    ,'どちらでもない'
                    ,'あまり自信がない'
                    ,'全く自信がない']}
## 図の軸に表示する質問の内容
labels = {'Q2':'スマートフォンへの自信 (Q2)'
         ,'Q4':'パソコンへの自信 (Q4)'}

crosstab = pd.crosstab(df[target_row],df[target_col])
crosstab = crosstab.reindex(order[target_row],axis='index')
crosstab = crosstab.reindex(order[target_col],axis='columns')

sns.heatmap( crosstab
           , cmap=plt.get_cmap('Reds')
           , yticklabels=1
           , xticklabels=1
           , linewidths=.5
           , annot=True
           , fmt='d')

plt.title('同時度数分布表 (人数)')
plt.ylabel(labels[target_row])
plt.xlabel(labels[target_col])
plt.xticks(rotation=30)
plt.tight_layout()
plt.savefig(dir_fig + 'cross_table_'+ target_row + '_' + target_col +'.png')
plt.close()

#独立性の検定
x2, p, dof, e = st.chi2_contingency(crosstab,correction=False)

print(p)
# >>> 2.2871064452593352e-84

if p <= 0.05:
    print('有意')
else:
    print('有意でない')
# >>> 有意

#有意であれば,残差分析を行う
##期待値表の計算
expect = pd.DataFrame(e)
expect.columns = crosstab.columns
expect.index   = crosstab.index

#残差の計算
residual = crosstab - expect

#標準化残差の計算
col_sum = crosstab.sum()
row_sum = crosstab.transpose().sum()
total   = col_sum.sum()
stdres  = residual.copy()

for c in stdres.columns:
    for r in stdres.index:
        stdres.at[r,c]    = stdres.at[r,c] \
                         / np.sqrt( expect.at[r,c]
                                  * (1 - col_sum[c] / total)
                                  * (1 - row_sum[r] / total))

#ヒートマップの作成
## 色の範囲は絶対値の最大に合わせ, 絶対値が 1.96 を超えるセルに * を付ける
lim   = np.abs(stdres.values).max()
annot = stdres.apply(lambda col: col.map(lambda v: f'{v:.1f}' + ('*' if abs(v) > 1.96 else '')))
sns.heatmap(stdres
           ,cmap=plt.get_cmap('RdBu_r')
           ,vmax=lim
           ,vmin=-lim
           ,center=0
           ,xticklabels=1
           ,linewidths=.5
           ,annot=annot
           ,fmt='')
plt.title('標準化残差 (* は絶対値が 1.96 を超えるセル)')
plt.ylabel(labels[target_row])
plt.xlabel(labels[target_col])
plt.xticks(rotation=30)
plt.tight_layout()
plt.savefig(dir_fig + 'stdres_' + target_row + '_' + target_col + '.png')
plt.close()

#次元削減
#------------------------------------------------------------------
#興味,自信,知識を統合して3変数にする

interest    = df[['Q' + str(x) + '_num' for x in [1,3,5,7]]]
confidence  = df[['Q' + str(x) + '_num' for x in [2,4,6,8]]]
knowledge   =  df[['Q' + str(x) + '_num' for x in range(9,18,1)]]

"""
#普通に平均値を出して,正規化するパターン
def min_max_scaling(x):
    return (x - x.min())/(x.max() - x.min())
interest = min_max_scaling(interest.mean(axis=1))
confidence = min_max_scaling(confidence.mean(axis=1))
knowledge = min_max_scaling(knowledge.mean(axis=1))

print(interest)
"""
#PCAやt-sneでやってみるパターン(ただし,解釈が難しい)
from sklearn.decomposition  import PCA

pca = PCA()

pca.fit(interest)
interest = pca.transform(interest)[:,0]
print('interest:', pca.explained_variance_ratio_[0])

pca.fit(confidence)
confidence = pca.transform(confidence)[:,0]
print('confidence:', pca.explained_variance_ratio_[0])

pca.fit(knowledge)
knowledge = pca.transform(knowledge)[:,0]
print('knowledge:', pca.explained_variance_ratio_[0])


#図示してみる

## interest-confidence
plt.scatter(y=interest,x=confidence)
plt.title('interest-confidence')
plt.ylabel('interest')
plt.xlabel('confidence')
plt.savefig(dir_fig + 'interest-confidence.png')
plt.close()

## interest-knowledge
plt.scatter(y=interest,x=knowledge)
plt.title('interest-knowledge')
plt.ylabel('interest')
plt.xlabel('knowledge')
plt.savefig(dir_fig + 'interest-knowledge.png')
plt.close()

## confidence-knowledge
plt.scatter(y=confidence,x=knowledge)
plt.title('confidence-knowledge')
plt.ylabel('confidence')
plt.xlabel('knowledge')
plt.savefig(dir_fig + 'confidence-knowledge.png')
plt.close()

##3変数
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(interest,confidence,knowledge)
plt.savefig(dir_fig + 'pca3d.png')
plt.close(fig)

#それぞれの相関係数を確認する
df['interest']      = interest
df['confidence']    = confidence
df['knowledge']     = knowledge

sns.heatmap(df[['interest','confidence','knowledge']].corr(method='spearman')
           ,cmap=plt.get_cmap('bwr')
           ,vmax=1.0
           ,center =0
           ,vmin=-1.0
           ,xticklabels=1
           ,yticklabels=1
           ,linewidths=.5
           ,annot=True)

plt.title('pca_spearman')
plt.xticks(rotation=45)
plt.savefig(dir_fig + 'pca_spearman.png')
plt.close()

#このあとクラスタリングしてみることで,学生の分類ができるかもしれない
#(けど今回はやらない)

import numpy as np

# 1 因子 4 指標のモデルのパラメータ
lam = np.array([1.0, 0.8, 0.6, 0.7])  # 因子負荷量 (1 つ目は 1 に固定)
phi = 0.5                              # 因子の分散
psi = np.array([0.3, 0.4, 0.5, 0.4])   # 独自因子の分散

# モデルが予言する共分散行列
Sigma = phi * np.outer(lam, lam) + np.diag(psi)
print(Sigma.round(3))
"""
[[0.8   0.4   0.3   0.35 ]
 [0.4   0.72  0.24  0.28 ]
 [0.3   0.24  0.68  0.21 ]
 [0.35  0.28  0.21  0.645]]
"""

# このモデルに従う仮想の回答者を 1000 人つくり, 共分散行列を計算する
rng = np.random.default_rng(0)
n = 1000
f = rng.normal(0, np.sqrt(phi), n)            # 因子の値 (観測できない)
eps = rng.normal(0, np.sqrt(psi), (n, 4))     # 独自因子
X = f[:, None] * lam + eps                    # 観測される 4 項目
print(np.cov(X, rowvar=False).round(3))
"""
[[0.735 0.348 0.288 0.324]
 [0.348 0.694 0.25  0.282]
 [0.288 0.25  0.683 0.195]
 [0.324 0.282 0.195 0.66 ]]
"""

#共分散構造分析
#------------------------------------------------------------------
#変数間の関係の仮説をモデルにして, データと矛盾しないかを確かめる.
import semopy
from semopy import Model

#データ型の変換 (すべて float でないといけない)
df_sem = df[num_columns].astype('float64')

#標準化する
df_sem = (df_sem - df_sem.mean()) / df_sem.std(ddof=0)

#わかりにくいから名前を変える
df_sem.rename({'Q2_num':'Smartphone'
              ,'Q4_num':'PC'
              ,'Q13_num':'Word'
              ,'Q14_num':'Excel'
              ,'Q15_num':'PowerPoint'
              ,'Q10_num':'Windows'
              ,'Q17_num':'Info'
              ,'Q9_num':'IT'
              ,'Q8_num':'Communication'
              ,'Q6_num':'SNS'}
              ,inplace=True
              ,axis='columns')

# 仮説 (文字列で指定する)
desc = '''
    # (測定モデル) 潜在変数 =~ 観測変数
    creation =~ Word + Excel + PowerPoint + Windows
    knowledge =~ Info + IT
    consumption =~ Communication + SNS

    # (構造モデル) 目的変数 ~ 説明変数
    PC ~ creation + knowledge + consumption
    Smartphone ~ creation + knowledge + consumption

    # (誤差の共分散) 変数 ~~ 変数
    PC ~~ Smartphone
    '''

# モデルを用意して推定する
mod = Model(desc)
res = mod.fit(df_sem)
print(res)
"""
Name of objective: MLW
Optimization method: SLSQP
Optimization successful.
Optimization terminated successfully
Objective value: 0.254
Number of iterations: 43
Params: 0.959 0.808 0.907 1.276 1.054 0.244 0.495 0.273 -0.102 0.341 0.800 0.144 0.645 0.365 0.290 0.371 0.614 0.496 0.295 0.525 0.364 0.228 0.635 0.134 0.772 0.152 0.462 0.386
"""

# 推定されたパラメータの一覧 (std_est=True で標準化推定値の列が付く)
estimates = mod.inspect(std_est=True)
print(estimates[['lval', 'op', 'rval', 'Estimate', 'Est. Std', 'p-value']])
"""
             lval  op           rval  Estimate  Est. Std   p-value
0            Word   ~       creation  1.000000  0.878675         -
1           Excel   ~       creation  0.959228  0.842957       0.0
2      PowerPoint   ~       creation  0.808264  0.710333       0.0
3         Windows   ~       creation  0.907300  0.797326       0.0
4            Info   ~      knowledge  1.000000  0.621641         -
5              IT   ~      knowledge  1.276020  0.793006       0.0
6   Communication   ~    consumption  1.000000  0.796807         -
7             SNS   ~    consumption  1.054072  0.839821       0.0
8              PC   ~       creation  0.243879  0.214321  0.002827
9              PC   ~      knowledge  0.494677  0.307485  0.000087
10             PC   ~    consumption  0.272566  0.217144       0.0
11     Smartphone   ~       creation -0.101796 -0.089464  0.182332
12     Smartphone   ~      knowledge  0.340858  0.211889  0.003632
13     Smartphone   ~    consumption  0.799502  0.636981       0.0
14    consumption  ~~    consumption  0.634833  1.000000       0.0
15    consumption  ~~      knowledge  0.133671  0.269868       0.0
16       creation  ~~       creation  0.772483  1.000000       0.0
17       creation  ~~    consumption  0.151833  0.216817       0.0
18       creation  ~~      knowledge  0.462311  0.846120       0.0
19      knowledge  ~~      knowledge  0.386469  1.000000       0.0
20             PC  ~~     Smartphone  0.143877  0.143851       0.0
21             PC  ~~             PC  0.644794  0.644631       0.0
22  Communication  ~~  Communication  0.365059  0.365098       0.0
23          Excel  ~~          Excel  0.289505  0.289424       0.0
24             IT  ~~             IT  0.371377  0.371141       0.0
25           Info  ~~           Info  0.613612  0.613562       0.0
26     PowerPoint  ~~     PowerPoint  0.495509  0.495427       0.0
27            SNS  ~~            SNS  0.294719  0.294701       0.0
28     Smartphone  ~~     Smartphone  0.525353  0.525298       0.0
29        Windows  ~~        Windows  0.364371  0.364271       0.0
30           Word  ~~           Word  0.228052  0.227930       0.0
"""

# 適合度の指標 (本章で説明した指標だけを表示する)
cols = ['DoF', 'chi2', 'chi2 p-value', 'CFI', 'TLI', 'RMSEA']
fit_stats = semopy.calc_stats(mod)
print(fit_stats.T.loc[cols])
"""
                   Value
DoF            27.000000
chi2          343.794134
chi2 p-value    0.000000
CFI             0.951725
TLI             0.919542
RMSEA           0.093192
"""

#パス図を描く (標準化推定値を表示する)
#plot_covs=True で共分散も表示する
semopy.semplot(mod
              ,dir_fig + 'SEM.png'
              ,plot_covs=True
              ,std_ests=True
              ,engine='dot')

#残差を確かめる
#------------------------------------------------------------------
#データの共分散行列とモデルが予言した共分散行列の差を, 相関の尺度で見る
import numpy as np

obs = mod.vars['observed']                               # モデルに含まれる観測変数
S = np.cov(df_sem[obs].values, rowvar=False, bias=True)  # データの共分散行列
Sigma = mod.calc_sigma()[0]                              # モデルが予言した共分散行列
sd = np.sqrt(np.diag(S))
resid = pd.DataFrame((S - Sigma) / np.outer(sd, sd), index=obs, columns=obs)
print(resid.round(2))
"""
               Communication  Excel    IT  Info    PC  PowerPoint   SNS  Smartphone  Windows  Word
Communication           0.00  -0.03  0.03  0.07  0.02        0.06  0.00       -0.00     0.09  0.02
Excel                  -0.03  -0.00  0.01 -0.04 -0.02       -0.03 -0.08       -0.06    -0.04  0.03
IT                      0.03   0.01 -0.00 -0.00  0.01       -0.07 -0.05       -0.00     0.09 -0.04
Info                    0.07  -0.04 -0.00 -0.00 -0.03        0.06  0.01        0.01     0.08 -0.03
PC                      0.02  -0.02  0.01 -0.03 -0.00       -0.05 -0.01       -0.00     0.11 -0.03
PowerPoint              0.06  -0.03 -0.07  0.06 -0.05       -0.00  0.02        0.01     0.02  0.02
SNS                     0.00  -0.08 -0.05  0.01 -0.01        0.02 -0.00        0.00     0.03 -0.01
Smartphone             -0.00  -0.06 -0.00  0.01 -0.00        0.01  0.00       -0.00     0.09 -0.01
Windows                 0.09  -0.04  0.09  0.08  0.11        0.02  0.03        0.09    -0.00 -0.03
Word                    0.02   0.03 -0.04 -0.03 -0.03        0.02 -0.01       -0.01    -0.03 -0.00
"""

#SRMR (残差の二乗平均の平方根. 対角と下三角の要素を使う)
srmr = np.sqrt((resid.values[np.tril_indices(len(obs))] ** 2).mean())
print('SRMR:', srmr)
"""
SRMR: 0.04042660130885509
"""

#修正したモデル
#------------------------------------------------------------------
desc2 = '''
    creation =~ Word + Excel + PowerPoint + Windows
    knowledge =~ Info + IT + Windows
    consumption =~ Communication + SNS

    PC ~ creation + knowledge + consumption
    Smartphone ~ creation + knowledge + consumption

    PC ~~ Smartphone
    '''
mod2 = Model(desc2)
mod2.fit(df_sem)

estimates2 = mod2.inspect(std_est=True)
print(estimates2[['lval', 'op', 'rval', 'Estimate', 'Est. Std', 'p-value']])
"""
             lval  op           rval  Estimate  Est. Std   p-value
0            Word   ~       creation  1.000000  0.902817         -
1           Excel   ~       creation  0.943159  0.851513       0.0
2      PowerPoint   ~       creation  0.787335  0.710805       0.0
3         Windows   ~       creation  0.226744  0.204725  0.000354
4         Windows   ~      knowledge  1.106513  0.670206       0.0
5            Info   ~      knowledge  1.000000  0.605709         -
6              IT   ~      knowledge  1.255615  0.760449       0.0
7   Communication   ~    consumption  1.000000  0.797009         -
8             SNS   ~    consumption  1.053600  0.839756       0.0
9              PC   ~       creation  0.001609  0.001452  0.981641
10             PC   ~      knowledge  0.890435  0.539311       0.0
11             PC   ~    consumption  0.239828  0.191169       0.0
12     Smartphone   ~       creation -0.168782 -0.152384  0.007119
13     Smartphone   ~      knowledge  0.458509  0.277699  0.000006
14     Smartphone   ~    consumption  0.783506  0.624525       0.0
15    consumption  ~~    consumption  0.635361  1.000000       0.0
16       creation  ~~       creation  0.815135  1.000000       0.0
17       creation  ~~    consumption  0.139921  0.194428       0.0
18      knowledge  ~~      knowledge  0.366825  1.000000       0.0
19      knowledge  ~~    consumption  0.138980  0.287881       0.0
20      knowledge  ~~       creation  0.448629  0.820432       0.0
21             PC  ~~     Smartphone  0.125249  0.125250       0.0
22             PC  ~~             PC  0.611823  0.611842       0.0
23  Communication  ~~  Communication  0.364856  0.364777       0.0
24          Excel  ~~          Excel  0.274935  0.274925       0.0
25             IT  ~~             IT  0.421749  0.421717       0.0
26           Info  ~~           Info  0.633018  0.633117       0.0
27     PowerPoint  ~~     PowerPoint  0.494811  0.494757       0.0
28            SNS  ~~            SNS  0.294856  0.294811       0.0
29     Smartphone  ~~     Smartphone  0.516227  0.516219       0.0
30        Windows  ~~        Windows  0.283743  0.283772       0.0
31           Word  ~~           Word  0.184933  0.184921       0.0
"""
print(semopy.calc_stats(mod2).T.loc[cols])
"""
                   Value
DoF            26.000000
chi2          174.089440
chi2 p-value    0.000000
CFI             0.977433
TLI             0.960942
RMSEA           0.064930
"""

semopy.semplot(mod2
              ,dir_fig + 'SEM2.png'
              ,plot_covs=True
              ,std_ests=True
              ,engine='dot')