1 ファイルの読み込み

library(readr)
dat <- read.csv("data/data_ch4-1.csv", header = TRUE)
head(dat)
student class sex score
S001 A M 82
S002 A M 60
S003 A M 66
S004 A F 84
S005 A M 78
S006 A F 82 
# 行数と列数の確認
str(dat)
## 'data.frame':    71 obs. of  4 variables:
##  $ student: chr  "S001" "S002" "S003" "S004" ...
##  $ class  : chr  "A" "A" "A" "A" ...
##  $ sex    : chr  "M" "M" "M" "F" ...
##  $ score  : int  82 60 66 84 78 82 46 70 48 70 ...
# 欠損値（NA）の有無を確認
complete.cases(dat) 
##  [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [13]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [25]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [37]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [49]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE
## [61]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
subset(dat, complete.cases(dat) == FALSE)
student class sex score
59 S059 B F NA 
# 読み込んだ4列のデータのうちNAのある行を削除
dat.2 <- na.omit(dat)
# NAのある行を削除したデータの冒頭を確認
str(dat.2)
## 'data.frame':    70 obs. of  4 variables:
##  $ student: chr  "S001" "S002" "S003" "S004" ...
##  $ class  : chr  "A" "A" "A" "A" ...
##  $ sex    : chr  "M" "M" "M" "F" ...
##  $ score  : int  82 60 66 84 78 82 46 70 48 70 ...
##  - attr(*, "na.action")= 'omit' Named int 59
##   ..- attr(*, "names")= chr "59"
# クラス別の学習者数
table(dat.2$class)
## 
##  A  B 
## 33 37
# クラス別の男女の学習者数
table(dat.2$class, dat.2$sex)
##    
##      F  M
##   A 10 23
##   B 22 15

2 クラス別得点分布

# パッケージの読み込み
library("lattice")
library("beeswarm")
#ヒストグラム
histogram(~ score | class, data = dat.2)

#箱ひげ図と蜂群図
boxplot(dat.2$score ~ dat.2$class, ylim = c(0, 100), main = "Result of the Exam", xlab = "class", ylab = "score")
beeswarm(dat.2$score ~ dat.2$class, ylim = c(0, 100), pch = 16, add = TRUE)

# パッケージの読み込み
library("psych")
# クラスごとの記述統計量の計算
describeBy(dat.2$score, dat.2$class)
## 
##  Descriptive statistics by group 
## group: A
##    vars  n  mean    sd median trimmed   mad min max range  skew kurtosis   se
## X1    1 33 60.24 15.81     60   60.67 17.79  24  86    62 -0.21    -0.87 2.75
## ------------------------------------------------------------ 
## group: B
##    vars  n  mean    sd median trimmed   mad min max range  skew kurtosis   se
## X1    1 37 72.27 16.11     76   72.97 20.76  40  94    54 -0.33    -1.15 2.65

3 ２標本t検定

3.1 ルービンテスト

library("car")
# 等分散性の検定
leveneTest(dat.2$score, dat.2$class, center = mean)
## Warning in leveneTest.default(dat.2$score, dat.2$class, center = mean):
## dat.2$class coerced to factor.
Df F value Pr(>F)
group 1 0.1953908 0.6598699
68 NA NA 
# 独立したt検定
t.test(dat.2$score ~ dat.2$class)
## 
##  Welch Two Sample t-test
## 
## data:  dat.2$score by dat.2$class
## t = -3.1489, df = 67.363, p-value = 0.002443
## alternative hypothesis: true difference in means between group A and group B is not equal to 0
## 95 percent confidence interval:
##  -19.65119  -4.40450
## sample estimates:
## mean in group A mean in group B 
##        60.24242        72.27027

3.2 対応のあるt検定

4 ファイルの読み込み

dat.3 <- read.csv("data/data_ch4-2.csv", header = TRUE)
head(dat.3)
student sex pre post
S001 M 78 95
S002 F 76 69
S003 M 79 72
S004 F 62 75
S005 F 56 70
S006 M 76 79 
# preとpostの列に欠損値がないか確認（TRUEがあれば，欠損値がある）
is.na(dat.3[3 : 4])
##         pre  post
##  [1,] FALSE FALSE
##  [2,] FALSE FALSE
##  [3,] FALSE FALSE
##  [4,] FALSE FALSE
##  [5,] FALSE FALSE
##  [6,] FALSE FALSE
##  [7,] FALSE FALSE
##  [8,] FALSE FALSE
##  [9,] FALSE FALSE
## [10,] FALSE FALSE
## [11,] FALSE FALSE
## [12,] FALSE FALSE
## [13,] FALSE FALSE
## [14,] FALSE FALSE
## [15,] FALSE FALSE
## [16,] FALSE FALSE
## [17,] FALSE FALSE
## [18,] FALSE FALSE
## [19,] FALSE FALSE
## [20,] FALSE FALSE
## [21,] FALSE FALSE
## [22,] FALSE FALSE
## [23,] FALSE FALSE
## [24,] FALSE FALSE
## [25,] FALSE FALSE
## [26,] FALSE FALSE
## [27,] FALSE FALSE
## [28,] FALSE FALSE
## [29,] FALSE FALSE
## [30,] FALSE FALSE
table(is.na(dat.3[3 : 4]))
## 
## FALSE 
##    60
# 男女の学習者数
table(dat.3$sex)
## 
##  F  M 
## 18 12
# 事前・事後テストの結果を重ねたヒストグラム（col = rgbで半透明の指定）
hist(dat.3$pre, col = rgb(1, 0, 0, 0.5), xlim = c(30, 100), ylim = c(0, 15), main = "Overlapping Histogram", xlab = "score")
hist(dat.3$post, col = rgb(0, 0, 1, 0.5), add = TRUE)

