多変量での検定

こんな感じのデータで検定できるかどうか聞かれた。

データの意味としては、O群とC群で画像のズレがx軸とy軸でとれるらしく、2群で差がありそうなことを言いたいらしい。
　
多変量に拡張された平均値の差の検定としてWilks検定があるらしい。RではMANOVAでできるらしい。

library(robust)
library(GGally)

## Create
dat <- read.table(header = TRUE, text = "
group   x       y
o      -2       1
o       0       4
o       1       1
o       2       2
o       3      -5
o       3       2
o       3       4
o       4       2
o       5       3
o       6       2
c      -5      -2
c      -4       0
c      -4       0
c      -1      -2
c      -1       0
c      -1      -2
c      -1       0
c       0      -1
c       0       0
c       0       2
")

plot(dat$x, dat$y, col=dat$group, pch=as.numeric(dat$group))
abline(h=0, v=0, lty=2)
legend("topleft", legend=paste("Group", c("O", "C")), bty="n", col=2:1, pch=2:1, cex=1.5)
# 外れ値の影響が大きそうなので、ロバスト回帰にしている
lm0 <- lmrob(y ~ 0 + x, data=dat)
abline(0, lm0$coefficients)

で、なんとなく2群とも、上の図で描いたような直線に乗っているっぽくて、適当に回転してしまってその直線上での新しいx座標でも考えて無難にt検定でもしたら?というのが案その1。

# 線形回帰の傾きから回転角を出す。
phi <- atan(lm0$coefficients)
# 回転行列
rot <- matrix(c(cos(-phi), sin(-phi), -sin(-phi), cos(-phi)), nr=2)
# 回転した座標
dat1 <- dat
for(i in seq(nrow(dat1))){
	dat1[i, 2:3] <- rot %*% t(dat[i, 2:3])
}

plot(dat1$x, dat1$y, col=dat$group, pch=as.numeric(dat$group))
abline(h=0, v=0, lty=2)
legend("topleft", legend=paste("Group", c("O", "C")), bty="n", col=2:1, pch=2:1, cex=1.5)

t.test(dat1$x ~ dat$group)

	Welch Two Sample t-test

data:  dat1$x by dat$group
t = -4.7934, df = 17.175, p-value = 0.0001646
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -6.736221 -2.620910
sample estimates:
mean in group c mean in group o 
      -1.764505        2.914060 

　
見た感じ、O群とC群の分布は重なっていなさそうだから、2群の重心というか中心点的なものを求めて、2点の距離をとしておくと、リサンプリングを使って、2群をひとまとめにした集合から10点ずつ偽サンプルを作成して、同じように2群の中心点の距離を求め、になる確率が求めるp値となる、というのが案その2。

ggpairs(dat,
	upper = list(continuous = "density", combo = "box"),
	lower = list(continuous = "smooth", combo = "dot"),
	color = "group"
)

# 中心点
G <- apply(dat[, 2:3], 2, tapply, dat$group, mean)
# 2点間の距離
R <- sqrt(sum(apply(G, 2, diff)^2))

niter <- 10000
res <- rep(0, niter)
for(i in seq(niter)){
	idx <- sample(seq(nrow(dat)), size=nrow(dat), replace=TRUE)
	dat_new <- dat[idx, 2:3]
	G_sim <- apply(dat_new, 2, tapply, dat$group, mean)
	R_sim <- sqrt(sum(apply(G_sim, 2, diff)^2))
	res[i] <- R_sim
}

hist(res, nclass=50)
abline(v=R, lty=2)
mean(res > R) # p value

[1] 0.0015