get the data into R

If you already know how to read your raw data into R, you can skip this section. We assume that you have your ratings organized in a table, in which each line corresponds to a rating event/trial. There is at least a column for the stimulus that was preferred by the rater and the one that was not. Additional columns are likely to be present and the table below provides an example. Most likely you will have organized this table in a spreadsheet software like Excel or OpenOffice. The perhaps simplest way of reading a data set into R is to first save your data table from the spreadsheet software as tab-delimited text file (File>Save as…> and then choose ‘Tab Delimited Text (.txt)’).

Having such a text-file, it is then easy to read that data set into R and save it as an object named xdata:³

# Windows
xdata <- read.table(file = "c:\\datafiles\\myfile.txt", sep = "\t", header = TRUE) 

# Mac
xdata <- read.table(file = "/Volumes/mydrive/myfile.txt", sep = "\t", header = TRUE)

str(xdata)

Note that you have to specify the full path to the place where the data file is saved.⁴ The sep = "\t" argument tells R that you used a tab to separate entries within a line. The header = TRUE argument tells R that the first line in your table are column headers. Finally, the str(xdata) command gives a brief overview over the data set, which you can use to check whether the import went smoothly (for example, as indicated by the number of lines (obs. in the output of str()), which should correspond to the number of trials in your data set).

If the import was successful, the only thing you need to know is that for the functions of the package to work we need to access the columns of the data table individually. This can be achieved by using the dollar character. For example xdata$preferred.stimulus returns the column with the preferred stimuli. If you are unfamiliar with this, hopefully it will become clear while going through the examples in the next section.

random data

Throughout this first part of the tutorial, we will use randomly generated data sets. We begin by creating such a random data set with the function randompairs(), which we name xdata.

set.seed(123)
xdata <- randompairs(nstim = 7, nint = 700, reverse = 0.1)
head(xdata)
#>   index winner loser
#> 1     1      o     s
#> 2     2      n     r
#> 3     3      j     r
#> 4     4      s     r
#> 5     5      j     r
#> 6     6      k     o

The command set.seed(123) simply ensures that each time you run this it will create the same data set as it shown in this tutorial. If you want to create a truly random data set, just leave out this line.

With this we created a data set with 7 different stimuli (nstim = 7), which were presented in presentations/trials (nint = 700). The head() command displays the first six lines of the newly created data set.⁵ You can see the stimulus IDs: the winner-column refers to the preferred stimulus (the winner of the trial) and the loser-column to the second/unpreferred stimulus (the loser). The index-column simply displays the original order in which the stimuli were presented. Note also that the data set is generated in a way such that there is always a preference for the stimulus that comes first in alphanumeric order, which is reversed in 10% of trials (reverse = 0.1). This ensures that a hierarchy of stimuli actually arises.

calculating the ratings

We then can go on to calculate the actual ratings from this sequence with the core function of this package, elochoice(). Again, to enable you to obtain the exact results as presented in this tutorial, I use the set.seed() function.

set.seed(123)
res <- elochoice(winner = xdata$winner, loser = xdata$loser, runs = 1000)
summary(res)
#> Elo ratings from 7 stimuli
#> total (mean/median) number of rating events: 700 (200/201)
#> range of rating events per stimulus: 187 - 214 
#> startvalue: 0 
#> k: 100 
#> randomizations: 1000

The two code pieces winner = xdata$winner and loser = xdata$loser specify the two columns in our data table that represent the winning (preferred) and losing (not preferred) stimuli, respectively. The runs = 1000 bit indicates how many random sequences of presentation order we want to generate.

We saved the results in an object named res, from which can we can get a brief summary with the summary(res) command. From this, we can see that there were 7 stimuli in the data set, each appearing between 187 and 214 times. The summary also notes that we randomized the original sequence 1000 times. In case you have larger data sets, you may want to reduce the number of randomizations to reduce the computation time and to explore whether everything runs as it should.⁶

Next, we obviously want to see what the actual Elo-ratings are for the stimuli used in the data set. For this, we use the function ratings():

ratings(res, show = "original", drawplot = FALSE)
#>    c    j    k    n    o    r    s 
#>  432  224   82    5 -169 -202 -372

The show = "original" argument specifies that we wish to see the ratings as obtained from the initial (original) data sequence that is present in our data set. If we want to have the ratings averaged across all randomizations (mElo), we change the argument to show = "mean":

ratings(res, show = "mean", drawplot = FALSE)
#>        c        j        k        n        o        r        s 
#>  459.205  247.165  141.170   16.752 -159.900 -235.153 -469.239

