立命館大学言語学研究会 (Day2)

 以下は、立命館大学言語学研究会のR講習会に関するメモです*1

記述統計

 まずは、平均値、中央値、分散、標準偏差などの記述統計を求めてみましょう。

x <- c(1, 2, 3, 4, 5)
# 平均値
mean(x)
# 中央値
median(x)

# グループAの5人の年収の調査(単位は万円)
a <- c(100, 200, 300, 400, 500)
mean(a)
median(a)
# グループBの5人の年収の調査(単位は万円)
b <- c(100, 200, 300, 400, 5000)
mean(b)
median(b)

x <- c(1, 2, 3, 4, 5)
# 最大値
max(x)
# 最小値
min(x)
# 分散
var(x)
# 標準偏差
sd(x)
# 要約統計量
summary(x)

# 行列の作成(復習を含む)
mat <- matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9), nrow = 3, ncol = 3, byrow = TRUE)
mat
# 行列の総和
sum(mat)
# 行列の平均値
mean(mat)
# 行列の1行目の総和
sum(mat[1, ])
# 行列の2〜3列目の平均値
mean(mat[, 2 : 3])
# 行ごとの総和
rowSums(mat)
# 列ごとの総和
colSums(mat)
# 行ごとの平均値
rowMeans(mat)
# 列ごとの平均値
colMeans(mat)

# 行ごとの総和(rowSums(mat)と同じ)
apply(mat, 1, sum)
# 列ごとの平均値(colMeans(mat)と同じ)
apply(mat, 2, mean)
# 行ごとの最大値
apply(mat, 1, max)
# 列ごとの要約統計量
apply(mat, 2, summary)

検定・効果量

 次に、検定をしてみましょう。

# クロス集計表の準備
tab <- matrix(c(96, 54, 52, 48), nrow = 2, ncol = 2, byrow = TRUE)
rownames(tab) <- c("Male", "Female")
colnames(tab) <- c("Jotai", "Keitai")
# クロス集計表の確認
tab
# クロス集計表の可視化
mosaicplot(tab)

# カイ自乗検定
chisq.test(tab, correct = FALSE)
# カイ自乗検定(イェーツの補正あり)
chisq.test(tab, correct = TRUE)

 2つ以上の頻度の差を検定する場合は、以下のようにします。多重比較の結果として得られるp値を解釈する際は、ボンフェローニの補正などをします。

# 2×3のクロス集計表の準備
tab.2 <- matrix(c(805, 414, 226, 99, 38, 12), nrow = 2, ncol = 3, byrow = TRUE)
rownames(tab.2) <- c("Correct", "Error")
colnames(tab.2) <- c("Level 1", "Level 2", "Level 3")
# クロス集計表の確認
tab.2

# 多重比較
# 1列目と2列目を検定
chisq.test(tab.2[, c(1, 2)], correct = FALSE)
# 1列目と3列目を検定
chisq.test(ab.2[, c(1, 3)], correct = FALSE)
# 2列目と3列目を検定
chisq.test(tab.2[, c(2, 3)], correct = FALSE)

 検定は非常に便利なものですが、検定結果がサンプルサイズに影響されることに注意する必要があります。そのため、p値だけでなく、オッズ比などの効果量も確認するとよいでしょう。

# 表中の数値を全て10倍
tab.3 <- tab * 10
# 10倍したデータの確認
tab.3
# カイ自乗検定
chisq.test(tab.3, correct = FALSE)

# オッズ比の計算
(tab[1, 1] / tab[2, 1]) / (tab[1, 2] / tab[2, 2])
# 10倍したデータでオッズ比を計算
(tab.3[1, 1] / tab.3[2, 1]) / (tab.3[1, 2] / tab.3[2, 2])

