R で e-typing のスコアを推定する
以前、こちらの記事で「e-typing の腕だめしテストでの成績を自動で取得し、毎回 csv に保存する」ということを紹介しました。
データが溜まってきたので、とりあえず今回は手軽なお遊びとして
スコアの計算方法を探る
ということをやってみます。
データの内容は以下の通りです。また、かっこがついているものは本記事ではこの表記を使うという意味です。
- スコア：タイピング能力を定量化したものであり、本記事での予測対象
- レベル：スコアによって一意に決まる称号
- 入力時間（タイム）：開始から全ての文章を打ち終えるまでにかかった時間
- 入力文字数（正入力数）：正しく入力された文字数
- ミス入力数（誤入力数）：ミス入力数
- WPM：Word Per Minute の略、1分間に何文字正しく打てるか
- 正解率：正入力数から誤入力数を引いて、正入力数で割ったもの
「レベル」はスコアによって一意に決まるため、今回は使用しません。また、一応特徴量の数は5個ですが、WPM と正解率は他の特徴量から計算されるものなので、最も細かい粒度の特徴量は入力時間、正入力数、ミス入力数の3つです。
WPM、正解率の計算方法は以下の通りです。
$$
\text{WPM} = \text{正入力数} \times \frac{60}{\text{タイム}} \\
\text{正解率} = \frac{\text{正入力数} – \text{誤入力数}}{\text{正入力数}}
$$
（正直、\(\frac{\text{正入力数}}{\text{正入力数} + \text{誤入力数}}\)にすべきだと思うんですが、サイトの方針のようです）
なお、答え（本当の計算方法）は最後に掲載してます、これはデベロッパーツールで e-typing のサイトのソースコードを辿って見つけました。
ここで、当然のことではありますが、線形回帰の場合アルゴリズム上ほぼ確実に切片項が式に現れますし、回帰係数も連続値なので「係数情報まで当てる」は無理です。なので、今回でのクリア条件は「スコア計算に使われている変数、変数同士の掛け算 or 割り算のパターンを特定する」です。
ということで、本当の計算方法を知らない体で色々やっていきます。
スコアの計算方法を探る
以下の方針で進めて行きます
- 分布確認
- 特徴量同士をかけたり割ったりしていろんな特徴量を作成
- 情報量基準やスパース推定法を使ったモデル選択
- テストデータに対する当てはまりをみて代表モデルを決定
3でクロスバリデーションを使ってハイパーパラメータを決定しますが、このときに時系列上の制約は考慮しません。例えば、「2021年12月中のデータで予測モデルを学習し、2021年11月中のデータを予測する」ということを許容します（本件ではなんら問題はない）。
また、これは今回のような問題に限った話なので要注意ですが、手元のデータに誤差はないので過学習ウェルカムなモデリングでいきます。そもそも解析的に求まっているスコアに対して誤差を仮定した確率モデルを適用するのは問題設定としてどうかみたいなのはありますが、これはお遊びである。
計算式推定の方針ですが、タイピングスコアがそれほど複雑な計算式だとは思えないので、線形重回帰分析をベースに計算式を探って行きます。
まず、サンプルサイズ \(N\) のとき、ベースとなるモデルは以下の通りです。
$$
\text{スコア}_i = \beta_1 \text{タイム}_i + \beta_2 \text{正入力数}_i + \beta_3 \text{誤入力数}_i + \\
\beta_4 \text{WPM}_i + \beta_5 \text{正解率}_i + \epsilon_i, \\
\epsilon_i \sim \mathcal{N}(0, \sigma^2), i = 1, \ldots, N.
$$
この線形重回帰をベースに、情報量基準を使ってどんな特徴量を使うのが正解に近いのかを探って行きます。
1. 分布確認
2行目でエンコーディングを 'CP932’ にしていますが、UTF-8 で保存しているのであれば locale = の部分は削除（デフォルト）で大丈夫です。
library(dplyr); library(tidyverse)
data_original = read_csv('typing_score.csv', locale = locale(encoding = 'cp932'))
print(data_original, n = 5, width = Inf)# A tibble: 133 x 8
date スコア レベル 入力時間 入力文字数 ミス入力数 WPM 正確率
<chr> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 2021/11/30 416 Professor 51.5 378 7 440. 98.1
2 2021/11/30 445 Professor 50 380 3 456. 99.2
3 2021/11/30 416 Professor 54.9 384 1 420. 99.7
4 2021/11/30 448 Professor 52.5 393 0 449. 100
5 2021/11/30 434 Professor 51.2 383 4 448. 99.0
# … with 128 more rows使うデータだけ取り出し、変数名を本記事と合わせてから分割します。また、最後に推定されたモデルが学習に使っていないデータ（テストデータ）に対してどれほどの精度で予測できるのかを確認するために、学習データには120サンプルだけ使用します。あと、正解率を小数点での表記にしておきます。
data = data_original[1:120,c(2,4:8)]
colnames(data)[2:4] = c('タイム', '正入力数', '誤入力数')
data = data %>% mutate(`正確率` = `正確率` / 100)
どんなデータが相手だろうが分布確認はします。
まず、分布の確認とピアソンの積率相関係数、散布図をみてみます。
library(GGally)
ggpairs(data = data,
upper = list(continuous = wrap("cor", size = 3, alpha = 0.7)),
diag = list(continuous = wrap('densityDiag', alpha = 0.7)),
lower = list(continuous = wrap('smooth', size = 0.5, alpha = 0.7))) +
theme_gray(base_family = "HiraKakuPro-W3")
一応特徴量間の関係について説明すると、タイムと正入力数が正相関なのは、このタイピングテストでは「正入力数はテスト終了までに打つ必要がある文字数」であるため、入力者のタイピング速度が大きくぶれなければ必然的に相関します。WPM は定義から、正入力数が増えたら増加、タイムが増えたら減少しますが、正入力数とタイムが正相関するため、タイムとの単相間は見かけ上では正になっています。
他の変数の影響も考慮するために、偏相関係数を計算します。偏相関係数は、例えば \(X \in \mathbb{R}^p\) において、\(X_1\) と \(X_2\) の相関係数を、\(X_3, \ldots, X_p\) の影響を取り除いた上で計算する量です。ここでいう「取り除く」というのは、「\(X_3, \ldots, X_p\) が与えられたという条件付きの相関係数」と定義されます。詳しくは[1]、[2]をご参照ください。
単相関、偏相関の順にヒートマップで可視化すると以下のようになります。
library(ppcor); library(reshape2)
library(gridExtra)
cor_gg = cor(data) %>% melt()
pcor_gg = pcor(data)$estimate %>% melt()
cor = ggplot(cor_gg, aes(x = Var1, y = Var2, fill = value)) +
geom_tile(show.legend = FALSE) +
geom_text(aes(label = round(value, 3))) +
theme_gray(base_family = "HiraKakuPro-W3") +
labs(x = '', y = '')
pcor = ggplot(pcor_gg, aes(x = Var1, y = Var2, fill = value)) +
geom_tile(show.legend = FALSE) +
geom_text(aes(label = round(value, 3))) +
theme_gray(base_family = "HiraKakuPro-W3") +
labs(x = '', y = '')
layout = matrix(1:2, nrow = 1)
fig_y_all = grid.arrange(cor, pcor, layout_matrix = layout) 
タイムと WPM の偏相関係数は負になっていることがわかります。ただし、この可視化の場合スコアを含めてしまっているため、WPM 正入力数、タイムとの偏相関係数の絶対値は小さくなっています。ちなみに、スコアを除いた偏相関係数の図は以下の通りです。
WPM とタイム、正入力数間の関係がより定義と整合する値になっていると思います。
以上の結果から、スコアの推定のポイントとなりそうなところをまとめると、
- スコアは WPM、正解率と偏相関が高い
- スコアはタイム、正入力数との偏相関の絶対値は小さいが、これは WPM や正解率に代替されている可能性がある
- スコアは誤入力数との単相関は小さいが、偏相関係数の絶対値でみると WPM、正解率の次に大きい
という感じです。簡単な可視化しかしていませんが、WPM と 正解率、誤入力数の影響が大きそうだということを念頭におきながら、モデリングしていきます。
2. 特徴量エンジニアリング
このセクションは色々と難しい部分があります。例えば、5この特徴量を使って二次の項を計算して加えたら、特徴量の総数は \(5 + 5 + {}_5 \mathrm{C}_2\ = 20\) 個になりますし、三角関数や対数、指数関数など、変換の種類とその組み合わせを考えると特徴量は膨大な数になります。特徴量を増やしすぎると、そもそも特徴量の解釈が難しくなることや、サンプルサイズとの比から推定結果が不安定になるといった問題が生じます。
今回は「タイピングのスコア」ということなので、sin や 対数といった変換をわざわざ行うのは考えにくく、また逆数は0割が起こることから計算式には入っていないと考え、二次の項までを作成してモデル選択を行なって行きます。
ちなみに、本来こういった特徴量作成はドメイン知識を活かした仮説ベースで行うことが多いです
ということで、各変数の二次の項を作ってみます。
data_quad = data
for (i in 1:(ncol(data) - 1)) {
for (j in i:(ncol(data) - 1)) {
data_quad = cbind(data_quad, data[,i + 1] * data[,j + 1])
colnames(data_quad)[ncol(data_quad)] = str_c(colnames(data)[i + 1], ':', colnames(data)[j + 1])
}
}力づくで二次形式を作っております。
このデータをベースにモデル選択をやっていきます。
3. 情報量基準を使ったモデル選択
情報量基準として特によく使われるものに AIC と BIC があると思います。今回は「正解の計算式を当てる」という目的なので、モデル選択の一致性を持つ BIC を使います。定義式はサンプルサイズ \(N\)、モデルに使った特徴量の数 \(p\)、正規分布に基づく対数尤度 \(\hat{L}\) を使って、以下で表されます。
$$
\text{BIC} = -2 \log \hat{L} + p \log N
$$
\(\hat{L}\) は正規分布を仮定した尤度関数に対数をとったものであり、以下のようにかけます。
$$
\hat{L} = -\frac{N}{2}\log {2\pi} – \frac{N}{2}\log {\sigma^2} – \frac{1}{2\sigma^2}\sum^N_i{(y_i – \hat{y}_i)^2}
$$
よって、BIC は
$$
BIC = N\log {2\pi\sigma^2} + \frac{1}{\sigma^2}\sum^N_{i = 1}{(y_i – \hat{y}_i)^2} + p\log {N}
$$
とかけます。\(\sigma^2\) は誤差分散ですが、予測残差の不偏分散で代用します。
$$
\sigma^2 = \frac{1}{N – p}\sum^N_i{(y_i – \hat{y}_i)^2}
$$
より詳しくは scikit-learn 公式サイト とその中で紹介されている元論文をご参照ください。
通常の線形重回帰モデルでの変数選択においては、leaps パッケージの regsubsets 関数を使います。
こちらは、総当たり法で変数の組み合わせを探していますが、nvmax を指定することでうまい具合にサボっています（サボらないと \(2^p\) 通り試すことになる）。
library(leaps)
regsubsets = regsubsets(`スコア` ~ ., data = data_quad, nvmax = 10)
reg_summary = summary(regsubsets)
index_best = which(reg_summary$bic == min(reg_summary$bic))
(coef_bic = coef(regsubsets, index_best) %>% round(digits = 2)) (Intercept) WPM `誤入力数:誤入力数` `WPM:正確率`
-0.71 -1.99 0.01 2.99 WPM、誤入力の2乗、WPM×正解率の3つが BIC 基準で選ばれました。
続いて、Lasso によるスパースモデリングを行い、罰則パラメータを BIC 基準で決めることによりモデル選択を行います。
まず、推定されたモデルを使って BIC を計算する関数を定義します。
CalcBIC = function(y, pred, df) {
num_sample = length(y)
var = sum((y - pred)^2) / (num_sample - df)
BIC = num_sample * log(2 * pi * var) + sum((y - pred)^2) / var + df * log(n_sample)
return(BIC)
}罰則パラメータの候補を適当に作って、それぞれに置いて BIC を計算し、最も BIC が小さくなる罰則パラメータでモデルを学習します。
library(glmnet)
lambda = seq(1e-4, 0.1, length.out = 1000)
bic_lambda = data.frame()
for (lambda in lambda) {
model = glmnet(x = as.matrix(data_quad[,-1]), y = data_quad$スコア, alpha = 1, lambda = lambda)
pred = predict(model, as.matrix(data_quad[,-1]))
bic = CalcBIC(y = data_quad$スコア, pred = pred, df = model$df)
bic_lambda = rbind(bic_lambda, data.frame(bic = bic, lambda = lambda))
}
ggplot(bic_lambda, aes(x = lambda, y = bic)) + geom_point()
bestlambda = bic_lambda$lambda[which(bic_lambda$bic == min(bic_lambda$bic))]
model = glmnet(x = as.matrix(data_quad[,-1]), y = data_quad$スコア, alpha = 1, lambda = bestlambda)
model$beta
20 x 1 sparse Matrix of class "dgCMatrix"
s0
タイム .
正入力数 .
誤入力数 .
WPM .
正確率 2.617408e+02
タイム:タイム .
タイム:正入力数 .
タイム:誤入力数 .
タイム:WPM .
タイム:正確率 .
正入力数:正入力数 .
正入力数:誤入力数 .
正入力数:WPM .
正入力数:正確率 2.185526e-03
誤入力数:誤入力数 .
誤入力数:WPM -4.173916e-04
誤入力数:正確率 .
WPM:WPM .
WPM:正確率 9.826801e-01
正確率:正確率 2.649470e+02正解率、正入力数×正解率、誤入力数× WPM、WPM ×正解率、正解率の2乗が選択されました。
次に、BIC 基準ではなく 10-folds Cross-validation で罰則パラメータを決定してみます。
set.seed(0)
bestlambda_cv = cv.glmnet(x = as.matrix(data_quad[,-1]), y = data_quad$スコア, alpha = 1, nfolds = 10)$lambda.1se
model_cv = glmnet(x = as.matrix(data_quad[,-1]), y = data_quad$スコア, alpha = 1, lambda = bestlambda_cv)
model_cv$beta20 x 1 sparse Matrix of class "dgCMatrix"
s0
タイム .
正入力数 .
誤入力数 .
WPM .
正確率 231.5852277
タイム:タイム .
タイム:正入力数 .
タイム:誤入力数 .
タイム:WPM .
タイム:正確率 .
正入力数:正入力数 .
正入力数:誤入力数 .
正入力数:WPM .
正入力数:正確率 .
誤入力数:誤入力数 .
誤入力数:WPM .
誤入力数:正確率 .
WPM:WPM .
WPM:正確率 0.9642199
正確率:正確率 290.6902518だいぶ絞られました！とにかく正解率にこだわった式ですね、もはやこんな採点方式もアリなのではと思います。
4. テストデータを使って代表モデル（回答）の決定
現状、データから最もそれっぽい計算式だと推された結果がこちら。
総当たり法（regsubsets）
$$
\text{スコア}_i = -0.71 – 1.99 \times \text{WPM}_i + 0.01 \times
\text{誤入力数}^2_i + 2.99 \times \text{WPM}_i \times \text{正解率}_i
$$
BIC 基準 Lasso
$$
\text{スコア}_i = -520.59 + 261.74 \times \text{正解率}_i + 0.002 \times
\text{正入力数}_i \times \text{正解率}_i – 0.0004 \times \text{誤入力数} \times \text{WPM}_i \\
+ 0.98 \times \text{WPM}_i \times \text{正確率}_i + 264.95 \times \text{正確率}_i^2
$$
CV Lasso
$$
\text{スコア}_i = -507.96 + 231.59 \times \text{正解率}_i + 0.96 \times \text{WPM}_i \times \text{正確率}_i + 290.69 \times \text{正確率}_i^2
$$
3つの候補がありますが、最初に取り置きしておいた24サンプルのテストデータを使って、最も精度が良いものを最終結果とします。ちなみに制度は MAPE、自由度調整済み決定係数で決めます。
コードは長いので割愛しますが、学習データと同様、テストデータも二次形式を作成し、推定したモデルの入力とします。指標の中身がわかるように、評価関数と結果だけ載せておきます。
metrics = function(y, pred, df) {
num_sample = length(y)
MAPE = sum(abs((pred - y) / y)) / num_sample * 100
R2_DFadj1 = sum((y - pred)^2) / (num_sample - df - 1)
R2_DFadj2 = sum((y - mean(y))^2) / (num_sample - 1)
R2_DFadj = 1 - R2_DFadj1 / R2_DFadj2
return(tibble(MAPE = MAPE, R2_DFadj = R2_DFadj))
}pred_bic = coef_bic[1] + as.matrix(test_quad[,c('WPM', '誤入力数:誤入力数', 'WPM:正確率')]) %*% coef_bic[-1]
pred_lasso_bic = predict(model, as.matrix(test_quad[,-1]))
pred_lasso_cv = predict(model_cv, as.matrix(test_quad[,-1]))
metrics(y = test$スコア, pred = pred_bic, df = 3)
metrics(y = test$スコア, pred = pred_lasso_bic, df = 5)
metrics(y = test$スコア, pred = pred_lasso_cv, df = 3)結果、MAPE が0.095%、自由度調整済み決定係数が1.00という脅威の記録を出した総当たり法によるモデルを採用します！さすが解析的にスコアが計算されているだけあってテストデータへの当てはまりも凄まじいものですね。（一応言っておきますが、もし実際の分析の場でこんなことが起きようものならほぼ確実にデータリーケージです）
答え合わせ
最終的な回答がこちら：
$$
\text{スコア}_i = -0.71 – 1.99 \times \text{WPM}_i + 0.01 \times
\text{誤入力数}^2_i + 2.99 \times \text{WPM}_i \times \text{正解率}_i
$$
冒頭で述べた通り、使われている変数や演算に着目した回答にすると、
「スコア計算には”WPM”、”誤入力数”、”正確率”の3つが使われており、”誤入力数”の2乗、”WPM”と”正確率”の積、WPM の3つの量に定数をかけて足し合わせたものがスコアとなっている」
です。
これに対し、サイトのソースコードを漁って見つけた答えがこちら。
$$
\text{スコア} = 60 \times \frac{\text{正入力数} – \text{誤入力数}}{\text{タイム}} \times
(\text{正解率}) ^ 2
$$
まあ遠すぎず近すぎずって感じでしょうか。ちなみに式変形してどうにか近づけようとしましたが諦めました（これ以上こんなことに時間を割くべきではない）。
皆さんはこんなことしてないでしっかりタイピング練習したり重回帰分析を学んでください。
参考
[1] Pythonで多変量解析、3変数以上の偏相関係数を算出してみた
