/*------------------------------------------------------------------------------*
 * File Name:	SurvivalAnalysis.c: Contains OriginC source code for all		*
 * 	Survival Analysis features.													*
 * Creation: GJL 1/28/2002														*
 * Purpose: Origin C file														*
 * Copyright ( c ) OriginLab Corp.	2000-2002									*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:															*
 *	GJL 1/28/2002: Split off source code from OStat.c							*
 *	CP/GJL 7/21/02 v7.0348 OSTAT_COMPILE_PROBLEM_WITH_TYPE_ERROR_MSG			*
 *		I have added Type_ErrorMsg1 function to replace all those that need		*
 *		second string arguments													*
 *	GJL 11/14/02 v7.0434 USE_OC_WKS_TYPE_NOT_LT									*
 *------------------------------------------------------------------------------*/

////////////////////////////////////////////////////////////////////////////////////
//
// System includes
#include <Origin.h>
//#include <OC_const.h>
//#include <string.h>
//#include <data.h>
//#include <math.h>
//#include <utilities.h>
//#include <stdio.h>
//#include <page.h>

// Includes definitions of SurvivalAnalysis constants, structures, and non-localized strings.
//  Also includes NAG definitions and prototypes of NAG functions used by SurvivalAnalysis.
#include "SurvivalAnalysis.h"

// Includes prototypes of Application Utilities.
#include "App_Utils.h"

// Includes definitions of all SurvivalAnalysis.c localized strings. $ causes CodeBuilder to look for 
// correct version of Local.h depending on Windows Regional settings
#include "$Local.h"
//
////////////////////////////////////////////////////////////////////////////////////

/**
		Implements the Kaplan-Meier (Product-Limit) Estimator. The computational
		engine is the NAG function nag_prod_limit_surviv_fn (g12aac). Additional
		reference is Chapter 2 of Applied Survival Analysis: Regression Modeling
		of Time to Event Data (1998), Hosmer, D.W. and Lemeshow, S. (1998) John
		Wiley and Sons Inc., New York, NY.
	Example:
		See [KaplanMeier] section of KaplanMeier.OGS for sample call.
	Parameters:
		nPts=Number of data points in Time and Censor datasets
		nResults=Bitwise integer specifying output options
			1=Output Survivorship Function in Resuls Log
			2=Send output to temporary worksheet in addition to Results Log
			4=Create plot of Survivorship Function (implemented in KaplanMeier.OGS)
			8=Plot Errors or Confidence Limits along with Survivorship Function  
		strSurvivalWKS=Name of temporary worksheet holding sorted copies of all
			required input (Time and Censor) and output datasets   
	Return:
		Depending on nResults, outputs Summary of Event and Censored Times, the 
		Survivorship Function, and Quartile Estimates to the Results Log and to
		a worksheet. Returns SA_NO_ERROR on successful exit or an SA_ error code.
*/
int saKaplanMeier( int nPts, int nResults, double dConLevel, string strSurvivalWKS )
{
	int iSize, iErr, nd;
	string strDatasetName;
	BOOL bErr;

	// Read Time dataset (including failure and censored times) into vector to pass to NAG
	strDatasetName = strSurvivalWKS + SA_TIME_DATASET_NAME;	
	Dataset dsTIME( strDatasetName );
	iSize = dsTIME.GetSize();
	vector<double> vTIME;
	vTIME.SetSize( iSize );
	vTIME = dsTIME;	
	dsTIME.Detach(); // No longer need input dataset

	// Read Censor dataset into vector to pass to NAG	
	strDatasetName = strSurvivalWKS + SA_CENSOR_DATASET_NAME;
	Dataset dsCENSOR( strDatasetName );
	iSize = dsCENSOR.GetSize();
	vector<int> vCENSOR;
	vCENSOR.SetSize( iSize );
	vCENSOR = dsCENSOR;
	
	// Perform basic stats on Censor dataset (all 0's and 1's) for summary output
	Dataset *pdsCENSOR = &dsCENSOR;
	BasicStats bsStat;
	BasicStats *pbsStat = &bsStat;
	bErr = Data_sum( pdsCENSOR, pbsStat );
	if( bErr == FALSE )
		return SA_DATA_SUM_ERROR1;
	
	dsCENSOR.Detach(); // No longer need input dataset

	// Contains distinct ordered failure times on output from NAG (does not include duplicates or censored times)
	vector<double> vTP;
	vTP.SetSize( iSize );

	// Contains Kaplan-Meier estimate of Survival Probability for each time in vTP on output from NAG 
	vector<double> vP;
	vP.SetSize( iSize );

	// Contains Kaplan-Meier estimate of error (standard deviation) of each probability in vP on output from NAG
	vector<double> vPSIG;
	vPSIG.SetSize( iSize );

	// From <NAG\OCN_g12.h>: g12aac nag_prod_limit_surviv_fn: Compute Kaplan-Meier estimate of
	// Survivorship Function. Frequency is NULL (count each point once; nd is the number of distinct
	// failure times output in vTP.
	iErr = nag_prod_limit_surviv_fn( nPts, vTIME, vCENSOR, NULL, &nd, vTP, vP, vPSIG );
	if( iErr != NE_NOERROR )
		return SA_PROD_LIMIT_ERROR;

	// Compute Quartile Estimates (Percentiles)
	double ardPercentile[3];
	iErr = saKaplanMeier_Compute_Percentiles( nd, vTP, vP, ardPercentile );
	if( iErr != SA_NO_ERROR )
		return iErr;
	
	// Compute Errors or Confidence Limits (upper and lower) for Percentiles
	double ardPerLLim[3];
	double ardPerULim[3];
	iErr = saKaplanMeier_Compute_Con_Limits( dConLevel, nd, vTP, vP, vPSIG, ardPerLLim, ardPerULim );
	if( iErr != SA_NO_ERROR )
		return iErr;

	// Output results in Results Log
	iErr = saKaplanMeier_Output_to_ResultsLog( nResults, nPts, vTIME, vCENSOR, pbsStat->total, vP, vPSIG, dConLevel,
		 ardPercentile, ardPerLLim, ardPerULim );
	if( iErr != SA_NO_ERROR )
		return iErr;

	// Output results in worksheet? If user wants worksheet or plot write everything out to worksheet
	if( nResults > 1 )
	{
		// Function to output results in worksheet
		iErr = saKaplanMeier_Output_to_Worksheet( strSurvivalWKS, nPts, vTIME, vCENSOR, nd , vP, vPSIG, dConLevel,
			 ardPercentile, ardPerLLim, ardPerULim );
		if( iErr != SA_NO_ERROR )
			return iErr;
	}

	return SA_NO_ERROR;
}

