5 おまけ

たしかに、ex-Gaussian分布を用いたフィッティングは、平均値や標準偏差の算出と比べると、 必要とされる計算の手順が多く複雑です。 また先述のとおり、この手法自体はそれほど新しいものでもないのですが、 研究者のプログラミング能力の不足やデータの歪曲に対する問題意識の欠如のため、 それほど一般に浸透しているわけではありません。 そのため、たとえば学会発表や論文などで唐突に ex-Gaussian分布を利用した反応時間解析をもちだすと、 聴衆や読者に不必要に怪しまれてしまうことも充分考えられます。 「平均値では結果が出ないから、ややこしい解析もちだしてきやがったな」 「平均値のような基本的な指標でちゃんと出てないってことは、 信頼できる結果ではなさそうだ」と。 よってそのような反応を受けるくらいなら、 反応時間分布の歪曲の問題は認識しながらも、 戦略的に平均と標準偏差のみを用いて話をするという態度もあり得ると思います。

しかしそれでは、平均値や標準偏差、 あるいはノンパラメトリクスなら中央値や四分位範囲といった指標は、 理論分布のフィッティングによる推定値よりも 「基本的」で「信頼できる」指標なのでしょうか。 たとえば平均値は、すべての観測値を足し合わせ、 それをサンプルサイズで割った値です。 しかし、通常の平均では単純に全サンプルを足し合わせますが、 それぞれのサンプルに別々の係数をかけたのちに足し合わせる、 重み付き平均というものも存在します。

こう考えると、われわれがいつも使っている平均値というものは、 すべての観測値の重みを1にした、 重み付き平均の非常に特殊な場合ということができます。

中央値についてはどうでしょうか。 中央値とは、データ全体を二分する （上側50%と下側50%に分ける）境となる値のことです。 しかしデータを切り分ける方法は、当然ながら二等分だけではありません。 データを四分割する境の点は四分位点と呼ばれ、 第1四分位点はデータを下側25%と上側75%に分ける点、 第3四分位点は逆に下側75%と上側25%に分ける点です。 百等分する点を考えれば百分位点となります。 百分位点という視点からみれば、 中央値は第50百分位点（下から数えて50番目の百分位点）と呼ぶことができるでしょう。 ではなぜ、第50百分位点はデータ要約において重要な代表値として特別扱いされ、 第49百分位点や第51百分位点は使えないのでしょうか。

もちろん多くのひとは、これらのいちゃもんに対して、 「重み付け平均や第49・51百分位点よりも、普通の平均値や中央値のほうが、 データを平等に利用していて恣意的でない」という反論をするでしょう。 しかし、その「平等」とか「恣意的」という理由は、 あくまで人間の意味的判断によるものであって、 「俺は納得できる」という主張以上のものではないのです。 不偏性や頑健性などの数理的特徴を用いて議論することもできるでしょうが、 すべての統計ユーザが、 不偏性という数理的性質のために平均値を選択し、 中心極限定理のためにパラメトリクスを選択しているとは到底思えません。 多くのひとは、単に「自分が自然に納得できる」というだけの理由をもって、 平均値や中央値、 標準偏差や四分位範囲などの指標を選択しているように思えます。 もしそうだとすれば、データの要約にどんな指標を用いるのがよいかは、 その人物の理解力や数学的センスに制約されることになります。

