// anova2.C Version 7.44  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-28 10:17:50 tac>
//
// 対応のある多時期間で二元配置分散分析を行う。
// 繰り返し測定型二元配置分散分析も含む。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -w anova2 anova2.C
//             -I/usr/local/include -L/usr/local/lib -ldans -lfun

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

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

const vecstr USAGE =            // プログラム使用法表示用
{
  "anova2  Version 7.44         Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:anova2 [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目を群項目(要因A)、以後の項目を時期項目(要因B)とした",
  "     二元配置分散分析をするんや。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定するんや。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -f=<ファイル名>:指定したファイルからオプションを読み込むんや。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示するんや。(-h:-h=man)",
  "  -res=[-][{def|all|a2|icc|m[{a|b|c|d|e|f|m|x}][{b|d|f|s|t}]|rm",
  "       |rp[<比較回数>]|sd|stat}][,…]:解析結果の出力項目を指定するんや。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    a2:一般型二元配置分散分析  icc:級内相関係数",
  "    m[{a|b|c|d|e|f|m|x}][{b|d|f|s|t}]:多重比較(ma:全群の多重比較",
  "    mb:時期1と他時期の多重比較  mc:全時期の多重比較  md:群1と他群の多重比較",
  "    me:群1と他群の時期別多重比較  mf:群1と他群の交互作用の多重比較",
  "    mm:全群の時期別多重比較  mx:全群の交互作用の多重比較)",
  "      b:Bonferroni型  d:Dunnett型  f:Fisher型  s:Scheffe型  t:Tukey型",
  "    rm:繰り返し測定型二元配置分散分析",
  "    rp[<比較回数>]:対応のある二元配置分散分析  sd:各要因のSD  stat:基礎統計量",
  "    オプション省略時:データに合わせて自動判定"
};

const string Title = "=== 多群・多時期の平均値の比較 ==="; // 表題

using vecstatv = vector<statv<double>>; // 統計量ベクトル配列の別名定義

enum atabtype                   // 分散分析表の種類
{
  atabnon,                      // 出力なし
  at2,                          // 繰り返しのない二元配置分散分析表
  at2r,                         // 繰り返しのある二元配置分散分析表
  at2rm,                        // 繰り返し測定型二元配置分散分析表
  at2rp,                        // 対応のある二元配置分散分析表
  atabetc                       // その他
};

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  mctype v_ma;                  // 群の多重比較法(全群比較)
  mctype v_mb;                  // 時期の多重比較法(時期1との比較)
  mctype v_mc;                  // 時期の多重比較法(全時期比較)
  mctype v_md;                  // 群の多重比較法(群1との比較)
  mctype v_me;                  // 時期別平均値の多重比較法(群1との比較)
  mctype v_mf;                  // 交互作用の多重比較法(群1との比較)
  mctype v_mm;                  // 時期別平均値の多重比法較(全群比較)
  mctype v_mx;                  // 交互作用の多重比較法(全群比較)
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
  int v_rp;                     // 対応のある二元配置分散分析の比較回数
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; }
  mctype ma() const             // 群の多重比較法を返す(全群比較)
  { return v_ma; }
  mctype mb() const           // 時期の多重比較法を返す(時期1との比較)
  { return v_mb; }
  mctype mc() const             // 時期の多重比較法を返す(全時期比較)
  { return v_mc; }
  mctype md() const             // 群の多重比較法を返す(群1との比較)
  { return v_md; }
  mctype me() const     // 時期別平均値の多重比較法を返す(群1との比較)
  { return v_me; }
  mctype mf() const         // 交互作用の多重比較法を返す(群1との比較)
  { return v_mf; }
  mctype mm() const        // 時期別平均値の多重比較法を返す(全群比較)
  { return v_mm; }
  mctype mx() const            // 交互作用の多重比較法を返す(全群比較)
  { return v_mx; }
  bool res(const string& code); // 分析結果出力フラッグを返す
  int rp() const         // 対応のある二元配置分散分析の比較回数を返す
  { return v_rp; }
  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_ma(mctuk), v_mb(mcdun), v_mc(mctuk), v_md(mcdun), v_me(mcbon), v_mf(mcbon), v_mm(mcbon), v_mx(mcbon), reslst(), resmap(), v_rp(1)
// コンストラクタ
{
  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 respos = resmap.find(code);
  if (respos == 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] == "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 < 3) return disperr(pname(), "項目が3つ以上必要なんや！");
  if (set_res() < 0) return -3; // -resオプションのセット
  return 0;                     // エラーなし
}

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["stat"] = true;
  resmap["a2"] = resmap["icc"] = resmap["ma"] = resmap["mb"] = resmap["mc"] = resmap["md"] = resmap["me"] = resmap["mf"] = resmap["mm"] = resmap["mx"] = resmap["rm"] = resmap["rp"] = resmap["sd"] = false;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["stat"] = false;
    for (const string& resstr : reslst)
    {
      mapsb_citr respos = resmap.find(resstr);
      if (resstr == "-")        // 全て出力せず
      { resmap["res"] = false; }
      else if (resstr[0] == '-' && resstr.size() > 1) // 特定の項目を出力せず
      {
        respos = resmap.find(resstr.substr(1));
        if (respos == resmap.end())
        { return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！"); }
        else
        { resmap[resstr.substr(1)] = false; }
      }
      else if (resstr == "def") // 省略時出力
      { resmap["def"] = resmap["stat"] = true; }
      else if (resstr == "all") // 全て出力
      {
        resmap["a2"] = resmap["icc"] = resmap["ma"] = resmap["mb"] = resmap["mc"] = resmap["md"] = resmap["me"] = resmap["mf"] = resmap["mm"] = resmap["mx"] = resmap["rm"] = resmap["sd"] = resmap["stat"] = true;
        // rp以外は全てtrue
      }
      else if (resstr == "m")   // mと指定された時
      { resmap["ma"] = resmap["mx"] = 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] == 'b' && resstr.size() > 2)
      {
        switch (resstr[2])
        {
        case 'b' : v_mb = mcbon;
          break;
        case 'd' : v_mb = mcdun;
          break;
        case 'f' : v_mb = mcfis;
          break;
        case 's' : v_mb = mcsch;
          break;
        default : return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！");
        }
        resmap["mb"] = true;
      }
      else if (resstr[0] == 'm' && resstr[1] == 'c' && resstr.size() > 2)
      {
        switch (resstr[2])
        {
        case 'b' : v_mc = mcbon;
          break;
        case 'f' : v_mc = mcfis;
          break;
        case 's' : v_mc = mcsch;
          break;
        case 't' : v_mc = mctuk;
          break;
        default : return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！");
        }
        resmap["mc"] = true;
      }
      else if (resstr[0] == 'm' && resstr[1] == 'd' && resstr.size() > 2)
      {
        switch (resstr[2])
        {
        case 'b' : v_md = mcbon;
          break;
        case 'd' : v_md = mcdun;
          break;
        case 'f' : v_md = mcfis;
          break;
        case 's' : v_md = mcsch;
          break;
        default : return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！");
        }
        resmap["md"] = true;
      }
      else if (resstr[0] == 'm' && resstr[1] == 'e' && resstr.size() > 2)
      {
        switch (resstr[2])
        {
        case 'b' : v_me = mcbon;
          break;
        case 'f' : v_me = mcfis;
          break;
        case 's' : v_me = mcsch;
          break;
        default : return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！");
        }
        resmap["me"] = true;
      }
      else if (resstr[0] == 'm' && resstr[1] == 'f' && resstr.size() > 2)
      {
        switch (resstr[2])
        {
        case 'b' : v_mf = mcbon;
          break;
        case 'f' : v_mf = mcfis;
          break;
        case 's' : v_mf = mcsch;
          break;
        default : return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！");
        }
        resmap["mf"] = true;
      }
      else if (resstr[0] == 'm' && resstr[1] == 'm' && resstr.size() > 2)
      {
        switch (resstr[2])
        {
        case 'b' : v_mm = mcbon;
          break;
        case 'f' : v_mm = mcfis;
          break;
        case 's' : v_mm = mcsch;
          break;
        default : return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！");
        }
        resmap["mm"] = true;
      }
      else if (resstr[0] == 'm' && resstr[1] == 'x' && resstr.size() > 2)
      {
        switch (resstr[2])
        {
        case 'b' : v_mx = mcbon;
          break;
        case 'f' : v_mx = mcfis;
          break;
        case 's' : v_mx = mcsch;
          break;
        default : return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！");
        }
        resmap["mx"] = true;
      }
      else if (resstr[0] == 'r' && resstr[1] == 'p' && resstr.size() > 2)
      {
        v_rp = int(strtodf(resstr.substr(2)));
        if (v_rp < 1) return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！");
        resmap["rp"] = true;
      }
      else if (respos == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！"); }
      else
      { resmap[resstr] = true; }
    }
    if ((resmap["def"] || resmap["a2"]) && resmap["rp"]) return disperr(pname(), "-res=a2とrpオプションは同時に指定できない！");
  }
  return 0;                     // エラーなし
}

