ResultsSect.Rmd

---
title: "ResultsSect"
author: "Freya Acar"
# header-includes:
output: pdf_document
editor_options: 
  chunk_output_type: inline
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(cache=TRUE)
```

## Descriptives
### Participants
Thirty CW (mean age = 25.80 years, SD = 3.84), 30 CM (mean age = 26.03 years, SD = 5.26), 40 TM (mean age = 24.38 years, SD = 5.35), and 41 TW (mean age = 24.88 years, SD = 6.20) participated in the study. One TW participant was excluded because no results could be extracted from FreeSurfer. Demographics can be observed in Table 1. The sample did not differ significantly in age [F(3, 136) = 0.74, p = 0.528].


### Table with demographic information
```{r read data, warning = FALSE, echo = FALSE}
# Read in data
  data.all <- read.csv("../1.Data/Behzad_all.csv", sep=";", dec=",")
  data.hyp <- read.csv("../1.Data/Behzad_hyp.csv", sep=";", dec=",")
```

```{r table1, warning = FALSE, echo = FALSE}
library(knitr)

dig <- 2 # number of digits after comma
groupn <- c("Group", "CW", "CM", "TM", "TW")
ager <- c("Age", paste(round(mean(data.all[data.all[,2]==1,4]), digits = dig), " ± ", 
                round(sd(data.all[data.all[,2]==1,4]), digits = dig), sep=""),
          paste(round(mean(data.all[data.all[,2]==2,4]), digits = dig), " ± ", 
                round(sd(data.all[data.all[,2]==2,4]), digits = dig), sep=""),
          paste(round(mean(data.all[data.all[,2]==3,4]), digits = dig), " ± ", 
                round(sd(data.all[data.all[,2]==3,4]), digits = dig), sep=""),
          paste(round(mean(data.all[data.all[,2]==4,4]), digits = dig), " ± ", 
                round(sd(data.all[data.all[,2]==4,4]), digits = dig), sep=""))
SESr <- c("SES", paste(round(mean(data.all[data.all[,2]==1,5]), digits = dig), " ± ", 
                round(sd(data.all[data.all[,2]==1,5]), digits = dig), sep=""),
          paste(round(mean(data.all[data.all[,2]==2,5]), digits = dig), " ± ", 
                round(sd(data.all[data.all[,2]==2,5]), digits = dig), sep=""),
          paste(round(mean(data.all[data.all[,2]==3,5]), digits = dig), " ± ", 
                round(sd(data.all[data.all[,2]==3,5]), digits = dig), sep=""),
          paste(round(mean(data.all[data.all[,2]==4,5]), digits = dig), " ± ", 
                round(sd(data.all[data.all[,2]==4,5]), digits = dig), sep=""))
Edr <- c("Education", paste(round(mean(data.all[data.all[,2]==1,6]), digits = dig), " ± ", 
                round(sd(data.all[data.all[,2]==1,6]), digits = dig), sep=""),
          paste(round(mean(data.all[data.all[,2]==2,6]), digits = dig), " ± ", 
                round(sd(data.all[data.all[,2]==2,6]), digits = dig), sep=""),
          paste(round(mean(data.all[data.all[,2]==3,6]), digits = dig), " ± ", 
                round(sd(data.all[data.all[,2]==3,6]), digits = dig), sep=""),
          paste(round(mean(data.all[data.all[,2]==4,6]), digits = dig), " ± ", 
                round(sd(data.all[data.all[,2]==4,6]), digits = dig), sep=""))
Hdr <- c("Handedness", paste(round(mean(data.all[data.all[,2]==1,7]), digits = dig), " ± ", 
                round(sd(data.all[data.all[,2]==1,7]), digits = dig), sep=""),
          paste(round(mean(data.all[data.all[,2]==2,7]), digits = dig), " ± ", 
                round(sd(data.all[data.all[,2]==2,7]), digits = dig), sep=""),
          paste(round(mean(data.all[data.all[,2]==3,7]), digits = dig), " ± ", 
                round(sd(data.all[data.all[,2]==3,7]), digits = dig), sep=""),
          paste(round(mean(data.all[data.all[,2]==4,7]), digits = dig), " ± ", 
                round(sd(data.all[data.all[,2]==4,7]), digits = dig), sep=""))
kable(rbind(groupn, ager, SESr, Edr, Hdr), row.names = FALSE, caption = "Demographics")
```

### Code for demographics
The code can be found in markdown version of this file, it is not printed in the PDF
```{r demogr, warning = FALSE, echo = FALSE, include = FALSE}
# Groups
  ng1 <- sum(data.all[,2]==1)
  ng2 <- sum(data.all[,2]==2)
  ng3 <- sum(data.all[,2]==3)
  ng4 <- sum(data.all[,2]==4)
  
  paste("Cisgender women = ", ng1, ", cisgender men = ", ng2, 
        ", transgender men = ", ng3, ", transgender women = ",ng4, sep="")

# Age
  # Cisgender women
  mean(data.all[data.all[,2]==1,4])
  sd(data.all[data.all[,2]==1,4])
  
  # Cisgender men
  mean(data.all[data.all[,2]==2,4])
  sd(data.all[data.all[,2]==2,4])
  
  # Transgender men
  mean(data.all[data.all[,2]==3,4])
  sd(data.all[data.all[,2]==3,4])
  
  # Transgender women
  mean(data.all[data.all[,2]==4,4])
  sd(data.all[data.all[,2]==4,4])
  
  # Difference between groups?
  fit <- lm(formula = data.all[,4] ~ as.factor(data.all[,2]))
  anova(fit)
  