結局わたしがいいたいのは、 皆さんが信じてやまない平均値などの「一般的」な指標も、 じつはそれほど確固たる理由をもって使われているわけではない、ということです。 これらの指標が万人に受け入れられている理由は、 万人がこの指標の定義になんとなく納得できるから、 ぐらいのものなのではないでしょうか。 意地悪な言いかたをすれば、 堪え性も数学的教養もない一般聴衆の「最大公約数」を想定したとき、 そのセンスでついてこられる計算のレベルが平均値であり、 標準偏差であるのでしょう。 そして筆者は、これがひとつの真理ではないかと思っています。 データ解析は、数値要約にしろフィッティングにしろ、 生データに何らかの処理を加えて人間に解釈しやすくするという意味ではどれも同じです。 そして解析の結果は、 「処理の結果どれだけ解釈しやすくなったか」と 「処理は妥当なものだったか」というふたつの天秤によってはかられます。 さらに究極的には、このふたつの論点は、 いずれも「結果をみる側にとって、それをどれだけ認められるか」 という問題に集約されます。 解析に使った手法は認めてもらえても、 それによりでてきた結果が学術的にまったく面白くなければ、 意味はありません。 逆にどんなに興味深い結果を示そうとも、 解析に使った手法があまりにも複雑で理解に苦しむものでは、 やはり誰も振り向いてくれないでしょう。 解析結果そのもののもつ意味と、 それを抽出するために使った方法の理解しやすさというふたつのバランスのなかで、 聞き手に理解され受け入れられる折衷点を探さなければならないのです。 戦略的にex-Gaussian分布解析を用いないということも考えられる、と最初に述べたのは、 このような思想が筆者のあたまのなかにあったからでした。

実験によって得られたデータは、 いつだって「そのデータ」たったひとつです。 しかしそれに対して可能な解析は無数に存在し、 「必ずこうすればよい」というような万能解は存在しません。 筆者個人としては、まとめの項でも書いたとおり、 そこに存在する差を見逃さないということがもっとも重要だと考えており、 そのためにex-Gaussian分布解析は有用なツールだと思います。 しかしときと場合によっては、あえてそれを用いないほうがよいことだってあるでしょう。 導き出される情報とその抽出のための手法の分かりやすさというふたつの観点のあいだで、 あなたのデータのもつ意味をもっとも効率的に伝えられる解析方法を選び出してください。

rexGauss <- function(n, mu=300, sigma=50, tau=50) {
  return(rnorm(n, mu, sigma) + rexp(n, 1/tau))
}

dexGauss <- function(x, mu=300, sigma=50, tau=50, log=FALSE, limit=10^-4) {
  if (tau < limit) {
    d <- dnorm(x, mean=mu, sd=sigma, log=TRUE)
  } else {
    d <- log(1 / tau) - (x - mu) / tau + (sigma^2) / (2 * (tau^2)) +
      pnorm((x - mu) / sigma - sigma / tau, 0, 1, log.p=TRUE)
  }
  if (!log) {
    d <- exp(d)
  }
  return(d)
}

