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
- １行目はタイトル行で読み飛ばされます。
- ２行目は観測値の数(n)、データ行列のカラム数(ncol)、固定として処理される独立変数の数(nfv)、ランダムとして処理される独立変数の数(nrv)、分散成分の数(nvpr)を指定しています。
- ３行目は各変数のレベル数(levels)を指定しています。
- ４行目は従属変数をもつデータ行列DATのカラム(yvid)、固定独立変数をもつDATのカラム(fvid)、ランダム独立変数をもつDATのカラム(rvid)、subject（個体）変数をもつDATのカラム(svid)、case weightをもつDATのカラム(cwid)、固定切片が含まれるかを示すフラグ(fint)、ランダム切片が含まれるかどうかを示すフラグ(rint)を指定しています。
- ５行目はランダム変数の分散を示すフラグ(vpr)を指定しています。
- ６〜２９行目はデータ行列(dat)を指定しています。
- ３０行目はガンマの初期値(gamma)を指定しています
- ３１行目は反復の最大値(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
- ２〜８行目に固定効果が出力されています。
- ３行目に固定切片と標準誤差が出力されています。
- ４〜８行目に固定独立変数のレベルごとにパラメータ推定と標準偏差が出力されています。
- １０〜２６行目にランダム効果が出力されています。
- １１行目にsubject変数のレベル１のランダム切片と標準誤差が出力されています。
- １２〜１４行目にランダム独立変数のレベルごとにパラメータ推定と標準偏差が出力されています。
- １５行目にsubject変数のレベル２のランダム切片と標準誤差が出力されています。
- １６〜１８行目にランダム独立変数のレベルごとにパラメータ推定と標準偏差が出力されています。
- １９行目にsubject変数のレベル３のランダム切片と標準誤差が出力されています。
- ２０〜２２行目にランダム独立変数のレベルごとにパラメータ推定と標準偏差が出力されています。
- ２３行目にsubject変数のレベル４のランダム切片と標準誤差が出力されています。
- ２４〜２６行目にランダム独立変数のレベルごとにパラメータ推定と標準偏差が出力されています。
- ２８〜３０行目に分散成分が出力されています。
- ３１行目にシグマの２乗が出力されています。
- ３２行目に対数尤度（log-likelihood）が出力されています。
- ３３行目に自由度が出力されています。
ソースコード
（本関数の詳細は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");
}
}
}
}