Test-retest reliability - Statistics#

Test-retest metrics analysis#

suppressWarnings(suppressMessages(library("car")))
suppressWarnings(suppressMessages(library("sjstats")))
suppressWarnings(suppressMessages(library("broom")))
suppressWarnings(suppressMessages(library("multcomp")))

library("car")
library("broom")
library("multcomp")
library("sjstats")
library("xtable")

path_resources <- file.path(getwd(), "..", "..", "resources", "data")
df_fs <- read.csv(file.path(path_resources, "FREESURFER_r2values.csv"))
df_cat <- read.csv(file.path(path_resources, "ACPC_CAT12_r2values.csv"))

data_col <- rbind(df_fs, df_cat)
data_col$software <- as.factor(data_col$software)

map_aov <- list("r_square"=aov(r_square ~ software + age + snr_total_mean, data = data_col), 
                "ICC"=aov(ICC ~ software + age + snr_total_mean, data = data_col))


for (rep_measure in c("r_square", "ICC")){
  
  print(sprintf("## #####################      %s      ################################", rep_measure))  
  
  my_anova <- map_aov[[rep_measure]]
  print(Anova(my_anova, type = "III"))
    
  latex_table_name = sprintf("test_retest_a2009s_%s_OASIS_first_level.txt", rep_measure)
  df = anova_stats(my_anova)[, c("df", "sumsq", "statistic", "p.value", "cohens.f")]
  print(xtable(df, caption = sprintf("ABIDEI III %s", rep_measure)), file=latex_table_name)

}
[1] "## #####################      r_square      ################################"
Anova Table (Type III tests)

Response: r_square
                Sum Sq  Df  F value Pr(>F)    
(Intercept)    0.31275   1  90.5660 <2e-16 ***
software       0.01503   1   4.3519 0.0374 *  
age            0.00843   1   2.4424 0.1186    
snr_total_mean 0.53121   1 153.8303 <2e-16 ***
Residuals      2.03050 588                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
[1] "## #####################      ICC      ################################"
Anova Table (Type III tests)

Response: ICC
                Sum Sq  Df  F value Pr(>F)    
(Intercept)    0.69861   1 576.1004 <2e-16 ***
software       0.00404   1   3.3333 0.0684 .  
age            0.00317   1   2.6123 0.1066    
snr_total_mean 0.16963   1 139.8833 <2e-16 ***
Residuals      0.71304 588                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Analysis of age effect on test-retest metrics#

for (rep_measure in c("r_square", "ICC")){
  print(sprintf("############### %s #################", rep_measure))
  equation = sprintf("%s ~ age + snr_total_mean", rep_measure)

  model_fs <- lm(equation, data = df_fs)
  model_cat <- lm(equation, data = df_cat)
  
  
  print("FREESURFER")
  print(summary(model_fs))

  par(mfrow = c(2, 2))
  #plot(model_fs)
    
    
  print("CAT12")
  print(summary(model_cat))

  #par(mfrow = c(2, 2))
  #plot(model_cat)
  for (var_name in c("age", "snr_total_mean")){

      
      n1 = nrow(df_fs)
      n2 = nrow(df_cat)
      degrees_freedom = n1+n1-2

      m1 = summary(model_fs)[["coefficients"]][var_name,"Estimate"]
      m2 = summary(model_cat)[["coefficients"]][var_name,"Estimate"]

      s1 =  summary(model_fs)[["coefficients"]][var_name, "Std. Error"] 
      s2 =  summary(model_cat)[["coefficients"]][var_name, "Std. Error"] 

      mean_diff = m1 - m2

      sp = (((n1-1)*(s1**2) + (n2-1)*(s2**2))/(n1+n2-2))**0.5
      cohen_d <- mean_diff/sp

      dem_samples = (1/n1+1/n2)**0.5
      t = mean_diff/(sp*dem_samples)

      p_value <- pt(-abs(t),degrees_freedom)

      t_critical_value <- qt(1-.05/2, degrees_freedom)

      print(sprintf("\n----------- %s  Slope: Independent t-test -----------\n", var_name))  
      print(sprintf("%s-%s -> ndegree: %d diff: %f t-value: %f, t-critical-value:%f, p-value: %e", "FS", "CAT12", n1+n2-2, mean_diff, t, t_critical_value, p_value))
      cat("Cohen's d effect %d", cohen_d)
      cat("\n\n")
  }
}
[1] "############### r_square #################"
[1] "FREESURFER"