/**
		Computes Quartile Estimates (the 25th, 50th, and 75th Percentiles) for the Kaplan-Meier (Product-Limit)
		Estimator function.
	Example:
		See function saKaplanMeier above for sample call.
	Parameters:
		nd=Number of distinct failure times in vTP (nd==size of vTP)
		vTP=Vector containing ordered failure times output by NAG (does not include duplicate or censored times)
		vP=Vector containing Kaplan-Meier estimate of Survival Probability for each time in vTP output by NAG
		ardPercentile=Three element array containing the 25th, 50th and 75th Percentiles
	Return:
		Array of doubles holding the 25th, 50th, and 75th percentiles of the Kaplan-Meier estimate
		of the Survivorship Function. Returns SA_NO_ERROR on successful exit.
*/
int saKaplanMeier_Compute_Percentiles( int nd, vector<double> &vTP, vector<double> &vP, double *ardPercentile )
{
	int ii, jj, iPercInterp, ivPofii, iP;
	double dPercInterp;

	// Get LabTalk system.graph.percinterp setting (from graph tab in Tools:Options DB).  ==1 when
	//  percentiles falling between two elements should be averaged. 
	LT_evaluate( "system.graph.percinterp", &dPercInterp );
	iPercInterp = dPercInterp;

	// Initialize values to negative time as flag for missing value
	ardPercentile[0] = -1;
	ardPercentile[1] = -1;
	ardPercentile[2] = -1;	
	jj = 0;
	iP = 75000; // Integerized percentile threshold for 25th percentile
	
	// Loop on three percentiles (25th, 50th, and 75th) as long as all non-missing value
	//  surviviorship function values have not been exausted
	for( ii = 0; ( ii < nd ) && ( jj < 3 ); ii++ )
	{
		// Allow for equality within tolerance by rounding to 10,000'sth place and integerizing
		ivPofii = vP[ii] * 100000 + 0.5;

		// If survivorship function value is less than or equal to percentile threshold...
		if( ivPofii <= iP )
		{
			// If equal to percentile threshold and interpolate option is enabled then try to interpolate...
			if( ivPofii == iP && iPercInterp )
			{
				// If next element exists to interpolate with then interpolate...
				if( ii + 1 < nd )  
					ardPercentile[jj] = ( vTP[ii] + vTP[ii + 1] ) / 2; // Interpolate
				// Else leave as missing value by not setting
			}
			else
				ardPercentile[jj] = vTP[ii]; // Don't interpolate
			ii--; // Decrement loop index to start looking for next percentile with current time value
			jj++; // Got current percentile...now get next one
			iP = 75000 - jj * 25000; // Compute integerized threshold for next percentile
		}
	}

	return SA_NO_ERROR;
}

/**
		Computes Upper and Lower Confidence Limits for the Percentiles of the Kaplan-Meier (Product-Limit)
		Estimator function.
	Example:
		See function saKaplanMeier above for sample call.
	Parameters:
		dConLevel=Confidence level
		nd=Number of distinct failure times in vTP (nd==size of vTP)
		vTP=Vector containing distinct ordered failure times output by NAG (does not include duplicate or 
			censored times)
		vP=Vector containing Kaplan-Meier estimate of Survival Probability for each time in vTP output by NAG
		vPSIG=Vector containing standard deviation of each probability in vP output by NAG	
		ardPerLLim=Three element array containing a lower confidence limit for the 25th, 50th, and 75th
			percentiles
		ardPerULim=Three element array containing an upper confidence limit for the 25th, 50th, and 75th
			percentiles
	Return:
		Two arrays of doubles holding the upper and lower confidence limits for the 25th, 50th, and 75th
		percentiles of the Kaplan-Meier estimate of the Survivorship Function. Returns SA_NO_ERROR on
		successful exit.
*/
int saKaplanMeier_Compute_Con_Limits( double dConLevel, int nd, vector<double> &vTP, vector<double> &vP,
	 vector<double> &vPSIG, double *ardPerLLim, double *ardPerULim )
{
	int ii, jj;
	double dP, dChi;
	NagError neErr;
	NagError *pneFail = &neErr;		// NAG error structure

	// Initialize all to -1 for later output as NANUM if not successfully computed
	ardPerLLim[0] = -1;
	ardPerLLim[1] = -1;
	ardPerLLim[2] = -1;	
	ardPerULim[0] = -1;
	ardPerULim[1] = -1;
	ardPerULim[2] = -1;

	// From <NAG\OCN_g01.h>: g01fcc nag_deviates_chi_sq: Compute Chi-Square deviate 
	dChi = nag_deviates_chi_sq( dConLevel, 1, pneFail );
	if( pneFail->code != NE_NOERROR )
		return SA_CHI_SQ_ERROR1;
	
	// Loop on three percentiles 25th, 50th, and 75th
	for( jj = 0; jj < 3; jj++ )
	{
		// Compute non-integerized threshold for current percentile
		dP = 0.75 - jj * 0.25;
		
		// Loop through all distict failure times
		for( ii = 0; ii < nd; ii++ )
		{
			// If condition is TRUE then time value may be upper or lower limit 
			if( fabs( vP[ii] - dP ) <= sqrt( dChi * pow( vPSIG[ii],2 ) ) )
			{
				if( ardPerLLim[jj] < 0 ) // Lower confidence limit is the first time the condition is TRUE
					ardPerLLim[jj] = vTP[ii];
				if( ii + 1 < nd ) // Upper confidence limit is the next time value the last time the condition is TRUE 
					ardPerULim[jj] = vTP[ii + 1];				
			} // End if condition TRUE
		} // End loop through all failure times
	} // End loop on three percentiles
	
	return SA_NO_ERROR;
}