class varclass                  // 変数クラス
{
private:
  map<string, vecstatv> gsmap;  // 群統計量マップ
  vector<string> gcodevec;      // 群コードベクトル
  vecstatv tstatvec;            // 群統計量ベクトル
  double v_sa;                  // 要因A平方和
  int v_dfa;                    // 要因A自由度
  double v_ssr;                 // 個体残差
  int v_dfsr;                   // 個体残差自由度
  double v_ss;                  // 個体間平方和
  int v_dfs;                    // 個体自由度
  double v_sb;                  // 要因B平方和
  int v_dfb;                    // 要因B自由度
  double v_sab;                 // 要因AxB平方和
  int v_dfab;                   // 要因AxB自由度
  double v_sr;                  // 残差平方和
  int v_dfr;                    // 残差自由度
  double v_st;                  // 全体平方和
  int v_dft;                    // 全体自由度
  double statv_sigma(vecstatv& statvec, const vecdbl& fdatavec);
  // データを累積する
  int set_gsmap(option& opt);  // 群統計量マップを整理して群数を返す
public:
  varclass(size_t ni) : gsmap(), gcodevec(), tstatvec(ni), v_sa(0.0), v_dfa(0), v_ssr(0.0), v_dfsr(0), v_ss(0.0), v_dfs(0), v_sb(0.0), v_dfb(0), v_sab(0.0), v_dfab(0), v_sr(0.0), v_dfr(0), v_st(0.0), v_dft(0)
  // デフォールトコンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  int g() const                 // 要因A(群)の数を返す
  { return (int)(gsmap.size()); }
  int b() const                 // 要因B(時期)の数を返す
  { return (int)(tstatvec.size()) - 1; }
  int ns() const                // 個体数を返す
  { return tstatvec[0].n(); }
  const string& gcode(size_t idx) const // 群コードを返す
  { return (idx < gcodevec.size()) ? gcodevec[idx] : BLANKs; }
  int n(size_t iidx, size_t jidx); // 群の例数を返す
  double mean(size_t iidx, size_t jidx); // 群の平均値を返す
  double vr(size_t iidx, size_t jidx);   // 群の不偏分散を返す
  double ss(size_t iidx, size_t jidx);   // 群の平方和を返す
  double se(size_t iidx, size_t jidx);   // 群の標準誤差を返す
  double sa() const             // 要因A平方和を返す
  { return v_sa; }
  int dfa() const               // 要因A自由度を返す
  { return v_dfa; }
  double ssr() const            // 個体残差平方和を返す
  { return v_ssr; }
  double ssr(double ssr0)       // 個体残差平方和を代入する
  { return (v_ssr = ssr0); }
  int dfsr() const              // 個体残差自由度を返す
  { return v_dfsr; }
  int dfsr(int dfsr0)           // 個体残差由度を代入する
  { return (v_dfsr = dfsr0); }
  double ss() const             // 個体平方和を返す
  { return v_ss; }
  int dfs() const               // 個体自由度を返す
  { return v_dfs; }
  double sb() const             // 要因B平方和を返す
  { return v_sb; }
  int dfb() const               // 要因B由度を返す
  { return v_dfb; }
  double sab() const            // 要因AxB平方和を返す
  { return v_sab; }
  int dfab() const              // 要因AxB由度を返す
  { return v_dfab; }
  double sr() const             // 残差平方和を返す
  { return v_sr; }
  double sr(double sr0)         // 残差平方和を代入する
  { return (v_sr = sr0); }
  int dfr() const               // 残差自由度を返す
  { return v_dfr; }
  int dfr(int dfr0)             // 残差自由度を代入する
  { return (v_dfr = dfr0); }
  double st() const             // 全体平方和を返す
  { return v_st; }
  int dft() const               // 全体自由度を返す
  { return v_dft; }
  void cal(option& opt, ldata& datalst); // データを集計して変数をセットする
};

double varclass::statv_sigma(vecstatv& statvec, const vecdbl& fdatavec)
// データを累積する
{
  double sum = 0.0;
  for (size_t j = 0; j < fdatavec.size(); ++j)
  {
    statvec[j].sigma(fdatavec[j]);
    sum += fdatavec[j];
  }
  return sum * sum / fdatavec.size(); // 個体平方和用の値を返す
}

int varclass::set_gsmap(option& opt)  // 群統計量マップを整理して群数を返す
{
  int g = (int)gsmap.size(), b = (int)tstatvec.size() - 1;
  gcodevec.resize(g);
  v_sa = 0.0;  v_sab = 0.0;
  int i = 0;
  for (auto& gstati : gsmap)
  {
    gcodevec[i] = gstati.first; // 群コードをべクルトに格納する
    for (int j = 0; j < b; ++j)
    {
      gstati.second[b] += gstati.second[j];
      tstatvec[j] += gstati.second[j];
      double m = gstati.second[j].mean();
      v_sab += (double)(gstati.second[j].n()) * m * m;
    }
    double m = gstati.second[b].mean();
    v_sa += (double)(gstati.second[b].n()) * m * m;
    ++i;
  }
  for (int j = 0; j < b; ++j)
  {
    tstatvec[b] += tstatvec[j];
    double m = tstatvec[j].mean();
    v_sb += (double)(tstatvec[j].n()) * m * m;
  }
  double m = tstatvec[b].mean();
  double cf = (double)(tstatvec[b].n()) * m * m;
  int ns = tstatvec[0].n();
  v_dfs = ns - 1;  v_dfb = b - 1;  v_ss -= cf;
  if (ns == g && (opt.res("def") == false && opt.res("a2") == false && opt.res("rp")))
  {                   // 対応のある二元配置分散分析の場合(-res=a2優先)
    v_st = tstatvec[b].sum2();  v_dft = tstatvec[b].n();
    v_sab = v_sb - cf;  v_dfab = b - 1;
    v_ssr = v_sa - cf;  v_dfsr = g - 1;
    v_sa = cf;  v_dfa = 1;  v_sb = v_sab;
    v_sr = v_st - v_sa - v_ssr - v_sab;
    v_dfr = v_dft - v_dfa - v_dfsr - v_dfab;
  }
  else
  {
    v_st = tstatvec[b].ss();  v_dft = tstatvec[b].n() - 1;
    v_sa -= cf;  v_dfa = g - 1;  v_sb -= cf;
    v_ssr = v_ss - v_sa;  v_dfsr = v_dfs - v_dfa;
    v_sab -= cf + v_sa + v_sb;  v_dfab = v_dfa * v_dfb;
    v_sr = v_st - v_ss - v_sb - v_sab;
    v_dfr = v_dft - v_dfs - v_dfb - v_dfab;
  }
  if (v_dfa < 1 || v_sa < FDeps) v_sa = 0.0;
  if (v_dfb < 1 || v_sb < FDeps) v_sb = 0.0;
  if (v_dfs < 1 || v_ss < FDeps) v_ss = 0.0;
  if (v_dfsr < 1 || v_ssr < FDeps) v_ssr = 0.0;
  if (v_dfr < 1 || v_sr < FDeps) v_sr = 0.0;
  return g;                     // 群数を返す
}

void varclass::cal(option& opt, ldata& datalst)
// データを集計して変数をセットする
{
  size_t ni = tstatvec.size();
  vecdbl fdatavec(ni - 1);
  v_ss = 0.0;
  for (auto& cdatai : datalst)
  {
    bool blank = (cdatai.sdata(0) == BLANKs) ? true : false;
    size_t j = 1;
    while (blank == false && j < ni)
    {
      if (cdatai.fdata(j) == BLANKf)
      { blank = true; }
      else
      {
        fdatavec[j - 1] = cdatai.fdata(j);
        ++j;
      }
    }
    if (blank == false)
    {
      map<string, vecstatv >::iterator gitr = gsmap.find(cdatai.sdata(0));
      if (gitr == gsmap.end())  // 新規群コード
      {
        vecstatv tstatvec(ni);
        v_ss += statv_sigma(tstatvec, fdatavec);
        gsmap.insert(map<string, vecstatv>::value_type(cdatai.sdata(0), tstatvec));
      }
      else                      // すでにある群
      { v_ss += statv_sigma(gitr->second, fdatavec); }
    }
  }
  set_gsmap(opt);               // 群統計量マップを整理する
}

