Description of the procedures and analysis present in Manuscript 2, Establishing a baseline for human cortical folding morphological variables: a multicenter study, at the Doctorate Thesis presented to the Programa de Pós-Graduação em Ciências Médicas at the Instituto D’Or de Pesquisa e Ensino as a partial requirement to obtain the Doctorate Degree.

Part of the data used here cannot be shared due to restrictions of the Ethic Committee. Data can be shared upon reasonable request to the corresponding author. To fulfill these limitation, we will generate random data to simulate the results.

Get in touch with us () in case any help is needed, our aim is to improve the code as needed!

1 RMARKDOWN AND R SETUP

setwd("~/GitHub/Typical-values")
## define functions

# test angular coeficinet versus theoretical value
test_coef <- function(reg, coefnum, val){
  co <- coef(summary(reg))
  tstat <- (co[coefnum,1] - val)/co[coefnum,2]
  2 * pt(abs(tstat), reg$df.residual, lower.tail = FALSE)
}

# wrap text
wrapper <- function(x, ...) paste(strwrap(x, ...), collapse = "\n")

is_outlier <- function(x, coef = 10) {
  stopifnot(is.numeric(x))
  y <- grDevices::boxplot.stats(x, coef = coef)$out
  x <- x %in% y
  return(x)
}
library(readr)
library(tidyverse)
library(lubridate)
library(ggpubr)
library(kableExtra)
library(broom)
library(lme4)
library(lsmeans)
library(MuMIn)
library(arm)
library(effects)
# COLOR BLIND PALETTE WITH BLACK
cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
cbbPalette2 <- c("#D55E00", "#E69F00", "#56B4E9", "#0072B2", "#CC79A7", "#009E73", "#F0E442")

2 set seed for random process

set.seed(1)
# dados_datasetscomp <- read_csv("dados_datasetscomp2.csv")
dados_datasetscomp <- read_csv("data_typical_values.csv")

dados_datasetscomp <- dados_datasetscomp %>%
  filter(Sample != "HCP500r", Diagnostic != "MCI")

3 DATA PREPARATION

# estimate cortical folding variables
dados_datasetscomp <- dados_datasetscomp %>%
  mutate(
    # create new variables
    GMvolume = ifelse(!is.na(GMvolume),GMvolume,AvgThickness*TotalArea),
    logAvgThickness = log10(AvgThickness),
    logTotalArea = log10(TotalArea),
    logExposedArea = log10(ExposedArea),
    localGI = TotalArea / ExposedArea,
    k = sqrt(AvgThickness) * TotalArea / (ExposedArea ^ 1.25),
    K = 1 / 4 * log10(AvgThickness ^ 2)  + log10(TotalArea) - 5 / 4 * log10(ExposedArea),
    S = 3 / 2 * log10(TotalArea) + 3 / 4 * log10(ExposedArea) - 9 /  4 * log10(AvgThickness ^
                                                                                 2) ,
    I = log10(TotalArea) + log10(ExposedArea) + log10(AvgThickness ^ 2),
    # c = as.double(ifelse(
    #   ROI == "hemisphere", NA, 4 * pi / GaussianCurvature
    # )),
    TotalArea_corrected = ifelse(ROI == "hemisphere", TotalArea, TotalArea * c),
    ExposedArea_corrected = ifelse(ROI == "hemisphere", ExposedArea, ExposedArea * c),
    logTotalArea_corrected = log10(TotalArea_corrected),
    logExposedArea_corrected = log10(ExposedArea_corrected),
    localGI_corrected = ifelse(
      ROI == "hemisphere",
      TotalArea / ExposedArea,
      TotalArea_corrected / ExposedArea_corrected
    ),
    k_corrected = ifelse(
      ROI == "hemisphere",
      sqrt(AvgThickness) * log10(TotalArea) / (log10(ExposedArea) ^ 1.25),
      sqrt(AvgThickness) * log10(TotalArea_corrected) / (log10(ExposedArea_corrected ^
                                                                 1.25))
    ),
    K_corrected =  ifelse(
      ROI == "hemisphere",
      1 / 4 * log10(AvgThickness ^ 2) + log10(TotalArea) - 5 / 4 * log10(ExposedArea),
      1 / 4 * log10(AvgThickness ^ 2) + log10(TotalArea_corrected) - 5 / 4 * log10(ExposedArea_corrected)
    ),
    I_corrected = ifelse(
      ROI == "hemisphere",
      log10(TotalArea) + log10(ExposedArea) + log10(AvgThickness ^ 2) ,
      log10(TotalArea_corrected) + log10(ExposedArea_corrected) + log10(AvgThickness ^ 2)
    ),
    S_corrected = ifelse(
      ROI == "hemisphere",
      3 / 2 * log10(TotalArea) + 3 / 4 * log10(ExposedArea) - 9 /  4 * log10(AvgThickness ^ 2) ,
      3 / 2 * log10(TotalArea_corrected) + 3 / 4 * log10(ExposedArea_corrected) - 9 /  4 * log10(AvgThickness ^ 2)
    ),
    Knorm = K_corrected / sqrt(1 + (1 / 4) ^ 2 + (5 / 4) ^ 2),
    Snorm = S_corrected / sqrt((3 / 2) ^ 2 + (3 / 4) ^ 2 + (9 / 4) ^ 2),
    Inorm = I_corrected / sqrt(1 ^ 2 + 1 ^ 2 + 1 ^ 1)
  )

