餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

人文系（哲学）大学教員のためのRによる成績評価方法（α版）

R言語

最近は理工系はもちろん、社会科学系の大学教員も研究で統計解析をしていることが多いので、Rで成績評価する方法を解説しても需要はあまり無さそうな気もするのですが、人文系でニーズがあるようなので。
名簿と試験の成績とレポートの評価の3つのCSVファイルがあるとして、それをつないで総合得点と成績を出し、基本統計量を出すところまでやってみましょう。なお、データは完全に架空のものです。

ファイルの読み込み

今回はウェブサイトに置いてあるファイルを読み込みますが、実際の場合はURLではなくファイル名を指定してください。なお、ファイルが読み込めないよと言うときは、`getwd()`で作業フォルダを確認し、RGuiのメニューからディレクトリの変更を行うか、`setwd("フォルダ名")`で移動することができるのを思い出しましょう。

# header = TRUEを指定しているので、1行目からデータフレームの列名が設定される
# ファイルがタブ区切りにするときは、sep = "\t"
# Windows環境で古い文字コード（ShiftJIS/CP932/Windows-31J）を使っている場合は、fileEncodingはcp932
df01 <- read.table("http://wh.anlyznews.com/R/dataset/grading/list.csv", 
    header = TRUE, sep = ",", fileEncoding = "utf-8")
df02 <- read.table("http://wh.anlyznews.com/R/dataset/grading/exam.csv", 
    header = TRUE, sep = ",", fileEncoding = "utf-8")
df03 <- read.table("http://wh.anlyznews.com/R/dataset/grading/report.csv", 
    header = TRUE, sep = ",", fileEncoding = "utf-8")

ファイルをデータフレームに読み込んだら、以下のように`summary`や`head`などで中身を確認したほうがよいです。

summary(df01)
summaryの例

テキスト型のところは長さ、数値型のところは基本統計量が出ます。

以下では`head`は頭から3行みていますが、3行でなくてもよいです。むしろ10行ぐらいは見たほうがよいと思います。

head(df01, 3)

なお、list.csvは、

学部 学籍番号 出席回数
文学部 L0001 10
文学部 L0002 13
現代社会学部 S0001 12

exam.csvは、

学部 学籍番号 期末試験
文学部 L0001 85
文学部 L0002 28
現代社会学部 S0001 27

report.csvは、

学部 学籍番号 レポート評価
文学部 L0001 D
現代社会学部 S0001 D
現代社会学部 S0002 A

と言うような構造になっています。
出席回数や点数が数字ではなく文字列になっている場合は、CSVファイルの該当行にNA以外の文字が混じっていることが多いです。文字化けは文字コードの不整合が原因になります。最近のアプリケーションはutf-8で出力してくれるので、もう見なくなりつつありますが。

データフレームの結合

データフレームが3つばらばらだと扱いづらいので、`merge`を使って結合します。

df04 <- merge(
     merge(
         df01,
         df02,
         by = c("学部", "学籍番号"),
         all.x = TRUE),
     df03,
     by = c("学部", "学籍番号"),
     all.x = TRUE)

内側の`merge`でデータフレームdf01とdf02をつないでいます。`by = c("学部", "学籍番号")`は学部と学籍番号をつなぎあわせるキーにするという指定です。`all.x`は結合される側に欠損値があっても、その行を残すという指定です。外側の`merge`で、それにdf03をつないでいます。つないだ結果はdf04に入ります。

この時点ではdf04の頭3行は以下のようになっています。

学部 学籍番号 出席回数 期末試験 レポート評価
現代社会学部 S0001 12 27 D
現代社会学部 S0002 11 NA A
現代社会学部 S0003 10 NA D

成績をつける

試験の点数とレポート評価から総合点をつけ、それを評価に換算します。

総合点をつける

点数のつけ方は色々とあると思いますが、レポート評価ごとに点数を設定し、それを単純に期末試験の点数に合算、上限は100とします。

