7章_相関分析
1 相関関係
1.1 ファイルの読み込み
dat <- read.csv("data/data_ch7-1.csv", header = TRUE)
# 1～6列目のうちNAのある行を削除
dat.2 <- na.omit(dat[, 1 : 6])
# NAのある行を削除したデータの冒頭を確認
head(dat.2)| student | class | sex | faculty | mid | end |
|---|---|---|---|---|---|
| S001 | A | M | F01 | 86 | 99 |
| S002 | A | F | F01 | 96 | 82 |
| S003 | A | M | F01 | 52 | 77 |
| S004 | A | F | F01 | 72 | 47 |
| S005 | A | F | F01 | 74 | 67 |
| S006 | A | M | F01 | 90 | 84 |
# midとend（5～6列目）の記述統計
summary(dat.2[, 5 : 6])## mid end
## Min. : 24.00 Min. : 8.00
## 1st Qu.: 60.00 1st Qu.: 59.00
## Median : 72.00 Median : 78.00
## Mean : 70.63 Mean : 73.58
## 3rd Qu.: 82.00 3rd Qu.: 94.00
## Max. :100.00 Max. :100.001.2 中間試験と期末試験の積率相関係数
cor(dat.2$mid, dat.2$end, method = "pearson")## [1] 0.6186017## 無相関検定と95%信頼区間
cor.test(dat.2$mid, dat.2$end, method = "pearson")##
## Pearson's product-moment correlation
##
## data: dat.2$mid and dat.2$end
## t = 21.79, df = 766, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5729303 0.6604433
## sample estimates:
## cor
## 0.61860171.3 対角線ありの散布図の作成
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(dat.2$mid, dat.2$end,
# ドットの大きさを変更
cex = 1.2,
# x軸・y軸のラベルと範囲を変更，タイトルを表示
xlab = "中間テスト得点", xlim = c(0, 100),
ylab = "期末テスト得点", ylim = c(0, 100))
# 散布図に切片a = 0，傾きb = 1の直線を表示
abline(a = 0, b = 1, lty = 1)
1.4 中間試験と期末試験の散布図
pairs.panels(dat.2[, 5 : 6],
# 回帰直線（黒色）とその95%信頼区間を表示
lm = TRUE, ci = TRUE, col = "black",
# ヒストグラムの色を変更
hist.col = "grey",
# ドットのスタイルを変更
pch = 21,
# 有意な相関係数にアスタリスクを付ける
stars = TRUE,
# 平滑線（smooth）と相関円（ellipses）を描写しない
smooth = FALSE, ellipses = FALSE)
1.4.1 学科F01だけで相関分析
cor(dat.2$mid[dat.2$faculty == "F01"], dat.2$end[dat.2$faculty == "F01"], method = "pearson")## [1] 0.59997171.4.2 学科別にプロット
xyplot(mid ~ end | faculty, data = dat.2,
# 中間の点数が80点以上か否かで色分け，異なるマーカーを使用
groups = mid >= 80, col = c("grey20", "grey40"),
pch = c(1, 16), cex = c(1.0, 1.0),
# 中間試験と期末試験の平均点ラインを描写
abline = list(h = mean(dat.2$mid),v = mean(dat.2$end),
col = "grey"))
1.4.3 F01~F09の並び替え
dat.2$faculty <- factor(dat.2$faculty,
levels = c("F07", "F08", "F09", "F04", "F05", "F06",
"F01", "F02", "F03"))
xyplot(mid ~ end | faculty, data = dat.2,
groups = mid >= 80, col = c("grey20", "grey40"),
pch = c(1, 16), cex = c(1.0, 1.0),
abline = list(h = mean(dat.2$mid),v = mean(dat.2$end),
col = "grey"))
2 授業評価アンケートの分析
2.1 Spearmanの相関係数
# 7～9列目のうちNAのある行を削除
dat.3 <- na.omit(dat[, 7 : 9])
# Spearmanの相関係数を計算
cor(dat.3, method = "spearman")## suki manzoku rikai
## suki 1.0000000 0.1754352 0.4985831
## manzoku 0.1754352 1.0000000 0.5177924
## rikai 0.4985831 0.5177924 1.00000002.2 相関係数の有意確率と信頼区間
res <- corr.test(dat.3, method = "spearman")
print(res, short = FALSE)## Call:corr.test(x = dat.3, method = "spearman")
## Correlation matrix
## suki manzoku rikai
## suki 1.00 0.18 0.50
## manzoku 0.18 1.00 0.52
## rikai 0.50 0.52 1.00
## Sample Size
## [1] 790
## Probability values (Entries above the diagonal are adjusted for multiple tests.)
## suki manzoku rikai
## suki 0 0 0
## manzoku 0 0 0
## rikai 0 0 0
##
## Confidence intervals based upon normal theory. To get bootstrapped values, try cor.ci
## raw.lower raw.r raw.upper raw.p lower.adj upper.adj
## suki-manzk 0.11 0.18 0.24 0 0.11 0.24
## suki-rikai 0.44 0.50 0.55 0 0.44 0.56
## manzk-rikai 0.46 0.52 0.57 0 0.45 0.582.3 散布図
plot(dat.3)
2.4 クロス集計表の作成
ques <- xtabs( ~ manzoku + rikai, data = dat.3)
ques## rikai
## manzoku 1 2 3 4 5
## 1 28 18 26 25 0
## 2 21 37 29 27 0
## 3 25 21 35 18 0
## 4 76 65 71 75 0
## 5 0 0 0 97 96# パッケージの読み込み
library("gplots")##
## 次のパッケージを付け加えます: 'gplots'## 以下のオブジェクトは 'package:stats' からマスクされています:
##
## lowess# 気球グラフのプロット
balloonplot(ques)
3 テストの妥当性と信頼性
3.1 ファイルの読み込み
dat.4 <- read.csv("data/data_ch7-2.csv", header = TRUE)
head(dat.4)| student | item1 | item2 | item3 | item4 | item5 |
|---|---|---|---|---|---|
| S001 | 1 | 1 | 1 | 1 | 0 |
| S002 | 1 | 0 | 0 | 1 | 1 |
| S003 | 0 | 1 | 0 | 1 | 1 |
| S004 | 1 | 1 | 1 | 1 | 1 |
| S005 | 0 | 0 | 0 | 0 | 0 |
| S006 | 1 | 1 | 0 | 0 | 1 |
library("tidyr")
d.gathered <- tidyr::gather(data = dat.4,key="item",value="score",2:6) # 列数指定
str(d.gathered)## 'data.frame': 50 obs. of 3 variables:
## $ student: chr "S001" "S002" "S003" "S004" ...
## $ item : chr "item1" "item1" "item1" "item1" ...
## $ score : int 1 1 0 1 0 1 1 0 1 1 ...attach(d.gathered)
res <- aov(score ~ factor(student) + item)
summary(res)## Df Sum Sq Mean Sq F value Pr(>F)
## factor(student) 9 5.6 0.6222 3.733 0.00215 **
## item 4 0.4 0.1000 0.600 0.66503
## Residuals 36 6.0 0.1667
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1str(summary(res))## List of 1
## $ :Classes 'anova' and 'data.frame': 3 obs. of 5 variables:
## ..$ Df : num [1:3] 9 4 36
## ..$ Sum Sq : num [1:3] 5.6 0.4 6
## ..$ Mean Sq: num [1:3] 0.622 0.1 0.167
## ..$ F value: num [1:3] 3.73 0.6 NA
## ..$ Pr(>F) : num [1:3] 0.00215 0.66503 NA
## - attr(*, "class")= chr [1:2] "summary.aov" "listof"summary(res)[[1]]["Df"]#自由度の取り出し| Df | |
|---|---|
| factor(student) | 9 |
| item | 4 |
| Residuals | 36 |
ff <- summary(res)[[1]]["F value"][1,]#F値の取り出し
1-1/ff## [1] 0.7321429生徒要因のF値であるffから算出した 1-1/ff
はクロンバックのα係数に等しい。
詳細はクロンバックの
α 信頼性係数について(訂正第2版) を参照
このことは、生徒による違いを識別することができることを意味している。
aovのsummaryからF値を抽出する方法は、分散分析の関数aov()が返す結果オブジェクトの中味
を参照。
3.2 信頼性係数
信頼性係数の理解については次のサイトが役に立った。
クロンバックのα係数を手計算してみる。
attach(dat.4)## 以下のオブジェクトは d.gathered からマスクされています:
##
## studentdat.5 <- dat.4 %>%
mutate( T = item1+item2+item3+item4+item5 )
str(dat.5)## 'data.frame': 10 obs. of 7 variables:
## $ student: chr "S001" "S002" "S003" "S004" ...
## $ item1 : int 1 1 0 1 0 1 1 0 1 1
## $ item2 : int 1 0 1 1 0 1 1 0 1 1
## $ item3 : int 1 0 0 1 0 0 1 0 1 1
## $ item4 : int 1 1 1 1 0 0 1 0 0 1
## $ item5 : int 0 1 1 1 0 1 0 0 0 1
## $ T : int 4 3 3 5 0 3 4 0 3 5var(dat.4[,-1])## item1 item2 item3 item4 item5
## item1 0.23333333 0.12222222 0.16666667 0.08888889 0.05555556
## item2 0.12222222 0.23333333 0.16666667 0.08888889 0.05555556
## item3 0.16666667 0.16666667 0.27777778 0.11111111 -0.05555556
## item4 0.08888889 0.08888889 0.11111111 0.26666667 0.11111111
## item5 0.05555556 0.05555556 -0.05555556 0.11111111 0.27777778sum(var(dat.4[,-1]))## [1] 3.111111var(dat.5$T)## [1] 3.111111b <- var(dat.4[,2])+var(dat.4[,3])+var(dat.4[,4])+var(dat.4[,5])+var(dat.4[,6])
a <- 5/4*(1-(b/var(dat.5$T)))
a## [1] 0.7321429分散共分散行列の要素の合計がTの分散になっている。
3.3 Cronbachのα係数を計算
alpha(dat.4[, -1])##
## Reliability analysis
## Call: alpha(x = dat.4[, -1])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.73 0.74 0.84 0.36 2.8 0.14 0.6 0.35 0.38
##
## lower alpha upper 95% confidence boundaries
## 0.46 0.73 1
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## item1 0.63 0.64 0.76 0.31 1.8 0.20 0.082 0.38
## item2 0.63 0.64 0.76 0.31 1.8 0.20 0.082 0.38
## item3 0.68 0.68 0.64 0.35 2.1 0.17 0.014 0.36
## item4 0.67 0.68 0.78 0.34 2.1 0.18 0.111 0.37
## item5 0.79 0.80 0.76 0.49 3.9 0.11 0.020 0.47
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## item1 10 0.78 0.79 0.73 0.63 0.7 0.48
## item2 10 0.78 0.79 0.73 0.63 0.7 0.48
## item3 10 0.72 0.72 0.72 0.51 0.5 0.53
## item4 10 0.73 0.72 0.62 0.54 0.6 0.52
## item5 10 0.48 0.47 0.39 0.20 0.5 0.53
##
## Non missing response frequency for each item
## 0 1 miss
## item1 0.3 0.7 0
## item2 0.3 0.7 0
## item3 0.5 0.5 0
## item4 0.4 0.6 0
## item5 0.5 0.5 0信頼性係数の詳細が出力される。
# パッケージの読み込み
library("RMeCab")
# 分析データの読み込み
# data_ch7-3.csvを選択
dat.5 <- read.csv("data/data_ch7-3.csv", header = TRUE)
#dat.5 <- read.csv(file.choose(), header = FALSE)
# 形態素解析
RMC <- RMeCabDF(dat.5)
# 100名分の自由回答を結合
RMC.2 <- unlist(RMC)
# 結合したデータの冒頭を確認
head(RMC.2)## 名詞 助詞 名詞 助詞 名詞 動詞
## "レポート" "の" "テーマ" "が" "複数" "ある"# 形態素解析の結果をデータフレームに変換
RMC.3 <- data.frame(RMC.2, names(RMC.2))
# 列名を編集
colnames(RMC.3) <- c("Word", "POS")
# データフレームの冒頭を確認
head(RMC.3)| Word | POS |
|---|---|
| レポート | 名詞 |
| の | 助詞 |
| テーマ | 名詞 |
| が | 助詞 |
| 複数 | 名詞 |
| ある | 動詞 |
# データフレームから名詞の行だけを抽出
noun <- RMC.3[RMC.3$POS == "名詞", ]
# 抽出した結果（冒頭の3つ）を確認
head(noun, 3)| Word | POS | |
|---|---|---|
| 1 | レポート | 名詞 |
| 3 | テーマ | 名詞 |
| 5 | 複数 | 名詞 |
# 名詞の頻度を集計
noun.list <- table(noun[, 1])
# 頻度の高い順に並び替え
noun.list.2 <- sort(noun.list, decreasing = TRUE)
# 並び替え結果をデータフレームに変換
noun.list.3 <- data.frame(noun.list.2)
# 列名を編集
colnames(noun.list.3) <- c("Word", "Freq")
# 頻度上位の名詞（冒頭の3つ）を確認
head(noun.list.3, 3)| Word | Freq |
|---|---|
| 授業 | 26 |
| 先生 | 22 |
| 教室 | 19 |