R を使った Resampling と作図

1. Resampling Methods

サンプルサイズが小さく，外れ値が含まれているデータを作って，ヒストグラムで分布を確認

``````# サンプルのデータを作る
sample <- c(8.72, 17.65, 13.53, 26.45, 24.34,
8.64, 12.06, 9.53, 9.42, 11, 10.18, 34.61, 12.41, 13.37)

# ヒストグラム
hist(sample, breaks = "FD", main = paste("Histogram of Sample Data"),
xlab = "Response Time", xlim = c(0, 40))

# ヒストグラム中の平均値に縦線を入れる場合
abline(v = mean(sample), col = "red", lwd = 2)
text(mean(sample), 5.5, "Mean = 15.13", cex = 2)
``````

ヒストグラムと箱ひげ図の両方を1つの図で見る

``````simple.hist.and.boxplot <- function(x) {
on.exit(par(op))
layout(matrix(c(1, 2), 2, 1), heights = c(3, 1))
par(mai = c(1.6, 2, 1, 1)/2)
hist(x, col = gray(0.95), xlim = c(0, 40), , xlab = "Response Time",
main = paste("Histogram of Sample Data"))
abline(v = mean(sample), col = "red", lwd = 2)
text(mean(sample), 5.5, "Mean = 15.13", cex = 1)
rug(x)
boxplot(x, horizontal = TRUE, ylim = range(c(0, 40)))
}
simple.hist.and.boxplot(sample)
``````

QQプロットで正規性の確認

``````# Q-Q plot
library(car)  # car パッケージを使用
qqPlot(sample, dist = "norm", labels = T)
``````

正規分布のサンプル

``````hist(rnorm(500, 0, 1), xlab = "", main = "", col = "lightgray",
border = "darkgray")
``````

``````norm.x = rnorm(500, 0, 1)
x = seq(-3.5, 3.5, 0.1)
dn = dnorm(x)
hist(norm.x, xlab = "", main = "Histogram of Normal Distribution",
col = "lightgray", border = "darkgray", prob = T)
lines(x, dn, col = "red", lwd = 2)
``````

Resampling のイメージを確認（スライド13）

``````dat <- c(1, 2, 3, 4, 5, 6, 7, 8, 9)
sample(dat, 9, replace = T)
``````
``````## [1] 5 2 8 8 5 9 5 7 8
``````
``````# ヒストグラムを確認する場合は hist(a) と書く
(a <- sample(dat, 9, replace = T))
``````
``````## [1] 7 5 6 8 3 5 2 8 5
``````
``````(b <- sample(dat, 9, replace = T))
``````
``````## [1] 3 9 3 8 9 9 2 3 9
``````
``````(c <- sample(dat, 9, replace = T))
``````
``````## [1] 6 2 9 3 6 6 5 6 5
``````
``````(d <- sample(dat, 9, replace = T))
``````
``````## [1] 9 9 9 1 2 2 5 3 7
``````
``````(e <- sample(dat, 9, replace = T))
``````
``````## [1] 5 5 2 1 9 7 1 7 9
``````

Resampling の例としてブートストラップ（10,000回）

``````x.boot <- numeric(10000)
for (i in 1:10000) {
x.boot[i] <- mean(sample(dat, 9, replace = T))
}

# ブートストラップのヒストグラム
hist(x.boot, col = gray(0.95))

# 平均値に線を入れる
abline(v = mean(x.boot), col = "red", lwd = 2)
``````

はじめの sample データをブートストラップ（2,000回）

``````sample <- c(8.72, 17.65, 13.53, 26.45, 24.34,
8.64, 12.06, 9.53, 9.42, 11, 10.18, 34.61, 12.41, 13.37)
x.boot <- numeric(2000)
for (i in 1:2000) {
x.boot[i] <- mean(sample(sample, 14, replace = T))
}

# ブートストラップのヒストグラム
hist(x.boot, breaks = "FD", ylim = c(0, 250),
main = paste("Histogram of Bootstrapped Mean"), xlab = "Response Time")

# 平均値に線を入れる
abline(v = mean(x.boot), col = "blue", lwd = 2)
text(mean(x.boot), 225, "Mean = 15.07", cex = 2)
``````

スライド29