# 追加パッケージのインストール(初回のみ)
install.packages("vcd", dependencies = TRUE)
# 追加パッケージの読み込み(Rを起動するごとに毎回)
library(vcd)
# オッズ比の計算
oddsratio(tab, log = FALSE)
# オッズ比の信頼区間(下限値,上限値)の計算
confint(oddsratio(tab, log = FALSE))
  • 練習問題
    • 以下のデータを使って、モザイクプロットを描きなさい
    • 同じデータに対して、カイ自乗検定(イェーツの補正なし)を実行しなさい
    • 同じデータに対して、オッズ比と、オッズ比の信頼区間を計算しなさい
A新聞 B新聞
内閣支持 225 292
内閣不支持 275 208

相関・回帰

 続いて、相関回帰をしてみます。

# 追加パッケージのインストール(初回のみ)
install.packages("corpora", dependencies = TRUE)
# 追加パッケージの読み込み(Rを起動するごとに毎回)
library(corpora)
# データセットの準備
data(BNCbiber)
# データの冒頭の5行のみを表示
head(BNCbiber, 5)
# 相関係数の計算
cor(BNCbiber[, 2], BNCbiber[, 4], method = "pearson")

# 散布図の描画
plot(BNCbiber[, 2], BNCbiber[, 4], xlab = "past tense", ylab = "present tense")

# 無相関検定
cor.test(BNCbiber[, 2], BNCbiber[, 4])

# スピアマンの順位相関係数の計算
cor(BNCbiber[, 2], BNCbiber[, 4], method = "spearman")

# ピアソンの積率相関係数
cor(BNCbiber[, 2 : 4], method = "pearson")
# スピアマンの順位相関係数の計算
cor(BNCbiber[, 2 : 4], method = "spearman")

# 追加パッケージのインストール(初回のみ)
install.packages("psych", dependencies = TRUE)
# 追加パッケージの読み込み(Rを起動するごとに毎回)
library(psych)
# 相関係数が表示された散布図行列の作成
pairs.panels(BNCbiber[, 2 : 4])

# 単回帰分析
lm.result <- lm(BNCbiber[, 2] ~ BNCbiber[, 4])
# 結果の確認
lm.result

# 回帰式の可視化
plot(BNCbiber[, 4], BNCbiber[, 2], xlab = "present tense", ylab = "past tense", pch = 16, col = "grey")
abline(lm.result)

# 重回帰分析
lm.result.2 <- lm(BNCbiber[, 2] ~ BNCbiber[, 3] + BNCbiber[, 4])
# 結果の確認
lm.result.2
  • 練習問題
    • carsデータセットを使って、散布図を描きなさい
    • 同じデータを使って、自動車の速度とブレーキを踏んでから止まるまでの距離の積率相関係数と順位相関係数を求めなさい
    • 同じデータを使って、説明変数を速度、目的変数を距離とする回帰分析を行い、その結果を可視化しなさい

多変量解析

 多変量解析と呼ばれる手法はたくさんありますが、ここでは対応分析クラスター分析を扱います。

# 追加パッケージのインストール(初回のみ)
install.packages("ca", dependencies = TRUE)
# 追加パッケージの読み込み(Rを起動するごとに毎回)
library(ca)
# データセットの準備
data(author)
# データの冒頭の5行のみを表示
head(author, 5)

# 対応分析
ca.result <- ca(author)
# 結果の可視化
plot(ca.result)

# 行データ(テキスト)のみを表示
plot(ca.result, what = c("all", "none"))
# 列データ(変数)のみを表示
plot(ca.result, what = c("none", "all"))

# 対応分析から得られた詳しい結果の確認
ca.result
# 行データの表示(第1~2次元のみ)
ca.result$rowcoord[, 1 : 2]
# 行データの第1次元の得点を並び換え
sort(ca.result$rowcoord[, 1], decreasing = TRUE)
# 行データの第2次元の得点を並び換え
sort(ca.result$rowcoord[, 2], decreasing = TRUE)
# 列データの表示(第1~2次元のみ)
ca.result$colcoord[, 1 : 2]
# 列データの第1次元の得点を並び換え
sort(ca.result$colcoord[, 1], decreasing = TRUE)
# 列データの第2次元の得点を並び換え
sort(ca.result$colcoord[, 2], decreasing = TRUE)