# 箱ひげ図に個人の得点分布を重ねた蜂群図
(score <- c(dat.3$pre, dat.3$post))
##  [1] 78 76 79 62 56 76 68 59 65 43 71 72 50 64 64 65 66 75 66 67 57 68 75 56 80
## [26] 68 62 67 90 75 95 69 72 75 70 79 74 68 75 56 91 77 63 69 79 67 75 67 71 73
## [51] 55 81 87 66 79 77 83 77 84 79
(group <- factor(c(rep("pre", 30), rep("post", 30)), levels = c("pre","post")))
##  [1] pre  pre  pre  pre  pre  pre  pre  pre  pre  pre  pre  pre  pre  pre  pre 
## [16] pre  pre  pre  pre  pre  pre  pre  pre  pre  pre  pre  pre  pre  pre  pre 
## [31] post post post post post post post post post post post post post post post
## [46] post post post post post post post post post post post post post post post
## Levels: pre post
boxplot(score ~ group, ylim = c(0, 100), main = "Result of the Pre-Post Test", xlab = "test", ylab = "score")
beeswarm(score ~ group, ylim = c(0, 100), pch = 16, add = TRUE)

# 記述統計量の確認
describe(dat.3[, 3 : 4])
vars n mean sd median trimmed mad min max range skew kurtosis se
pre 1 30 67.33333 9.660918 67 67.58333 9.6369 43 90 47 -0.1971041 0.1970595 1.763834
post 2 30 74.43333 8.977686 75 74.41667 8.1543 55 95 40 0.0092125 0.0121618 1.639094 
# p値の指数表示を回避
options(scipen = 10)
2.813e-05
## [1] 0.00002813
# 事前テストと事後テストの相関係数
cor(dat.3$pre, dat.3$post)
## [1] 0.6487101

4.1 # 個別推移図（スパゲティ・プロット）

# パッケージの読み込み
library("lattice")
library("latticeExtra")
## 
## Attaching package: 'latticeExtra'
## The following object is masked from 'package:ggplot2':
## 
##     layer
# 個別推移図（スパゲティ・プロット）
dat.3
student sex pre post
S001 M 78 95
S002 F 76 69
S003 M 79 72
S004 F 62 75
S005 F 56 70
S006 M 76 79
S007 M 68 74
S008 M 59 68
S009 F 65 75
S010 F 43 56
S011 F 71 91
S012 M 72 77
S013 F 50 63
S014 F 64 69
S015 M 64 79
S016 F 65 67
S017 F 66 75
S018 M 75 67
S019 F 66 71
S020 M 67 73
S021 F 57 55
S022 F 68 81
S023 F 75 87
S024 F 56 66
S025 F 80 79
S026 M 68 77
S027 M 62 83
S028 F 67 77
S029 M 90 84
S030 F 75 79 
df <- data.frame(score, group);df
score group
78 pre
76 pre
79 pre
62 pre
56 pre
76 pre
68 pre
59 pre
65 pre
43 pre
71 pre
72 pre
50 pre
64 pre
64 pre
65 pre
66 pre
75 pre
66 pre
67 pre
57 pre
68 pre
75 pre
56 pre
80 pre
68 pre
62 pre
67 pre
90 pre
75 pre
95 post
69 post
72 post
75 post
70 post
79 post
74 post
68 post
75 post
56 post
91 post
77 post
63 post
69 post
79 post
67 post
75 post
67 post
71 post
73 post
55 post
81 post
87 post
66 post
79 post
77 post
83 post
77 post
84 post
79 post 
1 : nrow(dat.3)
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30
c(rep(1 : nrow(dat.3)), rep(1 : nrow(dat.3)))
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
## [51] 21 22 23 24 25 26 27 28 29 30
factor(c(rep(1 : nrow(dat.3)), rep(1 : nrow(dat.3))))
##  [1] 1  2  3  4  5  6  7  8  9  10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 1  2  3  4  5  6  7  8  9  10 11 12 13 14 15 16 17 18 19 20
## [51] 21 22 23 24 25 26 27 28 29 30
## 30 Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 ... 30
df$indiv <- factor(c(rep(1 : nrow(dat.3)), rep(1 : nrow(dat.3))))
#対応を示す識別番号を列に追加
df$indiv
##  [1] 1  2  3  4  5  6  7  8  9  10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 1  2  3  4  5  6  7  8  9  10 11 12 13 14 15 16 17 18 19 20
## [51] 21 22 23 24 25 26 27 28 29 30
## 30 Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 ... 30
dim(df)
## [1] 60  3
each <- xyplot(score ~ group, group = indiv, type = c("l"), data = df, xlab = "test", ylab = "score")
xyplot(score ~ group,data=df,type="l")

each

all_mean <- c(mean(dat.3$pre), mean(dat.3$post))
fact <- factor(c("pre", "post"), levels = c("pre","post"))
all <- xyplot(all_mean ~ fact, col = "black", lwd = 5, type = c("l"), data = df)
each + as.layer(all, axes = NULL)

# 横軸に事前テスト，縦軸に事後テストをプロットした散布図
plot(dat.3$pre, dat.3$post, las = 1, pch = 16, xlab = "pretest", ylab = "posttest", main = NA, xlim = c(0, 100), ylim = c(0, 100))
lines(par()$usr[1 : 2], par()$usr[3 : 4], lty = 3)