metaboliteUtility:

#########################################################################################################################
### Handle original data
filePath <- "C:\\Users\\LiuShouye\\OneDrive\\UAMS_work\\longitudal_data_analysis\\"
datafile <- paste0(filePath,"Delay_discounting.txt")
data_l <- read.table(datafile,header = T,sep = "\t")


#In get_power_para function names(data) is FamHistory_Id (1 for positive and 2 for negative) variable1, variable2, variable3

get_power_para <- function(data)
{
  data_sub_pos <- data[which(data$FamHistory_id==1),]
  data_sub_neg <- data[which(data$FamHistory_id==2),]
  data_sub_mu <- aggregate(data,by=list(data$FamHistory_id),FUN=mean,na.rm = T)
  data_sub_sd <- aggregate(data,by=list(data$FamHistory_id),FUN=sd,na.rm = T)
  data_pos_var <- cov(data_sub_pos[2:length(names(data))])
  data_neg_var <- cov(data_sub_neg[2:length(names(data))])
  n1 = dim(data_sub_pos)[1]
  n2 = dim(data_sub_neg)[1]
  
  ## calculate the pooled variance
  data_pool_matrix <- ((n1-1)*data_pos_var + (n2-1)*data_neg_var)/(n1+n2-2)
  #pool_sd <-c((sqrt(data_pool_matrix))[1,1],(sqrt(data_pool_matrix))[2,2],(sqrt(data_pool_matrix))[3,3])
  pool_sd <-diag(sqrt(data_pool_matrix))
  #data_cor <- cor(data[2:4])
  pool_cor <- solve(diag(pool_sd)) %*%  data_pool_matrix %*% solve(diag(pool_sd))
  
  output <- list (data_mu = data_sub_mu,data_sd = data_sub_sd,n1=n1,n2=n2,
                  pooled_matrix=data_pool_matrix,pool_sd=pool_sd,pool_cor= pool_cor)
  
  return(output)
}



###################################################################################################################
##### This part is used to calculate t-test data without considering different times
orig_data <- data_l[c(4,5:7)]
log_data <- data_l[c(4,8:11)]
orig_data_omit_miss <- na.omit(orig_data)
log_data_omit_miss <- na.omit(log_data)
# This is used to perform power analysis by using original data
orig_power_para <- get_power_para(orig_data_omit_miss)
log_power_para <- get_power_para(log_data_omit_miss)
orig_power_para
log_power_para