fitexGauss <- function(
  rt, # A vector of RTs to analyze distribution
  type = c("ml", "chisq", "ks"), # Fitting procedure
  probs = seq(0.1, 0.9, by=0.2), # Quantiles for chisq procedure
  warn = FALSE, # Whether to print warnings in optim
  maxit = 1000, # Maximal number of iteration in optim
  qmax = 100, # Maximal number of quantiles in Q-Q plot
  ... # Additional arguments for optim
) {
  type <- type[1]
  rt <- rt[!is.na(rt)]
  nrt <- length(rt)
  arg <- list(...)
  arg$control$maxit <- maxit

  # 1st - 3rd cumulants
  k <- rep(NA, 3)
  k[1] <- mean(rt)
  k[2] <- var(rt)
  k[3] <- mean((rt - k[1])^3)

  # Estimates of ex-Gaussian parameters by moment statistics
  momfit <- c(mu=NA, sigma=NA, tau=NA)
  momfit[3] <- ifelse(k[3] > 0, (k[3]/2)^(1/3), 10^(-5))
  momfit[1] <- k[1] - momfit[3]
  v <- k[2] - momfit[3]^2
  momfit[2] <- ifelse(v > 0, sqrt(v), 10^(-5))

  # Definition of function to minimize in optim
  if (type == "ml") {
    # Sign-inverted sum of prob density (in log scale)
    fn <- function(THETA) {
      THETA <- exp(THETA)
      return(-sum(dexGauss(rt, THETA[1], THETA[2], THETA[3], TRUE)))
    }
  } else if (type == "chisq") {
    # Chi-squared statistics in contingency table
    n_lev <- 1:(length(probs) + 1)
    n_pred <- diff(c(0, probs, 1)) * nrt
    fn <- function(THETA) {
      THETA <- exp(THETA)
      q <- c(-Inf, qexGauss(probs, THETA[1], THETA[2], THETA[3]), Inf)
      n_obs <- table(cut(rt, breaks=q, labels=n_lev))
      return(sum(((n_obs - n_pred) ^ 2) / n_pred))
    }
  } else if (type == "ks") {
    # Y-axis discrepancy in cumulative prob dist
    rt_sorted <- sort(rt)
    p <- seq(1 / nrt, by=1 / nrt, length.out=nrt)
    fn <- function(THETA) {
      THETA <- exp(THETA)
      return(max(abs(pexGauss(rt_sorted, THETA[1], THETA[2], THETA[3]) - p)))
    }
  }
  optold <- options(warn=ifelse(warn, 0, -1))
  init <- momfit
  init <- log(init)
  loglikefit <- optim(par=init, fn=fn, arg)
  param <- loglikefit$par
  param <- exp(param)
  options(optold)

  if (FALSE) { # To check the failure in moment estimates
    if (any(momfit == 10^(-5))) {
      cat("fitexGauss: Moment estimate for sigma or tau became 0\n")
      print(matrix(c(momfit, param), byrow=TRUE, nrow=2,
        dimnames=list(c("mom","loglike"), c("mu","sigma","tau"))))
    }
  }

  rtnvalue <- list(mom=momfit, loglike=param,
    optimrslt=loglikefit[-1], convergence=loglikefit$convergence, rt=rt)
  class(rtnvalue) <- "htest.exGauss"

  # Calculate r^2 in Q-Q plot
  if (nrt > (qmax+1)) {
    # qmax-tiles should be separated by qmax+1 points
    # Require qmax+3 seps to eliminate lower and upper limits
    #   which are -Inf and Inf in the theoretical CDF
    obs <- quantile(rt, probs=seq(0, 1, length.out=qmax+3))
  } else {
    obs <- sort(rt)
  }
  l <- length(obs)
  pred <- qexGauss(seq(0, 1, length.out=l),
    param["mu"], param["sigma"], param["tau"])
  obs <- obs[-c(1, l)]
  pred <- pred[-c(1, l)]

  rtnvalue$lmrslt <- lm(obs~pred)
  s <- summary(rtnvalue$lmrslt)
  rtnvalue$r.squared <- s$r.squared

  # Kolmogorov-Smirnov Test
  rtnvalue$ksrslt <- suppressWarnings(ks.test(rt, pexGauss,
    mu=param["mu"], sigma=param["sigma"], tau=param["tau"]))

  return(rtnvalue)
}

pexGauss <- function(q, mu=300, sigma=50, tau=50,
  lower.tail=TRUE, log.p=FALSE) {
  mayinf <- -(q-mu)/tau
  if (any(is.infinite(exp(mayinf)))) {
    # Practically equal to Gaussian distribution
    #   when tau is very small
    p <- pnorm((q-mu)/sigma)
  } else {
    p <- pnorm((q-mu)/sigma) - exp(mayinf +
      (sigma^2)/(2*(tau^2)) + log(pnorm((q-mu)/sigma - sigma/tau)))
  }
  if (!lower.tail) {
    p <- 1 - p
  }
  if (log.p) {
    p <- log(p)
  }
  return(p)
}

qexGauss <- function(p, mu=300, sigma=50, tau=50,
  lower.tail=TRUE, log.p=FALSE) {
  if (log.p) {
    p <- exp(p)
  }
  if (!lower.tail) {
    p <- 1 - p
  }
  if (any(c(p < 0, p > 1))) {
    stop("qexGauss: <p> must be from 0 to 1")
  }
  pfunc <- function(x) {
    pexGauss(x, mu=mu, sigma=sigma, tau=tau)
  }
  m <- mu + tau
  bin <- ifelse(p > pfunc(m), 1, -1) * mean(c(sigma, tau))
  q <- mapply(FUN=function(p_tmp, b_tmp) {
    if (p_tmp == 0) {
      return(-Inf)
    } else if (p_tmp == 1) {
      return(Inf)
    }
    cri <- m
    while (sign(pfunc(cri) - p_tmp) == sign(pfunc(cri + b_tmp) - p_tmp)) {
      cri <- cri + b_tmp
    }
    func <- function(r){ pfunc(r) - p_tmp }
    return(uniroot(func, sort(c(cri, cri + b_tmp)))$root)
  }, p, bin)
  return(q)
}