``````# 2つのグループを作る
A <- c(20, 40, 30, 40, 40, 50)
B <- c(10, 70, 80, 90, 70, 80)

# 等分散を仮定の t 検定
t.test(A, B, var.equal = TRUE)
``````
``````##
##  Two Sample t-test
##
## data:  A and B
## t = -2.405, df = 10, p-value = 0.03698
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -57.79  -2.21
## sample estimates:
## mean of x mean of y
##     36.67     66.67
##
``````
``````# Welch の検定（等分散を仮定しない t 検定）
t.test(A, B)
``````
``````##
##  Welch Two Sample t-test
##
## data:  A and B
## t = -2.405, df = 6.269, p-value = 0.05115
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -60.2035   0.2035
## sample estimates:
## mean of x mean of y
##     36.67     66.67
##
``````
``````# マン・ホイットニーの U 検定
# 通常の wilcox.test(A,
#   B)はタイがあるため正確な値は求められない
# そのため，青木繁伸先生（群馬大学）のコードを使用
source("http://aoki2.si.gunma-u.ac.jp/R/src/U_test.R",
encoding = "euc-jp")
U.test(A, B, correct = FALSE)  # Wilcoxon検定
``````
``````##
##  マン・ホイットニーの U 検定
##
## data:  A and B
## U = 6.000, E(U) = 18.000, V(U) = 38.182, Z-value = 1.942, p-value
## = 0.05214
##
``````
``````# マン・ホイットニーの U 検定（正確確率検定版)
# 青木繁伸先生（群馬大学）のコードを使用
source("http://aoki2.si.gunma-u.ac.jp/R/src/exact-mw.R",
encoding = "euc-jp")
# 関数の定義が必要なので以下を入力
printf <- function(fmt, ...) {
cat(sprintf(fmt, ...))
}
# ----------
exact.mw(A, B, exact = TRUE)  # 正確確率検定
``````
``````## U = 6, Z = 1.94202, P 値 = 0.0521351
## 正確な P 値 = 0.05844155844
## 査察した分割表の個数は 226 個
``````
``````# 並べ替え検定（Permutation/Randomization Test）
library(exactRankTests)  # exactRankTests パッケージを使用
perm.test(A, B, conf.int = TRUE, exact = TRUE)
``````
``````##
##  2-sample Permutation Test
##
## data:  A and B
## T = 220, p-value = 0.05844
## alternative hypothesis: true mu is not equal to 0
## 95 percent confidence interval:
##  -53.33  -5.00
##
``````

2. Effect Size（効果量）

効果量が必要な理由をシミュレーションで証明

``````# 青木繁伸先生（群馬大学）のコードを参照
#   (http://aoki2.si.gunma-u.ac.jp/R/misc2.html)

# A群 n=20 で，M=50.00, SD=10.00 の正規乱数を作る
gendat1 <- function(n, mean, sd) return(scale(rnorm(n)) *
sd + mean)
A <- gendat1(20, 50, 10)
cat("mean =", mean(A), "   S.D.=", sd(A), "\n")  #平均と標準偏差の確認
``````
``````## mean = 50    S.D.= 10
``````
``````# B群 n=20 で，M=50.00, SD=10.00 の正規乱数を作る
B <- gendat1(20, 52, 10)
cat("mean =", mean(B), "   S.D.=", sd(A), "\n")  #平均と標準偏差の確認
``````
``````## mean = 52    S.D.= 10
``````
``````t.test(A, B)  # t 検定で確認
``````
``````##
##  Welch Two Sample t-test
##
## data:  A and B
## t = -0.6325, df = 38, p-value = 0.5309
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -8.402  4.402
## sample estimates:
## mean of x mean of y
##        50        52
##
``````
``````# X群 n=200 で，M=50.00, SD=10.00 の正規乱数を作る
gendat1 <- function(n, mean, sd) return(scale(rnorm(n)) *
sd + mean)
X <- gendat1(200, 50, 10)
cat("mean =", mean(X), "   S.D.=", sd(X), "\n")  #平均と標準偏差の確認
``````
``````## mean = 50    S.D.= 10
``````
``````# Y群 n=200 で，M=50.00, SD=10.00 の正規乱数を作る
Y <- gendat1(200, 52, 10)
cat("mean =", mean(Y), "   S.D.=", sd(Y), "\n")  #平均と標準偏差の確認
``````
``````## mean = 52    S.D.= 10
``````
``````t.test(X, Y)  # t検定で確認
``````
``````##
##  Welch Two Sample t-test
##
## data:  X and Y
## t = -2, df = 398, p-value = 0.04618
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.96594 -0.03406
## sample estimates:
## mean of x mean of y
##        50        52
##
``````

