# 注意！Colab用です！PC上のRでは実行しないこと
# ダウンロード
system("curl -L -o library.tar.gz https://github.com/kawa5902/LSAdata/raw/refs/heads/main/202502library.tar.gz")
# 解凍
system("tar -xzf /content/library.tar.gz -C /content")
.libPaths("library")12 intsvyを使う
理論編ではPISAを例に，LSAで一般的に利用されている技術 （標本ウェイト・Replication Method・項目反応理論・推算値法）について 解説してきました。 これらの説明を聞いて面倒だと思った人も多いでしょう。 幸い，こうした処理はPISAやTIMSSを分析するだけならば， それほど大きな壁ではありません。 基本的な分析（平均や相関，あるいは回帰）であれば， Replication MethodsやPVsの統合を自動化してくれるソフトウェアがあるからです。 ここでは，RでPISAやTIMSSといったLSAを分析するためのパッケージである intsvyを紹介します。
12.1 intsvyのインストール
以下のコマンドを入力すると， colab上でintsvyが使えるようになります。
なお，この方法は一般的なRのパッケージのインストール方法ではありません。 通常のRでは，パッケージを追加する際install.packages関数を使います。 たとえばintsvyをインストールする際は，以下のように行います。
# インストール
install.packages("intsvy")Colab上でもinstall.packages関数を利用してパッケージをインストールすることは可能です。 ただColabを終了すると，インストールしたパッケージはすべて消えてしまいます。 そのため再びColabを起動した際は，あらためてinstall.packagesを実行する必要があります。 intsvyのインストールには10分程度必要なので，1〜2回ならともかく 分析のたびにinstall.packagesを実行するのは苦痛です。
そこで本書では，あらかじめColab上にインストールしたパッケージの塊（library.tar.gz）を 別の場所（github）1に保存しておき， それをダウンロードしてきて使うという方法を採用しています。 この方法なら1分弱でintsvyを利用可能になります。 ただし裏技に近い使い方なので，Colab以外では使用しないでください。
12.2 平均や分散の計算（PVを使わない場合）
ここではPISA2012の日本のデータを例に，いくつかintsvyを使った分析をしてみましょう。
# データの読み込み
url <- "https://raw.githubusercontent.com/kawa5902/LSAdata/refs/heads/main/pisa2012stuJPN.csv"
jpn2012 <- read.csv(url)最初に，ESCSの平均値を計算することを考えましょう。 まず，ESCSは欠測値として9999が設定されているので， 事前にNAに変換しておく必要があります。 続いてpisa.mean関数を使うことで，平均を求めることができます。 この関数は，標本ウェイトやreplication weightを考慮した計算を行ってくれるので， 引数を設定するだけで適切な推定値や標準誤差が算出されます。 なお，本書ではpisa.mean関数がintsvyパッケージに含まれていることを明示するため， intsvy::pisa.meanという具合に，関数の前にパッケージ名を示す記法を採用しています。 引数のvariableに与えるのは文字なので，ダブルクオーテーション（""）で囲むことを忘れないでください。
# 欠測をNAに変換
jpn2012$ESCS[jpn2012$ESCS == 9999] <- NA
intsvy::pisa.mean(variable = "ESCS", data = jpn2012) Freq Mean s.e. SD s.e
1 6185 -0.07151303 0.01524119 0.7036702 0.007597166度数分布はpisa.table関数です。 この関数は度数（Freq）と標本ウェイトを考慮した割合（Percentage）， および割合の標準誤差（Std.err）を算出します。 pisa.mean関数の場合と同じく，intsvy::pisa.tableとして利用します。
intsvy::pisa.table(variable = "ST03Q01", data = jpn2012) ST03Q01 Freq Percentage Std.err.
1 1 521 8.20 0.36
2 2 456 7.27 0.38
3 3 501 7.89 0.38
4 4 532 8.30 0.36
5 5 533 8.43 0.34
6 6 521 8.07 0.30
7 7 574 9.00 0.40
8 8 560 8.94 0.36
9 9 541 8.56 0.43
10 10 547 8.64 0.32
11 11 510 7.98 0.38
12 12 555 8.72 0.38いずれの関数も引数byを設定することで， 下位集団ごとの推定値を得ることが可能です。 引数byに与えるのも文字なので，ダブルクオーテーション（""）で囲んでください。 ここでは， 男女（ST04Q01）別のESCSの平均値と， 男女（ST04Q01）別の生まれ月（ST03Q01）の度数分布を示しています。
intsvy::pisa.mean(variable = "ESCS", by = "ST04Q01", data = jpn2012) ST04Q01 Freq Mean s.e. SD s.e
1 1 2974 -0.07228828 0.02079859 0.7075644 0.010526624
2 2 3211 -0.07079883 0.01925711 0.7001372 0.008863656intsvy::pisa.table(variable = "ST03Q01", by = "ST04Q01", data = jpn2012) ST04Q01 ST03Q01 Freq Percentage Std.err.
1 1 1 263 8.73 0.51
2 1 2 225 7.60 0.51
3 1 3 226 7.48 0.49
4 1 4 249 8.11 0.48
5 1 5 268 9.00 0.49
6 1 6 247 7.97 0.42
7 1 7 263 8.57 0.58
8 1 8 264 8.78 0.49
9 1 9 249 8.40 0.62
10 1 10 275 9.16 0.50
11 1 11 226 7.43 0.47
12 1 12 266 8.76 0.51
13 2 1 258 7.73 0.49
14 2 2 231 6.97 0.46
15 2 3 275 8.25 0.49
16 2 4 283 8.47 0.46
17 2 5 265 7.92 0.48
18 2 6 274 8.15 0.45
19 2 7 311 9.38 0.51
20 2 8 296 9.08 0.50
21 2 9 292 8.70 0.50
22 2 10 272 8.17 0.51
23 2 11 284 8.49 0.49
24 2 12 289 8.69 0.52pisa.rho関数を使うと，指定した変数間の相関係数（およびその標準誤差）を算出できます。 ここでは，PV1READ・PV1MATH・ESCSの3つの変数間の相関係数を出力しています。
intsvy::pisa.rho(variables = c("PV1READ", "PV1MATH", "ESCS"), data = jpn2012) PV1READ Rho PV1READ s.e. PV1MATH Rho PV1MATH s.e. ESCS Rho ESCS s.e.
PV1READ 1.000 0.000 0.861 0.006 0.289 0.024
PV1MATH 0.861 0.006 1.000 0.000 0.320 0.025
ESCS 0.289 0.024 0.320 0.025 1.000 0.000
attr(,"class")
[1] "intsvy.rho" "matrix" "array" 12.3 平均や分散の計算（PVを使う場合）
前の節で取り上げた関数のうち，pisa.mean関数にはPVを考慮した pisa.mean.pvという関数が存在します。 この関数では，読解リテラシーや数学リテラシーの平均値と標準偏差を算出することが可能です。 PVの設定にはやや癖があり，paste0関数を使うことになっています。 また，引数byを設定することで，男女別の平均値などを出力することもできます。 ここでは，数学リテラシー（paste0("PV", 1:5, "MATH")）の平均値の出力と， それをさらに男女（ST04Q01）別に計算する方法を示します。 paste0内の数値（1:5）は，PVの数に対応しています。 PISA2012ではPV1からPV5までが存在するので，paste0("PV", 1:5, "MATH")というわけです。 PISA2015以降ではPVが10個になりますので，paste0("PV", 1:10, "MATH")とします。
intsvy::pisa.mean.pv(pvlabel = paste0("PV", 1:5, "MATH"), data = jpn2012) Freq Mean s.e. SD s.e
1 6351 536.41 3.59 93.52 2.19intsvy::pisa.mean.pv(
pvlabel = paste0("PV", 1:5, "MATH"), by = "ST04Q01", data = jpn2012
) ST04Q01 Freq Mean s.e. SD s.e
1 1 3021 527.01 3.59 88.06 2.63
2 2 3330 544.88 4.62 97.41 2.6512.4 回帰分析
回帰分析を行う関数として，pisa.reg，pisa.reg.pv関数が用意されています。 前者がPVを使わない場合，後者はPVを使う場合の関数です。 ここでは，PV1READを従属変数，ESCSを独立変数とした場合の回帰分析の例を示します。 式で言うと，
intsvy::pisa.reg(y = "PV1READ", x = "ESCS", data = jpn2012) Estimate Std. Error t value
(Intercept) 543.78 3.32 163.92
ESCS 39.24 3.69 10.63
R-squared 0.08 0.01 5.95出力は，Estimateの列が推定値，Std.Errの列が標準誤差を示します。 t valueの列やはR-squaredの行は気にしなくて構いません2。 ここでは，(Intercept)がESCSがESCSが1単位上昇するとPV1READが39.24ポイント上昇するという関係があることになります3。 intsvyを使うと，回帰係数にも標準誤差を算出してくれます。 標準誤差の解釈は，平均値の標準誤差とほぼ同様で 「母集団でPV1READを従属変数・ESCSを独立変数とする回帰分析を行った場合， ESCSの係数は39.24±3.69のどこかにある」（と信じよう） という意味になります4。
PVを使う場合は，以下のようになります。 PVの書き方は，pisa.mean.pvの場合と同じです。
intsvy::pisa.reg.pv(
x = "ESCS", pvlabel = paste0("PV", 1:5, "READ"), data = jpn2012
) Estimate Std. Error t value
(Intercept) 543.59 3.34 162.58
ESCS 37.99 3.91 9.71
R-squared 0.08 0.01 5.42出力の読み方は，先ほどとほぼ同じです。 PV1からPV5までを考慮すると， ESCSが1単位上昇すると読解リテラシーが37.99ポイント上昇するという関係があることになります。 標準誤差もPV1のみの場合は3.69だったのに対し，PV1からPV5を使うと3.91になっています。 推算値の使い方でも説明しましたが， 一般的な傾向として，PVをすべて使うと，PVが1つの場合に比べて標準誤差が大きくなります。
pisa.reg関数・pisa.reg.pv関数のいずれについても，重回帰分析を行うことも可能です。 重回帰分析を行う場合は，引数xに独立変数（仮にA・B・Cとします）を c("A", "B", "C")の形で指定します。 ここでは，pisa.reg.pv関数で読解リテラシーを従属変数， 性別（ST04Q01）とESCS（ESCS）を独立変数とした場合の重回帰分析の例を示します。 式で言うと，
intsvy::pisa.reg.pv(
x = c("ST04Q01", "ESCS"), pvlabel = paste0("PV", 1:5, "READ"), data = jpn2012
) Estimate Std. Error t value
(Intercept) 577.85 5.10 113.30
ST04Q01 -22.53 3.54 -6.36
ESCS 38.00 3.94 9.65
R-squared 0.09 0.01 6.51出力の読み方は，pisa.reg関数と同じです。 この場合，「性別（ST04Q01）が同じ集団で比べると， ESCSが1単位上昇すると読解リテラシーが38.00上昇する」という関係があることになります5。
12.5 相関係数（PVを使う場合）
2025年1月時点では，intsvyでPVを使った場合の相関係数は計算できません。 理由は不明ですが，pisa.rho.pv関数が実装されていないからです。 PV1からPV5を使って相関係数を算出する場合は， pisa.rho関数を使って自分で推定値を計算する必要があります。 たとえば，読解リテラシー（PV1READからPV5READ）とESCSの相関係数が計算したい場合は， 次のように入力し，ESCS Rhoの列の出力を見る必要があります。
intsvy::pisa.rho(
variables = c("PV1READ", "PV2READ", "PV3READ", "PV4READ", "PV5READ", "ESCS"),
data = jpn2012
) PV1READ Rho PV1READ s.e. PV2READ Rho PV2READ s.e. PV3READ Rho
PV1READ 1.000 0.000 0.885 0.005 0.887
PV2READ 0.885 0.005 1.000 0.000 0.883
PV3READ 0.887 0.005 0.883 0.005 1.000
PV4READ 0.888 0.005 0.883 0.005 0.884
PV5READ 0.886 0.005 0.882 0.005 0.885
ESCS 0.289 0.024 0.273 0.024 0.277
PV3READ s.e. PV4READ Rho PV4READ s.e. PV5READ Rho PV5READ s.e. ESCS Rho
PV1READ 0.005 0.888 0.005 0.886 0.005 0.289
PV2READ 0.005 0.883 0.005 0.882 0.005 0.273
PV3READ 0.000 0.884 0.005 0.885 0.006 0.277
PV4READ 0.005 1.000 0.000 0.883 0.005 0.279
PV5READ 0.006 0.883 0.005 1.000 0.000 0.284
ESCS 0.026 0.279 0.025 0.284 0.025 1.000
ESCS s.e.
PV1READ 0.024
PV2READ 0.024
PV3READ 0.026
PV4READ 0.025
PV5READ 0.025
ESCS 0.000
attr(,"class")
[1] "intsvy.rho" "matrix" "array" ESCS Rhoの列を見ると，PV1READからPV5READとESCSの相関係数が出力されています。 この5つの値をRubinのルールに従って平均したものが読解リテラシーとESCSの相関係数になります6。 同様に，ESCS s.e.の列を使えば，相関係数の標準誤差を算出することも可能です。