# 未提出者はE評価/レポート未提出による未受は無し（以下の1行をコメントアウトするとNAで未受扱い）
df04$レポート評価[is.na(df04$レポート評価)] <- "E" 
# レポートの評価を点数に換算
df04$レポート点数[!is.na(df04$レポート評価)] <- unlist(sapply(df04$レポート評価, 
    switch, "A"=20, "B"=15, "C"=10, "D"=5, "E"=0))
# 期末試験の点数とレポート点数を上限100で合算。
df04$合計点数 <- pmin(with(df04, 期末試験 + レポート点数), 100)

df04に列が追加されたので、頭10行を見てみましょう。

学部 学籍番号 出席回数 期末試験 レポート評価 レポート点数 合計点数
現代社会学部 S0001 12 27 D 5 32
現代社会学部 S0002 11 NA A 20 NA
現代社会学部 S0003 10 NA D 5 NA
現代社会学部 S0004 11 90 D 5 95
現代社会学部 S0005 13 97 B 15 100
現代社会学部 S0006 12 37 C 10 47
現代社会学部 S0007 11 85 D 5 90
現代社会学部 S0008 9 34 B 15 49
現代社会学部 S0009 9 87 E 0 87
現代社会学部 S0010 10 NA C 10 NA

評価をつける

評価のつけ方も色々とあると思うのですが、例として思いつく方法を3つ計算します。

## 絶対評価 ##
grade <- list("D"=0, "C"=50, "B"=65, "A"=80)
label_absolute <- "絶対評価"
df04[, label_absolute] <- rep(NA, nrow(df04))
for(g in names(grade)){
    df04[, label_absolute][df04$合計点数 >= grade[[g]]] <- g
}

まず、全員にNAを設定します。次に、0点以上の人にDを設定します。次に、50点以上の人にCをします。次に、65点以上の人にBを設定します。最後に、80点以上の人にAを設定します。成績85点の人のレコードはNA→D→C→B→Aと5回代入され、60点の人はNA→D→Cと3回代入されます。
なお、同じ日本語の列名を何回も書くのが冗長なので、変数を通して指定しています。列名に括弧を入れやすくなる御利益もあります。

## 相対評価（偏差値）##
mu <- mean(df04$合計点数, na.rm = TRUE)
sigma <- sd(df04$合計点数, na.rm = TRUE)
label_sigma <- "相対評価（偏差値）"
df04[, label_sigma] <- rep(NA, nrow(df04))
# 偏差値50がC、55がB、60がA
grade <- list("D"=0, "C"=mu, "B"=mu + sigma/2, "A"=mu + sigma) 
for(g in names(grade)){
    df04[, label_sigma][df04$合計点数 >= grade[[g]]] <- g
}

絶対評価と基本的には同じですが、平均値と標準偏差から評価の閾値を計算しています。

## 相対評価（順位）##
rank <- rank(-df04$合計点数, ties.method = "min", na.last = "keep")
label_rank <- "相対評価（順位）"
df04[, label_rank] <- rep(NA, nrow(df04))
 # 10位までがA、30位までがB、50位までがC、他はD
grade <- list("D"=Inf, "C"=50, "B"=30, "A"=10)
for(g in names(grade)){
    df04[, label_rank][rank <= grade[[g]]] <- g
}

順位の場合は、合計点数に-1をかけた値を`rank`に計算させています。`ties.method = "min",`はゴルフのような順位のつけ方にしろと言う指定で、`na.last`は欠損値を除外しないようにする指定です。また、点数と違って低いほうがよいので、代入時の比較演算子が以上から以下に変更しています。
`tail(df04, 3)`で末尾3行をみると、以下のようになっています。

学部 学籍番号 出席回数 期末試験 レポート評価 レポート点数 合計点数 絶対評価 相対評価（偏差値） 相対評価（順位）
文学部 L0051 12 86 E 0 86 A B B
文学部 L0052 12 93 D 5 98 A A A
文学部 L0053 11 26 B 15 41 D D D

成績の保存

df04を保存します。

write.table(df04, file = "成績.csv", sep = ",", row.names = FALSE, fileEncoding = "utf-8")

列を絞るときは、こんな感じで。

