Balanced Paired-Comparison Design

Description

Creates a (completely) balanced paired-comparison design.

Usage

balanced.pcdesign(nstimuli)

Arguments

Details

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.

Value

References

David, H. (1988). The method of paired comparisons. London: Griffin.

See Also

pcX, eba.

Examples

## 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)))

Bootstrap for Elimination-by-Aspects (EBA) Models

Description

Performs a bootstrap by resampling the individual data matrices.

Usage

boot(D, R = 100, A = 1:I, s = rep(1/J, J), constrained = TRUE)

Arguments

Details

The bootstrap functions eba.boot.constrained and eba.boot are automatically called by boot.

The code is experimental and may change in the future.

Value

See Also

eba, simulate.eba, paircomp.

Examples

data(pork)  # pork tasting data, 10 individual paired comparison matrices
m <- eba(apply(pork, 1:2, sum))  # fit Bradley-Terry-Luce model
b <- boot(pork, R = 200)         # resample 200 times

plot(coef(m), b$stat[, "mean"], log = "xy")
abline(0, 1, lty = 2)

Choice among Celebrities

Description

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.

Usage

data(celebrities)

Format

A square data frame containing the absolute choice frequencies and a diagonal of zeros; row stimuli are chosen over column stimuli.

Source

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

Examples

data(celebrities)
celebrities["LBJ", "HW"]  # 159 participants chose Johnson over Wilson

Circular Triads (Intransitive Cycles)

Description

Number of circular triads and coefficient of consistency.

Usage

circular(mat, alternative = c("two.sided", "less", "greater"),
         exact = NULL, correct = TRUE, simulate.p.value = FALSE,
         nsim = 2000)

Arguments

Details

Kendall's coefficient of consistency,

zeta = 1 - T/T_{max},

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 T_{max} = n(n^2 - 4)/24, if n is even, and T_{max} = n(n^2 - 1)/24 else, and n 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 n, 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 n < 8 and is only available for n > 4.

Value

References

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

See Also

eba, strans, kendall.u.

Examples

# 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

Covariance Matrix of the EBA Utility Scale

Description

Computes the (normalized) covariance matrix of the utility scale from the covariance matrix of elimination-by-aspects (EBA) model parameters.

Usage

cov.u(object, norm = "sum")

Arguments

Details

The additivity rule for covariances cov(x + y, z) = cov(x, z) + cov(y, z) is used for the computations.

If norm is not NULL, the unnormalized covariance matrix is transformed using a^2 cov(u), where the constant a results from the type of normalization applied.

Value

The (normalized) covariance matrix of the utility scale.

See Also

uscale, eba, wald.test.

Perceived Health Risk of Drugs

Description

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.

Usage

data(drugrisk)

Format

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.

Source

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.

Examples

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)

Elimination-by-Aspects (EBA) Models

Description

Fits a (multi-attribute) probabilistic choice model by maximum likelihood.

Usage

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"))

Arguments

Details

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.

Value

Author(s)

Florian Wickelmaier

References

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

See Also

strans, uscale, cov.u, group.test, wald.test, plot.eba, residuals.eba, logLik.eba, simulate.eba, kendall.u, circular, trineq, thurstone, nlm.

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.

Elimination-by-Aspects (EBA) Models with Order-Effect

Description

Fits a (multi-attribute) probabilistic choice model that accounts for the effect of the presentation order within a pair.

Usage

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, ...)

Arguments

Details

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.

Value

Author(s)

Florian Wickelmaier

References

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.

See Also

eba, group.test, plot.eba, residuals.eba, logLik.eba.

Examples

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 parameters

Auditory Unpleasantness of Environmental Sound

Description

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.

Usage

data(envirosound)

Format

A data frame containing 74 observations on 2 variables:

unpleasantness

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.

rt

median response time.

Details

Details of the recordings, including psychoacoustic metrics, are available as an attribute of the unpleasantness variable (see Examples).

Source

Zimmer, K., Ellermeier, W., & Schmid, C. (2004). Using probabilistic choice models to investigate auditory unpleasantness. Acta Acustica united with Acustica, 90(6), 1019–1028.

References

Johannsen, K., & Prante, H.U. (2001). Environmental sounds for psychoacoustic testing. Acta Acustica united with Acustica, 87(2), 290–293.

See Also

eba, strans, paircomp.

Examples

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
)

Group Effects in Elimination-by-Aspects (EBA) Models

Description

Tests for group effects in elimination-by-aspects (EBA) models.

Usage

group.test(groups, A = 1:I, s = rep(1/J, J), constrained = TRUE)

Arguments

Details

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.

Value

References

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

See Also

eba, wald.test.

Examples

## 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.

Weights Judging Data

Description

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.

Usage

data(heaviness)

Format

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).

Source

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

See Also

eba.order for a model that includes a within-pair order effect.

Examples

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)

Inclusion Rule

Description

Checks if a family of sets fulfills the inclusion rule.

Usage

inclusion.rule(A)

