// ancova.C Version 7.44  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-30 10:41:07 tac>
//
// 最初の項目を群項目、次の項目を集計項目Y、以後の項目を共変数Xとした共分散分析を行う。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W ancova.C
//              -I/usr/local/include -L/usr/local/lib -ldans -lfun

#include <dans.H>               // dans用ヘッダファイル(fun.H,mat.H含む)

using namespace dans;           // dans名前空間使用(名前空間matvec含む)

const vecstr USAGE =            // プログラム使用法表示用
{
  "ancova  Version 7.44         Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:ancova [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目を群項目、次の項目を集計項目Y、以後の項目を共変数Xとした",
  "     共分散分析をする。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定する。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -dnsname[=<出力ファイル名>]:DANSデータを出力するファイル名を指定する。",
  "    オプション省略時:DANSデータ出力せず  -dnsname:-dnsname=ancova_dns.d",
  "    出力項目:群別回帰係数",
  "  -f=<ファイル>:指定したファイルからオプションを読み込む。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示する。(-h:-h=man)",
  "  -res=[-][{def|all|atab|bt|gstat|m[{a|d}]{b|f|s|t}|pd",
  "      |tstat}][,…]:結果の出力項目を指定する。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    atab:共分散分析表  bt:回帰係数の検定  gstat:群別統計量",
  "    m[{a|d}]{b|f|s|t}:修正群差と非平行性の多重比較",
  "    pd:非平行性を残差にプールした分散分析表  tstat:全体の統計量",
  "    オプション省略時:自動判定",
  "    多重比較 b:Bonferroni型  f:Fisher型  s:Scheffe型  t:Tukey型"
};

const string Title = "=== 共分散分析(analysis of covariance, ANCOVA) ==="; // 表題
const string DNSNAMEdef = "ancova_dns.d"; // 省略時DANSファイル名

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  string v_dnsname;             // DANSファイル名
  mctype v_ma;                  // 修正群差の多重比較法
  mctype v_md;                  // 非平行性の多重比較法
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  double cc() const             // 信頼係数を返す
  { return v_cc; }
  const string& dnsname() const // DANSデータファイル名を返す
  { return v_dnsname; }
  mctype ma() const             // 修正群差の多重比較法を返す
  { return v_ma; }
  mctype md() const             // 非平行性の多重比較法を返す
  { return v_md; }
  bool res(const string& code); // 結果出力フラッグを返す
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check(size_t ni);         // オプションのチェックとセット
  int set_res();                // -resオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_cc(CCdef), v_dnsname(), v_ma(mctuk), v_md(mcbon), reslst(), resmap()
// コンストラクタ
{
  if (c_pname == nullptr) return; // 初期設定せず
  v_pname = getname(c_pname);     // パス名なしのプログラム名取得
  v_path = c_pname;
  strsize_t idx = v_path.find(v_pname);
  if (idx != string::npos) v_path.erase(idx); // プログラム名削除
}

bool option::res(const string& code) // 分析結果出力フラッグを返す
{
  mapsb_citr resitr = resmap.find(code);
  if (resitr == resmap.end()) return false; // code無し！
  return resmap[code];                      // true/falseを返す
}

