Agreement between methods with multiple observations per individual

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