survivalanalysis.c

图像处理的压缩算法

		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 );