# Intracranial Volume
  # Cisgender women
  mean(c(data.all[data.all[,2]==1,108], data.all[data.all[,2]==1,182]))
  sd(c(data.all[data.all[,2]==1,108], data.all[data.all[,2]==1,182]))
  
  # Cisgender men
  mean(c(data.all[data.all[,2]==2,108], data.all[data.all[,2]==2,182]))
  sd(c(data.all[data.all[,2]==2,108], data.all[data.all[,2]==2,182]))
  
  # Transgender men
  mean(c(data.all[data.all[,2]==3,108], data.all[data.all[,2]==3,182]))
  sd(c(data.all[data.all[,2]==3,108], data.all[data.all[,2]==3,182]))
  
  # Transgender women
  mean(c(data.all[data.all[,2]==4,108], data.all[data.all[,2]==4,182]))
  sd(c(data.all[data.all[,2]==4,108], data.all[data.all[,2]==4,182]))
  
  # Difference between groups?
  fit <- lm(formula = c(data.all[,108], data.all[,182]) ~ as.factor(c(data.all[,2],data.all[,2])))
  anova(fit)
```

```{r ICV}
  # Difference between groups for ICV
  fit <- lm(formula = c(data.all[,108], data.all[,182]) ~ as.factor(c(data.all[,2],data.all[,2])))
  anova(fit)

  # Average (sd) of ICV for every group
    print(paste("ICV in transgender women is ", mean(c(data.all[data.all[,2]==1,108], data.all[data.all[,2]==1,182])), " (",sd(c(data.all[data.all[,2]==1,108], data.all[data.all[,2]==1,182])),"), ",sep = ""))
    print(paste("ICV in transgender women is ", mean(c(data.all[data.all[,2]==2,108], data.all[data.all[,2]==2,182])), " (",sd(c(data.all[data.all[,2]==2,108], data.all[data.all[,2]==2,182])),"), ",sep = ""))
    print(paste("ICV in transgender women is ", mean(c(data.all[data.all[,2]==3,108], data.all[data.all[,2]==3,182])), " (",sd(c(data.all[data.all[,2]==3,108], data.all[data.all[,2]==3,182])),"), ",sep = ""))
    print(paste("ICV in transgender women is ", mean(c(data.all[data.all[,2]==4,108], data.all[data.all[,2]==4,182])), " (",sd(c(data.all[data.all[,2]==4,108], data.all[data.all[,2]==4,182])),"), ",sep = ""))


```

```{r plot age, echo = FALSE, fig.width=6, fig.height=4}
hist(data.all[,4], xlab = "Age", col = "seagreen3", border = "white", 
     main = "Histogram of age distribution", breaks = length(unique(data.all[,4])))
```

## Results
### Repeated measures
Questions: do we want (to account for) a correlation between regions of the same participant?  

First we look at the results for volume. Alle code used to compute this can be found in the markdown version of this document, but is not printed in the pdf.  
We fitted two mixed models. In both models a random intercept for every subject was added, in the second model age and total intracranial volume are added as covariates. An ANOVA is conducted on the results of the mixed models and an FDR-correction is applied over regions. For the regions where a statistically significant difference was found between the groups post-hoc paired comparisons were conducted that are bonferroni-corrected. The results for both models are very similar. A summary of the results is displayed in the tables below. All code used to obtain these results can be found in the RmD file.  
If no statistically significant difference is found for a regions "NA" is printed for that region.  


```{r repmes volume, warning = FALSE, include = FALSE, echo = FALSE, results='hide'}
# Libraries
  library(lme4)
  library(lmerTest)
  library(nlme)

# Construct data
  nt <- 2      # number of time points
  nr <- 22     # number of regions
  nsu <- dim(data.hyp)[1] # number of subjects
  
  suc <- rep(1:nsu, each = nt*nr)
  tc <- rep(1:nt, nsu*nr)
  rc <- rep(rep(1:nr, each = nt), nsu)
  gc <- c(rep(c("CW", "CM"), each = (30*nt*nr)), rep(c("TM", "TW"), each = (40*nt*nr)))
  
  j <- 0       # counter for for-loop
  pasteregions <- array(data=NA, dim = nt*nr*nsu)
  for(i in 1:nsu){
    for(r in 1:nr){
      for(t in 1:nt){
        j <- j + 1
        pasteregions[j] <- data.hyp[i,(2+r+(nr*(t-1)))]
      }
    }
  }
  data.lmer <- data.frame(cbind(suc, tc, rc, gc, pasteregions))
  data.lmer$pasteregions <- as.numeric(as.character(data.lmer$pasteregions))

# Fit model
  p.lme <- array(data = NA, dim = nr)
  for(r in 1:nr){
    mf <- lme(pasteregions ~ 1 + as.factor(gc), random = ~ 1 | suc, data = data.lmer[data.lmer$rc == r,])
    p.lme[r] <- anova(mf)[2,4]
  }
  
# FDR correction
  p.corr <- p.adjust(p.lme, method = "fdr")
  
# Data with age and total volume
  ac <- rep(data.all[,4], each = nt*nr)
  vc.1 <- rep(data.all[,108], each = nt*nr)
  vc.2 <- rep(data.all[,182], each = nt*nr)
  vc <- vc.1
  for(i in 1:length(vc.1)){
    vc[i] <- ifelse(i%%2==1, vc.1[i], vc.2[i])
  }
  
  data.lmer.cov <- data.frame(cbind(suc, tc, rc, gc, ac, vc, pasteregions))
  data.lmer.cov$pasteregions <- as.numeric(as.character(data.lmer.cov$pasteregions))
  data.lmer.cov$ac <- as.numeric(as.character(data.lmer.cov$ac))
  data.lmer.cov$vc <- as.numeric(as.character(data.lmer.cov$vc))
  
# Fit model
  p.lme.cov <- array(data = NA, dim = nr)
  for(r in 1:nr){
    mf <- lme(pasteregions ~ 1 + as.factor(gc) + ac + vc, random = ~ 1 | suc, data = data.lmer.cov[data.lmer.cov$rc == r,])
    p.lme.cov[r] <- anova(mf)[2,4]
  }
  
# FDR correction
  p.cov.corr <- p.adjust(p.lme.cov, method = "fdr")