箱ひげ図で確認

``````# 2つの箱ひげ図を並べる
par(family = "HiraKakuProN-W3")  # 日本語表示のため
par(mfrow = c(1, 2))
bp <- boxplot(A[, 1], B[, 1], ylim = c(0, 100),
names = c("A", "B"), main = "N = 40", ylab = "点数",
xlab = "Group", message = FALSE)
text(1, 5, paste("n = 20"))
text(2, 5, paste("n = 20"))
boxplot(X[, 1], Y[, 1], ylim = c(0, 100), names = c("X",
"Y"), main = "N = 400", ylab = "点数", xlab = "Group")
text(1, 5, paste("n = 200"))
text(2, 5, paste("n = 200"))
``````

Group X と Y (N = 400) の分布の重なり具合をヒストグラムで確認

``````library(ggplot2)  # ggplot2 パッケージを使用
df <- data.frame(cond = factor(rep(c("X", "Y"),
each = 200)), score = c(X, Y))
# 各群の平均値を割り当て
library(plyr)  # plyr パッケージを使用
cdf <- ddply(df, .(cond), summarise, rating.mean = mean(score))
# ヒストグラムの重ね合わせ（平均に線を引く）
ggplot(df, aes(x = score, fill = cond)) + geom_histogram(binwidth = 0.5,
alpha = 0.5, position = "identity") + geom_vline(data = cdf,
aes(xintercept = rating.mean, colour = cond), linetype = "dashed",
size = 1)
``````

効果量の計算

``````library(compute.es)  # 効果量を計算するパッケージ
mes(mean(A[, 1]), mean(B[, 1]), sd(A[, 1]), sd(B[,
1]), length(A[, 1]), length(B[, 1]))
``````
``````## \$MeanDifference
##        d    var.d        g    var.g
## -0.20000  0.10050 -0.19603  0.09655
##
## \$Correlation
##        r    var.r
## -0.09950  0.02439
##
## \$Log_Odds
##     log_odds var.log_odds
##      -0.3628       0.3306
##
## \$Fishers_z
##        z    var.z
## -0.09983  0.02703
##
## \$TotalSample
##  n
## 40
##
``````
``````mes(mean(X[, 1]), mean(Y[, 1]), sd(X[, 1]), sd(Y[,
1]), length(X[, 1]), length(Y[, 1]))
``````
``````## \$MeanDifference
##        d    var.d        g    var.g
## -0.20000  0.01005 -0.19962  0.01001
##
## \$Correlation
##         r     var.r
## -0.099504  0.002439
##
## \$Log_Odds
##     log_odds var.log_odds
##     -0.36276      0.03306
##
## \$Fishers_z
##         z     var.z
## -0.099834  0.002519
##
## \$TotalSample
##   n
## 400
##
``````

3. Visualization（図示方法）

“A picture is worth a thousand p values.” (Loftus, 1993)

棒グラフ（Barplot）

``````# データの読み込み
dat\$クラスサイズ <- as.factor(dat\$クラスサイズ)
dat\$教室 <- as.factor(dat\$教室)

# 2要因の組み合わせの群ごとの平均, SD, n, SE, 95%CIを計算
mean.group <- tapply(dat\$点数, list(dat\$教室, dat\$クラスサイズ), mean)
sd.group <- tapply(dat\$点数, list(dat\$教室, dat\$クラスサイズ), sd)
n.group <- tapply(dat\$点数, list(dat\$教室, dat\$クラスサイズ), length)
stderr <- function(x) sqrt(var(x)/length(x)) # 標準誤差を求める関数
stderr.group <- tapply(dat\$点数, list(dat\$教室, dat\$クラスサイズ), stderr)ci <- function(x) qt(0.975,length(x)-1)*sd(x)/sqrt(length(x)) # 95%CIを求める関数
ci.group <- tapply(dat\$点数, list(dat\$教室, dat\$クラスサイズ), ci)
``````