/**
		Outputs the results of the Kaplan-Meier (Product-Limit) Estimator to the Results Log.
	Example:
		See the function saKaplanMeier above for sample call.
	Parameters:
		nResults=Bitwise flag (bit 1) indicating whether to output the Survivorship Function
			to the Results Log
		nPts=Number of data points in Time and Censor datasets
		vTIME=Vector of input Time values
		vCENSOR=Vector of input Censor values
		dCensoredPts=Number of censored data points
		vP=Vector containing Survivorship Function values
		vPSIG=Vector containing Survivorship Function Error (standard deviation) values
		dConLevel=Confidence Level
		ardPercentile=Array containing 25th, 50th, and 75th Percentile values
		ardPerLLim=Array containing Lower Limits for Percentile values
		ardPerULim=Array containing Upper Limits for Percentile values
	Return:
		Outputs the results of the Kaplan-Meier Estimator (including a Main Header and Summary
		table, a Survivorship Function with Errors table, and a Quartile Estimates table) to the
		Results Log and returns SA_NO_ERROR on successful exit or an SA_ error code.
*/
int saKaplanMeier_Output_to_ResultsLog( int nResults, int nPts, vector<double> &vTIME, vector<int> &vCENSOR, double dCensoredPts,
	 vector<double> &vP, vector<double> &vPSIG, double dConLevel, double *ardPercentile, double *ardPerLLim,
	 double *ardPerULim )
{
	int ii, jj, iErr;
	string strDB;
	string strTimeData;
	string strCensorData;
	string strCensorValue;
	char szTemp[40];
	double dPrevUncensoredTime;

	string strOut;	
	int iRow, iCol;
	string strTablePath;
	string strTableName;

	// *** Output Survival Analysis Main Headers and Summary Table ***
	// Get name of Time dataset from dialog box	
	strDB = "%A=KaplanMeier!TimeDataEBX$";
	LT_execute( strDB );
	LT_get_str( "%A", szTemp, 40 );
	strTimeData = szTemp;

	// Get name of Censor dataset from dialog box
	strDB = "%A=KaplanMeier!CensorDataEBX$";
	LT_execute( strDB );
	LT_get_str( "%A", szTemp, 40 );
	strCensorData = szTemp;

	// Get Censor Value from dialog box
	strDB = "%A=KaplanMeier!CensorValueEBX$";
	LT_execute( strDB );	
	LT_get_str( "%A", szTemp, 40 );
	strCensorValue = szTemp;

	// Get name of dialog box (or feature) from Local.h
	strDB.Format( KM_MAIN_HDR );
	
	// Output survival analysis main headers and summary table for Kaplan-Meier Estimator
	iErr = saSurvival_Analysis_ResultsLog_Summary( strDB, strTimeData, strCensorData, strCensorValue,
		 nPts, dCensoredPts );
	if( iErr != SA_NO_ERROR )
		return iErr;

	// *** Output Survivorship Function and Errors ***
	if( nResults & 1)
	{
		Worksheet wksSurvFuncTable;
		
		// Create and attach to a temporary worksheet 
		strTablePath = GetFullPath( "saKMSurvFunc.OGW", SA_OGW_SUBFOLDER, TRUE );
		if( !strTablePath.IsFile() )
		{
			Type_ErrorMsg1( SA_FILE_NOT_FOUND_ERROR_MSG, strTablePath );
			return SA_ABORT_NO_ERROR_MSG;
		}
		wksSurvFuncTable.Open( strTablePath, CREATE_TEMP );
		strTableName = wksSurvFuncTable.GetPage().GetName();
		
		// Rows and columns are indexed from 0
		
		// Insert a blank line in survivorship function table for each data point + 1
		// (+ 1 for time 0) starting at row 5
		iRow = 4;
		Type_Insert_Blank_Lines_in_Wks( wksSurvFuncTable, iRow, nPts + 1 );
		
		// Output time value = 0 in column 2
		iCol = 1;
		wksSurvFuncTable.SetCell( iRow, iCol, 0.0 );
			
		// Output survivorship function value = 1 in column 3
		iCol++;
		wksSurvFuncTable.SetCell( iRow, iCol, 1.0 );
		
		// Output survivorship function error value = 0 in column 4
		iCol++;
		wksSurvFuncTable.SetCell( iRow, iCol, 0.0 );

		// Output non-zero time values starting in row 5
		iRow++;
		
		// Loop on all non-zero time values
		jj = -1;
		dPrevUncensoredTime = -1;
		for( ii = 0;ii < nPts; iRow++, ii++ )
		{
			// Output time, survivorship function, and survivorship function error values
			// Output time value in column 2
			iCol = 1;
			wksSurvFuncTable.SetCell( iRow, iCol, vTIME[ii] );
			
			if( vCENSOR[ii] == 1 ) // If censored value...
			{
				// Output missing value for survivorship function value in column 3
				iCol++;
				wksSurvFuncTable.SetCell( iRow, iCol, NANUM );
				
				// Output missing value for survivorship function error value in column 4
				iCol++;
				wksSurvFuncTable.SetCell( iRow, iCol, NANUM );
			}
			else // Else, not censored...
			{
				// Output survivorship function and survivorship function error value
				// Must index into vP and vPSIG differently than vTIME because only non-repeated and uncensored
				// survival/error values are in vP/vPSIG while all times are in vTIME 
				if( jj < 0 || vTIME[ii] != dPrevUncensoredTime )
					jj++;
				dPrevUncensoredTime = vTIME[ii];
				
				// Output survivorship function value in column 3
				iCol++;				
				wksSurvFuncTable.SetCell( iRow, iCol, vP[jj] );

				// Output survivorship function error value in column in column 4
				iCol++;				
				wksSurvFuncTable.SetCell( iRow, iCol, vPSIG[jj] );
			}
		}
		
	// GJL 11/14/02 v7.0434 USE_OC_WKS_TYPE_NOT_LT
	// Send output suvivorship function worksheet to Results Log
	wksSurvFuncTable.Type( TYPETARGET_OUTPUTLOG, NULL, 5 );
	
	}
	
	// *** Output Quartile Estimates ***	
	Worksheet wksQuartEstTable;
		
	// Create and attach to a temporary worksheet 
	strTablePath = GetFullPath( "saKMQuartEst.OGW", SA_OGW_SUBFOLDER, TRUE );
	if( !strTablePath.IsFile() )
	{
		Type_ErrorMsg1( SA_FILE_NOT_FOUND_ERROR_MSG, strTablePath );
		return SA_ABORT_NO_ERROR_MSG;
	}
	wksQuartEstTable.Open( strTablePath, CREATE_TEMP );
	strTableName = wksQuartEstTable.GetPage().GetName();
	
	// Rows and columns are indexed from 0
	
	// Output confidence level in row 3 and column 4
	iRow = 2;
	iCol = 3;
	strOut = LocalizeDouble( dConLevel*100, SA_PERCENT_FORMAT ) + KM_QUARTILE_HDR2;
	wksQuartEstTable.SetCell( iRow, iCol, strOut ); 

	// Loop on three percentiles (0.25, 0.50, and 0.75) starting in row 6
	iRow = 5;
	for( ii = 0; ii < 3; iRow++, ii++ )
	{
		// Output the percentile in column 3
		iCol = 2;
		if( ardPercentile[ii] >= 0 )
			wksQuartEstTable.SetCell( iRow, iCol, ardPercentile[ii] ); 

		// Output the percentile lower limit in column 4
		iCol++;
		if( ardPerLLim[ii] >= 0 )
			wksQuartEstTable.SetCell( iRow, iCol, ardPerLLim[ii] ); 
		
		// Output the percentile upper limit in column 5
		iCol++;
		if( ardPerULim[ii] >= 0 )
			wksQuartEstTable.SetCell( iRow, iCol, ardPerULim[ii] ); 
	}

	// GJL 11/14/02 v7.0434 USE_OC_WKS_TYPE_NOT_LT
	// Send output quartile estimates worksheet to Results Log
	wksQuartEstTable.Type( TYPETARGET_OUTPUTLOG, NULL, 6 );
	
	return SA_NO_ERROR;
}

