一対比較法を用いた刺激評価には，呈示した2つの刺激のうちいずれか一方を選ばせる方法と，どちらにより近いと思うかを段階評価させる方法があり，前者の方法はサーストン法，後者の方法はシェフェ法と呼ばれます。最近，シェフェ法によるデータを分析する機会があったので，その際に作ったスクリプトなどを公開しておきます。
じつはシェフェ法にはいろいろとバリエーションがあるのですが，ここでは各評価者がすべての刺激対について評価した場合を扱います。なお，刺激の呈示位置（左右）による影響はない（あるいは少なくともそれが相殺されている）ものとします。
lm()を用いる方法
このようなデータを分析する方法として，手軽なのは青木先生のサイトで紹介されている方法でしょう。この方法では，それぞれの組み合わせごとにどの評定段階の回答がいく使ったかを集計したうえで，ダミー変数を用いて回帰分析を行います。
この手法については青木先生がスクリプトを公開されており，それを使用すれば分析できるわけですが，このスクリプトを実行して分析するためにはデータを集計表として与える必要があります。しかし，入力されている実験・調査のデータはたいてい一人分が1行に入力されているので，まずはデータを集計表の形に変換しなければなりません。
比較的規模が大きくて，集計表にした時に度数が０になるセルがないようなデータであればtable()を使って比較的簡単に集計表を作成できますが，集計表に度数０のセルがあるようなデータの場合にはそうも行きません。そこで，まずは入力データから集計表を作成するスクリプトを作成しました。
試しに、ここで公開されている資料の第5表のデータを使用して集計表を作成してみます。
X<-matrix(c(2,2,-3,3,0,-2,3,0,-1,2,1,-1,0,1,-3,1,0,-2),ncol=3,byrow=T)
通常，実験・調査のデータの場合，次のような形で入力されていることでしょう。行が評価者，列は刺激の組み合わせです。
> X
[,1] [,2] [,3]
[1,] 2 2 -3
[2,] 3 0 -2
[3,] 3 0 -1
[4,] 2 1 -1
[5,] 0 1 -3
[6,] 1 0 -2
ここから集計表を作成するには，次のスクリプトを使用します。
### 集計表の作成
scheffe.table<-function(Y,c,n,lv=c(-2,-1,0,1,2)){
# Y: データフレーム
# c: 評価対象の種類
# n: 評価者の人数
# lv: 評価値
## データを回答者ごとに分割
mat <- split(Y,sample(n))
## 集計表
nr=choose(c,2)
freq <- matrix(0,nr,length(lv))
colnames(freq) <- lv
tmp<-sapply(mat, function(x){
for(i in 1:nr){
freq[i,grep(paste('^',x[i],'$',sep=''),colnames(freq))]<<-freq[i,grep(paste('^',x[i],'$',sep=''),colnames(freq))]+1}
})
freq
}
これを次のように実行します1。
res <- scheffe.table(X,3,6,lv=c(3:-3))
するとこのようになります。
> res
3 2 1 0 -1 -2 -3
[1,] 2 2 1 1 0 0 0
[2,] 0 1 2 3 0 0 0
[3,] 0 0 0 0 2 2 2
青木先生のスクリプトで分析を行うには，この行列を「一対比較の結果を表す行列(A)」として用います。
lm()の結果を得る
集計表を作成して，そのままlm()の結果まで出したければ，先ほどのスクリプトをワークスペースに読み込んでおいた上で次のスクリプトを実行します2。
### lmの実行
scheffe.lm<-function(Y,c,n,lv,dir=1){
## 集計表を作成
freq <- scheffe.table(Y,c,n,lv)
x <- combn(c, 2)
nc <- ncol(x)
X <- matrix(0, nc, c)
X[cbind(1:nc, x[1,])] <- 1*dir
X[cbind(1:nc, x[2,])] <- -1*dir
df <- data.frame(y= freq%*%lv, X[,2:nc])
lm(y~.,data=df)
}
すると，次のような結果が得られます。
> scheffe.lm(X,3,6,lv=c(3:-3))
Call:
lm(formula = y ~ ., data = df)
Coefficients:
(Intercept) X1 X2
5 16 9
この結果の切片(Intercept)は１番目の刺激の係数で，変数名と刺激番号が1つずつずれてしまっていますが，この結果から１番目の刺激（Intercept:5）がもっとも値が小さく，２番目の刺激（X1:16）がもっとも大きな値ですので，好まれている方から「刺激2」→「刺激3」→「刺激1」となります。
平均選好度の算出と分散分析
先ほどの方法は比較的簡単にできるのですが，得られる結果の数値の幅が評価尺度の範囲（+3 ~ -3）とは全然違ったものになっていて，いまひとつ結果の解釈が難しいものになっています。そこで，別の方法でも分析をしてみます。
先ほどサンプルデータの参照に使用したこの資料では，一対比較のデータを元に分散分析表を作成する方法が紹介されています。この方法では，上で説明したlm()を使う方法とは違い，リーグ戦の勝敗表のような形(行数，列数の両方が刺激個数と同じ正方行列)で集計表を作る必要があります。
ここで活躍するのがupper.tri()という関数です。先ほどと同じデータの１行目は，評価者1が「刺激1-刺激2」，「刺激1-刺激3」，「刺激2-刺激3」のペアについて評価した結果です。これを必要な形の行列に変形してみましょう。
まずは「刺激数(3個) 」と同じ行数・列数の正方行列を作成します。行列の中身はすべて0としておきます。
m <- matrix(0,3,3)
この行列に，一人目の評価者のデータ(X[1,])を設定します。ここで，upper.tri()を使用します。
このupper.tri()は，そのセルが行列の上三角の成分の場合にはTRUEそうでなければFALSEを返します。これでデータを入力するセルを指定することで，上三角の成分にだけデータを入力することができます。
m[upper.tri(m)] <- X[1,]
実行した結果，次のような行列ができました3。
> m
[,1] [,2] [,3]
[1,] 0 2 2
[2,] 0 0 -3
[3,] 0 0 0
下三角部分が空欄のままですので，行列を転置して下半分の値を設定します。
m <- m - t(m)
評価者一人分についてだけですが，分析用の行列が完成しました。
> m
[,1] [,2] [,3]
[1,] 0 2 2
[2,] -2 0 -3
[3,] -2 3 0
この表の列の平均は，その評価者の各評価対象に対する選好度となります。この評価者の場合，[,1]の列の平均は(0+(-2)+(-2))/3 = -1.33，[,2]の平均は(2+0+3)/3 = 1.67，[,3]の平均は(2+(-3)+0)/3= -0.33で，2番目の対象を最も好み，1番目の対象を最も好まないということがわかります。そして，それぞれの評価対象についてこれを評価者全員で平均したものが，各評価対象の平均選好度になります。そして，それらの値に統計的な差があるかどうかの検定には，分散分析が用いられます。
評価者全体での各対象の平均選好度の算出と，それについての分散分析，および多重比較を行うスクリプトは以下の通りです。scheffe.anova(データフレーム, 刺激の数, 評価者の人数)の形で実行すると，各対象の平均選好度，分散分析表，そして多重比較の結果が出力されます。
サンプルのデータ(X)を使って実行した結果は次の通りです。
> scheffe.anova(X,3,6)
$`Mean Preferences`
1 2 3
-0.8333333 1.277778 -0.4444444
$`Analysis of Variance`
Df Sum Sq Mean Sq F Value PR(>F) Sig.
Main Effect 2 45.444444 22.7222222 12.7018634 0.01096751 *
Main * Subj 10 5.222222 0.5222222 0.2919255 0.95386007
Combination 1 1.388889 1.3888889 0.7763975 0.41858023
Error 5 8.944444 1.7888889 NA NA
Total 18 61.000000 NA NA NA
$`Significance Threshold`
0.10 0.05 0.01 0.001
criterion 1.17182 1.450695 2.1991 3.679442
$`Paired Comparisons`
Contrast Difference LCI.95 UCI.95 LCI.99 UCI.99 Sig.
1 1 vs 2 -2.1111111 -3.561806 -0.660416 -4.310211 0.08798862 *
2 1 vs 3 -0.3888889 -3.561806 -0.660416 -4.310211 0.08798862
3 2 vs 3 1.7222222 -1.839584 1.061806 -2.587989 1.81021085 *
この結果から，2番目の刺激が最も好まれ(1.27)，1番目の刺激が最も好まれない(-0.83)ことがわかります。また，分散分析の結果，主効果（刺激の違い）が有意で，多重比較の結果からは，刺激1と2，刺激2と3の間に5%水準で差が見られ，刺激1と3には有意な差がないことがわかります。
#### ベクトルデータから上側の三角行列を作成
up.tri.mat <- function(x, n){
stopifnot(length(x) == (n-1)*n/2)
mat <- matrix(0,n,n)
mat[upper.tri(mat,diag=F)]<- x
mat
}
### 各刺激評定の個人ごとの合計値
sum.scheffe <- function(x, n){
mat <- up.tri.mat(x,n)
mat <- mat - t(mat)
apply(mat,2,sum)
}
### 各刺激評定の個人ごとの下三角の平方和
sumsq.scheffe.lowsq <- function(x, n){
mat <- low.tri.mat(x,n)
apply(mat^2,2,sum)
}
### 刺激対ごとの平均値の差の算出
calc.pw.diff <- function(x,y,m){
m[x] - m[y]
}
### 各刺激の平均評定値
scheffe.mean<-function(Y,c,n){
z <- apply(Y,1,sum.scheffe,c)
z.sum <- apply(z,1,sum)
z.sum/(c*n)
}
### 評定値に対するANOVA
scheffe.anova<-function(Y,c,n){
z <- apply(Y,1,sum.scheffe,c) # sum.scheffe(x,c) 評定値(x)を行列に格納し，各刺激に対する個人ごとの評定合計を求める
z.mean <- as.matrix(scheffe.mean(Y,c,n)) # scheffe.mean(x,c,n) 回答者全体での各刺激の評定平均を求める
z.sum <- apply(z,1,sum) # 回答者ごとの評定の総和
z.sumsq <- apply(z^2,1,sum) # 回答者ごとの評定の平方和
z.lowsq <- apply(Y,1,sumsq.scheffe.lowsq,c) # 回答者ごとの下三角の平方和
z.low <- apply(Y,1,low.tri.mat,c) # 全員分を多し合わせた下三角行列の作成
### 平方和と自由度の算出
S.a <- sum(z.sum^2)/(c*n); df.a <- c-1
S.aB <- sum(z.sumsq)/c - S.a; df.aB <- (c-1)*(n-1)
S.g <- sum(apply(z.low,1,sum)^2)/n - S.a; df.g <- (c-1)*(c-2)/2
S.T <- sum(z.lowsq); df.T <- c*n*(c-1)/2;
S.e <- S.T - S.a - S.aB - S.g; df.e <- (c-1)*(c-2)*(n-1)/2
### 平均平方とF，p値の算出
SS <- c(S.a, S.aB, S.g, S.e, S.T);
df <- c(df.a, df.aB, df.g, df.e, df.T)
MS <- (SS/df)[1:4]
Fv <- MS[1:3]/MS[4]
p <- mapply(pf,Fv,df[1:3],df.e, lower.tail=FALSE)
sig <- sapply(p, function(x) {
if(x < .001) '***'
else if(x < .01) '** '
else if(x < .05) '* '
else if(x < .10) '. '
else ' '
})
### ここから先は結果の出力用
# 評定平均
colnames(z.mean)<-''; rownames(z.mean)<-c(1:c)
# 分散分析表
# 行数の足りない部分はとりあえず'NA'とする
length(MS)<-length(SS); length(Fv)<-length(SS); length(p)<-length(SS)
# 数値だけを連結してデータフレームに変換
avtab <- as.data.frame(cbind(df,SS,MS,Fv,p));
# 文字を含む列をデータフレームに追加
avtab$sig <-c(sig,'','')
colnames(avtab)<-c('Df','Sum Sq','Mean Sq', 'F Value', 'PR(>F)', 'Sig.')
rownames(avtab)<-c('Main Effect','Main * Subj','Combination', 'Error', 'Total')
# 差の臨界値
thres <- matrix(c(qtukey(.10,c,df.e,lower.tail = F),
qtukey(.05,c,df.e,lower.tail = F),
qtukey(.01,c,df.e,lower.tail = F),
qtukey(.001,c,df.e,lower.tail = F)), ncol=4) * sqrt(MS[4]/(c*n))
colnames(thres) <- c('0.10','0.05','0.01','0.001'); rownames(thres)<-c('criterion')
# 一対比較
cb <- combn(c(1:c),2)
pw.cont <- sapply(c(1:ncol(cb)), function(x){paste(cb[1,x],cb[2,x],sep=' vs ')})
pw.diff <- sapply(c(1:ncol(cb)), function(x){z.mean[cb[1,x]] - z.mean[cb[2,x]]})
pw.95cl <- sapply(c(1:ncol(cb)), function(x){pw.diff[cb[1,x]] - thres[2]})
pw.95cu <- sapply(c(1:ncol(cb)), function(x){pw.diff[cb[1,x]] + thres[2]})
pw.99cl <- sapply(c(1:ncol(cb)), function(x){pw.diff[cb[1,x]] - thres[3]})
pw.99cu <- sapply(c(1:ncol(cb)), function(x){pw.diff[cb[1,x]] + thres[3]})
pw.sig <- sapply(pw.diff, function(x) {
if(abs(x) > thres[4]) '***'
else if(abs(x) > thres[3]) '** '
else if(abs(x) > thres[2]) '* '
else if(abs(x) > thres[1]) '. '
else ' '
})
pw.diff<-data.frame('Contrast'=pw.cont, 'Difference'=pw.diff,
'LCI.95'=pw.95cl, 'UCI.95'=pw.95cu,
'LCI.99'=pw.99cl, 'UCI.99'=pw.99cu,
'Sig.'=pw.sig)
# 結果の表示（構造体としてまとめた上で出力）
return(structure(list( "Mean Preferences"=t(z.mean),
"Analysis of Variance"=avtab,
"Significance Threshold"=thres,
"Paired Comparisons"=pw.diff
)))
}
-
青木先生のサイトにあるスクリプトでは，左側が正の値，右側が負の値の向きになっていますので，それに揃えます。 ↩
-
左側が負の値，右側が正の値になっているデータでは，そのまま実行すると結果の符号が逆になってしまいます。その場合には，データの符号を逆転させてから分析を実行するか，スクリプト実行時に
dir=-1オプションを指定してください。 ↩ -
資料とは表の向きが異なっていますが，これで問題ありません。資料の第1図では尺度の左側に負の値が置かれており，より好まれるものが負の値になるような集計がなされています。しかし，より好まれるものがプラスの値になる方が直感的にはわかりやすいため，ここでは好まれる方が正の値になるように集計しています。 ↩