``````par(family = "HiraKakuProN-W3")  # 日本語表示のため
barplot(mean.group, beside = T, col = c("gray",
"white"), las = 1, legend = c("普通教室", "CALL教室"),
ylim = c(0, 100), xlab = "クラスサイズ", ylab = "点数",
names = c("10人クラス", "20人クラス", "40人クラス"),
main = "Mean", yaxp = c(0, 100, 10))
``````

エラーバーつき棒グラフ（SDの場合）

``````par(family = "HiraKakuProN-W3")  # 日本語表示のため
mids <- barplot(mean.group, beside = T, col = c("gray",
"white"), las = 1, legend = c("普通教室", "CALL教室"),
ylim = c(0, 100), xlab = "クラスサイズ", ylab = "点数",
names = c("10人クラス", "20人クラス", "40人クラス"),
main = "Mean and SD", yaxp = c(0, 100, 10))
arrows(mids, mean.group - sd.group, mids, mean.group +
sd.group, code = 3, angle = 90, length = 0.1)
``````

エラーバーつき棒グラフ（SEの場合）

``````par(family = "HiraKakuProN-W3")  # 日本語表示のため
mids <- barplot(mean.group, beside = T, col = c("gray",
"white"), las = 1, legend = c("普通教室", "CALL教室"),
ylim = c(0, 100), xlab = "クラスサイズ", ylab = "点数",
names = c("10人クラス", "20人クラス", "40人クラス"),
main = "Mean and SE", yaxp = c(0, 100, 10))
arrows(mids, mean.group - stderr.group, mids,
mean.group + stderr.group, code = 3, angle = 90, length = 0.1)
``````

エラーバーつき棒グラフ（95%信頼区間の場合）

``````par(family = "HiraKakuProN-W3")  # 日本語表示のため
mids <- barplot(mean.group, beside = T, col = c("gray",
"white"), las = 1, legend = c("普通教室", "CALL教室"),
ylim = c(0, 100), xlab = "クラスサイズ", ylab = "点数",
names = c("10人クラス", "20人クラス", "40人クラス"),
main = "Mean and 95% CI", yaxp = c(0, 100, 10))
arrows(mids, mean.group - ci.group, mids, mean.group +
ci.group, code = 3, angle = 90, length = 0.1)
``````

``````par(family = "HiraKakuProN-W3")  # 日本語表示のため
mids <- barplot(mean.group, beside = T, col = c("gray",
"white"), las = 1, legend = c("普通教室", "CALL教室"),
ylim = c(0, 100), xlab = "クラスサイズ", ylab = "点数",
names = c("10人クラス", "20人クラス", "40人クラス"),
main = "Mean and 95% CI", yaxp = c(0, 100, 10))
arrows(mids, mean.group - ci.group, mids, mean.group +
ci.group, code = 3, angle = 90, length = 0.1)
text(mids, 3, paste("n =", n.group))
``````

``````par(family = "HiraKakuProN-W3")  # 日本語表示のため
mids <- barplot(mean.group, beside = T, col = c("gray",
"white"), las = 1, legend = c("普通教室", "CALL教室"),
ylim = c(0, 100), xlab = "クラスサイズ", ylab = "点数",
names = c("10人クラス", "20人クラス", "40人クラス"),
main = "Mean and 95% CI", yaxp = c(0, 100, 10))
arrows(mids, mean.group - ci.group, mids, mean.group +
ci.group, code = 3, angle = 90, length = 0.1)
text(mids, 3, paste("n =", n.group))

# 線を加える
segments(x0 = 1.5, y0 = 80, x1 = 1.5, y1 = 83)  #xとyの数値で始点と終点を指定
segments(x0 = 1.5, y0 = 83, x1 = 2.5, y1 = 83)
segments(x0 = 2.5, y0 = 83, x1 = 2.5, y1 = 70)

segments(x0 = 4.5, y0 = 70, x1 = 4.5, y1 = 73)  #xとyの数値で始点と終点を指定
segments(x0 = 4.5, y0 = 73, x1 = 5.5, y1 = 73)
segments(x0 = 5.5, y0 = 70, x1 = 5.5, y1 = 73)