write.table(df04[, c("学部", "学籍番号", "合計点数", label_absolute, label_sigma, label_rank)], 
    file = "成績.csv", sep = ",", row.names = FALSE, fileEncoding = "utf-8")

成績の検索

特定の人や集団を絞り込みたいときは`subset`や`grep`を使うと便利です。なおデータフレームと括弧で絞り込むときには、括弧の中のカンマの位置に注意してください。カンマの左側が行番号/行名、右側が列番号/列名になるのですが、カンマは必ず必要になります。

# 学籍番号で絞る
subset(df04, 学籍番号=="L0021")
# 正規表現で学籍番号の下2桁が21で絞る
df04[grep("21$", df04$学籍番号), ]
# 絶対評価Aで絞る
subset(df04, 絶対評価=="A")
# 文学部のB評価を抽出する
subset(df04, 絶対評価=="B" & 学部 == "文学部")
# 期末試験が未受で絞る
df04[is.na(df04$期末試験), ]

集計

データフレームを整備しておけば、集計は簡単です。`summary(df04)`とすると全体の集計が出たりしますが、大学の講義は他学部の学生の履修などもあるので、学生の属性ごとの平均点を見たくなったりするかも知れません。これも`aggregate`や`xtabs`で簡単に集計できます。

学部ごとの平均点

# 学部ごとの平均点（meanをminやmaxに置き換えると、最低点や最高点になる）
aggregate(期末試験 ~ 学部, data = df04, mean, na.rm = TRUE)
学部 期末試験
現代社会学部 50.05
文学部 56.96

複数の列にそれぞれ学部別の集計をとることもできます。

aggregate(cbind(期末試験, レポート点数, 合計点数) ~ 学部, data = df04, mean, na.rm = TRUE)
学部 期末試験 レポート点数 合計点数
現代社会学部 50.05 9.5 59.25
文学部 56.96 8.2 64.56

学部ごとの最低点、平均点、最高点

複数の`aggregate`が戻してくるデータフレームの列を`cbind`で連結し、一つの表にすることもできます。

with(df04, {
    model <- 期末試験 ~ 学部
    a <- cbind(
        aggregate(model, data = df04, min, na.rm = TRUE), 
        aggregate(model, data = df04, mean, na.rm = TRUE)[, 2], 
        aggregate(model, data = df04, max, na.rm = TRUE)[, 2])
    colnames(a) <- c("学部", "最低点", "平均点", "最高点")
    a
})
学部 最低点 平均点 最高点
現代社会学部 14 50.05 97
文学部 12 56.96 97

学部ごとのレポート評価の頻度

# 学部ごとのレポート評価の頻度
xtabs(~ レポート評価 + 学部, data = df04)
レポート評価 現代社会学部 文学部
A 7 6
B 13 16
C 5 4
D 16 13
E 6 14

プロット

中央値と第1、第3四分位点と最大値と最小値とはずれ値の図示であるboxplotで視覚的に比較する方法もあります。

boxplot(合計点数 ~ 学部, data=df04, varwidth = TRUE)

`varwidth = TRUE`はサブサンプルのサイズごとに棒の太さを変える指定です。

barplotの例

ヒストグラムはこんな感じで。

# 10点刻み
hist(df04$合計点数, breaks=seq(0, 100, 10), main="絶対評価の分布")
# D, C, B, Aの閾値刻み
# hist(df04$合計点数, breaks=c(0, 50, 65, 80, 100), main="絶対評価の分布")
histの例


なお、画像ファイルとして保存する場合は、例えば以下のようにします。

png(filename = "学部別合計得点.png",  width = 640, height = 640, bg="white", type="cairo")
# Windowsの場合
windowsFonts(Meiryo = windowsFont("Meiryo"))
par(family="Meiryo")
# MacOSの場合
# par(family= "HiraKakuProN-W3")
boxplot(合計点数 ~ 学部, data=df04, varwidth = TRUE)
dev.off()

フォント指定が煩雑ですが、日本語フォント指定をしないと文字化けします。

計量分析

差の検定

統制変数が入れられないので、社会科学方面よりも理工系のラボで親しまれているものですが、試験の点数の比較などには有益です。

序数などの頻度の検定

点数ではなくて、A、B、C、D、Eの評価なのでクロス表の検定を行います。

# 学部ごとのレポート評価に差があるか比較しましょう
ct <- xtabs(~ レポート評価 + 学部, data = df04)
chisq.test(ct, simulate.p.value = TRUE)
# セルごとにp値を計算
round(2*pnorm(abs(chisq.test(ct, simulate.p.value = TRUE)$stdres), lower.tail=FALSE), 4)

なお、今回のデータでは統計的に有意ではありませんでした。

基数の差の検定

学部ごとの平均点に差があるかを比較しましょう。統計学の教科書的にはウェルチのt検定を用いたくなりますが、3つ以上の部分集合を比較する場合は別に多重比較検定が必要になり、また図示の方法を考えなければなりません。理工系のラボではウェルチのt検定と信頼区間の図が定番なのですが、二つのサブサンプルの信頼区間のエラーバーが重なっていても、ウェルチのt検定で差が統計的に有意であることがあり、あまりよい慣習とは言えません。広く知られた方法ではないですが、ガブリエル比較区間を用います。

library(rgabriel)

# 以下の4行以外はだいたいコピペで使えると思います
r_rgabriel <- rgabriel(df04$合計点数, df04$学部, 0.025)
r_aggregate <- aggregate(合計点数 ~ 学部, data = df04, mean, na.rm = TRUE)
clst <- unique(df04$学部)
noc <- length(clst)

means <- r_aggregate[ ,2]

gci <- sapply(1:noc, function(i){
    c(means[i] - r_rgabriel[i], means[i] + r_rgabriel[i], use.names = FALSE)
})
dimnames(gci) <- list(c("lower", "upper"), clst)

dotchart(means, xlim = c(0.9, 1.1)*range(gci),
    main = "95%ガブリエル比較区間", cex = 1.25, labels = colnames(gci),
    pch = 21, pt.cex = 2, xlab = "点数",  
    bg = "white", color = "black", lcolor = "darkgray", panel.first = {
        arrows(gci[1, ], 1:noc, gci[2, ], 1:noc, 
            length = 0.1, angle = 90, code = 3)
    })
rgabrielの例

ヒゲ部分が重なっているので、二つのサブサンプルに統計的に有意な差が無いことが分かります。

回帰分析

図書館利用頻度やアルバイトの時間などの学生ごとのデータが取れていれば、成績との関係を回帰分析などで調べることもできます。アンケート結果のファイルが別にあるとした分析例を示します。

被説明変数が基数データ
# ファイルロード後、即結合している
df05 <- merge(df04, read.table("http://wh.anlyznews.com/R/dataset/grading/lifestyle.csv", 
    header = TRUE, sep = ",", fileEncoding = "utf-8"), by = c("学部", "学籍番号"))
r_lm <- lm(期末試験 ~ 出席回数 + 図書館利用頻度 + アルバイト時間, data = df05)
summary(r_lm)
lmの例

アルバイト時間の影響が大きそうです。プロットもしてみましょう。

# 出席回数と図書館利用頻度が平均であった場合に補正した試験点数を計算
補正期末試験 <- with(df05, 期末試験 - coef(r_lm)[2]*(出席回数 - mean(出席回数)) - coef(r_lm)[3]*(図書館利用頻度 - mean(図書館利用頻度)))
# プロットする範囲
xlim <- c(0, max(df05$アルバイト時間))
# 予測値を計算
px <- seq(xlim[1], xlim[2], length.out=100)
prdct <- predict(r_lm, data.frame(
    アルバイト時間 = px,
    出席回数 = mean(df05$出席回数),
    図書館利用頻度 = mean(df05$図書館利用頻度)), interval = "prediction")
# プロットする
plot(df05$アルバイト時間, 補正期末試験, xlim = xlim, panel.first = {
    # 95%予測区間をピンク斜線で描く
    px <- c(px, rev(px))
    py <- c(prdct[,2], rev(prdct[,3]))
    polygon(px, py, density=25, col="pink", border=NA)
}, pch = 21, bg = "white")
# 予測線を描く
lines(xlim, c(prdct[1,1], prdct[nrow(prdct),1]), lty = 3)
# 凡例を書く
legend("topright", c("観測値", "予測値"), 
    lty = c(NA, 3), 
    seg.len = 3, 
    col = "black",
    pch = c(21, NA), 
    bty = "n")
lmからのplotの例
被説明変数が序数データ

レポート評価など基数とみなせない場合は、順序ロジットを使います。この推定はMASSパッケージを使うのが手っ取り早いです。

# "A"から"D"の評価に1から4の数字を割り当て、順序属性を与える
# Eは0人なので除外
順序化レポート評価 <- ordered(df05$レポート評価, c("D", "C", "B", "A"))

# 順序ロジット回帰をする
library(MASS)
r_polr <- polr(順序化レポート評価 ~ 出席回数 + 図書館利用頻度 + アルバイト時間, data = df05, Hess = TRUE)
summary(r_polr) # 推定結果
polrの例

なお推定された係数に観測値をかけて足したものは、BからEの学生 iの評価 jの生起確率 \pi_{ij}自体ではなく、 log(\pi_{ij}/(1-\pi_{ij}))になるので注意してください。

予測値は以下のように求められます。

# 各学生のレポート評価を予測
predict(r_polr)
# 各学生のレポート評価の確率を"A"から"D"まで予測
predict(r_polr, type="probs")

解釈がしづらいので、限界効果を見ましょう。ererパッケージを使うので、入っていない場合は`install.packages(c("erer"))`しておいてください。

library(erer)
ocME(r_polr) # 限界効果の計算
ocMEの例

説明変数が平均値と同じ学生の図書館利用頻度を1単位増やすと、Aの可能性が2.2%、Bの可能性が25.9%ぐらい増え、他の評価の可能性が減ります。

ところで、ここで用いたererパッケージは本日時点でCRAN checks: erer results [issues need fixing before 2024-06-29]と出ていて、来週あたりにインストールしづらくなるかも知れません。追加のパッケージを使わず限界効果を計算する方法を別エントリーに書きました。追記(2024/07/02 10:31):他を試してみたのですが、代替手段としてはmarginaleffectsパッケージが良さそうでした。marginaleffects::avg_slopesで平均限界効果が簡単に出ます。ただし平均におけるmarginalな効果は、ダミー変数があるとうまく出せそうになかったです。

図が欲しい場合は、effectsパッケージを使うと手ごろです。

library(effects)
r_pe <- predictorEffect("図書館利用頻度", r_polr)
plot(r_pe)
effects::predictorEffectの例

図書館利用頻度以外を平均値だと仮定したときの、図書館利用頻度が評価A、B、C、Dのそれぞれの確率に与える影響を見ています。限界効果と言うのは、この図の線の傾きのことになります。

おまけ（tidyverse）

Rの標準構成だと以上の操作でよいのですが、最近はtidyverseと呼ばれるパッケージ群を使ってデータ処理を行うことが一般化しています。こちらのコードも記しておきます。

# tidyverseが入っているか確認してロード
loadpkg <- function(pkgs){
	for(pkg in pkgs){
		if(!any(suppressWarnings(library(quietly=TRUE, verbose=FALSE)$results[,"Package"] == pkg))){
			stop(sprintf("Do install.packages(\"%s\") before runnning this script.", pkg))
		}
		library(pkg, character.only = TRUE)
	}
}
loadpkg("tidyverse")

df01 <- read_csv('http://wh.anlyznews.com/R/dataset/grading/list.csv')
df02 <- read_csv('http://wh.anlyznews.com/R/dataset/grading/exam.csv')
df03 <- read_csv('http://wh.anlyznews.com/R/dataset/grading/report.csv')

df04 <- df01 %>% left_join(df02) %>% left_join(df03) %>% mutate(
    レポート点数 = unlist(sapply(with(., {
# 未提出者はE評価/レポート未提出による未受は無し
        tmp <- レポート評価
        tmp[is.na(レポート評価)] <- "E"
        tmp
    }), switch, "A"=20, "B"=15, "C"=10, "D"=5, "E"=0))) %>% mutate(
        合計点数 = pmin(with(., 期末試験 + レポート点数), 100)
    ) %>% mutate("絶対評価" = with(., {
    grade <- list("D"=0, "C"=50, "B"=65, "A"=80)
    v <- rep(NA, nrow(.))
    for(g in names(grade)){
        v <- case_when(合計点数 >= grade[[g]] ~ g, .default = v)
    }
    v
})) %>% mutate("相対評価（偏差値）" = with(., {
    mu <- mean(合計点数, na.rm = TRUE)
    sigma <- sd(合計点数, na.rm = TRUE)
    grade <- list("D"=0, "C"=mu, "B"=mu + sigma/2, "A"=mu + sigma) 
    v <- rep(NA, nrow(.))
    for(g in names(grade)){
        v <- case_when(合計点数 >= grade[[g]] ~ g, .default = v)
    }
    v
})) %>% mutate("相対評価（順位）" = with(., {
    rank <- rank(-合計点数, ties.method = "min", na.last = "keep")
    # 10位までがA、30位までがB、50位までがC、他はD
    grade <- list("D"=Inf, "C"=50, "B"=30, "A"=10)
    v <- rep(NA, nrow(.))
    for(g in names(grade)){
        v <- case_when(rank <= grade[[g]] ~ g, .default = v)
    }
    v
}))

# 集計結果を保存
write_csv(df04, "成績.csv")

# 集計はaggregateを使って変わらないので略

# 学籍番号で検索
filter(df04, 学籍番号=="L0021")
# 正規表現で検索
filter(df04, str_detect(学籍番号, "21$"))
# 絶対評価Aで絞る
filter(df04, 絶対評価=="A")
# 文学部のB評価を抽出する
filter(df04, 絶対評価=="B" & 学部 == "文学部")
# 期末試験が未受で絞る
filter(df04, is.na(期末試験))

# プロットはバイオリンを用いる
windowsFonts(Meiryo = windowsFont("Meiryo"))
ggplot(df04, aes(x=学部, y=合計点数, fill=学部)) + geom_violin() + theme(text = element_text(family="Meiryo"))
# 保存するときは以下
# ggsave(file = "学部別分布.png", type = "cairo", dpi = 150, width = 6, height = 3)
violin plotの例

回帰分析とプロット

ggplot2だとggeffectsパッケージを使うと手軽に信頼区間や予測区間が描けます。

loadpkg("ggeffects")

# アンケートデータをダウンロードして結合
df05 <- df04 %>% right_join(read_csv('http://wh.anlyznews.com/R/dataset/grading/lifestyle.csv'))
# 回帰分析
r_lm <- lm(期末試験 ~ 出席回数 + 図書館利用頻度 + アルバイト時間, data = df05)

# 予測区間を求める
gprdct <- ggpredict(r_lm, terms = c("アルバイト時間"), interval = "prediction")

# 内部コードがUTF-8ではない古いOSもしくはRをお使いで、
# gsub("~", "", trim_ws(unlist(strsplit(split = pattern, x = .x,  でエラー: input string 3 is invalid in this locale
# などとエラーが出る場合は、上1行をコメントアウトして、下2行コメントイン
# a <- df05$アルバイト時間
# gprdct <- ggpredict(lm(期末試験 ~ 出席回数 + 図書館利用頻度 + a, data = df05), terms = "a", interval = "prediction")

# プロットする
# windowsFonts(Meiryo = windowsFont("Meiryo"))
ggplot() + 
    geom_point(data = df05, aes(x = アルバイト時間, y = 期末試験)) + 
    geom_line(data = gprdct, aes(x = x, y = predicted)) + 
    geom_ribbon(data = gprdct, aes(x = x, ymin = conf.low, ymax = conf.high), alpha = 0.1) # + theme(text = element_text(family="Meiryo"))
# ggsave(file = "r_lm_ggeffects.png", type = "cairo", dpi = 150, width = 6, height = 3)
ggeffectsの利用例

まとめ

とりあえずありそうな例を並べてみました。抜けているタスクがある場合は、御指摘いただけると、加筆するかも知れません。