# create age intervals
dados_datasetscomp$Age_interval <- cut(dados_datasetscomp$Age,
                                       breaks = c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100),
                                       right = FALSE,
                                       include.lowest = TRUE)

dados_datasetscomp$Age_interval10 <- cut(dados_datasetscomp$Age,
                                         breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100),
                                         right = FALSE,
                                         include.lowest = TRUE)
dados_all <- dados_datasetscomp %>% filter(!is.na(logAvgThickness), ExposedArea != 0 | !is.na(localGI), !is.infinite(logExposedArea)) %>% 
  droplevels()

dados_datasetscomp <- dados_all

dados_datasetscomp$Diagnostic <- as.factor(dados_datasetscomp$Diagnostic)
dados_datasetscomp$Diagnostic <- relevel(dados_datasetscomp$Diagnostic, ref = "CTL")

4 Deaging

# define age for deaging
Age.cor = 25

# DEAGING + HARMONIZATION FULL DATA ----

dados_datasetscomp_rate <-
  filter(dados_datasetscomp, Diagnostic == "CTL")

dados_datasetscomp_rate$Sample <-
  as.factor(dados_datasetscomp_rate$Sample)

dados_datasetscomp_rate$ROI <-
  factor(dados_datasetscomp_rate$ROI,
         levels = c("hemisphere", "F", "O", "P", "T"))

4.1 GM volume

m.1 <-
  lme4::lmer(GMvolume ~ Age * ROI + (1|Sample:ROI), data = dados_datasetscomp_rate)

re <- as_tibble(ranef(m.1)) %>%
  filter(grpvar == "Sample:ROI") %>%
  mutate(
    GM_shift = condval,
    sd_GM_shift = condsd,
    Sample = str_split(grp, pattern = ":", simplify = TRUE)[, 1],
    ROI = str_split(grp, pattern = ":", simplify = TRUE)[, 2]
  ) %>%
  dplyr::select(-c(condval, grpvar, term, condsd, grp))

Age.trend <- as_tibble(lstrends(m.1, ~ ROI, var = "Age")) %>%
  mutate(Age.trend_GM = Age.trend) %>%
  dplyr::select(c(ROI, Age.trend_GM))

dados_datasetscomp <- full_join(dados_datasetscomp, Age.trend) %>%
  full_join(re) %>%
  mutate(
    GMvolume_shiftc = GMvolume - GM_shift,
    GMvolume_age_decay = GMvolume - Age.trend_GM * (Age - Age.cor),
    GMvolume_age_decay_shiftc = GMvolume - GM_shift - Age.trend_GM *
      (Age - Age.cor)
  )

4.2 AvgThickness

m.1 <-
  lme4::lmer(AvgThickness ~ Age * ROI + (1|Sample:ROI), data = dados_datasetscomp_rate)

re <- as_tibble(ranef(m.1)) %>%
  filter(grpvar == "Sample:ROI") %>%
  mutate(
    T_shift = condval,
    sd_T_shift = condsd,
    Sample = str_split(grp, pattern = ":", simplify = TRUE)[, 1],
    ROI = str_split(grp, pattern = ":", simplify = TRUE)[, 2]
  ) %>%
  dplyr::select(-c(condval, grpvar, term, condsd, grp))