/**
		Outputs the results of the Kaplan-Meier (Product-Limit) Estimator to a worksheet.
	Example:
		See the function saKaplanMeier above for sample call.
	Parameters:
		strSurvivalWKS=Name of output worksheet
		nPts=Number of data points in Time and Censor datasets
		vTIME=Vector of input Time values
		vCENSOR=Vector of input Censor values
		nd=Number of distinct failure times (nd==size of vTP)
		vP=Vector containing Survivorship Function values
		vPSIG=Vector containing Survivorship Function Error (standard deviation) values
		dConLevel=Confidence Level
		ardPercentile=Array containing 25th, 50th, and 75th Percentile values
		ardPerLLim=Array containing Lower Limits for Percentile values
		ardPerULim=Array containing Upper Limits for Percentile values
	Return:
		Outputs the results of the Kaplan-Meier Estimator (the Survivorship Function with Errors table and Upper
		and Lower Confidence Limits, and a Quartile Estimates table) to a worksheet. Returns SA_NO_ERROR on
		successful exit or an SA_ error code.
*/
int saKaplanMeier_Output_to_Worksheet( string strSurvivalWKS, int nPts, vector<double> &vTIME, vector<int> &vCENSOR,
	 int nd, vector<double> &vP, vector<double> &vPSIG, double dConLevel, double *ardPercentile, double *ardPerLLim,
	 double *ardPerULim )
{
	int ii, jj;
	double dZ, dPrevUncensoredTime;
	string strDatasetName;
	NagError neErr;
	NagError *pneFail = &neErr;		// NAG error structure

	// Output Time dataset	
	strDatasetName = strSurvivalWKS + SA_TIME_DATASET_NAME;
	Dataset dsTIME( strDatasetName );
	dsTIME.SetSize( nPts + 1 );
	
	// Output Survivorship Function dataset
	strDatasetName = strSurvivalWKS + SA_SURVIVAL_DATASET_NAME;
	Dataset dsSURVIVAL( strDatasetName );
	dsSURVIVAL.SetSize( nPts + 1 );

	// Output Survivorship Function Error (Standard Deviation) dataset
	strDatasetName = strSurvivalWKS + SA_STDERROR_DATASET_NAME;	
	Dataset dsERROR( strDatasetName );
	dsERROR.SetSize( nPts + 1 );

	// Output Survivorship Function Lower Confidence Limit dataset
	strDatasetName = strSurvivalWKS + KM_SURLLIM_DATASET_NAME;	
	Dataset dsSURLLIM( strDatasetName );
	dsSURLLIM.SetSize( nPts + 1 );
	
	// Output Survivorship Function Upper Confidence Limit dataset
	strDatasetName = strSurvivalWKS + KM_SURULIM_DATASET_NAME;
	Dataset dsSURULIM( strDatasetName );
	dsSURULIM.SetSize( nPts + 1 );
	
	// Output Percentile dataset
	strDatasetName = strSurvivalWKS + KM_PERCENT_DATASET_NAME;
	Dataset dsPERCENT( strDatasetName );
	dsPERCENT.SetSize( 3 );
	
	// Output Survivorship Function Percentile Estimate dataset
	strDatasetName = strSurvivalWKS + KM_ESTIMATE_DATASET_NAME;
	Dataset dsESTIMATE( strDatasetName );
	dsESTIMATE.SetSize( 3 );
	
	// Output Survivorship Function Percentile Estimate Lower Limit dataset
	strDatasetName = strSurvivalWKS + KM_PERLLIM_DATASET_NAME;
	Dataset dsPERLLIM( strDatasetName );
	dsPERLLIM.SetSize( 3 );

	// Output Survivorship Function Percentile Estimate Upper Limit dataset
	strDatasetName = strSurvivalWKS + KM_PERULIM_DATASET_NAME;	
	Dataset dsPERULIM( strDatasetName );
	dsPERULIM.SetSize( 3 );

	// *** Add time=0 to survivorship function, compute confidence limits ( upper and lower limits ) ***
	// *** for survivorship function, and copy all to output datasets in Survival worksheet          ***
	
	// From <NAG\OCN_g01.h>: g01fac nag_deviates_normal: Compute Normal Deviate for Survivorship Function Upper and
	//  Lower Limits
	dZ = nag_deviates_normal( Nag_LowerTail, ( dConLevel+1 )/2, pneFail );
	if( pneFail->code != NE_NOERROR )
		return SA_NORMAL_ERROR;

	dPrevUncensoredTime = -1;
	jj = nd;
	for( ii = nPts; ii > 0; ii-- )      // Loop backward through all input time values
	{
		if( vCENSOR[ii-1] == 1 ) 		// If time value ii-1 was censored and not a failure time...
		{
			dsTIME[ii] = vTIME[ii - 1]; // Output time value anyway 
			dsSURVIVAL[ii] = NANUM;     //  and NANUM for rest
			dsERROR[ii] = NANUM;
			dsSURLLIM[ii] = NANUM;
			dsSURULIM[ii] = NANUM;
		}
		else                            // Else time value not censored (is a failure time)... 
		{
			if( jj == nd || vTIME[ii - 1] != dPrevUncensoredTime ) // If new failure time get index of failure time in vTP (and index of associated  vectors vP and vPSIG)
				jj--;                                              //  vTP only contains distinct failure times (not censored or duplicate
			                                                       //  times) so skip decrementing index when times are censored or repeated 
			dPrevUncensoredTime = vTIME[ii - 1];                   // Set previous uncensored time 
			dsTIME[ii] = vTIME[ii - 1];                            // Output time
			dsSURVIVAL[ii] = vP[jj];                               // Output Survivorship Function value
			dsERROR[ii] = vPSIG[jj];                               // Output Survivorship Function Error (standard deviation) value
			dsSURLLIM[ii] = dsSURVIVAL[ii] - dZ * dsERROR[ii];     // Output Lower Limit of Survivorship Function value
			if( dsSURLLIM[ii] < 0 )                                // Lower Limit has min of 0
				dsSURLLIM[ii] = 0;
			dsSURULIM[ii] = dsSURVIVAL[ii] + dZ * dsERROR[ii];     // Output Upper Limit of Survivorship Function value
			if( dsSURULIM[ii] > 1 )                                // Upper Limit has max of 1
				dsSURULIM[ii] = 1;
		}
	}
	dsTIME[ii] = 0.0;                                              // Output values for time 0 (not computed...just known)
	dsSURVIVAL[ii] = 1.0;
	dsERROR[ii] = 0.0;
	dsSURLLIM[ii] = 1.0;
	dsSURULIM[ii] = 1.0;

	// *** Output Quartile Estimates and Confidence Limits to output datasets in Survival worksheet *** 
	for( jj = 0; jj < 3; jj++ )
	{
		dsPERCENT[jj] = ( jj + 1 ) * 0.25;        // Compute and output decimal percent

		if( ardPercentile[jj] < 0 )               // Output Percentile
			dsESTIMATE[jj] = NANUM;
		else
			dsESTIMATE[jj] = ardPercentile[jj];

		if( ardPerLLim[jj] < 0 )                  // Output Lower Limit of Percentile
			dsPERLLIM[jj] = NANUM;
		else
			dsPERLLIM[jj] = ardPerLLim[jj];

		if( ardPerULim[jj] < 0 )                  // Output Upper Limit of Percentile
			dsPERULIM = NANUM;
		else
			dsPERULIM[jj] = ardPerULim[jj];				
	}

	return SA_NO_ERROR;
}

