BrainAGEsimulation.Rmd

---
title: "BrainAGESim"
output:
  pdf_document: default
  html_document: default
---

```{r}
library(kernlab)
library(ggplot2)
library(e1071)

remove(list = ls())
set.seed(12345)


today <- Sys.Date()
plotABFs <- F
saveplot = FALSE; ntraining = 500; ntesting = 500; nvars = 100; errorEps = 10; plot_prefix = ''; errorEta = 0.5;
beta0 = 100; beta1 = -1; 
errorEpsPct = 2; #2
errorEtaPct = 0.4;  #1.2
# nonDecaysd = 5;
decayPct = 0.75;
scaleFac = 1/80; #1/120; 1/2500
# scaleFac = 1/4500;

nruns <- 100

w1Sd <- 1
w2Sd <- 1
w3Sd <- 1
w1MeanMean <- 20 #2
w2MeanMean <- 20 #2
w3MeanMean <- 20 #2

cbPaletteGrey <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
cbPaletteBlack <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

ntotal <- ntraining + ntesting
#generate random subject ages
age <- runif(ntotal, 20, 56)
ageTrain <- age[1:ntraining]
ageTest <- age[(ntraining+1):ntotal]
ABF1 <- beta0 + age * beta1 # age basis function 1
ABF2 <- 1/424*sqrt(90^5-age^5) - 100 # age basis function 2
ABF3 <- 2.6*100^2/age^2 + 16 # age basis function 3
```