Age.trend <- as_tibble(lstrends(m.1, ~ ROI, var = "Age")) %>%
  mutate(Age.trend_T = Age.trend) %>%
  dplyr::select(c(ROI, Age.trend_T))

dados_datasetscomp <- full_join(dados_datasetscomp, Age.trend) %>%
  full_join(re) %>%
  mutate(
    AvgThickness_shiftc = AvgThickness - T_shift,
    AvgThickness_age_decay = AvgThickness - Age.trend_T * (Age - Age.cor),
    AvgThickness_age_decay_shiftc = AvgThickness - T_shift - Age.trend_T *
      (Age - Age.cor),
    logAvgThickness_shiftc = log10(AvgThickness_shiftc),
    logAvgThickness_age_decay = log10(AvgThickness_age_decay),
    logAvgThickness_age_decay_shiftc = log10(AvgThickness_age_decay_shiftc)
  )

4.3 TotalArea

m.1 <-
  lme4::lmer(TotalArea_corrected ~ Age * ROI + (1 | Sample:ROI), data = dados_datasetscomp_rate)

re <- as_tibble(ranef(m.1)) %>%
  filter(grpvar == "Sample:ROI") %>%
  mutate(
    AT_shift = condval,
    sd_AT_shift = condsd,
    Sample = str_split(grp, pattern = ":", simplify = TRUE)[, 1],
    ROI = str_split(grp, pattern = ":", simplify = TRUE)[, 2]
  ) %>%
  dplyr::select(-c(condval, grpvar, term, condsd, grp))

Age.trend <-
  as_tibble(lstrends(m.1, ~ ROI, var = "Age")) %>%
  mutate(Age.trend_AT = Age.trend) %>%
  dplyr::select(c(ROI, Age.trend_AT))

dados_datasetscomp <- full_join(dados_datasetscomp, re) %>%
  full_join(Age.trend) %>%
  mutate(
    TotalArea_shiftc = TotalArea_corrected - AT_shift,
    TotalArea_age_decay = TotalArea_corrected - Age.trend_AT * (Age - Age.cor),
    TotalArea_age_decay_shiftc = TotalArea_corrected - AT_shift - Age.trend_AT *
      (Age - Age.cor),
    logTotalArea_shiftc = log10(TotalArea_shiftc),
    logTotalArea_age_decay = log10(TotalArea_age_decay),
    logTotalArea_age_decay_shiftc = log10(TotalArea_age_decay_shiftc)
  )

4.4 ExposedArea

m.1 <-
  lme4::lmer(
    ExposedArea_corrected  ~ Age * ROI + (1 | Sample:ROI), data = dados_datasetscomp_rate)

re <- as_tibble(ranef(m.1)) %>%
  filter(grpvar == "Sample:ROI") %>%
  mutate(
    AE_shift = condval,
    sd_AE_shift = condsd,
    Sample = str_split(grp, pattern = ":", simplify = TRUE)[, 1],
    ROI = str_split(grp, pattern = ":", simplify = TRUE)[, 2]
  ) %>%
  dplyr::select(-c(condval, grpvar, term, condsd, grp))

Age.trend <-
  as_tibble(lstrends(m.1, ~ ROI, var = "Age")) %>%
  mutate(Age.trend_AE = Age.trend) %>%
  dplyr::select(c(ROI, Age.trend_AE))

dados_datasetscomp <- full_join(dados_datasetscomp, re) %>%
  full_join(Age.trend) %>%
  mutate(
    ExposedArea_shiftc = ExposedArea_corrected - AE_shift,
    ExposedArea_age_decay = ExposedArea_corrected - Age.trend_AE * (Age - Age.cor),
    ExposedArea_age_decay_shiftc = ExposedArea_corrected - AE_shift - Age.trend_AE * (Age - Age.cor),
    logExposedArea_shiftc = log10(ExposedArea_shiftc),
    logExposedArea_age_decay = log10(ExposedArea_age_decay),
    logExposedArea_age_decay_shiftc = log10(ExposedArea_age_decay_shiftc)
  )

4.5 Lobes correction factor

m.1 <-
  lme4::lmer(c ~ Age * ROI + (1|Sample:ROI), data = dados_datasetscomp_rate)