////////////////////////////////////////////////////////////////////////////////////
/**
		Implements the Cox Proportional Hazards Model. The computational engine is
		the NAG function nag_surviv_cox_model (g12bac). Additional reference is
		Chapter 3 of Applied Survival Analysis: Regression Modeling of Time to Event
		Data (1998), Hosmer, D.W. and Lemeshow, S. (1998) John Wiley and Sons Inc.,
		New York, NY.
	Example:
		See [CoxPHM] section of CoxPHM.OGS for sample call.
	Parameters:
		nPts=Number of data points in Time, Censor, and Covariate datasets
		nCovars=Number of Covariate datasets
		nResults=Bitwise integer specifying output options
			1=Output Survivorship Function in Resuls Log
			2=Send output to temporary worksheet in addition to Results Log
			4=Create plot of Survivorship Function (implemented in CoxPHM.OGS)
		strSurvivalWKS=Name of temporary worksheet holding sorted copies of all
			required input datasets (Time, Censor, and Covariate) and output datasets   
	Return:
		Depending on nResults, outputs a Summary of Event and Censored values, a table of
		Covariate Parameter Estimates, and the Survivorship Function to the Results Log and
		to a temporary worksheet. Returns SA_NO_ERROR on successful exit or an SA_ error
		code.
*/
int saCoxPHM( int nPts, int nCovars, int nResults, string strSurvivalWKS )
{
	int iSize, iErr, ii, jj, nd;
	string strDB, strCovar, strDatasetName;
	char szTemp[40];
	double ddev;
	BOOL bErr;

	// Attach to text dataset in temporary worksheet that will hold names of Covariate datasets 
	strDatasetName = strSurvivalWKS + COX_COVARIATE_DATASET_NAME;
	Dataset dsCovarNAMES( strDatasetName );
	dsCovarNAMES.SetSize( nCovars );

	// Create dataset and matrix objects to hold covariate data
	Dataset dsCOVAR;		// Holds one covariate dataset at a time
	matrix<double> mZ;		// Holds all covariate datasets
	mZ.SetSize( nPts, nCovars );

	// vSZ[jj]==1 tells NAG that jjth Covariate is to be used in model...all are used so all == 1
	vector<int> vSZ;
	vSZ.SetSize( nCovars );
	vSZ = 1;

	// vB[jj] is jjth parameter estimate on entry to NAG and jjth parameter value on exit from NAG
	vector<double> vB;
	vB.SetSize( nCovars );
	vB = 0;

	// vSE[jj] is asymptotic standard error of jjth parameter estimate on exit from NAG 
	vector<double> vSE;
	vSE.SetSize( nCovars );
	vSE = 0;
	
	// Loop on all Covariate datasets, write actual name of raw Covariate datasets out to text dataset
	//  and read temporary sorted Covariate datasets into Covariate matrix mZ  
	for( jj = 0; jj < nCovars; jj++ )
	{
		// Read actual name of jjth raw Covariate dataset from Dialog Builder listbox 
		strDB.Format( "%%A=CoxPHM!CovariatesLBX.v%d$", jj + 1 );
		LT_execute( strDB );
		LT_get_str( "%A", szTemp, 40 );
		strCovar = szTemp;
		
		// Write out actual name of jjth raw Covariate dataset to temporary worksheet
		bErr = dsCovarNAMES.SetText( jj, strCovar );
		if( bErr == FALSE )
			return SA_SET_TEXT_ERROR;
		
		// Constuct name of the jjth temporary sorted Covariate dataset to use in computation
		strCovar.Format( COX_SORT_COVAR_DATA_NAME, strSurvivalWKS, jj + 1 );
		
		// Attach to the jjth temporary sorted Covariate dataset and copy into jjth column in Covariate matrix
		dsCOVAR.Attach( strCovar );
		for( ii = 0; ii < nPts; ii++ )
		{
			mZ[ii][jj] = dsCOVAR[ii];
		}
		
		// Finished with jjth Covariate
		dsCOVAR.Detach();
	}
		
	// Get Time dataset and store in vector
	strDatasetName = strSurvivalWKS + SA_TIME_DATASET_NAME;
	Dataset dsTIME( strDatasetName );
	iSize = dsTIME.GetSize();
	vector<double> vTIME;
	vTIME.SetSize( iSize );
	vTIME = dsTIME;
	dsTIME.Detach(); // No longer need raw Time dataset	
	
	// Get Censor dataset and store in vector
	strDatasetName = strSurvivalWKS + SA_CENSOR_DATASET_NAME;
	Dataset dsCENSOR( strDatasetName );
	iSize = dsCENSOR.GetSize();
	vector<int> vCENSOR;
	vCENSOR.SetSize( iSize );
	vCENSOR = dsCENSOR;

	// Compute basic descriptive statistics on Censor dataset...use Total in Summary of Event and Censored values
	Dataset *pdsCENSOR = &dsCENSOR;
	BasicStats bsStat;
	BasicStats *pbsStat = &bsStat;
	bErr = Data_sum( pdsCENSOR, pbsStat );
	if( bErr == FALSE )
		return SA_DATA_SUM_ERROR2;
	
	dsCENSOR.Detach(); // No longer need raw Censor dataset	

	// Output failure times...not used by CoxPHM but NAG requires
	vector<double> vTP;
	vTP.SetSize( iSize );

	// Output Survivorship function...NAG expects matrix but second dimension is 1 so vector is used
	vector<double> vSUR;
	vSUR.SetSize( iSize );

	// Score Function not used by CoxPHM but NAG requires
	vector<double> vSC;
	vSC.SetSize( nCovars );

	// Covariate matrix not used by CoxPHM but NAG requires
	vector<double> vCOV;
	vCOV.SetSize( nCovars*( nCovars + 1 )/2 );

	// Residuals not used by CoxPHM but NAG requires
	vector<double> vRES;
	vRES.SetSize( iSize );
	
	// From <NAG\OCN_g12.h>: g12bac nag_surviv_cox_model: Compute estimates of parameters of Cox model
	//  and Survivorship function values. CoxPHM does not use strata (ns=0), Offsets and Stratum are not used
	//  and ==NULL, second dimension of vSUR is 1 (tdsur==1), do not print (iprint=0), and do not use Outfile
	//  which also ==NULL. 
	iErr = nag_surviv_cox_model( nPts, nCovars, 0, mZ, nCovars, vSZ, nCovars, vTIME, vCENSOR, NULL, NULL, &ddev,
		 vB, vSE, vSC, vCOV, vRES, &nd, vTP, vSUR, 1, nPts, COX_MODEL_TOLERANCE, COX_MODEL_ITERATIONS, 0, NULL );
	if( iErr != NE_NOERROR )
		return SA_SURVIV_COX_ERROR;

	// Vector to hold output Chi-Square values
	vector<double> vCHISQ;
	vCHISQ.SetSize( nCovars );
	
	// Vector to hold output P values
	vector<double> vPVAL;
	vPVAL.SetSize( nCovars );
	
	// Vector to hold output Hazard values
	vector<double> vHAZARD;
	vHAZARD.SetSize( nCovars );
	
	// Compute Chi-Square, P(Chi-Square), and Hazard values
	iErr = saCoxPHM_Compute_ChiSq_Prob_Hazard( nCovars, vB, vSE, vCHISQ, vPVAL, vHAZARD );
	if( iErr != SA_NO_ERROR )
		return iErr;

	// Output results in Results Log
	iErr = saCoxPHM_Output_to_ResultsLog( nResults, strSurvivalWKS, nPts, nCovars, vTIME, vCENSOR, pbsStat->total, ddev, 
		 vSUR, vB, vSE, vCHISQ, vPVAL, vHAZARD );
	if( iErr != SA_NO_ERROR )
		return iErr;

	// Output results in worksheet? If user wants worksheet or plot write everything out to worksheet 
	if( nResults > 1 )
	{
		// Function to send output to worksheet
		iErr = saCoxPHM_Output_to_Worksheet ( strSurvivalWKS, nPts, nCovars, vTIME, vCENSOR, nd, vSUR, vB, vSE, vCHISQ,
			 vPVAL, vHAZARD );
		if( iErr != SA_NO_ERROR )
			return iErr;
	}
	
	return SA_NO_ERROR;
}