segments(x0 = 7.5, y0 = 51, x1 = 7.5, y1 = 68)  #xとyの数値で始点と終点を指定
segments(x0 = 7.5, y0 = 68, x1 = 8.5, y1 = 68)
segments(x0 = 8.5, y0 = 68, x1 = 8.5, y1 = 64)

# アステリスクを加える
text(x = 2, y = 85, labels = "***", family = "mono",
font = 1, ps = 8)
text(x = 8, y = 70, labels = "***", family = "mono",
font = 1, ps = 8)
``````

箱ひげ図（Boxplot）

そのため最近では，棒グラフよりも箱ひげ図のほうがデータの提示に適しているとされている。

Larson-hall, J., & Herrington, R. (2010). Improving data analysis in second language acquisition by utilizing modern developments in applied statistics. Applied Linguistics, 31, 368–390. doi:10.1093/applin/amp038

``````par(family = "HiraKakuProN-W3")  # 日本語表示のため
dat\$クラスサイズ <- as.factor(dat\$クラスサイズ)
dat\$教室 <- as.factor(dat\$教室)
dat2 <- split(dat, dat\$クラスサイズ)
par(mfrow=c(1,3)) # 縦に1つ，横に3つのグラフを並べる
boxplot(dat2[[1]]\$点数~dat2[[1]]\$教室, names=c("普通教室", "CALL教室"), ylim=c(0,100), main="10人クラス", ylab="点数")
boxplot(dat2[[2]]\$点数~dat2[[2]]\$教室, names=c("普通教室", "CALL教室"), ylim=c(0,100), main="20人クラス", ylab="点数")
boxplot(dat2[[3]]\$点数~dat2[[3]]\$教室, names=c("普通教室", "CALL教室"), ylim=c(0,100), main="40人クラス", ylab="点数")
``````