re <- as_tibble(ranef(m.1)) %>%
  filter(grpvar == "Sample:ROI") %>%
  mutate(
    c_shift = condval,
    sd_c_shift = condsd,
    Sample = str_split(grp, pattern = ":", simplify = TRUE)[, 1],
    ROI = str_split(grp, pattern = ":", simplify = TRUE)[, 2]
  ) %>%
  dplyr::select(-c(condval, grpvar, term, condsd, grp))

Age.trend <- as_tibble(lstrends(m.1, ~ ROI, var = "Age")) %>%
  mutate(Age.trend_c = Age.trend) %>%
  dplyr::select(c(ROI, Age.trend_c))

dados_datasetscomp <- full_join(dados_datasetscomp, Age.trend) %>%
  full_join(re) %>%
  mutate(
    c_shiftc = c - c_shift,
    c_age_decay = c - Age.trend_c * (Age - Age.cor),
    c_age_decay_shiftc = c - c_shift - Age.trend_c *
      (Age - Age.cor)
  )
# ----
dados_datasetscomp <- dados_datasetscomp %>%
  mutate(
    localGI_age_decay = TotalArea_age_decay/ExposedArea_age_decay,
    localGI_shiftc = TotalArea_shiftc/ExposedArea_shiftc,
    localGI_age_decay_shiftc = TotalArea_age_decay_shiftc/ExposedArea_age_decay_shiftc,
    K_age_decay = log10(TotalArea_age_decay) + 1/4*log10(AvgThickness_age_decay^2) - 5/4*log10(ExposedArea_age_decay),
    K_shiftc = log10(TotalArea_shiftc) + 1/4*log10(AvgThickness_shiftc^2) - 5/4*log10(ExposedArea_shiftc),
    K_age_decay_shiftc = log10(TotalArea_age_decay_shiftc) + 1/4*log10(AvgThickness_age_decay_shiftc^2) - 5/4*log10(ExposedArea_age_decay_shiftc),
    I_age_decay = log10(TotalArea_age_decay) + log10(ExposedArea_age_decay) + log10(AvgThickness_age_decay^2),
    I_shiftc = log10(TotalArea_shiftc) + log10(ExposedArea_shiftc) + log10(AvgThickness_shiftc^2),
    I_age_decay_shiftc = log10(TotalArea_age_decay_shiftc) + log10(ExposedArea_age_decay_shiftc) + log10(AvgThickness_age_decay_shiftc^2),
    S_age_decay = 3/2*log10(TotalArea_age_decay) + 3/4*log10(ExposedArea_age_decay) - 9/4*log10(AvgThickness_age_decay^2),
    S_shiftc = 3/2*log10(TotalArea_shiftc) + 3/4*log10(ExposedArea_shiftc) - 9/4*log10(AvgThickness_shiftc^2),
    S_age_decay_shiftc = 3/2*log10(TotalArea_age_decay_shiftc) + 3/4*log10(ExposedArea_age_decay_shiftc) - 9/4*log10(AvgThickness_age_decay_shiftc^2),
    Knorm_shiftc = K_shiftc / sqrt(1 + (1 / 4) ^ 2 + (5 / 2) ^ 2),
    Snorm_shiftc = S_shiftc / sqrt((3 / 2) ^ 2 + (3 / 4) ^ 2 + (9 / 4) ^ 2),
    Inorm_shiftc = I_shiftc / sqrt(1 ^ 2 + 1 ^ 2 + 1 ^ 2),
    Knorm_age_decay = K_age_decay / sqrt(1 + (1 / 4) ^ 2 + (5 / 2) ^ 2),
    Snorm_age_decay = S_age_decay / sqrt((3 / 2) ^ 2 + (3 / 4) ^ 2 + (9 / 4) ^ 2),
    Inorm_age_decay = I_age_decay / sqrt(1 ^ 2 + 1 ^ 2 + 1 ^ 2),
    Knorm_age_decay_shiftc = K_age_decay_shiftc / sqrt(1 + (1 / 4) ^ 2 + (5 / 2) ^ 2),
    Snorm_age_decay_shiftc = S_age_decay_shiftc / sqrt((3 / 2) ^ 2 +  (3 / 4) ^ 2 + (9 / 4) ^ 2),
    Inorm_age_decay_shiftc = I_age_decay_shiftc / sqrt(1 ^ 2 + 1 ^ 2 + 1 ^ 2)
  )

5 RESULTS

5.1 Data description

