ワインに関する興味ある実験は色々あるが、下記の論文は、困難な実験を周到な計画かつ高度な手法によって成功した例である。自然な考えかたでモデルを作り、最尤法を用いてパラメータを決定する（p.824に手法の簡単な説明がある）。

Oberfeld D, Hecht H, Allendorf U and Wickelmaier F. (2009). Ambient lighting modifies the flavor of wine. Jornal of Sensory Studies, 24 797-832. http://www.staff.uni-mainz.de/oberfeld/downloads/oberfeld_et_al_2009_JSS_wine_color.pdf

この実験の趣旨は、ワインの味が、それを飲む場の照明の色彩によって影響をうけるか、ということにある。完全な成果を出すのに難航し、３度の実験を行っている。なお、第3著者は協力したReingau ワイナリー Fritz Allendorf の跡継ぎであろうか、所属がワイナリーになっている。

実験では、白ワインの色調が解らないように黒のグラスで提供された。ワインはグラスが大変重要であるが、この実験ではZwiesel を用いている。Riedel と並ぶ著名な業者であり、この実験においてグラスの問題はない。但し私は黒のグラスワインというのは見たことがない。ワインは QbAクラスを用いている。部屋の照明を Red, Blue, Green, White の４つに設定し、そこで被験者はグラスのワインを飲み、その印象を記憶する。

第一実験はワイナリーで行われた。意味のある結果は味や香りについては出ず、総合評価に関して有意な差が出た。つまり、「好きである」の１０段階評価（0~9)」と「最高この値段で買ってもよい」において、Red, Blue がWhite, Greenに優越する。残念ながらデータが明記されないのでここで詳細は紹介できないが、t検定の結果では効果量r（論文には記載がない）も .21~.38 あるので認めてよい。

第二と第三の実験は同じく４色の照明を施す実験室で行われた。同じ部屋にプロジェクタで色を設定する。第一実験では出なかった、味と香りについての結果を得るためである。第二実験では、spiciness で Blue, GreenがWhite に対して優越、fruitiness において Green がRed Blueに対して優越、bitterness においてBlueがGreen Whiteに対して優越、となった。つまり、「味」に関しての結果はでたが、香りについては有意と認められる差は出なかった。

そこで、最後の第三実験となる。同じワインであるにもかかわらず、部屋の照明によって「香り」の印象が変化するか、を測定する目的で行う。この処理が圧巻である。どうしても差を出すために、同じ部屋で照明を２通りに変化させて１グラスづつ提供して、一対比較法で行う。モデルはBTL(Bradley-Terry-Luce)を採用する。ところが、重大な問題が発生する。これは官能実験でしばしば起こることなのだが、順序効果(Order-Effect)が強く生じている。つまり、ワイン提示の順序が、照明という要因を上回って影響している。この場合、同じワインであるが、初めに提供した方が照明にかかわらず選択されることが多くなっているのである。

そこで、BTLモデルのDavidson and Beaverによる拡張を用いることとする。この著者達は手法にも細心の注意を払っている。ANOVA の実行においては、アンバランスかつ異分散でも頑健なBrown-Forsythe プロセス、いわゆるF* test（残余平方和と第二自由度の調整）で行っている。t検定はもちろんWelchである。

順序効果としては例えば次の通り。第一グラスをBlue照明で、第二グラスをRed照明で飲んだ場合、26被験者がBlueの元での方がfruitierな香りであるとし、19被験者がRedの元での方がfruitierとする。照明が逆順の場合には、10被験者がBlue（第二グラス）を、35被験者がRed（第一グラス）をfruitierとした。このデータは原論文の Table 5.に提供されているので、追処理をおこなうことができる。このFruitiness 以外に、Preference, Spiciness, Sweetness が提示されている。それらも同様に入力する。完全のためにここに提供する。

preference <-
structure(c(0, 32, 22, 27, 0, 25, 23, 24, 0, 0, 13, 13, 18, 0, 
10, 12, 11, 0), .Dim = c(3L, 3L, 2L), .Dimnames = structure(list(
    `>` = c("Blue", "Red", "White"), `<` = c("Blue", "Red", "White"
    ), order = c("1", "2")), .Names = c(">", "<", "order")))
fruitiness <-
structure(c(0, 35, 21, 26, 0, 26, 22, 28, 0, 0, 10, 14, 19, 0, 
9, 13, 7, 0), .Dim = c(3L, 3L, 2L), .Dimnames = structure(list(
    `>` = c("Blue", "Red", "White"), `<` = c("Blue", "Red", "White"
    ), order = c("1", "2")), .Names = c(">", "<", "order")))
spiciness <-
structure(c(0, 21, 16, 24, 0, 18, 20, 17, 0, 0, 24, 19, 21, 0, 
17, 15, 18, 0), .Dim = c(3L, 3L, 2L), .Dimnames = structure(list(
    `>` = c("Blue", "Red", "White"), `<` = c("Blue", "Red", "White"
    ), order = c("1", "2")), .Names = c(">", "<", "order")))
sweetness <-
structure(c(0, 32, 18, 22, 0, 19, 21, 23, 0, 0, 13, 17, 23, 0, 
16, 14, 12, 0), .Dim = c(3L, 3L, 2L), .Dimnames = structure(list(
    `>` = c("Blue", "Red", "White"), `<` = c("Blue", "Red", "White"
    ), order = c("1", "2")), .Names = c(">", "<", "order")))