int varclass::n(size_t iidx, size_t jidx) // 群の例数を返す
{
  int n = 0, b = (int)tstatvec.size() - 1;
  if (iidx < gsmap.size())     // 群別
  {
    if (jidx < tstatvec.size()) n = (gsmap[gcodevec[iidx]])[jidx].n();
    else n = (gsmap[gcodevec[iidx]])[b].n();
  }
  else                          // 群合計
  {
    if (jidx < tstatvec.size()) n = tstatvec[jidx].n();
    else n = tstatvec[b].n();
  }
  return n;
}

double varclass::mean(size_t iidx, size_t jidx) // 群の平均値を返す
{
  int b = (int)tstatvec.size() - 1;
  double m = 0.0;
  if (iidx < gsmap.size())      // 群別
  {
    if (jidx < tstatvec.size()) m = (gsmap[gcodevec[iidx]])[jidx].mean();
    else m = (gsmap[gcodevec[iidx]])[b].mean();
  }
  else                          // 群合計
  {
    if (jidx < tstatvec.size()) m = tstatvec[jidx].mean();
    else m = tstatvec[b].mean();
  }
  return m;
}

double varclass::vr(size_t iidx, size_t jidx) // 群の不偏分散を返す
{
  int b = (int)tstatvec.size() - 1;
  double vr = 0.0;
  if (iidx < gsmap.size())      // 群別
  {
    if (jidx < tstatvec.size()) vr = (gsmap[gcodevec[iidx]])[jidx].vr();
    else vr = (gsmap[gcodevec[iidx]])[b].vr();
  }
  else                          // 群合計
  {
    if (jidx < tstatvec.size()) vr = tstatvec[jidx].vr();
    else vr = tstatvec[b].vr();
  }
  return vr;
}

double varclass::ss(size_t iidx, size_t jidx) // 群の平方和を返す
{
  int b = (int)tstatvec.size() - 1;
  double ss = 0.0;
  if (iidx < gsmap.size())      // 群別
  {
    if (jidx < tstatvec.size()) ss = (gsmap[gcodevec[iidx]])[jidx].ss();
    else ss = (gsmap[gcodevec[iidx]])[b].ss();
  }
  else                          // 群合計
  {
    if (jidx < tstatvec.size()) ss = tstatvec[jidx].ss();
    else ss = tstatvec[b].ss();
  }
  return ss;
}

double varclass::se(size_t iidx, size_t jidx) // 群の標準誤差を返す
{
  int b = (int)tstatvec.size() - 1;
  double se = 0.0;
  if (iidx < gsmap.size())      // 群別
  {
    if (jidx < tstatvec.size()) se = (gsmap[gcodevec[iidx]])[jidx].se();
    else se = (gsmap[gcodevec[iidx]])[b].se();
  }
  else                          // 群合計
  {
    if (jidx < tstatvec.size()) se = tstatvec[jidx].se();
    else se = tstatvec[b].se();
  }
  return se;
}

class anova2class               // anova2メインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_id(ostream& fout);   // タイトルとデータ識別情報を出力する
  void put_stat(ostream& fout, option& opt, varclass& var);
  // 群別基礎統計量を出力する
  void put_icc(ostream& fout, option& opt, varclass& var);
  // 級内相関係数(信頼性係数)を出力する
  void put_iccci(ostream& fout, const string& kicc, double cc, double v1, int df1, double v2, int df2, double b) const;
  // 級内相関係数の信頼区間を出力する
  void put_icc2ci(ostream& fout, double cc, double icc21, double va, double vb, double vr, int dfa, int dfb) const;
  // ICC(2,1)とICC(2,b)の信頼区間を出力する
  void put_rmanova2(ostream& fout, option& opt, varclass& var);
  // 繰り返し測定型二元配置分散分析結果を出力する
  void put_ranova2(ostream& fout, option& opt, varclass& var);
  // 繰り返しのある二元配置分散分析結果を出力する
  void put_nranova2(ostream& fout, option& opt, varclass& var);
  // 繰り返しのない二元配置分散分析結果を出力する
  void put_panova2(ostream& fout, option& opt, varclass& var);
  // 対応のある二元配置分散分析結果を出力する
  void put_anova2(ostream& fout, atabtype attype, varclass& var) const;
  // 二元配置分散分析表を出力する
  void put_sd(ostream& fout, atabtype attype, varclass& var) const;
  // 各要因のSDを出力する
  void put_mcmpad(ostream& fout, mctype mc, bool ma, double sr, int dfr, varclass& var) const;
  // 群の多重比較を出力する(全群比較／群1との比較)
  void put_mdciad(ostream& fout, option& opt, mctype mc, bool ma, double sr, int dfr, varclass& var) const;
  // 2群間の平均値の差の信頼区間を出力する(全群比較／群1との比較)
  void put_mcmpbc(ostream& fout, mctype mc, bool mb, double sr, int dfr, varclass& var) const;
  // 時期の多重比較を出力する(時期1との比較／全時期比較)
  void put_mdcibc(ostream& fout, option& opt, mctype mc, bool mb, double sr, int dfr, varclass& var) const;
  // 2時期間の平均値の差の信頼区間を出力する(時期1との比較／全時期比較)
  void put_mcmpxf(ostream& fout, mctype mc, bool mcx, double sr, int dfr, varclass& var) const;
  // 交互作用の多重比較を出力する(全群比較／群1との比較)
  void put_mcmpme(ostream& fout, mctype mc, bool mm, varclass& var) const;
  // 時期別群別平均値の多重検定を出力する(全群比較／群1との比較)
  void put_mdcime(ostream& fout, option& opt, mctype mc, bool mm, varclass& var) const;
  // 時期別群別平均値の差の信頼区間を出力する(全群比較／群1との比較)
  void put_ttest1(ostream& fout, option& opt, double ssr, int dfsr, varclass& var);
  // 総平均値の1標本t検定(対応のある二元配置の時)
  void put_mcmppm(ostream& fout, option& opt, varclass& var) const;
  // 時期別平均値の多重検定を出力する(対応のある二元配置の時、全群比較)
  void put_mdcipm(ostream& fout, option& opt, varclass& var) const;
  // 時期別平均値の信頼区間を出力する(対応のある二元配置の時、時期1との比較)
public:
  anova2class() : fname(), id(), datalst() // コンストラクタ
  {}
  ~anova2class()                // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み保存・累積する
  void put(option& opt);                     // 集計結果を出力する
};

void anova2class::put_id(ostream& fout) // タイトルとデータ識別情報を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
  fout << "要因A  :" << id.item(0) << "\n";
  for (size_t i = 1; i < id.ni(); ++i)
  {
    fout << "要因B" << setw(2) << i;
    fout << ":" << id.item(i) << "\n";
  }
}

