こんな感じのデータで検定できるかどうか聞かれた。
データの意味としては、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