```

```{r posthoc, warning = FALSE, include = FALSE, echo = FALSE, results='hide'}
# create an object with all possible combinations
  allcomb <- combn(c(1:4),2)
  allcomb.txt <- array(data=NA, dim = dim(allcomb)[2])
  labels.g <- c("CW", "CM", "TM", "TW")
  for(i in 1:dim(allcomb)[2]){
    allcomb.txt[i] <- paste(labels.g[allcomb[1,i]], " vs ", labels.g[allcomb[2,i]], sep="")
  }


# create object to save results
  pt <- array(data = NA, dim = c(length(p.corr),dim(allcomb)[2]))
  pt.cov <- array(data = NA, dim = c(length(p.cov.corr),dim(allcomb)[2]))
  pt.corr <- array(data=NA, dim=dim(pt))
  pt.cov.corr <- array(data=NA, dim=dim(pt.cov))
  
# Pair-wise comparisons and correct results
  for(r in 1:nr){
    # general model
    if (p.corr[r] > 0.05) {
        pt[r,] <- rep(NA,dim(allcomb)[2])
    }else{
      for(i in 1:dim(allcomb)[2]){
        pt[r,i] <- unlist(summary(lme(pasteregions ~ 1 + as.factor(gc), random = ~ 1 | suc, data = data.lmer[data.lmer$rc == r & (data.lmer$gc == labels.g[allcomb[1,i]] | data.lmer$gc == labels.g[allcomb[2,i]]),]))[20])[10]
      }
      pt.corr[r,] <- p.adjust(pt[r,], method = "bonferroni")
    }
    # model with covariates
    if (p.cov.corr[r] > 0.05) {
        pt.cov[r,] <- rep(NA,dim(allcomb)[2])
    }else{
      for(i in 1:dim(allcomb)[2]){
        pt.cov[r,i] <- summary(lme(pasteregions ~ 1 + as.factor(gc) + ac + vc, random = ~ 1 | suc, data = data.lmer.cov[data.lmer.cov$rc == r & (data.lmer.cov$gc == labels.g[allcomb[1,i]] | data.lmer.cov$gc == labels.g[allcomb[2,i]]),]))$tTable[2,5]
      }
      pt.cov.corr[r,] <- p.adjust(pt.cov[r,], method = "bonferroni")
    }
  }
```

```{r tables, echo = FALSE}
# Model without covariates
  kable(rbind(c("region", "No cov", "With cov"), cbind(substring(names(data.hyp[,3:24]), 4), round(p.corr, 3), round(p.cov.corr, 3))), caption = "ANOVA for volume with FDR correction for model with and without covariates")
  kable(cbind(c(" ",substring(names(data.hyp[,3:24]), 4)),rbind(c("CW vs CM","CW vs TM","CW vs TW","CM vs TM","CM vs TW","TM vs TW"),round(pt.corr,3))), caption = "Group-wise comparison for volume in model with no covariates")
  kable(cbind(c(" ",substring(names(data.hyp[,3:24]), 4)),rbind(c("CW vs CM","CW vs TM","CW vs TW","CM vs TM","CM vs TW","TM vs TW"),round(pt.cov.corr,3))), caption = "Group-wise comparison for volume in model with covariates")
```

Then we do the same computations for thickness. However, for thickness we leave out intracranial volume as a covariate (cf. e-mail Meredith Braskie).   

```{r repmes thickn, warning = FALSE, include = FALSE, echo = FALSE, results='hide'}
# Construct data
  nt <- 2      # number of time points
  nr <- 10     # number of regions
  nsu <- dim(data.hyp)[1]
  
  suc <- rep(1:nsu, each = nt*nr)
  tc <- rep(1:nt, nsu*nr)
  rc <- rep(rep(1:nr, each = nt), nsu)
  gc <- c(rep(c("CW", "CM"), each = (30*nt*nr)), rep(c("TM", "TW"), each = (40*nt*nr)))
  
  j <- 0       # counter for for-loop
  pasteregions <- array(data=NA, dim = nt*nr*nsu)
  for(i in 1:nsu){
    for(r in 1:nr){
      for(t in 1:nt){
        j <- j + 1
        pasteregions[j] <- data.hyp[i,(46+r+(nr*(t-1)))]
      }
    }
  }
  data.lmer.th <- data.frame(cbind(suc, tc, rc, gc, pasteregions))
  data.lmer.th$pasteregions <- as.numeric(as.character(data.lmer.th$pasteregions))

# Fit model
  p.lme.th <- array(data = NA, dim = nr)
  for(r in 1:nr){
    mf <- lme(pasteregions ~ 1 + as.factor(gc), random = ~ 1 | suc, data = data.lmer.th[data.lmer.th$rc == r,])
    p.lme.th[r] <- anova(mf)[2,4]
  }
  
# FDR correction
  p.corr.th <- p.adjust(p.lme.th, method = "fdr")
  
# Data with age and total volume
  ac <- rep(data.all[,4], each = nt*nr)

  data.lmer.th.cov <- data.frame(cbind(suc, tc, rc, gc, ac, pasteregions))
  data.lmer.th.cov$pasteregions <- as.numeric(as.character(data.lmer.th.cov$pasteregions))
  data.lmer.th.cov$ac <- as.numeric(as.character(data.lmer.th.cov$ac))
  
# Fit model
  p.lme.cov.th <- array(data = NA, dim = nr)
  for(r in 1:nr){
    mf <- lme(pasteregions ~ 1 + as.factor(gc) + ac, random = ~ 1 | suc, data = data.lmer.th.cov[data.lmer.th.cov$rc == r,])
    p.lme.cov.th[r] <- anova(mf)[2,4]
  }
  
# FDR correction
  p.cov.corr.th <- p.adjust(p.lme.cov.th, method = "fdr")
