Agreement between methods with multiple observations per individual

Question

I am having trouble interpreting agreement analysis in R.

I have a dataset similar to that below with multiple paired observations per subject:

x <- c(4,6,3,2,6,7,8,4,3,2,6,7,8,3,3,6,8,2,5,1,6,8,7,1,4)
y <- c(7,7,3,6,3,7,7,2,4,1,3,5,0,6,3,2,1,2,8,7,3,3,4,6,3)
id <- x <- c(1,1,1,1,2,2,3,3,3,3,4,4,4,4,4,4,5,6,7,8,9,9,9,10,10)

I want to assess the agreement between method x and y using Bland and Altman stats. x is the reference/gold standard method. id represents the subjects individual id number. According to Bland and Altman publication (https://www-users.york.ac.uk/~mb55/meas/bland2007.pdf) I should examine the mean variance within subjects and between subjects using one-way ANOVA:

res.aov <- aov((x-y) ~ factor(id), data = df)
summary(res.aov)
            Df Sum Sq Mean Sq F value   Pr(>F)    
factor(id)   9 296.41   32.93   7.292 0.000435 ***
Residuals   15  67.75    4.52                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

I interpret this output as there is statistical difference in variance of difference between measurements using method x and y within each subject. So next, I have to examine the mean bias and limits-of-agreement (LOA). First, if I ignore the subjects and calculating means and LOA for all paired observations, as if they were from different subjects I get the following:

ba <- bland.altman.stats(df$x,df$y)
ba$lines
lower.limit  mean.diffs upper.limit 
-7.074781    0.560000    8.194781 

When I do not ignore the subjects I first calculate the means of measurements within each subject before before calculating mean bias and LOA:

df2 <- aggregate(df[,2:3], list(df$id), mean, na.rm=TRUE, na.action=NULL)
# And then calculating mean bias and LOA: 
ba2 <- bland.altman.stats(df2$x,df2$y)
ba2$lines
lower.limit  mean.diffs upper.limit 
-5.833102    1.175000    8.183102 

So, the LOA are wider if I ignore the subjects and treat the observations as different subjects. Does that mean that I should do this?

Moreover, in one parameter (not visible here) the p-value from one-way ANOVA was >0.05, but when comparing the methods mean bias was different and LOA was wider when observations were treated as different subjects compared to if I did not ignore subjects. What should I do in this case?

Answer 1

Your question is somewhat unclear, you should have given more detailed information, such as

• Are the observations for same id really paired, that is, each pair taken at the same time, say, or not? This has consequences for how the data could be analyzed.

• Are the observations for same id given in time order (or some other relevant order)?

• You should have given the name of the R package you are using. There are at least two packages on CRAN for Bland-Altman method comparisons! blandr and BlandAltmanLeh. It seems (from function names) that you are using the last one, but you should not let us guessing!

So, some comments: After the first anova you present, you say I interpret this output as there is statistical difference in variance of difference between measurements using method x and y within each subject. I'm not sure what you try to say (what is statistical difference in variance of difference?) but, what it says is that there is indeed difference between the subjects. That is not so much, it would be shocking in this context if variance within subjects where not (much) smaller than the total variance. Anyhow, in this context the F test given by the anova is not of much interest, the anova is just used as an algorithm giving some numbers used in the further analysis (see the paper you linked.)

Before continuing, you must be clear about the modeling assumptions I asked about above! You use the function bland.altman.stats, but from its help page it does not seem that it can be used for your situation with multiple observations per id, so commenting on the results will not be helpful. But, not, So, the LOA are wider if I ignore the subjects and treat the observations as different subjects. Does that mean that I should do this? You should not use the analysis method that gives wider limits, you should use the one appropriate for your data structure! And, it is not even clear that the R package you have chosen can offer you that.

Finally, I will give some plots you could have done with your data as a help, clarifying some questions one should make:

x <- c(4,6,3,2,6,7,8,4,3,2,6,7,8,3,3,6,8,2,5,1,6,8,7,1,4)
y <- c(7,7,3,6,3,7,7,2,4,1,3,5,0,6,3,2,1,2,8,7,3,3,4,6,3)
id  <- c(1,1,1,1,2,2,3,3,3,3,4,4,4,4,4,4,5,6,7,8,9,9,9,10,10)

library(tidyverse)

mydf <- tibble(x, y, id=as.factor(id))

ggplot(mydf, aes( (x+y)/2, y-x, color=id)) + geom_point() +
    geom_hline(yintercept=0, color="red") + ggtitle("Bland-Altman plot for each pair")

Observe that there is quite a lot of variation in the x coordinates within same id, so maybe what we are measuring is not quite constant?

If you are willing to assume that we are really measuring a constant, in which case it is reasonable to use the mean of the (x+y)/2 as x-coordinate, the resulting Bland_Altman plot is

mean_by_id <- mydf %>% group_by(id) %>% summarize(meanmean=mean( (x+y)/2))

### Then joining it to mydf:

mydf <- mydf %>% left_join(mean_by_id, by="id")

ggplot(mydf, aes(meanmean, y-x, color=id))+geom_point() +
    geom_hline(yintercept=0, color="red")+ggtitle("Bland-Altman plot wit a common mean measurement within id")

Finally, a plot to see if there could be some time order (within id):

mydf <- mydf %>% mutate(row=as.numeric(row.names(mydf)))
minrow_by_id <- mydf %>% group_by(id) %>% summarize(minrow=min(row))
mydf <- mydf %>% left_join(minrow_by_id, by="id")
mydf <- mydf %>% mutate(order=row+1-minrow, row=NULL)

ggplot(mydf, aes(order, y-x, color=id, group=id)) + geom_point()  +geom_line() +
    ggtitle("Is there some time order within id?")