tpreference <- preference
tpreference[,,2] <- t(preference[,,2])
tfruitiness <- fruitiness
tfruitiness[,,2] <- t(fruitiness[,,2])
tspiciness  <- spiciness
tspiciness[,,2]  <- t(spiciness[,,2])
tsweetness  <- sweetness
tsweetness[,,2]  <- t(sweetness[,,2])

Fruitiness の全体を Table 5.にしたがって3dim array fruitinessに収めると次の通りである。ただし、のちの分析には　fruitier[,,2] を転置した、tfruitiness を用いる。tfruitierでいえば、orderはワイングラスの順番を意味しており、数字は行列の行を選択した人数を表している。

fruitiness
## , , order = 1
## 
##        <
## >       Blue Red White
##   Blue     0  26    22
##   Red     35   0    28
##   White   21  26     0
## 
## , , order = 2
## 
##        <
## >       Blue Red White
##   Blue     0  19    13
##   Red     10   0     7
##   White   14   9     0
tfruitiness
## , , order = 1
## 
##        <
## >       Blue Red White
##   Blue     0  26    22
##   Red     35   0    28
##   White   21  26     0
## 
## , , order = 2
## 
##        <
## >       Blue Red White
##   Blue     0  10    14
##   Red     19   0     9
##   White   13   7     0

BTL model の拡張は、もちろんRで行うことができる。パッケージ eba は EBA(Elimination-by-Aspects) model を対象とするが、拡張BTL modelを取り扱える。そもそもこのパッケージの著者 F.Wickelmaierは本論文で言及されている次の論文の著者である。かつては Matlabで提供していたが、2012からは他の関数も R で提供している。現在としては当然の措置であろう。

WICKELMAIER, F. and SCHMID, C. 2004. A Matlab function to estimate choice model parameters from paired-comparison data. Behav. Res. Methods Instrum. Comput. 36, 29-40. http://homepages.uni-tuebingen.de/florian.wickelmaier（著者のHP）

具体的なプログラムと結果は次の通りである。このグラフはTable 5. にそっくりになることを目指したので、細部の処理に手間がかかっているが、本質的なところはわずかである。原論文はRのlatticeを使ったかと思われる(Matlab でもlatticeのimplementはあったのかなとも思うが）。

このように、論文の中核部分でRならでは最新の処理ができている。グラフでは信頼区間(95%)を表示している。検定ではなく、信頼区間が重なるかどうかで判断れば、Fruitiness, Sweetness では Redが優越、Spicinessは他と傾向が異なることがわかる。

require(eba)
# Table 5.
vscol <- c("preference","fruitiness","spiciness","sweetness")
Vscol <- c("Preference","Fruitiness","Spiciness","Sweetness")
tvscol <- paste("t",vscol,sep="")
Norm.Estm <- array(0, c(3,4,4))
dimnames(Norm.Estm) <- list(c("Blue","Red","White"),
        c("Norm.Estimate", "Std.Error", "lower","upper"),　vscol)
names(dimnames(Norm.Estm)) <- c("ambient","value","response")

plotbtl <- function(EBA, vscol, tvscol)
{
  for (item in tvscol){
    command <- paste("ueos <- summary(",EBA,"(",item,"))")
    eval(parse(text=command))　　　# essential
    dat <- coef(ueos)/coef(ueos)[3,1]
    dat[,3] <- dat[,1]-dat[,2]
    dat[,4] <- dat[,1]+dat[,2]
    Norm.Estm[,,match(item, tvscol)] <- dat                   
    }
  
  ymin <- min(Norm.Estm[,3,])
  ymax <- max(Norm.Estm[,4,])
  opar <- par(mar=c(3,3,2,1), mfrow=c(2,2)) 
  # usually: par(mar=c(5,4,3,1), mfrow=c(1,1))
  
  xval <- 1:3
  for (cnt in 1:4){
    dname <- tvscol[cnt]
    plot(Norm.Estm[,1,cnt] ~ xval, type="n", ylim=c(ymin,ymax), 
         xlim=c(0.5, 3.5), xaxt="n", main=Vscol[cnt])
    abline(h=1, col="gray")
    arrows(xval, Norm.Estm[,3,cnt],
           xval, Norm.Estm[,4,cnt], .05, 90, 3, "black", lwd=2)
    points(Norm.Estm[,1,cnt]~ xval, type="p", pch=21, 
           bg="white", lwd=2,cex=1.5)
    axis(1, at=1:3, dimnames(Norm.Estm)[[1]])
    }
}

plotbtl("eba.order", vscol, tvscol)

この実例を通常の BTLmodel で処理するとどうなるであろう。それには、tpreference のorder=1,2を加えた行列 spreference を作成し、eba で処理すればよい。その結果は、全く同じに見える。

実はこの例の場合、ordereffect = 0.9854 であって、殆ど１である(H0: ordereffect = 1 の検定では p=.915)。よって普通のBTLモデルで十分。つまり、現実の order effect はかなり出ているのだが、相互作用というか、加え合わせると影響が薄まってしまった。パッケージにあるデータ heaviness では、ordereffect = 1.33734 であり、検定は p=.0035 で棄却となる。

spreference <- apply(tpreference, c(1,2), sum)
sfruitiness <- apply(tfruitiness, c(1,2), sum)
sspiciness <- apply(tspiciness, c(1,2), sum)
ssweetness <- apply(tsweetness, c(1,2), sum)

svscol <- paste("s",vscol,sep="")
plotbtl("eba", vscol, svscol)