```

```{r posthoc thick, warning = FALSE, include = FALSE, echo = FALSE, results='hide'}
# create object to save results
  pt.th <- array(data = NA, dim = c(length(p.corr.th),dim(allcomb)[2]))
  pt.cov.th <- array(data = NA, dim = c(length(p.cov.corr.th),dim(allcomb)[2]))
  pt.corr.th <- array(data=NA, dim=dim(pt.th))
  pt.cov.corr.th <- array(data=NA, dim=dim(pt.cov.th))
  
# Pair-wise comparisons and correct results
  for(r in 1:nr){
    # general model
    if (p.corr.th[r] > 0.05) {
        pt.th[r,] <- rep(NA,dim(allcomb)[2])
    }else{
      for(i in 1:dim(allcomb)[2]){
        pt.th[r,i] <- unlist(summary(lme(pasteregions ~ 1 + as.factor(gc), random = ~ 1 | suc, data = data.lmer.th[data.lmer.th$rc == r & (data.lmer.th$gc == labels.g[allcomb[1,i]] | data.lmer.th$gc == labels.g[allcomb[2,i]]),]))[20])[10]
      }
      pt.corr.th[r,] <- p.adjust(pt.th[r,], method = "bonferroni")
    }
    # model with covariates
    if (p.cov.corr.th[r] > 0.05) {
        pt.cov.th[r,] <- rep(NA,dim(allcomb)[2])
    }else{
      for(i in 1:dim(allcomb)[2]){
        pt.cov.th[r,i] <- summary(lme(pasteregions ~ 1 + as.factor(gc) + ac, random = ~ 1 | suc, data = data.lmer.th.cov[data.lmer.th.cov$rc == r & (data.lmer.th.cov$gc == labels.g[allcomb[1,i]] | data.lmer.th.cov$gc == labels.g[allcomb[2,i]]),]))$tTable[2,5]
      }
      pt.cov.corr.th[r,] <- p.adjust(pt.cov.th[r,], method = "bonferroni")
    }
  }
```


```{r tables thick, echo = FALSE}
# Model without covariates
  kable(rbind(c("region", "No cov", "With cov"), cbind(substring(names(data.hyp[,47:56]), 4), round(p.corr.th, 3), round(p.cov.corr.th, 3))), caption = "ANOVA for thickness with FDR correction for model with and without covariates")
  kable(cbind(c(" ",substring(names(data.hyp[,47:56]), 4)),rbind(c("CW vs CM","CW vs TM","CW vs TW","CM vs TM","CM vs TW","TM vs TW"),round(pt.corr.th,3))), caption = "Group-wise comparison for thickness in model with no covariates")
    kable(cbind(c(" ",substring(names(data.hyp[,47:56]), 4)),rbind(c("CW vs CM","CW vs TM","CW vs TW","CM vs TM","CM vs TW","TM vs TW"),round(pt.cov.corr.th,3))), caption = "Group-wise comparison for thickness in model with covariates")
```

And surface area:  

```{r repmes surfavg, warning = FALSE, include = FALSE, echo = FALSE, results='hide'}
# Construct data
  nt <- 2      # number of time points
  nr <- 10     # number of regions
  nsu <- dim(data.hyp)[1]
  
  suc <- rep(1:nsu, each = nt*nr)
  tc <- rep(1:nt, nsu*nr)
  rc <- rep(rep(1:nr, each = nt), nsu)
  gc <- c(rep(c("CW", "CM"), each = (30*nt*nr)), rep(c("TM", "TW"), each = (40*nt*nr)))
  
  j <- 0       # counter for for-loop
  pasteregions <- array(data=NA, dim = nt*nr*nsu)
  for(i in 1:nsu){
    for(r in 1:nr){
      for(t in 1:nt){
        j <- j + 1
        pasteregions[j] <- data.hyp[i,(66+r+(nr*(t-1)))]
      }
    }
  }
  data.lmer.sa <- data.frame(cbind(suc, tc, rc, gc, pasteregions))
  data.lmer.sa$pasteregions <- as.numeric(as.character(data.lmer.sa$pasteregions))

# Fit model
  p.lme.sa <- array(data = NA, dim = nr)
  for(r in 1:nr){
    mf <- lme(pasteregions ~ 1 + as.factor(gc), random = ~ 1 | suc, data = data.lmer.sa[data.lmer.sa$rc == r,])
    p.lme.sa[r] <- anova(mf)[2,4]
  }
  
# FDR correction
  p.corr.sa <- p.adjust(p.lme.sa, method = "fdr")
  
# Data with age and total volume
  ac <- rep(data.all[,4], each = nt*nr)
  vc.1 <- rep(data.all[,108], each = nt*nr)
  vc.2 <- rep(data.all[,182], each = nt*nr)
  vc <- vc.1
  for(i in 1:length(vc.1)){
    vc[i] <- ifelse(i%%2==1, vc.1[i], vc.2[i])
  }
  
  data.lmer.sa.cov <- data.frame(cbind(suc, tc, rc, gc, ac, vc, pasteregions))
  data.lmer.sa.cov$pasteregions <- as.numeric(as.character(data.lmer.sa.cov$pasteregions))
  data.lmer.sa.cov$ac <- as.numeric(as.character(data.lmer.sa.cov$ac))
  data.lmer.sa.cov$vc <- as.numeric(as.character(data.lmer.sa.cov$vc))
  
# Fit model
  p.lme.cov.sa <- array(data = NA, dim = nr)
  for(r in 1:nr){
    mf <- lme(pasteregions ~ 1 + as.factor(gc) + ac + vc, random = ~ 1 | suc, data = data.lmer.sa.cov[data.lmer.sa.cov$rc == r,])
    p.lme.cov.sa[r] <- anova(mf)[2,4]
  }
  
# FDR correction
  p.cov.corr.sa <- p.adjust(p.lme.cov.sa, method = "fdr")
