PISA 2022データを読む

はじめに

PISA 2018データを読むほぼそのままである。

2023年12月5日、PISA 2022 resultsが公開された。データはPISA 2022 Databaseからダウンロードできる。ここではSPSS用データの生徒アンケート（STU_QQQ_SPSS.zipを展開して得られるCY08MSP_STU_QQQ.savという2Gバイトほどのファイル）をPythonで読み込んでみる。

SPSSのsavデータを読むにはpandasの read_spss() でもできるが，これはpyreadstatパッケージを呼び出しているだけである。以下ではpyreadstatを直接使う（pip install pyreadstat しておく）：

import pyreadstat

全部一度に読む（これがエラーになるなら次に進んでください）

メモリの多いマシンなら全部読み込んでしまうのが手っ取り早い：

df, meta = pyreadstat.read_sav('CY08MSP_STU_QQQ.sav')

dfは

import sys
sys.getsizeof(df)

とすればわかるように，約7Gのメモリを占める。

df.shape

と打ち込めば613744行1279列であることがわかる。

これ以降は次と同じである。

一部分のデータだけ読み込む（メモリの少ないマシン用）

まずはメタデータだけ読んでみよう：

df, meta = pyreadstat.read_sav('CY08MSP_STU_QQQ.sav', metadataonly=True)

メタデータに含まれる変数名とその説明を表示する（長い出力に注意！）：

meta.column_names_to_labels

出力は次のようになる：

{'CNT': 'Country code 3-character',
 'CNTRYID': 'Country Identifier',
 'CNTSCHID': 'Intl. School ID',
 'CNTSTUID': 'Intl. Student ID',
・・・（以下略）・・・
}

これでは見にくいので，JSON形式で保存しておくことにする：

import json

with open("column_names_to_labels.json", "w") as f:
    json.dump(meta.column_names_to_labels, f, ensure_ascii=False, indent=2)

with open("variable_value_labels.json", "w") as f:
    json.dump(meta.variable_value_labels, f, ensure_ascii=False, indent=2)

これらはそれぞれ変数ラベルとその説明，変数ラベルとその値の説明である。例えば "IC182Q02JA" というラベルについては前者が

  "IC182Q02JA": "Agree/disagree: I am interested in learning [computer programming].",

後者が

  "IC182Q02JA": {
    "1.0": "Strongly disagree",
    "2.0": "Disagree",
    "3.0": "Agree",
    "4.0": "Strongly agree",
    "95.0": "Valid Skip",
    "97.0": "Not Applicable",
    "98.0": "Invalid",
    "99.0": "No Response"
  },

となっている。国名（3文字）とこの項目だけを改めて読み込んでみよう：

df, meta = pyreadstat.read_sav('CY08MSP_STU_QQQ.sav', usecols=['CNT', 'IC182Q02JA'])

"IC182Q02JA" の回答の度数分布を表示する：

df["IC182Q02JA"].value_counts()

結果は

IC182Q02JA
3.0    134464
2.0     93118
4.0     50305
1.0     44096
Name: count, dtype: int64

となって順番通りにならない。value_counts(sort=False) としてもランダムな順序になる。そこで

df[['CNT', 'IC182Q02JA']].groupby('IC182Q02JA').count()

とすると

               CNT
IC182Q02JA        
1.0          44096
2.0          93118
3.0         134464
4.0          50305

のように表らしくなってくれる。あるいは別解として

import pandas as pd

pd.crosstab(df["IC182Q02JA"], "count")  # "count" はダミー

とすると

col_0        count
IC182Q02JA        
1.0          44096
2.0          93118
3.0         134464
4.0          50305

になってくれるようだ。crosstab() の本来の使い方は

pd.crosstab(df["CNT"], df["IC182Q02JA"])

のようにして国名と回答の度数分布を

IC182Q02JA   1.0   2.0    3.0   4.0
CNT                                
ALB          207   640   1812   859
ARG          856  1697   3224  1299
AUS         2586  4376   3631  1002
AUT         1365  1569   1416   969
BEL         1180  2285   2200   759
BGR          391  1112   1986   704
・・・後略・・・

のように出力することである。日本に限れば

df["IC182Q02JA"][df["CNT"] == "JPN"].value_counts()
# または df[df["CNT"] == "JPN"].groupby('IC182Q02JA').count()

とすればよい。

国ごとに回答3＋回答4の割合を求めるには

s = pd.crosstab(df["CNT"], df["IC182Q02JA"]).apply(lambda x: (x.iloc[2]+x.iloc[3])/sum(x), axis=1)
s.sort_values(inplace=True)

手抜き棒グラフにする：

import matplotlib.pyplot as plt

plt.barh(s.index, s)
plt.savefig('231206a.svg', bbox_inches="tight")

国名を日本語にし，ドットプロットにしてみる：

co = pd.read_csv("https://okumuralab.org/~okumura/stat/data/countries.csv")
dic = { r['三字']: r['国名'] for i,r in co.iterrows() }
cj = [ dic[x] for x in s.index ]

plt.plot(s, cj, "o")
plt.grid()
plt.savefig('231206b.svg', bbox_inches="tight")

Last modified: