最近は理工系はもちろん、社会科学系の大学教員も研究で統計解析をしていることが多いので、Rで成績評価する方法を解説しても需要はあまり無さそうな気もするのですが、人文系でニーズがあるようなので。
名簿と試験の成績とレポートの評価の3つのCSVファイルがあるとして、それをつないで総合得点と成績を出し、基本統計量を出すところまでやってみましょう。なお、データは完全に架空のものです。
ファイルの読み込み
今回はウェブサイトに置いてあるファイルを読み込みますが、実際の場合はURLではなくファイル名を指定してください。なお、ファイルが読み込めないよと言うときは、`getwd()`で作業フォルダを確認し、RGuiのメニューからディレクトリの変更を行うか、`setwd("フォルダ名")`で移動することができるのを思い出しましょう。
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)
テキスト型のところは長さ、数値型のところは基本統計量が出ます。
以下では`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`はサブサンプルのサイズごとに棒の太さを変える指定です。
ヒストグラムはこんな感じで。
# 10点刻み hist(df04$合計点数, breaks=seq(0, 100, 10), main="絶対評価の分布") # D, C, B, Aの閾値刻み # hist(df04$合計点数, breaks=c(0, 50, 65, 80, 100), main="絶対評価の分布")
なお、画像ファイルとして保存する場合は、例えば以下のようにします。
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()
フォント指定が煩雑ですが、日本語フォント指定をしないと文字化けします。
計量分析
図書館利用頻度やアルバイトの時間などの学生ごとのデータが取れていれば、成績との関係を回帰分析などで調べることもできます。アンケート結果のファイルが別にあるとした分析例を示します。
被説明変数が基数データ
# ファイルロード後、即結合している 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)
アルバイト時間の影響が大きそうです。プロットもしてみましょう。
# 出席回数と図書館利用頻度が平均であった場合に補正した試験点数を計算 補正期末試験 <- 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")
被説明変数が序数データ
レポート評価など基数とみなせない場合は、順序ロジットを使います。この推定はMASSパッケージを使うのが手っ取り早いです。
# "A"から"E"の評価に1から5の数字を割り当て、順序属性を与える 順序化レポート評価 <- with(df05, { n <- c("A", "B", "C", "D", "E") r <- length(n):1 names(r) <- n as.ordered(r[レポート評価]) }) # 順序ロジット回帰をする library(MASS) r_polr <- polr(順序化レポート評価 ~ 出席回数 + 図書館利用頻度 + アルバイト時間, data = df05, Hess = TRUE) summary(r_polr) # 推定結果
なお推定された係数に観測値をかけて足したものは、BからEの学生の評価
の生起確率
自体ではなく、
になるので注意してください。
予測値は以下のように求められます。
# 各学生のレポート評価を予測 predict(r_polr) # 各学生のレポート評価の確率を"A"から"D"まで予測 predict(r_polr, type="probs")
解釈がしづらいので、限界効果を見ましょう。
library(erer) ocME(r_polr) # 限界効果の計算
説明変数が平均値と同じ学生の図書館利用頻度を1単位増やすと、Aの可能性が2.2%、Bの可能性が25.9%ぐらい増え、他の評価の可能性が減ります。
ところで、ここで用いたererパッケージは本日時点でCRAN checks: erer results [issues need fixing before 2024-06-29]と出ていて、来週あたりにインストールしづらくなるかも知れません。追加のパッケージを使わず限界効果を計算する方法を別エントリーに書きました。
おまけ（tidyverse）
Rの標準構成だと以上の操作でよいのですが、最近はtidyverseと呼ばれるパッケージ群を使ってデータ処理を行うことが一般化しています。こちらのコードも記しておきます。
# tidyverseが入っているか確認 if(!any(suppressWarnings(library(quietly=TRUE, verbose=FALSE)$results[,"Package"] == "tidyverse"))){ stop("Do install.packages(\"tidyverse\") before runnning this script.") } library(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)
回帰分析とプロット
ggplot2だとggeffectsパッケージを使うと手軽に信頼区間や予測区間が描けます。
# ggeffectsが入っているか確認 if(!any(suppressWarnings(library(quietly=TRUE, verbose=FALSE)$results[,"Package"] == "ggeffects"))){ stop("Do install.packages(\"ggeffects\") before runnning this script.") } library(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)
まとめ
とりあえずありそうな例を並べてみました。抜けているタスクがある場合は、御指摘いただけると、加筆するかも知れません。