```

```{r posthoc surfavg, warning = FALSE, include = FALSE, echo = FALSE, results='hide'}
# create object to save results
  pt.sa <- array(data = NA, dim = c(length(p.corr.sa),dim(allcomb)[2]))
  pt.cov.sa <- array(data = NA, dim = c(length(p.cov.corr.sa),dim(allcomb)[2]))
  pt.corr.sa <- array(data=NA, dim=dim(pt.sa))
  pt.cov.corr.sa <- array(data=NA, dim=dim(pt.cov.sa))
  
# Pair-wise comparisons and correct results
  for(r in 1:nr){
    # general model
    if (p.corr.sa[r] > 0.05) {
        pt.sa[r,] <- rep(NA,dim(allcomb)[2])
    }else{
      for(i in 1:dim(allcomb)[2]){
        pt.sa[r,i] <- unlist(summary(lme(pasteregions ~ 1 + as.factor(gc), random = ~ 1 | suc, data = data.lmer.sa[data.lmer.sa$rc == r & (data.lmer.sa$gc == labels.g[allcomb[1,i]] | data.lmer.sa$gc == labels.g[allcomb[2,i]]),]))[20])[10]
      }
      pt.corr.sa[r,] <- p.adjust(pt.sa[r,], method = "bonferroni")
    }
    # model with covariates
    if (p.cov.corr.sa[r] > 0.05) {
        pt.cov.sa[r,] <- rep(NA,dim(allcomb)[2])
    }else{
      for(i in 1:dim(allcomb)[2]){
        pt.cov.sa[r,i] <- summary(lme(pasteregions ~ 1 + as.factor(gc) + ac + vc, random = ~ 1 | suc, data = data.lmer.sa.cov[data.lmer.sa.cov$rc == r & (data.lmer.sa.cov$gc == labels.g[allcomb[1,i]] | data.lmer.sa.cov$gc == labels.g[allcomb[2,i]]),]))$tTable[2,5]
      }
      pt.cov.corr.sa[r,] <- p.adjust(pt.cov.sa[r,], method = "bonferroni")
    }
  }
```


```{r tables surfavg, echo = FALSE, warning = FALSE}
# Model without covariates
  kable(rbind(c("region", "No cov", "With cov"), cbind(substring(names(data.hyp[,67:(67+nr-1)]), 4), round(p.corr.sa, 3), round(p.cov.corr.sa, 3))), caption = "ANOVA for surface area with FDR correction for model with and without covariates")
  kable(cbind(c(" ",substring(names(data.hyp[,67:86]), 4)),rbind(c("CW vs CM","CW vs TM","CW vs TW","CM vs TM","CM vs TW","TM vs TW"),round(pt.corr.sa,3))), caption = "Group-wise comparison for surface area in model with no covariates")
    kable(cbind(c(" ",substring(names(data.hyp[,67:86]), 4)),rbind(c("CW vs CM","CW vs TM","CW vs TW","CM vs TM","CM vs TW","TM vs TW"),round(pt.cov.corr.sa,3))), caption = "Group-wise comparison for surface area in model with covariates")
```

\newpage

### Compare to results from 1 timepoint  
In this section we explore what the results would have been if we would have been limited to one scan for the analysis. Tables 11, 14 and 17 show the p-values resulting from the anova's conducted on T1, T2 and both respectively. You can observe the (slight) increase in power here.  
```{r T1, echo = FALSE, warning = FALSE, include = FALSE}
# Volume
  # Dataset
    T1.vol <- data.hyp[,3:24]

  # FDR
    nr <- 22
    p.T1.vol <- array(data=NA, dim = nr)
    for(r in 1:nr){
      p.T1.vol[r] <- unlist(anova(lm(T1.vol[,r] ~ as.factor(data.hyp[,2]))))[9]
    }
    p.T1.vol.corr <- p.adjust(p.T1.vol, method = "fdr")
    
  # Bonferroni
    pt.T1.vol <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    pt.T1.vol.corr <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    for(r in 1:nr){
      # general model
      if (p.T1.vol.corr[r] > 0.05) {
          pt.T1.vol[r,] <- rep(NA,dim(allcomb)[2])
      }else{
        for(i in 1:dim(allcomb)[2]){
          pt.T1.vol[r,i] <- unlist(t.test(T1.vol[data.hyp[,2]==allcomb[1,i],r], T1.vol[data.hyp[,2]==allcomb[2,i],r])[3])
        }
        pt.T1.vol.corr[r,] <- p.adjust(pt.T1.vol[r,], method = "bonferroni")
      }
    }

    
# Thickness
  # Dataset
    T1.th <- data.hyp[,47:56]

  # FDR
    nr <- 10
    p.T1.th <- array(data=NA, dim = nr)
    for(r in 1:nr){
      p.T1.th[r] <- unlist(anova(lm(T1.th[,r] ~ as.factor(data.hyp[,2]))))[9]
    }
    p.T1.th.corr <- p.adjust(p.T1.th, method = "fdr")
    
  # Bonferroni
    pt.T1.th <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    pt.T1.th.corr <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    for(r in 1:nr){
      # general model
      if (p.T1.th.corr[r] > 0.05) {
          pt.T1.th[r,] <- rep(NA,dim(allcomb)[2])
      }else{
        for(i in 1:dim(allcomb)[2]){
          pt.T1.th[r,i] <- unlist(t.test(T1.th[data.hyp[,2]==allcomb[1,i],r], T1.th[data.hyp[,2]==allcomb[2,i],r])[3])
        }
        pt.T1.th.corr[r,] <- p.adjust(pt.T1.th[r,], method = "bonferroni")
      }
    }

    