void anova2class::put_stat(ostream& fout, option& opt, varclass& var)
// 群別基礎統計量を出力する
{
  constexpr int LENline = 77;   // 区切り線の長さ
  constexpr int LENhed = 9;     // ヘッダの長さ
  int g = var.g(), b = var.b(), ns = var.ns();
  if (ns > g || opt.res("a2"))  // 繰り返しあり/一般型二元配置分散分析
  {
    if (ns == g)                // 繰り返しのない二元配置分散分析
    { fout << "個体     :個体(要因A)別・全時期(要因B)合計基礎統計量\n"; }
    else
    { fout << "群       :群(要因A)別・全時期(要因B)合計基礎統計量\n"; }
    putline(fout, LENline);
    for (int i = 0; i < g; ++i)
    {
      string gcodei(var.gcode(i));
      int npad = gcodei.size() - strclen(gcodei, CHARwide);
      fout << left << setw(LENhed + npad) << gcodei << right << ":";
      putmsd(fout, var.n(i, b), var.mean(i, b), var.vr(i, b));
    }
    putline(fout, LENline);
    fout << "\n";
  }
  if (ns == g)                  // 繰り返しのない二元配置分散分析
  { fout << "時期     :時期(要因B)別・全個体(要因A)合計基礎統計量\n"; }
  else
  { fout << "時期     :時期(要因B)別・全群(要因A)合計基礎統計量\n"; }
  putline(fout, LENline);
  for (int j = 0; j < b; ++j)
  {
    fout << left << setw(9) << j + 1 << right << ":";
    putmsd(fout, var.n(g, j), var.mean(g, j), var.vr(g, j));
  }
  putline(fout, LENline);
  if (ns > g)                   // 繰り返しあり
  {
    fout << "\n群   時期:群(要因A)別・時期(要因B)別基礎統計量\n";
    putline(fout, LENline);
    for (int i = 0; i < g; ++i)
    {
      for (int j = 0; j < b; ++j)
      {
        string gcodei(var.gcode(i));
        int npad = gcodei.size() - strclen(gcodei, CHARwide);
        fout << left << setw(LENhed - 4 + npad) << gcodei << right;
        fout << setw(4) << j + 1 << ":";
        putmsd(fout, var.n(i, j), var.mean(i, j), var.vr(i, j));
      }
    }
    putline(fout, LENline);
  }
  fout << "全体     :";
  putmsd(fout, var.n(g, b), var.mean(g, b), var.vr(g, b));
}

void anova2class::put_icc(ostream& fout, option& opt, varclass& var)
// 級内相関係数(信頼性係数)を出力する
{
  double a = (double)var.g(), b = (double)var.b(), sw = var.st() - var.sa();
  int dfw = var.dft() - var.dfa();
  double sr = var.sab();
  int dfr = var.dfab();
  double va = (var.dfa() > 0) ? var.sa() / (double)var.dfa() : 0.0;
  double vb = (var.dfb() > 0) ? var.sb() / (double)var.dfb() : 0.0;
  double vw = (dfw > 0) ? sw / (double)dfw : 0.0;
  double vr = (dfr > 0) ? sr / (double)dfr : 0.0;
  double vtmp = va + (double)var.dfb() * vw;
  double icc11 = (vtmp > 0.0) ? (va - vw) / vtmp : 0.0;
  double icc1b = (va > 0.0) ? (va - vw) / va : 0.0;
  vtmp = va + (double)var.dfb() * vr + b * (vb - vr) / a;
  double icc21 = (vtmp > 0.0) ? (va - vr) / vtmp : 0.0;
  vtmp = va + (vb - vr) / a;
  double icc2b = (vtmp > 0.0) ? (va - vr) / vtmp : 0.0;
  vtmp = va + (double)var.dfb() * vr;
  double icc31 = (vtmp > 0.0) ? (va - vr) / vtmp : 0.0;
  double icc3b = (va > 0.0) ? (va - vr) / va : 0.0;
  fout << "・級内相関係数(intraclass correlation coefficient)\n";
  double sda = (va - vw) / b, sdr = sqrt(vw);
  sda = (sda > 0.0) ? sqrt(sda) : 0.0;
  fout << "  ICC(1,1)=" << icc11 << " ICC(1," << b << ")=" << icc1b;
  fout << "  σA=" << sda << " σR=" << sdr << "\n";
  put_iccci(fout, "ICC(1", opt.cc(), va, var.dfa(), vw, dfw, b);
  sda = (va - vr) / b;  sdr = sqrt(vr);
  sda = (sda > 0.0) ? sqrt(sda) : 0.0;
  double sdb = (vb - vr) / a;
  sdb = (sdb > 0.0) ? sqrt(sdb) : 0.0;
  fout << "  ICC(2,1)=" << icc21 << " ICC(2," << b << ")=" << icc2b;
  fout << "  σA=" << sda << " σB=" << sdb << " σR=" << sdr << "\n";
  put_icc2ci(fout, opt.cc(), icc21, va, vb, vr, var.dfa(), var.dfb());
  fout << "  ICC(3,1)=r11=" << icc31 << " ICC(3," << b << ")=α係数=" << icc3b;
  fout << "  σA=" << sda << " σR=" << sdr << "\n";
  put_iccci(fout, "ICC(3", opt.cc(), va, var.dfa(), vr, dfr, b);
}

void anova2class::put_iccci(ostream& fout, const string& kicc, double cc, double v1, int df1, double v2, int df2, double b) const
// 級内相関係数の信頼区間を出力する
{
  double f = (v2 > 0.0) ? v1 / v2 : 0.0;
  double ah = (1.0 - cc) * 0.5;
  double fl = f / ptof(ah, df1, df2), fu = f / ptof(1.0 - ah, df1, df2);
  double iccl = (fl - 1.0) / (fl + (b - 1.0));
  double iccu = (fu - 1.0) / (fu + (b - 1.0));
  fout << "  " << kicc << ",1)の" << cc * 100.0 << "%信頼区間：下限=" << iccl << " 上限=" << iccu << "\n";
  iccl = (fl > 0.0) ? (fl - 1.0) / fl : 0.0;
  iccu = (fu > 0.0) ? (fu - 1.0) / fu : 0.0;
  fout << "  " << kicc << "," << b << ")の" << cc * 100.0 << "%信頼区間：下限=" << iccl << " 上限=" << iccu << "\n";
}

void anova2class::put_icc2ci(ostream& fout, double cc, double icc21, double va, double vb, double vr, int dfa, int dfb) const
// ICC(2,1)とICC(2,b)の信頼区間を出力する
{
  double f = (vr > 0.0) ? vb / vr : 0.0;
  double a = (double)(dfa + 1), b = (double)(dfb + 1);
  double bicc = b * icc21;
  double t1 = a * (1.0 + (double)dfb * icc21) - bicc;
  double t2 = (double)dfa * (double)dfb * (bicc * f + t1) * (bicc * f + t1);
  int df2 = int(t2 / ((double)dfa * bicc * bicc * f * f + t1 * t1) + 0.5);
  double ah = (1.0 - cc) * 0.5;
  double fl = ptof(ah, dfa, df2), fu = ptof(ah, df2, dfa);
  double iccl = (a * (va - fl * vr)) / (fl * (b * vb + (a * b - a - b) * vr) + a * va);
  double iccu = (a * (fu * va - vr)) / (b * vb + (a * b - a - b) * vr + a * fu * va);
  fout << "  ICC(2,1)の" << cc * 100.0 << "%信頼区間：下限=" << iccl << " 上限=" << iccu << "\n";
  double iccbl = b * iccl / (1.0 + (double)dfb * iccl);
  double iccbu = b * iccu / (1.0 + (double)dfb * iccu);
  fout << "  ICC(2," << b << ")の" << cc * 100.0 << "%信頼区間：下限=" << iccbl << " 上限=" << iccbu << "\n";
}

void anova2class::put_rmanova2(ostream& fout, option& opt, varclass& var)
// 繰り返し測定型二元配置分散分析結果を出力する
{
  put_anova2(fout, at2rm, var);
  if (opt.res("sd"))
  {
    fout << "\n";
    put_sd(fout, at2rm, var);
  }
  if (opt.res("ma"))            // 全群の多重比較
  {
    fout << "\n";
    put_mcmpad(fout, opt.ma(), true, var.ssr(), var.dfsr(), var);
    fout << "\n";
    put_mdciad(fout, opt, opt.ma(), true, var.ssr(), var.dfsr(), var);
  }
  if (opt.res("md"))            // 群1と他群の多重比較
  {
    fout << "\n";
    put_mcmpad(fout, opt.md(), false, var.ssr(), var.dfsr(), var);
    fout << "\n";
    put_mdciad(fout, opt, opt.md(), false, var.ssr(), var.dfsr(), var);
  }
  if (opt.res("mb"))            // 時期1と他時期の多重比較
  {
    fout << "\n";
    put_mcmpbc(fout, opt.mb(), true, var.sr(), var.dfr(), var);
    fout << "\n";
    put_mdcibc(fout, opt, opt.mb(), true, var.sr(), var.dfr(), var);
  }
  if (opt.res("mc"))            // 全時期の多重比較
  {
    fout << "\n";
    put_mcmpbc(fout, opt.mc(), false, var.sr(), var.dfr(), var);
    fout << "\n";
    put_mdcibc(fout, opt, opt.mc(), false, var.sr(), var.dfr(), var);
  }
  if (opt.res("mx"))            // 全群の交互作用の多重比較
  {
    fout << "\n";
    put_mcmpxf(fout, opt.mx(), true, var.sr(), var.dfr(), var);
  }
  if (opt.res("mf"))            // 群1と他群の交互作用の多重比較
  {
    fout << "\n";
    put_mcmpxf(fout, opt.mf(), false, var.sr(), var.dfr(), var);
  }
  if (opt.res("mm"))            // 全群の時期別多重比較
  {
    fout << "\n";
    put_mcmpme(fout, opt.mm(), true, var);
    fout << "\n";
    put_mdcime(fout, opt, opt.mm(), true, var);
  }
  if (opt.res("me"))            // 群1と他群の時期別多重比較
  {
    fout << "\n";
    put_mcmpme(fout, opt.me(), false, var);
    fout << "\n";
    put_mdcime(fout, opt, opt.me(), false, var);
  }
}