/**
		Computes the Chi-Square, P(Chi-Square) and Hazard values
	Example:
		See function saCoxPHM above for sample call.
	Parameters:
		nCovars=Number of Covariate datasets
		vB=Vector holding Covariate parameter estimates
		vSE=Vector holding standard errors of Covariate parameter estimates
		vCHISQ=Output vector holding Chi-Square values 
		vPVAL=Output vector holding P values
		vHAZARD=Output vector holding Hazard values
	Return:
		Returns Chi-Square, P(Chi-Square), and Hazard values and SA_NO_ERROR on successful exit
		or an SA_ error code.
*/
int saCoxPHM_Compute_ChiSq_Prob_Hazard( int nCovars, vector<double> &vB, vector<double> &vSE, vector<double> &vCHISQ,
	 vector<double> &vPVAL, vector<double> &vHAZARD )
{
	int jj;
	double dChiSq;
	NagError neErr;
	NagError *pneFail = &neErr;		// NAG error structure

	// Loop on all Covariate datasets
	for( jj = 0; jj < nCovars; jj++ )
	{
		vCHISQ[jj] = pow( vB[jj] / vSE[jj], 2 ); // Compute Chi-Square values
		// From <NAG\OCN_g01.h>: g01ecc nag_prob_chi_sq: Compute P Values	
		vPVAL[jj] = 1 - nag_prob_chi_sq( Nag_LowerTail, vCHISQ[jj], 1, pneFail ); // Compute P(Chi-Square) values
		if( pneFail->code != NE_NOERROR )
			return SA_CHI_SQ_ERROR2;
		vHAZARD[jj] = exp( vB[jj] );  // Compute Hazard values
	}

	return SA_NO_ERROR;
}