``````# 箱ひげ図に個別データをプロットする
Normal <- c(0.2, 0.29, 0.051, 0.38, 0.067, 0.04,
0.111, 0.202, 0.109, 0.124, -0.03, 0.033)
Scanning <- c(0.043, 0.374, 0.637, 0.579, 0.477,
0.098, 0.394, 0.376, -0.018, 0.648, 0.118, 0.175)
score <- c(Normal, Scanning)
condition <- factor(c(rep("Normal", 12), rep("Scanning",
12)))
boxplot(score ~ condition, xlab = "Condition",
ylab = "Oxy-Hb (mM-mm)", ylim = c(-0.1, 0.7), yaxp = c(-0.1,
0.7, 8))
stripchart(score ~ condition, pch = 16, vert = TRUE,
``````

``````Normal <- c(0.2, 0.29, 0.051, 0.38, 0.067, 0.04,
0.111, 0.202, 0.109, 0.124, -0.03, 0.033)
Scanning <- c(0.043, 0.374, 0.637, 0.579, 0.477,
0.098, 0.394, 0.376, -0.018, 0.648, 0.118, 0.175)
score <- c(Normal, Scanning)
condition <- factor(c(rep("Normal", 12), rep("Scanning",
12)))
boxplot(score ~ condition, xlab = "Condition",
ylab = "Oxy-Hb (mM-mm)", ylim = c(-0.1, 0.7), yaxp = c(-0.1,
0.7, 8))
stripchart(score ~ condition, pch = 16, vert = TRUE,
x <- c(1, 2)
a <- mean(Normal)
b <- mean(Scanning)
means <- c(a, b)
c <- sd(Normal)
d <- sd(Scanning)
sds <- c(c, d)
points(x + 0.2, means, pch = 18, col = "red",
cex = 2)
arrows(x + 0.2, means, x + 0.2, means + sds, length = 0.1,
angle = 45, col = "red")
arrows(x + 0.2, means, x + 0.2, means - sds, length = 0.1,
angle = 45, col = "red")
``````

個人の変化を示す図

``````Normal <- c(0.2, 0.29, 0.051, 0.38, 0.067, 0.04,
0.111, 0.202, 0.109, 0.124, -0.03, 0.033)
Scanning <- c(0.043, 0.374, 0.637, 0.579, 0.477,
0.098, 0.394, 0.376, -0.018, 0.648, 0.118, 0.175)
Value <- c(Normal, Scanning)
Condition <- factor(c(rep("Normal", 12), rep("Scanning",
12)))
dat <- data.frame(Value, Condition)
dat\$Subject <- factor(c(rep(1:12), rep(1:12)))  # id番号付与

library(lattice)  # lattice パッケージを使用
(each <- xyplot(Value ~ Condition, group = Subject,
type = c("l"), data = dat, ylab = "Oxy-Hb (mM-mm)", ylim = c(-0.1,
0.7), scales = list(y = list(at = c(-0.1, 0, 0.1,
0.2, 0.3, 0.4, 0.5, 0.6, 0.7)))))
``````

（参考）http://ito-hi.blog.so-net.ne.jp/2009-03-02-1

``````library(lattice)  # lattice パッケージを使用
each <- xyplot(Value ~ Condition, group = Subject,
type = c("l"), data = dat, ylab = "Oxy-Hb (mM-mm)", ylim = c(-0.1,
0.7), scales = list(y = list(at = c(-0.1, 0, 0.1,
0.2, 0.3, 0.4, 0.5, 0.6, 0.7))))
library(latticeExtra)  # lattice パッケージの図を重ね合わせるパッケージ
normal <- subset(dat, Condition == "Normal")
a <- mean(normal\$Value)
scanning <- subset(dat, Condition == "Scanning")
b <- mean(scanning\$Value)
value <- c(a, b)
condition <- factor(c(rep("Normal", 1), rep("Scanning",
1)))
dat2 <- data.frame(value, condition)
all <- xyplot(value ~ condition, col = "black",
lwd = 5, type = c("l"), data = dat, ylab = "Oxy-Hb (mM-mm)",
ylim = c(-0.1, 0.7), scales = list(y = list(at = c(-0.1,
0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7))))

each + as.layer(all, axes = NULL)
``````

事前・事後データ比較における散布図の活用

``````par(family = "HiraKakuProN-W3")  # 日本語表示のため
# 青木繁伸先生（群馬大学）のコードを参照
source("http://aoki2.si.gunma-u.ac.jp/R/src/mix2.R",
encoding = "euc-jp")
# 前田（2008, p. 257）の人数，平均，標準偏差を再現
mu <- matrix(c(538, 565.1, 421.2, 424), byrow = TRUE,
ncol = 2)
sigma <- matrix(c(35.2, 71.4, 95.5, 111.2), byrow = TRUE,
ncol = 2)
r <- c(0.8, 0.8)  # 事前・事後の相関係数がわからないため0.8としている
n <- 393  # 合計人数
d <- mix2(n, mu, sigma, r, prob = c(1.1, 8.9))  # 人数の割合
``````
``````## 538 565.1 35.2 71.4 0.8
## 421.2 424 95.5 111.2 0.8
``````
``````dat <- as.data.frame(d)  # データフレームの形に変更
a <- subset(dat, which == "1")  # グループごとに分割
b <- subset(dat, which == "2")  # グループごとに分割

# データを描画
par(mar = c(5, 5, 2, 2))  #余白の指定．下，左，上，右の順
par(xaxs = "i", yaxs = "i")
plot(a\$d.1, a\$d.2, pch = 1, cex = 1.3, xlim = c(0,
1000), ylim = c(0, 1000), las = 1, xlab = "事前", ylab = "事後",
cex.lab = 1.2, xaxp = c(0, 1000, 10), yaxp = c(0, 1000,
10))
points(b\$d.1, b\$d.2, pch = 20, cex = 0.25)  #2グループ目のデータを図に入れる
lines(par()\$usr[1:2], par()\$usr[3:4], lty = 2)  #対角線をダッシュで引く
lines(c(mean(a\$d.1), mean(a\$d.1)), c(0, mean(a\$d.2)),
lty = 2)  #x1, x2, y1, y2 の順
lines(c(0, mean(a\$d.1)), c(mean(a\$d.2), mean(a\$d.1)),
lty = 2)  #x1, x2, y1, y2 の順
lines(c(mean(b\$d.1), mean(b\$d.1)), c(0, mean(b\$d.2)),
lty = 1)  #x1, x2, y1, y2 の順
lines(c(0, mean(b\$d.1)), c(mean(b\$d.2), mean(b\$d.1)),
lty = 1)  #x1, x2, y1, y2 の順
``````

シミュレーションでデータを生成しているので，結果は毎回異なることに注意する。