plot.htest.exGauss <- function(
  x,
  y = NA,
  main = "",
  xlab = "Time [ms]",
  ylab = ifelse(freq, "Frequency", "Density"),
  border = "black",
  ylim = NULL,
  freq = TRUE,
  breaks = 30,
  plot = TRUE,
  ...
) {
  rt  <- x$rt
  fit <- x$loglike

  # Draw histogram
  histrslt <- hist(rt, breaks=breaks, plot=plot,
    border=border,
    main=main, xlab=xlab, ylab=ylab,
    ylim=ylim, freq=freq)

  if (x$optimrslt$convergence == 0) {
    # Density function with fitted exGauss parameters
    # xlim <- range(histrslt$breaks)
    xlim <- par("usr")[1:2]
    x <- seq(xlim[1], xlim[2], length.out=100)
    y <- dexGauss(x, fit["mu"], fit["sigma"], fit["tau"])

    if (!freq) { # Density plot
      h <- 1
    } else { # Frequency plot
      h <- max(histrslt$counts) / max(histrslt$density)
    }
    lines(x, y * h, ...)
  }

  invisible()
}

nsub <- 15

source("exGauss.R")
set.seed(100)

word  <- c(rep(rep(c("red", "green"), 2), each=50), rep(c("xxx", "xxxxx"), each=50))
color <- rep(c("red", "green", "green", "red", "red", "green"), each=50)
cond  <- rep(c("cng", "incng", "neut"), each=100)
proto <- data.frame(word=word, color=color, cond=cond)

err <- function(
  eprob # Probability of error
) {
  sample(c(0,1), 100, replace=TRUE, prob=c(eprob,100-eprob))
}

dir.create("./data", showWarnings=FALSE)
options(warn=-1)
makesample <- function(sub) {

  mu    <- c(500, 600, 520) + rnorm(3, sd=10) + rnorm(1, sd=50)
  sigma <- c(50,  110, 50)  + rnorm(3, sd=5)  + rnorm(1, sd=10)
  tau   <- c(110, 115, 90)  + rnorm(3, sd=5)  + rnorm(1, sd=10)

  RT  <- round(as.vector(mapply(FUN=rexGauss, 100, mu, sigma, tau)))
  cor <- as.vector(sapply(FUN=err, X=c(3,15,5)))
  tmp <- proto
  tmp$RT  <- RT
  tmp$cor <- cor
  tmp <- tmp[sample(1:300),]
  tmp$trial <- 1:300
  tmp <- tmp[,c("trial", "word", "color", "cond", "RT", "cor")]

  filename <- paste("./data/sub_", sub, ".txt", sep="")
  con <- file(filename, open="w")
  textout <- c("Pseudo-data of Stroop task for ex-Gaussian RT analysis",
      "Experimenter\tSatou Tanaka",
      paste("Experimental date\t", Sys.Date()-nsub+sub, sep=""))
  writeLines(text=textout, con=con, sep="\n")
  cat("\n", file=con)
  close(con)
  write.table(tmp, file=filename, append=TRUE, sep="\t", quote=FALSE,
    col.names=TRUE, row.names=FALSE)

}
lapply(X=1:nsub, FUN=makesample)
options(warn=0)

nsub <- 15 # Number of subjects

source("exGauss.R")

