配置环境

• 基础编程环境

conda create -n icc conda-forge::r-tidyverse conda-forge::r-irkernel conda-forge::r-lme4 conda-forge::r-psych conda-forge::r-blandaltmanleh -y
conda run -n icc Rscript -e "IRkernel::installspec(name='icc', displayname='r-icc')"

使用示例

require(psych)
require(BlandAltmanLeh)
require(tidyverse)
require(RColorBrewer)
options(repr.plot.width=6, repr.plot.height=12)
RColorBrewer::display.brewer.all(type = "all")

R1 = readr::read_tsv(na = '/', locale = locale(encoding = 'GBK'), show_col_types = FALSE,
                                 file = './data/R1.tsv')
R2 = readr::read_tsv(na = '/', locale = locale(encoding = 'GBK'), show_col_types = FALSE,
                                 file = './data/R2.tsv')
R1 %>% head

# Single 衡量一次测量的可靠性 Average 最终结果是多次评分的平均值
# raters_absolute 每个对象由随机抽取的观察者评分
# random_raters 所有对象由同一组观察者的评分
# fixed_raters 所有对象由固定的观察者评分
f_check_icc_df <- function(R1_, R2_, vars_, type='Single_fixed_raters'){
    res <- data.frame(var = vars_, ICC = rep(0, length(vars_)), F=rep(0, length(vars_)), 
                     p.value = rep(1, length(vars_)), lower=rep(0, length(vars_)), upper=rep(1, length(vars_)))
    for (i in 1:nrow(res)) {
        var = res$var[i]
        data = data.frame(R1_[[var]], R2_[[var]])
        result <- psych::ICC(data)
        result <- result$results[type,]
        res[i, c('ICC', 'F', 'p.value', 'lower', 'upper')] <- result[c('ICC', 'F', 'p', 'lower bound', 'upper bound')]
    }
    return(res)
}

vars = colnames(R1)[-1]
ICC_result <- f_check_icc_df(R1, R2, vars)
ICC_result %>% head
readr::write_tsv(x = ICC_result, file = './data/icc.tsv')

bland.altman.ggplot2 <- function (group1, group2, var_n, two = 1.96, mode = 1, conf.int = 0, geom_count = FALSE) 
{
    ba <- BlandAltmanLeh::bland.altman.stats(group1 = group1, group2 = group2, 
        two = two, mode = mode, conf.int = conf.int)
    values <- data.frame(m = ba$means, d = ba$diffs)
    geom_my <- if (geom_count) 
        ggplot2::geom_count
    else ggplot2::geom_point
    m <- NULL
    d <- NULL
    p <- ggplot2::ggplot(values, ggplot2::aes(x = m, y = d)) + geom_my() + 
    geom_point(color = RColorBrewer::brewer.pal(8, 'Set1')[5], size = 2) +
    theme_minimal() + labs(title = paste0("Inter-observer Agreement ", var_n)) +
    ggplot2::xlab("Mean of measurements") + ggplot2::ylab("Difference") +
    theme(
        axis.title.x = element_text(size = 16),  # X轴标题字体大小
        axis.title.y = element_text(size = 16),  # Y轴标题字体大小
        axis.text.x = element_text(size = 16),   # X轴刻度字体大小
        axis.text.y = element_text(size = 16),   # Y轴刻度字体大小
        panel.border = element_rect(color = "black", fill = NA, linewidth = 1),  # 方形边框
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5)  # 居中加粗标题
    ) + 
    ggplot2::geom_hline(yintercept = ba$lines['mean.diffs'], linetype = "dashed", linewidth = 1, color = RColorBrewer::brewer.pal(8, 'Set1')[8]) + 
    ggplot2::geom_hline(yintercept = ba$lines['lower.limit'], linetype = "dashed", linewidth = 1, color = RColorBrewer::brewer.pal(8, 'Set1')[3]) + 
    ggplot2::geom_hline(yintercept = ba$lines['upper.limit'], linetype = "dashed", linewidth = 1, color = RColorBrewer::brewer.pal(8, 'Set1')[1]) +
    ggplot2::annotate("rect", xmin = -Inf, xmax = Inf, ymin = ba$lines['lower.limit'], ymax = ba$lines['upper.limit'],
           alpha = 0.2, fill = "lightblue")
    
    return(p)
}

options(repr.plot.width=6, repr.plot.height=4)
vars = colnames(R1)[-1]
for (var in vars){
    bap <- bland.altman.ggplot2(
        R1[[var]], R2[[var]], var_n=var
    )
    print(bap)
}