void anova2class::put_ranova2(ostream& fout, option& opt, varclass& var)
// 繰り返しのある二元配置分散分析結果を出力する
{
  int g = var.g(), ns = var.ns();
  atabtype attype = at2r;
  var.sr(var.ssr() + var.sr());  var.dfr(var.dfsr() + var.dfr());
  var.ssr(0.0);  var.dfsr(0);
  put_anova2(fout, attype, var);
  if (opt.res("sd"))
  {
    fout << "\n";
    put_sd(fout, attype, var);
  }
  if (opt.res("ma"))            // 全群の多重比較
  {
    fout << "\n";
    put_mcmpad(fout, opt.ma(), true, var.sr(), var.dfr(), var);
    fout << "\n";
    put_mdciad(fout, opt, opt.ma(), true, var.sr(), var.dfr(), var);
  }
  if (opt.res("md"))            // 群1と他群の多重比較
  {
    fout << "\n";
    put_mcmpad(fout, opt.md(), false, var.sr(), var.dfr(), var);
    fout << "\n";
    put_mdciad(fout, opt, opt.md(), false, var.sr(), var.dfr(), var);
  }
  if (opt.res("mb"))            // 時期1と他時期の多重比較
  {
    fout << "\n";
    put_mcmpbc(fout, opt.mb(), true, var.sr(), var.dfr(), var);
    fout << "\n";
    put_mdcibc(fout, opt, opt.mb(), true, var.sr(), var.dfr(), var);
  }
  if (opt.res("mc"))            // 全時期の多重比較
  {
    fout << "\n";
    put_mcmpbc(fout, opt.mc(), false, var.sr(), var.dfr(), var);
    fout << "\n";
    put_mdcibc(fout, opt, opt.mc(), false, var.sr(), var.dfr(), var);
  }
  if (opt.res("mx"))            // 全群の交互作用の多重比較
  {
    fout << "\n";
    put_mcmpxf(fout, opt.mx(), true, var.sr(), var.dfr(), var);
  }
  if (opt.res("mf"))            // 群1と他群の交互作用の多重比較
  {
    fout << "\n";
    put_mcmpxf(fout, opt.mf(), false, var.sr(), var.dfr(), var);
  }
  if (opt.res("mm"))            // 全群の時期別多重比較
  {
    fout << "\n";
    put_mcmpme(fout, opt.mm(), true, var);
    fout << "\n";
    put_mdcime(fout, opt, opt.mm(), true, var);
  }
  if (opt.res("me"))            // 群1と他群の時期別多重比較
  {
    fout << "\n";
    put_mcmpme(fout, opt.me(), false, var);
    fout << "\n";
    put_mdcime(fout, opt, opt.me(), false, var);
  }
}

void anova2class::put_nranova2(ostream& fout, option& opt, varclass& var)
// 繰り返しのない二元配置分散分析結果を出力する
{
  int g = var.g(), ns = var.ns();
  atabtype attype = at2;
  var.sr(var.sab());  var.dfr(var.dfab());
  put_anova2(fout, attype, var);
  if (opt.res("sd")) // 群の多重比較はせず、交互作用と時期別多重比較は不可能
  {
    fout << "\n";
    put_sd(fout, attype, var);
  }
  if (opt.res("mb"))            // 時期1と他時期の多重比較
  {
    fout << "\n";
    put_mcmpbc(fout, opt.mb(), true, var.sr(), var.dfr(), var);
    fout << "\n";
    put_mdcibc(fout, opt, opt.mb(), true, var.sr(), var.dfr(), var);
  }
  if (opt.res("mc"))            // 全時期の多重比較
  {
    fout << "\n";
    put_mcmpbc(fout, opt.mc(), false, var.sr(), var.dfr(), var);
    fout << "\n";
    put_mdcibc(fout, opt, opt.mc(), false, var.sr(), var.dfr(), var);
  }
}

void anova2class::put_panova2(ostream& fout, option& opt, varclass& var)
// 対応のある二元配置分散分析結果を出力する
{
  int g = var.g(), ns = var.ns();
  atabtype attype = at2rp;
  put_anova2(fout, attype, var);
  if (opt.res("sd"))
  {
    fout << "\n";
    put_sd(fout, attype, var);
  }
  fout << "\n";
  put_ttest1(fout, opt, var.ssr(), var.dfsr(), var);
  if (opt.res("mm"))            // 時期別多重比較
  {
    fout << "\n";
    put_mcmppm(fout, opt, var);
    fout << "\n";
    put_mdcipm(fout, opt, var);
  }
}

void anova2class::put_anova2(ostream& fout, atabtype attype, varclass& var) const
// 二元配置分散分析表を出力する
{
  constexpr int LENline = 78;   // 区切り線の長さ
  fout << "・";
  if (attype == at2rm) fout << "繰り返し測定型";
  else if (attype == at2rp) fout << "対応のある";
  fout << "二元配置分散分析(two-way layout analysis of variance)\n";
  string tstr("分散分析表(ANOVA table)");
  int npad = tstr.size() - strclen(tstr, CHARwide);
  fout << setw(51 + npad) << tstr;
  fout << "\n";
  fout << "要因              平方和    自由度    平均平方和         F値   有意確率p値\n";
  putline(fout, LENline);
  double v = 0.0, f = 0.0, pr = 0.0, vr = 0.0;
  if (attype == at2rm)          // 繰り返し測定型
  {
    calvfp(v, f, pr, var.sa(), var.dfa(), var.ssr(), var.dfsr());
    tstr = "A:群";
    npad = tstr.size() - strclen(tstr, CHARwide);
    fout << left << setw(10 + npad) << tstr << right << setw(14) << var.sa();
    fout << setw(10) << var.dfa() << setw(14) << v << setw(12) << f;
    fout << setw(14) << pr << pind(pr) << "\n";
    calvfp(v, f, pr, var.ssr(), var.dfsr(), var.sr(), var.dfr());
    tstr = "個体残差";
    npad = tstr.size() - strclen(tstr, CHARwide);
    fout << left << setw(10 + npad) << tstr << right << setw(14) << var.ssr();
    fout << setw(10) << var.dfsr() << setw(14) << v << "\n";
    putline(fout, LENline);
    calvfp(v, f, pr, var.ss(), var.dfs(), var.sr(), var.dfr());
    tstr = "個体";
    npad = tstr.size() - strclen(tstr, CHARwide);
    fout << left << setw(10 + npad) << tstr << right << setw(14) << var.ss();
    fout << setw(10) << var.dfs() << setw(14) << v << setw(12) << f;
    fout << setw(14) << pr << pind(pr) << "\n";
  }
  else                          // 繰り返し測定型以外
  {
    calvfp(v, f, pr, var.sa(), var.dfa(), var.sr(), var.dfr());
    if (attype == at2) tstr = "A:個体";
    else if (attype == at2rp) tstr = "C:総平均";
    else tstr = "A:群";
    npad = tstr.size() - strclen(tstr, CHARwide);
    fout << left << setw(10 + npad) << tstr << right << setw(14) << var.sa();
    fout << setw(10) << var.dfa() << setw(14) << v << setw(12) << f;
    fout << setw(14) << pr << pind(pr) << "\n";
    if (attype == at2rp)        // 対応のある時
    {
      calvfp(v, f, pr, var.ssr(), var.dfsr(), var.sr(), var.dfr());
      tstr = "AxC:個体";
      npad = tstr.size() - strclen(tstr, CHARwide);
      fout << left << setw(10 + npad) << tstr << right << setw(14) << var.ssr();
      fout << setw(10) << var.dfsr() << setw(14) << v << setw(12) << f;
      fout << setw(14) << pr << "\n";
    }
  }
  if (attype != at2rp)         // 対応のある二元配置分散分析ではない時
  {
    vr = calvfp(v, f, pr, var.sb(), var.dfb(), var.sr(), var.dfr());
    tstr = "B:時期";
    npad = tstr.size() - strclen(tstr, CHARwide);
    fout << left << setw(10 + npad) << tstr << right << setw(14) << var.sb();
    fout << setw(10) << var.dfb() << setw(14) << v << setw(12) << f;
    fout << setw(14) << pr << pind(pr) << "\n";
  }
  if (attype == at2r || attype == at2rm || attype == at2rp) // 繰り返しのある/繰り返し測定型/対応のある時
  {
    vr = calvfp(v, f, pr, var.sab(), var.dfab(), var.sr(), var.dfr());
    if (attype == at2rp) tstr = "BxC:時期";
    else tstr = "AxB";
    npad = tstr.size() - strclen(tstr, CHARwide);
    fout << left << setw(10 + npad) << tstr << right << setw(14) << var.sab();
    fout << setw(10) << var.dfab() << 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) << var.sr();
  fout << setw(10) << var.dfr() << setw(14) << vr << "\n";
  putline(fout, LENline);
  tstr = "全体";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << var.st();
  fout << setw(10) << var.dft() << "\n";
}