Table 1
Sample Diagnostic N age age_range
ADNI CTL 868 75 ± 6.5 56 ; 96
ADNI AD 542 75 ± 8 56 ; 92
AHEAD CTL 100 42 ± 19 24 ; 76
AOMICPIOP1 CTL 208 22 ± 1.8 18 ; 26
AOMICPIOP2 CTL 224 22 ± 1.8 18 ; 26
HCPr900 CTL 881 29 ± 3.6 24 ; 37
IDOR CTL 77 66 ± 8.4 43 ; 80
IDOR AD 13 77 ± 6.1 63 ; 86
IXI-Guys CTL 314 51 ± 16 20 ; 86
IXI-HH CTL 181 47 ± 17 20 ; 82
IXI-IOP CTL 68 42 ± 17 20 ; 86
NKI CTL 168 34 ± 19 4 ; 85
OASIS CTL 312 45 ± 24 18 ; 94
Diagnostic N age age_range
CTL 3095 46 ± 23 4 ; 96
AD 555 75 ± 8 56 ; 92

5.2 Uncertainties estimation for baseline

5.2.1 K ~ Age

m.1 <- lme4::lmer(K ~ Age * Diagnostic +(1|Sample) + (1|Sample:Diagnostic) , data = filter(dados_datasetscomp, ROI == "hemisphere"))

Model summary and R squared:

summary(m.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: K ~ Age * Diagnostic + (1 | Sample) + (1 | Sample:Diagnostic)
##    Data: filter(dados_datasetscomp, ROI == "hemisphere")
## 
## REML criterion at convergence: -42555
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.9733 -0.6065  0.0381  0.6446  5.2725 
## 
## Random effects:
##  Groups            Name        Variance  Std.Dev.
##  Sample:Diagnostic (Intercept) 4.607e-06 0.002146
##  Sample            (Intercept) 4.410e-04 0.021001
##  Residual                      2.658e-04 0.016304
## Number of obs: 7912, groups:  Sample:Diagnostic, 13; Sample, 11
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)      -4.920e-01  6.412e-03  -76.73
## Age              -8.612e-04  1.675e-05  -51.40
## DiagnosticAD     -8.439e-02  5.444e-03  -15.50
## Age:DiagnosticAD  8.855e-04  6.366e-05   13.91
## 
## Correlation of Fixed Effects:
##             (Intr) Age    DgnsAD
## Age         -0.113              
## DiagnostcAD -0.039  0.221       
## Ag:DgnstcAD  0.030 -0.264 -0.882
r.squaredGLMM(m.1)
##            R2m       R2c
## [1,] 0.4372324 0.7897408

Natural fluctuation error:

signif(sigma.hat(m.1)$sigma$data,2)
## [1] 0.016

Type B error:

signif(sigma.hat(m.1)$sigma$Sample,2)
## (Intercept) 
##       0.021

5.2.2 S ~ Age

m.1 <- lme4::lmer(S ~ Age * Diagnostic +(1|Sample) + (1|Sample:Diagnostic) , data = filter(dados_datasetscomp, ROI == "hemisphere"))

Model summary and R squared:

summary(m.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S ~ Age * Diagnostic + (1 | Sample) + (1 | Sample:Diagnostic)
##    Data: filter(dados_datasetscomp, ROI == "hemisphere")
## 
## REML criterion at convergence: -10789.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5837 -0.6903 -0.0315  0.6406  5.5519 
## 
## Random effects:
##  Groups            Name        Variance  Std.Dev.
##  Sample:Diagnostic (Intercept) 0.0003168 0.01780 
##  Sample            (Intercept) 0.0057953 0.07613 
##  Residual                      0.0147855 0.12160 
## Number of obs: 7912, groups:  Sample:Diagnostic, 13; Sample, 11
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)       9.0840260  0.0242645 374.375
## Age               0.0017030  0.0001246  13.664
## DiagnosticAD      0.1969053  0.0413152   4.766
## Age:DiagnosticAD -0.0013879  0.0004747  -2.924
## 
## Correlation of Fixed Effects:
##             (Intr) Age    DgnsAD
## Age         -0.222              
## DiagnostcAD -0.079  0.212       
## Ag:DgnstcAD  0.058 -0.263 -0.866
r.squaredGLMM(m.1)
##            R2m      R2c
## [1,] 0.1515722 0.399721

Natural fluctuation error:

signif(sigma.hat(m.1)$sigma$data,2)
## [1] 0.12