Call:
lm(formula = equation, data = df_fs)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.48284 -0.00963  0.01062  0.03161  0.09143 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    0.4619635  0.0739683   6.245 1.48e-09 ***
age            0.0003258  0.0004657   0.700    0.485    
snr_total_mean 0.0472921  0.0056473   8.374 2.33e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06254 on 293 degrees of freedom
Multiple R-squared:  0.222,	Adjusted R-squared:  0.2167 
F-statistic:  41.8 on 2 and 293 DF,  p-value: < 2.2e-16

[1] "CAT12"

Call:
lm(formula = equation, data = df_cat)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.45921 -0.01085  0.00608  0.02597  0.08161 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    0.4644845  0.0649180   7.155 6.74e-12 ***
age            0.0006414  0.0004088   1.569    0.118    
snr_total_mean 0.0457822  0.0049563   9.237  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.05489 on 293 degrees of freedom
Multiple R-squared:  0.2453,	Adjusted R-squared:  0.2401 
F-statistic:  47.6 on 2 and 293 DF,  p-value: < 2.2e-16

[1] "\n----------- age  Slope: Independent t-test -----------\n"
[1] "FS-CAT12 -> ndegree: 590 diff: -0.000316 t-value: -8.762571, t-critical-value:1.963993, p-value: 1.011985e-17"
Cohen's d effect %d -0.7202789

[1] "\n----------- snr_total_mean  Slope: Independent t-test -----------\n"
[1] "FS-CAT12 -> ndegree: 590 diff: 0.001510 t-value: 3.457183, t-critical-value:1.963993, p-value: 2.925716e-04"
Cohen's d effect %d 0.2841787

[1] "############### ICC #################"
[1] "FREESURFER"

Call:
lm(formula = equation, data = df_fs)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.307774 -0.005504  0.005264  0.017648  0.051610 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    0.6993503  0.0429895  16.268  < 2e-16 ***
age            0.0002070  0.0002707   0.765    0.445    
snr_total_mean 0.0264509  0.0032822   8.059 1.97e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.03635 on 293 degrees of freedom
Multiple R-squared:  0.2075,	Adjusted R-squared:  0.2021 
F-statistic: 38.35 on 2 and 293 DF,  p-value: 1.609e-15

[1] "CAT12"

Call:
lm(formula = equation, data = df_cat)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.298704 -0.006643  0.003960  0.014671  0.045904 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    0.6951402  0.0394216  17.634  < 2e-16 ***
age            0.0003857  0.0002482   1.554    0.121    
snr_total_mean 0.0261443  0.0030097   8.687 2.68e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.03333 on 293 degrees of freedom
Multiple R-squared:  0.2222,	Adjusted R-squared:  0.2168 
F-statistic: 41.84 on 2 and 293 DF,  p-value: < 2.2e-16

[1] "\n----------- age  Slope: Independent t-test -----------\n"
[1] "FS-CAT12 -> ndegree: 590 diff: -0.000179 t-value: -8.371640, t-critical-value:1.963993, p-value: 2.069315e-16"
Cohen's d effect %d -0.6881446

[1] "\n----------- snr_total_mean  Slope: Independent t-test -----------\n"
[1] "FS-CAT12 -> ndegree: 590 diff: 0.000307 t-value: 1.184747, t-critical-value:1.963993, p-value: 1.182972e-01"
Cohen's d effect %d 0.09738563