Useful links

https://benwhalley.github.io/just-enough-r/icc-and-vpc.html

https://neoacademic.com/2011/11/16/computing-intraclass-correlations-icc-as-estimates-of-interrater-reliability-in-spss/

http://cran.nexr.com/web/packages/sjstats/vignettes/mixedmodels-statistics.html

https://cran.r-project.org/web/packages/SimplyAgree/vignettes/intro_vignette.html

Data

#data_all <- read_excel("/Users/ferrec14/Library/Mobile Documents/com~apple~CloudDocs/Documents/H-LAND/t1-t2-flair/anatacomp_stats.xlsx", sheet = "newstats2")
data_volume <- read_excel("/Users/ferrec14/Library/Mobile Documents/com~apple~CloudDocs/Documents/H-LAND/t1-t2-flair/anatacomp_stats.xlsx", sheet = "data_volume")
## New names:
## * `` -> ...9

List of volume estimates

Kruskal-wallis, dunn test, bonferroni

kw <- lapply(roi_volume, #roi_volume or roi_thickness list
        function(x)
         kruskal.test(as.formula(paste0(x,"~ factor(seq)")), data=data_volume)) 

do.call(rbind.data.frame, kw) %>% 
  write.csv("results_kruskal.csv")


dunn <- lapply(roi_volume, #roi_volume or roi_thickness list
        function(x)
        FSA::dunnTest(as.formula(paste0(x,"~ as.factor(seq)")), data=data_volume, method="bonferroni"))

dunn_all <- list()
for (i in 1:length(dunn)){
    print(dunn[[i]])
  dunn_all[[i]] <- dunn[[i]]
}

Mixed model + Manual ICC ICC is the variance explained by the grouping (subject) in relation to the total variance. “How strongly units in the same group (within each subject) resemble each other”. ICC explains the variation between values (t1, t2, flair) within the same group (subject) Reliability: extent at which measures can be replicated, takes into account degree of correlation and agreement between measures (ratio of true variance over total variance) - Pearson measures correlation - Bland altman, t test measures agreement

map_dfr(roi_volume,
                  function(x){
                    formula_mlm = as.formula(paste0(x,"~ factor(seq) +  factor(session) +(1|subnum) ")); #can modify the model
                    #https://rpsychologist.com/r-guide-longitudinal-lme-lmer
                    model_fit = lmer(formula_mlm,data=data_volume);
                    re_variances = VarCorr(model_fit,comp="Variance") %>% 
                      data.frame() %>% 
                      dplyr::mutate(variable = x);
                    return(re_variances)
                    
                  }) %>% 
  dplyr::select(variable,grp,vcov) %>%  #data frame with a column with all rois, columns for each factor from the model with variances for each roi, and a final column of resiaul variance
  pivot_wider(names_from="grp",values_from="vcov") %>%
  dplyr::mutate(icc = subnum/(subnum +Residual))  #creates new columns with ICC for each roi 
  #write.csv("results_data_volume_onseq.csv")

ICC using SimpleAgree package

library(SimplyAgree)
#https://cran.r-project.org/web/packages/SimplyAgree/vignettes/intro_vignette.html
data_volume_reli <- data_volume %>% dplyr::select(subnum, roi_volume)
test1 = reli_stats(
  data = data_volume_reli,
  wide = TRUE,
  col.names = roi_volume
)
test1 #results are the same regardless of what the data.frame included (roi + seq or subnum or both or etc )

ICC using irr package

library("irr")
#data_volume_long <- data_volume %>% pivot_longer(cols = rh_G_S_frontomargin_volume:WM_hypointensities ,names_to= "region")
data_volume_icc <- read_excel("/Users/ferrec14/Library/Mobile Documents/com~apple~CloudDocs/Documents/H-LAND/t1-t2-flair/dat_volume_long.xlsx")

#data_volume_icc$seq_t1 = apply(data_volume_icc, 1, function(row) { if (row['seq'] == "1") row['value'] else NA })
#data_volume_icc$seq_t2 = apply(data_volume_icc, 1, function(row) { if (row['seq'] == "2") row['value'] else NA })
#data_volume_icc$seq_flair = apply(data_volume_icc, 1, function(row) { if (row['seq'] == "3") row['value'] else NA })

data_volume_icc_limited <- data_volume_icc[c("seq_t1", "seq_t2", "seq_flair")]
irr::icc(data_volume_icc_limited,  model="twoway",type = "agreement", unit = "single")

ICC using performance package

model_icc_performance <- list()
for (i in roi_volume) {
    f <- formula(paste(i,"~  factor(seq) +  factor(session) +(1|session) +(1|subnum) "))
    model_icc_performance[[i]] <- lmer(f, data=data_volume)  #runs a lmer model for each regio
    }
iccs <- lapply(model_icc_performance,
       function(x){
         performance::icc(x, by_group = TRUE) #applies icc to each model, output is the ICC for each group (random intercept)
       }) 
do.call(rbind.data.frame, iccs) 
#write.csv("icc_performance_volume.csv") 

ICC across sessions #same as within subjects but taking “session” as group

data_volume_t1 <- subset(data_volume, data_volume$for_t1_atleast2 =="1") 
data_volume_t2 <- subset(data_volume, data_volume$for_t2_atleast2 =="1") 
data_volume_flair <- subset(data_volume, data_volume$for_flair_atleast2 =="1")