histRT <- function(
  rt, # List[3] containing RTs for each condition
  leg, # List[3] containing 5 stats & params for each condition
  nbin=30 # Number of bins in the histograms
) {

  ylim <- c(min(unlist(rt)), max(unlist(rt))+1)
  br   <- seq(ylim[1], ylim[2], length.out=nbin)
  leg  <- round(leg)
  prefix <- c("Mean:", "SD:", "Mu:", "Sigma:", "Tau:")

  hist(rt$cng, breaks=br, main="Congruent")
  legend("topright", legend=paste(prefix, leg[1,]), bty="n")

  hist(rt$neut, breaks=br, main="Neutral")
  legend("topright", legend=paste(prefix, leg[2,]), bty="n")

  hist(rt$incng, breaks=br, main="Incongruent")
  legend("topright", legend=paste(prefix, leg[3,]), bty="n")

}
analyzeSub <- function(filename) {

  sub <- matrix(NA, nrow=3, ncol=5) # Container for statistics and estimates
  tmp <- read.table(filename, header=TRUE, skip=4)
  cor <- tmp[as.logical(tmp$cor),] # Correct trials only

  cond <- cor$cond
  RT   <- cor$RT
  RT   <- split(RT, cond)

  sub[,1]    <- sapply(X=RT, FUN=mean)[c("cng","neut","incng")]
  sub[,2]    <- sapply(X=RT, FUN=sd)[c("cng","neut","incng")]
  sub[1,3:5] <- fitexGauss(RT$cng)$loglike
  sub[2,3:5] <- fitexGauss(RT$neut)$loglike
  sub[3,3:5] <- fitexGauss(RT$incng)$loglike

  histRT(RT, sub)

  return(sub)

}
statRT <- function(
  param, # Matrix[ncond, nsub] containing estimates of one parameter
  title # Name of the parameter
) {

  param <- data.frame(RT=as.vector(param), cond=rep(c("cng","neut","incng"), nsub),
      sub=rep(paste("sub", 1:nsub, sep=""), each=3))
  rsltaov <- aov(RT ~ cond + sub, data=param)

  title <- paste("### ", title, " ###", sep="")
  cat(paste(c(rep("#", nchar(title)), "\n"), collapse=""))
  cat(title, "\n", sep="")
  cat(paste(c(rep("#", nchar(title)), "\n\n"), collapse=""))

  cat("### ANOVA ###\n\n")
  print(summary(rsltaov))
  cat("\n")

  cat("### TukeyHSD ###\n\n")
  print(TukeyHSD(rsltaov, "cond"))

}

postscript("histRT.eps", width=3, height=6, paper="special", horizontal=FALSE)
par(mfrow=c(3,1), mar=c(3,4,3,1), las=1)
rslt <- lapply(FUN=analyzeSub, X=paste("./data/sub_", 1:nsub, ".txt", sep=""))
dev.off()
rslt <- array(unlist(rslt), dim=c(3,5,nsub))

mstats <- apply(X=rslt, MARGIN=c(1,2), FUN=mean)
sstats <- apply(X=rslt, MARGIN=c(1,2), FUN=sd)

postscript("paramRT.eps", width=5, height=5, paper="special", horizontal=FALSE)
x <- barplot(mstats, beside=TRUE, names.arg=c("mean","SD","mu","sigma","tau"),
  legend.text=c("congruent","neutral","incongruent"), las=1)
x <- as.vector(x)
mmax <- as.vector(mstats + sstats)
mmin <- as.vector(mstats - sstats)
# arrows(x, mmax, x, mmin, code=3, angle=90, length=0.05) # For error bars
dev.off()

sink("statrslt.txt")
mapply(FUN=statRT, list(rslt[,1,], rslt[,2,], rslt[,3,], rslt[,4,], rslt[,5,]),
  c("Mean", "SD", "Mu", "Sigma", "Tau"))
sink()

vincentize <- function(
  rt, # List[N] containing RTs for N subjects
  n=20, # Number of quantile
  xlim=NA,
  ...
) {

vin <- function(x) {
  return(quantile(x, prob=seq(0, 1, length.out=n+1)))
}
quan <- lapply(X=rt, FUN=vin)
quan <- do.call(rbind, quan)

mquan <- colMeans(quan)
uarea <- 1 / n # Unit area
h <- uarea / diff(mquan)

if (is.na(xlim[1])) {
  xlim <- range(mquan[-length(mquan)])
}
plot(0, 0, type="n", xlim=xlim, ylim=c(0, max(h)),
  xlab="", ylab="Probability Density", ...)
rect(mquan[-length(mquan)], 0, mquan[-1], h)

invisible(list(quan=quan, mquan=mquan))

}