/**
		Outputs the results of the Cox Proportional Hazards Model to the Results Log.
	Example:
		See the function saCoxPHM above for sample call.
	Parameters:
		nResults=Bitwise flag (bit 1) indicating whether to output the Survivorship function
			to Results Log
		strSurvivalWKS=Name of current Survival Analysis worksheet
		nPts=Number of data points in Time, Censor, and Covariate datasets 
		nCovars=Number of Covariate datasets
		vTIME=Vector of Time values
		vCENSOR=Vector of Censor values
		dCensoredPts=Number of censored data points
		ddev=-2 ln L
		vSUR=Vector containing Survivorship function values
		vB=Vector containing Covariate Parameter Estimates
		vSE=Vector containing Standard Errors
		vCHISQ=Vector containing Chi-Square Statistics
		vPVAL=Vector containing P Values
		vHAZARD=Vector containing Hazard Ratios
	Return:
		Outputs the results of the Cox Proportional Hazards Model (including a Summary of Event and Censored
		values, a Covariate Parameters Estimate table, and a Survivorship Function table) to the Results
		Log and returns SA_NO_ERROR on successful exit or an SA_ error code.
*/
int saCoxPHM_Output_to_ResultsLog( int nResults, string strSurvivalWKS, int nPts, int nCovars, vector<double> &vTIME,
	 vector<int> &vCENSOR, double dCensoredPts, double ddev, vector<double> &vSUR, vector<double> &vB,
	 vector<double> &vSE, vector<double> &vCHISQ, vector<double> &vPVAL, vector<double> &vHAZARD )
{
	int ii, jj, iErr;
	double dPrevUncensoredTime;
	string str;
	string strDB;
	string strTimeData;
	string strCensorData;
	string strCensorValue;
	string strCovar;
	string strDatasetName;
	char szTemp[40];
	BOOL bErr;

	int iRow, iCol;
	string strTablePath;
	string strTableName;
	Worksheet wksParamTable;
	
	// *** Output Survival Analysis Main Headers and Summary Table ***
	// Get name of Time dataset from dialog box
	strDB = "%A=CoxPHM!TimeDataEBX$";
	LT_execute( strDB );
	LT_get_str( "%A", szTemp, 40 );
	strTimeData = szTemp;

	// Get name of Censor dataset from dialog box
	strDB = "%A=CoxPHM!CensorDataEBX$";
	LT_execute( strDB );
	LT_get_str( "%A", szTemp, 40 );
	strCensorData = szTemp;
	
	// Get Censor Value from dialog box
	strDB = "%A=CoxPHM!CensorValueEBX$";
	LT_execute( strDB );	
	LT_get_str( "%A", szTemp, 40 );
	strCensorValue = szTemp;
	
	// Get name of dialog box (or feature) from Local.h
	strDB.Format( COX_MAIN_HDR );

	// Output survival analysis main headers and summary table for Cox Proportional Hazards Model
	iErr = saSurvival_Analysis_ResultsLog_Summary( strDB, strTimeData, strCensorData, strCensorValue,
		 nPts, dCensoredPts );
	if( iErr != SA_NO_ERROR )
		return iErr;

	// *** Output Cox Parameter Estimates Table ***
	// Attach to dataset holding names of covariate datasets
	strDatasetName = strSurvivalWKS + COX_COVARIATE_DATASET_NAME;
	Dataset dsCovarNAMES( strDatasetName );
	dsCovarNAMES.SetSize( nCovars );
	
	// Create and attach to a temporary output worksheet 
	strTablePath = GetFullPath( "saCoxParamEst.OGW", SA_OGW_SUBFOLDER, TRUE );
	if( !strTablePath.IsFile() )
	{
		Type_ErrorMsg1( SA_FILE_NOT_FOUND_ERROR_MSG, strTablePath );
		return SA_ABORT_NO_ERROR_MSG;
	}
	wksParamTable.Open( strTablePath, CREATE_TEMP );
	strTableName = wksParamTable.GetPage().GetName();
	
	// Rows and columns are indexed from 0
	
	// Insert a blank line in parameter estimates table for each covariate starting at row 6
	iRow = 5;
	Type_Insert_Blank_Lines_in_Wks( wksParamTable, iRow, nCovars );
 		
	// Loop on all covariate datasets counting covariates with jj and current row with iRow
	for( jj = 0; jj < nCovars; iRow++, jj++ )
	{
		// Output name of covariate dataset in column 2
		bErr = dsCovarNAMES.GetText( jj, strCovar );
		if( bErr == FALSE )
			return SA_GET_TEXT_ERROR;
		iCol = 1;
		wksParamTable.SetCell( iRow, iCol, strCovar );
		
		// Output parameter estimate in column 3
		iCol++;
		wksParamTable.SetCell( iRow, iCol, vB[jj] );
		
		// Output standard error in column 4		
		iCol++;
		wksParamTable.SetCell( iRow, iCol, vSE[jj] );
		
		// Output chi-square statistic in column 5		
		iCol++;
		wksParamTable.SetCell( iRow, iCol, vCHISQ[jj] );
		
		// Output p value in column 6
		iCol++;
		wksParamTable.SetCell( iRow, iCol, vPVAL[jj] );
		
		// Output hazard ratio in column 7
		iCol++;
		wksParamTable.SetCell( iRow, iCol, vHAZARD[jj] );
	}
	
	// Output -2 ln L after separater row and in column 3
	iRow++;
	iCol = 2;
	wksParamTable.SetCell( iRow, iCol, ddev );

	// GJL 11/14/02 v7.0434 USE_OC_WKS_TYPE_NOT_LT
	// Send output parameter estimates worksheet to Results Log
	wksParamTable.Type( TYPETARGET_OUTPUTLOG, NULL, 8 );
	
	// *** Output Survivorship Function ***	
	if( nResults & 1 )
	{
		Worksheet wksSurvFuncTable;

		// Create and attach to a temporary worksheet 
		strTablePath = GetFullPath( "saCoxSurvFunc.OGW", SA_OGW_SUBFOLDER, TRUE );
		if( !strTablePath.IsFile() )
		{
			Type_ErrorMsg1( SA_FILE_NOT_FOUND_ERROR_MSG, strTablePath );
			return SA_ABORT_NO_ERROR_MSG;
		}
		wksSurvFuncTable.Open( strTablePath, CREATE_TEMP );
		strTableName = wksSurvFuncTable.GetPage().GetName();
		
		// Rows and columns are indexed from 0
		
		// Insert a blank line in survivorship function table for each data point + 1
		// (+ 1 for time 0) starting at row 5
		iRow = 4;
		Type_Insert_Blank_Lines_in_Wks( wksSurvFuncTable, iRow, nPts + 1 );
		
		// Output time value = 0 in column 2
		iCol = 1;
		wksSurvFuncTable.SetCell( iRow, iCol, 0.0 );
			
		// Output survivorship function value = 1 for time = 0 in column 3
		iCol++;
		wksSurvFuncTable.SetCell( iRow, iCol, 1.0 );

		// Output non-zero time values starting in row 5
		iRow++;
		
		// Loop on all non-zero time values
		jj = -1;
		dPrevUncensoredTime = -1;
		for( ii = 0;ii < nPts; iRow++, ii++ )
		{
			// Output time value in column 2
			iCol = 1;
			wksSurvFuncTable.SetCell( iRow, iCol, vTIME[ii] );
			
			// Output survivorship function value in column 3
			iCol++;
			if( vCENSOR[ii] == 1 ) // If censored...
			{
				// Output missing value for survivorship function value
				wksSurvFuncTable.SetCell( iRow, iCol, NANUM );
			}
			else // Else if not censored...
			{
				// Output survivorship function value 
				// Must index into vSUR differently than vTIME because only non-repeated and uncensored survival
				//  values are in vSUR while all times are in vTIME 
				if( jj < 0 || vTIME[ii] != dPrevUncensoredTime )
					jj++;
				dPrevUncensoredTime = vTIME[ii];
				wksSurvFuncTable.SetCell( iRow, iCol, vSUR[jj] );
			}
		}			

	// GJL 11/14/02 v7.0434 USE_OC_WKS_TYPE_NOT_LT
	// Send output survivorship function worksheet to Results Log
	wksSurvFuncTable.Type( TYPETARGET_OUTPUTLOG, NULL, 4 );

	}

	return SA_NO_ERROR;
}