Arguments

Details

The inclusion rule is necessary and sufficient for a tree structure on the aspect sets:

Structure theorem. A family \{x' | x \in T\} of aspect sets is representable by a tree iff either x' \cap y' \supset x' \cap z' or x' \cap z' \supset x' \cap y' for all x, y, z in T. (Tversky and Sattath, 1979, p. 546)

Value

Either TRUE if the inclusion rule holds for A, or FALSE otherwise.

References

Tversky, A., & Sattath, S. (1979). Preference trees. Psychological Review, 86, 542–573. doi:10.1037/0033-295X.86.6.542

See Also

eba, trineq, strans.

Examples

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 Coefficient of Agreement

Description

Kendall's u coefficient of agreement between judges.

Usage

kendall.u(M, correct = TRUE)

Arguments

Details

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 -1/(m - 1), if m is even, and -1/m, if m is odd, where m 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.

Value

References

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

See Also

schoolsubjects, eba, strans, circular.

Examples

data(schoolsubjects)
lapply(schoolsubjects, kendall.u)  # better-than-chance agreement

Linear Coefficients to Bradley-Terry-Luce (BTL) Estimates

Description

Transforms linear model coefficients to Bradley-Terry-Luce (BTL) model parameter estimates.

Usage

linear2btl(object, order = FALSE)

Arguments

Details

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.

Value

References

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

See Also

eba, eba.order, glm, pcX.

Examples

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)

Log-Likelihood of an eba Object

Description

Returns the log-likelihood value of the (multi-attribute) probabilistic choice model represented by object evaluated at the estimated parameters.

Usage

## S3 method for class 'eba'
logLik(object, ...)

Arguments

Value

The log-likelihood of the model represented by object evaluated at the estimated parameters.

See Also

eba, logLik.lm, AIC.

Examples

data(heaviness)
btl1 <- eba(heaviness[, , order=1])
logLik(btl1)
AIC(btl1)
BIC(btl1)

Mallows-Bradley-Terry Model

Description

Fits a Mallows-Bradley-Terry (MBT) model by maximum likelihood.

Usage

mbt(data, bootstrap = FALSE, nsim = 1000, ...)

Arguments

Details

mbt provides a front end for glm. See Critchlow and Fligner (1991) and Mallows (1957) for details.

Value

Author(s)

Florian Wickelmaier

References

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

See Also

tartness, glm.

Examples

data(tartness)        # tartness rankings of salad dressings (Vargo, 1989)
mbt(tartness, bootstrap=TRUE, nsim=500)  # fit Mallows-Bradley-Terry model

Paired-Comparison Design Matrix

Description

Computes a paired-comparison design matrix.

Usage

pcX(nstimuli, omitRef = TRUE)

Arguments

Details

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.

Value

A matrix having (nstimuli - 1)*nstimuli/2 rows and nstimuli - 1 columns (if the reference category is omitted).

References

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

See Also

eba, thurstone, glm, balanced.pcdesign, linear2btl.

Examples

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)

Diagnostic Plot for EBA Models

Description

Plots elimination-by-aspects (EBA) model residuals against fitted values.

Usage

## S3 method for class 'eba'
plot(x, xlab = "Predicted choice probabilities",
     ylab = "Deviance residuals", ...)

Arguments

Details

The deviance residuals are plotted against the predicted choice probabilities for the upper triangle of the paired-comparison matrix.

See Also

eba, residuals.eba.

Examples

## 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

Pork Tasting Data

Description

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.

Usage

data(pork)

Format

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.

Source

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

Examples

data(pork)
apply(pork, 1:2, sum)  # aggregate choice frequencies

Residuals for EBA Models

Description

Computes deviance and Pearson residuals for eba objects.

Usage

## S3 method for class 'eba'
residuals(object, type = c("deviance", "pearson"), ...)

Arguments

Details

Residuals are computed from the upper triangle of the paired-comparison matrix. See residuals.glm for details.

Value

A vector of residuals having as many elements as pairs of stimuli.

See Also

eba, residuals.glm, plot.eba.

Examples

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 X2

Preference for School Subjects

Description

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).

Usage

data(schoolsubjects)

Format

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.

Source

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

Examples

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)

Simulate Responses from EBA Models

Description

Simulates responses from the distribution corresponding to a fitted eba model object.

Usage

## S3 method for class 'eba'
simulate(object, nsim, seed, pool = TRUE, ...)

Arguments

Details

Responses are simulated by rbinom with sizes taken from the original sample and probabilities computed from the model object.

Value

A paired-comparison matrix of (pooled) responses.

See Also

eba, boot.

Examples

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)

Quality of Multichannel Reproduced Sound

Description

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.

Usage

data("soundquality")

Format

SQpreference A data frame containing 783 observations on 6 variables:

id

factor, listener ID.

time

factor, listening experiment before or after elicitation and scaling of more specific auditory attributes.

progmat

factor, the program material: Beethoven, Rachmaninov, Steely Dan, Sting.