void option::set(const string& optstr, char devc)
// オプション指定によってオプションを設定する
{
  string opt[2];                // オプションのopt[0]:項目,opt[1]:内容
  if (getopt(opt, optstr, devc) < 1) return; // オプション指定なし
  int err = 0;
  if (opt[0] == "cc")           // 信頼係数の指定
  {
    v_cc = opt[1].empty() ? CCdef : strtodf(opt[1]);
    if (v_cc < 0.0 || 1.0 <= v_cc) err = badopt(pname(), optstr, USAGE, 0.0, 1.0);
  }
  else if (opt[0] == "dnsname") // DANSファイル名指定
  {
    if (opt[1].empty()) v_dnsname = DNSNAMEdef;
    else v_dnsname = opt[1];
    if (v_dnsname.size() < 3) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "f")       // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else if (opt[0] == "res")     // 結果の出力指定
  {
    if (strsplit(reslst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else                          // 間違ったオプション指定
  { err = nonopt(pname(), optstr, USAGE); }
  if (err > 0) throw help_except("ヘルプ表示");
  if (err < 0) throw option_error(string("オプション指定違反→") + optstr);
}

int option::check(size_t ni)    // オプションのセットとチェック
{
  if (ni < 2) return disperr(pname(), "項目が2つ以上必要です！");
  if (set_res() < 0) return -3; // -resオプションのエラー
  return 0;                     // エラーなし
}

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["atab"] = resmap["gstat"] = resmap["tstat"] = true;
  resmap["bt"] = resmap["ma"] = resmap["md"] = resmap["pd"] = false;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["atab"] = resmap["gstat"] = resmap["tstat"] = false;
    for (const string& resstr : reslst)
    {
      mapsb_citr resitr = resmap.find(resstr);
      if (resstr == "-")        // 全て出力せず
      { resmap["res"] = false; }
      else if (resstr[0] == '-' && resstr.size() > 1) // 特定の項目を出力せず
      {
        resitr = resmap.find(resstr.substr(1));
        if (resitr == resmap.end())
        { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
        else
        { resmap[resstr.substr(1)] = false; }
      }
      else if (resstr == "def") // 省略時出力
      {
        resmap["def"] = resmap["atab"] = resmap["cc"] = resmap["gstat"] = resmap["tstat"] = true;
      }
      else if (resstr == "all") // 全ての結果を出力する
      {
        resmap["def"] = resmap["atab"] = resmap["gstat"] = resmap["tstat"] = resmap["bt"] = resmap["ma"] = resmap["md"] = resmap["pd"] = true;
      }
      else if (resstr == "m")      // -res=mと指定された時
      { resmap["ma"] = resmap["md"] = true; }
      else if (resstr[0] == 'm' && resstr[1] == 'a' && resstr.size() > 2)
      {
        switch (resstr[2])
        {
        case 'b' : v_ma = mcbon;
          break;
        case 'f' : v_ma = mcfis;
          break;
        case 's' : v_ma = mcsch;
          break;
        case 't' : v_ma = mctuk;
          break;
        default : return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！");
        }
        resmap["ma"] = true;
      }
      else if (resstr[0] == 'm' && resstr[1] == 'd' && resstr.size() > 2)
      {
        switch (resstr[2])
        {
        case 'b' : v_md = mcbon;
          break;
        case 'f' : v_md = mcfis;
          break;
        case 's' : v_md = mcsch;
          break;
        default : return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！");
        }
        resmap["md"] = true;
      }
      else if (resitr == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class gstat                     // 群別統計量クラス
{
private:
  mat<double> smat;             // 積和行列
  int v_n;                      // 例数
  vecdbl mevec;                 // xとyの平均ベクトル
  vecdbl ssvec;                 // xとyの平方和ベクトル
public:
  gstat() : smat(), v_n(0), mevec(), ssvec() // デフォールトコンストラクタ
  {}
  explicit gstat(size_t nd) : smat(nd), v_n(0), mevec(nd), ssvec(nd)
  // コンストラクタ
  {}
  gstat(const gstat& gstat0) : smat(gstat0.smat), v_n(gstat0.v_n), mevec(gstat0.mevec), ssvec(gstat0.ssvec)
  // コピーコンストラクタ
  {}
  virtual ~gstat()              // デストラクタ
  {}
  const gstat& operator=(const gstat& gstat0); // コピー代入演算子
  const gstat& operator+=(const gstat& gstat0); // 加算代入演算子
  void sigma(const vecdbl& fdata);              // データを累積する
  int setsmat(bool sets);       // 平均と積和行列を設定して例数を返す
  void calreg();                // 回帰式を計算する
  int n() const                 // 例数を返す
  { return v_n; }
  double me(size_t idx) const   // 平均値を返す
  { return (idx < mevec.size()) ? mevec[idx] : 0.0; }
  double my() const             // yの平均値を返す
  { return mevec[mevec.size() - 1]; }
  double ss(size_t idx) const   // 平方和を返す
  { return (idx < ssvec.size()) ? ssvec[idx] : 0.0; }
  double vr(size_t idx) const   // 分散を返す
  { return (idx < ssvec.size() && n() > 1) ? ssvec[idx] / (n() - 1) : 0.0; }
  double syy() const            // yの平方和を返す
  { return ssvec[ssvec.size() - 1]; }
  double vyy() const            // yの分散を返す
  { return (n() > 1) ? ssvec[ssvec.size() - 1] / (n() - 1) : 0.0; }
  double sr() const             // 回帰残差を返す
  { return smat(smat.r() - 1, smat.c() - 1); }
  double b(size_t idx) const    // 回帰係数を返す
  { return (idx < smat.r() - 1) ? smat(idx, smat.c() - 1) : 0.0; }
  double vbinv(size_t idx) const // 回帰係数の分散要素を返す
  { return (idx < smat.r() - 1) ? smat(idx, idx) : 0.0; }
  double vxr(const vec<double>& dxv) const; // 多重比較用dV^(-1)dを返す
};

const gstat& gstat::operator=(const gstat& gstat0) // コピー代入演算子
{
  if (&gstat0 == this) return *this; // 自分自身への代入を防ぐ
  smat = gstat0.smat;  v_n = gstat0.v_n;
  mevec = gstat0.mevec;  ssvec = gstat0.ssvec;
  return *this;
}

const gstat& gstat::operator+=(const gstat& gstat0) // 加算代入演算子
{
  smat += gstat0.smat;  v_n += gstat0.v_n;
  return *this;
}

void gstat::sigma(const vecdbl& fdatavec) // データを累積する
{
  mevec[0] = 1.0;
  for (size_t i = 0; i < fdatavec.size(); ++i) mevec[i + 1] = fdatavec[i];
  for (size_t i = 0; i <= fdatavec.size(); ++i)
  {
    double x = mevec[i];
    for (size_t j = i; j <= fdatavec.size(); ++j) smat.add(i, j, x * mevec[j]);
  }
}

int gstat::setsmat(bool sets)  // 平均と積和行列をセットして例数を返す
{
  if (sets)                     // 単純積和行列を積和行列にする
  {
    v_n = (int)(smat(0, 0));
    for (size_t i = 0; i < smat.r() - 1; ++i)
    { for (size_t j = i + 1; j < smat.c(); ++j) smat.assign(j, i, smat(i, j)); }
    smat.sweep(0);
  }
  for (size_t i = 0; i < smat.r(); ++i)
  { mevec[i] = smat(0, i);  ssvec[i] = smat(i, i); }
  return v_n;                   // 例数を返す
}

void gstat::calreg()            // 回帰式を計算する
{
  for (size_t i = 1; i < smat.r() - 1; ++i) smat.sweep(i);
}

double gstat::vxr(const vec<double>& dxv) const // 多重比較用dV^(-1)dを返す
{
  int p = dxv.dim();
  vec<double> tv = smat.sub(1, 1, p, p) * dxv;
  return inp(dxv, tv);
}

using mapsg = map<string, gstat>;  // 群別統計量マップの別名定義
using mapsg_itr = mapsg::iterator; // 反復子の別名定義
using mapsg_citr = mapsg::const_iterator; // const反復子の別名定義

class varclass                 // 変数クラス
{
private:
  int v_ni;                    // 項目数
  mapsg v_gsmap;               // 群統計量マップ
  gstat v_gstatall;            // 全体統計量オブジェクト
  gstat v_gstatcom;            // 共通統計量オブジェクト
  double v_sbc;                // 共通回帰平方和
  double v_sa;                 // 群差平方和
  double v_sbt;                // 全体回帰平方和
  double v_saa;                // 修正群差平方和
  double v_sd;                 // 非平行性平方和
  double v_sr;                 // 残差平方和
  void cal_ancova();           // 共分散分析用各種パラメターを計算する
public:
  varclass(size_t ni) : v_ni(ni), v_gsmap(), v_gstatall(ni), v_gstatcom(ni), v_sbc(0.0), v_sa(0.0), v_sbt(0.0), v_saa(0.0), v_sd(0.0), v_sr(0.0)
  // コンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  int ni() const                // 項目数を返す
  { return v_ni; }
  int nt() const                // 全例数を返す
  { return v_gstatall.n(); }
  size_t ng() const             // 群数を返す
  { return v_gsmap.size(); }
  const mapsg& gsmap() const    // 群統計量マップを返す
  { return v_gsmap; }
  const gstat& gstatall() const // 全体統計量オブジェクトを返す
  { return v_gstatall; }
  const gstat& gstatcom() const // 共通統計量オブジェクトを返す
  { return v_gstatcom; }
  int dfb() const               // 回帰自由度を返す
  { return ni() - 2; }
  int dfa() const               // 群自由度を返す
  { return ng() - 1; }
  int dfd() const               // 非平行性自由度を返す
  { return dfb() * dfa(); }
  int dfr() const               // 残差自由度を返す
  { return dft() - dfa() - ng() * dfb(); }
  int dft() const               // 全体自由度を返す
  { return v_gstatall.n() - 1; }
  double sbc() const            // 共通回帰平方和を返す
  { return v_sbc; }
  double sa() const             // 群差平方和を返す
  { return v_sa; }
  double sbt() const            // 全体回帰平方和を返す
  { return v_sbt; }
  double saa() const            // 修正群差平方和を返す
  { return v_saa; }
  double sd() const             // 非平行性平方和を返す
  { return v_sd; }
  double sr() const             // 残差平方和を返す
  { return v_sr; }
  double st() const             // 全体平方和を返す
  { return v_gstatall.syy(); }
  void cal(option& opt, ldata& datalst); // データを集計して各種統計量を求める
};

void varclass::cal_ancova()    // 共分散分析用各種パラメターを計算する
{
  int g = v_gsmap.size();
  v_sa = 0.0;
  double ssb = 0.0;
  for (auto& gstati : v_gsmap)
  {
    v_sa += gstati.second.n() * gstati.second.my() * gstati.second.my();
    ssb += gstati.second.syy() - gstati.second.sr();
  }
  v_sa -= v_gstatall.n() * v_gstatall.my() * v_gstatall.my();
  if (v_sa < FDeps) v_sa = 0.0;
  v_sbt = st() - v_gstatall.sr();
  if (v_sbt < FDeps) v_sbt = 0.0;
  v_sbc = v_gstatcom.syy() - v_gstatcom.sr();
  if (v_sbc < FDeps) v_sbc = 0.0;
  v_saa = v_sa + v_sbc - v_sbt;
  if (v_saa < FDeps) v_saa = 0.0;
  v_sd = ssb - v_sbc;
  if (v_sd < FDeps) v_sd = 0.0;
  v_sr = st() - v_sa - ssb;
  if (dfr() < 1 || v_sr < FDeps) v_sr = 0.0;
}

void varclass::cal(option& opt, ldata& datalst)
// データを集計して各種統計量を求める
{
  int p = v_ni - 2;
  vecdbl fdatavec(p + 1);
  for (auto& cdatai : datalst)
  {
    bool blank = false;
    fdatavec[p] = cdatai.fdata(1);
    for (int j = 0; blank == false && j < p; ++j)
    {
      fdatavec[j] = cdatai.fdata(j + 2);
      if (fdatavec[j] == BLANKf) blank = true;
    }
    if (cdatai.sdata(0) != BLANKs && blank == false && fdatavec[p] != BLANKf)
    {
      mapsg_itr gitr = v_gsmap.find(cdatai.sdata(0));
      if (gitr == v_gsmap.end()) // 新規群コード
      {
        gstat gstati(v_ni);
        gstati.sigma(fdatavec);
        v_gsmap.insert(make_pair(cdatai.sdata(0), gstati));
      }
      else                      // すでにある群
      { gitr->second.sigma(fdatavec); }
    }
  }
  for (auto& gstati : v_gsmap)
  {
    v_gstatall += gstati.second;
    gstati.second.setsmat(true);
    v_gstatcom += gstati.second;
    gstati.second.calreg();
  }
  v_gstatall.setsmat(true);  v_gstatcom.setsmat(false);
  v_gstatall.calreg();  v_gstatcom.calreg();
  cal_ancova();
}

class ancovaclass               // ancovaメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout); // 結果のタイトルとデータ識別情報を出力する
  void put_mreg(ostream& fout, option& opt, const varclass& var);
  // 群別平均と回帰式を出力する
  void put_mregtab(ostream& fout, const gstat& gstati, const gstat& gstatcom) const;
  // 平均と回帰式の表を出力する
  void put_bg(ostream& fout, option& opt, const varclass& var) const;
  // 群別回帰係数の検定結果を出力する
  void put_bc(ostream& fout, option& opt, const varclass& var) const;
  // 共通回帰係数の検定結果を出力する
  void put_ancovatab(ostream& fout,double sbc, int dfb, double sa, int dfa, double sbt, double saa, double sd, int dfd, double sr, int dfr, double st, int dft) const;
  // 共分散分析表を出力する
  void put_macmp(ostream& fout, option& opt, double sr, int dfr, const varclass& var) const;
  // 修正群差の多重比較を出力する
  void put_maci(ostream& fout, option& opt, double sr, int dfr, const varclass& var) const;
  // 修正群差の信頼区間を出力する
  void put_mdcmp(ostream& fout, option& opt, double sr, int dfr, const varclass& var) const;
  // 非平行性の多重比較を出力する
  void put_dnsid(ofstream& fout);
  // DANSデータのタイトルとデータ識別情報を出力する
  void put_dns(option& opt, const varclass& var);
  // 群別回帰係数をデータファイル形式で出力する
public:
  ancovaclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~ancovaclass()                // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

void ancovaclass::put_resid(ostream& fout)
// 結果のタイトルとデータ識別情報を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
  fout << "群項目    :" << id.item(0) << "\n";
  fout << "集計項目y :" << id.item(1) << "\n";
  for (size_t i = 2; i < id.ni(); ++i)
  {
    fout << "共変数x" << setw(3) << i - 1;
    fout << ":" << id.item(i) << "\n";
  }
}

void ancovaclass::put_mreg(ostream& fout, option& opt, const varclass& var)
// 群別平均と回帰式を出力する
{
  if (opt.res("gstat"))
  {
    int gno = 1;
    for (auto& gstati : var.gsmap())
    {
      fout << "\n・群" << setw(5) << gno;
      fout << ":" << id.item(0) << "=" << gstati.first << "\n";
      put_mregtab(fout, gstati.second, var.gstatcom());
      ++gno;
    }
  }
  if (opt.res("tstat"))
  {
    fout << "\n・全体\n";
    put_mregtab(fout, var.gstatall(), var.gstatcom());
  }
}

void ancovaclass::put_mregtab(ostream& fout, const gstat& gstati, const gstat& gstatcom) const
// 平均と回帰式の表を出力する
{
  constexpr int LENline = 77;   // 区切り線の長さ
  constexpr int Nterm = 4;      // 1行あたりの項数
  putline(fout, LENline);
  int p = id.ni() - 2;
  double bc0 = gstati.my(), syy = gstati.syy();
  fout << "y        :";
  putmsd(fout, gstati.n(), bc0, gstati.vyy());
  for (int j = 1; j <= p; ++j)
  {
    fout << "x" << left << setw(8) << j << right << ":";
    putmsd(fout, gstati.n(), gstati.me(j), gstati.vr(j));
  }
  putline(fout, LENline);
  double sb = syy - gstati.sr();
  double r2 = (syy > 0.0) ? sb / syy : 0.0;
  fout << "群別回帰式:y=" << setw(13) << gstati.b(0);
  for (int j = 1; j <= p; ++j)
  {
    if ((j + 1) % Nterm == 1) fout << "\n" << setw(13) << "";
    fout << showpos << setw(10) << gstati.b(j) << noshowpos;
    fout << "x" << left << setw(2) << j << right;
    bc0 -= gstatcom.b(j) * gstati.me(j);
  }
  fout << "\n共通回帰式:y=" << setw(13) << bc0;
  for (int j = 1; j <= p; ++j)
  {
    if ((j + 1) % Nterm == 1) fout << "\n" << setw(13) << "";
    fout << showpos << setw(10) << gstatcom.b(j) << noshowpos;
    fout << "x" << left << setw(2) << j << right;
  }
  double r = sqrt(r2);
  if (p == 1) r = sgn(gstati.b(1)) * r;
  int dfr = gstati.n() - 1 - p;
  double f = (gstati.sr() > 0.0) ? sb / (double)p / gstati.sr() * (double)dfr : 0.0;
  ldbl pr = ftop(f, p, dfr);
  fout << "\n群別回帰式の寄与率:r2乗=" << r2 << "  r=" << r << "  有意確率:p=" << pr << pind(pr) << endl;
}

void ancovaclass::put_bg(ostream& fout, option& opt, const varclass& var) const
// 群別回帰係数の検定結果を出力する
{
  const int LENline = 78;       // 区切り線の長さ
  int p = id.ni() - 2, g = (var.gsmap()).size();
  double vr = (var.dfr() > 0) ? var.sr() / (double)(var.dfr()) : 0.0;
  int df1 = 1, dc = g * (g - 1) / 2;
  string method = "Bonferroni";
  switch (opt.md())
  {
  case mcfis : method = "Fisher";  dc = 1;
    break;
  case mcsch : method = "Scheffe";  df1 = g - 1;  dc = 1;
  }
  fout << "\n";
  fout << "・群別回帰係数の" << method << "型多重比較\n";
  fout << "群  -群    変数  回帰係数差" << setw(5) << opt.cc() * 100.0 << "%CI下限        上限         F値 有意確率p値\n";
  putline(fout, LENline);
  int gno1 = 1;
  for (mapsg_citr gitr1 = (var.gsmap()).begin(); gitr1 != (var.gsmap()).end(); ++gitr1, ++gno1)
  {
    mapsg_citr gitr2 = gitr1;
    ++gitr2;
    if (gitr2 != (var.gsmap()).end())
    {
      int gno2 = gno1 + 1;
      while (gitr2 != (var.gsmap()).end())
      {
        for (int i = 1; i <= p; ++i)
        {
          double bdi = (gitr2->second).b(i) - (gitr1->second).b(i);
          double vbdi = vr * ((gitr1->second).vbinv(i) + (gitr2->second).vbinv(i));
          double fi = (sgn(vbdi) > 0) ? bdi * bdi / vbdi / (double)df1 : 0.0;
          double pr = calmcpr(opt.md(), fi, g, df1, var.dfr(), dc);
          double ci = calmcci(opt.md(), opt.cc(), g, df1, var.dfr(), dc, vbdi);
          fout << left << setw(4) << gno2 << "-" << setw(4) << gno1;
          fout << "  x" << setw(3) << i << right << setw(12) << bdi;
          fout << setw(12) << bdi - ci << setw(12) << bdi + ci;
          fout << setw(12) << fi << setw(12) << pr << pind(pr) << "\n";
        }
        ++gitr2;  ++gno2;
      }
    }
  }
  putline(fout, LENline);
}

void ancovaclass::put_bc(ostream& fout, option& opt, const varclass& var) const
// 共通回帰係数の検定結果を出力する
{
  constexpr int LENline = 68;   // 区切り線の長さ
  int p = id.ni() - 2;
  double vr = (var.dfr() > 0) ? var.sr() / (double)(var.dfr()) : 0.0;
  fout << "\n";
  fout << "・共通回帰係数の区間推定と検定\n";
  fout << "変数     回帰係数" << setw(5) << opt.cc() * 100.0 << "%CI下限        上限         F値 有意確率p値\n";
  putline(fout, LENline);
  for (int i = 1; i <= p; ++i)
  {
    double bi = (var.gstatcom()).b(i), vbi = vr * (var.gstatcom()).vbinv(i);
    double fi = (sgn(vbi) > 0) ? bi * bi / vbi : 0.0;
    double pr = ftop(fi, 1, var.dfr());
    double ci = calmcci(mcfis, opt.cc(), 1, 1, var.dfr(), 1, vbi);
    fout << "x" << left << setw(4) << i << right << setw(12) << bi;
    fout << setw(12) << bi - ci << setw(12) << bi + ci;
    fout << setw(12) << fi << setw(12) << pr << pind(pr) << "\n";
  }
  putline(fout, LENline);
}

void ancovaclass::put_ancovatab(ostream& fout,double sbc, int dfb, double sa, int dfa, double sbt, double saa, double sd, int dfd, double sr, int dfr, double st, int dft) const
// 共分散分析表を出力する
{
  fout << "\n";
  constexpr int LENline = 78;   // 区切り線の長さ
  string tstr;
  if (dfd > 0) tstr = "共分散分析表(ANCOVA table)";
  else tstr = "共分散分析表(非平行性プール)";
  int w = (LENline - strclen(tstr, CHARwide)) / 2 + tstr.size();
  fout << setw(w) << tstr;
  fout << "\n";
  fout << "要因              平方和    自由度    平均平方和         F値   有意確率p値\n";
  putline(fout, LENline);
  double v = 0.0, f = 0.0, pr = 0.0;
  double vr = calvfp(v, f, pr, sa, dfa, sr, dfr);
  tstr = "群差";
  int npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << sa;
  fout << setw(10) << dfa << setw(14) << v << setw(12) << f;
  fout << setw(14) << pr << pind(pr) << "\n";
  calvfp(v, f, pr, sbc, dfb, sr, dfr);
  tstr = "共通回帰";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << sbc;
  fout << setw(10) << dfb << setw(14) << v << setw(12) << f;
  fout << setw(14) << pr << pind(pr) << "\n";
  putline(fout, LENline);
  calvfp(v, f, pr, saa, dfa, sr, dfr);
  tstr = "修正群差";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << saa;
  fout << setw(10) << dfa << setw(14) << v << setw(12) << f;
  fout << setw(14) << pr << pind(pr) << "\n";
  calvfp(v, f, pr, sbt, dfb, sr, dfr);
  tstr = "全体回帰";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << sbt;
  fout << setw(10) << dfb << setw(14) << v << setw(12) << f;
  fout << setw(14) << pr << pind(pr) << "\n";
  putline(fout, LENline);
  if (dfd > 0)
  {
    calvfp(v, f, pr, sd, dfd, sr, dfr);
    tstr = "非平行性";
    npad = tstr.size() - strclen(tstr, CHARwide);
    fout << left << setw(10 + npad) << tstr << right << setw(14) << sd;
    fout << setw(10) << dfd << setw(14) << v << setw(12) << f;
    fout << setw(14) << pr << pind(pr) << "\n";
  }
  tstr = "残差";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << sr;
  fout << setw(10) << dfr << setw(14) << vr << "\n";
  putline(fout, LENline);
  tstr = "全体";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << st;
  fout << setw(10) << dft << endl;
}

void ancovaclass::put_macmp(ostream& fout, option& opt, double sr, int dfr, const varclass& var) const
// 修正群差の多重比較を出力する
{
  constexpr int LENline = 77;   // 区切り線の長さ
  if (sr <= 0.0) return;        // 出力なし
  int g = (var.gsmap()).size(), p = id.ni() - 2;
  int df1 = 1, dc = 1;
  string method = "Tukey";
  switch (opt.ma())
  {
  case mcbon : method = "Bonferroni";  dc = g * (g - 1) / 2;
    break;
  case mcfis : method = "Fisher";
    break;
  case mcsch : method = "Scheffe";  df1 = g - 1;
  }
  fout << "\n・修正群差の" << method << "型多重比較:共変数=全体平均値  ";
  if (opt.ma() == mctuk)
  {
    fout << "群数=" << g << "  自由度=" << dfr << "\n";
    fout << "群 :  修正平均値 - 群 :  修正平均値    修正群差         q値   有意確率p値\n";
  }
  else
  {
    fout << "第1自由度=" << df1 << "  第2自由度=" << dfr << "\n";
    fout << "群 :  修正平均値 - 群 :  修正平均値    修正群差         F値   有意確率p値\n";
  }
  putline(fout, LENline);
  double vr = sr / dfr;
  vec<double> dxv(p);
  double yc = 0.0;
  for (int k = 1; k <= p; ++k) yc += (var.gstatcom()).b(k) * (var.gstatall()).me(k);
  mapsg_citr gitri = (var.gsmap()).begin();
  for (int i = 0; i < g - 1; ++i, ++gitri)
  {
    double ni = (double)(gitri->second.n());
    double bc0i = gitri->second.my();
    for (int k = 1; k <= p; ++k) bc0i -= (var.gstatcom()).b(k) * gitri->second.me(k);
    double yi = bc0i + yc;
    mapsg_citr gitrj = gitri;
    ++gitrj;
    int j = i + 1;
    while (gitrj != (var.gsmap()).end())
    {
      double bc0j = gitrj->second.my();
      for (int k = 1; k <= p; ++k)
      {
        bc0j -= (var.gstatcom()).b(k) * gitrj->second.me(k);
        dxv.assign(k - 1, gitri->second.me(k) - gitrj->second.me(k));
      }
      double yj = bc0j + yc;
      double dy = bc0j - bc0i;
      double vd = (1.0 / ni + 1.0 / (double)(gitrj->second.n()) + (var.gstatcom()).vxr(dxv)) * vr;
      double f = (sgn(vd) > 0) ? dy * dy / vd / (double)df1 : 0.0;
      double pr = calmcpr(opt.ma(), f, g, df1, dfr, dc);
      fout << left << setw(3) << j + 1 << ":" << right << setw(12) << yj << " - ";
      fout << left << setw(3) << i + 1 << ":" << right << setw(12) << yi << setw(12) << dy;
      fout << setw(12) << f << setw(14) << pr << pind(pr) << "\n";
      ++gitrj;  ++j;
    }
  }
  putline(fout, LENline);
}

void ancovaclass::put_maci(ostream& fout, option& opt, double sr, int dfr, const varclass& var) const
// 修正群差の信頼区間を出力する
{
  constexpr int LENline = 57;   // 区切り線の長さ
  if (sr <= 0.0) return;        // 出力なし
  int g = (var.gsmap()).size(), p = id.ni() - 2;
  int df1 = 1, dc = 1;
  string method = "Tukey";
  switch (opt.ma())
  {
  case mcbon : method = "Bonferroni";  dc = g * (g - 1) / 2;
    break;
  case mcfis : method = "Fisher";
    break;
  case mcsch : method = "Scheffe";  df1 = g - 1;
  }
  fout << "\n・修正群差の" << method << "型" << opt.cc() * 100.0 << "%同時信頼区間(simultaneous confidence interval)\n";
  fout << "群  - 群     修正群差    絶対精度        下限        上限\n";
  putline(fout, LENline);
  double vr = sr / dfr;
  vec<double> dxv(p);
  mapsg_citr gitri = (var.gsmap()).begin();
  for (int i = 0; i < g - 1; ++i, ++gitri)
  {
    double ni = (double)(gitri->second.n());
    double bc0i = gitri->second.my();
    for (int k = 1; k <= p; ++k) bc0i -= (var.gstatcom()).b(k) * gitri->second.me(k);
    mapsg_citr gitrj = gitri;
    ++gitrj;
    int j = i + 1;
    while (gitrj != (var.gsmap()).end())
    {
      double bc0j = gitrj->second.my();
      for (int k = 1; k <= p; ++k)
      {
        bc0j -= (var.gstatcom()).b(k) * gitrj->second.me(k);
        dxv.assign(k - 1,  gitri->second.me(k) - gitrj->second.me(k));
      }
      double dy = bc0j - bc0i;
      double vd = (1.0 / ni + 1.0 / (double)(gitrj->second.n()) + (var.gstatcom()).vxr(dxv)) * vr;
      double ci = calmcci(opt.ma(), opt.cc(), g, df1, dfr, dc, vd);
      fout << left << setw(4) << j + 1 << "- " << setw(3) << i + 1 << right;
      fout << setw(12) << dy << setw(12) << ci << setw(12) << dy - ci << setw(12) << dy + ci << "\n";
      ++gitrj;  ++j;
    }
  }
  putline(fout, LENline);
}

void ancovaclass::put_mdcmp(ostream& fout, option& opt, double sr, int dfr, const varclass& var) const
// 非平行性の多重比較を出力する
{
  constexpr int LENline = 57;   // 区切り線の長さ
  if (sr <= 0.0) return;        // 出力なし
  int g = (var.gsmap()).size(), p = id.ni() - 2;
  int df1 = p, dc = g * (g - 1) / 2;
  string method = "Bonferroni";
  switch (opt.md())
  {
  case mcfis : method = "Fisher";  dc = 1;
    break;
  case mcsch : method = "Scheffe";  df1 = (g - 1) * p;  dc = 1;
  }
  fout << "\n・非平行性の" << method << "型多重比較:第1自由度=" << df1 << "  第2自由度=" << var.dfr() << "\n";
  fout << "群  - 群   回帰係数差平方和         F値   有意確率p値\n";
  putline(fout, LENline);
  double vr = sr / dfr;
  vec<double> dbv(p);
  mapsg_citr gitri = (var.gsmap()).begin();
  for (int i = 0; i < g - 1; ++i, ++gitri)
  {
    mapsg_citr gitrj = gitri;
    ++gitrj;
    int j = i + 1;
    while (gitrj != (var.gsmap()).end())
    {
      gstat gstatc(p + 2);
      gstatc += gitri->second;  gstatc += gitrj->second;
      gstatc.calreg();
      for (int k = 1; k <= p; ++k)
      { dbv.assign(k - 1, gitri->second.b(k) - gitrj->second.b(k)); }
      double sdb = gstatc.vxr(dbv);
      double f = (sgn(vr) > 0) ? sdb / vr / df1 : 0.0;
      double pr = calmcpr(opt.md(), f, g, df1, dfr, dc);
      fout << left << setw(4) << j + 1 << "- " << setw(3) << i + 1 << right;
      fout << setw(18) << sdb << setw(12) << f;
      fout << setw(14) << pr << pind(pr) << "\n";
      ++gitrj;  ++j;
    }
  }
  putline(fout, LENline);
}

void ancovaclass::put_dnsid(ofstream& fout)
// DANSデータのタイトルとデータ識別情報を出力する
{
  idd tid(id.title(), id.ni());
  tid.assign_item(0, id.item(0));
  ostringstream ostr;
  ostr << "群別回帰式の切片";
  tid.assign_item(1, ostr.str());
  for (size_t i = 2; i < id.ni(); ++i)
  {
    ostr.str("");
    ostr << "群別回帰式の" << id.item(i) << "の回帰係数";
    tid.assign_item(i, ostr.str());
  }
  id.put_title(fout, formtag);
  id.put_stra(fout, formtag);
  tid.put_item(fout, formtag);
}

void ancovaclass::put_dns(option& opt, const varclass& var)
// 群別回帰係数をデータファイル形式で出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.dnsname()).c_str()) == false) return; // ファイルオープン失敗
  put_dnsid(fout);
  datac cdata(id.ni());
  for (auto& gstati : var.gsmap())
  {
    cdata.assign(0, gstati.first);
    cdata.assign(1, (gstati.second).b(0));
    for (size_t i = 2; i < id.ni(); ++i) cdata.assign(i, (gstati.second).b(i - 1));
    cdata.put(fout);
  }
}

int ancovaclass::get(const char* c_fname, option& opt)
// データを読み込み保存・累積する
{
  if (c_fname == nullptr) fname.erase();
  else fname = c_fname;
  int nc = fgetdata(fname, id, datalst);
  if (nc < 0) throw runtime_error("実行時エラー！");
  if (opt.check(id.ni()) != 0) throw option_error("実行時オプションエラー！");
  return nc;
}

void ancovaclass::put(option& opt) // 集計結果を出力する
{
  varclass var(id.ni());          // 変数オブジェクト
  var.cal(opt, datalst);
  if (var.nt() < 1)             // 有効なデータなし
  {
    if (opt.res("res"))
    {
      put_resid(cout);
      cout << "※有効なデータなし！\n\n" << endl;
    }
    return;
  }
  double vd = 0.0, fd = 0.0, pr = 0.0;
  calvfp(vd, fd, pr, var.sd(), var.dfd(), var.sr(), var.dfr());
  if (opt.res("res"))
  {
    put_resid(cout);
    put_mreg(cout, opt, var);
    if (opt.res("bt"))
    {
      put_bg(cout, opt, var);
      put_bc(cout, opt, var);
    }
    if (var.ng() < 2)
    { cout << "※共分散分析不可能！\n"; }
    else
    {
      if (opt.res("atab")) put_ancovatab(cout, var.sbc(), var.dfb(), var.sa(), var.dfa(), var.sbt(), var.saa(), var.sd(), var.dfd(), var.sr(), var.dfr(), var.st(), var.dft());
      if (opt.res("ma"))
      {
        put_macmp(cout, opt, var.sr(), var.dfr(), var);
        put_maci(cout, opt, var.sr(), var.dfr(), var);
      }
      if (opt.res("md")) put_mdcmp(cout, opt, var.sr(), var.dfr(), var);
      if (var.dfd() > 0)          // 非平行性をプールした結果出力
      {
        if ((opt.res("def") && fd < 1.0) || opt.res("pd"))
        {
          double srt = var.sd() + var.sr();
          int dfrt = var.dfd() + var.dfr();
          if (opt.res("atab")) put_ancovatab(cout, var.sbc(), var.dfb(), var.sa(), var.dfa(), var.sbt(), var.saa(), 0.0, 0, srt, dfrt, var.st(), var.dft());
          if (opt.res("ma"))
          {
            put_macmp(cout, opt, srt, dfrt, var);
            put_maci(cout, opt, srt, dfrt, var);
          }
        }
      }
    }
    cout << "\n" << endl;
  }
  if (opt.dnsname() != "") put_dns(opt, var);
}

int main(int argc, char *argv[]) // メインルーチン
{
  ancovaclass ancova;           // メインクラス変数
  option opt(argv[0]);          // オプション
  try
  {
    opt.setdef();               // オプションの初期設定
    bool fok = true;            // ファイル処理フラッグ
    for (int i = 1; i <= argc; ++i)
    {
      if (i >= argc)
      {
        if (fok) ancova.get(nullptr, opt); // 標準入力から読み込む
      }
      else if (*(argv[i]) == CHRopt)
      { opt.set(argv[i] + 1);  fok = true; }
      else
      { ancova.get(argv[i], opt);  fok = false; }
    }
    ancova.put(opt);
  }
  catch (help_except& helpex)   // ヘルプファイル表示
  { return 0; }
  catch (option_error& opterr)  // オプション指定エラー
  { return 1; }
  catch (runtime_error& rterr)  // 実行時エラー
  {
    disperr(opt.pname(), rterr.what());
    return 2;
  }
  catch (...)                   // 予期せぬ例外
  {
    disperr(opt.pname(), "内容不明なエラー発生！");
    return 255;
  }
  return 0;
}
