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 (fernandahmoraes@gmail.com) in case any help is needed, our aim is to improve the code as needed!
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")
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")
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")
# 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")
# 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"))
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)
)
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)
)
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)
)
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)
)
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)
)
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 |
figS6a <- ggplot(filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL"),
aes(x = Age , y = K)) +
geom_point(aes(color = Sample, fill = Sample, alpha = 0.4)) +
geom_smooth(color = "black", method = "lm") +
# stat_cor() +
theme_pubr() +
guides(alpha = "none")+
labs(x = "Age [years]") +
scale_x_continuous(limits = c(0,100))
cor.test(filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL")$K, filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL")$Age)
##
## Pearson's product-moment correlation
##
## data: filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL")$K and filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL")$Age
## t = -100.83, df = 6800, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.7834543 -0.7643992
## sample estimates:
## cor
## -0.7741021
figS6b <- ggplot(filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL"),
aes(x = Age, y = S)) +
geom_point(aes(color = Sample, fill = Sample, alpha = 0.4)) +
geom_smooth(color = "black", method = "lm") +
# stat_cor() +
theme_pubr() +
guides(alpha = "none", color = FALSE,fill = FALSE)+
labs(x = "Age [years]") +
theme(legend.position= "none") +
scale_x_continuous(limits = c(0,100))
cor.test(filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL")$S, filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL")$Age)
##
## Pearson's product-moment correlation
##
## data: filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL")$S and filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL")$Age
## t = 37.406, df = 6800, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3931982 0.4326213
## sample estimates:
## cor
## 0.4131033
figS6c <- ggplot(filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL"),
aes(x = Age, y = I)) +
geom_point(aes(color = Sample, fill = Sample, alpha = 0.4)) +
geom_smooth(color = "black", method = "lm") +
# stat_cor() +
theme_pubr() +
guides(alpha = "none", color = FALSE,fill = FALSE)+
labs(x = "Age [years]") +
theme(legend.position= "none") +
scale_x_continuous(limits = c(0,100))
cor.test(filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL")$I, filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL")$Age)
##
## Pearson's product-moment correlation
##
## data: filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL")$I and filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL")$Age
## t = -70.698, df = 6800, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6643648 -0.6369635
## sample estimates:
## cor
## -0.6508761
figS6 <- ggarrange(figS6a, ggarrange(figS6b, figS6c,labels = c("B","C"), nrow = 1, ncol = 2, font.label = list(size = 11)),labels = c("A"), nrow = 2, ncol = 1, font.label = list(size = 11), common.legend = TRUE, legend = "top")
ggsave("figS6.pdf", plot = figS6, dpi=1200, width = 18, height = 22, units = "cm", device = "pdf")
FIGURE 1
figS6
Through age, every Healthy Control subject for the independent morphological component, K, S, and I. For HCPr900 and AHEAD, ages are determined by an interval of years, thereby considering the mean age. The solid line represents a linear regression with a 95% confidence interval.
summary(lm(
1 / 2 * log10(AvgThickness) + log10(TotalArea) ~ log10(ExposedArea),
data = filter(dados_datasetscomp,
Sample != "Mota&Houzel2015",
ROI == "hemisphere"),
na.action = na.omit
))
##
## Call:
## lm(formula = 1/2 * log10(AvgThickness) + log10(TotalArea) ~ log10(ExposedArea),
## data = filter(dados_datasetscomp, Sample != "Mota&Houzel2015",
## ROI == "hemisphere"), na.action = na.omit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.147267 -0.026113 -0.001921 0.029974 0.105636
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.19727 0.05515 -21.71 <2e-16 ***
## log10(ExposedArea) 1.39352 0.01201 115.99 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03537 on 7910 degrees of freedom
## Multiple R-squared: 0.6298, Adjusted R-squared: 0.6297
## F-statistic: 1.345e+04 on 1 and 7910 DF, p-value: < 2.2e-16
summary(lm(
1 / 2 * log10(AvgThickness) + log10(TotalArea) ~ log10(ExposedArea),
data = filter(dados_datasetscomp,
Sample != "Mota&Houzel2015",
ROI == "hemisphere", Age > 18 & Age < 40),
na.action = na.omit
))
##
## Call:
## lm(formula = 1/2 * log10(AvgThickness) + log10(TotalArea) ~ log10(ExposedArea),
## data = filter(dados_datasetscomp, Sample != "Mota&Houzel2015",
## ROI == "hemisphere", Age > 18 & Age < 40), na.action = na.omit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.075674 -0.012016 0.002633 0.013807 0.065150
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.55119 0.05033 -10.95 <2e-16 ***
## log10(ExposedArea) 1.25912 0.01095 115.02 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02063 on 3534 degrees of freedom
## Multiple R-squared: 0.7892, Adjusted R-squared: 0.7891
## F-statistic: 1.323e+04 on 1 and 3534 DF, p-value: < 2.2e-16
## Displays confidence interval
# tidy(lm(
# 1 / 2 * log10(AvgThickness) + log10(TotalArea) ~ log10(ExposedArea),
# data = dados_hemi_v1,
# na.action = na.omit
# ), conf.int = TRUE)
paste(
"Student's t test comapring slope with theoretical value 1.25. t = ",
signif(abs(coef(summary(
lm(
1 / 2 * log10(AvgThickness) + log10(TotalArea) ~ log10(ExposedArea),
data = filter(dados_datasetscomp,
Sample != "Mota&Houzel2015",
ROI == "hemisphere"),
na.action = na.omit
)
))[2, 1] - 1.25) / coef(summary(
lm(
1 / 2 * log10(AvgThickness) + log10(TotalArea) ~ log10(ExposedArea),
data = filter(dados_datasetscomp,
Sample != "Mota&Houzel2015",
ROI == "hemisphere"),
na.action = na.omit
)
))[2, 2], 3)
)
## [1] "Student's t test comapring slope with theoretical value 1.25. t = 11.9"
paste(
"Student's t test comapring slope with theoretical value 1.25. p value = ",
signif(test_coef(
lm(
1 / 2 * log10(AvgThickness) + log10(TotalArea) ~ log10(ExposedArea),
data = filter(dados_datasetscomp,
Sample != "Mota&Houzel2015",
ROI == "hemisphere"),
na.action = na.omit
),
2,
1.25
),3)
)
## [1] "Student's t test comapring slope with theoretical value 1.25. p value = 1.3e-32"
lm_Age <-
filter(
dados_datasetscomp,
Sample != "Mota&Houzel2015",
ROI == "hemisphere",
Diagnostic == "CTL",
Age_interval != "[0,5)",
Age_interval != "[5,10)",
Age_interval != "[10,15)"
) %>%
group_by(Age_interval) %>%
do(fit_Age = tidy(
lm(
1 / 2 * log10(AvgThickness) + log10(TotalArea) ~ log10(ExposedArea),
data = .,
na.action = na.omit
),
conf.int = TRUE
)) %>%
unnest(cols = c(fit_Age))
lm_Age <- lm_Age %>% mutate(Age_interval = as.double((str_sub(Age_interval,2,3))))
cor.test(filter(lm_Age, term == "log10(ExposedArea)")$estimate, filter(lm_Age, term == "log10(ExposedArea)")$Age_interval)
##
## Pearson's product-moment correlation
##
## data: filter(lm_Age, term == "log10(ExposedArea)")$estimate and filter(lm_Age, term == "log10(ExposedArea)")$Age_interval
## t = -4.2916, df = 15, p-value = 0.0006426
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9013924 -0.4069511
## sample estimates:
## cor
## -0.742386
Part of Table 2 - Natural fluctuation and Type B errors
m.1 <- lme4::lmer(GMvolume ~ 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: GMvolume ~ Age * Diagnostic + (1 | Sample) + (1 | Sample:Diagnostic)
## Data: filter(dados_datasetscomp, ROI == "hemisphere")
##
## REML criterion at convergence: 183968.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8395 -0.6727 -0.0320 0.6298 4.5351
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample:Diagnostic (Intercept) 27722513 5265
## Sample (Intercept) 366049150 19132
## Residual 733577651 27085
## Number of obs: 7912, groups: Sample:Diagnostic, 13; Sample, 11
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 303365.74 6119.08 49.577
## Age -1061.86 27.79 -38.216
## DiagnosticAD -67730.51 9855.71 -6.872
## Age:DiagnosticAD 725.19 105.75 6.858
##
## Correlation of Fixed Effects:
## (Intr) Age DgnsAD
## Age -0.196
## DiagnostcAD -0.083 0.195
## Ag:DgnstcAD 0.052 -0.263 -0.809
r.squaredGLMM(m.1)
## R2m R2c
## [1,] 0.3957362 0.6067994
Mean value:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
as_tibble(mean)[1,2]
## # A tibble: 1 x 1
## lsmean
## <dbl>
## 1 250218.
Natural fluctuation error:
signif(sigma.hat(m.1)$sigma$data,2)
## [1] 27000
Type B error:
signif(sigma.hat(m.1)$sigma$Sample,2)
## (Intercept)
## 19000
m.1 <- lme4::lmer(AvgThickness ~ 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:
## AvgThickness ~ Age * Diagnostic + (1 | Sample) + (1 | Sample:Diagnostic)
## Data: filter(dados_datasetscomp, ROI == "hemisphere")
##
## REML criterion at convergence: -13111.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.8038 -0.6167 0.0104 0.6495 3.4140
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample:Diagnostic (Intercept) 0.0006026 0.02455
## Sample (Intercept) 0.0071880 0.08478
## Residual 0.0110146 0.10495
## Number of obs: 7912, groups: Sample:Diagnostic, 13; Sample, 11
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.7038037 0.0270734 99.869
## Age -0.0044384 0.0001077 -41.194
## DiagnosticAD -0.3421007 0.0405587 -8.435
## Age:DiagnosticAD 0.0029937 0.0004098 7.305
##
## Correlation of Fixed Effects:
## (Intr) Age DgnsAD
## Age -0.172
## DiagnostcAD -0.084 0.183
## Ag:DgnstcAD 0.045 -0.263 -0.762
r.squaredGLMM(m.1)
## R2m R2c
## [1,] 0.4615061 0.6845916
Mean value:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
as_tibble(mean)[1,2]
## # A tibble: 1 x 1
## lsmean
## <dbl>
## 1 2.48
Natural fluctuation error:
signif(sigma.hat(m.1)$sigma$data,2)
## [1] 0.1
Type B error:
signif(sigma.hat(m.1)$sigma$Sample,2)
## (Intercept)
## 0.085
m.1 <- lme4::lmer(TotalArea ~ 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: TotalArea ~ Age * Diagnostic + (1 | Sample) + (1 | Sample:Diagnostic)
## Data: filter(dados_datasetscomp, ROI == "hemisphere")
##
## REML criterion at convergence: 167642
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9066 -0.7162 -0.0476 0.6633 4.2777
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample:Diagnostic (Intercept) 0 0
## Sample (Intercept) 23918246 4891
## Residual 93180143 9653
## Number of obs: 7912, groups: Sample:Diagnostic, 13; Sample, 11
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 113025.061 1543.436 73.229
## Age -244.189 9.866 -24.750
## DiagnosticAD -13184.422 2852.379 -4.622
## Age:DiagnosticAD 152.485 37.675 4.047
##
## Correlation of Fixed Effects:
## (Intr) Age DgnsAD
## Age -0.277
## DiagnostcAD -0.073 0.259
## Ag:DgnstcAD 0.072 -0.263 -0.992
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
r.squaredGLMM(m.1)
## R2m R2c
## [1,] 0.2353489 0.3915348
Mean value:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
as_tibble(mean)[1,2]
## # A tibble: 1 x 1
## lsmean
## <dbl>
## 1 100803.
Natural fluctuation error:
signif(sigma.hat(m.1)$sigma$data,2)
## [1] 9700
Type B error:
signif(sigma.hat(m.1)$sigma$Sample,2)
## (Intercept)
## 4900
m.1 <- lme4::lmer(ExposedArea ~ 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:
## ExposedArea ~ Age * Diagnostic + (1 | Sample) + (1 | Sample:Diagnostic)
## Data: filter(dados_datasetscomp, ROI == "hemisphere")
##
## REML criterion at convergence: 148576.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8868 -0.7093 -0.0618 0.6706 3.8432
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample:Diagnostic (Intercept) 0 0.0
## Sample (Intercept) 279141 528.3
## Residual 8382585 2895.3
## Number of obs: 7912, groups: Sample:Diagnostic, 13; Sample, 11
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 40714.770 207.022 196.669
## Age -40.602 2.859 -14.204
## DiagnosticAD -262.931 853.269 -0.308
## Age:DiagnosticAD 1.535 11.275 0.136
##
## Correlation of Fixed Effects:
## (Intr) Age DgnsAD
## Age -0.597
## DiagnostcAD -0.151 0.249
## Ag:DgnstcAD 0.152 -0.255 -0.992
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
r.squaredGLMM(m.1)
## R2m R2c
## [1,] 0.09982676 0.1288366
Mean value:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
as_tibble(mean)[1,2]
## # A tibble: 1 x 1
## lsmean
## <dbl>
## 1 38683.
Natural fluctuation error:
signif(sigma.hat(m.1)$sigma$data,2)
## [1] 0.095
Type B error:
signif(sigma.hat(m.1)$sigma$Sample,2)
## (Intercept)
## 0.11
m.1 <- lme4::lmer(localGI ~ 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
Mean value:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
as_tibble(mean)[1,2]
## # A tibble: 1 x 1
## lsmean
## <dbl>
## 1 2.60
Natural fluctuation error:
signif(sigma.hat(m.1)$sigma$data,2)
## [1] 0.095
Type B error:
signif(sigma.hat(m.1)$sigma$Sample,2)
## (Intercept)
## 0.11
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
Mean value:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
as_tibble(mean)[1,2]
## # A tibble: 1 x 1
## lsmean
## <dbl>
## 1 -0.535
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
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
Mean value:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
as_tibble(mean)[1,2]
## # A tibble: 1 x 1
## lsmean
## <dbl>
## 1 9.17
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
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
Mean value:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
as_tibble(mean)[1,2]
## # A tibble: 1 x 1
## lsmean
## <dbl>
## 1 10.4
Natural fluctuation error:
signif(sigma.hat(m.1)$sigma$data,2)
## [1] 27000
Type B error:
signif(sigma.hat(m.1)$sigma$Sample,2)
## (Intercept)
## 0
FIGURE 2
K_lifespan
K through age for Healthy Controls (A) raw data, (B) after harmonizing (removing the estimated residual from the Linear Mixed Model) from T, AT, and AE resulting in the harmonized values of K.
# dados_datasetscomp$Diagnostic <- as.factor(dados_datasetscomp$Diagnostic)
# dados_datasetscomp$Diagnostic <- relevel(dados_datasetscomp$Diagnostic, ref = "CTL")
summary(lm(
1 / 2 * log10(AvgThickness_shiftc) + log10(TotalArea_shiftc) ~ log10(ExposedArea_shiftc),
data = filter(dados_datasetscomp,
Sample != "Mota&Houzel2015",
ROI == "hemisphere"),
na.action = na.omit
))
##
## Call:
## lm(formula = 1/2 * log10(AvgThickness_shiftc) + log10(TotalArea_shiftc) ~
## log10(ExposedArea_shiftc), data = filter(dados_datasetscomp,
## Sample != "Mota&Houzel2015", ROI == "hemisphere"), na.action = na.omit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.135371 -0.017810 0.003617 0.020761 0.100026
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.392350 0.042559 -32.72 <2e-16 ***
## log10(ExposedArea_shiftc) 1.436420 0.009281 154.78 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02835 on 7910 degrees of freedom
## Multiple R-squared: 0.7518, Adjusted R-squared: 0.7517
## F-statistic: 2.396e+04 on 1 and 7910 DF, p-value: < 2.2e-16
## Displays confidence interval
# tidy(lm(
# 1 / 2 * log10(AvgThickness) + log10(TotalArea) ~ log10(ExposedArea),
# data = dados_hemi_v1,
# na.action = na.omit
# ), conf.int = TRUE)
lm_Age <-
filter(
dados_datasetscomp,
Sample != "Mota&Houzel2015",
ROI == "hemisphere",
Diagnostic == "CTL",
Age_interval != "[0,5)",
Age_interval != "[5,10)",
Age_interval != "[10,15)"
) %>%
group_by(Age_interval) %>%
do(fit_Age = tidy(
lm(
1 / 2 * log10(AvgThickness_shiftc) + log10(TotalArea_shiftc) ~ log10(ExposedArea_shiftc),
data = .,
na.action = na.omit
),
conf.int = TRUE
)) %>%
unnest(cols = c(fit_Age))
lm_Age %>%
kable(digits = 2) %>%
kable_styling()
Age_interval | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
---|---|---|---|---|---|---|---|
[15,20) | (Intercept) | -0.89 | 0.19 | -4.76 | 0.00 | -1.26 | -0.52 |
[15,20) | log10(ExposedArea_shiftc) | 1.33 | 0.04 | 32.78 | 0.00 | 1.25 | 1.41 |
[20,25) | (Intercept) | -0.59 | 0.05 | -11.25 | 0.00 | -0.70 | -0.49 |
[20,25) | log10(ExposedArea_shiftc) | 1.27 | 0.01 | 110.57 | 0.00 | 1.25 | 1.29 |
[25,30) | (Intercept) | -0.51 | 0.06 | -8.70 | 0.00 | -0.62 | -0.39 |
[25,30) | log10(ExposedArea_shiftc) | 1.25 | 0.01 | 98.73 | 0.00 | 1.22 | 1.27 |
[30,35) | (Intercept) | -0.51 | 0.07 | -7.65 | 0.00 | -0.64 | -0.38 |
[30,35) | log10(ExposedArea_shiftc) | 1.25 | 0.01 | 86.34 | 0.00 | 1.22 | 1.28 |
[35,40) | (Intercept) | -0.69 | 0.14 | -4.94 | 0.00 | -0.97 | -0.41 |
[35,40) | log10(ExposedArea_shiftc) | 1.29 | 0.03 | 42.30 | 0.00 | 1.23 | 1.35 |
[40,45) | (Intercept) | -0.14 | 0.15 | -0.94 | 0.35 | -0.44 | 0.16 |
[40,45) | log10(ExposedArea_shiftc) | 1.17 | 0.03 | 35.56 | 0.00 | 1.10 | 1.23 |
[45,50) | (Intercept) | -0.26 | 0.14 | -1.81 | 0.07 | -0.54 | 0.02 |
[45,50) | log10(ExposedArea_shiftc) | 1.19 | 0.03 | 38.09 | 0.00 | 1.13 | 1.25 |
[50,55) | (Intercept) | -0.48 | 0.14 | -3.56 | 0.00 | -0.75 | -0.21 |
[50,55) | log10(ExposedArea_shiftc) | 1.24 | 0.03 | 41.84 | 0.00 | 1.18 | 1.30 |
[55,60) | (Intercept) | -0.25 | 0.14 | -1.79 | 0.07 | -0.53 | 0.03 |
[55,60) | log10(ExposedArea_shiftc) | 1.19 | 0.03 | 38.40 | 0.00 | 1.13 | 1.25 |
[60,65) | (Intercept) | 0.13 | 0.10 | 1.28 | 0.20 | -0.07 | 0.32 |
[60,65) | log10(ExposedArea_shiftc) | 1.10 | 0.02 | 51.31 | 0.00 | 1.06 | 1.15 |
[65,70) | (Intercept) | 0.20 | 0.12 | 1.66 | 0.10 | -0.04 | 0.44 |
[65,70) | log10(ExposedArea_shiftc) | 1.09 | 0.03 | 41.05 | 0.00 | 1.03 | 1.14 |
[70,75) | (Intercept) | -0.27 | 0.10 | -2.74 | 0.01 | -0.46 | -0.08 |
[70,75) | log10(ExposedArea_shiftc) | 1.19 | 0.02 | 55.27 | 0.00 | 1.15 | 1.23 |
[75,80) | (Intercept) | 0.30 | 0.10 | 3.02 | 0.00 | 0.11 | 0.50 |
[75,80) | log10(ExposedArea_shiftc) | 1.06 | 0.02 | 48.67 | 0.00 | 1.02 | 1.10 |
[80,85) | (Intercept) | -0.43 | 0.14 | -3.00 | 0.00 | -0.71 | -0.15 |
[80,85) | log10(ExposedArea_shiftc) | 1.22 | 0.03 | 39.02 | 0.00 | 1.16 | 1.28 |
[85,90) | (Intercept) | 0.25 | 0.21 | 1.18 | 0.24 | -0.17 | 0.67 |
[85,90) | log10(ExposedArea_shiftc) | 1.07 | 0.05 | 23.26 | 0.00 | 0.98 | 1.16 |
[90,95) | (Intercept) | 0.19 | 0.38 | 0.48 | 0.63 | -0.60 | 0.98 |
[90,95) | log10(ExposedArea_shiftc) | 1.09 | 0.08 | 12.89 | 0.00 | 0.91 | 1.26 |
[95,100] | (Intercept) | 0.86 | 4.17 | 0.21 | 0.86 | -17.08 | 18.80 |
[95,100] | log10(ExposedArea_shiftc) | 0.93 | 0.91 | 1.02 | 0.41 | -3.00 | 4.86 |
lm_Age <- lm_Age %>% mutate(Age_interval = as.double((str_sub(Age_interval,2,3))))
cor.test(filter(lm_Age, term == "log10(ExposedArea_shiftc)")$estimate, filter(lm_Age, term == "log10(ExposedArea_shiftc)")$Age_interval)
##
## Pearson's product-moment correlation
##
## data: filter(lm_Age, term == "log10(ExposedArea_shiftc)")$estimate and filter(lm_Age, term == "log10(ExposedArea_shiftc)")$Age_interval
## t = -5.7812, df = 15, p-value = 3.624e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9372002 -0.5829227
## sample estimates:
## cor
## -0.8307961
amostra_Coef <- filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL", Age_interval != "[0,5)", Age_interval != "[95,100]") %>%
group_by(Age_interval) %>%
do(fit_amostra = tidy(lm(1/2 * log10(AvgThickness_shiftc) + log10(TotalArea_shiftc) ~ log10(ExposedArea_shiftc), data = ., na.action = na.omit), conf.int = TRUE, conf.level = 0.95)) %>%
unnest(fit_amostra)
amostras_Coef_age <- filter(dados_datasetscomp, ROI == "hemisphere", Diagnostic == "CTL", Age_interval != "[0,5)", Age_interval != "[95,100]") %>%
group_by(Age_interval) %>%
summarise(
N = n_distinct(SUBJ),
min_age = min(Age),
max_age = max(Age))
amostra_Coef <- full_join(amostra_Coef, amostras_Coef_age) %>% filter(term == "log10(ExposedArea_shiftc)")
fig_slope_ageinterval <- ggplot(amostra_Coef, aes(y = estimate, x = Age_interval)) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
geom_point() +
geom_line(group = 1) +
geom_hline(yintercept = 1.25, linetype = "dashed") +
theme_pubr() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size=9)
) + labs(y ="Slope (after harmonization)", x = "Age interval") + geom_text(aes(label = N), nudge_y = 0.2)
ggsave("fig_slope_ageinterval.pdf", plot = fig_slope_ageinterval, dpi=1200, width = 18, height = 14, units = "cm", device = "pdf")
FIGURE 3
fig_slope_ageinterval
Slope alpha derived from the Cortical Folding model from Mota & Houzel through age for Healthy Controls. Data points from (0, 5] and [95, 100] years old were omitted due to the reduced data. The numbers on top of each point indicates the number of subjects at each interval. Mean value for all ages were 1.38+-0.012 and correlation of alpha and Age is Pearson’s r=-0.69, p=0.0031 from [15,20) years old.
m.1 <- lme4::lmer(GMvolume_shiftc ~ Age * Diagnostic +(1|Sample) + (1|Sample:Diagnostic) , data = filter(dados_datasetscomp, ROI == "hemisphere"))
summary(m.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## GMvolume_shiftc ~ Age * Diagnostic + (1 | Sample) + (1 | Sample:Diagnostic)
## Data: filter(dados_datasetscomp, ROI == "hemisphere")
##
## REML criterion at convergence: 183907.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8431 -0.6738 -0.0306 0.6297 4.5365
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample:Diagnostic (Intercept) 0.00e+00 0
## Sample (Intercept) 0.00e+00 0
## Residual 7.33e+08 27074
## Number of obs: 7912, groups: Sample:Diagnostic, 13; Sample, 11
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 303346.92 740.16 409.837
## Age -1062.07 14.43 -73.583
## DiagnosticAD -72189.07 7734.48 -9.333
## Age:DiagnosticAD 732.24 102.94 7.113
##
## Correlation of Fixed Effects:
## (Intr) Age DgnsAD
## Age -0.896
## DiagnostcAD -0.096 0.086
## Ag:DgnstcAD 0.126 -0.140 -0.992
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
m.lstk <- lstrends(m.1, ~ Diagnostic, var="Age")
pairs <- pairs(m.lstk)
mean <- as_tibble(mean)
m.lstk <- as_tibble(m.lstk)
lst <- m.lstk %>%
mutate(TX = paste(signif(Age.trend, 2), "±", signif(SE, 2)))
tableGM <- full_join(mean, lst, by = c("Diagnostic")) %>%
mutate(percent = Age.trend*100/abs(lsmean), Variable = "GM_shiftc")
coefs <- data.frame(coef(summary(m.1)))
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
coefs
## Estimate Std..Error t.value p.z
## (Intercept) 303346.9208 740.16472 409.837041 0.000000e+00
## Age -1062.0694 14.43361 -73.583075 0.000000e+00
## DiagnosticAD -72189.0707 7734.47570 -9.333415 0.000000e+00
## Age:DiagnosticAD 732.2418 102.94214 7.113140 1.134426e-12
pairs %>%
kable() %>%
kable_styling()
contrast | estimate | SE | df | z.ratio | p.value |
---|---|---|---|---|---|
CTL - AD | -732.2418 | 102.9421 | Inf | -7.11314 | 0 |
ggplot(data = as_tibble(pairs), aes(
x = contrast,
y = estimate,
ymin = estimate - SE,
ymax = estimate + SE)) +
geom_hline(yintercept = 0,
linetype = "11",
colour = "grey60") +
geom_pointrange( position = position_dodge(width = 0.2)) +
geom_text(aes(label = str_c("p adj ", signif(p.value, digits = 2))), nudge_x = 0.3, nudge_y = 0.0003) +
coord_flip() +
labs(y = expression('Differences in slopes ('*Delta*'GM/'*Delta*'Age)'), x = "Contrast") +
theme_pubr() +
theme(axis.title = element_text(size = 11),
axis.text = element_text(size = 10), text = element_text(size = 10))
e <- effect("Age:Diagnostic", m.1)
e <- as.data.frame(e)
ggplot(e, aes(x = Age, y = fit,ymin=lower, ymax=upper, shape = Diagnostic)) + geom_pointrange() + geom_line() + theme_pubr() + labs(x = "Age [years]", y = "K", shape = "Diagnostic") +
theme(axis.title = element_text(size = 11),
axis.text = element_text(size = 10), text = element_text(size = 10))
Mean value:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
as_tibble(mean)[1,2]
## # A tibble: 1 x 1
## lsmean
## <dbl>
## 1 250189.
Natural fluctuation error:
signif(sigma.hat(m.1)$sigma$data,2)
## [1] 27000
Type B error:
signif(sigma.hat(m.1)$sigma$Sample,2)
## (Intercept)
## 0
m.1 <- lme4::lmer(AvgThickness_shiftc ~ Age * Diagnostic +(1|Sample) + (1|Sample:Diagnostic) , data = filter(dados_datasetscomp, ROI == "hemisphere"))
summary(m.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: AvgThickness_shiftc ~ Age * Diagnostic + (1 | Sample) + (1 |
## Sample:Diagnostic)
## Data: filter(dados_datasetscomp, ROI == "hemisphere")
##
## REML criterion at convergence: -13173.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.8059 -0.6167 0.0103 0.6500 3.4058
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample:Diagnostic (Intercept) 0.00000 0.0000
## Sample (Intercept) 0.00000 0.0000
## Residual 0.01101 0.1049
## Number of obs: 7912, groups: Sample:Diagnostic, 13; Sample, 11
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.704e+00 2.869e-03 942.562
## Age -4.440e-03 5.594e-05 -79.376
## DiagnosticAD -3.695e-01 2.998e-02 -12.325
## Age:DiagnosticAD 3.029e-03 3.990e-04 7.592
##
## Correlation of Fixed Effects:
## (Intr) Age DgnsAD
## Age -0.896
## DiagnostcAD -0.096 0.086
## Ag:DgnstcAD 0.126 -0.140 -0.992
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
m.lstk <- lstrends(m.1, ~ Diagnostic, var="Age")
mean <- as_tibble(mean)
m.lstk <- as_tibble(m.lstk)
lst <- m.lstk %>%
mutate(TX = paste(signif(Age.trend, 2), "±", signif(SE, 2)))
tableT <- full_join(mean, lst, by = c("Diagnostic")) %>%
mutate(percent = Age.trend*100/abs(lsmean), Variable = "T_shiftc")
Mean value:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
as_tibble(mean)[1,2]
## # A tibble: 1 x 1
## lsmean
## <dbl>
## 1 2.48
Natural fluctuation error:
signif(sigma.hat(m.1)$sigma$data,2)
## [1] 0.1
Type B error:
signif(sigma.hat(m.1)$sigma$Sample,2)
## (Intercept)
## 0
m.1 <- lme4::lmer(TotalArea_shiftc ~ Age * Diagnostic +(1|Sample) + (1|Sample:Diagnostic) , data = filter(dados_datasetscomp, ROI == "hemisphere"))
summary(m.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## TotalArea_shiftc ~ Age * Diagnostic + (1 | Sample) + (1 | Sample:Diagnostic)
## Data: filter(dados_datasetscomp, ROI == "hemisphere")
##
## REML criterion at convergence: 167586.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9088 -0.7169 -0.0470 0.6641 4.2805
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample:Diagnostic (Intercept) 0 0
## Sample (Intercept) 0 0
## Residual 93066204 9647
## Number of obs: 7912, groups: Sample:Diagnostic, 13; Sample, 11
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 113025.783 263.740 428.550
## Age -244.375 5.143 -47.515
## DiagnosticAD -13212.654 2755.994 -4.794
## Age:DiagnosticAD 152.779 36.681 4.165
##
## Correlation of Fixed Effects:
## (Intr) Age DgnsAD
## Age -0.896
## DiagnostcAD -0.096 0.086
## Ag:DgnstcAD 0.126 -0.140 -0.992
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
sigma.hat(m.1)
## $sigma
## $sigma$data
## [1] 9647.083
##
## $sigma$`Sample:Diagnostic`
## (Intercept)
## 0
##
## $sigma$Sample
## (Intercept)
## 0
##
##
## $cors
## $cors$data
## [1] NA
##
## $cors$`Sample:Diagnostic`
## [1] NA
##
## $cors$Sample
## [1] NA
r.squaredGLMM(m.1)
## R2m R2c
## [1,] 0.2795186 0.2795186
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
m.lstk <- lstrends(m.1, ~ Diagnostic, var="Age")
mean <- as_tibble(mean)
m.lstk <- as_tibble(m.lstk)
lst <- m.lstk %>%
mutate(
TX = paste(signif(Age.trend, 2), "±", signif(SE, 2)))
tableAT <- full_join(mean, lst, by = c("Diagnostic")) %>%
mutate(percent = Age.trend*100/abs(lsmean), Variable = "AT_shiftc")
Mean value:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
as_tibble(mean)[1,2]
## # A tibble: 1 x 1
## lsmean
## <dbl>
## 1 100795.
Natural fluctuation error:
signif(sigma.hat(m.1)$sigma$data,2)
## [1] 2900
Type B error:
signif(sigma.hat(m.1)$sigma$Sample,2)
## (Intercept)
## 6e-07
m.1 <- lme4::lmer(ExposedArea_shiftc ~ Age * Diagnostic +(1|Sample) + (1|Sample:Diagnostic) , data = filter(dados_datasetscomp, ROI == "hemisphere"))
summary(m.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## ExposedArea_shiftc ~ Age * Diagnostic + (1 | Sample) + (1 | Sample:Diagnostic)
## Data: filter(dados_datasetscomp, ROI == "hemisphere")
##
## REML criterion at convergence: 148540.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8881 -0.7093 -0.0605 0.6726 3.8448
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample:Diagnostic (Intercept) 1.753e-13 4.187e-07
## Sample (Intercept) 3.654e-13 6.045e-07
## Residual 8.372e+06 2.893e+03
## Number of obs: 7912, groups: Sample:Diagnostic, 13; Sample, 11
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 40729.774 79.104 514.890
## Age -41.214 1.543 -26.718
## DiagnosticAD -321.787 826.609 -0.389
## Age:DiagnosticAD 2.208 11.002 0.201
##
## Correlation of Fixed Effects:
## (Intr) Age DgnsAD
## Age -0.896
## DiagnostcAD -0.096 0.086
## Ag:DgnstcAD 0.126 -0.140 -0.992
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
m.lstk <- lstrends(m.1, ~ Diagnostic, var="Age")
mean <- as_tibble(mean)
m.lstk <- as_tibble(m.lstk)
lst <- m.lstk %>% mutate(
TX = paste(signif(Age.trend, 2), "±", signif(SE, 2)))
tableAE <- full_join(mean, lst, by = c("Diagnostic")) %>% mutate(percent = Age.trend*100/abs(lsmean), Variable = "AE_shiftc")
Mean value:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
as_tibble(mean)[1,2]
## # A tibble: 1 x 1
## lsmean
## <dbl>
## 1 38667.
Natural fluctuation error:
signif(sigma.hat(m.1)$sigma$data,2)
## [1] 2900
Type B error:
signif(sigma.hat(m.1)$sigma$Sample,2)
## (Intercept)
## 6e-07
m.1 <- lme4::lmer(localGI_shiftc ~ Age * Diagnostic +(1|Sample) + (1|Sample:Diagnostic) , data = filter(dados_datasetscomp, ROI == "hemisphere"))
summary(m.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## localGI_shiftc ~ Age * Diagnostic + (1 | Sample) + (1 | Sample:Diagnostic)
## Data: filter(dados_datasetscomp, ROI == "hemisphere")
##
## REML criterion at convergence: -14560.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7703 -0.6798 -0.0048 0.6759 6.4327
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample:Diagnostic (Intercept) 0.000000 0.00000
## Sample (Intercept) 0.000000 0.00000
## Residual 0.009238 0.09611
## Number of obs: 7912, groups: Sample:Diagnostic, 13; Sample, 11
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.778e+00 2.628e-03 1057.14
## Age -3.509e-03 5.124e-05 -68.49
## DiagnosticAD -3.243e-01 2.746e-02 -11.81
## Age:DiagnosticAD 3.811e-03 3.655e-04 10.43
##
## Correlation of Fixed Effects:
## (Intr) Age DgnsAD
## Age -0.896
## DiagnostcAD -0.096 0.086
## Ag:DgnstcAD 0.126 -0.140 -0.992
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
m.lstk <- lstrends(m.1, ~ Diagnostic, var="Age")
mean <- as_tibble(mean)
m.lstk <- as_tibble(m.lstk)
lst <- m.lstk %>%
mutate(
TX = paste(signif(Age.trend, 2), "±", signif(SE, 2)))
tableGI <- full_join(mean, lst, by = c("Diagnostic")) %>%
mutate(percent = Age.trend*100/abs(lsmean), Variable = "localGI_shiftc")
Mean value:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
as_tibble(mean)[1,2]
## # A tibble: 1 x 1
## lsmean
## <dbl>
## 1 2.60
Natural fluctuation error:
signif(sigma.hat(m.1)$sigma$data,2)
## [1] 0.096
Type B error:
signif(sigma.hat(m.1)$sigma$Sample,2)
## (Intercept)
## 0
m.1 <- lme4::lmer(K_shiftc ~ Age * Diagnostic +(1|Sample) + (1|Sample:Diagnostic) , data = filter(dados_datasetscomp, ROI == "hemisphere"))
summary(m.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: K_shiftc ~ Age * Diagnostic + (1 | Sample) + (1 | Sample:Diagnostic)
## Data: filter(dados_datasetscomp, ROI == "hemisphere")
##
## REML criterion at convergence: -42446
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.5322 -0.6041 0.0391 0.6431 5.2736
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample:Diagnostic (Intercept) 0.0000000 0.00000
## Sample (Intercept) 0.0000000 0.00000
## Residual 0.0002717 0.01648
## Number of obs: 7912, groups: Sample:Diagnostic, 13; Sample, 11
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -4.916e-01 4.507e-04 -1090.92
## Age -8.600e-04 8.788e-06 -97.86
## DiagnosticAD -8.706e-02 4.709e-03 -18.49
## Age:DiagnosticAD 8.940e-04 6.268e-05 14.26
##
## Correlation of Fixed Effects:
## (Intr) Age DgnsAD
## Age -0.896
## DiagnostcAD -0.096 0.086
## Ag:DgnstcAD 0.126 -0.140 -0.992
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
coefs <- data.frame(coef(summary(m.1)))
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
coefs
## Estimate Std..Error t.value p.z
## (Intercept) -0.4916413892 4.506654e-04 -1090.92322 0
## Age -0.0008599766 8.788218e-06 -97.85563 0
## DiagnosticAD -0.0870584281 4.709304e-03 -18.48647 0
## Age:DiagnosticAD 0.0008939612 6.267857e-05 14.26263 0
Mean value:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
as_tibble(mean)[1,2]
## # A tibble: 1 x 1
## lsmean
## <dbl>
## 1 -0.535
Mean values:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
pairs_means <- pairs(mean)
mean <- as_tibble(mean)
pairs_means %>%
kable() %>%
kable_styling()
contrast | estimate | SE | df | z.ratio | p.value |
---|---|---|---|---|---|
CTL - AD | -0.1403122 | 0.0121857 | Inf | -11.51445 | 0 |
paste("Difference in pct ", as_tibble(pairs_means)[2]*100/mean$lsmean[1])
## [1] "Difference in pct -1.53044891669301"
Trends:
m.lstk <- lstrends(m.1, ~ Diagnostic, var="Age")
pairs_trends <- pairs(m.lstk)
m.lstk <- as_tibble(m.lstk)
pairs_trends %>%
kable() %>%
kable_styling()
contrast | estimate | SE | df | z.ratio | p.value |
---|---|---|---|---|---|
CTL - AD | -0.000894 | 6.27e-05 | Inf | -14.26263 | 0 |
lst <- m.lstk %>%
mutate(
TX = paste(signif(Age.trend, 2), "±", signif(SE, 2)))
tableK <- full_join(mean, lst, by = c("Diagnostic")) %>%
mutate(percent = Age.trend*100/abs(lsmean), Variable = "K_shiftc")
Figures:
fig_TX_K_age_c <- ggplot(data = as_tibble(pairs), aes(
x = contrast,
y = estimate,
ymin = estimate - SE,
ymax = estimate + SE)) +
geom_hline(yintercept = 0,
linetype = "11",
colour = "grey60") +
geom_pointrange( position = position_dodge(width = 0.2)) +
geom_text(aes(label = str_c("p adj ", signif(p.value, digits = 2))), nudge_x = 0.3, nudge_y = 0.0003) +
coord_flip() +
labs(y = expression('Differences in slopes ('*Delta*'K/'*Delta*'Age)'), x = "Contrast") +
theme_pubr() +
theme(axis.title = element_text(size = 11),
axis.text = element_text(size = 10), text = element_text(size = 10))
e <- effect("Age:Diagnostic", m.1)
e <- as.data.frame(e)
fig_K_age_c <- ggplot(e, aes(x = Age, y = fit,ymin=lower, ymax=upper, shape = Diagnostic)) + geom_pointrange() + geom_line() + theme_pubr() + labs(x = "Age [years]", y = "K", shape = "Diagnostic") +
theme(axis.title = element_text(size = 11),
axis.text = element_text(size = 10), text = element_text(size = 10))
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
m.1 <- lme4::lmer(S_shiftc ~ Age * Diagnostic +(1|Sample) + (1|Sample:Diagnostic) , data = filter(dados_datasetscomp, ROI == "hemisphere"))
summary(m.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_shiftc ~ Age * Diagnostic + (1 | Sample) + (1 | Sample:Diagnostic)
## Data: filter(dados_datasetscomp, ROI == "hemisphere")
##
## REML criterion at convergence: -10771.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4726 -0.6906 -0.0334 0.6362 5.8254
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample:Diagnostic (Intercept) 0.00000 0.0000
## Sample (Intercept) 0.00000 0.0000
## Residual 0.01492 0.1221
## Number of obs: 7912, groups: Sample:Diagnostic, 13; Sample, 11
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 9.088e+00 3.339e-03 2721.724
## Age 1.606e-03 6.511e-05 24.671
## DiagnosticAD 2.073e-01 3.489e-02 5.943
## Age:DiagnosticAD -1.339e-03 4.644e-04 -2.884
##
## Correlation of Fixed Effects:
## (Intr) Age DgnsAD
## Age -0.896
## DiagnostcAD -0.096 0.086
## Ag:DgnstcAD 0.126 -0.140 -0.992
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
coefs <- data.frame(coef(summary(m.1)))
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
coefs
## Estimate Std..Error t.value p.z
## (Intercept) 9.087638690 3.338928e-03 2721.723676 0.000000e+00
## Age 0.001606360 6.511088e-05 24.671142 0.000000e+00
## DiagnosticAD 0.207341356 3.489069e-02 5.942599 2.805383e-09
## Age:DiagnosticAD -0.001339213 4.643782e-04 -2.883885 3.928021e-03
Mean value:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
as_tibble(mean)[1,2]
## # A tibble: 1 x 1
## lsmean
## <dbl>
## 1 9.17
Mean values:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
pairs_means <- pairs(mean)
mean <- as_tibble(mean)
pairs_means %>%
kable() %>%
kable_styling()
contrast | estimate | SE | df | z.ratio | p.value |
---|---|---|---|---|---|
CTL - AD | -0.1403122 | 0.0121857 | Inf | -11.51445 | 0 |
paste("Difference in pct ", as_tibble(pairs_means)[2]*100/mean$lsmean[1])
## [1] "Difference in pct -1.53044891669301"
Trends:
m.lstk <- lstrends(m.1, ~ Diagnostic, var="Age")
pairs_trends <- pairs(m.lstk)
m.lstk <- as_tibble(m.lstk)
pairs_trends %>%
kable() %>%
kable_styling()
contrast | estimate | SE | df | z.ratio | p.value |
---|---|---|---|---|---|
CTL - AD | 0.0013392 | 0.0004644 | Inf | 2.883885 | 0.003928 |
lst <- m.lstk %>% mutate(
TX = paste(signif(Age.trend, 2), "±", signif(SE, 2)))
tableS <- full_join(mean, lst, by = c("Diagnostic")) %>%
mutate(percent = Age.trend*100/abs(lsmean), Variable = "S_shiftc")
Figures:
fig_TX_S_age_c <- ggplot(data = as_tibble(pairs), aes(
x = contrast,
y = estimate,
ymin = estimate - SE,
ymax = estimate + SE)) +
geom_hline(yintercept = 0,
linetype = "11",
colour = "grey60") +
geom_pointrange( position = position_dodge(width = 0.2)) +
geom_text(aes(label = str_c("p adj ", signif(p.value, digits = 2))), nudge_x = 0.3, nudge_y = 0.0003) +
coord_flip() +
labs(y = expression('Differences in slopes ('*Delta*'S/'*Delta*'Age)'), x = "Contrast") +
theme_pubr() +
theme(axis.title = element_text(size = 11),
axis.text = element_text(size = 10), text = element_text(size = 10))
e <- effect("Age:Diagnostic", m.1)
e <- as.data.frame(e)
fig_S_age_c <- ggplot(e, aes(x = Age, y = fit,ymin=lower, ymax=upper, shape = Diagnostic)) + geom_pointrange() + geom_line() + theme_pubr() + labs(x = "Age [years]", y = "S", shape = "Diagnostic") +
theme(axis.title = element_text(size = 11),
axis.text = element_text(size = 10),
text = element_text(size = 10), legend.position = "none")
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
m.1 <- lme4::lmer(I_shiftc ~ Age * Diagnostic +(1|Sample) + (1|Sample:Diagnostic) , data = filter(dados_datasetscomp, ROI == "hemisphere"))
summary(m.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: I_shiftc ~ Age * Diagnostic + (1 | Sample) + (1 | Sample:Diagnostic)
## Data: filter(dados_datasetscomp, ROI == "hemisphere")
##
## REML criterion at convergence: -16982.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.7921 -0.6605 0.0035 0.6753 3.2170
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample:Diagnostic (Intercept) 0.000000 0.00000
## Sample (Intercept) 0.000000 0.00000
## Residual 0.006801 0.08247
## Number of obs: 7912, groups: Sample:Diagnostic, 13; Sample, 11
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.053e+01 2.255e-03 4669.930
## Age -3.074e-03 4.397e-05 -69.929
## DiagnosticAD -1.938e-01 2.356e-02 -8.226
## Age:DiagnosticAD 1.702e-03 3.136e-04 5.426
##
## Correlation of Fixed Effects:
## (Intr) Age DgnsAD
## Age -0.896
## DiagnostcAD -0.096 0.086
## Ag:DgnstcAD 0.126 -0.140 -0.992
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
coefs <- data.frame(coef(summary(m.1)))
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
coefs
## Estimate Std..Error t.value p.z
## (Intercept) 10.528806334 2.254596e-03 4669.930202 0.000000e+00
## Age -0.003074501 4.396583e-05 -69.929319 0.000000e+00
## DiagnosticAD -0.193800114 2.355978e-02 -8.225888 2.220446e-16
## Age:DiagnosticAD 0.001701526 3.135693e-04 5.426315 5.752943e-08
Mean value:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
as_tibble(mean)[1,2]
## # A tibble: 1 x 1
## lsmean
## <dbl>
## 1 10.4
Mean values:
mean <- lsmeans(m.1, ~ Diagnostic, var="Age")
pairs_means <- pairs(mean)
mean <- as_tibble(mean)
pairs_means %>%
kable() %>%
kable_styling()
contrast | estimate | SE | df | z.ratio | p.value |
---|---|---|---|---|---|
CTL - AD | 0.1086367 | 0.0082284 | Inf | 13.20271 | 0 |
paste("Difference in pct ", as_tibble(pairs_means)[2]*100/mean$lsmean[1])
## [1] "Difference in pct 1.04710882243449"
Trends:
m.lstk <- lstrends(m.1, ~ Diagnostic, var="Age")
pairs_trends <- pairs(m.lstk)
m.lstk <- as_tibble(m.lstk)
pairs_trends %>%
kable() %>%
kable_styling()
contrast | estimate | SE | df | z.ratio | p.value |
---|---|---|---|---|---|
CTL - AD | -0.0017015 | 0.0003136 | Inf | -5.426315 | 1e-07 |
lst <- m.lstk %>% mutate(
TX = paste(signif(Age.trend, 2), "±", signif(SE, 2)))
tableI <- full_join(mean, lst, by = c("Diagnostic")) %>%
mutate(percent = Age.trend*100/abs(lsmean), Variable = "I_shiftc")
Figures:
fig_TX_I_age_c <- ggplot(data = as_tibble(pairs), aes(
x = contrast,
y = estimate,
ymin = estimate - SE,
ymax = estimate + SE)) +
geom_hline(yintercept = 0,
linetype = "11",
colour = "grey60") +
geom_pointrange( position = position_dodge(width = 0.2)) +
geom_text(aes(label = str_c("p adj ", signif(p.value, digits = 2))), nudge_x = 0.3, nudge_y = 0.0003) +
coord_flip() +
labs(y = expression('Differences in slopes ('*Delta*'I/'*Delta*'Age)'), x = "Contrast") +
theme_pubr() +
theme(axis.title = element_text(size = 11),
axis.text = element_text(size = 10), text = element_text(size = 10))
e <- effect("Age:Diagnostic", m.1)
e <- as.data.frame(e)
fig_I_age_c <- ggplot(e, aes(x = Age, y = fit,ymin=lower, ymax=upper, shape = Diagnostic)) + geom_pointrange() + geom_line() + theme_pubr() + labs(x = "Age [years]", y = "I", shape = "Diagnostic") +
theme(axis.title = element_text(size = 11),
axis.text = element_text(size = 10),
text = element_text(size = 10), legend.position = "none")
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
TABLE 3
rate %>%
kable(digits = 2) %>%
kable_styling()
Variable | CTL | AD |
---|---|---|
T_shiftc | -0.0044 ± 5.6e-05 ( -0.18 ) | -0.0014 ± 4e-04 ( -0.062 ) |
AT_shiftc | -240 ± 5.1 ( -0.24 ) | -92 ± 36 ( -0.096 ) |
AE_shiftc | -41 ± 1.5 ( -0.11 ) | -39 ± 11 ( -0.1 ) |
K_shiftc | -0.00086 ± 8.8e-06 ( -0.0094 ) | 3.4e-05 ± 6.2e-05 ( 0.00037 ) |
S_shiftc | 0.0016 ± 6.5e-05 ( 0.018 ) | 0.00027 ± 0.00046 ( 0.0029 ) |
I_shiftc | -0.0031 ± 4.4e-05 ( -0.03 ) | -0.0014 ± 0.00031 ( -0.013 ) |
localGI_shiftc | -0.0035 ± 5.1e-05 ( -0.13 ) | 3e-04 ± 0.00036 ( 0.012 ) |
GM_shiftc | -1100 ± 14 ( -0.42 ) | -330 ± 100 ( -0.15 ) |
fig_age_c <- ggarrange(fig_K_age_c, ggarrange(fig_S_age_c, fig_I_age_c,labels = c("B","C"), nrow = 1, ncol = 2, font.label = list(size = 11)),labels = c("A"), nrow = 2, ncol = 1, font.label = list(size = 11), common.legend = TRUE, legend = "top")
ggsave("fig_age_c.pdf", plot = fig_age_c, dpi=1200, width = 18, height = 18, units = "cm", device = "pdf")
FIGURE4
fig_age_c
Fitted values extracted from the linear regression model after data harmonization. Bars represent the 95 confidence interval. Complete summary of linear models in Supplementary Material. (A) Concerning K, Alzheimer’s Disease (AD) has a shallow slope, meaning small changes with Age. Compared to healthy controls, AD values of K are almost constant and similar to older subjects (AD slope p-value<0.0001 and CTL slope p-value<0.0001; pairwise comparison estimate -0.000903, p<0.0001). (B) For S, the AD and CTL patterns have similar intercepts, but statistically different slopes (AD slope p-value=0.004 and CTL slope p-value<0.0001; pairwise comparison estimate 0.0013, p=0.004). (C) For I, that reflects brain volume, CTL has a decreasing volume with aging, while AD has a smaller slope (AD slope p-value<0.0001 and CTL slope p-value<0.0001; pairwise comparison estimate -0.00168, p<0.0001).
dados_datasetscomp_lobes <- filter(dados_datasetscomp, ROI != "hemisphere", Diagnostic == "CTL" | Diagnostic == "AD", Sample != "IDOR-CCD-Control")
dados_datasetscomp_lobes$Diagnostic <- factor(dados_datasetscomp_lobes$Diagnostic, levels = c("CTL", "MCI","AD"))
dados_datasetscomp_lobes$Sample <- as.factor(dados_datasetscomp_lobes$Sample)
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 |
## Linear mixed model fit by REML ['lmerMod']
## Formula: K_shiftc ~ Age * Diagnostic * ROI + (1 | Sample:Diagnostic:ROI)
## Data: dados_datasetscomp_lobes
##
## REML criterion at convergence: -139166.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.3529 -0.6169 0.0203 0.6414 5.3896
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample:Diagnostic:ROI (Intercept) 0.0000000 0.00000
## Residual 0.0004934 0.02221
## Number of obs: 29188, groups: Sample:Diagnostic:ROI, 48
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -5.062e-01 6.403e-04 -790.494
## Age -9.198e-04 1.247e-05 -73.752
## DiagnosticAD -9.008e-02 6.343e-03 -14.201
## ROIO 4.222e-02 9.055e-04 46.621
## ROIP 4.204e-02 9.055e-04 46.429
## ROIT 2.639e-02 9.055e-04 29.142
## Age:DiagnosticAD 9.744e-04 8.447e-05 11.536
## Age:ROIO 2.355e-04 1.764e-05 13.352
## Age:ROIP -2.060e-04 1.764e-05 -11.680
## Age:ROIT -3.828e-05 1.764e-05 -2.170
## DiagnosticAD:ROIO 3.362e-02 8.970e-03 3.748
## DiagnosticAD:ROIP -4.009e-02 8.970e-03 -4.469
## DiagnosticAD:ROIT 6.636e-03 8.970e-03 0.740
## Age:DiagnosticAD:ROIO -4.070e-04 1.195e-04 -3.407
## Age:DiagnosticAD:ROIP 4.587e-04 1.195e-04 3.840
## Age:DiagnosticAD:ROIT -2.046e-04 1.195e-04 -1.713
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
Trends:
pairs %>%
kable() %>%
kable_styling()
contrast | estimate | SE | df | z.ratio | p.value | Contrast.1 | Contrast.2 | Diagnostic.1 | ROI.1 | Diagnostic.2 | ROI.2 | contrast_ROI |
---|---|---|---|---|---|---|---|---|---|---|---|---|
CTL-AD | -0.0009744 | 8.45e-05 | Inf | -11.535618 | 0 | CTL F | AD F | CTL | F | AD | F | F-F |
CTL-AD | -0.0005674 | 8.45e-05 | Inf | -6.716977 | 0 | CTL O | AD O | CTL | O | AD | O | O-O |
CTL-AD | -0.0014331 | 8.45e-05 | Inf | -16.966067 | 0 | CTL P | AD P | CTL | P | AD | P | P-P |
CTL-AD | -0.0007698 | 8.45e-05 | Inf | -9.113069 | 0 | CTL T | AD T | CTL | T | AD | T | T-T |
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_shiftc ~ Age * Diagnostic * ROI + (1 | Sample:Diagnostic:ROI)
## Data: dados_datasetscomp_lobes
##
## REML criterion at convergence: -29380.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2192 -0.6718 -0.0240 0.6368 5.8473
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample:Diagnostic:ROI (Intercept) 0.00000 0.0000
## Residual 0.02126 0.1458
## Number of obs: 29188, groups: Sample:Diagnostic:ROI, 48
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 8.731e+00 4.203e-03 2076.999
## Age 2.681e-03 8.187e-05 32.752
## DiagnosticAD 2.785e-01 4.164e-02 6.689
## ROIO -2.149e-01 5.945e-03 -36.156
## ROIP 6.513e-01 5.945e-03 109.566
## ROIT -4.981e-02 5.945e-03 -8.379
## Age:DiagnosticAD -2.415e-03 5.545e-04 -4.355
## Age:ROIO -1.212e-03 1.158e-04 -10.464
## Age:ROIP -1.198e-04 1.158e-04 -1.035
## Age:ROIT -5.236e-04 1.158e-04 -4.522
## DiagnosticAD:ROIO -1.154e-01 5.889e-02 -1.960
## DiagnosticAD:ROIP 2.092e-01 5.889e-02 3.552
## DiagnosticAD:ROIT -2.767e-02 5.889e-02 -0.470
## Age:DiagnosticAD:ROIO 9.338e-04 7.842e-04 1.191
## Age:DiagnosticAD:ROIP -2.442e-03 7.842e-04 -3.114
## Age:DiagnosticAD:ROIT 9.578e-04 7.842e-04 1.221
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
Trends:
pairs %>%
kable() %>%
kable_styling()
contrast | estimate | SE | df | z.ratio | p.value | Contrast.1 | Contrast.2 | Diagnostic.1 | ROI.1 | Diagnostic.2 | ROI.2 | contrast_ROI |
---|---|---|---|---|---|---|---|---|---|---|---|---|
CTL-AD | 0.0024148 | 0.0005545 | Inf | 4.354852 | 0.0003550 | CTL F | AD F | CTL | F | AD | F | F-F |
CTL-AD | 0.0014810 | 0.0005545 | Inf | 2.670780 | 0.1315406 | CTL O | AD O | CTL | O | AD | O | O-O |
CTL-AD | 0.0048567 | 0.0005545 | Inf | 8.758395 | 0.0000000 | CTL P | AD P | CTL | P | AD | P | P-P |
CTL-AD | 0.0014570 | 0.0005545 | Inf | 2.627543 | 0.1459513 | CTL T | AD T | CTL | T | AD | T | T-T |
## Linear mixed model fit by REML ['lmerMod']
## Formula: I_shiftc ~ Age * Diagnostic * ROI + (1 | Sample:Diagnostic:ROI)
## Data: dados_datasetscomp_lobes
##
## REML criterion at convergence: -47911.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5924 -0.6391 0.0025 0.6538 4.0292
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample:Diagnostic:ROI (Intercept) 0.00000 0.0000
## Residual 0.01127 0.1061
## Number of obs: 29188, groups: Sample:Diagnostic:ROI, 48
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.029e+01 3.060e-03 3363.538
## Age -3.634e-03 5.959e-05 -60.975
## DiagnosticAD -2.590e-01 3.031e-02 -8.546
## ROIO -7.612e-01 4.327e-03 -175.911
## ROIP 3.175e-01 4.327e-03 73.373
## ROIT 7.542e-02 4.327e-03 17.429
## Age:DiagnosticAD 2.767e-03 4.036e-04 6.855
## Age:ROIO 7.472e-04 8.428e-05 8.866
## Age:ROIP 5.915e-04 8.428e-05 7.018
## Age:ROIT 6.682e-04 8.428e-05 7.928
## DiagnosticAD:ROIO 1.045e-01 4.286e-02 2.437
## DiagnosticAD:ROIP -2.042e-03 4.286e-02 -0.048
## DiagnosticAD:ROIT 8.302e-03 4.286e-02 0.194
## Age:DiagnosticAD:ROIO -1.215e-03 5.708e-04 -2.129
## Age:DiagnosticAD:ROIP 5.997e-05 5.708e-04 0.105
## Age:DiagnosticAD:ROIT -7.663e-04 5.708e-04 -1.342
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
Trends:
pairs %>%
kable() %>%
kable_styling()
contrast | estimate | SE | df | z.ratio | p.value | Contrast.1 | Contrast.2 | Diagnostic.1 | ROI.1 | Diagnostic.2 | ROI.2 | contrast_ROI |
---|---|---|---|---|---|---|---|---|---|---|---|---|
CTL-AD | -0.0027668 | 0.0004036 | Inf | -6.854835 | 0.0000000 | CTL F | AD F | CTL | F | AD | F | F-F |
CTL-AD | -0.0015513 | 0.0004036 | Inf | -3.843424 | 0.0030554 | CTL O | AD O | CTL | O | AD | O | O-O |
CTL-AD | -0.0028267 | 0.0004036 | Inf | -7.003407 | 0.0000000 | CTL P | AD P | CTL | P | AD | P | P-P |
CTL-AD | -0.0020005 | 0.0004036 | Inf | -4.956275 | 0.0000198 | CTL T | AD T | CTL | T | AD | T | T-T |