```{r}
# set.seed(1234)
alphas <- c(0, runif(1, -2, -1))
betas <- c(0, runif(1, -2, -1))
gammas <- c(0, runif(1, -2, -1))
albega1 <- expand.grid(alphas, betas, gammas)
albega1[4:5,] <- albega1[c(5,4),]
albega1[c(5,7),] <- albega1[c(7,5),]

alphas <- c(0, runif(1, -2, -1))
betas <- c(0, runif(1, -2, -1))
gammas <- c(0, runif(1, -2, -1))
albega2 <- expand.grid(alphas, betas, gammas)
albega2[4:5,] <- albega2[c(5,4),]
albega2[c(5,7),] <- albega2[c(7,5),]

albega <- rbind(albega1, albega2)
ncombs <- nrow(albega)
albegaSums <- rowSums(abs(albega))
albega <- albega/ifelse(albegaSums != 0, albegaSums, 1)

pCombWiAge.df <- matrix(NA, nrow = ncombs, ncol = nruns)
pCombNoAge.df <- matrix(NA, nrow = ncombs, ncol = nruns)
r2vec <- vector(mode = "numeric", length = nruns)
namesResults <- c("AlBeGa", "CorWithAge", "pNoAge", "pWiAge", "pCombNoAge", "pCombWiAge", "CorWithA")
resultss <- array(NA, c(ncombs, length(namesResults) - 1, nruns))
funcInd <- list() # functional indices
funcInd$ABF1 <- 1:8
funcInd$ABF2 <- 1:8
funcInd$ABF3 <- 1:8
allFuncInd <- unlist(funcInd)
effectSizes <- matrix(NA, ncol = 3, nrow = nruns)

system.time(
for (myrun in 1:nruns){
  w1Mean <- rnorm(nvars, mean = w1MeanMean, sd = w1Sd)
  w2Mean <- rnorm(nvars, mean = w2MeanMean, sd = w2Sd)
  w3Mean <- rnorm(nvars, mean = w3MeanMean, sd = w3Sd)
  
  results <- data.frame(matrix(NA, nrow = ncombs, ncol = length(namesResults)))
  colnames(results) <- namesResults
  covsNoNoise <- t(as.matrix(albega) %*% as.matrix(rbind(ABF1, ABF2, ABF3)))
  errorEtaMat <- matrix(rnorm(ncombs*ntotal, 0, errorEtaPct*abs(covsNoNoise)), 
                        ncol = ncombs, nrow = ntotal, byrow = T)
  rawCovs <- covsNoNoise + errorEtaMat
  covs <- rawCovs
  covsTest <- covs[(ntraining+1):ntotal, ]
  x_total <- matrix(NA, ncol = nvars, nrow = ntotal)
  
  D1 <- rowSums(covs[,funcInd$ABF1, drop = F])*scaleFac
  D2 <- rowSums(covs[,funcInd$ABF2, drop = F])*scaleFac
  D3 <- rowSums(covs[,funcInd$ABF3, drop = F])*scaleFac
  
  backgroundInd <- floor(decayPct*nvars) 
  cohen1 <- vector(mode = "numeric", length = backgroundInd)
  cohen2 <- vector(mode = "numeric", length = backgroundInd)
  cohen3 <- vector(mode = "numeric", length = backgroundInd)
  for (k in 1:nvars){
    if (k <= backgroundInd){
      xNoNoise <- 
        (w1Mean[k] - D1)*ABF1 + 
        (w2Mean[k] - D2)*ABF2 +
        (w3Mean[k] - D3)*ABF3
      x_total[,k] <- xNoNoise + rnorm(ntotal, mean = 0, sd = errorEpsPct*abs(xNoNoise))
    } 
  }
  effectSize <- c(mean(abs((w1Mean-mean(D1))/w1Mean)),
                  mean(abs((w2Mean-mean(D2))/w2Mean)),
                  mean(abs((w3Mean-mean(D3))/w3Mean)))
  effectSizes[myrun, ] <- effectSize
  sdRemain <- sd(x_total[,1:backgroundInd])
  x_total[, (backgroundInd+1):nvars] <- rnorm(ntotal*(nvars - backgroundInd), 
                                              mean = 0, sd = sdRemain)
  names(x_total) <- paste('x', 1:nvars, sep = '')
  x_training <- x_total[1:ntraining, ]
  x_testing <- x_total[(ntraining + 1): ntotal, ]
  
  #fit a model on the trianing data
  train_y <- as.matrix(ageTrain)
  train_x <- as.matrix(x_training)
  this_model <- svm(x = train_x, y = train_y)
  
  names(x_testing) <- paste('x', 1:nvars, sep = '')
  predictions <- predict(this_model, x_testing) # get predictions on the testing set 
  BrainAGE <- predictions - ageTest     #BrainAGE is yhat - y
  
  mae <- mean(abs(BrainAGE))
  ss_total <- sum((ageTest - mean(ageTest))^2)
  ss_resid <- sum((BrainAGE - mean(BrainAGE))^2)
  (r2 <- 1 - ss_resid / ss_total)
  results$CorWithA <- abs(cor(covsTest, ageTest))
  
  for (i in 1:ncombs){
    results[i,"AlBeGa"] <- paste0(albega[i,], collapse = ";")
    results[i, "CorWithAge"] <- cor(covsTest[,i], ageTest)
  }
  
  for (i in 1:ncombs){
    model_i <- lm(BrainAGE ~ covsTest[,i])
    results[i, "pNoAge"] <- summary(model_i)$coefficients["covsTest[, i]", "Pr(>|t|)"]
    modeliAge <- lm(BrainAGE ~ covsTest[,i] + ageTest + ageTest^2 + 1/ageTest)
    results[i, "pWiAge"] <- summary(modeliAge)$coefficients["covsTest[, i]", "Pr(>|t|)"]
  }
  
  myFrame <- data.frame(covsTest, BrainAGE, ageTest)
  myFrameNonlin <- myFrame
  myFrameNonlin$ageSq <- myFrameNonlin$ageTest^2
  myFrameNonlin$ageInv <- 1/myFrameNonlin$ageTest
  model_55 <- lm(BrainAGE ~ . - ageTest, data = myFrame)
  results[, "pCombNoAge"] <- summary(model_55)$coefficients[2:(ncombs + 1), "Pr(>|t|)"]
  modelAge56 <- lm(BrainAGE ~ . , data = myFrame)

  results[, "pCombWiAge"] <- summary(modelAge56)$coefficients[2:(ncombs + 1), "Pr(>|t|)"]
  resultsRound <- results
  resultsRound[,-1] <- round(resultsRound[,-1], 5)  
  predAgedf <- data.frame(ageTest, predictions, BrainAGE)
  plotPred <- F
  if ((myrun == nruns) & (plotPred == T)){
    linreg <- lm(predictions ~ ageTest, data = predAgedf)
    b1 <- linreg$coefficients[2]
    b0 <- linreg$coefficients[1]
    predAgedf$BrainAGER <- linreg$residuals
    a1 <- ggplot(predAgedf, aes(x = ageTest, y = predictions)) + 
      geom_point(alpha = 0.6, shape = 1, size = 0.5) + 
      geom_abline(slope = 1) + labs(x = "Age", y = "Predicted Age") + theme_bw() + 
      geom_smooth(method = "lm", color = cbPaletteBlack[3], se = FALSE)
    a2 <- ggplot(predAgedf, aes(x = ageTest, y = BrainAGE)) + 
      geom_point(shape = 1, size = 0.5) + 
      geom_smooth(method = "lm", color = cbPaletteBlack[3], se = FALSE) + 
      labs(x = "Age", y = "BrainAGE") + theme_bw()
    a3 <- ggplot(predAgedf, aes(x = ageTest, y = BrainAGER)) + 
      geom_point(shape = 1, size = 0.5) + 
      geom_smooth(method = "lm", color = cbPaletteBlack[3], se = FALSE) + 
      labs(x = "Age", y = "BrainAGER") + theme_bw()
    
    # ggsave(file = "predictedVsReal.pdf", plot = a1, height = 3, width = 3)
    # ggsave(file = "resVsReal.pdf", plot = a2, height = 3, width = 3)
    # ggsave(file = "BrainAGERvsAge.pdf", plot = a3, height = 3, width = 3)
  }
  
  pCombWiAge.df[, myrun] <- results$pCombWiAge
  pCombNoAge.df[, myrun] <- results$pCombNoAge
  r2vec[myrun] <- r2
  resultss[ , , myrun] <- as.matrix(results[,-1])
}
)

xCorrs <- rowMeans(resultss[,1,])

pCombWiAge.df <- data.frame(pCombWiAge.df)
pCombNoAge.df <- data.frame(pCombNoAge.df)
pCombWiAge.df$cov <- 1:ncombs
pCombWiAge.df$relateV <- pCombWiAge.df$cov %in% allFuncInd
pCombWiAge.df$corWiAge <- xCorrs
pCombWiAgeMelt <- data.frame(reshape2::melt(pCombWiAge.df, id.vars = c("cov", "relateV", "corWiAge")))
colnames(pCombWiAgeMelt)[5] <- "pval"
pCombWiAgeMelt$nlogpval <- -log10(pCombWiAgeMelt$pval)
pCombWiAgeMelt$ageAsCov <- T

pCombNoAge.df$cov <- 1:ncombs
pCombNoAge.df$relateV <- pCombNoAge.df$cov %in% allFuncInd
pCombNoAge.df$corWiAge <- xCorrs
pCombNoAgeMelt <- data.frame(reshape2::melt(pCombNoAge.df, id.vars = c("cov", "relateV", "corWiAge")))
colnames(pCombNoAgeMelt)[5] <- "pval"
pCombNoAgeMelt$nlogpval <- -log10(pCombNoAgeMelt$pval)
pCombNoAgeMelt$ageAsCov <- F

pCombAgeMelt <- rbind(pCombWiAgeMelt, pCombNoAgeMelt)
pCombAgeMelt$ageAsCov <- as.factor(pCombAgeMelt$ageAsCov)
levels(pCombAgeMelt$ageAsCov) <- c("BrainAGE ~ Covariates", "BrainAGE ~ Covariates + Age")
```



```{r}
myeffectSize <- formatC(mean(effectSizes), digit = 3, format="f")

p2 <- ggplot(pCombAgeMelt, 
            aes(group = cov, x = cov, y = nlogpval, fill = relateV)) + 
  facet_grid(~ ageAsCov) + 
    scale_fill_manual(values = cbPaletteGrey[-1]) +
  theme_bw() + labs(y = "-log(p)", x = "Covariates") +
  theme(legend.position = c(0.85, 0.85)) + 
  guides(
    fill = guide_legend(title="Influences imaging features", reverse = T),
    alpha = F) + 
  scale_x_continuous(breaks = seq(0, 16, 2))+
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(title = paste("Fold Change =", myeffectSize)) +
  geom_boxplot(alpha = 0.6) + geom_hline(yintercept = -log10(0.05))

print(p2)

# ggsave(file = paste0("Random", myeffectSize, "FoldChange", nruns, "runs.pdf"), plot = p2, width = 8, height = 5)
save(myeffectSize, pCombAgeMelt, nruns, file = paste0("Random", myeffectSize, "FoldChange", nruns, "runs.RData"))
```

```{r}
sessionInfo()
```