---
title: "Spacekime Analytics (Time Complexity and Inferential Uncertainty)"
subtitle: "Structured Big Data Analytics Case-Study (UKBB)"
author: "SOCR Team "
date: "`r format(Sys.time(),'%m/%d/%Y')`"
output: 
  html_document:
    theme: spacelab
    highlight: tango
    includes:
      before_body: TCIU_header.html
    toc: true
    number_sections: true
    toc_depth: 2
    toc_float:
      collapsed: false
      smooth_scroll: true
    code_folding: hide
---

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


```{r kime, message=F, warning=F}

# Get this figure: fig <- get_figure("MattSundquist", 4064)
# Get this figure's data: data <- get_figure("MattSundquist", 4064)$data
# Add data to this figure: p <- add_trace(p, x=c(4, 5), y=c(4, 5), kwargs=list(filename="Klein bottle", fileopt="extend"))
# Get y data of first trace: y1 <- get_figure("MattSundquist", 4064)$data[[1]]$y

library(plotly)
```

# Structured Big Data Analytics Case-Study

Next, we will look at another interesting example of a large structured tabular dataset. The goal remains the same - examine the effects of indexing complex data only using kime-order (time) and comparing the data representations as well as the subsequent data analytics. In this case-study, we will use the [UK Biobank (UKBB) data](https://www.ukbiobank.ac.uk/data-showcase).

A previous investigation [Predictive Big Data Analytics using the UK Biobank Data](https://doi.org/10.1038/s41598-019-41634-y), based on $7,614$ imaging, clinical, and phenotypic features and neuroimaging data of $9,914$ UKBB subjects reported the twenty most salient derived imaging biomarkers. By jointly representing and modeling the significant clinical and demographic variables along with specific salient neuroimaging features, the researchers predicted the presence and progression of depression and mental health of participating volunteers. We will explore the effects of kime-direction on the findings based on the same data and methods. For ease of demonstration, efficient calculations, and direct interpretation, we start by transforming the data into a tighter computable object of dimensions $9,914\times 107$.

```{r message=F, warning=F}
UKBB_data <- get(load("E:/Ivo.dir/Research/UMichigan/Publications_Books/2019/DataScience_Book_Value_Uncertainty_Kime_2019/other/UKBB_data_cluster_label.Rdata")) 
# str(UKBB_data)
UKBB_Colnames <- colnames(UKBB_data); View(UKBB_Colnames); dim(UKBB_data)   # 9914 7615

# Extract the top-50 derived NI biomarkers (data_summary_cluster_2.xlsx), per
# https://drive.google.com/drive/folders/1SdAtefp_taabNL70JvwJZSexTkzXiEKD 
top50_NI_Biomarkers <- c("lh_BA_exvivo_area__lh_WhiteSurfArea_area", "rh_BA_exvivo_area__rh_WhiteSurfArea_area", "rh_aparc.a2009s_area__rh_WhiteSurfArea_area", "rh_aparc_area__rh_WhiteSurfArea_area", "lh_aparc_area__lh_WhiteSurfArea_area", "lh_aparc.a2009s_area__lh_WhiteSurfArea_area", "aseg__SupraTentorialVol", "aseg__SupraTentorialVolNotVent", "aseg__SupraTentorialVolNotVentVox", "aseg__BrainSegVol", "aseg__BrainSegVolNotVentSurf", "aseg__BrainSegVolNotVent", "aseg__CortexVol", "aseg__rhCortexVol", "aseg__lhCortexVol", "aseg__TotalGrayVol", "aseg__MaskVol", "rh_aparc.DKTatlas_area__rh_superiortemporal_area", "rh_aparc.DKTatlas_area__rh_superiorfrontal_area", "lh_aparc.DKTatlas_area__lh_superiorfrontal_area", "lh_aparc.DKTatlas_area__lh_lateralorbitofrontal_area", "lh_aparc.DKTatlas_area__lh_superiortemporal_area", "aseg__EstimatedTotalIntraCranialVol", "lh_aparc_area__lh_lateralorbitofrontal_area", "rh_aparc.DKTatlas_area__rh_lateralorbitofrontal_area", "lh_aparc_area__lh_superiorfrontal_area", "rh_aparc_area__rh_superiortemporal_area", "rh_aparc.a2009s_area__rh_G.S_cingul.Ant_area", "rh_aparc_area__rh_superiorfrontal_area", "lh_aparc_area__lh_rostralmiddlefrontal_area", "wmparc__wm.lh.lateralorbitofrontal", "wmparc__wm.lh.insula", "rh_aparc_area__rh_medialorbitofrontal_area", "lh_BA_exvivo_area__lh_BA3b_exvivo_area", "lh_aparc.DKTatlas_area__lh_postcentral_area", "lh_aparc.DKTatlas_volume__lh_lateralorbitofrontal_volume", "lh_aparc.DKTatlas_area__lh_insula_area", "aseg__SubCortGrayVol", "lh_aparc.a2009s_area__lh_G_orbital_area", "lh_aparc_area__lh_superiortemporal_area", "rh_aparc.DKTatlas_area__rh_insula_area", "lh_aparc.DKTatlas_area__lh_precentral_area", "lh_aparc.pial_area__lh_lateralorbitofrontal_area", "lh_aparc.DKTatlas_area__lh_rostralmiddlefrontal_area", "lh_aparc_area__lh_postcentral_area", "lh_aparc.pial_area__lh_superiorfrontal_area", "rh_aparc_area__rh_rostralmiddlefrontal_area", "wmparc__wm.lh.superiortemporal", "lh_aparc.pial_area__lh_rostralmiddlefrontal_area", "rh_aparc.DKTatlas_volume__rh_lateralorbitofrontal_volume")

# Extract the main clinical features (binary/dichotomous and categorical/polytomous)
#### binary
top25_BinaryClinical_Biomarkers <- c("X1200.0.0", "X1200.2.0", "X1170.0.0", "X1190.2.0", "X1170.2.0", "X2080.0.0", "X6138.2.2", "X20117.0.0", "X6138.0.2", "X2877.0.0", "X20117.2.0", "X2877.2.0", "X1190.0.0", "X4968.2.0", "X1249.2.0", "X1190.1.0", "X1170.1.0", "X2080.2.0", "X4292.2.0", "X2050.0.0", "X1628.0.0", "X1200.1.0", "X20018.2.0", "X4292.0.0", "X3446.0.0")
#### polytomous
top31_PolytomousClinical_Biomarkers <- c("X31.0.0", "X22001.0.0", "X1950.0.0", "X1950.2.0", "X1980.0.0", "X2040.2.0", "X1980.2.0", "X2030.0.0", "X2090.0.0", "X2040.0.0", "X1618.2.0", "X1618.0.0", "X1210.0.0", "X2030.2.0", "X2000.0.0", "X1930.0.0", "X2090.2.0", "X2000.2.0", "X1210.2.0", "X1618.1.0", "X4653.2.0", "X1970.2.0", "X1970.0.0", "X1980.1.0", "X1930.2.0", "X4598.2.0", "X4598.0.0", "X4653.0.0", "X2090.1.0", "X2040.1.0", "X4631.2.0")

# Extract derived computed phenotype 
derivedComputedPhenotype <- UKBB_Colnames[length(UKBB_Colnames)]

# Construct the Computable data object including all salient predictors and derived cluster phenotype
ColNameList <- c(top50_NI_Biomarkers, top25_BinaryClinical_Biomarkers, 
                     top31_PolytomousClinical_Biomarkers, derivedComputedPhenotype); length(ColNameList)
col.index <- which(colnames(UKBB_data) %in% ColNameList); length(col.index)
tight107_UKBB_data <- UKBB_data[ , col.index]; dim(tight107_UKBB_data)
## View(tight107_UKBB_data[1:10, ])   # Confirm the tight data object organization 
```

We can investigate the effects of the *kime-phase* on the resulting data analytic inference obtained using the UKBB data.

```{r}
library(caret)
################################
#  1. First deal with the missing values  # summary(tight107_UKBB_data)
##### "mi" imputation - LONG
#library(mi) # use pmm (predictive mean matching) imputation method for the missing variables.
#mdf <- missing_data.frame(tight107_UKBB_data)
#show(mdf)   # ; head(mdf); dim(mdf)
#options(mc.cores = 4)
#imputations <- mi(mdf, n.iter=4, n.chains=1, verbose=T)
#imp.data.frames <- complete(imputations, 1)
#imp_tight107_UKBB_data <- imp.data.frames[[1]]

# split-off the computed phenotype (avoid manipulationg hte outcome that will be predicted during analysis phase
colnames(tight107_UKBB_data)[length(tight107_UKBB_data)] <- "cluster_2_cluster"
colnames(tight107_UKBB_data)[length(tight107_UKBB_data)]   # fix the problem with $ in variable name
y_pheno  <- tight107_UKBB_data[ , length(tight107_UKBB_data)]
tight106_UKBB_data <- tight107_UKBB_data[, -length(tight107_UKBB_data)]; dim(tight106_UKBB_data)

# MICE imputation (With parallel core computing)
# all predictors with absolute correlation over 0.4 AND at least 30% usable cases
library(mice)
library(foreach)
library(doParallel)
# set-up local parallel cluster
number_cores <- detectCores() - 2
clustUKBB <- makeCluster(number_cores)
clusterSetRNGStream(clustUKBB, 1234)
registerDoParallel(clustUKBB)

imp_tight106_UKBB_data <-
  foreach(no = 1:number_cores, 
          .combine = ibind, 
          .export = "tight106_UKBB_data",
          .packages = "mice") %dopar%
{
  mice(tight106_UKBB_data, m=2,maxit=3, printFlag=T, seed=1234, method = 'cart')
}

str(imp_tight106_UKBB_data)

# save(imp_tight106_UKBB_data, file = "C:/Users/Dinov/Desktop/imp_tight106_UKBB_data.Rdata")
# test <- load(file = "C:/Users/Dinov/Desktop/imp_tight106_UKBB_data.Rdata")

# imp_tight106_UKBB_data <- mice(tight106_UKBB_data, m=2,maxit=3, printFlag=T, seed=1234, method = 'cart')
### MICE NOTE: When the column features have a number of unbalanced factors, these categorical variables are transformed into dummy indicator-variables.
### There is a high probability that the resulting columns may be linear combinations of one-another. 
### The default MICE imputation methods, linear regression, may not be able to solve the linear matrix equations due to matrix low rank, i.e., matrices may cannot be inverted causing errors like "system is computationally singular". RTo solve that problem, we can change the method to "cart". 
### Also use seed before/during imputation for reproducible results.
comp_imp_tight106_UKBB_data <- as.matrix(complete(imp_tight106_UKBB_data), 
                                         dimnames = list(NULL, colnames(tight106_UKBB_data)))

################################################
##  2. Split the UKBB into 9 epochs
# First normalize all 50 derived NI biomarkers *using scale* as the NI biomarkers have vastly different distributions
for (i in 1:50) {
  comp_imp_tight106_UKBB_data[, i] <- scale(comp_imp_tight106_UKBB_data[, i])
}

# Next configure the 11 epochs
comp_imp_tight106_UKBB_data <- comp_imp_tight106_UKBB_data[1:9900, ] # remove the last 14 cases to make the 11 epochs of size 900 observations yeach
is.matrix(comp_imp_tight106_UKBB_data); dim(comp_imp_tight106_UKBB_data)
dim(comp_imp_tight106_UKBB_data) <- c(11, 900, dim(comp_imp_tight106_UKBB_data)[2])
dim(comp_imp_tight106_UKBB_data)
# double check epoch split worked 
# identical(comp_imp_tight106_UKBB_data[10, 900, 12], as.matrix(complete(imp_tight106_UKBB_data))[10*900, 12])
epochs_tight106_UKBB_data_1 <- comp_imp_tight106_UKBB_data[1, , ]; dim(epochs_tight106_UKBB_data_1)

# 3. Transform all 9 epochs (Big datasets/signals) to k-space (Fourier domain)
x1 <- c(1:900)
FT_epochs_tight106_UKBB <- array(complex(), c(11, 900, dim(comp_imp_tight106_UKBB_data)[3]))
mag_FT_epochs_tight106_UKBB <- array(complex(), c(11, 900, dim(comp_imp_tight106_UKBB_data)[3]))
phase_FT_epochs_tight106_UKBB <- array(complex(), c(11, 900, dim(comp_imp_tight106_UKBB_data)[3]))
for (i in 1:11) {
  FT_epochs_tight106_UKBB[i, , ] <- fft(comp_imp_tight106_UKBB_data[i, , ])
  X2 <- FT_epochs_tight106_UKBB[i, , ]
  # plot(fftshift1D(log(Re(X2)+2)), main = "log(fftshift1D(Re(FFT(tight106_UKBB))))") 
  mag_FT_epochs_tight106_UKBB[i, , ] <- sqrt(Re(X2)^2+Im(X2)^2); 
  # plot(log(fftshift1D(Re(X2_mag))), main = "log(Magnitude(FFT(tight106_UKBB)))") 
  phase_FT_epochs_tight106_UKBB[i, , ] <- atan2(Im(X2), Re(X2)); 
  # plot(fftshift1D(X2_phase), main = "Shift(Phase(FFT(tight106_UKBB)))")
}

# Compute the Average Phase of all 11 epochs (this will be needed later to confirm better data analytics)
avgPhase_FT_epochs_tight106_UKBB <- apply(phase_FT_epochs_tight106_UKBB, c(2,3), mean)
dim(avgPhase_FT_epochs_tight106_UKBB)

### Test the process to confirm calculations
# X2<-FT_epochs_tight106_UKBB[1,,];X2_mag<-mag_FT_epochs_tight106_UKBB[1,,];X2_phase<-phase_FT_epochs_tight106_UKBB[1,,]
# Real2 = X2_mag * cos(X2_phase)
# Imaginary2 = X2_mag * sin(X2_phase)
# man_hat_X2 = Re(fft(Real2 + 1i*Imaginary2, inverse = T)/length(X2))
# ifelse(abs(man_hat_X2[5,10] - comp_imp_tight106_UKBB_data[1, 5, 10]) < 0.001, "Perfect Syntesis", "Problems!!!")
#######

##### 4. Invert back to spacetime the epochs_tight106_UKBB_data_1 signal with NIL and Average phase
# Start with Nil Phase
Real = mag_FT_epochs_tight106_UKBB[1, , ] * cos(0)  # cos(mag_FT_epochs_tight106_UKBB[1, , ])
Imaginary = mag_FT_epochs_tight106_UKBB[1, , ] * sin(0)   # sin(mag_FT_epochs_tight106_UKBB[1, , ])
ift_NilPhase_X2mag = Re(fft(Real+1i*Imaginary, inverse = T)/length(mag_FT_epochs_tight106_UKBB[1,,]))
# display(ift_NilPhase_X2mag, method = "raster")
# dim(ift_NilPhase_X2mag); View(ift_NilPhase_X2mag); # compare to View(comp_imp_tight106_UKBB_data[1, , ])

# Next, synthesize the data using Average Phase
Real_Avg = mag_FT_epochs_tight106_UKBB[1, , ] * cos(avgPhase_FT_epochs_tight106_UKBB)
Imaginary_Avg = mag_FT_epochs_tight106_UKBB[1, , ] * sin(avgPhase_FT_epochs_tight106_UKBB)
ift_AvgPhase_X2mag = Re(fft(Real_Avg+1i*Imaginary_Avg, inverse = T)/length(mag_FT_epochs_tight106_UKBB[1,,]))
# is.complex(ift_AvgPhase_X2mag); dim(ift_AvgPhase_X2mag)

# To Transform the entire UKBB data to k-space (Fourier domain)
# library(EBImage)
#FT_UKBB_data <- fft(comp_imp_tight106_UKBB_data)
#X2 <- FT_UKBB_data  # display(FT_UKBB_data, method = "raster") 
#mag_FT_UKBB_data <- sqrt(Re(X2)^2+Im(X2)^2) 
###  # plot(log(fftshift1D(Re(X2_mag))), main = "log(Magnitude(FFT(timeseries)))") 
#phase_FT_UKBB_data <- atan2(Im(X2), Re(X2)) 
### Test the process to confirm calculations
# X2<-FT_UKBB_data; X2_mag <- mag_FT_UKBB_data; X2_phase<-phase_FT_UKBB_data
# Real2 = X2_mag * cos(X2_phase)
# Imaginary2 = X2_mag * sin(X2_phase)
# man_hat_X2 = Re(fft(Real2 + 1i*Imaginary2, inverse = T)/length(X2))
# ifelse(abs(man_hat_X2[5,10] - comp_imp_tight106_UKBB_data[5, 10]) < 0.001, "Perfect Syntesis", "Problems!!!")
#######
# Then we can Invert back the complete UKBB FT data into spacetime using nil phase
#Real = mag_FT_UKBB_data * cos(0)  # cos(phase_FT_UKBB_data)
#Imaginary = mag_FT_UKBB_data * sin(0)   # sin(phase_FT_UKBB_data)
#ift_NilPhase_X2mag = Re(fft(Real+1i*Imaginary, inverse = T)/length(FT_UKBB_data))
# display(ift_NilPhase_X2mag, method = "raster")
# dim(ift_NilPhase_X2mag); View(ift_NilPhase_X2mag); # compare to View(aqi_data1)
#summary(comp_imp_tight106_UKBB_data); summary(ift_NilPhase_X2mag)

# 5. Epoch 1: Perform Random Forest prediction (based on ift_TruePhase_X2mag==Original==epochs_tight106_UKBB_data_1) of:
##### Ever depressed for a whole week 1 #########################################
library("randomForest")
# y_pheno <- comp_imp_tight106_UKBB_data[,"X4598.2.0"] ### Ever depressed for a whole week 1
# y_pheno <- as.factor(y_pheno)
colnames(epochs_tight106_UKBB_data_1) <- colnames(tight106_UKBB_data)
set.seed(1234)
rf_depressed <- randomForest(as.factor(epochs_tight106_UKBB_data_1[ ,"X4598.2.0"]) ~ . , 
          data=epochs_tight106_UKBB_data_1[ , !(colnames(epochs_tight106_UKBB_data_1) %in% c("X4598.0.0", "X4598.2.0"))])
rf_depressed
pred1 = predict(rf_depressed, type="class")
confusionMatrix(pred1, as.factor(epochs_tight106_UKBB_data_1[ ,"X4598.2.0"]))
############## Accuracy : 0.7911 #############################

##### plot a simple decision tree
library(rpart); library(rpart.plot)
## there is an error in installing the rattle package
library(rattle)  ## this is for the fancyRpartPlot
label <- as.factor(epochs_tight106_UKBB_data_1[ ,"X4598.2.0"])
data1 <- as.data.frame(epochs_tight106_UKBB_data_1[ , !(colnames(epochs_tight106_UKBB_data_1) %in%
                      c("X4598.0.0", "X4598.2.0", "cluster_2_cluster"))])
data1$label <- label
depress_tree <- rpart(label ~ ., control=rpart.control(minsplit=30, cp=0.001, maxdepth=30), data=data1) 
rpart.plot(depress_tree, type = 4, extra = 1, clip.right.labs = F, tweak=4)
pred1 = predict(depress_tree, type="class")
# table(pred1, label)
library(caret)
confusionMatrix(pred1, label)
########################### Accuracy : 0.8511  ##############################

### Prune decision tree
prune_depress_tree <- prune(depress_tree,
                            cp=depress_tree$cptable[which.min(depress_tree$cptable[,"xerror"]),"CP"])

# plot the pruned tree
plot(prune_depress_tree, uniform=TRUE, main="Pruned UKBB Decision Tree")
text(prune_depress_tree, use.n=TRUE, all=TRUE, cex=.8)
rpart.plot(prune_depress_tree, type = 4, extra = 1, clip.right.labs = F, tweak=2)
pred2 = predict(prune_depress_tree, type="class")
confusionMatrix(pred2, label)
######################### Accuracy : 0.7867 ###############################


###### 6. FOR the COMPLETE UKBB data: plot a simple decision tree
#library(rpart); library(rpart.plot)
#library(rattle)  ## this is for the fancyRpartPlot
dim(comp_imp_tight106_UKBB_data) <- c(11*900, dim(comp_imp_tight106_UKBB_data)[3])
colnames(comp_imp_tight106_UKBB_data) <- colnames(tight106_UKBB_data)
label3 <- as.factor(comp_imp_tight106_UKBB_data[ ,"X4598.2.0"])
data3 <- as.data.frame(comp_imp_tight106_UKBB_data[ , !(colnames(comp_imp_tight106_UKBB_data) %in%
                      c("X4598.0.0", "X4598.2.0", "cluster_2_cluster"))])
depress_tree3 <- rpart(label3 ~ ., control=rpart.control(minsplit=30, cp=0.001, maxdepth=30), data=data3) 
rpart.plot(depress_tree3, type = 4, extra = 1, clip.right.labs = F, tweak=2)
pred3 = predict(depress_tree3, type="class")
confusionMatrix(pred3, label3)
#################### Accuracy : 0.8238 ###########################

### Prune decision tree
prune_depress_tree3 <- prune(depress_tree3,
                            cp=depress_tree3$cptable[which.min(depress_tree3$cptable[,"xerror"]),"CP"])
rpart.plot(prune_depress_tree3, type = 4, extra = 1, clip.right.labs = F, tweak=2,  main="Pruned UKBB Decision Tree")
pred3 = predict(prune_depress_tree3, type="class")
confusionMatrix(pred3, label3)
#################### Accuracy : 0.8132 ###########################

# RF
rf_depressed_Complete <- randomForest(label3 ~ . ,  data=data3)
rf_depressed_Complete
pred3 = predict(rf_depressed_Complete, type="class")
confusionMatrix(pred3, label3)
# In random forests, there is no need for cross-validation or a separate test set to get an unbiased estimate of the test set error. The out-of-bag (oob) error estimate is internally computed during the run... In particular, when newdata is not provided, the *predict.randomForest()==predict()* method automatically returns the out-of-bag prediction.
##################### Accuracy : 0.808 #########################


# 7. Perform Random Forest prediction (based on ift_AvgPhase_X2mag) of:
##### Ever depressed for a whole week 1 #########################################
# library("randomForest")
# y_pheno <- comp_imp_tight106_UKBB_data[,"X4598.2.0"] ### Ever depressed for a whole week 1
# y_pheno <- as.factor(y_pheno)
label2 <- as.factor(epochs_tight106_UKBB_data_1[ ,"X4598.2.0"]) 
# use the real outcome not synthesized (as.factor(epochs_tight106_UKBB_data_1[ ,"X4598.2.0"]))
colnames(ift_AvgPhase_X2mag) <- colnames(tight106_UKBB_data)
data2 <- as.data.frame(ift_AvgPhase_X2mag[ , !(colnames(ift_AvgPhase_X2mag) %in%
                      c("X4598.0.0", "X4598.2.0", "cluster_2_cluster"))])
set.seed(1234)
rf_depressed_AvgPhase <- randomForest(label2 ~ . , importance=T, nodesize=30, mtry=100, ntree= 1000, data=data2)
rf_depressed_AvgPhase
pred1_AvgPhase <- predict(rf_depressed_AvgPhase, type="class")
confusionMatrix(pred1_AvgPhase, label2)

##### plot a simple decision tree
# library(rpart); library(rpart.plot)
## there is an error in installing the rattle package
# library(rattle)  ## this is for the fancyRpartPlot
set.seed(1234)
depress_tree_AvgPhase <- rpart(label2 ~ ., control=rpart.control(minsplit=30,cp=0.001,maxdepth=30),data=data2) 
pred1_AvgPhase <- predict(depress_tree_AvgPhase, type="class")
# library(caret)
confusionMatrix(pred1_AvgPhase, label2)
rpart.plot(depress_tree_AvgPhase, type = 2, extra = 1, clip.right.labs = F, varlen=5, faclen=5, tweak=4)
################### Accuracy : 0.7956 #############################

### Prune decision tree
prune_depress_tree_AvgPhase <- prune(depress_tree_AvgPhase,
        cp=depress_tree_AvgPhase$cptable[which.min(depress_tree_AvgPhase$cptable[,"xerror"]),"CP"]/(1.5))
# plot the pruned tree
#plot(prune_depress_tree_AvgPhase, uniform=TRUE, main="Pruned UKBB Decision Tree (Avg-Phase Synthesis)")
#text(prune_depress_tree_AvgPhase, use.n=TRUE, all=TRUE, cex=.8)
rpart.plot(prune_depress_tree_AvgPhase, type = 4, extra = 1, clip.right.labs = F, varlen=5, faclen=5, tweak=1.5)
pred2_AvgPhase = predict(prune_depress_tree_AvgPhase, type="class")
confusionMatrix(pred2_AvgPhase, label2)



# 8. Perform Random Forest prediction (based on ift_NilPhase_X2mag) of:
##### Ever depressed for a whole week 1 #########################################
# library("randomForest")
# y_pheno <- comp_imp_tight106_UKBB_data[,"X4598.2.0"] ### Ever depressed for a whole week 1
# y_pheno <- as.factor(y_pheno)
label2 <- as.factor(epochs_tight106_UKBB_data_1[ ,"X4598.2.0"]) 
# use the real outcome not synthesized (as.factor(epochs_tight106_UKBB_data_1[ ,"X4598.2.0"]))
colnames(ift_NilPhase_X2mag) <- colnames(tight106_UKBB_data)
data2 <- as.data.frame(ift_NilPhase_X2mag[ , !(colnames(ift_NilPhase_X2mag) %in%
                      c("X4598.0.0", "X4598.2.0", "cluster_2_cluster"))])
set.seed(1234)
rf_depressed_NilPhase <- randomForest(label2 ~ . , importance=T, nodesize=30, mtry=100, ntree= 1000, data=data2)
rf_depressed_NilPhase
pred1_NilPhase <- predict(rf_depressed_NilPhase, type="class")
confusionMatrix(pred1_NilPhase, label2)

##### plot a simple decision tree
# library(rpart); library(rpart.plot)
## there is an error in installing the rattle package
# library(rattle)  ## this is for the fancyRpartPlot
set.seed(1234)
depress_tree_NilPhase <- rpart(label2 ~ ., control=rpart.control(minsplit=30,cp=0.001,maxdepth=30),data=data2) 
pred1_NilPhase <- predict(depress_tree_NilPhase, type="class")
# library(caret)
confusionMatrix(pred1_NilPhase, label2)
rpart.plot(depress_tree_NilPhase, type = 2, extra = 1, clip.right.labs = F, varlen=5, faclen=5, tweak=2)
################### Accuracy : 0.79  #############################

### Prune decision tree
prune_depress_tree_NilPhase <- prune(depress_tree_NilPhase,
            cp=depress_tree_NilPhase$cptable[which.min(depress_tree_NilPhase$cptable[,"xerror"]),"CP"]/(1.5))
# plot the pruned tree
#plot(prune_depress_tree_NilPhase, uniform=TRUE, main="Pruned UKBB Decision Tree (Nil-Phase Synthesis)")
#text(prune_depress_tree_NilPhase, use.n=TRUE, all=TRUE, cex=.8)
rpart.plot(prune_depress_tree_NilPhase, type = 4, extra = 1, clip.right.labs = F, varlen=5, faclen=5, tweak=2)
pred2_NilPhase = predict(prune_depress_tree_AvgPhase, type="class")
confusionMatrix(pred2_NilPhase, label2)


# 9. Compare the analytics results from #3, ..., #8!
# Compare the original against nil-phase and avg-phase synthesized data
origNilNil_6rows_Compare <- rbind(head(epochs_tight106_UKBB_data_1), head(ift_NilPhase_X2mag), head(ift_AvgPhase_X2mag))
x1 <- 1:dim(epochs_tight106_UKBB_data_1)[2]
orig.mean <- apply(epochs_tight106_UKBB_data_1, 2, mean)
nilPhase.mean <- apply(ift_NilPhase_X2mag, 2, mean)
avgPhase.mean <- apply(ift_AvgPhase_X2mag, 2, mean)

plot(x1, orig.mean, main = "Comparing original UKBB against nil-phase and avg-phase synthesized data",
     col="green", lwd = 3, type="l", lty=1, xlab = "Features", ylab = "Averages across cases")
lines(x1, nilPhase.mean, col = "red", lwd = 3, lty=1)
lines(x1, avgPhase.mean, col = "blue", lwd = 3, lty=1)
legend("top", bty="n", legend=c(
  sprintf("Original"), sprintf("Nil-Phase Reconstruction"), 
  sprintf("Average-Phase Reconstruction")), 
  col=c("green", "red", "blue"), lty=c(1,1,1), lwd=c(3,3,3), cex=0.9)

stopCluster(clustUKBB)

################### Plot_ly Figure 6.15
time <- x1[1:length(orig.mean)]
p <- plot_ly()
p <- add_lines(p, x=~time,  y=~orig.mean, 
               name = "Original.Mean", type = 'scatter', mode = 'lines', 
               hoverinfo = 'name', line=list(color='green', width=4))
p <- add_lines(p, x=~time, y=~nilPhase.mean, 
               name = "ift_NilPhase.Mean", type = 'scatter', mode = 'lines', 
               hoverinfo = 'name', 
               line=list(color='red', width = 2))
p <- add_lines(p, x=~time, y=~avgPhase.mean, 
               name = "ift_AvgPhase.Mean", type = 'scatter', mode = 'lines', 
               hoverinfo = 'name', 
               line=list(color='blue', width = 2)) %>%
  layout(xaxis = list(range = c(0,108), title="Features"),
    yaxis = list(range = c(0,6), title="Averages across cases"),
    title=sprintf('Comparison:  Corr(Real, NilPhase) = %s; Comparison:  Corr(Real, AvgPhase) = %s',
                  format(cor(orig.mean, nilPhase.mean), digits=3),
                  format(cor(orig.mean, avgPhase.mean), digits=3)),
    legend = list(orientation = "h",   # show entries horizontally
                     xanchor = "center",  # use center of legend as anchor
                     x = 0.5, y=-0.1,
                     size = 30))      
p

```

As we have the benefit of hind-side, the prior study [Predictive Big Data Analytics using the UK Biobank Data](https://doi.org/10.1038/s41598-019-41634-y) had identified the most salient derived neuroimaging and clinical biomarkers ($k=107$) in this study, including the physician identified clinical outcomes and the computed phenotypes using unsupervised machine learning methods. Therefore, for simplicity, this demonstration only focuses on these features, without a previous feature-selection preprocessing step. We are examining the effects of the kime phase on the scientific inference.

Note that averaging the phases in the Fourier domain assumes that either we have a very large number of samples that effectively span the range of kime-angles, or that the sampling is uniform on the kime-angles space. When these assumptions are violated, other phase-aggregation or phase-ensembling methods (e.g., weighted mean, non-parametric measures of centrality, or Bayesian strategies) may need to utilized to ensure the reliability of the final inference.



<!--html_preserve-->
<div>
    	<footer><center>
			<a href="https://www.socr.umich.edu/">SOCR Resource</a>
				Visitor number <img class="statcounter" src="https://c.statcounter.com/5714596/0/038e9ac4/0/" alt="Web Analytics" align="middle" border="0">
				<script type="text/javascript">
					var d = new Date();
					document.write(" | " + d.getFullYear() + " | ");
				</script> 
				<a href="https://socr.umich.edu/img/SOCR_Email.png"><img alt="SOCR Email"
	 			title="SOCR Email" src="https://socr.umich.edu/img/SOCR_Email.png"
	 			style="border: 0px solid ;"></a>
	 		 </center>
	 	</footer>

	<!-- Start of StatCounter Code -->
		<script type="text/javascript">
			var sc_project=5714596; 
			var sc_invisible=1; 
			var sc_partition=71; 
			var sc_click_stat=1; 
			var sc_security="038e9ac4"; 
		</script>
		
		<script type="text/javascript" src="https://www.statcounter.com/counter/counter.js"></script>
	<!-- End of StatCounter Code -->
	
	<!-- GoogleAnalytics -->
		<script src="https://www.google-analytics.com/urchin.js" type="text/javascript"> </script>
		<script type="text/javascript"> _uacct = "UA-676559-1"; urchinTracker(); </script>
	<!-- End of GoogleAnalytics Code -->
</div>
<!--/html_preserve-->