repet

the repetition number.

session

the experimental session coding the presentation order of the program material.

preference

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:

id

factor, listener ID.

progmat

factor, the program material.

width, spaciousness, envelopment, distance, clarity, brightness, elevation, naturalness

Paired comparison of class paircomp.

SQsubjects A data frame containing 78 observations on 18 variables:

id

factor, listener ID

status

factor, selection status; 40 listeners were selected.

HLmax

maximum hearing level between 250 and 8000 Hz

stereowidth

stereo-width discrimination threshold

fluency

word fluency score

age

subject age

gender

factor, subject gender

education

factor, education class

background, experience, listenmusic, concerts, instrument, critical, cinema, hifi, surround, earliertests

indicators of prior experience

Details

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").

Note

One listener (ID 62) dropped out after contributing the first set of preference judgments.

References

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

See Also

eba, strans, paircomp.

Examples

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")

Stochastic Transitivity

Description

Checks the weak, moderate, and strong stochastic transitivity.

Usage

strans(M)

Arguments

Details

Weak (WST), moderate (MST), and strong (SST) stochastic transitivity hold for a set of choice probabilities P, whenever if P_{ij} \ge 0.5 and P_{jk} \ge 0.5, then

P_{ik} \ge 0.5 (WST),

P_{ik} \ge \min(P_{ij}, P_{jk}) (MST),

P_{ik} \ge \max(P_{ij}, P_{jk}) (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 \ge 0.5. 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).

Value

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.

References

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

See Also

eba, circular, kendall.u, trineq.

Examples

data(celebrities)           # absolute choice frequencies
strans(celebrities)         # WST and MST hold, but not SST
strans(celebrities)$pcm     # reordered relative frequencies
strans(celebrities)$violdf  # transitivity violations

Tartness Rankings of Salad Dressings

Description

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.

Usage

data(tartness)

Format

A data frame consisting the rankings and their frequencies.

Source

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

References

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.

Examples

data(tartness)

Thurstone-Mosteller Model (Case V)

Description

Fits a Thurstone-Mosteller model (Case V) by maximum likelihood.

Usage

thurstone(M)

Arguments

Details

thurstone provides a front end for glm. See Critchlow and Fligner (1991) for more details.

Value

References

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

See Also

eba, strans, pcX, kendall.u, circular, glm.

Examples

## 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

Trinary Inequality

Description

Checks if binary choice probabilities fulfill the trinary inequality.

Usage

trineq(M, A = 1:I)

Arguments

Details

For any triple of stimuli x, y, z, the trinary inequality states that, if P(x, y) > 1/2 and (xy)z, then

R(x, y, z) > 1,

where R(x, y, z) = R(x, y) R(y, z) R(z, x), R(x, y) = P(x, y)/P(y, x), and (xy)z denotes that x and y share at least one aspect that z does not have (Tversky and Sattath, 1979, p. 554).

inclusion.rule checks if a family of aspect sets is representable by a tree.

Value

Results checking the trinary inequality.

References

Tversky, A., & Sattath, S. (1979). Preference trees. Psychological Review, 86, 542–573. doi:10.1037/0033-295X.86.6.542

See Also

eba, inclusion.rule, strans.

Examples

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 triple

Utility Scale of an EBA Choice Model

Description

Extract the (normalized) utility scale for an elimination-by-aspects (EBA) model.

Usage

uscale(object, norm = "sum", log = FALSE)

Arguments

Details

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.

Value

The (normalized) utility scale of the stimuli.

See Also

eba, cov.u, wald.test.

Examples

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) := 0

Testing Linear Hypotheses in Elimination-by-Aspects (EBA) Models

Description

Tests linear hypotheses of the form Cp = 0 in elimination-by-aspects (EBA) models using the Wald test.

Usage

wald.test(object, C, u.scale = TRUE)

Arguments

Details

The Wald test statistic,

W = (Cp)' [C cov(p) C']^{-1} (Cp),

is approximately chi-square distributed with rk(C) degrees of freedom.

C is usually of full rank and must have as many columns as there are parameters in p.

Value

See Also

eba, group.test, uscale, cov.u.

Examples

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)

Wine Tasting Data

Description

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.

Usage

data("winetaste")

Format

ambilight A data frame containing 230 observations on 10 variables:

preference, fruitiness, spiciness, sweetness

Paired comparison of class paircomp; judgments for one of the 6 ordered pairs of the blue, red, and white lighting conditions.

age

subject age

gender

factor, subject gender

sensesmell

self-rating of sense of smell and taste.

likewine

self-rating of general liking of wine.

drinkwine

factor, frequency of wine consumption.

redwhite

factor, preference for red or white wine.

redwines A data frame containing 11 observations on 7 variables:

bitterness, fruitiness, sourness, roundness, preference

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.

best

factor, Which of the five wines did you like best?

worst

factor, Which of the five wines did you like worst?

Details

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).

References

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

See Also

eba, eba.order, paircomp.

Examples

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)