1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110
|
/*=========================================================================
Program: Visualization Toolkit
Module: vtkContingencyStatisticsGnuR.cxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
/*-------------------------------------------------------------------------
Copyright 2011 Sandia Corporation.
Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
the U.S. Government retains certain rights in this software.
-------------------------------------------------------------------------*/
#include "vtkContingencyStatisticsGnuR.h"
#include "vtkRInterface.h"
#include "vtkToolkits.h"
#include "vtkStatisticsAlgorithmPrivate.h"
#include "vtkDoubleArray.h"
#include "vtkIdTypeArray.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkObjectFactory.h"
#include "vtkStringArray.h"
#include "vtkStdString.h"
#include "vtkTable.h"
#include "vtkVariantArray.h"
#include <map>
#include <vector>
#include <sstream>
vtkStandardNewMacro(vtkContingencyStatisticsGnuR);
vtkContingencyStatisticsGnuR::vtkContingencyStatisticsGnuR()
{
}
// ----------------------------------------------------------------------
vtkContingencyStatisticsGnuR::~vtkContingencyStatisticsGnuR()
{
}
// ----------------------------------------------------------------------
void vtkContingencyStatisticsGnuR::PrintSelf( ostream &os, vtkIndent indent )
{
this->Superclass::PrintSelf( os, indent );
}
// ----------------------------------------------------------------------
void vtkContingencyStatisticsGnuR::CalculatePValues( vtkTable* outTab)
{
vtkIdTypeArray* dimCol = vtkArrayDownCast<vtkIdTypeArray>( outTab->GetColumn(0) );
vtkDoubleArray* chi2Col = vtkArrayDownCast<vtkDoubleArray>( outTab->GetColumn(1));
vtkDoubleArray* chi2yCol = vtkArrayDownCast<vtkDoubleArray>( outTab->GetColumn(2));
// Prepare VTK - R interface
vtkRInterface* ri = vtkRInterface::New();
// Use the calculated DOFs and Chi square statistics as inputs to the Chi square function
ri->AssignVTKDataArrayToRVariable( dimCol, "d" );
ri->AssignVTKDataArrayToRVariable( chi2Col, "chi2" );
ri->AssignVTKDataArrayToRVariable( chi2yCol, "chi2y" );
// Now prepare R script and calculate the p-values (in a single R script evaluation for efficiency)
std::ostringstream rs;
rs << "p<-c();"
<< "py<-c();"
<< "for(i in 1:"
<< dimCol->GetNumberOfTuples()
<< "){"
<< "p<-c(p,1-pchisq(chi2[i],d[i]));"
<< "py<-c(py,1-pchisq(chi2y[i],d[i]))"
<< "}";
ri->EvalRscript( rs.str().c_str() );
// Retrieve the p-values
vtkDoubleArray* testChi2Col = vtkArrayDownCast<vtkDoubleArray>( ri->AssignRVariableToVTKDataArray( "p" ) );
vtkDoubleArray* testChi2yCol = vtkArrayDownCast<vtkDoubleArray>( ri->AssignRVariableToVTKDataArray( "py" ) );
if ( ! testChi2Col || ! testChi2yCol
|| testChi2Col->GetNumberOfTuples() != dimCol->GetNumberOfTuples()
|| testChi2yCol->GetNumberOfTuples() != dimCol->GetNumberOfTuples() )
{
vtkWarningMacro( "Something went wrong with the R calculations. Reported p-values will be invalid." );
this->Superclass::CalculatePValues( outTab );
}
else
{
outTab->AddColumn( testChi2Col );
outTab->AddColumn( testChi2yCol );
}
testChi2Col->SetName( "P" );
testChi2yCol->SetName( "P Yates" );
// Clean up
ri->Delete();
}
|