/**
		Outputs the results of the Cox Proportional Hazards Model to a worksheet.
	Example:
		See the function saCoxPHM above for sample call.
	Parameters:
		strSurvivalWKS=Name of current Survival Analysis worksheet
		nPts=Number of data points in Time, Censor, and Covariate datasets 
		nCovars=Number of Covariate datasets
		vTIME=Vector of Time values
		vCENSOR=Vector of Censor values
		nd=Number of distinct failure times
		vSUR=Vector containing Survivorship function values
		vB=Vector containing Covariate Parameter Estimates
		vSE=Vector containing Standard Errors
		vCHISQ=Vector containing Chi-Square Statistics
		vPVAL=Vector containing P Values
		vHAZARD=Vector containing Hazard Ratios
	Return:
		Outputs the results of the Cox Proportional Hazards Model (including a Covariate Parameters Estimate
		table and a Survivorship Function table) to a worksheet. Returns SA_NO_ERROR on successful exit.
*/
int saCoxPHM_Output_to_Worksheet( string strSurvivalWKS, int nPts, int nCovars, vector<double> &vTIME,
	 vector<int> &vCENSOR, int nd, vector<double> &vSUR, vector<double> &vB, vector<double> &vSE,
	 vector<double> &vCHISQ, vector<double> &vPVAL, vector<double> &vHAZARD )
{
	int ii, jj;
	double dPrevUncensoredTime;
	string strDatasetName;

	// Output Time dataset
	strDatasetName = strSurvivalWKS + SA_TIME_DATASET_NAME;
	Dataset dsTIME( strDatasetName );
	dsTIME.SetSize( nPts + 1 );

	// Output Survivorship Function dataset
	strDatasetName = strSurvivalWKS + SA_SURVIVAL_DATASET_NAME;
	Dataset dsSURVIVAL( strDatasetName );
	dsSURVIVAL.SetSize( nPts + 1 );

	// *** Output Survivorship Function to Worksheet ***
	dPrevUncensoredTime = -1;
	jj = nd;
	// Loop backward on all points in Time dataset
	for( ii = nPts; ii > 0; ii-- )
	{
		if( vCENSOR[ii-1] == 1 )           // If censored time...
		{
			dsTIME[ii] = vTIME[ii - 1];    // Output time value and NANUM for survival function value
			dsSURVIVAL[ii] = NANUM;
		}
		else                               // Else an event (failure) time...
		{
			if( jj == nd || vTIME[ii - 1] != dPrevUncensoredTime ) // If first non-censored time or non-repeated time... 
				jj--;                      // Decrement index into vSUR which only holds survival function values for non-censored and non-repeated times  
			dPrevUncensoredTime = vTIME[ii - 1]; // Remember previous non-censored time 
			dsTIME[ii] = vTIME[ii - 1];    // Output time value
			dsSURVIVAL[ii] = vSUR[jj];     // Output Survivorship Function value
			
		}
	}
	// Output survivorship function value for time==0; not computed but are known 
	dsTIME[ii] = 0.0;
	dsSURVIVAL[ii] = 1.0;

	// Output Covariate Parameter Estimate dataset name	
	strDatasetName = strSurvivalWKS + COX_PARAMEST_DATASET_NAME;
	Dataset dsPARAM( strDatasetName );
	dsPARAM.SetSize( nCovars );

	// Output Covariate Parameter Estimate Standard Error dataset name	
	strDatasetName = strSurvivalWKS + SA_STDERROR_DATASET_NAME;		
	Dataset dsERROR( strDatasetName );
	dsERROR.SetSize( nCovars );

	// Output Covariate Parameter Estimate Chi-Square dataset name	
	strDatasetName = strSurvivalWKS + COX_CHISQ_DATASET_NAME;		
	Dataset dsCHISQ( strDatasetName );
	dsCHISQ.SetSize( nCovars );

	// Output Covariate Parameter Estimate P Value dataset name	
	strDatasetName = strSurvivalWKS + COX_PVALUE_DATASET_NAME;		
	Dataset dsPVAL( strDatasetName );
	dsPVAL.SetSize( nCovars );

	// Output Covariate Parameter Estimate Hazard Ratio dataset name	
	strDatasetName = strSurvivalWKS + COX_HAZARD_DATASET_NAME;		
	Dataset dsHAZARD( strDatasetName );
	dsHAZARD.SetSize( nCovars );
		
	// Copy estimates of parameters et al for Cox model to output datasets in Survival worksheet  
	for( ii = 0; ii < nCovars; ii++ )
	{
		dsPARAM[ii] = vB[ii];
		dsERROR[ii] = vSE[ii];
		dsCHISQ[ii] = vCHISQ[ii];
		dsPVAL[ii] = vPVAL[ii];
		dsHAZARD[ii] = vHAZARD[ii];
	}
	
	return SA_NO_ERROR;
}

/**
		Outputs to the Results Log the survival analysis main headers and a summary of events
		and censored values for both the Kaplan-Meier Estimator and the Cox Proportional Hazards
		Model. 
	Example:
		See functions saKaplanMeier_Output_to_ResultsLog and saCoxPHM_Output_to_ResultsLog above
		for sample call.
	Parameters:
		strDB=Dialog Box or feature name
		strTimeData=Name of Time dataset
		strCensorData=Name of Censor dataset
		strCensorValue=Censor value
		nPts=Number of data points in Time and Censor datasets
		dCensoredPts=Number of data points that are censored 
	Return:
		Outputs the survival analysis main headers and a summary table and returns SA_NO_ERROR
		on successful exit.
*/
static int saSurvival_Analysis_ResultsLog_Summary( string strDB, string strTimeData, string strCensorData, string strCensorValue,
	 int nPts, double dCensoredPts )
{
	int iRow, iCol;
	int iEvents;
	int iCensored;
	double dnPts;
	double dPercentCensored;
	string strTablePath;
	string strOut;
	string strTableName;
	Worksheet wksSumTable;

	// Compute summary values for output
	iCensored = (int) dCensoredPts;
	iEvents = nPts - iCensored;
	dnPts = (double) nPts;
	dPercentCensored = 10000.0 * dCensoredPts / dnPts + 0.5; // Round to 100ths place
	dPercentCensored = (int) dPercentCensored;
	dPercentCensored = dPercentCensored / 100;

	// Create and attach to a temporary worksheet 
	strTablePath = GetFullPath( "saSAEventSum.OGW", SA_OGW_SUBFOLDER, TRUE );
	if( !strTablePath.IsFile() )
	{
		Type_ErrorMsg1( SA_FILE_NOT_FOUND_ERROR_MSG, strTablePath );
		return SA_ABORT_NO_ERROR_MSG;
	}
	wksSumTable.Open( strTablePath, CREATE_TEMP );
	strTableName = wksSumTable.GetPage().GetName();
	
	// Rows and columns are indexed from 0

	// *** Output Survival Analysis Main Headers (starting in row 1 and column 1) ***
	// Output the dialog box name (feature name)
	iRow = 0;
	iCol = 0;
	wksSumTable.SetCell( iRow, iCol, strDB );
	
	// Output the name of the Time dataset
	iRow = 2;
	iCol = 3;
	wksSumTable.SetCell( iRow, iCol, strTimeData );
	
	// Output the name of the Censor dataset
	iRow++;
	wksSumTable.SetCell( iRow, iCol, strCensorData );
	
	// Output the Censor Value
	iRow++;
	wksSumTable.SetCell( iRow, iCol, strCensorValue );

	// *** Output Summary Data (in row 12 and starting in column 2) ***
	// Output the number of data points in the Time and Censor datasets
	iRow = 11;
	iCol = 1;
	wksSumTable.SetCell( iRow, iCol, nPts );
	
	// Output the number of events in column 3
	iCol++;
	wksSumTable.SetCell( iRow, iCol, iEvents );
	
	// Output the number of censored data points in column 4
	iCol++;
	wksSumTable.SetCell( iRow, iCol, iCensored );
	
	// Output the percent of censored events in column 5
	iCol++;
	strOut = LocalizeDouble( dPercentCensored, SA_PERCENT_FORMAT );
	wksSumTable.SetCell( iRow, iCol, strOut );
	
	// GJL 11/14/02 v7.0434 USE_OC_WKS_TYPE_NOT_LT
	// Send output worksheet to Results Log
	wksSumTable.Type( TYPETARGET_OUTPUTLOG, NULL, 6 );

	return SA_NO_ERROR;
}