void anova2class::put_sd(ostream& fout, atabtype attype, varclass& var) const
// 各要因のSDを出力する
{
  double va = (var.dfa() > 0) ? var.sa() / (double)var.dfa() : 0.0;
  double vs = (var.dfs() > 0) ? var.ss() / (double)var.dfs() : 0.0;
  double vb = (var.dfb() > 0) ? var.sb() / (double)var.dfb() : 0.0;
  double vab = (var.dfab() > 0) ? var.sab() / (double)var.dfab() : 0.0;
  double vr = (var.dfr() > 0) ? var.sr() / (double)var.dfr() : 0.0;
  double vt = (var.dft() > 0) ? var.st() / (double)var.dft() : 0.0;
  double a = (double)(var.dfa() + 1), b = (double)(var.dfb() + 1);
  double ar = (double)(var.dfs() + 1);
  double r = ar / a;
  double s2a = (va > vr) ? va - vr : 0.0;
  if (attype == at2) s2a /= b;
  else s2a /= b * r;
  double s2s = (vs > vr) ? (vs - vr) / b : 0.0;
  double s2b = (vb > vr) ? vb - vr : 0.0;
  if (attype == at2) s2b /= a;
  else s2b /= ar;
  double s2ab = (vab > vr) ? (vab - vr) / r : 0.0;
  fout << "・各要因の分散(V)と標準偏差(SD)\n";
  fout << "  要因A:" << "V=" << s2a << " SD=" << sqrt(s2a);
  if (attype == at2rm)          // 繰り返し測定型
  {
    fout << "  個体:" << "V=" << s2s << " SD=" << sqrt(s2s);
  }
  fout << "\n  要因B:" << "V=" << s2b << " SD=" << sqrt(s2b);
  if (attype != at2)            // 繰り返しあり
  {
    fout << "  交互作用:" << "V=" << s2ab << " SD=" << sqrt(s2ab) << "\n";
  }
  fout << "  残差:" << "V=" << vr << " SD=" << sqrt(vr) << "\n";
}

void anova2class::put_mcmpad(ostream& fout, mctype mc, bool ma, double sr, int dfr, varclass& var) const
// 群の多重比較を出力する(全群比較／群1との比較)
{
  constexpr int LENline = 74;   // 区切り線の長さ
  constexpr int LENhed = 9;     // ヘッダの長さ
  int g = var.g(), b = var.b();
  if (g < 2 || sr <= 0.0) return; // 出力なし
  int df1 = 1, dc = 1;
  double vr = sr / (double)dfr;
  string method = "Tukey";
  switch (mc)
  {
  case mcbon : method = "Bonferroni";
    dc = ma ? g * (g - 1) / 2 : g - 1;
    break;
  case mcdun : method = "Dunnett";
    break;
  case mcfis : method = "Fisher";
    break;
  case mcsch : method = "Scheffe";  df1 = g - 1;
  }
  fout << "・群の" << method << "型多重比較(" << method << " type multiple comparison)\n";
  if (mc == mctuk)
  {
    fout << "群       - 群                q値        群数      自由度   有意確率p値\n";
  }
  else if (mc == mcdun)
  {
    fout << "群       - 群                d値        群数      自由度   有意確率p値\n";
  }
  else
  {
    fout << "群       - 群                F値   第1自由度   第2自由度   有意確率p値\n";
  }
  putline(fout, LENline);
  int iend = ma ? g - 1 : 1;
  for (int i = 0; i < iend; ++i)
  {
    string gcodei(var.gcode(i));
    int npadi = gcodei.size() - strclen(gcodei, CHARwide);
    double ni = (double)var.n(i, b);
    double mi = var.mean(i, b);
    for (int j = i + 1; j < g; ++j)
    {
      double nj = (double)var.n(j, b);
      double dm = var.mean(j, b) - mi;
      double f = dm * dm / df1 / (1.0 / ni + 1.0 / nj) / vr;
      double r = (ni == nj) ? 0.5 : sqrt(ni / (ni + nj) * nj / (ni + nj));
      double pr = calmcpr(mc, f, g, df1, dfr, dc, r);
      string gcodej(var.gcode(j));
      int npad = gcodej.size() - strclen(gcodej, CHARwide);
      fout << left << setw(LENhed + npad) << gcodej << "- ";
      fout << setw(LENhed + npadi) << gcodei << right << setw(12) << f;
      if (mc == mctuk || mc == mcdun) fout << setw(12) << g;
      else fout << setw(12) << df1;
      fout << setw(12) << dfr << setw(14) << pr << pind(pr) << "\n";
    }
  }
  putline(fout, LENline);
}

void anova2class::put_mdciad(ostream& fout, option& opt, mctype mc, bool ma, double sr, int dfr, varclass& var) const
// 2群間の平均値の差の信頼区間を出力する(全群比較／群1との比較)
{
  constexpr int LENline = 68;   // 区切り線の長さ
  constexpr int LENhed = 9;     // ヘッダの長さ
  int g = var.g(), b = var.b();
  if (g < 2 || sr <= 0.0) return; // 出力なし
  int df1 = 1, dc = 1;
  double vr = sr / (double)dfr;
  string method = "Tukey";
  switch (mc)
  {
  case mcbon : method = "Bonferroni";
    dc = ma ? g * (g - 1) / 2 : g - 1;
    break;
  case mcdun : method = "Dunnett";
    break;
  case mcfis : method = "Fisher";
    break;
  case mcsch : method = "Scheffe";  df1 = g - 1;
  }
  fout << "・群の" << method << "型" << opt.cc() * 100.0 << "%同時信頼区間(simultaneous confidence interval)\n";
  fout << "群       - 群         平均値の差    絶対精度        下限        上限\n";
  putline(fout, LENline);
  int iend = ma ? g - 1 : 1;
  for (int i = 0; i < iend; ++i)
  {
    string gcodei(var.gcode(i));
    int npadi = gcodei.size() - strclen(gcodei, CHARwide);
    int ni = var.n(i, b);
    double mi = var.mean(i, b);
    for (int j = i + 1; j < g; ++j)
    {
      int nj = var.n(j, b);
      double dm = var.mean(j, b) - mi;
      double vd = (1.0 / ni + 1.0 / nj) * vr;
      double r = (ni == nj) ? 0.5 : sqrt(ni / (ni + nj) * nj / (ni + nj));
      double ci = calmcci(mc, opt.cc(), g, df1, dfr, dc, vd, r);
      string gcodej(var.gcode(j));
      int npad = gcodej.size() - strclen(gcodej, CHARwide);
      fout << left << setw(LENhed + npad) << gcodej << "- ";
      fout << setw(LENhed + npadi) << gcodei << right << setw(12) << dm;
      fout << setw(12) << ci << setw(12) << dm - ci << setw(12) << dm + ci << "\n";
    }
  }
  putline(fout, LENline);
}

