R を使った Resampling と作図

2012年12月01日（土）外国語教育メディア学会（LET）中部支部第80回支部研究大会のシンポジウム

『統計手法を用いたデータ分析とその解釈—何が必要でどう利用すべきか—』での

発表「より好ましい統計解析と図示方法」で使用したコードや関連した内容です。

使用したスライドやリンク集はこちらにあります。

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) {
    op <- par(no.readonly = TRUE)
    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)

それぞれの検定で p 値の比較（Permutation/Randomization Test）

スライド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 <- read.csv("http://mizumot.com/handbook/wp-content/uploads/ch07-1.csv",header=T, fileEncoding="CP932")
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)

各群の n を入れる場合

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 <- read.csv("http://mizumot.com/handbook/wp-content/uploads/ch07-1.csv",header=T, fileEncoding="CP932")
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, 
    add = TRUE)

箱ひげ図に平均値と標準偏差を加える場合 （参照）http://www.youtube.com/watch?v=xzY3GnyLAuM

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, 
    add = 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)

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

前田啓朗（2008）「WBTを援用した授業で成功した学習者・成功しなかった学習者」で 使用されている，群ごとのデータを散布図で比較する方法をシミュレーションデータで再現。

特定の平均値，標準偏差，相関係数を持ったデータの生成 （参照）http://aoki2.si.gunma-u.ac.jp/R/mix2.html

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 の順

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