# 相対頻度の計算
author.2 <- author / apply(author, 1, sum)
# ユークリッド距離の計算
dist.result <- dist(author.2, method = "euclidean")
# 階層型クラスター分析(ウォード法)
hclust.result <- hclust(dist.result, method = "ward.D2")
# 結果の可視化
plot(hclust.result)

# データセットの転置
author.3 <- t(author.2)
# ユークリッド距離の計算
dist.result.2 <- dist(author.3, method = "euclidean")
# 階層型クラスター分析(ウォード法)
hclust.result.2 <- hclust(dist.result.2, method = "ward.D2")
# 結果の可視化
plot(hclust.result.2)

# 階層型クラスターつきのヒートマップ
heatmap(author.2)

テキストマイニング

 最後に、RMeCabを使って、簡単なテキストマイニングデモをします。RMeCabのインストール方法は、他のパッケージと異なるため、公式サイトをよく読んでください(テキストデータの扱い方は、WindowsMacで異なる部分もあります)。

# 追加パッケージの読み込み(Rを起動するごとに毎回)
library(RMeCab)
# 短い文章の形態素解析
RMeCabC("すもももももももものうち")
# 形態素解析結果の保存
RMeCabC.result <- RMeCabC("すもももももももものうち")
# データ形式の確認
class(RMeCabC.result)
# データ形式の変換
RMeCabC.result.2 <- unlist(RMeCabC.result)
RMeCabC.result.2
# データのクラスの確認
class(RMeCabC.result.2)

# 解析結果の一部のみを表示
RMeCabC.result.2[1]
RMeCabC.result.2[2]
RMeCabC.result.2[1 : 3]

# 品詞情報のみを表示
names(RMeCabC.result.2)

# 単語の原形を復元
RMeCabC.result.3 <- RMeCabC("私は焼肉を食べていない。", 1)
RMeCabC.result.4 <- unlist(RMeCabC.result.3)
RMeCabC.result.4

# 追加パッケージのインストール(初回のみ)
install.packages("wordcloud", dependencies = TRUE)
# 追加パッケージの読み込み(Rを起動するごとに毎回)
library(wordcloud)

# RMeCabText関数で形態素解析(解析対象のファイルを選択)
RMeCabText.result <- RMeCabText(file.choose())
# RMeCabText関数の結果の確認
head(RMeCabText.result, 5)
# 単語ベクトルの作成
RMeCabText.result.2 <- unlist(sapply(RMeCabText.result, "[[", 1))
# 単語ベクトルの確認
head(RMeCabText.result.2, 5)
# ワードクラウドを描画
wordcloud(RMeCabText.result.2, min.freq = 2, random.order = FALSE)

# 追加パッケージのインストール(初回のみ)
install.packages("igraph", dependencies = TRUE)
# 追加パッケージの読み込み(Rを起動するごとに毎回)
library(igraph)
# NgramDFによる共起語の集計(解析対象のファイルを選択)
NgramDF.result <- NgramDF(file.choose(), type = 1, N = 2, pos = "名詞")
# 共起頻度2以上のペアのみを抽出
NgramDF.result.2 <- subset(NgramDF.result, Freq > 1)
# ネットワークの描画
g <- graph.data.frame(NgramDF.result.2, directed = FALSE)
plot(g, vertex.label = V(g)$name, vertex.color = "grey")
# 共起頻度3以上のペアのみを抽出
NgramDF.result.3 <- subset(NgramDF.result, Freq > 2)
# ネットワークの描画
g.2 <- graph.data.frame(NgramDF.result.3, directed = FALSE)
plot(g.2, vertex.label = V(g.2)$name, vertex.color = "grey")

参考図書

Rによるやさしいテキストマイニング

Rによるやさしいテキストマイニング

*1:本講習会は、学内のクローズドなイベントです

*2:CSVファイルを読み込む際、 ISO <- read.csv(file(file.choose(), encoding = "cp932"), header = TRUE, row.names = 1) のように文字コードで cp932 を指定すること