Useful links
https://benwhalley.github.io/just-enough-r/icc-and-vpc.html
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")