R: MICをためしてみる (2)

前回のつづき。有意差が気になる向きには、ということでランダマイゼーションテストをやってみた。

こういうコード。

library(minerva)

## bootstrap function
mic.fun <- function(dat, ind) {
  mine(dat$x, dat[sample.int(nrow(dat)), ]$y)$MIC
}

set.seed(12345)
n <- 500
data <- vector("list", 3)
mic <- vector("list", 3)

data[[1]] <- data.frame(x = runif(n, -1, 1))
data[[1]]$y <- data[[1]]$x + rnorm(n, 0, 0.1)
data[[2]] <- data.frame(x = runif(n, -1, 1))
data[[2]]$y <- sqrt(1 - data[[2]]$x^2) + rnorm(n, 0, 0.1)
data[[3]] <- data.frame(x = runif(n, -1, 1))
data[[3]]$y <- rnorm(n, 0, 0.1)

## plot data
png(width = 520, height = 520)
lapply(data, function(dat) plot(dat$x, dat$y, xlab = "x", ylab = "y"))
dev.off()

## Pearson's correlation
lapply(data, function(dat) cor.test(dat$x, dat$y, method = "pearson"))

## randomization
mic[[1]] <- replicate(2000, mine(data[[1]]$x, data[[1]][sample.int(n), ]$y)$MIC)
mine(data[[1]]$x, data[[1]]$y)
quantile(mic[[1]], probs = c(0.025, 0.975))

mic[[2]] <- replicate(2000, mine(data[[2]]$x, data[[2]][sample.int(n), ]$y)$MIC)
mine(data[[2]]$x, data[[2]]$y)
quantile(mic[[2]], probs = c(0.025, 0.975))

mic[[3]] <- replicate(2000, mine(data[[3]]$x, data[[3]][sample.int(n), ]$y)$MIC)
mine(data[[3]]$x, data[[3]]$y)
quantile(mic[[3]], probs = c(0.025, 0.975))

結果

data[[1]]

data[[2]]

data[[3]]

> ## Pearson's correlation
> lapply(data, function(dat) cor.test(dat$x, dat$y, method = "pearson"))
[[1]]

	Pearson's product-moment correlation

data:  dat$x and dat$y 
t = 127.1799, df = 498, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0 
95 percent confidence interval:
 0.9820854 0.9873632 
sample estimates:
      cor 
0.9849522 


[[2]]

	Pearson's product-moment correlation

data:  dat$x and dat$y 
t = 0.4846, df = 498, p-value = 0.6281
alternative hypothesis: true correlation is not equal to 0 
95 percent confidence interval:
 -0.06610423  0.10919490 
sample estimates:
       cor 
0.02171222 


[[3]]

	Pearson's product-moment correlation

data:  dat$x and dat$y 
t = -0.208, df = 498, p-value = 0.8353
alternative hypothesis: true correlation is not equal to 0 
95 percent confidence interval:
 -0.09692998  0.07843608 
sample estimates:
         cor 
-0.009318603 


> 
> ## randomization
> mic[[1]] <- replicate(2000, mine(data[[1]]$x, data[[1]][sample.int(n), ]$y)$MIC)
> mine(data[[1]]$x, data[[1]]$y)
$MIC
[1] 0.9623742

$MAS
[1] 0.05576998

$MEV
[1] 0.9623742

$MCN
[1] 5.321928

$`MIC-R2`
[1] -0.0077567

> quantile(mic[[1]], probs = c(0.025, 0.975))
     2.5%     97.5% 
0.1418036 0.1891671 
> 
> mic[[2]] <- replicate(2000, mine(data[[2]]$x, data[[2]][sample.int(n), ]$y)$MIC)
> mine(data[[2]]$x, data[[2]]$y)
$MIC
[1] 0.6487575

$MAS
[1] 0.4096336

$MEV
[1] 0.6487575

$MCN
[1] 5.169925

$`MIC-R2`
[1] 0.6482861

> quantile(mic[[2]], probs = c(0.025, 0.975))
     2.5%     97.5% 
0.1414939 0.1872861 
> 
> mic[[3]] <- replicate(2000, mine(data[[3]]$x, data[[3]][sample.int(n), ]$y)$MIC)
> mine(data[[3]]$x, data[[3]]$y)
$MIC
[1] 0.159976

$MAS
[1] 0.01830776

$MEV
[1] 0.159976

$MCN
[1] 5.321928

$`MIC-R2`
[1] 0.1598891

> quantile(mic[[3]], probs = c(0.025, 0.975))
     2.5%     97.5% 
0.1417466 0.1869522 

それらしい結果だが、これでよいのだろうか。

タグ：R