# Surface Area
  # Dataset
    T1.sa <- data.hyp[,67:76]

  # FDR
    nr <- 10
    p.T1.sa <- array(data=NA, dim = nr)
    for(r in 1:nr){
      p.T1.sa[r] <- unlist(anova(lm(T1.sa[,r] ~ as.factor(data.hyp[,2]))))[9]
    }
    p.T1.sa.corr <- p.adjust(p.T1.sa, method = "fdr")
    
  # Bonferroni
    pt.T1.sa <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    pt.T1.sa.corr <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    for(r in 1:nr){
      # general model
      if (p.T1.sa.corr[r] > 0.05) {
          pt.T1.sa[r,] <- rep(NA,dim(allcomb)[2])
      }else{
        for(i in 1:dim(allcomb)[2]){
          pt.T1.sa[r,i] <- unlist(t.test(T1.sa[data.hyp[,2]==allcomb[1,i],r], T1.sa[data.hyp[,2]==allcomb[2,i],r])[3])
        }
        pt.T1.sa.corr[r,] <- p.adjust(pt.T1.sa[r,], method = "bonferroni")
      }
    }

```

```{r T2, echo = FALSE, warning = FALSE, include = FALSE}
# Volume
  # Dataset
    T2.vol <- data.hyp[,25:46]

  # FDR
    nr <- 22
    p.T2.vol <- array(data=NA, dim = nr)
    for(r in 1:nr){
      p.T2.vol[r] <- unlist(anova(lm(T2.vol[,r] ~ as.factor(data.hyp[,2]))))[9]
    }
    p.T2.vol.corr <- p.adjust(p.T2.vol, method = "fdr")
    
  # Bonferroni
    pt.T2.vol <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    pt.T2.vol.corr <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    for(r in 1:nr){
      # general model
      if (p.T2.vol.corr[r] > 0.05) {
          pt.T2.vol[r,] <- rep(NA,dim(allcomb)[2])
      }else{
        for(i in 1:dim(allcomb)[2]){
          pt.T2.vol[r,i] <- unlist(t.test(T2.vol[data.hyp[,2]==allcomb[1,i],r], T2.vol[data.hyp[,2]==allcomb[2,i],r])[3])
        }
        pt.T2.vol.corr[r,] <- p.adjust(pt.T2.vol[r,], method = "bonferroni")
      }
    }

    
# Thickness
  # Dataset
    T2.th <- data.hyp[,57:66]

  # FDR
    nr <- 10
    p.T2.th <- array(data=NA, dim = nr)
    for(r in 1:nr){
      p.T2.th[r] <- unlist(anova(lm(T2.th[,r] ~ as.factor(data.hyp[,2]))))[9]
    }
    p.T2.th.corr <- p.adjust(p.T2.th, method = "fdr")
    
  # Bonferroni
    pt.T2.th <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    pt.T2.th.corr <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    for(r in 1:nr){
      # general model
      if (p.T2.th.corr[r] > 0.05) {
          pt.T2.th[r,] <- rep(NA,dim(allcomb)[2])
      }else{
        for(i in 1:dim(allcomb)[2]){
          pt.T2.th[r,i] <- unlist(t.test(T2.th[data.hyp[,2]==allcomb[1,i],r], T2.th[data.hyp[,2]==allcomb[2,i],r])[3])
        }
        pt.T2.th.corr[r,] <- p.adjust(pt.T2.th[r,], method = "bonferroni")
      }
    }

    
# Surface Area
  # Dataset
    T2.sa <- data.hyp[,77:86]

  # FDR
    nr <- 10
    p.T2.sa <- array(data=NA, dim = nr)
    for(r in 1:nr){
      p.T2.sa[r] <- unlist(anova(lm(T2.sa[,r] ~ as.factor(data.hyp[,2]))))[9]
    }
    p.T2.sa.corr <- p.adjust(p.T2.sa, method = "fdr")
    
  # Bonferroni
    pt.T2.sa <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    pt.T2.sa.corr <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    for(r in 1:nr){
      # general model
      if (p.T2.sa.corr[r] > 0.05) {
          pt.T2.sa[r,] <- rep(NA,dim(allcomb)[2])
      }else{
        for(i in 1:dim(allcomb)[2]){
          pt.T2.sa[r,i] <- unlist(t.test(T2.sa[data.hyp[,2]==allcomb[1,i],r], T2.sa[data.hyp[,2]==allcomb[2,i],r])[3])
        }
        pt.T2.sa.corr[r,] <- p.adjust(pt.T2.sa[r,], method = "bonferroni")
      }
    }
```

### Volume  
```{r compare V, echo = FALSE, warning = FALSE}
# Volume
  kable(caption = "Comparison volume", rbind(c("Region", "T1", "T2", "Rep Meas"), cbind(substring(names(data.hyp[,3:24]), 4), round(p.T1.vol.corr, 3), round(p.T2.vol.corr, 3), round(p.corr, 3))))
  kable(caption = "Volume: T1", cbind(c(" ",substring(names(data.hyp[,3:24]), 4)),rbind(c("CW vs CM","CW vs TM","CW vs TW","CM vs TM","CM vs TW","TM vs TW"),round(pt.T1.vol.corr,3))))
  kable(caption = "Volume: T2", cbind(c(" ",substring(names(data.hyp[,25:46]), 4)),rbind(c("CW vs CM","CW vs TM","CW vs TW","CM vs TM","CM vs TW","TM vs TW"),round(pt.T2.vol.corr,3))))
```

### Thickness  
```{r compare T, echo = FALSE, warning = FALSE}
# Thickness
  kable(caption = "Comparison thickness", rbind(c("Region", "T1", "T2", "Rep Meas"), cbind(substring(names(data.hyp[,47:56]), 4), round(p.T1.th.corr, 3), round(p.T2.th.corr, 3), round(p.corr.th, 3))))
  kable(caption = "Thickness: T1", cbind(c(" ",substring(names(data.hyp[,47:56]), 4)),rbind(c("CW vs CM","CW vs TM","CW vs TW","CM vs TM","CM vs TW","TM vs TW"),round(pt.T1.th.corr,3))))
  kable(caption = "Thickness: T2", cbind(c(" ",substring(names(data.hyp[,57:66]), 4)),rbind(c("CW vs CM","CW vs TM","CW vs TW","CM vs TM","CM vs TW","TM vs TW"),round(pt.T1.th.corr,3))))
