using System; using System.Text; using System.Runtime.InteropServices; using System.Text.RegularExpressions; using System.IO; using NagLibrary; public class NagG02Functions { public static void Main() { int n = 0, ncol = 0, nfv = 0, nrv = 0, nvpr = 0; int[] levels = null, rvid = null, fvid = null, vpr = null; int svid = 0, cwid = 0, fint = 0, rint = 0; int fl = 0, fnlevel = 0, rnlevel = 0, lgamma = 0, lb = 0, tddat = 0; int maxit = 0, warnp = 0; double tol = 0; int nff = 0, nrf = 0, df = 0; double like = 0; double[] dat = null; double[] gamma = null, b = null, se = null; int i = 0, j = 0, k = 0, l = 0; int yvid = 0; nag_declarations.NagError fail = new nag_declarations.NagError(); fail.message = new byte[512]; Encoding enc = Encoding.ASCII; String line; Regex spaces = new Regex(@"\s+"); FileStream g02jaceDataHandle = new FileStream("g02jace.d", FileMode.Open, FileAccess.Read); StreamReader g02jaceDataReader = new StreamReader(g02jaceDataHandle); // skip heading in data file line = g02jaceDataReader.ReadLine(); // Read in the problem size information line = g02jaceDataReader.ReadLine(); string[] field = spaces.Split(line.Trim()); n = int.Parse(field[0]); ncol = int.Parse(field[1]); nfv = int.Parse(field[2]); nrv = int.Parse(field[3]); nvpr = int.Parse(field[4]); // Check problem size if (n < 0 || ncol < 0 || nfv < 0 || nrv < 0 || nvpr < 0) { Console.WriteLine("Invalid problem size, at least one of n, ncol, nfv, nrv or nvpr is < 0"); Environment.Exit(1); } // allocate some memory levels = new int[ncol]; fvid = new int[nfv]; rvid = new int[nrv]; vpr = new int[nrv]; // read the number of levels for each variable line = g02jaceDataReader.ReadLine(); field = spaces.Split(line.Trim()); for (i = 0; i < ncol; ++i) levels[i] = int.Parse(field[i]); // read in model information line = g02jaceDataReader.ReadLine(); field = spaces.Split(line.Trim()); yvid = int.Parse(field[0]); for (i = 0; i < nfv; ++i) fvid[i] = int.Parse(field[i + 1]); for (i = 0; i < nrv; ++i) rvid[i] = int.Parse(field[nfv + i + 1]); svid = int.Parse(field[nfv + nrv + 1]); cwid = int.Parse(field[nfv + nrv + 2]); fint = int.Parse(field[nfv + nrv + 3]); rint = int.Parse(field[nfv + nrv + 4]); // read in the variance component flag line = g02jaceDataReader.ReadLine(); field = spaces.Split(line.Trim()); for (i = 0; i < nrv; ++i) vpr[i] = int.Parse(field[i]); // if no subject specified, ignore rint if (svid == 0) rint = 0; // count the number of levels in the fixed parameters for (i = 0, fnlevel = 0; i < nfv; ++i) { fl = levels[fvid[i] - 1] - 1; fnlevel += (fl < 1) ? 1 : fl; } if (fint == 1) fnlevel++; // count the number of levels in the random parameters for (i = 0, rnlevel = 0; i < nrv; ++i) { rnlevel += levels[rvid[i] - 1]; } if (rint != 0) rnlevel++; // calculate the sizes of the output arrays if (rint == 1) { lgamma = nvpr + 2; } else { lgamma = nvpr + 1; } if (svid != 0) { lb = fnlevel + levels[svid - 1] * rnlevel; } else { lb = fnlevel + rnlevel; } tddat = ncol; // allocate remaining memory dat = new double[n * ncol]; gamma = new double[lgamma]; b = new double[lb]; se = new double[lb]; // read the data matrix for (i = 0; i < n; ++i) { line = g02jaceDataReader.ReadLine(); field = spaces.Split(line.Trim()); for (j = 0; j < ncol; ++j) dat[i * ncol + j] = double.Parse(field[j]); } // read the initial values for gamma line = g02jaceDataReader.ReadLine(); field = spaces.Split(line.Trim()); for (i = 0; i < lgamma - 1; ++i) gamma[i] = double.Parse(field[i]); // read the maximum number of iterations line = g02jaceDataReader.ReadLine(); field = spaces.Split(line.Trim()); maxit = int.Parse(field[0]); g02jaceDataReader.Close(); g02jaceDataHandle.Close(); /* nag_reml_mixed_regsn (g02jac). * Linear mixed effects regression using Restricted Maximum Likelihood (REML) */ nag_declarations.g02jac(n, ncol, dat, tddat, levels, yvid, cwid, nfv, fvid, fint, nrv, rvid, nvpr, vpr, rint, svid, gamma, ref nff, ref nrf, ref df, ref like, lb, b, se, maxit, tol, ref warnp, ref fail); if (fail.code == 0) { if (warnp != 0) { Console.WriteLine("Warning: At least one variance component was " + "estimated to be negative and then reset to zero"); } Console.WriteLine("Fixed effects (Estimate and Standard Deviation)"); Console.WriteLine(""); k = 0; if (fint == 1) { Console.WriteLine("Intercept " + b[k].ToString(" #####.0000") + se[k].ToString(" #####.0000")); ++k; } for (i = 0; i < nfv; ++i) { for (j = 1; j <= levels[fvid[i] - 1]; ++j) { if (levels[fvid[i] - 1] != 1 && j == 1) continue; int ii = i + 1; Console.WriteLine("Variable " + ii.ToString() + " Level " + j.ToString("#### ") + b[k].ToString("#####.0000 ") + se[k].ToString("#####.0000")); ++k; } } Console.WriteLine(""); Console.WriteLine("Random Effects (Estimate and Standard Deviation)"); if (svid == 0) { for (i = 0; i < nrv; ++i) { for (j = 1; j <= levels[rvid[i] - 1]; ++j) { int ii = i + 1; Console.WriteLine("Variable " + ii.ToString() + " Level " + j.ToString() + b[k].ToString("####.0000 ") + se[k].ToString("####.0000 ") ); ++k; } } } else { for (l = 1; l <= levels[svid - 1]; ++l) { if (rint == 1) { Console.WriteLine("Intercept for Subject Level " + l.ToString("## ") + " " + b[k].ToString("####.0000 ") + se[k].ToString("####.0000")); ++k; } for (i = 0; i < nrv; ++i) { for (j = 1; j <= levels[rvid[i] - 1]; ++j) { int ii = i + 1; Console.WriteLine("Subject Level " + l.ToString() + " Variable " + ii.ToString("## ") + " Level " + j.ToString("## ") + b[k].ToString("####.0000 ") + se[k].ToString("####.0000")); ++k; } } } } Console.WriteLine(""); Console.WriteLine("Variance Components"); for (i = 1; i <= nvpr + rint; ++i) { Console.WriteLine(i.ToString("#### ") + gamma[i - 1].ToString("#####.####")); } Console.WriteLine("SIGMA^2 = " + gamma[nvpr + rint].ToString("#####.####")); Console.WriteLine("-2LOG LIKE = " + like.ToString("#####.####")); Console.WriteLine("DF = " + df.ToString()); } else { Console.WriteLine("Error from nag_reml_mixed_regsn (g02jac):"); Console.WriteLine(enc.GetString(fail.message)); } Environment.Exit(fail.code); } }