| Title: | Elimination-by-Aspects Models |
|---|---|
| Description: | Fitting and testing multi-attribute probabilistic choice models, especially the Bradley-Terry-Luce (BTL) model (Bradley & Terry, 1952 <doi:10.1093/biomet/39.3-4.324>; Luce, 1959), elimination-by-aspects (EBA) models (Tversky, 1972 <doi:10.1037/h0032955>), and preference tree (Pretree) models (Tversky & Sattath, 1979 <doi:10.1037/0033-295X.86.6.542>). |
| Authors: | Florian Wickelmaier [aut, cre] |
| Maintainer: | Florian Wickelmaier <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 1.10-0 |
| Built: | 2024-11-26 06:24:40 UTC |
| Source: | CRAN |
Creates a (completely) balanced paired-comparison design.
balanced.pcdesign(nstimuli)balanced.pcdesign(nstimuli)
nstimuli |
number of stimuli in the paired-comparison design |
When nstimuli is odd, the presentation order is completely balanced,
that is any given stimulus appears an equal number of times as the first
and second member of a pair. When nstimuli is even, the presentation
order is balanced as much as possible. Subjects should be equally assigned
to listA and listB for the purpose of balancing the
within-pair presentation order across a sample of subjects. Pairs should be
re-randomized for each subject. See David (1988) for details.
pairs |
a character array holding the balanced pairs |
listA |
the vector pairs in the original within-pair order |
listB |
the vector of pairs in the inverted within-pair order |
David, H. (1988). The method of paired comparisons. London: Griffin.
## Create balanced design for 6 stimuli bp <- balanced.pcdesign(6) ## Replicate each within-pair order 10 times and re-randomize cbind(replicate(10, sample(bp$listA)), replicate(10, sample(bp$listB)))## Create balanced design for 6 stimuli bp <- balanced.pcdesign(6) ## Replicate each within-pair order 10 times and re-randomize cbind(replicate(10, sample(bp$listA)), replicate(10, sample(bp$listB)))
Performs a bootstrap by resampling the individual data matrices.
boot(D, R = 100, A = 1:I, s = rep(1/J, J), constrained = TRUE)boot(D, R = 100, A = 1:I, s = rep(1/J, J), constrained = TRUE)
D |
either a 3d array consisting of the individual paired
comparison matrices or an object of class |
R |
the number of bootstrap samples |
A |
a list of vectors consisting of the stimulus aspects;
the default is |
s |
the starting vector with default |
constrained |
logical, if TRUE (default), parameters are constrained to be positive |
The bootstrap functions eba.boot.constrained and eba.boot
are automatically called by boot.
The code is experimental and may change in the future.
p |
the matrix of bootstrap vectors |
stat |
the matrix of bootstrap statistics, including parameter means, standard errors, and confidence limits |
data(pork) # pork tasting data, 10 individual paired comparison matrices btl1 <- eba(apply(pork, 1:2, sum)) # fit Bradley-Terry-Luce model b <- boot(pork, R = 200) # resample 200 times plot(coef(btl1), b$stat[, "mean"], log = "xy") abline(0, 1, lty = 2)data(pork) # pork tasting data, 10 individual paired comparison matrices btl1 <- eba(apply(pork, 1:2, sum)) # fit Bradley-Terry-Luce model b <- boot(pork, R = 200) # resample 200 times plot(coef(btl1), b$stat[, "mean"], log = "xy") abline(0, 1, lty = 2)
Rumelhart and Greeno (1971) presented 234 participants with pairs of names of nine celebrities: the politicians L. B. Johnson (LBJ), Harold Wilson (HW), and Charles De Gaulle (CDG); the athletes Johnny Unitas (JU), Carl Yastrzemski (CY), and A. J. Foyt (AJF); the actresses Brigitte Bardot (BB), Elizabeth Taylor (ET), and Sophia Loren (SL). Participants were instructed to choose the person with whom they would rather spend an hour of discussion.
data(celebrities)data(celebrities)
A square data frame containing the absolute choice frequencies and a diagonal of zeros; row stimuli are chosen over column stimuli.
Rumelhart, D.L., & Greeno, J.G. (1971). Similarity between stimuli: An experimental test of the Luce and Restle choice models. Journal of Mathematical Psychology, 8, 370–381. doi:10.1016/0022-2496(71)90038-1
data(celebrities) celebrities["LBJ", "HW"] # 159 participants chose Johnson over Wilsondata(celebrities) celebrities["LBJ", "HW"] # 159 participants chose Johnson over Wilson
Number of circular triads and coefficient of consistency.
circular(mat, alternative = c("two.sided", "less", "greater"), exact = NULL, correct = TRUE, simulate.p.value = FALSE, nsim = 2000)circular(mat, alternative = c("two.sided", "less", "greater"), exact = NULL, correct = TRUE, simulate.p.value = FALSE, nsim = 2000)
mat |
a square matrix or a data frame consisting of (individual) binary choice data; row stimuli are chosen over column stimuli |
alternative |
a character string specifying the alternative hypothesis,
must be one of |
exact |
a logical indicating whether an exact p-value should be computed |
correct |
a logical indicating whether to apply continuity correction in the chi-square approximation for the p-value |
simulate.p.value |
a logical indicating whether to compute p-values by Monte Carlo simulation |
nsim |
an integer specifying the number of replicates used in the Monte Carlo test |
Kendall's coefficient of consistency,
lies between one (perfect consistency) and zero,
where T is the observed number of circular triads,
and the maximum possible number of circular triads is
, if is even, and
else, and is the
number of stimuli or objects being judged. For details see Kendall and
Babington Smith (1940) and David (1988).
Kendall (1962) discusses a test of the hypothesis that the number of
circular triads T is different (smaller or greater) than expected
when choosing randomly. For small , an exact p-value is computed,
based on the exact distributions listed in Alway (1962) and in Kendall
(1962). Otherwise, an approximate chi-square test is computed. In this
test, the sampling distribution is measured from lower to higher values of
T, so that the probability that T will be exceeded is the
complement of the probability for chi2. The chi-square approximation
may be incorrect if and is only available for .
T |
number of circular triads |
T.max |
maximum possible number of circular triads |
T.exp |
expected number of circular triads |
zeta |
Kendall's coefficient of consistency |
chi2, df, correct
|
the chi-square statistic and degrees of freedom for the approximate test, and whether continuity correction has been applied |
p.value |
the p-value for the test (see Details) |
simulate.p.value, nsim
|
whether the p-value is based on simulations, number of simulation runs |
Alway, G.G. (1962). The distribution of the number of circular triads in paired comparisons. Biometrika, 49, 265–269. doi:10.1093/biomet/49.1-2.265
David, H. (1988). The method of paired comparisons. London: Griffin.
Kendall, M.G. (1962). Rank correlation methods. London: Griffin.
Kendall, M.G., & Babington Smith, B. (1940). On the method of paired comparisons. Biometrika, 31, 324–345. doi:10.1093/biomet/31.3-4.324
# A dog's preferences for six samples of food # (Kendall and Babington Smith, 1940, p. 326) dog <- matrix(c(0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0), 6, 6, byrow=TRUE) dimnames(dog) <- setNames(rep(list(c("meat", "biscuit", "chocolate", "apple", "pear", "cheese")), 2), c(">", "<")) circular(dog, alternative="less") # moderate consistency subset(strans(dog)$violdf, !wst) # list circular triads# A dog's preferences for six samples of food # (Kendall and Babington Smith, 1940, p. 326) dog <- matrix(c(0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0), 6, 6, byrow=TRUE) dimnames(dog) <- setNames(rep(list(c("meat", "biscuit", "chocolate", "apple", "pear", "cheese")), 2), c(">", "<")) circular(dog, alternative="less") # moderate consistency subset(strans(dog)$violdf, !wst) # list circular triads
Computes the (normalized) covariance matrix of the utility scale from the covariance matrix of elimination-by-aspects (EBA) model parameters.
cov.u(object, norm = "sum")cov.u(object, norm = "sum")
object |
an object of class |
norm |
either |
The additivity rule for covariances
is used for the computations.
If norm is not NULL, the unnormalized covariance matrix is
transformed using , where the constant results from
the type of normalization applied.
The (normalized) covariance matrix of the utility scale.
In summer 2007, a survey was conducted at the Department of Psychology, University of Tuebingen. Hundred and ninety-two participants were presented with all 15 unordered pairs of the names of six drugs or substances and asked to choose the drug they judged as more dangerous for their health. The six drugs were alcohol (alc), tobacco (tob), cannabis (can), ecstasy (ecs), heroine (her), and cocaine (coc). Choice frequencies were aggregated in four groups defined by gender and age.
data(drugrisk)data(drugrisk)
A 3d array consisting of four square matrices of choice frequencies (row drugs are judged over column drugs):
drugrisk[, , group = "female30"]holds the choices of the 48 female participants up to 30 years of age.
drugrisk[, , group = "female31"]holds the choices of the 48 female participants from 31 years of age.
drugrisk[, , group = "male30"]holds the choices of the 48 male participants up to 30 years of age.
drugrisk[, , group = "male31"]holds the choices of the 48 male participants from 31 years of age.
Wickelmaier, F. (2008). Analyzing paired-comparison data in R using probabilistic choice models. Presented at the R User Conference 2008, August 12-14, Dortmund, Germany.
data(drugrisk) ## Bradley-Terry-Luce (BTL) model btl <- eba(drugrisk[, , group = "male30"]) ## Elimination-by-aspects (EBA) model, 1 additional aspect A1 <- list(c(1), c(2,7), c(3,7), c(4,7), c(5,7), c (6,7)) eba1 <- eba(drugrisk[, , group = "male30"], A1) ## EBA model, 2 additional aspects A2 <- list(c(1), c(2,7), c(3,7), c(4,7,8), c(5,7,8), c(6,7,8)) eba2 <- eba(drugrisk[, , group = "male30"], A2) ## Model selection anova(btl, eba1, eba2) ## Utility scale values (BTL for females, EBA for males) dotchart(cbind( apply(drugrisk[, , 1:2], 3, function(x) uscale(eba(x), norm = 1)), apply(drugrisk[, , 3:4], 3, function(x) uscale(eba(x, A2), norm = 1)) ), xlab="Utility scale value (BTL and EBA models)", main="Perceived health risk of drugs", panel.first = abline(v = 1, col = "gray"), log = "x") mtext("(Wickelmaier, 2008)", line = 0.5)data(drugrisk) ## Bradley-Terry-Luce (BTL) model btl <- eba(drugrisk[, , group = "male30"]) ## Elimination-by-aspects (EBA) model, 1 additional aspect A1 <- list(c(1), c(2,7), c(3,7), c(4,7), c(5,7), c (6,7)) eba1 <- eba(drugrisk[, , group = "male30"], A1) ## EBA model, 2 additional aspects A2 <- list(c(1), c(2,7), c(3,7), c(4,7,8), c(5,7,8), c(6,7,8)) eba2 <- eba(drugrisk[, , group = "male30"], A2) ## Model selection anova(btl, eba1, eba2) ## Utility scale values (BTL for females, EBA for males) dotchart(cbind( apply(drugrisk[, , 1:2], 3, function(x) uscale(eba(x), norm = 1)), apply(drugrisk[, , 3:4], 3, function(x) uscale(eba(x, A2), norm = 1)) ), xlab="Utility scale value (BTL and EBA models)", main="Perceived health risk of drugs", panel.first = abline(v = 1, col = "gray"), log = "x") mtext("(Wickelmaier, 2008)", line = 0.5)
Fits a (multi-attribute) probabilistic choice model by maximum likelihood.
eba(M, A = 1:I, s = rep(1/J, J), constrained = TRUE) OptiPt(M, A = 1:I, s = rep(1/J, J), constrained = TRUE) ## S3 method for class 'eba' summary(object, ...) ## S3 method for class 'eba' anova(object, ..., test = c("Chisq", "none"))eba(M, A = 1:I, s = rep(1/J, J), constrained = TRUE) OptiPt(M, A = 1:I, s = rep(1/J, J), constrained = TRUE) ## S3 method for class 'eba' summary(object, ...) ## S3 method for class 'eba' anova(object, ..., test = c("Chisq", "none"))
M |
a square matrix or a data frame consisting of absolute choice frequencies; row stimuli are chosen over column stimuli |
A |
a list of vectors consisting of the stimulus aspects;
the default is |
s |
the starting vector with default |
constrained |
logical, if TRUE (default), parameters are constrained to be positive |
object |
an object of class |
test |
should the p-values of the chi-square distributions be reported? |
... |
additional arguments; none are used in the summary method;
in the anova method they refer to additional objects of class |
eba is a wrapper function for OptiPt. Both functions can be
used interchangeably. See Wickelmaier and Schmid (2004) for further
details.
The probabilistic choice models that can be fitted to paired-comparison data are the Bradley-Terry-Luce (BTL) model (Bradley, 1984; Luce, 1959), preference tree (Pretree) models (Tversky and Sattath, 1979), and elimination-by-aspects (EBA) models (Tversky, 1972), the former being special cases of the latter.
A represents the family of aspect sets. It is usually a list of
vectors, the first element of each being a number from 1 to I;
additional elements specify the aspects shared by several stimuli. A
must have as many elements as there are stimuli. When fitting a BTL model,
A reduces to 1:I (the default), i.e. there is only one aspect
per stimulus.
The maximum likelihood estimation of the parameters is carried out by
nlm. The Hessian matrix, however, is approximated by
nlme::fdHess. The likelihood functions L.constrained and
L are called automatically.
See group.test for details on the likelihood ratio
tests reported by summary.eba.
coefficients |
a vector of parameter estimates |
estimate |
same as |
logL.eba |
the log-likelihood of the fitted model |
logL.sat |
the log-likelihood of the saturated (binomial) model |
goodness.of.fit |
the goodness of fit statistic including the likelihood ratio fitted vs. saturated model (-2logL), the degrees of freedom, and the p-value of the corresponding chi-square distribution |
u.scale |
the unnormalized utility scale of the stimuli; each utility scale value is defined as the sum of aspect values (parameters) that characterize a given stimulus |
hessian |
the Hessian matrix of the likelihood function |
cov.p |
the covariance matrix of the model parameters |
chi.alt |
the Pearson chi-square goodness of fit statistic |
fitted |
the fitted paired-comparison matrix |
y1 |
the data vector of the upper triangle matrix |
y0 |
the data vector of the lower triangle matrix |
n |
the number of observations per pair ( |
mu |
the predicted choice probabilities for the upper triangle |
nobs |
the number of pairs |
Florian Wickelmaier
Bradley, R.A. (1984). Paired comparisons: Some basic procedures and examples. In P.R. Krishnaiah & P.K. Sen (eds.), Handbook of Statistics, Volume 4. Amsterdam: Elsevier. doi:10.1016/S0169-7161(84)04016-5
Luce, R.D. (1959). Individual choice behavior: A theoretical analysis. New York: Wiley.
Tversky, A. (1972). Elimination by aspects: A theory of choice. Psychological Review, 79, 281–299. doi:10.1037/h0032955
Tversky, A., & Sattath, S. (1979). Preference trees. Psychological Review, 86, 542–573. doi:10.1037/0033-295X.86.6.542
Wickelmaier, F., & Schmid, C. (2004). A Matlab function to estimate choice model parameters from paired-comparison data. Behavior Research Methods, Instruments, and Computers, 36, 29–40. doi:10.3758/BF03195547
strans, uscale, cov.u,
group.test, wald.test, plot.eba,
residuals.eba, logLik.eba,
simulate.eba,
kendall.u, circular, trineq,
thurstone, nlm.
data(celebrities) # absolute choice frequencies btl1 <- eba(celebrities) # fit Bradley-Terry-Luce model A <- list(c(1,10), c(2,10), c(3,10), c(4,11), c(5,11), c(6,11), c(7,12), c(8,12), c(9,12)) # the structure of aspects eba1 <- eba(celebrities, A) # fit elimination-by-aspects model summary(eba1) # goodness of fit plot(eba1) # residuals versus predicted values anova(btl1, eba1) # model comparison based on likelihoods confint(eba1) # confidence intervals for parameters uscale(eba1) # utility scale ci <- 1.96 * sqrt(diag(cov.u(eba1))) # 95% CI for utility scale values dotchart(uscale(eba1), xlim=c(0, .3), main="Choice among celebrities", xlab="Utility scale value (EBA model)", pch=16) # plot the scale arrows(uscale(eba1)-ci, 1:9, uscale(eba1)+ci, 1:9, .05, 90, 3) # error bars abline(v=1/9, lty=2) # indifference line mtext("(Rumelhart and Greeno, 1971)", line=.5) ## See data(package = "eba") for application examples.data(celebrities) # absolute choice frequencies btl1 <- eba(celebrities) # fit Bradley-Terry-Luce model A <- list(c(1,10), c(2,10), c(3,10), c(4,11), c(5,11), c(6,11), c(7,12), c(8,12), c(9,12)) # the structure of aspects eba1 <- eba(celebrities, A) # fit elimination-by-aspects model summary(eba1) # goodness of fit plot(eba1) # residuals versus predicted values anova(btl1, eba1) # model comparison based on likelihoods confint(eba1) # confidence intervals for parameters uscale(eba1) # utility scale ci <- 1.96 * sqrt(diag(cov.u(eba1))) # 95% CI for utility scale values dotchart(uscale(eba1), xlim=c(0, .3), main="Choice among celebrities", xlab="Utility scale value (EBA model)", pch=16) # plot the scale arrows(uscale(eba1)-ci, 1:9, uscale(eba1)+ci, 1:9, .05, 90, 3) # error bars abline(v=1/9, lty=2) # indifference line mtext("(Rumelhart and Greeno, 1971)", line=.5) ## See data(package = "eba") for application examples.
Fits a (multi-attribute) probabilistic choice model that accounts for the effect of the presentation order within a pair.
eba.order(M1, M2 = NULL, A = 1:I, s = c(rep(1/J, J), 1), constrained = TRUE) ## S3 method for class 'eba.order' summary(object, ...)eba.order(M1, M2 = NULL, A = 1:I, s = c(rep(1/J, J), 1), constrained = TRUE) ## S3 method for class 'eba.order' summary(object, ...)
M1, M2
|
two square matrices or data frames consisting of absolute choice frequencies in both within-pair orders; row stimuli are chosen over column stimuli. If M2 is empty (default), M1 is assumed to be a 3d array containing both orders |
A |
see |
s |
the starting vector with default |
constrained |
see |
object |
an object of class |
... |
additional arguments |
The choice models include a single multiplicative order effect,
order, that is constant for all pairs (see Davidson and Beaver,
1977). An order effect < 1 (> 1) indicates a bias in favor of the first
(second) interval. See eba for choice models without order
effect.
Several likelihood ratio tests are performed (see also
summary.eba).
EBA.order tests an order-effect EBA model against a saturated
binomial model; this corresponds to a goodness of fit test of the former
model.
Order tests an EBA model with an order effect constrained to 1
against an unconstrained order-effect EBA model; this corresponds to a test
of the order effect.
Effect tests an order-effect indifference model (where all scale
values are equal, but the order effect is free) against the order-effect EBA
model; this corresponds to testing for a stimulus effect; order0 is
the estimate of the former model.
Wickelmaier and Choisel (2006) describe a model that generalizes the Davidson-Beaver model and allows for an order effect in Pretree and EBA models.
coefficients |
a vector of parameter estimates, the last component holds the order-effect estimate |
estimate |
same as |
logL.eba |
the log-likelihood of the fitted model |
logL.sat |
the log-likelihood of the saturated (binomial) model |
goodness.of.fit |
the goodness of fit statistic including the likelihood ratio fitted vs. saturated model (-2logL), the degrees of freedom, and the p-value of the corresponding chi-square distribution |
u.scale |
the unnormalized utility scale of the stimuli; each utility scale value is defined as the sum of aspect values (parameters) that characterize a given stimulus |
hessian |
the Hessian matrix of the likelihood function |
cov.p |
the covariance matrix of the model parameters |
chi.alt |
the Pearson chi-square goodness of fit statistic |
fitted |
3d array of the fitted paired-comparison matrices |
y1 |
the data vector of the upper triangle matrices |
y0 |
the data vector of the lower triangle matrices |
n |
the number of observations per pair ( |
mu |
the predicted choice probabilities for the upper triangles |
M1, M2
|
the data matrices |
Florian Wickelmaier
Davidson, R.R., & Beaver, R.J. (1977). On extending the Bradley-Terry model to incorporate within-pair order effects. Biometrics, 33, 693–702.
Wickelmaier, F., & Choisel, S. (2006). Modeling within-pair order effects in paired-comparison judgments. In D.E. Kornbrot, R.M. Msetfi, & A.W. MacRae (eds.), Fechner Day 2006. Proceedings of the 22nd Annual Meeting of the International Society for Psychophysics (p. 89–94). St. Albans, UK: The ISP.
eba, group.test, plot.eba,
residuals.eba, logLik.eba.
data(heaviness) # weights judging data ebao1 <- eba.order(heaviness) # Davidson-Beaver model summary(ebao1) # goodness of fit plot(ebao1) # residuals versus predicted values confint(ebao1) # confidence intervals for parametersdata(heaviness) # weights judging data ebao1 <- eba.order(heaviness) # Davidson-Beaver model summary(ebao1) # goodness of fit plot(ebao1) # residuals versus predicted values confint(ebao1) # confidence intervals for parameters
Zimmer et al. (2004) investigated the auditory unpleasantness of twelve short binaural recordings (Johannsen and Prante, 2001); recordings were presented via headphones to 74 participants.
data(envirosound)data(envirosound)
A data frame containing 74 observations on 2 variables:
paired comparison of class paircomp;
judgments for all 66 paired comparisons from 12 recordings:
circular saw, stadium, dentist's drill, waterfall, ship's horn, stone in
well, typewriter, hooves, fan, howling wind, tyre on gravel, wasp.
median response time.
Details of the recordings, including psychoacoustic metrics, are available
as an attribute of the unpleasantness variable (see Examples).
Zimmer, K., Ellermeier, W., & Schmid, C. (2004). Using probabilistic choice models to investigate auditory unpleasantness. Acta Acustica united with Acustica, 90(6), 1019–1028.
Johannsen, K., & Prante, H.U. (2001). Environmental sounds for psychoacoustic testing. Acta Acustica united with Acustica, 87(2), 290–293.
requireNamespace("psychotools") data(envirosound) set.seed(1019) ## Choice-model representation of unpleasantness mat <- summary(envirosound$unpleasantness, pcmatrix = TRUE) strans(mat) btl1 <- eba(mat) eba1 <- eba(mat, A = list(c(1, 13), c(2, 13), c(3, 13), c(4, 13), c(5, 13), c(6, 13), c(7, 13), c(8, 13), c(9, 13), c(10, 13), c(11, 13), 12)) eba2 <- eba(mat, A = list(c(1, 13), c(2, 13), c(3, 13), c(4, 13), c(5, 13), c(6, 13), c(7, 13, 14), c(8, 13, 14), c(9, 13, 14), c(10, 13, 14), c(11, 13, 14), 12), s = runif(14)) anova(btl1, eba1, eba2) sounds <- psychotools::covariates(envirosound$unpleasantness) sounds$u <- 10 * uscale(eba2, norm = 9) # u(fan) := 10 plot(magnitude ~ u, sounds, log = "x", type = "n", xlab = "Indirect scaling (EBA model)", ylab = "Direct magnitude estimation", main = "Auditory unpleasantness of environmental sound") mtext("(Zimmer et al., 2004)", line = 0.5) abline(lm(magnitude ~ log10(u), sounds)) text(magnitude ~ u, sounds, labels = abbreviate(rownames(sounds), 4)) ## Predicting unpleasantness from psychoacoustic metrics summary( lm(log(u) ~ scale(sharpness, scale = FALSE) + scale(roughness, scale = FALSE):I(loudness.5 > 27), sounds[-12, ]) # w/o wasp )requireNamespace("psychotools") data(envirosound) set.seed(1019) ## Choice-model representation of unpleasantness mat <- summary(envirosound$unpleasantness, pcmatrix = TRUE) strans(mat) btl1 <- eba(mat) eba1 <- eba(mat, A = list(c(1, 13), c(2, 13), c(3, 13), c(4, 13), c(5, 13), c(6, 13), c(7, 13), c(8, 13), c(9, 13), c(10, 13), c(11, 13), 12)) eba2 <- eba(mat, A = list(c(1, 13), c(2, 13), c(3, 13), c(4, 13), c(5, 13), c(6, 13), c(7, 13, 14), c(8, 13, 14), c(9, 13, 14), c(10, 13, 14), c(11, 13, 14), 12), s = runif(14)) anova(btl1, eba1, eba2) sounds <- psychotools::covariates(envirosound$unpleasantness) sounds$u <- 10 * uscale(eba2, norm = 9) # u(fan) := 10 plot(magnitude ~ u, sounds, log = "x", type = "n", xlab = "Indirect scaling (EBA model)", ylab = "Direct magnitude estimation", main = "Auditory unpleasantness of environmental sound") mtext("(Zimmer et al., 2004)", line = 0.5) abline(lm(magnitude ~ log10(u), sounds)) text(magnitude ~ u, sounds, labels = abbreviate(rownames(sounds), 4)) ## Predicting unpleasantness from psychoacoustic metrics summary( lm(log(u) ~ scale(sharpness, scale = FALSE) + scale(roughness, scale = FALSE):I(loudness.5 > 27), sounds[-12, ]) # w/o wasp )
Tests for group effects in elimination-by-aspects (EBA) models.
group.test(groups, A = 1:I, s = rep(1/J, J), constrained = TRUE)group.test(groups, A = 1:I, s = rep(1/J, J), constrained = TRUE)
groups |
a 3d array containing one aggregate choice matrix per group |
A |
a list of vectors consisting of the stimulus aspects;
the default is |
s |
the starting vector with default |
constrained |
logical, if TRUE (default), EBA parameters are constrained to be positive |
The five tests are all based on likelihood ratios.
Overall compares a 1-parameter Poisson model to a saturated Poisson
model, thereby testing the equality of the frequencies in each cell of the
array. This test corresponds to simultaneously testing for a null effect of
(1) the context induced by a given pair, (2) the grouping factor, (3) the
stimuli, and (4) the imbalance between pairs. The deviances of the
remaining tests sum to the total deviance associated with the overall test.
EBA.g tests an EBA group model against a saturated binomial group
model, which corresponds to a goodness of fit test of the EBA group model.
Group tests an EBA model having its parameters restricted to be equal
across groups (single set of parameters) against the EBA group model
allowing its parameters to vary freely across groups (one set of parameters
per group); this corresponds to testing for group differences.
Effect tests an indifference model (where all choice probabilities
are equal to 0.5) against the restricted EBA model (single set of
parameters), which corresponds to testing for a stimulus effect.
Imbalance tests for differences in the number of observations per
pair by comparing the average sample size (1-parameter Poisson model) to the
actual sample sizes (saturated Poisson model).
See Duineveld, Arents, and King (2000) for further details, and Choisel and Wickelmaier (2007) for an application.
tests |
a table displaying the likelihood ratio test statistics |
Choisel, S., & Wickelmaier, F. (2007). Evaluation of multichannel reproduced sound: Scaling auditory attributes underlying listener preference. Journal of the Acoustical Society of America, 121, 388–400. doi:10.1121/1.2385043
Duineveld, C.A.A., Arents, P., & King, B.M. (2000). Log-linear modelling of paired comparison data from consumer tests. Food Quality and Preference, 11, 63–70. doi:10.1016/s0950-3293(99)00040-3
## Bradley-Terry-Luce model data(pork) # Is there a difference between Judge 1 and Judge 2? groups <- simplify2array(list(apply(pork[, , 1:5], 1:2, sum), apply(pork[, , 6:10], 1:2, sum))) group.test(groups) # Yes, there is. ## Elimination-by-aspects model data(drugrisk) # Do younger and older males judge risk of drugs differently? A2 <- list(c(1), c(2,7), c(3,7), c(4,7,8), c(5,7,8), c(6,7,8)) group.test(drugrisk[, , 3:4], A2) # Yes.## Bradley-Terry-Luce model data(pork) # Is there a difference between Judge 1 and Judge 2? groups <- simplify2array(list(apply(pork[, , 1:5], 1:2, sum), apply(pork[, , 6:10], 1:2, sum))) group.test(groups) # Yes, there is. ## Elimination-by-aspects model data(drugrisk) # Do younger and older males judge risk of drugs differently? A2 <- list(c(1), c(2,7), c(3,7), c(4,7,8), c(5,7,8), c(6,7,8)) group.test(drugrisk[, , 3:4], A2) # Yes.
Beaver and Gokhale (1975) presented fifty subjects with all 20 ordered pairs of bottles filled with lead shot and asked them to choose the bottle that felt heavier. The mass of the bottles was 90, 95, 100, 105, and 110 grams, respectively. Choice frequencies were aggregated across subjects for the two within-pair presentation orders.
data(heaviness)data(heaviness)
A 3d array consisting of two square matrices:
heaviness[, , order = 1]holds the choices where the row stimulus was presented first for each pair (in the upper triangle, and vice versa in the lower triangle).
heaviness[, , order = 2]holds the choices where the column stimulus was presented first for each pair (in the upper triangle, and vice versa in the lower triangle).
Beaver, R.J., & Gokhale, D.V. (1975). A model to incorporate within-pair order effects in paired comparisons. Communications in Statistics, 4, 923–939. doi:10.1080/03610927308827302
eba.order for a model that includes a within-pair
order effect.
data(heaviness) ## 6 subjects chose 90g over 100g, when 90g was presented first. heaviness["90g", "100g", order=1] ## 44 subjects chose 100g over 90g, when 90g was presented first. heaviness["100g", "90g", order=1] ## 14 subjects chose 90g over 100g, when 90g was presented second. heaviness["90g", "100g", order=2] ## 36 subjects chose 100g over 90g, when 90g was presented second. heaviness["100g", "90g", order=2] ## Bradley-Terry-Luce (BTL) model for each within-pair order btl1 <- eba(heaviness[, , 1]) btl2 <- eba(heaviness[, , 2]) xval <- seq(90, 110, 5) plot(uscale(btl1) ~ xval, type="n", log="y", xlab="Mass of lead shot filled bottle (in g)", ylab="Perceived heaviness (BTL model)", main="Weights judging experiment") mtext("(Beaver and Gokhale, 1975)", line=.5) arrows(xval, uscale(btl1) - 1.96*sqrt(diag(cov.u(btl1))), xval, uscale(btl1) + 1.96*sqrt(diag(cov.u(btl1))), .05, 90, 3, "gray") arrows(xval, uscale(btl2) - 1.96*sqrt(diag(cov.u(btl2))), xval, uscale(btl2) + 1.96*sqrt(diag(cov.u(btl2))), .05, 90, 3, "gray") points(uscale(btl1) ~ xval, type="b", pch=16) points(uscale(btl2) ~ xval, type="b", lty=2, pch=21, bg="white") text(90.3, .01, "Lower weight judged first", adj=0) text(90.3, .08, "Higher weight judged first", adj=0)data(heaviness) ## 6 subjects chose 90g over 100g, when 90g was presented first. heaviness["90g", "100g", order=1] ## 44 subjects chose 100g over 90g, when 90g was presented first. heaviness["100g", "90g", order=1] ## 14 subjects chose 90g over 100g, when 90g was presented second. heaviness["90g", "100g", order=2] ## 36 subjects chose 100g over 90g, when 90g was presented second. heaviness["100g", "90g", order=2] ## Bradley-Terry-Luce (BTL) model for each within-pair order btl1 <- eba(heaviness[, , 1]) btl2 <- eba(heaviness[, , 2]) xval <- seq(90, 110, 5) plot(uscale(btl1) ~ xval, type="n", log="y", xlab="Mass of lead shot filled bottle (in g)", ylab="Perceived heaviness (BTL model)", main="Weights judging experiment") mtext("(Beaver and Gokhale, 1975)", line=.5) arrows(xval, uscale(btl1) - 1.96*sqrt(diag(cov.u(btl1))), xval, uscale(btl1) + 1.96*sqrt(diag(cov.u(btl1))), .05, 90, 3, "gray") arrows(xval, uscale(btl2) - 1.96*sqrt(diag(cov.u(btl2))), xval, uscale(btl2) + 1.96*sqrt(diag(cov.u(btl2))), .05, 90, 3, "gray") points(uscale(btl1) ~ xval, type="b", pch=16) points(uscale(btl2) ~ xval, type="b", lty=2, pch=21, bg="white") text(90.3, .01, "Lower weight judged first", adj=0) text(90.3, .08, "Higher weight judged first", adj=0)
Checks if a family of sets fulfills the inclusion rule.
inclusion.rule(A)inclusion.rule(A)
A |
a list of vectors consisting of the stimulus aspects of an elimination-by-aspects model |
The inclusion rule is necessary and sufficient for a tree structure on the aspect sets:
Structure theorem. A family of aspect sets is
representable by a tree iff either or
for all in .
(Tversky and Sattath, 1979, p. 546)
Either TRUE if the inclusion rule holds for A, or FALSE
otherwise.
Tversky, A., & Sattath, S. (1979). Preference trees. Psychological Review, 86, 542–573. doi:10.1037/0033-295X.86.6.542
A <- list(c(1, 5), c(2, 5), c(3, 6), c(4, 6)) # tree inclusion.rule(A) B <- list(c(1, 5), c(2, 5, 6), c(3, 6), c(4, 6)) # lattice inclusion.rule(B)A <- list(c(1, 5), c(2, 5), c(3, 6), c(4, 6)) # tree inclusion.rule(A) B <- list(c(1, 5), c(2, 5, 6), c(3, 6), c(4, 6)) # lattice inclusion.rule(B)
Kendall's u coefficient of agreement between judges.
kendall.u(M, correct = TRUE)kendall.u(M, correct = TRUE)
M |
a square matrix or a data frame consisting of absolute choice frequencies; row stimuli are chosen over column stimuli |
correct |
logical, if |
Kendall's u (Kendall and Babington Smith, 1940) takes on values between
min.u (minimum agreement) and 1 (maximum agreement).
The minimum min.u equals , if is even,
and , if is odd, where is the number of subjects
(judges).
The null hypothesis in the chi-square test is that the agreement between judges is by chance.
It is assumed that there is an equal number of observations per pair and that each subject judges each pair only once.
u |
Kendall's u coefficient of agreement |
min.u |
the minimum value for u |
chi2 |
the chi-square statistic for a test that the agreement is by chance |
df |
the degrees of freedom |
p.value |
the p-value of the test |
Kendall, M.G., & Babington Smith, B. (1940). On the method of paired comparisons. Biometrika, 31, 324–345. doi:10.1093/biomet/31.3-4.324
schoolsubjects, eba, strans,
circular.
data(schoolsubjects) lapply(schoolsubjects, kendall.u) # better-than-chance agreementdata(schoolsubjects) lapply(schoolsubjects, kendall.u) # better-than-chance agreement
Transforms linear model coefficients to Bradley-Terry-Luce (BTL) model parameter estimates.
linear2btl(object, order = FALSE)linear2btl(object, order = FALSE)
object |
an object of class |
order |
logical, does the model include an order effect? Defaults to FALSE |
The design matrix used by glm or lm usually results from
a call to pcX. It is assumed that the reference category is
the first level. The covariance matrix is estimated by employing the delta
method. See Imrey, Johnson, and Koch (1976) for more details.
btl.parameters |
a matrix; the first column holds the BTL parameter estimates, the second column the approximate standard errors |
cova |
the approximate covariance matrix of the BTL parameter estimates |
linear.coefs |
a vector of the original linear coefficients as returned
by |
Imrey, P.B., Johnson, W.D., & Koch, G.G. (1976). An incomplete contingency table approach to paired-comparison experiments. Journal of the American Statistical Association, 71, 614–623. doi:10.2307/2285591
data(drugrisk) y1 <- t(drugrisk[, , 1])[lower.tri(drugrisk[, , 1])] y0 <- drugrisk[, , 1][ lower.tri(drugrisk[, , 1])] ## Fit BTL model using glm (maximum likelihood) btl.glm <- glm(cbind(y1, y0) ~ 0 + pcX(6), binomial) linear2btl(btl.glm) ## Fit BTL model using lm (weighted least squares) btl.lm <- lm(log(y1/y0) ~ 0 + pcX(6), weights=y1*y0/(y1 + y0)) linear2btl(btl.lm)data(drugrisk) y1 <- t(drugrisk[, , 1])[lower.tri(drugrisk[, , 1])] y0 <- drugrisk[, , 1][ lower.tri(drugrisk[, , 1])] ## Fit BTL model using glm (maximum likelihood) btl.glm <- glm(cbind(y1, y0) ~ 0 + pcX(6), binomial) linear2btl(btl.glm) ## Fit BTL model using lm (weighted least squares) btl.lm <- lm(log(y1/y0) ~ 0 + pcX(6), weights=y1*y0/(y1 + y0)) linear2btl(btl.lm)
Returns the log-likelihood value of the (multi-attribute) probabilistic
choice model represented by object evaluated at the estimated
parameters.
## S3 method for class 'eba' logLik(object, ...)## S3 method for class 'eba' logLik(object, ...)
object |
an object inheriting from class |
... |
some methods for this generic require additional arguments; none are used in this method. |
The log-likelihood of the model represented by
object evaluated at the estimated parameters.
data(heaviness) btl1 <- eba(heaviness[, , order=1]) logLik(btl1) AIC(btl1) BIC(btl1)data(heaviness) btl1 <- eba(heaviness[, , order=1]) logLik(btl1) AIC(btl1) BIC(btl1)
Fits a Mallows-Bradley-Terry (MBT) model by maximum likelihood.
mbt(data, bootstrap = FALSE, nsim = 1000, ...)mbt(data, bootstrap = FALSE, nsim = 1000, ...)
data |
a data frame, the first t columns containing the ranks, the (t + 1)th column containing the frequencies |
bootstrap |
logical. Return a parametric bootstrap p-value? |
nsim |
number of bootstrap replicates |
... |
further aguments passed to |
mbt provides a front end for glm. See Critchlow and Fligner
(1991) and Mallows (1957) for details.
coefficients |
a vector of parameter estimates (scale values) constrained to sum to unity |
goodness.of.fit |
the goodness of fit statistic including the
likelihood ratio fitted vs. saturated model (-2logL), the degrees of
freedom, the p-value of the corresponding chi-square distribution, and
if |
perm.idx |
the names of the non-zero frequency ranks |
y |
the vector of rank frequencies including zeros |
mbt.glm |
the output from a call to |
Florian Wickelmaier
Critchlow, D.E., & Fligner, M.A. (1991). Paired comparison, triple comparison, and ranking experiments as generalized linear models, and their implementation in GLIM. Psychometrika, 56, 517–533. doi:10.1007/bf02294488
Mallows, C.L. (1957). Non-null ranking models. I. Biometrika, 44, 114–130. doi:10.1093/biomet/44.1-2.114
data(tartness) # tartness rankings of salad dressings (Vargo, 1989) mbt(tartness, bootstrap=TRUE, nsim=500) # fit Mallows-Bradley-Terry modeldata(tartness) # tartness rankings of salad dressings (Vargo, 1989) mbt(tartness, bootstrap=TRUE, nsim=500) # fit Mallows-Bradley-Terry model
Computes a paired-comparison design matrix.
pcX(nstimuli, omitRef = TRUE)pcX(nstimuli, omitRef = TRUE)
nstimuli |
number of stimuli in the paired-comparison design |
omitRef |
logical, if |
The design matrix can be used when fitting a Bradley-Terry-Luce (BTL)
model or a Thurstone-Mosteller (TM) model by means of glm
or lm. See Critchlow and Fligner (1991) for more details.
A matrix having (nstimuli - 1)*nstimuli/2 rows and
nstimuli - 1 columns (if the reference category is omitted).
Critchlow, D.E., & Fligner, M.A. (1991). Paired comparison, triple comparison, and ranking experiments as generalized linear models, and their implementation in GLIM. Psychometrika, 56, 517–533. doi:10.1007/bf02294488
eba, thurstone, glm,
balanced.pcdesign, linear2btl.
data(drugrisk) # absolute choice frequencies btl <- eba(drugrisk[, , 1]) # fit Bradley-Terry-Luce model using eba summary(btl) y1 <- t(drugrisk[, , 1])[lower.tri(drugrisk[, , 1])] y0 <- drugrisk[, , 1][ lower.tri(drugrisk[, , 1])] ## Fit Bradley-Terry-Luce model using glm btl.glm <- glm(cbind(y1, y0) ~ 0 + pcX(6), binomial) summary(btl.glm) ## Fit Thurstone Case V model using glm tm.glm <- glm(cbind(y1, y0) ~ 0 + pcX(6), binomial(probit)) summary(tm.glm)data(drugrisk) # absolute choice frequencies btl <- eba(drugrisk[, , 1]) # fit Bradley-Terry-Luce model using eba summary(btl) y1 <- t(drugrisk[, , 1])[lower.tri(drugrisk[, , 1])] y0 <- drugrisk[, , 1][ lower.tri(drugrisk[, , 1])] ## Fit Bradley-Terry-Luce model using glm btl.glm <- glm(cbind(y1, y0) ~ 0 + pcX(6), binomial) summary(btl.glm) ## Fit Thurstone Case V model using glm tm.glm <- glm(cbind(y1, y0) ~ 0 + pcX(6), binomial(probit)) summary(tm.glm)
Plots elimination-by-aspects (EBA) model residuals against fitted values.
## S3 method for class 'eba' plot(x, xlab = "Predicted choice probabilities", ylab = "Deviance residuals", ...)## S3 method for class 'eba' plot(x, xlab = "Predicted choice probabilities", ylab = "Deviance residuals", ...)
x |
an object of class |
xlab, ylab, ...
|
graphical parameters passed to plot. |
The deviance residuals are plotted against the predicted choice probabilities for the upper triangle of the paired-comparison matrix.
## Compare two choice models data(celebrities) # absolute choice frequencies btl1 <- eba(celebrities) # fit Bradley-Terry-Luce model A <- list(c(1,10), c(2,10), c(3,10), c(4,11), c(5,11), c(6,11), c(7,12), c(8,12), c(9,12)) # the structure of aspects eba1 <- eba(celebrities, A) # fit elimination-by-aspects model anova(btl1, eba1) # model comparison based on likelihoods par(mfrow = 1:2) # residuals versus fitted values plot(btl1, main = "BTL", ylim = c(-4, 4.5)) # BTL doesn't fit well plot(eba1, main = "EBA", ylim = c(-4, 4.5)) # EBA fits better## Compare two choice models data(celebrities) # absolute choice frequencies btl1 <- eba(celebrities) # fit Bradley-Terry-Luce model A <- list(c(1,10), c(2,10), c(3,10), c(4,11), c(5,11), c(6,11), c(7,12), c(8,12), c(9,12)) # the structure of aspects eba1 <- eba(celebrities, A) # fit elimination-by-aspects model anova(btl1, eba1) # model comparison based on likelihoods par(mfrow = 1:2) # residuals versus fitted values plot(btl1, main = "BTL", ylim = c(-4, 4.5)) # BTL doesn't fit well plot(eba1, main = "EBA", ylim = c(-4, 4.5)) # EBA fits better
Bradley and Terry (1952) provide the individual choice matrices of two judges choosing between pairs of three samples of pork meet. The pigs had been fed on either corn (C), corn plus peanut supplement (Cp), or corn plus a large peanut supplement (CP). Each judge does five repetitions.
data(pork)data(pork)
A 3d array consisting of ten square matrices. The first five matrices correspond to the five repetitions of the first judge, the last five to the repetitions of the second judge. Row stimuli are chosen (preferred) over column stimuli.
Bradley, R.A., & Terry, M.E. (1952). Rank analysis of incomplete block designs. I. The method of paired comparisons. Biometrika, 39, 324–345. doi:10.1093/biomet/39.3-4.324
data(pork) apply(pork, 1:2, sum) # aggregate choice frequenciesdata(pork) apply(pork, 1:2, sum) # aggregate choice frequencies
Computes deviance and Pearson residuals for eba objects.
## S3 method for class 'eba' residuals(object, type = c("deviance", "pearson"), ...)## S3 method for class 'eba' residuals(object, type = c("deviance", "pearson"), ...)
object |
an object of class |
type |
the type of residuals which should be returned; the
alternatives are: |
... |
further arguments passed to or from other methods; none are used in this method. |
Residuals are computed from the upper triangle of the paired-comparison
matrix. See residuals.glm for details.
A vector of residuals having as many elements as pairs of stimuli.
data(celebrities) # absolute choice frequencies btl1 <- eba(celebrities) # fit Bradley-Terry-Luce model sum( resid(btl1)^2 ) # deviance G2 deviance(btl1) sum( resid(btl1, "pearson")^2 ) # Pearson X2data(celebrities) # absolute choice frequencies btl1 <- eba(celebrities) # fit Bradley-Terry-Luce model sum( resid(btl1)^2 ) # deviance G2 deviance(btl1) sum( resid(btl1, "pearson")^2 ) # Pearson X2
Two classes of children (ages 11 to 13) were asked to state their preferences with respect to certain school subjects. Each child was given a sheet on which were written the possible pairs of subjects and asked to underline the one preferred in each case (Kendall and Babington Smith, 1940).
data(schoolsubjects)data(schoolsubjects)
A list containing two square matrices of aggregate choice frequencies (row entries are preferred over column entries):
schoolsubjects[["boys"]]holds the frequencies of 21 boys choosing among 13 school subjects: woodwork, gymnastics, art, science, history, geography, arithmetic, religion, English literature, commercial subjects, algebra, English grammar, geometry.
schoolsubjects[["girls"]]holds the frequencies of 25 girls choosing among 11 school subjects: gymnastics, science, art, domestic science, history, arithmetic, geography, English literature, religion, algebra, English grammar.
Kendall, M.G., & Babington Smith, B. (1940). On the method of paired comparisons. Biometrika, 31, 324–345. doi:10.1093/biomet/31.3-4.324
data(schoolsubjects) m <- lapply(schoolsubjects, eba) # Bradley-Terry-Luce (BTL) model par(mfrow = 1:2) dotchart(uscale(m$boys), main = "Boys' preferences", log = "x") dotchart(uscale(m$girls), main = "Girls' preferences", log = "x") mtext("Utility scale value (BTL model)", outer = TRUE, side = 1, line = -2) mtext("(Kendall and Babington Smith, 1940)", outer = TRUE, line = -4)data(schoolsubjects) m <- lapply(schoolsubjects, eba) # Bradley-Terry-Luce (BTL) model par(mfrow = 1:2) dotchart(uscale(m$boys), main = "Boys' preferences", log = "x") dotchart(uscale(m$girls), main = "Girls' preferences", log = "x") mtext("Utility scale value (BTL model)", outer = TRUE, side = 1, line = -2) mtext("(Kendall and Babington Smith, 1940)", outer = TRUE, line = -4)
Simulates responses from the distribution corresponding to a fitted
eba model object.
## S3 method for class 'eba' simulate(object, nsim, seed, pool = TRUE, ...)## S3 method for class 'eba' simulate(object, nsim, seed, pool = TRUE, ...)
object |
an object of class |
nsim |
currently not used |
seed |
currently not used |
pool |
logical, if TRUE (default), pooled responses (summed across respondents) are returned |
... |
further arguments passed to or from other methods; none are used in this method |
Responses are simulated by rbinom with sizes taken from the
original sample and probabilities computed from the model object.
A paired-comparison matrix of (pooled) responses.
data(celebrities) # absolute choice frequencies A <- list(c(1,10), c(2,10), c(3,10), c(4,11), c(5,11), c(6,11), c(7,12), c(8,12), c(9,12)) # the structure of aspects eba1 <- eba(celebrities, A) # fit elimination-by-aspects model ## Parametric bootstrap of goodness-of-fit test LR.stat <- replicate(200, deviance(eba(simulate(eba1), A))) hist(LR.stat, col="lightgray", border="white", freq=FALSE, breaks=20, xlim=c(0, 60), main="Parametric bootstrap") curve(dchisq(x, df=eba1$goodness.of.fit["df"]), add=TRUE) abline(v=deviance(eba1), lty=2)data(celebrities) # absolute choice frequencies A <- list(c(1,10), c(2,10), c(3,10), c(4,11), c(5,11), c(6,11), c(7,12), c(8,12), c(9,12)) # the structure of aspects eba1 <- eba(celebrities, A) # fit elimination-by-aspects model ## Parametric bootstrap of goodness-of-fit test LR.stat <- replicate(200, deviance(eba(simulate(eba1), A))) hist(LR.stat, col="lightgray", border="white", freq=FALSE, breaks=20, xlim=c(0, 60), main="Parametric bootstrap") curve(dchisq(x, df=eba1$goodness.of.fit["df"]), add=TRUE) abline(v=deviance(eba1), lty=2)
Paired comparison judgments of 40 selected listeners with respect to
eight audio reproduction modes and four types of music:
SQpreference includes judgments on overall preference;
SQattributes includes judgments on specific spatial and timbral
auditory attributes;
SQsubjects: includes information about the listeners involved.
data("soundquality")data("soundquality")
SQpreference A data frame containing 783 observations on 6 variables:
factor, listener ID.
factor, listening experiment before or after elicitation and scaling of more specific auditory attributes.
factor, the program material: Beethoven, Rachmaninov, Steely Dan, Sting.
the repetition number.
the experimental session coding the presentation order of the program material.
paired comparison of class paircomp;
preferences for all 28 paired comparisons from 8 audio reproduction modes:
Mono, Phantom Mono, Stereo, Wide-Angle Stereo, 4-channel Matrix,
5-channel Upmix 1, 5-channel Upmix 2, and 5-channel Original.
SQattributes A data frame containing 156 observations on 10 variables:
factor, listener ID.
factor, the program material.
Paired comparison of class
paircomp.
SQsubjects A data frame containing 78 observations on 18 variables:
factor, listener ID
factor, selection status; 40 listeners were selected.
maximum hearing level between 250 and 8000 Hz
stereo-width discrimination threshold
word fluency score
subject age
factor, subject gender
factor, education class
indicators of prior experience
The data were collected within a series of experiments conducted at the Sound Quality Research Unit (SQRU), Department of Acoustics, Aalborg University, Denmark, between September 2004 and March 2005.
The results of scaling listener preference and spatial and timbral auditory attributes are reported in Choisel and Wickelmaier (2007). See examples for replication code. Details about the loudspeaker setup and calibration are given in Choisel and Wickelmaier (2006). The attribute elicitation procedure is described in Wickelmaier and Ellermeier (2007) and in Choisel and Wickelmaier (2006). The selection of listeners for the experiments is described in Wickelmaier and Choisel (2005).
Portions of these data are also available via data("SoundQuality",
package = "psychotools").
One listener (ID 62) dropped out after contributing the first set of preference judgments.
Choisel, S., & Wickelmaier, F. (2006). Extraction of auditory features and elicitation of attributes for the assessment of multichannel reproduced sound. Journal of the Audio Engineering Society, 54(9), 815–826.
Choisel, S., & Wickelmaier, F. (2007). Evaluation of multichannel reproduced sound: Scaling auditory attributes underlying listener preference. Journal of the Acoustical Society of America, 121(1), 388–400. doi:10.1121/1.2385043
Wickelmaier, F., & Choisel, S. (2005). Selecting participants for listening tests of multichannel reproduced sound. Presented at the AES 118th Convention, May 28–31, Barcelona, Spain, convention paper 6483.
Wickelmaier, F., & Ellermeier, W. (2007). Deriving auditory features from triadic comparisons. Perception & Psychophysics, 69(2), 287–297. doi:10.3758/BF03193750
requireNamespace("psychotools") data(soundquality) ######### Replication code for Choisel and Wickelmaier (2007) ###### ### A. Scaling overall preference ## Participants summary(subset(SQsubjects, status == "selected")) ## Transitivity violations aggregate(preference ~ progmat + time, data = SQpreference, function(x) unlist(strans(summary(x, pcmatrix = TRUE))[ c("weak", "moderate", "strong")])) ## BTL modeling prefdf <- aggregate(preference ~ progmat + time, data = SQpreference, function(x) uscale(eba(summary(x, pcmatrix = TRUE)))) ## Preference scale p <- t(prefdf[prefdf$time == "before", 3]) colnames(p) <- levels(prefdf$progmat) dotchart(p, main = "Quality of multichannel reproduced sound", xlab = "Estimated preference (BTL model)", log = "x", panel.first = abline(v = 1/8, col = "gray")) points(x = t(prefdf[prefdf$time == "after", 3]), y = c(31:38, 21:28, 11:18, 1:8), pch = 3) legend("topleft", c("before", "after"), pch = c(1, 3)) mtext("(Choisel and Wickelmaier, 2007)", line = 0.5) ### B. Scaling specific auditory attributes ## Transitivity violations aggregate( x = SQattributes[-(1:2)], by = list(progmat = SQattributes$progmat), FUN = function(x) strans(summary(x, pcmatrix = TRUE))[ c("weak", "moderate", "strong")], simplify = FALSE ) ## BTL modeling uscal <- aggregate( x = SQattributes[-(1:2)], by = list(progmat = SQattributes$progmat), FUN = function(x) uscale(eba(summary(x, pcmatrix = TRUE))) ) uscal <- data.frame( progmat = rep(levels(SQattributes$progmat), each = 8), repmod = factor(1:8, labels = labels(SQattributes$width)), as.data.frame(sapply(uscal[-1], t)) ) ## EBA modeling: envelopment, width uscal$envelopment[uscal$progmat == "SteelyDan"] <- uscale(eba(summary(SQattributes[SQattributes$progmat == "SteelyDan", "envelopment"], pcmatrix = TRUE), A = list(c(1,9), c(2,9), c(3,9), c(4,9), 5, 6, c(7,9), 8))) uscal$width[uscal$progmat == "SteelyDan"] <- uscale(eba(summary(SQattributes[SQattributes$progmat == "SteelyDan", "width"], pcmatrix = TRUE), A = list(1, 2, c(3,9), c(4,9), c(5,9), 6, 7, c(8,9)))) ### C. Relating overall preference to specific attributes ## Principal components: classical music pca1 <- princomp( ~ width + spaciousness + envelopment + distance + clarity + brightness + elevation, data = uscal, subset = progmat %in% c("Beethoven", "Rachmaninov"), cor = TRUE ) ## Loadings on first two components L <- varimax(loadings(pca1) %*% diag(pca1$sdev)[, 1:2]) ## Factor scores f <- scale(predict(pca1)[, 1:2]) %*% L$rotmat dimnames(f) <- list( abbreviate(rep(labels(SQattributes$width), 2), 2), names(pca1$sdev)[1:2] ) biplot(f, 2*L$loadings, cex = 0.8, expand = 1.5, main = "Preference and auditory attributes: classical music", xlim = c(-1.5, 2), ylim = c(-2, 2)) ## Predicting preference classdf <- cbind( pref = as.vector(t(prefdf[prefdf$time == "after" & prefdf$progmat %in% c("Beethoven", "Rachmaninov"), 3])), as.data.frame(f) ) m1 <- lm(pref ~ Comp.1 + Comp.2 + I(Comp.1^2), classdf) c1 <- seq(-1.5, 2.0, length.out = 101) c2 <- seq(-2.0, 2.0, length.out = 101) z <- matrix(predict(m1, expand.grid(Comp.1 = c1, Comp.2 = c2)), 101, 101) contour(c1, c2, z, add = TRUE, col = "darkblue") ## Principal components: pop music pca2 <- princomp( ~ width + spaciousness + envelopment + distance + clarity + brightness + elevation, data = uscal, subset = progmat %in% c("SteelyDan", "Sting"), cor = TRUE ) L <- varimax(loadings(pca2) %*% diag(pca2$sdev)[, 1:2]) f[] <- scale(predict(pca2)[, 1:2]) %*% L$rotmat biplot(f, 2*L$loadings, cex = 0.8, main = "Preference and auditory attributes: pop music", xlim = c(-1.5, 2), ylim = c(-2.5, 1.5)) popdf <- cbind( pref = as.vector(t(prefdf[prefdf$time == "after" & prefdf$progmat %in% c("SteelyDan", "Sting"), 3])), as.data.frame(f) ) m2 <- lm(pref ~ Comp.1 + Comp.2 + I(Comp.2^2), popdf) c1 <- seq(-1.5, 2.0, length.out = 101) c2 <- seq(-2.5, 1.5, length.out = 101) z <- matrix(predict(m2, expand.grid(Comp.1 = c1, Comp.2 = c2)), 101, 101) contour(c1, c2, z, add = TRUE, col = "darkblue")requireNamespace("psychotools") data(soundquality) ######### Replication code for Choisel and Wickelmaier (2007) ###### ### A. Scaling overall preference ## Participants summary(subset(SQsubjects, status == "selected")) ## Transitivity violations aggregate(preference ~ progmat + time, data = SQpreference, function(x) unlist(strans(summary(x, pcmatrix = TRUE))[ c("weak", "moderate", "strong")])) ## BTL modeling prefdf <- aggregate(preference ~ progmat + time, data = SQpreference, function(x) uscale(eba(summary(x, pcmatrix = TRUE)))) ## Preference scale p <- t(prefdf[prefdf$time == "before", 3]) colnames(p) <- levels(prefdf$progmat) dotchart(p, main = "Quality of multichannel reproduced sound", xlab = "Estimated preference (BTL model)", log = "x", panel.first = abline(v = 1/8, col = "gray")) points(x = t(prefdf[prefdf$time == "after", 3]), y = c(31:38, 21:28, 11:18, 1:8), pch = 3) legend("topleft", c("before", "after"), pch = c(1, 3)) mtext("(Choisel and Wickelmaier, 2007)", line = 0.5) ### B. Scaling specific auditory attributes ## Transitivity violations aggregate( x = SQattributes[-(1:2)], by = list(progmat = SQattributes$progmat), FUN = function(x) strans(summary(x, pcmatrix = TRUE))[ c("weak", "moderate", "strong")], simplify = FALSE ) ## BTL modeling uscal <- aggregate( x = SQattributes[-(1:2)], by = list(progmat = SQattributes$progmat), FUN = function(x) uscale(eba(summary(x, pcmatrix = TRUE))) ) uscal <- data.frame( progmat = rep(levels(SQattributes$progmat), each = 8), repmod = factor(1:8, labels = labels(SQattributes$width)), as.data.frame(sapply(uscal[-1], t)) ) ## EBA modeling: envelopment, width uscal$envelopment[uscal$progmat == "SteelyDan"] <- uscale(eba(summary(SQattributes[SQattributes$progmat == "SteelyDan", "envelopment"], pcmatrix = TRUE), A = list(c(1,9), c(2,9), c(3,9), c(4,9), 5, 6, c(7,9), 8))) uscal$width[uscal$progmat == "SteelyDan"] <- uscale(eba(summary(SQattributes[SQattributes$progmat == "SteelyDan", "width"], pcmatrix = TRUE), A = list(1, 2, c(3,9), c(4,9), c(5,9), 6, 7, c(8,9)))) ### C. Relating overall preference to specific attributes ## Principal components: classical music pca1 <- princomp( ~ width + spaciousness + envelopment + distance + clarity + brightness + elevation, data = uscal, subset = progmat %in% c("Beethoven", "Rachmaninov"), cor = TRUE ) ## Loadings on first two components L <- varimax(loadings(pca1) %*% diag(pca1$sdev)[, 1:2]) ## Factor scores f <- scale(predict(pca1)[, 1:2]) %*% L$rotmat dimnames(f) <- list( abbreviate(rep(labels(SQattributes$width), 2), 2), names(pca1$sdev)[1:2] ) biplot(f, 2*L$loadings, cex = 0.8, expand = 1.5, main = "Preference and auditory attributes: classical music", xlim = c(-1.5, 2), ylim = c(-2, 2)) ## Predicting preference classdf <- cbind( pref = as.vector(t(prefdf[prefdf$time == "after" & prefdf$progmat %in% c("Beethoven", "Rachmaninov"), 3])), as.data.frame(f) ) m1 <- lm(pref ~ Comp.1 + Comp.2 + I(Comp.1^2), classdf) c1 <- seq(-1.5, 2.0, length.out = 101) c2 <- seq(-2.0, 2.0, length.out = 101) z <- matrix(predict(m1, expand.grid(Comp.1 = c1, Comp.2 = c2)), 101, 101) contour(c1, c2, z, add = TRUE, col = "darkblue") ## Principal components: pop music pca2 <- princomp( ~ width + spaciousness + envelopment + distance + clarity + brightness + elevation, data = uscal, subset = progmat %in% c("SteelyDan", "Sting"), cor = TRUE ) L <- varimax(loadings(pca2) %*% diag(pca2$sdev)[, 1:2]) f[] <- scale(predict(pca2)[, 1:2]) %*% L$rotmat biplot(f, 2*L$loadings, cex = 0.8, main = "Preference and auditory attributes: pop music", xlim = c(-1.5, 2), ylim = c(-2.5, 1.5)) popdf <- cbind( pref = as.vector(t(prefdf[prefdf$time == "after" & prefdf$progmat %in% c("SteelyDan", "Sting"), 3])), as.data.frame(f) ) m2 <- lm(pref ~ Comp.1 + Comp.2 + I(Comp.2^2), popdf) c1 <- seq(-1.5, 2.0, length.out = 101) c2 <- seq(-2.5, 1.5, length.out = 101) z <- matrix(predict(m2, expand.grid(Comp.1 = c1, Comp.2 = c2)), 101, 101) contour(c1, c2, z, add = TRUE, col = "darkblue")
Checks the weak, moderate, and strong stochastic transitivity.
strans(M)strans(M)
M |
a square matrix or a data frame consisting of absolute choice frequencies; row stimuli are chosen over column stimuli |
Weak (WST), moderate (MST), and strong (SST) stochastic transitivity hold
for a set of choice probabilities , whenever if
and , then
(WST),
(MST),
(SST).
See Suppes, Krantz, Luce, and Tversky (1989/2007, chap. 17) for an introduction to the representation of choice probabilities.
If WST holds, a permutation of the indices of the matrix exists such that
the proportions in the upper triangular matrix are . This
rearranged matrix is stored in pcm. If WST does not hold, cells in
the upper triangular matrix that are smaller than 0.5 are replaced by 0.5.
The deviance resulting from this restriction is reported in wst.fit.
The approximate likelihood ratio test for significance of the WST violations is according to Tversky (1969); for a more exact test of WST see Iverson and Falmagne (1985).
A table displaying the number of violations of the weak, moderate, and strong stochastic transitivity, the number of tests, the error ratio (violations/tests), and the mean and maximum deviation from the minimum probability for which the corresponding transitivity would hold.
weak |
number of violations of WST |
moderate |
number of violations of MST |
strong |
number of violations of SST |
n.tests |
number of transitivity tests performed |
wst.violations |
a vector containing
|
mst.violations |
a vector containing
|
sst.violations |
a vector containing
|
pcm |
the permuted square matrix of relative choice frequencies |
ranking |
the ranking of the objects, which corresponds to the colnames of pcm |
chkdf |
data frame reporting the choice proportions for each triple in each permutation |
violdf |
data frame reporting for each triple which type of transitivity holds or does not hold |
wst.fit |
likelihood ratio test of WST (see Details) |
wst.mat |
restricted matrix that satisfies WST |
Iverson, G., & Falmagne, J.-C. (1985). Statistical issues in measurement. Mathematical Social Sciences, 10, 131–153. doi:10.1016/0165-4896(85)90031-9
Suppes, P., Krantz, D.H., Luce, R.D., & Tversky, A. (1989/2007). Foundations of measurement. Volume II. Mineola, N.Y.: Dover Publications.
Tversky, A. (1969). Intransitivity of preferences. Psychological Review, 76, 31–48. doi:10.1037/h0026750
eba, circular, kendall.u,
trineq.
data(celebrities) # absolute choice frequencies strans(celebrities) # WST and MST hold, but not SST strans(celebrities)$pcm # reordered relative frequencies strans(celebrities)$violdf # transitivity violationsdata(celebrities) # absolute choice frequencies strans(celebrities) # WST and MST hold, but not SST strans(celebrities)$pcm # reordered relative frequencies strans(celebrities)$violdf # transitivity violations
The data were collected by Vargo (1989; cited in Critchlow and Fligner, 1991). Each of 32 judges is asked to rank four salad dressing preparations according to tartness, with a rank of 1 being assigned to the formulation judged to be the most tart.
data(tartness)data(tartness)
A data frame consisting the rankings and their frequencies.
Critchlow, D.E., & Fligner, M.A. (1991). Paired comparison, triple comparison, and ranking experiments as generalized linear models, and their implementation in GLIM. Psychometrika, 56, 517–533. doi:10.1007/bf02294488
Vargo, M.D. (1989). Microbiological spoilage of a moderate acid food system using a dairy-based salad dressing model. Unpublished masters thesis, Ohio State University, Department of Food Science and Nutrition, Columbus, OH.
data(tartness)data(tartness)
Fits a Thurstone-Mosteller model (Case V) by maximum likelihood.
thurstone(M)thurstone(M)
M |
a square matrix or a data frame consisting of absolute choice frequencies; row stimuli are chosen over column stimuli |
thurstone provides a front end for glm. See Critchlow and
Fligner (1991) for more details.
estimate |
a vector of parameter estimates (scale values), first element is set to zero |
goodness.of.fit |
the goodness of fit statistic including the likelihood ratio fitted vs. saturated model (-2logL), the degrees of freedom, and the p-value of the corresponding chi-square distribution |
tm.glm |
the output from a call to |
Critchlow, D.E., & Fligner, M.A. (1991). Paired comparison, triple comparison, and ranking experiments as generalized linear models, and their implementation in GLIM. Psychometrika, 56, 517–533. doi:10.1007/bf02294488
eba, strans, pcX,
kendall.u, circular, glm.
## Taste data (David, 1988, p. 116) taste <- matrix(c( 0, 3, 2, 2, 12, 0, 11, 3, 13, 4, 0, 5, 13, 12, 10, 0), 4, 4, byrow=TRUE) dimnames(taste) <- setNames(rep(list(c("A1", "A2", "A3", "A4")), 2), c(">", "<")) thurstone(taste) # Thurstone-Mosteller model fits OK## Taste data (David, 1988, p. 116) taste <- matrix(c( 0, 3, 2, 2, 12, 0, 11, 3, 13, 4, 0, 5, 13, 12, 10, 0), 4, 4, byrow=TRUE) dimnames(taste) <- setNames(rep(list(c("A1", "A2", "A3", "A4")), 2), c(">", "<")) thurstone(taste) # Thurstone-Mosteller model fits OK
Checks if binary choice probabilities fulfill the trinary inequality.
trineq(M, A = 1:I)trineq(M, A = 1:I)
M |
a square matrix or a data frame consisting of absolute choice frequencies; row stimuli are chosen over column stimuli |
A |
a list of vectors consisting of the stimulus aspects;
the default is |
For any triple of stimuli , the trinary inequality states
that, if and , then
where ,
, and denotes that and
share at least one aspect that does not have
(Tversky and Sattath, 1979, p. 554).
inclusion.rule checks if a family of aspect sets is
representable by a tree.
Results checking the trinary inequality.
n |
number of tests of the trinary inequality |
prop |
proportion of triples confirming the trinary inequality |
quant |
quantiles of |
n.tests |
number of transitivity tests performed |
chkdf |
data frame reporting |
Tversky, A., & Sattath, S. (1979). Preference trees. Psychological Review, 86, 542–573. doi:10.1037/0033-295X.86.6.542
data(celebrities) # absolute choice frequencies A <- list(c(1,10), c(2,10), c(3,10), c(4,11), c(5,11), c(6,11), c(7,12), c(8,12), c(9,12)) # the structure of aspects trineq(celebrities, A) # check trinary inequality for tree A trineq(celebrities, A)$chkdf # trinary inequality for each tripledata(celebrities) # absolute choice frequencies A <- list(c(1,10), c(2,10), c(3,10), c(4,11), c(5,11), c(6,11), c(7,12), c(8,12), c(9,12)) # the structure of aspects trineq(celebrities, A) # check trinary inequality for tree A trineq(celebrities, A)$chkdf # trinary inequality for each triple
Extract the (normalized) utility scale for an elimination-by-aspects (EBA) model.
uscale(object, norm = "sum", log = FALSE)uscale(object, norm = "sum", log = FALSE)
object |
an object of class |
norm |
either |
log |
should the log of the utility scale values be returned? Defaults to FALSE. |
Each utility scale value is defined as the sum of aspect values (EBA model parameters) that characterize a given stimulus. First these sums are computed for all stimuli, then normalization (if any) is applied. As each type of normalization corresponds to a multiplication by a positive real number, the ratio between scale values remains constant.
The (normalized) utility scale of the stimuli.
data(drugrisk) A <- list(c(1), c(2,7), c(3,7), c(4,7,8), c(5,7,8), c(6,7,8)) eba1 <- eba(drugrisk[, , group = "male30"], A) # EBA model uscale(eba1) # sum-to-unity normalization uscale(eba1, norm=1) # u(alcohol) := 1 uscale(eba1, norm=5) # u(heroine) := 1 uscale(eba1, norm=NULL) # no normalization uscale(eba1, norm=1, log=TRUE) # log utility scale, log u(alcohol) := 0data(drugrisk) A <- list(c(1), c(2,7), c(3,7), c(4,7,8), c(5,7,8), c(6,7,8)) eba1 <- eba(drugrisk[, , group = "male30"], A) # EBA model uscale(eba1) # sum-to-unity normalization uscale(eba1, norm=1) # u(alcohol) := 1 uscale(eba1, norm=5) # u(heroine) := 1 uscale(eba1, norm=NULL) # no normalization uscale(eba1, norm=1, log=TRUE) # log utility scale, log u(alcohol) := 0
Tests linear hypotheses of the form in elimination-by-aspects
(EBA) models using the Wald test.
wald.test(object, C, u.scale = TRUE)wald.test(object, C, u.scale = TRUE)
object |
an object of class |
C |
a matrix of contrasts, specifying the linear hypotheses |
u.scale |
logical, if TRUE the test is performed on the utility scale, if FALSE the test is performed on the EBA parameters directly |
The Wald test statistic,
is approximately chi-square distributed with degrees of
freedom.
C is usually of full rank and must have as many columns as there
are parameters in p.
C |
the matrix of contrasts, specifying the linear hypotheses |
W |
the Wald test statistic |
df |
the degrees of freedom ( |
pval |
the p-value of the test |
eba, group.test, uscale,
cov.u.
data(celebrities) # absolute choice frequencies A <- list(c(1,10), c(2,10), c(3,10), c(4,11), c(5,11), c(6,11), c(7,12), c(8,12), c(9,12)) # the structure of aspects eba1 <- eba(celebrities, A) # fit elimination-by-aspects model ## Test whether JU, CY, and AJF have equal utility scale values C1 <- rbind(c(0,0,0,1,-1, 0,0,0,0), c(0,0,0,1, 0,-1,0,0,0)) wald.test(eba1, C1) ## Test whether the three branch parameters are different C2 <- rbind(c(0,0,0,0,0,0,0,0,0,1,-1, 0), c(0,0,0,0,0,0,0,0,0,1, 0,-1)) wald.test(eba1, C2, u.scale = FALSE)data(celebrities) # absolute choice frequencies A <- list(c(1,10), c(2,10), c(3,10), c(4,11), c(5,11), c(6,11), c(7,12), c(8,12), c(9,12)) # the structure of aspects eba1 <- eba(celebrities, A) # fit elimination-by-aspects model ## Test whether JU, CY, and AJF have equal utility scale values C1 <- rbind(c(0,0,0,1,-1, 0,0,0,0), c(0,0,0,1, 0,-1,0,0,0)) wald.test(eba1, C1) ## Test whether the three branch parameters are different C2 <- rbind(c(0,0,0,0,0,0,0,0,0,1,-1, 0), c(0,0,0,0,0,0,0,0,0,1, 0,-1)) wald.test(eba1, C2, u.scale = FALSE)
Paired comparison judgments for two wine tasting studies:
ambilight includes the results of a study on the effect of ambient
lighting on the flavor of wine; redwines includes judgments on the
sensory quality of red wines.
data("winetaste")data("winetaste")
ambilight A data frame containing 230 observations on 10 variables:
Paired comparison of
class paircomp; judgments for one of the 6 ordered pairs
of the blue, red, and white lighting conditions.
subject age
factor, subject gender
self-rating of sense of smell and taste.
self-rating of general liking of wine.
factor, frequency of wine consumption.
factor, preference for red or white wine.
redwines A data frame containing 11 observations on 7 variables:
Paired
comparison of class paircomp; judgments for all 10 pairs
from 5 red wines: Primitivo di Manduria, Cotes du Rhone, Bourgogne,
Shiraz cuvee, Cabernet Sauvignon.
factor, Which of the five wines did you like best?
factor, Which of the five wines did you like worst?
The ambilight data are from Experiment 3 in Oberfeld et al. (2009).
The redwines data were collected among the members of the Sound
Quality Research Unit (SQRU), Department of Acoustics, Aalborg University,
Denmark, in 2004. Details of the red wines are available as an attribute of
the preference variable (see Examples).
Oberfeld, D., Hecht, H., Allendorf, U., & Wickelmaier, F. (2009). Ambient lighting modifies the flavor of wine. Journal of Sensory Studies, 24(6), 797–832. doi:10.1111/j.1745-459X.2009.00239.x
requireNamespace("psychotools") data(winetaste) ## No effect of ambient lighting on flavor (Oberfeld et al., 2009) m <- lapply(ambilight[, c("preference", "fruitiness", "spiciness", "sweetness")], function(x) eba.order(summary(x, pcmatrix = TRUE))) lapply(m, summary) u <- sapply(m, uscale, norm = 3) dotchart( u, xlim = c(0.5, 2), pch = 16, panel.first = abline(v = 1, col = "gray"), main = "Ambient lighting and the flavor of wine", xlab = "Utility scale value (Davidson-Beaver model)" ) ci <- sapply(m, function(x) 1.96 * sqrt(diag(cov.u(x)))) arrows( u - ci, c(16:18, 11:13, 6:8, 1:3), u + ci, c(16:18, 11:13, 6:8, 1:3), .05, 90, 3) mtext("Oberfeld et al. (2009)", line = 0.5) ## Sensory quality of red wines psychotools::covariates(redwines$preference) # details of the wines m <- lapply(redwines[, c("bitterness", "fruitiness", "sourness", "roundness", "preference")], function(x) eba(summary(x, pcmatrix = TRUE))) lapply(m, summary) u <- sapply(m, uscale) dotchart( u[order(u[, "preference"]), ], log = "x", panel.first = abline(v = 1/5, col = "gray"), main = "SQRU red wine tasting", xlab = "Utility scale value (BTL model), choice proportion (+)" ) points(as.vector( prop.table(table(redwines$best))[order(u[, "preference"])] ), 1:5, pch = 3)requireNamespace("psychotools") data(winetaste) ## No effect of ambient lighting on flavor (Oberfeld et al., 2009) m <- lapply(ambilight[, c("preference", "fruitiness", "spiciness", "sweetness")], function(x) eba.order(summary(x, pcmatrix = TRUE))) lapply(m, summary) u <- sapply(m, uscale, norm = 3) dotchart( u, xlim = c(0.5, 2), pch = 16, panel.first = abline(v = 1, col = "gray"), main = "Ambient lighting and the flavor of wine", xlab = "Utility scale value (Davidson-Beaver model)" ) ci <- sapply(m, function(x) 1.96 * sqrt(diag(cov.u(x)))) arrows( u - ci, c(16:18, 11:13, 6:8, 1:3), u + ci, c(16:18, 11:13, 6:8, 1:3), .05, 90, 3) mtext("Oberfeld et al. (2009)", line = 0.5) ## Sensory quality of red wines psychotools::covariates(redwines$preference) # details of the wines m <- lapply(redwines[, c("bitterness", "fruitiness", "sourness", "roundness", "preference")], function(x) eba(summary(x, pcmatrix = TRUE))) lapply(m, summary) u <- sapply(m, uscale) dotchart( u[order(u[, "preference"]), ], log = "x", panel.first = abline(v = 1/5, col = "gray"), main = "SQRU red wine tasting", xlab = "Utility scale value (BTL model), choice proportion (+)" ) points(as.vector( prop.table(table(redwines$best))[order(u[, "preference"])] ), 1:5, pch = 3)