```

### Surface Area  
```{r compare SA, echo = FALSE, warning = FALSE}
# Surface area
  kable(caption = "Comparison surface area", rbind(c("Region", "T1", "T2", "Rep Meas"), cbind(substring(names(data.hyp[,67:76]), 4), round(p.T1.sa.corr, 3), round(p.T2.sa.corr, 3), round(p.corr.sa, 3))))
  kable(caption = "Surface area: T1", cbind(c(" ",substring(names(data.hyp[,67:76]), 4)),rbind(c("CW vs CM","CW vs TM","CW vs TW","CM vs TM","CM vs TW","TM vs TW"),round(pt.T1.sa.corr,3))))
  kable(caption = "Surface area: T2", cbind(c(" ",substring(names(data.hyp[,77:86]), 4)),rbind(c("CW vs CM","CW vs TM","CW vs TW","CM vs TM","CM vs TW","TM vs TW"),round(pt.T2.sa.corr,3))))
```



## Simulations
First we need to define the parameters of our simulations. 
```{r Simulations parameters}
# variance/sd epsilon
seps <- 1  

# Number of simulations
asim <- 5000

# Effect size
delta <- 0.8

# Number of participants
n <- 30
n.1 <- n/2      # in the first group
n.2 <- n/2      # in the second group

# Level of statistical significance 
alpha <- 0.05

# Correlation between first and second measurement
rho <- seq(0.01,0.99,0.01)
```

Then we prepare objects to store our results
```{r Create objects}
# Number of simulations
pow.mean1<-vector("numeric",length(rho))
pow.mean2<-vector("numeric",length(rho))
pow.mean3<-vector("numeric",length(rho))
```

```{r Load libraries, warning=FALSE, echo=FALSE}
# Libraries
library(lme4)
library(lmerTest)
```

```{r Simulations code}
# Loop over preset correlations between measure 1 and measure 2
for(i in 1:length(rho)){
  # Create objects to store power in for every simulations
  pow.1<-vector("numeric",asim)
  pow.2<-vector("numeric",asim)
  pow.3<-vector("numeric",asim)
  
  for(k in 1:asim){
    # Scenario 1: lower bound of power curve
    # two groups with equal amount of subjects, groups differ with an effect size delta
    # Construct a vector that determines in which group each subject falls
    x<-c(rep(1,n.1),rep(0,n.2))
    
    # Vector with observations in the set of participants
    y<-rnorm(n,0,seps)       
    # Add an effect size to the first group
    y[1:n.1]<-y[1:n.1]+delta    
    
    # Boolean of whether an effect is detected, this is later used to compute the power
    pow.1[k]<-summary(lm(y~x))$coef[2,4]<alpha
    
    # Scenario 2: upper bound of power curve
    # two groups with equal amount of subjects, twice as many as scenario 1, groups differ with an effect size delta
    # Construct a vector that determines in which group each subject falls
    x2<-c(rep(1,(n.1*2)),rep(0,(n.2*2)))
    
    # Vector with observations in the set of participants
    y2<-rnorm(n*2,0,seps)
    # Add an effect size to the first group
    y2[1:(n.1*2)]<-y2[1:(n.1*2)]+delta
    
    # Boolean of whether an effect is detected, this is later used to compute the power
    pow.2[k]<-summary(lm(y2~x2))$coef[2,4]<alpha
    
    # Scenario 3: two measurements for every subject, same amount of subjects as in scenario 1
    # two groups with equal amount of subjects, correlation between measurements, groups differ with an effect size delta
    # Construct a vector that determines in which group each subject falls
    x3<-c(rep(1,n.1),rep(0,n.2),rep(1,n.1),rep(0,n.2))
    
    # Vector with first observation of every participant
    y3<-rnorm(n,0,seps)
    # Factor to multiply second set of observations with to obtain results in line with predefined correlation
    alpac<-sqrt(rho[i]^2/(1-rho[i]^2)*seps)
    # Construct second set of observations that are correlated with first set (y3)
    y3.2u<-alpac*y3+rnorm(n)
    y3.2<-y3.2u/sqrt(var(y3.2u))
    # Add effect size to the first group of participants
    y3[1:n.1]<-y3[1:n.1]+delta
    y3.2[1:n.1]<-y3.2[1:n.1]+delta
    # Combine both observations in 1 vector
    y3o<-c(y3,y3.2)
    
    # Define subject numbers
    subject<-rep(1:n,2)
    # Construct mixed model
    mm<-lmer(y3o ~ x3 + (1 | subject))   
    # Boolean of whether an effect is detected, this is later used to compute the power
    pow.3[k]<-summary(mm)$coef[2,5]<alpha
  }
  pow.mean1[i]<-mean(pow.1)
  pow.mean2[i]<-mean(pow.2)
  pow.mean3[i]<-mean(pow.3)
}
```


### Results
```{r Simulation results}
# power of taking both measures into account
plot(rho,pow.mean1,ylim=c(0,1), type="l")
par(new = TRUE)
plot(rho,pow.mean2,ylim=c(0,1), type="l")
par(new = TRUE)
plot(rho,pow.mean3,ylim=c(0,1), type="p", col = "goldenrod3", pch = 2)
```

```{r simulation and correlation}
# Correlation of hypothesized regions
  bg.one <- 3
  nd.one <- 24
  ln <- nd.one-bg.one+1
  
  # Object to save correlations
  corrhyp <- array(data=NA, dim = ln)
  
  # Compute ANOVA for every predictor and save p-value
  for(i in bg.one:nd.one){
    corrhyp[i-bg.one+1] <- cor(x = data.hyp[,i], y = data.hyp[,i + ln])
  }
  
  summary(corrhyp)
  
  kable(cbind(substring(names(data.hyp[,3:24]), 4), round(corrhyp,3)))
  
  
