上位ランキングの変動から人気の変遷を分析する方法
太郎丸 博
2015年8月11日
1 Beginning: きっかけ
2015年8月9日の第47回関西計量社会学研究会で
塚常健太「日本の新生児命名ランキングに関する基礎分析」
を聞いてたいへん面白かったので、分析法について思いつきでコメントしたのであるが、帰りの電車の中で色々と考えてみるに、私の提案した方法にはいささか難点があるように思えたので、この種のデータの分析法について、調べたことや思ったことを書き留めておこう。
2 Question: 上位ランキングだけのデータ
社会学では流行や人気の変化を調べたい場合がある。例えばどんな本や音楽の人気があるのか、それはどう変化したのか、といったことを調べた社会学者としては見田宗介が有名だろう。見田の場合は、そのような流行の変化を通して人々の社会意識の変化を明らかにしようとしたわけである。見田は統計的な分析はほとんど行っていないが、売上が上位1～10位までの本や楽曲の名前がベストテンといったかたちでわかることは多い。このようなデータを使ってうまく流行の変化を調べられないか、というのがここでの問題である。
この問題は、具体的には左側センサー・データをどう扱うか、という問題と、変化の激しさをどう比較するかという問題にわけて考えられる。
2.1 Left Censored Data: 左側センサー・データ
データは、毎年（あるいは毎月、毎週でもよい）のベストテンの書名や楽曲名である。これを使って本や楽曲の人気がどの程度変化したのか知りたいとする。このような変化の指標として順位相関係数を計算するという方法がある。相関係数が大きいほど順位は変化しておらず、小さいほど変化していると解釈できる。
しかし、相関係数を計算するためには、2つの時点の対象となる書名や楽曲の順位がわかっていないといけないのだが、11位以下の書名や楽曲の順位がわからないということはしばしばある。上であげた新生児命名ランキングに関しても同じ問題があり、正確な相関係数は計算できない。そこで、適当なデータ処理をして概算で相関係数を得る必要がある。
このような上位ランクだけがわかるデータは左側センサー・データの一種である。左側センサー・データとは、ある変数の値がある値以下であることはわかっているが、具体的にどんな値かはわからないようなデータのことである。以下では上記の塚常報告の資料から転載した1912～1921年に生まれた男女の名前、上位十位の順位を、左側センサー・データのサンプル・データとして分析していく。ただし、一年でもランクインした名前が別の年に11位以下になっている場合は一律に11という値を入力している（が、実際には何位なのかは不明である）。
setwd("C:/Users/Hiroshi/Dropbox/HP/documents")
NameData <- read.csv("NameData.csv") # データの読み込み
head(NameData, 18) # データの最初の18行だけ表示## name sex r1912 r1913 r1914 r1915 r1916 r1917 r1918 r1919 r1920 r1921
## 1 正一 male 1 7 6 10 11 8 11 11 11 11
## 2 清 male 2 5 2 1 2 2 1 2 1 1
## 3 正雄 male 3 3 3 4 8 6 10 7 10 10
## 4 正 male 4 4 5 5 11 9 8 8 9 9
## 5 茂 male 5 2 10 3 6 10 7 9 2 3
## 6 武雄 male 6 8 11 7 11 11 11 11 11 11
## 7 正治 male 7 11 11 11 11 11 11 11 11 11
## 8 三郎 male 8 6 4 2 3 1 2 1 3 2
## 9 正夫 male 9 10 11 11 11 11 11 11 11 11
## 10 一郎 male 10 11 11 8 5 4 4 5 6 6
## 11 正二 male 11 1 11 11 11 11 11 11 11 11
## 12 義雄 male 11 9 11 9 11 11 5 6 11 11
## 13 正三 male 11 11 1 11 11 11 11 11 11 11
## 14 勇 male 11 11 7 11 4 3 3 3 4 4
## 15 実 male 11 11 8 6 7 5 6 4 5 7
## 16 秀雄 male 11 11 9 11 9 11 11 11 11 11
## 17 辰雄 male 11 11 11 11 1 11 11 11 11 11
## 18 辰男 male 11 11 11 11 10 11 11 11 11 11summary(NameData[, 1:5]) # データセットの1～5列目の変数の要約統計量## name sex r1912 r1913
## キミ : 1 female:28 Min. : 1.000 Min. : 1.000
## きみ : 1 male :20 1st Qu.: 6.750 1st Qu.: 6.750
## キヨ : 1 Median :11.000 Median :11.000
## きよ : 1 Mean : 8.708 Mean : 8.708
## トミ : 1 3rd Qu.:11.000 3rd Qu.:11.000
## ハナ : 1 Max. :11.000 Max. :11.000
## (Other):42
## r1914
## Min. : 1.000
## 1st Qu.: 6.750
## Median :11.000
## Mean : 8.708
## 3rd Qu.:11.000
## Max. :11.000
## 例えば上のデータの下から 2番め（17行目）の「辰雄」という名前は 1916 年に急に 1位になったがその前後はランク外だったので、11 という値になっている。
とうぜん時間がたてば昔人気があった名前はどんどんトップテンから脱落し、新しい名前がトップテンを占めるようになる。このような変化があると、 11 という値がほとんどのセルをしめるようになってしまう。塚常報告では 1912～2013年の間の変化を調べようとしているので、このような左側センサーの影響は深刻である。これがここで論じる第一の問題である。
2.2 Comparison of Fluidity: 変化の激しさを比較する
塚常報告で検討されていた問題の1つは、男性の名前の人気と女性の名前の人気では、どちらのほうが変化が激しいか、というものであった。上記のように相関係数が変化の激しさを示すとすれば、男女で相関係数の大きさを比較すればこの問題に答えられる。
しかし、長期の時系列データが得られている場合、男性のほうが相関係数が大きいときもあれば、女性のほうが大きい時もある、という具合であまり一貫した結果が得られない場合も考えられる。例えば、塚常報告の 1912～1916年の名前順位の相関係数を、男女に関して計算した結果は以下の通りである。
n.years <- ncol(NameData) - 2 # 何年分のデータか
cor1 <- matrix(NA, n.years - 1, 2) # 相関係数を代入するための空の行列を用意
for(j in 1:2){ # 男女別に同じ計算をするので２回ループを回す
d0 <- subset(NameData,
subset = sex == levels(NameData$sex)[j], # 一方の性別だけ取って
select = 3: (n.years+2) # 1, 2列目は要らないので落とす
)
for(i in 1 : (n.years - 1)){ # 隣接する年次の相関係数を計算して cor1 に代入
cor1[i, j] <- cor(d0[,i], d0[, i + 1])
}
}
cor2 <- data.frame(sex=gl(2, n.years -1, labels=levels(NameData$sex) ), # cor1 の相関係数を
year=rep(1913:1921, 2), # ggplot で処理しやすいようにロング形式に変換
cor=c(cor1[,1], cor1[,2])
)
library(ggplot2) # 男女別に相関係数の折れ線グラフを描く
plot1 <- ggplot(cor2, aes(x=year, y = cor, group = sex, colour=sex))
plot1 + geom_line() + ggtitle("図1 前年の順位と当年の順位の相関係数") + ylab("Spearman の順位相関係数")
年による変化が激しく、男性の方が高くなったり、女性のほうが高くなったりで、どちらのほうが平均的に見て高いのか判断が難しい。それゆえ男女別にすべての時期を通算した平均的な相関係数を計算し、その差を検定したい。
また、相関係数が毎年上下動を繰り返す傾向が見られる。つまり、前年に相関係数が上がると翌年は下がる、という傾向である。これは前年に名前のランキンがあまり変化しないと、翌年は大きく変化しやすいということである。これは負の自己相関の一種といえる。このような傾向の存在は、単純な平均値の差の検定が満たすべき仮定に反するものである。それゆえ、分析に工夫が必要となる。
3 SEM アプローチ
私は最初、構造方程式モデリング (Structural Equation Modeling: SEM) がこの問題に対する最適なアプローチだと考えた。基本的には、1913年のランキングを 1912年のランキングに回帰させ、1914年のランキングは 1913年のランキングに回帰させ、というモデルを作り、これらの係数をすべて同じ値（下の図では b としている）になるよう制約をかける。下の図はこのようなモデルを図示したものである。
この b は、前年と当年のランクの平均的な関係の強さを示すと考えられる。 これを男女別に推定し、係数 b を男女で比較する。具体的には、男女で b が等しいと仮定したモデルと、異なると仮定したモデルを推定し、両者の適合度の違いを尤度比検定すればよい（あるいは BIC などをつかってもよい）。以下は実際の計算
library(lavaan)## This is lavaan 0.5-18
## lavaan is BETA software! Please report any bugs.model1 <- 'r1913 ~ c(b1, b2)*r1912
r1914 ~ c(b1, b2)*r1913
r1915 ~ c(b1, b2)*r1914
r1916 ~ c(b1, b2)*r1915
r1917 ~ c(b1, b2)*r1916
r1918 ~ c(b1, b2)*r1917
r1919 ~ c(b1, b2)*r1918
r1920 ~ c(b1, b2)*r1919
r1921 ~ c(b1, b2)*r1920' # 係数が男女で異なるモデル
model2 <- 'r1913 ~ c(b1, b1)*r1912
r1914 ~ c(b1, b1)*r1913
r1915 ~ c(b1, b1)*r1914
r1916 ~ c(b1, b1)*r1915
r1917 ~ c(b1, b1)*r1916
r1918 ~ c(b1, b1)*r1917
r1919 ~ c(b1, b1)*r1918
r1920 ~ c(b1, b1)*r1919
r1921 ~ c(b1, b1)*r1920' # 係数が男女で等しいモデル
sem1 <- sem(model1, data = NameData, group = "sex")
sem2 <- sem(model2, data = NameData, group = "sex")
coef(sem1)[c("b1", "b2")]## b1 b2
## 0.8791345 0.6526202anova(sem1, sem2)## Chi Square Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## sem1 88 2281.3 2352.4 164.39
## sem2 89 2290.8 2360.0 175.88 11.491 1 0.0006995 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fitMeasures(sem1)[c("cfi", "rmsea", "srmr")]## cfi rmsea srmr
## 0.7905302 0.1901817 0.5194765b1 が男性、b2 が女性なので、男性のほうが係数が大きい（つまりランキングの変化は小さい）。尤度比検定の結果も有意である。
ちなみに CFI や RMSEA を見ると、上のモデルはどちらも適合度が著しく低い。つまり名前のランクは2年以上前のランクからも影響を受けるか（以前ランクが高かった名前は後にはランクが下がりやすい？）、あるいは誤差相関がある（急激にランクの上がった名前は翌年には逆にランクが下がりやすい？）と考えられるが、この辺りは既存の知見を参考にしながら最適なモデルを探索すればよい。
この分析は非標準化解で行ったが、あらかじめ用いる変数を標準化すれば回帰係数は相関係数と一致するのでどうしても標準化したい人はすればよい（変数のスケールはすべて同じなので標準化は必要ないと思うが）。また、ランク外は機械的に 11位にしたまま分析したが、これも好きなように処理すればよい。例えば M-plus ならば左側センサーを仮定した SEM などできそうな気がする（が未確認）。また塚常報告のように、一定期間以上ランク外の場合、欠損値にするといった処理をして、リストワイズではなくペアワイズや ML、Imputation などの方法で欠損値処理してもよかろう。
3.1 Limitation of SEM
SEM アプローチは、上記のように短いスパンで変化を見る場合、非常に有効である。しかし、塚常報告のように 100年ものあいだの変化を見る場合にはいささか無理が生じる。相関係数（あるいは分散と共分散）の推定が怪しくなりやすいのである。
名前のランキングの場合、毎年 1～2 程度の新しい名前がランキングにはいってきた。それゆえ、毎年1つ新しい名前がランキングに入るとすると、100年間で約100の名前が一度はランクインすることになる。つまり、データとしては 100ケースあるが、毎年順位がわかるのはそのうち 10 だけで、残りの 90 の名前は 11位以下だということしかわからない。そのため、1912年と 1913年のランクの相関を見るとき、9割近くの名前はどちらの時点でも 11位以下ということになる。これが相関係数を過大に推定させはしないか、というのが第一の懸念である（この辺りは左側センサーをどう処理するかに強く依存する）。
また、1912年と 2013年のランキングの相関をとると、必ず負の相関になる。なぜなら、100年間にわたって人気を維持した名前は殆ど無いので、1912年のトップテンは 2013年には 11位以下、2013年のトップテンは 1912年には 11位以下になるからである。もしも塚常報告のように一定期間ランク外の名前には欠損値を与えるとすると、ペアワイズで処理しても 1912年と 2013年のランキングの相関は有効ケースがなくなり、計算出来ないだろう。
1912年と 2013年のランキングの相関は、関心の対象とはなりにくいと思う（塚常報告でも問題になっていない）が、構造方程式モデリングでは、用いるすべての変数の分散共分散行列からモデルのパラメータを推定するので、分散共分散行列の値が歪んでいれば、とうぜん得られる推定値も歪んでしまう。
このような歪みは致命的なものなのか、それともあまり実害ない程度のものなのか、シミュレーションなどして確かめてみないとよくわからないので、今のところ、長期的な変化の趨勢を見る際にお勧めとはいえないのである。個人的には SEM のエレガントさは好きなのだが、著しく多数の左側センサー・ケースが生じることがネックと言えよう。
4 Time Series Approach: 相関係数の時系列データ分析
SEM に代わるもう一つの方法として考えられるのは、図1 のような相関係数のデータを時系列データとみなして分析するということである。実際に観察されているのは相関係数そのものではなく、名前のランキングなのであるが、相関係数の変動が時系列モデルに従うと仮定できるのならばこの方法も有効と思われる。この方法のメリットは、隣接する時点（例えば 1912年と 1913年）の相関係数だけを計算すればよく、離れた時点の相関は無視できる点にある。それゆえ、塚常報告のように長期間ランク外の名前は欠損値として処理することもできる。
4.1 Trend: トレンドの有無のチェック
時系列分析では、右肩上がりの上昇傾向などがある場合、そのようなトレンドはあらかじめコントロールするか、除外するといった方法が取られる。ここでも上の 1912-1921年の名前を例にトレンドのチェックから始めよう。まず男女別に OLS で各年の相関係数を調査年に回帰させた結果が以下である。
lm0 <- vector("list", 2)
for(i in 1:2){
lm0[[i]] <- lm(cor ~ year, data = cor2, subset= sex==levels(sex)[i])
print(levels(cor2$sex)[i])
print(summary(lm0[[i]]))
}## [1] "female"
##
## Call:
## lm(formula = cor ~ year, data = cor2, subset = sex == levels(sex)[i])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.186406 -0.048181 0.007482 0.048221 0.177239
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -59.71206 29.24612 -2.042 0.0805 .
## year 0.03147 0.01526 2.063 0.0780 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1182 on 7 degrees of freedom
## Multiple R-squared: 0.3781, Adjusted R-squared: 0.2892
## F-statistic: 4.256 on 1 and 7 DF, p-value: 0.07803
##
## [1] "male"
##
## Call:
## lm(formula = cor ~ year, data = cor2, subset = sex == levels(sex)[i])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16128 -0.11679 -0.01911 0.09084 0.20136
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -128.62866 33.43465 -3.847 0.00632 **
## year 0.06745 0.01744 3.867 0.00615 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1351 on 7 degrees of freedom
## Multiple R-squared: 0.6812, Adjusted R-squared: 0.6356
## F-statistic: 14.96 on 1 and 7 DF, p-value: 0.006153男女とも年次が最近になるほど相関が大きくなる傾向が見られる。女性は .03, 男性は .07 程度の傾きで女性は 10% 水準、男性は 1% 水準で有意である。女性に関しては右肩上がりのトレンドがあるかどうかは不確かであるが、男女を合わせて1つのサンプルと考えるならば、性別と年次の交互作用を仮定するほど男女の年次の傾きに違いがあるわけではない。これらを考え合わせると、男女とも右肩上がりのトレンドがあると判断すべきだろう。
次に誤差相関の有無を調べる。上の回帰分析の結果得られた残差（誤差）に系列相関（自己相関）がないか相関係数と偏相関係数を男女別に計算した。その結果が下の図である。これらはコレログラムとも呼ばれるものである。
par(mfrow=c(2,2))
for(i in 1:2){
acf(resid(lm0[[2]]), main=paste(levels(cor2$sex)[i], "の残差の系列相関"),
ylab="自己相関係数")
pacf(resid(lm0[[2]]), main=paste(levels(cor2$sex)[i], "の残差の偏系列相関"),
ylab="自己偏相関係数")
}
この図は横軸の Lag が何時点前との相関かを示し、縦軸が相関係数、または偏相関係数の値を示している。相関係数に関しては Lag の値は 0 から始まるが、これは完全に同じ残差どうしの相関をとっているだけなので、必ず 1 になる。重要なのは、Lag が 1 以上の時の自己相関係数や偏自己相関係数である。青い点線は自己相関係数や自己偏相関係数の帰無仮説のもとでの両側 95% 水準の限界値であるので、この青線よりも 0 から離れていたら、自己（偏）相関係数は有意である。結果を見ると、一時点前の誤差と -0.3 程度の自己相関があるが、有意ではない。またそれ以上前の時点の誤差ともある程度は負の自己相関があるようであるが、いずれも有意ではないので、無視することにする。つまり先ほどの OLS の結果を採択すればよい。
4.2 Cross-Sexual Difference: 男女差の検討
まず、男女の残差の間に相関がないかチェックしてみた。相関がある場合それを考慮した分析が必要になるからである。結果を見ると 0.4 程度の相関があるがやはり有意ではないので、両者は独立であると仮定して分析を進める。つまり、最初の OLS の推定結果をやはりそのまま採択するということである。この分析結果から予測される回帰直線を示したのが、下の図である。
cor.test(resid(lm0[[1]]), resid(lm0[[2]]))##
## Pearson's product-moment correlation
##
## data: resid(lm0[[1]]) and resid(lm0[[2]])
## t = 1.2163, df = 7, p-value = 0.2633
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3410328 0.8468878
## sample estimates:
## cor
## 0.4176905plot1 + geom_point(aes(colour = sex)) + geom_smooth(aes(group = sex), method = "lm") +
ggtitle("図2 前年の順位と当年の順位の相関係数とその回帰直線") +
ylab("Spearman の順位相関係数")
赤い線が女性、青い線が男性の回帰直線である。回帰直線の周りの影は直線の 95% 信頼区間である。この図を見ればわかるように男性と女性の信頼区間は完全に重なっており、有意な違いがあるとは言えない。以下のように男女を合わせたサンプルで相関係数を性別と年次に回帰させてもみたが、やはり有意な男女差はない。
lm1 <- lm(cor ~ sex + year, data = cor2)
summary(lm1)##
## Call:
## lm(formula = cor ~ sex + year, data = cor2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.21525 -0.06089 -0.02230 0.10809 0.24920
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -94.19775 23.23147 -4.055 0.001037 **
## sexmale 0.05478 0.06258 0.875 0.395191
## year 0.04946 0.01212 4.081 0.000983 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1328 on 15 degrees of freedom
## Multiple R-squared: 0.5374, Adjusted R-squared: 0.4757
## F-statistic: 8.712 on 2 and 15 DF, p-value: 0.003084つまり、有意な男女差は無いという結果であった。
4.3 Comparison: SEM と時系列分析の比較
SEM と時系列分析ではかなり分析結果が違っていた。これは誤差がどのように発生すると仮定して、そのバラつきをどう推定するかがまったく違っているからである。データの性質から言うと SEM のほうが適切な推定であるが、すでに触れたように、ランク外の名前がたくさんある場合、必ずしも SEM のほうが適切とは言い切れない。現時点ではどちらがいいかは判断できないので、シミュレーションなどでバイアスや検定力などを検討してみる必要があろう。
5 Correlation between Left Censored Variables: 左側センサーされた変数どうしの順位相関係数
おまけ。R の場合、NADA パッケージに cenken() という関数があり、Kedall の tau を左側センサーを考慮して計算してくれる。使い方は
cenken(変数 y, y のセンサーを示す変数, 変数 x, x のセンサーを示す変数)
で、変数 y と x の順位相関係数 (Kendall’s tau) を計算してくれる。「y, x のセンサーを示す変数」とはセンサーされているとき（新生児の命名ランキングの例では 11位以下のとき） TRUE、されていないとき（新生児の命名ランキングの例では 10位以上のとき） FALSE となる変数のこと。以下は使用例。
library(NADA)## Loading required package: survival
##
## Attaching package: 'NADA'
##
## 以下のオブジェクトは 'package:stats' からマスクされています:
##
## corcenken(d0$r1912, d0$r1912==11, d0$r1913, d0$r1913==11)## $slope
## [1] 0.5251465
##
## $intercept
## [1] 1.899414
##
## $tau
## [1] 0.03157895
##
## $p
## [1] 0.832328
##
## attr(,"class")
## [1] "cenken"Kendall’s の tau だけでなく、回帰分析の切片と傾き、および傾きの検定の p 値も返す。計算の仕方などは cenken() のヘルプに参考文献が載っているので、それを見ると良いだろう。