void anova2class::put_mcmpbc(ostream& fout, mctype mc, bool mb, double sr, int dfr, varclass& var) const
// 時期の多重比較を出力する(時期1との比較／全時期比較)
{
  constexpr int LENline = 74;   // 区切り線の長さ
  constexpr int LENhed = 9;     // ヘッダの長さ
  int g = var.g(), b = var.b();
  if (b < 2 || sr <= 0.0) return; // 出力なし
  int df1 = 1, dc = 1;
  double vr = sr / (double)dfr;
  string method = "Dunnett";
  switch (mc)
  {
  case mcbon : method = "Bonferroni";
    dc = mb ? b - 1 : b * (b - 1) / 2;
    break;
  case mcfis : method = "Fisher";
    break;
  case mcsch : method = "Scheffe";  df1 = b - 1;
    break;
  case mctuk : method = "Tukey";
  }
  fout << "・時期の" << method << "型多重比較(" << method << " type multiple comparison)\n";
  if (mc == mcdun)
  { fout << "時期     - 時期              d値      時期数      自由度   有意確率p値\n"; }
  else if (mc == mctuk)
  { fout << "時期     - 時期              q値      時期数      自由度   有意確率p値\n"; }
  else
  { fout << "時期     - 時期              F値   第1自由度   第2自由度   有意確率p値\n"; }
  putline(fout, LENline);
  int iend = mb ? 1 : b - 1;
  for (int i = 0; i < iend; ++i)
  {
    double ni = (double)var.n(g, i);
    double mi = var.mean(g, i);
    for (int j = i + 1; j < b; ++j)
    {
      double nj = (double)var.n(g, j);
      double dm = var.mean(g, j) - mi;
      double f = dm * dm / df1 / (1.0 / ni + 1.0 / nj) / vr;
      double r = (ni == nj) ? 0.5 : sqrt(ni / (ni + nj) * nj / (ni + nj));
      double pr = calmcpr(mc, f, b, df1, dfr, dc, r);
      fout << left << setw(LENhed) << j + 1 << "- " << setw(9) << i + 1;
      fout << right << setw(12) << f;
      if (mc == mcdun || mc == mctuk) fout << setw(12) << b;
      else fout << setw(12) << df1;
      fout << setw(12) << dfr << setw(14) << pr << pind(pr) << "\n";
    }
  }
  putline(fout, LENline);
}

void anova2class::put_mdcibc(ostream& fout, option& opt, mctype mc, bool mb, double sr, int dfr, varclass& var) const
// 2時期間の平均値の差の信頼区間を出力する(時期1との比較／全時期比較)
{
  constexpr int LENline = 68;   // 区切り線の長さ
  constexpr int LENhed = 9;     // ヘッダの長さ
  int g = var.g(), b = var.b();
  if (b < 2 || sr <= 0.0) return; // 出力なし
  int df1 = 1, dc = 1;
  double vr = sr / (double)dfr;
  string method = "Dunnett";
  switch (mc)
  {
  case mcbon : method = "Bonferroni";
    dc = mb ? b - 1 : b * (b - 1) / 2;
    break;
  case mcfis : method = "Fisher";
    break;
  case mcsch : method = "Scheffe";
    df1 = b - 1;
    break;
  case mctuk : method = "Tukey";
  }
  fout << "・時期の" << method << "型" << opt.cc() * 100.0 << "%同時信頼区間(simultaneous confidence interval)\n";
  fout << "時期     - 時期       平均値の差    絶対精度        下限        上限\n";
  putline(fout, LENline);
  int iend = mb ? 1 : b - 1;
  for (int i = 0; i < iend; ++i)
  {
    double ni = (double)var.n(g, i);
    double mi = var.mean(g, i);
    for (int j = i + 1; j < b; ++j)
    {
      double nj = (double)var.n(g, j);
      double dm = var.mean(g, j) - mi;
      double vd = (1.0 / ni + 1.0 / nj) * vr;
      double r = (ni == nj) ? 0.5 : sqrt(ni / (ni + nj) * nj / (ni + nj));
      double ci = calmcci(mc, opt.cc(), b, df1, dfr, dc, vd, r);
      fout << left << setw(LENhed) << j + 1 << "- " << setw(LENhed) << i + 1;
      fout << right << setw(12) << dm << setw(12) << ci << setw(12) << dm - ci << setw(12) << dm + ci << "\n";
    }
  }
  putline(fout, LENline);
}

void anova2class::put_mcmpxf(ostream& fout, mctype mc, bool mcx, double sr, int dfr, varclass& var) const
// 交互作用の多重比較を出力する(全群比較／群1との比較)
{
  constexpr int LENline = 74;   // 区切り線の長さ
  constexpr int LENhed = 9;     // ヘッダの長さ
  int g = var.g(), b = var.b();
  int df1 = b - 1;
  if (g < 2 || sr <= 0.0) return; // 出力せず
  int dc = mcx ? g * (g - 1) / 2 : g - 1;
  double vr = sr / (double)dfr;
  string method = "Bonferroni";
  switch (mc)
  {
  case mcfis : method = "Fisher";  dc = 1;
    break;
  case mcsch : method = "Scheffe";  df1 = (g - 1) * (b - 1);  dc = 1;
  }
  fout << "・交互作用の" << method << "型多重比較(" << method << " type multiple comparison)\n";
  fout << "群       - 群                F値   第1自由度   第2自由度   有意確率p値\n";
  putline(fout, LENline);
  int iend = mcx ? g - 1 : 1;
  for (int i = 0; i < iend; ++i)
  {
    string gcodei(var.gcode(i));
    int npadi = gcodei.size() - strclen(gcodei, CHARwide);
    double ni = var.n(i, 0);
    double mi = var.mean(i, b);
    for (int j = i + 1; j < g; ++j)
    {
      double nj = (double)var.n(j, 0);
      double mj = var.mean(j, b);
      double m = -b * (mj - mi) * (mj - mi);
      for (int k = 0; k < b; ++k)
      {
        double dm = var.mean(j, k) - var.mean(i, k);
        m += dm * dm;
      }
      double f = m / df1 / (1.0 / ni + 1.0 / nj) / vr;
      double pr = calmcpr(mc, f, g, df1, dfr, dc);
      string gcodej(var.gcode(j));
      int npad = gcodej.size() - strclen(gcodej, CHARwide);
      fout << left << setw(LENhed + npad) << gcodej << "- ";
      fout << setw(LENhed + npadi) << gcodei << right << setw(12) << f;
      fout << setw(12) << df1 << setw(12) << dfr << setw(14) << pr << pind(pr) << "\n";
    }
  }
  putline(fout, LENline);
}

void anova2class::put_mcmpme(ostream& fout, mctype mc, bool mm, varclass& var) const
// 時期別群別平均値の多重検定を出力する(全群比較／群1との比較)
{
  constexpr int LENline = 70;   // 区切り線の長さ
  constexpr int LENhed = 9;     // ヘッダの長さ
  int g = var.g(), b = var.b(), df1 = 1;
  if (g < 2 || b < 1) return;   // 出力せず
  int dc = mm ? b * g * (g - 1) / 2 : b * (g - 1);
  if (mc == mcfis)
  {
    dc = 1;
    fout << "・時期別群別平均値の2標本t検定(two sample t test)\n";
  }
  else
  { fout << "・時期別群別平均値のBonferroni型多重検定(Bonferroni type multiple test)\n"; }
  fout << "時期    群       - 群                t値      自由度   有意確率p値\n";
  putline(fout, LENline);
  for (int k = 0; k < b; ++k)
  {
    int iend = mm ? g - 1 : 1;
    for (int i = 0; i < iend; ++i)
    {
      string gcodei(var.gcode(i));
      int npadi = gcodei.size() - strclen(gcodei, CHARwide);
      double ni = (double)var.n(i, k);
      double mi = var.mean(i, k);
      double ssi = var.ss(i, k);
      for (int j = i + 1; j < g; ++j)
      {
        double nj = (double)var.n(j, k);
        double dm = var.mean(j, k) - mi;
        double ssj = var.ss(j, k);
        int dfr = ni + nj - 2;
        double vr = (ssi + ssj) / (double)dfr * (ni + nj) / ni / nj;
        double t = (vr > 0.0) ? dm / sqrt(vr) : 0.0;
        double f = t * t;
        double pr = calmcpr(mc, f, 2, 1, dfr, dc);
        string gcodej(var.gcode(j));
        int npad = gcodej.size() - strclen(gcodej, CHARwide);
        fout << left << setw(8) << k + 1 << setw(LENhed + npad) << gcodej;
        fout << "- " << setw(LENhed + npadi) << gcodei << right;
        fout << setw(12) << t << setw(12) << dfr << setw(14) << pr << pind(pr) << "\n";
      }
    }
    putline(fout, LENline);
  }
}

