関連情報
ホーム > 製品 > NAG数値計算ライブラリ > C#向けNAGライブラリ > サンプルソースコード集 > 制限つき最尤法(REML)を用いた線形混合効果回帰

C#による 制限つき最尤法(REML)を用いた線形混合効果回帰

C#によるサンプルソースコード
使用関数名:g02ja

Keyword: 制限つき最尤法, REML, 線形混合効果回帰

概要

本サンプルは制限つき最尤法(REML)を用いた線形混合効果回帰のフィットを行うC#によるサンプルプログラムです。 本サンプルは以下に示されるデータについて線形混合効果回帰のフィットを行います。

線形混合効果回帰のデータ 

※本サンプルはNAG Library for .NETに含まれる関数 g02ja() のExampleコードです。本サンプル及び関数の詳細情報は g02ja のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで

入力データ

(本関数の詳細はg02ja のマニュアルページを参照)

このデータをダウンロード
g02ja Example Program Data
24 5 3 1 1
1 4 3 2 3
1 3 4 5 3 2 0 1 1
1
56 1 1 1 1
50 1 2 1 1
39 1 3 1 1
30 2 1 1 1
36 2 2 1 1
33 2 3 1 1
32 3 1 1 1
31 3 2 1 1
15 3 3 1 1
30 4 1 1 1
35 4 2 1 1
17 4 3 1 1
41 1 1 2 1
36 1 2 2 2
35 1 3 2 3
25 2 1 2 1
28 2 2 2 2
30 2 3 2 3
24 3 1 2 1
27 3 2 2 2
19 3 3 2 3
25 4 1 2 1
30 4 2 2 2
18 4 3 2 3
1.0 1.0
-1

  • 1行目はタイトル行で読み飛ばされます。
  • 2行目は観測値の数(n)、データ行列のカラム数(ncol)、固定として処理される独立変数の数(nfv)、ランダムとして処理される独立変数の数(nrv)、分散成分の数(nvpr)を指定しています。
  • 3行目は各変数のレベル数(levels)を指定しています。
  • 4行目は従属変数をもつデータ行列DATのカラム(yvid)、固定独立変数をもつDATのカラム(fvid)、ランダム独立変数をもつDATのカラム(rvid)、subject(個体)変数をもつDATのカラム(svid)、case weightをもつDATのカラム(cwid)、固定切片が含まれるかを示すフラグ(fint)、ランダム切片が含まれるかどうかを示すフラグ(rint)を指定しています。
  • 5行目はランダム変数の分散を示すフラグ(vpr)を指定しています。
  • 6〜29行目はデータ行列(dat)を指定しています。
  • 30行目はガンマの初期値(gamma)を指定しています
  • 31行目は反復の最大値(maxit)を指定しています。ゼロより小さい場合は100が使用されます。

出力結果

(本関数の詳細はg02ja のマニュアルページを参照)

この出力例をダウンロード
g02ja Example Program Results
 Fixed effects (Estimate and Standard Deviation) 
  Intercept                37.0000    4.6674
  Variable   1 Level   2    1.0000    3.5173
  Variable   1 Level   3  -11.0000    3.5173
  Variable   2 Level   2   -8.2500    2.1635
  Variable   3 Level   2    0.5000    3.0596
  Variable   3 Level   3    7.7500    3.0596

 Random Effects (Estimate and Standard  Deviation)
  Intercept for Subject Level   1     10.7631    4.4865
  Subject Level   1 Variable   1 Level   1    3.7276    3.0331
  Subject Level   1 Variable   1 Level   2   -1.4476    3.0331
  Subject Level   1 Variable   1 Level   3    0.3733    3.0331
  Intercept for Subject Level   2     -0.5269    4.4865
  Subject Level   2 Variable   1 Level   1   -3.7171    3.0331
  Subject Level   2 Variable   1 Level   2   -1.2253    3.0331
  Subject Level   2 Variable   1 Level   3    4.8125    3.0331
  Intercept for Subject Level   3     -5.6450    4.4865
  Subject Level   3 Variable   1 Level   1    0.5903    3.0331
  Subject Level   3 Variable   1 Level   2    0.3987    3.0331
  Subject Level   3 Variable   1 Level   3   -2.3806    3.0331
  Intercept for Subject Level   4     -4.5912    4.4865
  Subject Level   4 Variable   1 Level   1   -0.6009    3.0331
  Subject Level   4 Variable   1 Level   2    2.2742    3.0331
  Subject Level   4 Variable   1 Level   3   -2.8052    3.0331

  Variance Components
     1   62.3958
     2   15.3819
  SIGMA^2     =     9.3611
  -2LOG LIKE  =   119.7618
 DF          =  16

  • 2〜8行目に固定効果が出力されています。
  • 3行目に固定切片と標準誤差が出力されています。
  • 4〜8行目に固定独立変数のレベルごとにパラメータ推定と標準偏差が出力されています。
  • 10〜26行目にランダム効果が出力されています。
  • 11行目にsubject変数のレベル1のランダム切片と標準誤差が出力されています。
  • 12〜14行目にランダム独立変数のレベルごとにパラメータ推定と標準偏差が出力されています。
  • 15行目にsubject変数のレベル2のランダム切片と標準誤差が出力されています。
  • 16〜18行目にランダム独立変数のレベルごとにパラメータ推定と標準偏差が出力されています。
  • 19行目にsubject変数のレベル3のランダム切片と標準誤差が出力されています。
  • 20〜22行目にランダム独立変数のレベルごとにパラメータ推定と標準偏差が出力されています。
  • 23行目にsubject変数のレベル4のランダム切片と標準誤差が出力されています。
  • 24〜26行目にランダム独立変数のレベルごとにパラメータ推定と標準偏差が出力されています。
  • 28〜30行目に分散成分が出力されています。
  • 31行目にシグマの2乗が出力されています。
  • 32行目に対数尤度(log-likelihood)が出力されています。
  • 33行目に自由度が出力されています。

