text: #!/usr/bin/env python
  # encoding:utf-8
  """
  Author : Yuanqing Mei
  Date : 2021/4/4
  Time: 9:29
  File: AlvesMeta.py
  HomePage : http://github.com/yuanqingmei
  Email : dg1533019@smail.nju.edu.cn
  
  Deriving Alves's method threshold by meta-analysis.
  computes the variance of high level threshold: var(T)=var(cX)=c^2×var(X)=[E(T)/E(X)]^2×var(X)
  """
  
  
  def Alves_meta(t_dir="F:\\NJU\\MTmeta\\experiments\\unsupervised\\trainingData\\",
                m_dir="F:\\NJU\\MTmeta\\experiments\\unsupervised\\AlvesThresholdMinified\\"):
  
      import os
      import csv
      from scipy.stats import norm  # norm.cdf() the cumulative normal distribution function in Python
      from scipy import stats  # 根据卡方分布计算p值: p_value=1.0-stats.chi2.cdf(chisquare,freedom_degree)
      import numpy as np
      import pandas as pd
  
      # 输入：两个匿名数组，effect_size中存放每个study的effect size，variance存放对应的方差
      # 输出：fixed effect model固定效应元分析后的结果，包括
      #      (1)fixedMean：固定效应元分析后得到的效应平均值；(2) fixedStdError：固定效应元分析的效应平均值对应的标准错
      def fixed_effect_meta_analysis(effect_size, variance):
          fixed_weight = []
          sum_Wi = 0
          sum_WiYi = 0
          d = {}  # return a dict
          study_number = len(variance)
          for i in range(study_number):
              if variance[i] == 0:
                  continue
              fixed_weight.append(1 / variance[i])
              sum_Wi = sum_Wi + fixed_weight[i]
              sum_WiYi = sum_WiYi + effect_size[i] * fixed_weight[i]
          fixedMean = sum_WiYi / sum_Wi  # 固定模型元分析后得到的效应平均值
          fixedStdError = (1 / sum_Wi) ** 0.5  # 固定模型元分析的效应平均值对应的标准错
          d['fixedMean'] = fixedMean
          d['fixedStdError'] = fixedStdError
          return d
  
      # 输入：两个匿名数组，effect_size中存放每个study的effect size，variance存放对应的方差
      # 输出：random effect model随机效应元分析后的结果，包括：
      #      (1) randomMean：随机模型元分析后得到的效应平均值； (2) randomStdError：随机模型元分析的效应平均值对应的标准错
      def random_effect_meta_analysis(effect_size, variance):
  
          sum_Wi = 0
          sum_WiWi = 0
          sum_WiYi = 0  # Sum(Wi*Yi), where i ranges from 1 to k, and k is the number of studies
          sum_WiYiYi = 0  # Sum(Wi*Yi*Yi), where i ranges from 1 to k, and k is the number of studies
  
          sum_Wistar = 0
          sum_WistarYi = 0
          d = {}  # return a dict
  
          study_number = len(variance)
          fixed_weight = [0 for i in range(study_number)]  # 固定模型对应的权值
          random_weight = [0 for i in range(study_number)]  # 随机模型对应的权值
  
          for i in range(study_number):
              if variance[i] == 0:
                  continue
              fixed_weight[i] = 1 / variance[i]
              sum_Wi = sum_Wi + fixed_weight[i]
              sum_WiWi = sum_WiWi + fixed_weight[i] * fixed_weight[i]
              sum_WiYi = sum_WiYi + effect_size[i] * fixed_weight[i]
              sum_WiYiYi = sum_WiYiYi + fixed_weight[i] * effect_size[i] * effect_size[i]
  
          Q = sum_WiYiYi - sum_WiYi * sum_WiYi / sum_Wi
          df = study_number - 1
          C = sum_Wi - sum_WiWi / sum_Wi
  
          # 当元分析过程中只有一个study研究时，没有研究间效应，故研究间的方差为零
          if study_number == 1:
              T2 = 0
          else:
              T2 = (Q - df) / C  # sample estimate of tau square
  
          if T2 < 0:
              T2 = 0  # 20210411，若T2小于0，取0,   M.Borenstein[2009] P114
  
          for i in range(study_number):
              random_weight[i] = 1 / (variance[i] + T2)  # random_weight 随机模型对应的权值
  
          for i in range(study_number):
              sum_Wistar = sum_Wistar + random_weight[i]
              sum_WistarYi = sum_WistarYi + random_weight[i] * effect_size[i]
  
          randomMean = sum_WistarYi / sum_Wistar  # 随机模型元分析后得到的效应平均值
          randomStdError = (1 / sum_Wistar) ** 0.5  # 随机模型元分析的效应平均值对应的标准错
          # 当元分析过程中只有一个study研究时，没有研究间异质性，故异质性为零
          if study_number == 1:
              I2 = 0
          else:
              I2 = ((Q - df) / Q) * 100  # Higgins et al. (2003) proposed using a statistic, I2,
              # the proportion of the observed variance reflects real differences in effect size
  
          if I2 < 0:
              I2 = 0        # 20210418，若I2小于0，取0,   M.Borenstein[2009] P110
  
          pValue_Q = 1.0 - stats.chi2.cdf(Q, df)  # pValue_Q = 1.0 - stats.chi2.cdf(chisquare, freedom_degree)
  
          d["C"] = C
          d["mean"] = randomMean
          d["stdError"] = randomStdError
          d["LL_CI"] = randomMean - 1.96 * randomStdError  # The 95% lower limits for the summary effect
          d["UL_CI"] = randomMean + 1.96 * randomStdError  # The 95% upper limits for the summary effect
          d["ZValue"] = randomMean / randomStdError  # a Z-value to test the null hypothesis that the mean effect is zero
          # 20210419 双侧检验时需要增加绝对值符号np.abs
          d["pValue_Z"] = 2 * (1 - norm.cdf(np.abs(randomMean / randomStdError)))  # norm.cdf() 返回标准正态累积分布函数值
          d["Q"] = Q
          d["df"] = df
          d["pValue_Q"] = pValue_Q
          d["I2"] = I2
          d["tau"] = T2 ** 0.5
          d["LL_ndPred"] = randomMean - 1.96 * (T2 ** 0.5)  # tau、randomMean 已知情况下的新出现的study的effctsize所落的区间
          d["UL_ndPred"] = randomMean + 1.96 * (T2 ** 0.5)  # tau、randomMean 已知情况下的新出现的study的effctsize所落的区间
          # tau、randomMean 未知情况（估计）下的新出现的study的effctsize所落的区间
          # stats.t.ppf(0.975,df)返回学生t分布单尾alpha=0.025区间点(双尾是alpha=0.05)的函数，它是stats.t.cdf()累积分布函数的逆函数
          d["LL_tdPred"] = randomMean - stats.t.ppf(0.975, df) * ((T2 + randomStdError * randomStdError) ** 0.5)
          # tau、randomMean 未知情况（估计）下的新出现的study的effctsize所落的区间
          d["UL_tdPred"] = randomMean + stats.t.ppf(0.975, df) * ((T2 + randomStdError * randomStdError) ** 0.5)
  
          fixedMean = sum_WiYi / sum_Wi  # 固定模型元分析后得到的效应平均值
          fixedStdError = (1 / sum_Wi) ** 0.5  # 固定模型元分析的效应平均值对应的标准错
          d['fixedMean'] = fixedMean
          d['fixedStdError'] = fixedStdError
          return d
  
      metric_dir = t_dir
      meta_dir = m_dir
      os.chdir(metric_dir)
      print(os.getcwd())
  
      # 显示所有行
      pd.set_option('display.max_rows', None)
  
      # 显示所有列
      pd.set_option('display.max_columns', None)
  
      # the item of row of dataframe
      pd.set_option('display.width', 5000)
  
      # read_csv(path, keep_default_na=False, na_values=[""])  只有一个空字段将被识别为NaN
      df = pd.read_csv(meta_dir + "Alves_Threshold_minified.csv", keep_default_na=False, na_values=[""])
  
      metric_names = df.columns.values.tolist()
  
      k = 0
  
      for metric in metric_names:
  
          # chooses the high level value as the threshold
          if metric[-5:] != "_High":
              continue
          # There are no more than 10 values of DCAEC and DCMEC metric, which are greater than 0. drop it
          if (metric[:-5] == "DCAEC") or (metric[:-5] == "DCMEC"):
              continue
  
          print("the current metric is ", metric)
          print("the current metric is ", metric[:-5])
  
          threshold_effect_size = df.loc[:, metric]
          # SPD、ACAIC和ACMIC度量的high水平上阈值全为零，故取very_high水平上的阈值
          if metric[:-5] == "ACAIC" or metric[:-5] == "ACMIC" or metric[:-5] == "SPD":
              threshold_effect_size = df.loc[:, metric[:-5] + "_Very-High"]
  
          threshold_variance = []
          for i in range(len(df["fileName"])):
              file_name = df.loc[i, "fileName"]
              df_metric = pd.read_csv(metric_dir + file_name + ".csv", keep_default_na=False, na_values=[""])
              mean_metric = df_metric[metric[:-5]].mean()
              variance_metric = df_metric[metric[:-5]].var()
              threshold_metric = df.loc[i, metric]
              if metric[:-5] == "ACAIC" or metric[:-5] == "ACMIC" or metric[:-5] == "SPD":
                  threshold_metric = df.loc[i, metric[:-5] + "_Very-High"]
              # print(mean_metric,variance_metric,threshold_metric,(threshold_metric/mean_metric) ** 2 * variance_metric)
              threshold_variance.append((threshold_metric/mean_metric) ** 2 * variance_metric)
  
          metaThreshold = pd.DataFrame()
          metaThreshold['EffectSize'] = threshold_effect_size
          metaThreshold['Variance'] = threshold_variance
          try:
              resultMetaAnalysis = random_effect_meta_analysis(
                  np.array(metaThreshold[metaThreshold["EffectSize"] > 0].loc[:, "EffectSize"]),
                  np.array(metaThreshold[metaThreshold["EffectSize"] > 0].loc[:, "Variance"]))
  
              with open(meta_dir + "Alves_metaThresholds.csv", 'a+', encoding="utf-8", newline='') as f:
                  writer_f = csv.writer(f)
                  if os.path.getsize(meta_dir + "Alves_metaThresholds.csv") == 0:
                      writer_f.writerow(["metric", "Alves_metaThreshold", "Alves_metaThreshold_stdError", "LL_CI", "UL_CI",
                                         "ZValue", "pValue_Z", "Q", "df", "pValue_Q", "I2", "tau"])
                  writer_f.writerow([metric, resultMetaAnalysis["mean"], resultMetaAnalysis["stdError"],
                                     resultMetaAnalysis["LL_CI"], resultMetaAnalysis["UL_CI"],
                                     resultMetaAnalysis["ZValue"], resultMetaAnalysis["pValue_Z"],
                                     resultMetaAnalysis["Q"], resultMetaAnalysis["df"], resultMetaAnalysis["pValue_Q"],
                                     resultMetaAnalysis["I2"], resultMetaAnalysis["tau"]])
  
          except Exception as err1:
              print(err1)
  
          k += 1
          # if k == 1:
          #     break
  
  
  if __name__ == '__main__':
      Alves_meta()
  
      print("The __name__ is ", __name__, ". This is end of AlvesMeta.py!")
  