void anova2class::put_mdcime(ostream& fout, option& opt, mctype mc, bool mm, varclass& var) const
// 時期別群別平均値の差の信頼区間を出力する(全群比較／群1との比較)
{
  constexpr int LENline = 76;   // 区切り線の長さ
  constexpr int LENhed = 9;     // ヘッダの長さ
  int g = var.g(), b = var.b(), df1 = 1;
  if (g < 2 || b < 1) return;   // 出力せず
  int dc = mm ? b * g * (g - 1) / 2 : b * (g - 1);
  if (mc == mcfis)
  {
    dc = 1;
    fout << "・時期別群別平均値の" << opt.cc() * 100.0 << "%信頼区間(confidence interval)\n";
  }
  else
  { fout << "・時期別群別平均値のBonferroni型" << opt.cc() * 100.0 << "%同時信頼区間(simultaneous CI)\n"; }
  fout << "時期    群       - 群         平均値の差    絶対精度        下限        上限\n";
  putline(fout, LENline);
  for (int k = 0; k < b; ++k)
  {
    int iend = mm ? g - 1 : 1;
    for (int i = 0; i < iend; ++i)
    {
      string gcodei(var.gcode(i));
      int npadi = gcodei.size() - strclen(gcodei, CHARwide);
      double ni = (double)var.n(i, k);
      double mi = var.mean(i, k);
      double ssi = var.ss(i, k);
      for (int j = i + 1; j < g; ++j)
      {
        double nj = (double)var.n(j, k);
        double dm = var.mean(j, k) - mi;
        double ssj = var.ss(j, k);
        int dfr = ni + nj - 2;
        double vr = (ssi + ssj) / (double)dfr * (ni + nj) / ni / nj;
        double ci = calmcci(mc, opt.cc(), 2, 1, dfr, dc, vr);
        string gcodej(var.gcode(j));
        int npad = gcodej.size() - strclen(gcodej, CHARwide);
        fout << left << setw(8) << k + 1 << setw(LENhed + npad) << gcodej;
        fout << "- " << setw(LENhed + npadi) << gcodei << right;
        fout << setw(12) << dm << setw(12) << ci << setw(12) << dm - ci << setw(12) << dm + ci << "\n";
      }
    }
    putline(fout, LENline);
  }
}

void anova2class::put_ttest1(ostream& fout, option& opt, double ssr, int dfsr, varclass& var)
// 総平均値の1標本t検定(対応のある二元配置の時)
{
  int g = var.g(), b = var.b();
  double me = var.mean(g, b), vr = ssr / (double)dfsr / (double)var.n(g, b);
  double t = (vr > 0.0) ? me / sqrt(vr) : 0.0;
  double f = t * t;
  ldbl pr = calmcpr(mcbon, f, 1, 1, dfsr, opt.rp());
  if (opt.rp() > 1)
  { fout << "・総平均値のBonferroni型多重検定(Bonferroni type multiple test)\n"; }
  else
  { fout << "・総平均値の1標本t検定(one sample t test)\n"; }
  fout << "  t=" << t << "  自由度" << dfsr << "  有意確率p=" << pr << pind(pr) << "\n";
  double ci = calmcci(mcbon, opt.cc(), 1, 1, dfsr, opt.rp(), vr);
  if (opt.rp() > 1)
  { fout << "  総平均値のBonferroni型" << opt.cc() * 100.0 << "%同時信頼区間="; }
  else
  { fout << "  総平均値の" << opt.cc() * 100.0 << "%信頼区間="; }
  fout << me << "±" << ci << "(" << me - ci << "-" << me + ci << ")\n";
}

void anova2class::put_mcmppm(ostream& fout, option& opt, varclass& var) const
// 時期別平均値の多重検定を出力する(対応のある二元配置の時、時期1との比較)
{
  constexpr int LENline = 62;   // 区切り線の長さ
  int g = var.g(), b = var.b();
  if (b < 2) return;            // 出力せず
  int dc = b * opt.rp();
  if (opt.mm() == mcfis)
  {
    dc = 1;
    fout << "・時期別平均値の1標本t検定(one sample t test)\n";
  }
  else
  { fout << "・時期別平均値のBonferroni型多重検定(Bonferroni type multiple test)\n"; }
  fout << "時期          平均値         t値      自由度   有意確率p値\n";
  putline(fout, LENline);
  for (int j = 0; j < b; ++j)
  {
    int dfj = var.n(g, j) - 1;
    double mj = var.mean(g, j);
    double sej = var.se(g, j);
    double tj = (sej > 0.0) ? mj / sej : 0.0;
    double f = tj * tj;
    double pr = calmcpr(opt.mm(), f, 1, 1, dfj, dc);
    fout << left << setw(8) << j + 1 << right;
    fout << setw(12) << mj << setw(12) << tj;
    fout << setw(12) << dfj << setw(14) << pr << pind(pr) << "\n";
  }
  putline(fout, LENline);
}

void anova2class::put_mdcipm(ostream& fout, option& opt, varclass& var) const
// 時期平均値の信頼区間を出力する(対応のある二元配置の時、時期1との比較)
{
  int g = var.g(), b = var.b();
  if (b < 2) return;            // 出力せず
  int dc = b * opt.rp();
  if (opt.mm() == mcfis)
  {
    dc = 1;
    fout << "・時期平均値の" << opt.cc() * 100.0 << "%信頼区間(confidence interval)\n";
  }
  else
  {
    fout << "・時期平均値のBonferroni型" << opt.cc() * 100.0 << "%同時信頼区間(simultaneous CI)\n";
  }
  fout << "時期          平均値    絶対精度        下限        上限\n";
  constexpr int LENline = 56;   // 区切り線の長さ
  putline(fout, LENline);
  for (int j = 0; j < b; j++)
  {
    int dfj = var.n(g, j) - 1;
    double mj = var.mean(g, j);
    double vrj = var.vr(g, j) / (double)var.n(g, j);
    double ci = calmcci(opt.mm(), opt.cc(), 1, 1, dfj, dc, vrj);
    fout << setw(4) << j + 1 << setw(16) << mj;
    fout << setw(12) << ci << setw(12) << mj - ci << setw(12) << mj + ci;
    fout << "\n";
  }
  putline(fout, LENline);
}

int anova2class::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 anova2class::put(option& opt) // 集計結果を出力する
{
  if (opt.res("res") == false) return; // 出力なし
  put_id(cout);
  varclass var(id.ni());        // 計算用変数クラス
  var.cal(opt, datalst);
  if (var.g() < 1)
  {
    cout << "※有効なデータなし！\n\n" << endl;
    return;
  }
  if (opt.res("stat"))
  {
    cout << "\n";
    put_stat(cout, opt, var);
  }
  if (var.ns() < 2 || var.st() <= 0.0) // 分散分析不可能
  {
    cout << "※分散分析不可能！\n\n" << endl;
    return;
  }
  if (var.ns() > var.g())       // 繰り返しあり
  {
    if (opt.res("def") || opt.res("rm")) // 繰り返し測定型二元配置分散分析
    {
      cout << "\n";
      put_rmanova2(cout, opt, var);
    }
    if (opt.res("a2"))          // 繰り返しのある二元配置分散分析
    {
      cout << "\n";
      put_ranova2(cout, opt, var);
    }
  }
  else                          // 繰り返しなし
  {
    if (opt.res("icc"))
    {
      cout << "\n";
      put_icc(cout, opt, var);
    }
    if (opt.res("def") || opt.res("a2"))
    {
      cout << "\n";
      put_nranova2(cout, opt, var);
    }
    else if (opt.res("rp"))
    {
      cout << "\n";
      put_panova2(cout, opt, var);
    }
  }
  cout << "\n" << endl;
}

int main(int argc, char *argv[]) // メインルーチン
{
  anova2class anova2;           // メインクラス変数
  option opt(argv[0]);          // オプション
  try
  {
    opt.setdef();               // オプションの初期設定
    bool fok = true;            // ファイル処理フラッグ
    for (int i = 1; i <= argc; ++i)
    {
      if (i >= argc)
      {
        if (fok) anova2.get(nullptr, opt); // 標準入力から読み込む
      }
      else if (*(argv[i]) == CHRopt)
      { opt.set(argv[i] + 1);  fok = true; }
      else
      { anova2.get(argv[i], opt);  fok = false; }
    }
    anova2.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;
}