Similarly, you can also request all ratings from all randomizations, using show = "all" or return the ranges across all randomizations with show = "range".

If you wish to export ratings, you can use the write.table() function, which saves your results in a text file that can easily be opened in a spreadsheet program. First, we save the rating results into a new object, myratings, which we then export. Note that the resulting text file will be in a ‘long’ format, i.e. each stimulus along with its original or mean rating will appear as one row.

myratings <- ratings(res, show = "mean", drawplot = FALSE)
# Windows
xdata <- write.table(myratings, "c:\\datafiles\\myratings.txt", sep = "\t", header = TRUE)
# Mac
xdata <- write.table(myratings, "/Volumes/mydrive/myratings.txt", sep = "\t", header = TRUE)

If you want to export the ratings from each single randomization (i.e. show = "all") or the range of ratings across all randomizations (show="range"), the layout of the text file will be ‘wide’, i.e. each stimulus appears as its own column with each row representing ratings after one randomization or two rows representing the minimum and maximum rating values. By default, R appends row names to the text output, which is not convenient in this case so we turn this option off with row.names = FALSE}.

myratings <- ratings(res, show = "all", drawplot = FALSE)
# Windows
xdata <- write.table(myratings, "c:\\datafiles\\myratings.txt", sep = "\t", header = TRUE, row.names = FALSE)
# Mac
xdata <- write.table(myratings, "/Volumes/mydrive/myratings.txt", sep = "\t", header = TRUE, row.names = FALSE)

Finally, the ratings() function also allows you to take a first graphical glance at how the randomizations affected the ratings.

ratings(res, show = NULL, drawplot = TRUE)

Elo ratings of 7 stimuli after 700 rating events and 1000 randomizations of the sequence. The black circles represent the average rating at the end of the 1000 generated sequences for each stimulus, and the black lines represent their ranges. The grey circles show the final ratings from the original sequence.

This figure shows the mean ratings across the 1000 randomizations as black circles, while the ratings from the original/initial sequence are indicated by the smaller grey circles The vertical bars represent the ranges of Elo-ratings across the 1000 randomizations for each stimulus. The stimulus with the highest rating is c and the stimulus with the lowest rating is s. Recall that the data generation with randompairs uses an underlying ranking of stimuli that is based on an alphabetical order: hence c should have a higher rating than j, which should be higher rated than k and so on. The rating process recovers this nicely here because there were not too many exceptions to this rule during the data generation (reversals = 0.1).

background

We define an reliability-index as \(R = 1- \frac{\sum{u}}{N}\), where \(N\) is the total number of rating events/trials for which a binary expectation for the outcome of the trial existed⁷ and \(u\) is a vector containing 0’s and 1’s, in which a \(0\) indicates that the preference in this trial was according to the expectation (i.e. the stimulus with the higher Elo-rating before the trial was preferred), and a \(1\) indicates a trial in which the expectation was violated, i.e. the stimulus with the lower Elo-rating before the trial was preferred (an ‘upset’). In other words, \(R\) is the proportion of trials that went in accordance with the expectation. Note that trials without any expectation, i.e. those for which ratings for both stimuli are identical, are excluded from the calculation.⁸ Subtracting the proportion from \(1\) is done to ensure that if there are no upsets (\(\sum {u}=0\)), the index is \(1\), thereby indicating complete agreement between expectation and observed rating events. For example, in the table below there are ten rating events/trials, four of which go against the expectation (i.e. they are upsets), yielding an reliability-index of \(1-4/10=0.6\).

This approach can be extended to calculate a weighted reliability index, where the weight is given by the absolute Elo-rating difference, an index we denote \(R'\) and which is defined as \(R' = 1- \sum_{i=1}^{N} {\frac{u_i * w_i}{\sum{w}}}\), where \(u_i\) is the same vector of 0’s and 1’s as described above and \(w_i\) is the absolute Elo-rating difference, i.e. the weight. Following this logic, stronger violations of the expectation contribute stronger to the reliability index than smaller violations. For example, column 3 in the table below contains fictional rating differences (absolute values), which for illustrative purposes are assigned in a way such that the four largest rating differences (200, 200, 280, 300) correspond to the four upsets, which should lead to a smaller reliability index as compared to the simpler version described earlier. Applying this leads to \(R'=1-0.57=0.43\). In contrast, if we apply the smallest rating differences (90, 100, 120, 140, as per column 4 in the table) to the upsets, this should lead to a larger reliability-index, which it does: \(R'=1-0.26=0.74\).