Type B error:

signif(sigma.hat(m.1)$sigma$Sample,2)
## (Intercept) 
##       0.076

5.2.3 I ~ Age

m.1 <- lme4::lmer(I ~ Age * Diagnostic +(1|Sample) + (1|Sample:Diagnostic) , data = filter(dados_datasetscomp, ROI == "hemisphere"))

Model summary and R squared:

summary(m.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: I ~ Age * Diagnostic + (1 | Sample) + (1 | Sample:Diagnostic)
##    Data: filter(dados_datasetscomp, ROI == "hemisphere")
## 
## REML criterion at convergence: -17050
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.9306 -0.6593 -0.0005  0.6762  3.2113 
## 
## Random effects:
##  Groups            Name        Variance  Std.Dev.
##  Sample:Diagnostic (Intercept) 0.0001859 0.01364 
##  Sample            (Intercept) 0.0013945 0.03734 
##  Residual                      0.0067038 0.08188 
## Number of obs: 7912, groups:  Sample:Diagnostic, 13; Sample, 11
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)       1.053e+01  1.259e-02 836.069
## Age              -3.118e-03  8.371e-05 -37.254
## DiagnosticAD     -1.829e-01  2.848e-02  -6.421
## Age:DiagnosticAD  1.747e-03  3.196e-04   5.467
## 
## Correlation of Fixed Effects:
##             (Intr) Age    DgnsAD
## Age         -0.287              
## DiagnostcAD -0.109  0.200       
## Ag:DgnstcAD  0.075 -0.262 -0.844
r.squaredGLMM(m.1)
##            R2m       R2c
## [1,] 0.4505911 0.5554086

Natural fluctuation error:

signif(sigma.hat(m.1)$sigma$data,2)
## [1] 0.082

Type B error:

signif(sigma.hat(m.1)$sigma$Sample,2)
## (Intercept) 
##       0.037

5.3 Without outliers

filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL" | Diagnostic == "AD") %>%
  mutate(K_outlier = is_outlier(filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL" | Diagnostic == "AD")$K)) %>%
  filter(K_outlier == TRUE) %>%
  dplyr::select(c(SUBJ, Age, Gender, Diagnostic, Sample)) %>% unique()
## # A tibble: 0 x 5
## # ... with 5 variables: SUBJ <chr>, Age <dbl>, Gender <chr>, Diagnostic <fct>,
## #   Sample <chr>
outlier_subj <- filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL" | Diagnostic == "AD") %>%
  mutate(K_outlier = is_outlier(filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL" | Diagnostic == "AD")$K)) %>%
  filter(K_outlier == TRUE) %>%
  dplyr::select(c(SUBJ, Age, Gender, Diagnostic, Sample)) %>% unique()

paste("N outliers = ", n_distinct((filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL" | Diagnostic == "AD") %>%
  mutate(K_outlier = is_outlier(filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL" | Diagnostic == "AD")$K)) %>%
  filter(K_outlier == TRUE))$SUBJ))
## [1] "N outliers =  0"
dados_datasetscomp_rate <- filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL" | Diagnostic == "AD")  %>%
  anti_join(outlier_subj)

dados_datasetscomp_rate$Diagnostic <- factor(dados_datasetscomp_rate$Diagnostic, levels = c("CTL","AD"))
dados_datasetscomp_rate$Sample <- as.factor(dados_datasetscomp_rate$Sample)

5.3.1 K ~ Age

m.1 <- lme4::lmer(K ~ Age * Diagnostic +(1|Sample) + (1|Sample:Diagnostic) , data = filter(dados_datasetscomp_rate, ROI == "hemisphere"))

Model summary and R squared:

summary(m.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: K ~ Age * Diagnostic + (1 | Sample) + (1 | Sample:Diagnostic)
##    Data: filter(dados_datasetscomp_rate, ROI == "hemisphere")
## 
## REML criterion at convergence: -42555
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.9733 -0.6065  0.0381  0.6446  5.2725 
## 
## Random effects:
##  Groups            Name        Variance  Std.Dev.
##  Sample:Diagnostic (Intercept) 4.607e-06 0.002146
##  Sample            (Intercept) 4.410e-04 0.021001
##  Residual                      2.658e-04 0.016304
## Number of obs: 7912, groups:  Sample:Diagnostic, 13; Sample, 11
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)      -4.920e-01  6.412e-03  -76.73
## Age              -8.612e-04  1.675e-05  -51.40
## DiagnosticAD     -8.439e-02  5.444e-03  -15.50
## Age:DiagnosticAD  8.855e-04  6.366e-05   13.91
## 
## Correlation of Fixed Effects:
##             (Intr) Age    DgnsAD
## Age         -0.113              
## DiagnostcAD -0.039  0.221       
## Ag:DgnstcAD  0.030 -0.264 -0.882
r.squaredGLMM(m.1)
##            R2m       R2c
## [1,] 0.4372324 0.7897408

Natural fluctuation error:

signif(sigma.hat(m.1)$sigma$data,2)
## [1] 0.016

Type B error:

signif(sigma.hat(m.1)$sigma$Sample,2)
## (Intercept) 
##       0.021

5.3.2 S ~ Age

m.1 <- lme4::lmer(S ~ Age * Diagnostic +(1|Sample) + (1|Sample:Diagnostic) , data = filter(dados_datasetscomp_rate, ROI == "hemisphere"))

Model summary and R squared:

summary(m.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S ~ Age * Diagnostic + (1 | Sample) + (1 | Sample:Diagnostic)
##    Data: filter(dados_datasetscomp_rate, ROI == "hemisphere")
## 
## REML criterion at convergence: -10789.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5837 -0.6903 -0.0315  0.6406  5.5519 
## 
## Random effects:
##  Groups            Name        Variance  Std.Dev.
##  Sample:Diagnostic (Intercept) 0.0003168 0.01780 
##  Sample            (Intercept) 0.0057953 0.07613 
##  Residual                      0.0147855 0.12160 
## Number of obs: 7912, groups:  Sample:Diagnostic, 13; Sample, 11
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)       9.0840260  0.0242645 374.375
## Age               0.0017030  0.0001246  13.664
## DiagnosticAD      0.1969053  0.0413152   4.766
## Age:DiagnosticAD -0.0013879  0.0004747  -2.924
## 
## Correlation of Fixed Effects:
##             (Intr) Age    DgnsAD
## Age         -0.222              
## DiagnostcAD -0.079  0.212       
## Ag:DgnstcAD  0.058 -0.263 -0.866
r.squaredGLMM(m.1)
##            R2m      R2c
## [1,] 0.1515722 0.399721

Natural fluctuation error:

signif(sigma.hat(m.1)$sigma$data,2)
## [1] 0.12

Type B error:

signif(sigma.hat(m.1)$sigma$Sample,2)
## (Intercept) 
##       0.076

5.3.3 I ~ Age

m.1 <- lme4::lmer(I ~ Age * Diagnostic +(1|Sample) + (1|Sample:Diagnostic) , data = filter(dados_datasetscomp_rate, ROI == "hemisphere"))

Model summary and R squared:

summary(m.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: I ~ Age * Diagnostic + (1 | Sample) + (1 | Sample:Diagnostic)
##    Data: filter(dados_datasetscomp_rate, ROI == "hemisphere")
## 
## REML criterion at convergence: -17050
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.9306 -0.6593 -0.0005  0.6762  3.2113 
## 
## Random effects:
##  Groups            Name        Variance  Std.Dev.
##  Sample:Diagnostic (Intercept) 0.0001859 0.01364 
##  Sample            (Intercept) 0.0013945 0.03734 
##  Residual                      0.0067038 0.08188 
## Number of obs: 7912, groups:  Sample:Diagnostic, 13; Sample, 11
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)       1.053e+01  1.259e-02 836.069
## Age              -3.118e-03  8.371e-05 -37.254
## DiagnosticAD     -1.829e-01  2.848e-02  -6.421
## Age:DiagnosticAD  1.747e-03  3.196e-04   5.467
## 
## Correlation of Fixed Effects:
##             (Intr) Age    DgnsAD
## Age         -0.287              
## DiagnostcAD -0.109  0.200       
## Ag:DgnstcAD  0.075 -0.262 -0.844
r.squaredGLMM(m.1)
##            R2m       R2c
## [1,] 0.4505911 0.5554086

Natural fluctuation error:

signif(sigma.hat(m.1)$sigma$data,2)
## [1] 0.082

Type B error:

signif(sigma.hat(m.1)$sigma$Sample,2)
## (Intercept) 
##       0.037