ソースコード

(本関数の詳細はg02ja のマニュアルページを参照)

※本サンプルソースコードは .NET環境用の科学技術・統計計算ライブラリである「NAG Library for .NET」の関数を呼び出します。
サンプルのコンパイル及び実行方法


このソースコードをダウンロード
//      g02ja Example Program Text
//      C# version, NAG Copyright 2008
using System;
using NagLibrary;
using System.IO;
namespace NagDotNetExamples
{
  public class G02JAE
  {
    static bool defaultdata = true;
    static string datafile = "";
    static void Main(String[] args)
    {
      if (args.Length == 1)
      {
        defaultdata = false;
        datafile = args[0];
      }
      StartExample();
    }
    public static void StartExample()
    {
      try
      {
        DataReader sr = null;
        if (defaultdata)
        {
          sr = new DataReader("exampledata/g02jae.d");
        }
        else
        {
          sr = new DataReader(datafile);
        }
        double reml,   tol; int cwid,  df,  fint,  i,  j,  k,  l,  lb,  maxit,  n,  ncol,  nff,  nfv,  nrf,  nrv,  nv,  nvpr,  rint,  svid,  warn,  yvid;
        // 
        int ifail;
        Console.WriteLine("g02ja Example Program Results");
        // 
        // 
        //      Skip heading in data file
        sr.Reset();
        // 
        //      Read in the problem size information
        sr.Reset();
        n = int.Parse(sr.Next());
        ncol = int.Parse(sr.Next());
        nfv = int.Parse(sr.Next());
        nrv = int.Parse(sr.Next());
        nvpr = int.Parse(sr.Next());
        nv = nfv + nrv;
        double[,] dat = new double[n, ncol];
        double[] gamma = new double[nvpr+2];
        int[] fvid = new int[nfv];
        int[] levels = new int[ncol];
        int[] rvid = new int[nrv];
        int[] vpr = new int[nrv];
        // 
        //      Check array sizes
        if (ncol < 2 )
        {
          Console.Write(" {0}", "Problem too small");
          goto L200;
        }
        // 
        //      Read in number of levels for each variable
        sr.Reset();
        for (i = 1; i <= ncol; i++)
        {
          levels[i - 1] = int.Parse(sr.Next());
        }
        // 
        //      Read in model information
        sr.Reset();
        yvid = int.Parse(sr.Next());
        for (i = 1; i <= nfv; i++)
        {
          fvid[i - 1] = int.Parse(sr.Next());
        }
        for (i = 1; i <= nrv; i++)
        {
          rvid[i - 1] = int.Parse(sr.Next());
        }
        svid = int.Parse(sr.Next());
        cwid = int.Parse(sr.Next());
        fint = int.Parse(sr.Next());
        rint = int.Parse(sr.Next());
        // 
        //      If no subject specified, then ignore rint
        //      Check remaining array sizes
        if (n < 1 || (nvpr + rint) < 1)
        {
          Console.Write(" {0}", "Problem too small");
          goto L200;
        }
        // 
        //      Read in the variance component flag
        sr.Reset();
        for (i = 1; i <= nrv; i++)
        {
          vpr[i - 1] = int.Parse(sr.Next());
        }
        // 
        //      Read in the Data matrix
        for (i = 1; i <= n; i++)
        {
          sr.Reset();
          for (j = 1; j <= ncol; j++)
          {
            dat[i - 1, j - 1] = double.Parse(sr.Next());
          }
        }
        // 
        //      Read in the initial values for gamma
        sr.Reset();
        if (svid == 0)
        {
          rint = 0;
        }
        for (i = 1; i <= nvpr + rint; i++)
        {
          gamma[i - 1] = double.Parse(sr.Next());
        }
        // 
        //      Read in the maximum number of iterations
        sr.Reset();
        maxit = int.Parse(sr.Next());
        //      Calculate lb
        lb = rint;
        for (i=1; i<=nrv; i++)
        {
          lb += levels[rvid[i-1]-1];
        }
        if (svid != 0) 
        {
          lb = lb*levels[svid-1];
        }
        lb += fint;
        for (i = 1; i<=nfv; i++)
        {
          lb += Math.Max(levels[fvid[i-1]-1]-1,1);
        }
        double[] b = new double[lb];
        double[] se = new double[lb];
        // 
        //      Run the analysis
        tol = 0.00e0;
        warn = 0;
        G02.g02ja(n, ncol, dat, levels, yvid, cwid, nfv, fvid, fint, nrv, rvid, nvpr, vpr, rint, svid,
                 
        gamma, out nff, out nrf, out df, out reml, b, se, maxit, tol, out warn, out ifail);
        // 
        if (ifail == 0)
        {
          //         Output results
          if (warn != 0)
          {
            Console.Write(" {0} {1}", "Warning: At least one variance component was ", "estimated to be negative and then reset to zero");
          }
          Console.Write(" {0}", "Fixed effects (Estimate and Standard Deviation)");
          Console.WriteLine(" ");
          k = 1;
          if (fint == 1)
          {
            Console.WriteLine("  {0}{1,10:f4}{2,10:f4}", "Intercept             ", b[k - 1], se[k - 1]);
            k = k + 1;
          }
          for (i = 1; i <= nfv; i++)
          {
            for (j = 1; j <= levels[fvid[i - 1] - 1]; j++)
            {
              if ((levels[fvid[i - 1] - 1] != 1) && (j == 1))
              {
                goto L40;
              }
              Console.WriteLine("  {0}{1,4}{2}{3,4}{4,10:f4}{5,10:f4}", "Variable", i, " Level", j, b[k - 1], se[k - 1]);
              k = k + 1;
              L40: ;
            }
          }
          // 
          Console.WriteLine(" ");
          Console.WriteLine(" {0} {1}", "Random Effects (Estimate and Standard",
          " Deviation)");
          if (svid == 0)
          {
            for (i = 1; i <= nrv; i++)
            {
              for (j = 1; j <= levels[rvid[i - 1] - 1]; j++)
              {
                Console.WriteLine("  {0}{1,4}{2}{3,4}{4,10:f4}{5,10:f4}", "Variable", i, " Level", j, b[k - 1], se[k - 1]);
                k = k + 1;
              }
            }
          }
          else
          {
            for (l = 1; l <= levels[svid - 1]; l++)
            {
              Console.Write("  {0}{1,4}", "Intercept for Subject Level", l);
              if (rint == 1)
              {
                Console.WriteLine("  {0,10:f4}{1,10:f4}", b[k - 1], se[k - 1]);
                k = k + 1;
              }
              for (i = 1; i <= nrv; i++)
              {
                for (j = 1; j <= levels[rvid[i - 1] - 1]; j++)
                {
                  Console.WriteLine("  {0}{1,4}{2}{3,4}{4}{5,4}{6,10:f4}{7,10:f4}", "Subject Level", l, " Variable", i, " Level",
                  j, b[k - 1], se[k - 1]);
                  k = k + 1;
                }
              }
            }
          }
          // 
          Console.WriteLine(" ");
          Console.WriteLine(" {0}", " Variance Components");
          for (i = 1; i <= (nvpr + rint); i++)
          {
            Console.WriteLine("  {0,4}{1,10:f4}", i, gamma[i - 1]);
          }
          // 
          Console.WriteLine("  {0}{1,10:f4}", "SIGMA^2     = ", gamma[nvpr + rint + 1 - 1]);
          Console.WriteLine("  {0}{1,10:f4}", "-2LOG LIKE  = ", reml);
          Console.WriteLine(" {0} {1}", "DF          = ", df);
        }
        else
        {
          Console.WriteLine(" ");
          Console.WriteLine("** g02ja failed with ifail = {0,5}", ifail);
        }
        L200: ;
        // 
      }
      catch (Exception e)
      {
        Console.WriteLine(e.Message);
        Console.Write( "Exception Raised");
      }
    }
  }
}


Results matter. Trust NAG.

Privacy Policy | Trademarks