In sum, the reliability index represents how many rating events went against the expectations. The weighted version refines this measure by integrating the strength of deviations from expectations, such that many large deviations will lead to smaller index values, while many small deviations will lead to larger index values.

application

To calculate the reliability-index, we use the function reliability(). Note that we calculated our initial Elo-ratings based on 1000 randomizations, so to save space, I’ll display only the first six lines of the results (the head(...) function).

upsets <- reliability(res)
head(upsets)
#>       upset upset.wgt totIA
#> 1 0.8581662 0.8847880   698
#> 2 0.8610315 0.8885191   698
#> 3 0.8553009 0.8876968   698
#> 4 0.8569385 0.8913181   699
#> 5 0.8640916 0.8858861   699
#> 6 0.8647482 0.8842455   695

Each line in this table represents one randomization.⁹ The first column represents the unweighted and the second the weighted reliability index (\(R\) and \(R'\)), which is followed by the total number of trials that contributed to the calculation of the index. Note that this number cannot reach 1000 (our total number of trials in the data set) because at least for the very first trial we did not have an expectation for the outcome of that trial.

We then calculate the average values for both the unweighted and weighted upset indices.

mean(upsets$upset)
#> [1] 0.8548718
mean(upsets$upset.wgt)
#> [1] 0.886231

Remember that our data set contained a fairly low number of ‘reversals’, i.e. 10% of trials went against the predefined preference (e.g. ‘A’ is preferred over ‘K’). As such, the \(R\) and \(R'\) values are fairly high (\(\sim 0.8 - 0.9\)). If we create another data set, in which reversals are more common, we can see that the values go down.

set.seed(123)
xdata <- randompairs(nstim = 7, nint = 700, reverse = 0.3)
res <- elochoice(winner = xdata$winner, loser = xdata$loser, runs = 1000)
upsets <- reliability(res)
mean(upsets$upset)
#> [1] 0.6084496
mean(upsets$upset.wgt)
#> [1] 0.6462743

how many raters?¹⁰

A note of caution. The functions described in the following section have not been very thoroughly tested (yet). If they fail or produce confusing results, please let us know!

We may also wonder how many raters are needed to achieve stability in ratings. The function raterprog() allows assessing this. The idea here is to start with one rater, calculate reliability, add the second rater, re-calculate reliability, add the third rater, re-calculate reliability, etc. The idea is that as more raters are included, reliability should converge on some ‘final’ value. Please note that for now the function only returns the weighted reliability index \(R'\).

In order to calculate this progressive reliability, obviously you need a column in your data set that reflects rater IDs. For this demonstration we use a subset¹¹ of a data set that contains this information and is described further below.

data(physical)
# limit to 10 raters
physical <- subset(physical, raterID %in% c(1, 2, 8, 10, 11, 12, 23, 27, 31, 47))
set.seed(123)
res <- raterprog(physical$Winner, physical$Loser, physical$raterID, progbar = FALSE)

raterprogplot(res)

Reliability index \(R'\) as function of number of raters included in the rating process. In this example, it appears as if including 5 raters is sufficient, as additional raters improve rating reliability relatively little.

When we look at these results graphically, we may conclude that with this data set using the data from the first 5 raters may be sufficient to achieve high reliability (e.g. larger than 0.8).

It may be advisable, however, to not necessarily rely on the original sequence of raters involved. In the general spirit of the package, we can randomize the order with which raters are included. We can do this by increasing the ratershuffle= argument from its default value.¹² Note that doing so will increase computation time even further, so start with small values if applying this to your own data to test whether the functions work as intended.

set.seed(123)
res <- raterprog(physical$Winner, physical$Loser, physical$raterID, progbar = FALSE, ratershuffle = 10)

raterprogplot(res)

Reliability index \(R'\) as function of number of raters included in the rating process. Here, we used the original rater order and an additional 9 random orders. Grey bars reflect quartiles, grey points are the results from the original rater order (compare to figure~ ef{fig:raterprog1}), and black points are the average values of the 10 rater orders used. The interpretation would likely to be the same as for figure~ ef{fig:raterprog1}.

The results of this figure probably would be interpreted in a quite similar same way as those of the figure above: from four or five raters onwards, we achieve reliability larger than 0.8. One thing to note here is the original rater sequence (the grey circles fall outside the quartiles of the reshuffled orders for the low rater numbers). This is not a problem per se, but just illustrates that the original rater order might not be very representative. On other other thing to note here is that the variation in reliability becomes smaller with larger rater numbers.