# Increase in power
```



# Added later
Here we create table 11 and table 17 with covariates

## Table 11 - Volume
```{r T1T2cov vol}
# Volume
  # Dataset
    T1.vol <- data.hyp[,3:24]

  # FDR
    nr <- 22
    p.T1.vol.cov <- array(data=NA, dim = nr)
    for(r in 1:nr){
      p.T1.vol.cov[r] <- unlist(summary(aov(T1.vol[,r] ~ data.all[,4] + data.all[,108] + as.factor(data.hyp[,2]))))[19]
    }
    p.T1.vol.cov.corr <- p.adjust(p.T1.vol.cov, method = "fdr")
    
  # Bonferroni
    # pt.T1.vol.cov <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    # pt.T1.vol.cov.corr <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    # for(r in 1:nr){
    #   # general model
    #   if (p.T1.vol.cov.corr[r] > 0.05) {
    #       pt.T1.vol.cov[r,] <- rep(NA,dim(allcomb)[2])
    #   }else{
    #     for(i in 1:dim(allcomb)[2]){
    #       pt.T1.vol.cov[r,i] <- unlist(t.test(T1.vol[data.hyp[,2]==allcomb[1,i],r], T1.vol[data.hyp[,2]==allcomb[2,i],r])[3])
    #     }
    #     pt.T1.vol.cov.corr[r,] <- p.adjust(pt.T1.vol.cov[r,], method = "bonferroni")
    #   }
    # }

    

# Surface Area
  # Dataset
    T1.sa <- data.hyp[,67:76]

  # FDR
    nr <- 10
    p.T1.sa.cov <- array(data=NA, dim = nr)
    for(r in 1:nr){
      p.T1.sa.cov[r] <- unlist(summary(aov(T1.sa[,r] ~ data.all[,4] + data.all[,108] + as.factor(data.hyp[,2]))))[19]
    }
    p.T1.sa.cov.corr <- p.adjust(p.T1.sa.cov, method = "fdr")
    
  # Bonferroni
    # pt.T1.sa.cov <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    # pt.T1.sa.cov.corr <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    # for(r in 1:nr){
    #   # general model
    #   if (p.T1.sa.cov.corr[r] > 0.05) {
    #       pt.T1.sa.cov[r,] <- rep(NA,dim(allcomb)[2])
    #   }else{
    #     for(i in 1:dim(allcomb)[2]){
    #       pt.T1.sa.cov[r,i] <- unlist(t.test(T1.sa[data.hyp[,2]==allcomb[1,i],r], T1.sa[data.hyp[,2]==allcomb[2,i],r])[3])
    #     }
    #     pt.T1.sa.cov.corr[r,] <- p.adjust(pt.T1.sa[r,], method = "bonferroni")
    #   }
    # }

# Volume
  # Dataset
    T2.vol <- data.hyp[,25:46]

  # FDR
    nr <- 22
    p.T2.vol.cov <- array(data=NA, dim = nr)
    for(r in 1:nr){
      p.T2.vol.cov[r] <- unlist(summary(aov(T2.vol[,r] ~ data.all[,4] + data.all[,108] + as.factor(data.hyp[,2]))))[19]
    }
    p.T2.vol.cov.corr <- p.adjust(p.T2.vol.cov, method = "fdr")
    
  # Bonferroni
    # pt.T2.vol.cov <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    # pt.T2.vol.cov.corr <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    # for(r in 1:nr){
    #   # general model
    #   if (p.T2.vol.cov.corr[r] > 0.05) {
    #       pt.T2.vol.cov[r,] <- rep(NA,dim(allcomb)[2])
    #   }else{
    #     for(i in 1:dim(allcomb)[2]){
    #       pt.T2.vol.cov[r,i] <- unlist(t.test(T2.vol[data.hyp[,2]==allcomb[1,i],r], T2.vol[data.hyp[,2]==allcomb[2,i],r])[3])
    #     }
    #     pt.T2.vol.cov.corr[r,] <- p.adjust(pt.T2.vol.cov[r,], method = "bonferroni")
    #   }
    # }

    
# Surface Area
  # Dataset
    T2.sa <- data.hyp[,77:86]

  # FDR
    nr <- 10
    p.T2.sa.cov <- array(data=NA, dim = nr)
    for(r in 1:nr){
      p.T2.sa.cov[r] <- unlist(summary(aov(T2.sa[,r] ~ data.all[,4] + data.all[,108] + as.factor(data.hyp[,2]))))[19]
    }
    p.T2.sa.cov.corr <- p.adjust(p.T2.sa.cov, method = "fdr")
    
  # Bonferroni
    # pt.T2.sa.cov <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    # pt.T2.sa.cov.corr <- array(data = NA, dim = c(nr,dim(allcomb)[2]))
    # for(r in 1:nr){
    #   # general model
    #   if (p.T2.sa.cov.corr[r] > 0.05) {
    #       pt.T2.sa.cov[r,] <- rep(NA,dim(allcomb)[2])
    #   }else{
    #     for(i in 1:dim(allcomb)[2]){
    #       pt.T2.sa.cov[r,i] <- unlist(t.test(T2.sa[data.hyp[,2]==allcomb[1,i],r], T2.sa[data.hyp[,2]==allcomb[2,i],r])[3])
    #     }
    #     pt.T2.sa.cov.corr[r,] <- p.adjust(pt.T2.sa.cov[r,], method = "bonferroni")
    #   }
    # }
```
```{r table cov}
# Volume
  kable(caption = "Comparison volume", rbind(c("Region", "T1", "T2", "Rep Meas"), cbind(substring(names(data.hyp[,3:24]), 4), round(p.T1.vol.cov.corr, 3), round(p.T2.vol.cov.corr, 3), round(p.cov.corr, 3))))
# Surface area
  kable(caption = "Comparison surface area", rbind(c("Region", "T1", "T2", "Rep Meas"), cbind(substring(names(data.hyp[,67:76]), 4), round(p.T1.sa.cov.corr, 3), round(p.T2.sa.cov.corr, 3), round(p.cov.corr.sa, 3))))
```