# Best stats test for comparing this count data in R

So I've come across a tricky dataset which comprises of several participants who read a vision testing chart using first their right eye, then the left eye. Each line of the chart has 5 letters and as they read down the letters on each subsequent line get smaller (just like a normal opticians letter chart). I have a dataset like below which assigns each of the letters to a column. So the first letter on each line is A, the second column is B etc.. upto the fifth column which is E.

We then counted how many errors (error scores) each participant made in each column:

ID <- c("harry","Ben","sheila","Rocky","Tom","Anna","Lucy","Bob","Derek")
A <- c(0,0,1,1,2,1,0,2,0)
B <- c(2,1,1,2,0,2,0,0,0)
C <- c(0,0,0,0,1,2,1,0,0)
D <- c(0,2,1,1,3,0,0,0,0)
E <- c(1,2,2,3,1,0,0,0,2)
data <- cbind.data.frame(ID,A,B,C,D,E)


My dataset has quite a few more participants than this and as you can see the error scores are generally very low (0,1,2,3), although they could theoretically be higher.
I want to compare column A to column E to see if there is a significant difference. I figured the values are so low that comparing means is not the way to go so I'm looking at comparing the distributions. I thought a McNemar test but as its not dichotomous I'm looking at GLMM or CLMM with ID as a random variable to control for the fact that the same participant per row is reading the chart for each column value.

So far I have converted to long format changed error scores to factors:

library(tidyverse)
data_long <- data %>%
pivot_longer(!ID, names_to = "column", values_to = "Error_score")
data_long$ID = factor(data_long$ID, levels=unique(data_long$ID)) data_long$column = factor(data_long$column, levels=unique(data_long$column))


Then fitted the CLMM using the ordinal package and using artistic licence to say that the error score is ordinal in nature:

library(ordinal)
model_clmm = clmm(factor(Error_score) ~ column + (1|ID),
data = data_long,
threshold = "equidistant")


Then used a post-hoc test using estimated marginal means to compare columns:

library(emmeans)
marginal = emmeans(model_clmm, ~ column)


Now I'm just not sure if this is the correct way to go or whether I should be using GLMM as below:

library(lme4)
data_long$Error_score <- as.numeric((data_long$Error_score))
model_glmm <- glmer(Error_score ~ column + (1 | ID), family = poisson(link = "log"),
data = data_long)
summary(model_glmm)
library(emmeans)
marginal = emmeans(model_glmm, ~ column)


These two methods come up with very different results (p-values) for comparing column error scores. Any help from the community would be so incredibly appreciated.

Just to try and clarify, essentially we are not distinguishing between lines but have tallied up the errors per column to see if participants make more errors on the left side of the chart compared to the right side of the chart. So comparing the total error score for all participants in column A versus column E. Basically like an ANOVA (but something appropriate for this data) to compare error scores in A vs error scores in E.

There are 20 lines on the chart but some could only read the top 4 or 5 lines whereas some could read 15-16 lines. This is why we went for error scores rather than absolute number of letters read, since each person read vastly different number of letters, but all tended to make errors at or near the last line they could read.

So the hypothesis we want to check is that the total error score in column A (left side of chart) is statistically significantly lower than column E (right side of chart)

Thanks all

• Am I correct that you aren't distinguishing among the lines on the chart with decreasing letter sizes, just summing errors over all of the lines? How many lines were there on the chart (i.e., chances to make an error for each letter)? This might be formulated as a binomial model, although I suspect that would end up pretty similar to the Poisson model. It would also help to edit the question to show the discrepancy in model results. Please provide that information by editing the question, as comments are easy to overlook and can be deleted.
– EdM
Commented Oct 12, 2022 at 15:51
• I can't execute the code. emmeans cannot compute the clm vcov matrix. But at a high level, why are you even conducting a statistical test? What's even the hypothesis? Commented Oct 12, 2022 at 15:59
• thanks both, I've edited the code to hopefully actually run this time. Also added comments to hopefully clarify the experiment. Commented Oct 12, 2022 at 16:36

Putting aside for a moment whether the "error score" makes sense, and whether we should be worrying too much about p-values in any event, it's not clear that you are getting such big differences between the models in terms of p-values as you seem to fear.

marginal1 <- emmeans(model_clmm, ~ column)
marginal2 <- emmeans(model_glmm, ~ column)
summary(pairs1)$$p.value; summary(pairs2)$$p.value
#  [1] 0.9993745 0.8933498 0.9992225 0.8886367 0.7996809 0.9905929 0.9571348 0.9668942 0.4109463 0.7942288
#  [1] 0.9990259 0.8996124 1.0000000 0.8833376 0.7895380 0.9990261 0.9596780 0.8996074 0.4136298 0.8833436


Scroll across to see that the corresponding p-values in the two rows (not labeled in this simple display) agree quite well given how different the models are. Maybe you are having more problems with the real data. The coefficients and their contrasts per se will necessarily differ between the two models, as they are modeling different things.

In terms of these 2 models, note that you would presumably have a fixed upper limit to the "error score" values, which might call into question the Poisson assumption. On that basis the ordinal model would seem to be better suited to your data insofar as the assumptions of the clmm model hold. But, in this sample at least, you also seem to have very few high "error score" values, so the Poisson model might not be too bad.

What bothers me is how this "error score" ignores the total number of lines that each individual could read at all. That would seem to be an important aspect of the study to control for. I might, for example, expect to see this modeled instead as a binary model of successes/failures taking all of the rows into account, allowing for inter-individual differences in overall success with the random-effect term. If this score is well established in your field of interest, OK, but I suspect that otherwise you might have trouble convincing a reviewer of your work that this type of "error score" is useful.

• Thank you for this comprehensive answer. The error score data given above was made up but the actual data also contains only 0,1,2 or 3 errors per column. When I run the script with the actual data I get a sig. difference (p value < 0.05) between A and E but no others. The main question I had was really about the type of method I should be using for data like this. Commented Oct 12, 2022 at 18:46
• You are definitely correct that error score is not an established term and the data is not ideal. Unfortunately the raw data is unavailable to go back and assess overall numbers of letters read. It would indeed have been better to control for overall performance differences using a random effect. I also share concerns about the overuse of p values but in the ophthalmology field I'm in, reviewers tend to insist on them. Commented Oct 12, 2022 at 18:51
• I should say with the actual data I get a sig diff between col A and col E when using clmm but not when using glmm Commented Oct 12, 2022 at 19:29
• @holmes the ordinal model would be the easiest to justify given the upper limit of 3 for the outcome. With a Poisson mean of 1 you would expect about 1% of cases to have a value of 4 or more. You don't have any, arguing against a simple Poisson model. If you had hypothesized a difference specifically between colA and colE prior to seeing the data, you don't need to do all 10 pairwise tests (with Tukey adjustment) and thus save a lot pf power. You just compare those two columns in a specified contrast (using the results of the full model, however).
– EdM
Commented Oct 12, 2022 at 19:32
• Thank you EdM. Appreciate your help. Needed some reassurance that I was on roughly the right track. Commented Oct 15, 2022 at 9:28