https://github.com/cran/BDgraph
Tip revision: 5de6bccb9fb5a330b645447ff3957bb280817926 authored by Reza Mohammadi on 11 May 2018, 10:38:18 UTC
version 2.49
version 2.49
Tip revision: 5de6bcc
compare.R
## ------------------------------------------------------------------------------------------------|
# Copyright (C) 2012 - 2018 Reza Mohammadi |
# |
# This file is part of BDgraph package. |
# |
# BDgraph is free software: you can redistribute it and/or modify it under |
# the terms of the GNU General Public License as published by the Free |
# Software Foundation; see <https://cran.r-project.org/web/licenses/GPL-3>. |
# |
# Maintainer: Reza Mohammadi <a.mohammadi@uva.nl> |
## ------------------------------------------------------------------------------------------------|
# This function reports below measures to assess the performance of estimated graphs: |
## ------------------------------------------------------------------------------------------------|
# True positive: number of correctly estimated links. |
# True negative: number of true non-existing links which is correctly estimated. |
# False positive: number of links which are not in the true graph, but are incorrectly estimated. |
# False negative: number of links which they are in the true graph, but are not estimated. |
# F1-score: weighted average of the positive predictive and true positive rate. |
# Specificity: Specificity value reaches. |
# Sensitivity: Sensitivity value reaches. |
# MCC: Matthews Correlation Coefficients (MCC). |
## ------------------------------------------------------------------------------------------------|
compare = function( sim.obj, bdgraph.obj, bdgraph.obj2 = NULL, bdgraph.obj3 = NULL,
colnames = NULL, vis = FALSE )
{
if( is.matrix( sim.obj ) )
{
if( ( sum( sim.obj == 0 ) + sum( sim.obj == 1 ) ) != ( p ^ 2 ) ) stop( "Element of 'sim.obj' must be 0 or 1" )
G = sim.obj
}
if( is.matrix( bdgraph.obj ) )
{
if( ( sum( bdgraph.obj == 0 ) + sum( bdgraph.obj == 1 ) ) != ( p ^ 2 ) ) stop( "Element of 'bdgraph.obj' must be 0 or 1" )
est = bdgraph.obj
}
if( !is.matrix( sim.obj ) && ( class( sim.obj ) == "sim" ) ) G <- sim.obj $ G
if( !is.matrix( bdgraph.obj ) && ( class( bdgraph.obj ) == "bdgraph" ) ) est <- BDgraph::select( bdgraph.obj )
if( !is.matrix( bdgraph.obj ) && ( class( bdgraph.obj ) == "ssgraph" ) ) est <- BDgraph::select( bdgraph.obj )
if( class( bdgraph.obj ) == "select" ) est <- bdgraph.obj $ refit
G = as.matrix( G ) # G is the adjacency matrix of true graph
est = as.matrix( est ) # est is the adjacency matrix of estimated graph
p = nrow( G )
if( sum( dim( G ) == dim( est ) ) != 2 ) stop( "'sim.obj' and 'bdgraph.obj' have non-conforming size" )
G[ lower.tri( G, diag = TRUE ) ] = 0
est[ lower.tri( est, diag = TRUE ) ] = 0
result = matrix( 1, 8, 2 )
result[ , 2 ] = compute_measures( G = G, est_G = est )
if( !is.null( bdgraph.obj2 ) )
{
if( is.matrix( bdgraph.obj2 ) )
{
if( ( sum( bdgraph.obj2 == 0 ) + sum( bdgraph.obj2 == 1 ) ) != ( p ^ 2 ) ) stop( "Element of 'bdgraph.obj2' must be 0 or 1" )
est2 = bdgraph.obj2
}
if( !is.matrix( bdgraph.obj2 ) && ( class( bdgraph.obj2 ) == "bdgraph" ) ) est2 <- BDgraph::select( bdgraph.obj2 )
if( !is.matrix( bdgraph.obj2 ) && ( class( bdgraph.obj2 ) == "ssgraph" ) ) est2 <- BDgraph::select( bdgraph.obj2 )
if( class( bdgraph.obj2 ) == "select" ) est2 <- bdgraph.obj2 $ refit
est2 = as.matrix( est2 )
if( sum( dim( G ) == dim( est2 ) ) != 2 ) stop( "'sim.obj' and 'bdgraph.obj2' have non-conforming size" )
est2[ lower.tri( est2, diag = TRUE ) ] = 0
result = cbind( result, compute_measures( G = G, est_G = est2 ) )
}
if( !is.null( bdgraph.obj3 ) )
{
if( is.matrix( bdgraph.obj3 ) )
{
if( ( sum( bdgraph.obj3 == 0 ) + sum( bdgraph.obj3 == 1 ) ) != ( p ^ 2 ) ) stop( "Element of 'bdgraph.obj3' must be 0 or 1" )
est3 = bdgraph.obj3
}
if( !is.matrix( bdgraph.obj3 ) && ( class( bdgraph.obj3 ) == "bdgraph" ) ) est3 <- BDgraph::select( bdgraph.obj3 )
if( !is.matrix( bdgraph.obj3 ) && ( class( bdgraph.obj3 ) == "ssgraph" ) ) est3 <- BDgraph::select( bdgraph.obj3 )
if( class( bdgraph.obj3 ) == "select" ) est3 <- bdgraph.obj3 $ refit
est3 = as.matrix( est3 )
if( sum( dim( G ) == dim( est3 ) ) != 2 ) stop( "'sim.obj' and 'bdgraph.obj3' have non-conforming size" )
est3[ lower.tri( est3, diag = TRUE ) ] = 0
result = cbind( result, compute_measures( G = G, est_G = est3 ) )
}
result[ c( 3, 4 ), 1 ] = 0
result[ 1, 1 ] = sum( G )
result[ 2, 1 ] = p * ( p - 1 ) / 2 - result[ 1, 1 ]
result[ is.na( result ) ] = 0
if( is.null( colnames ) )
{
colnames = c( "True", "estimate" )
if( !is.null( bdgraph.obj2 ) ) colnames = c( colnames, "estimate2" )
if( !is.null( bdgraph.obj3 ) ) colnames = c( colnames, "estimate3" )
}
colnames( result ) <- colnames
rownames( result ) <- c( "true positive", "true negative", "false positive", "false negative",
"F1-score", "specificity", "sensitivity", "MCC" )
if( vis )
{
G_igraph <- graph.adjacency( G, mode = "undirected", diag = FALSE )
est_igraph <- graph.adjacency( est, mode = "undirected", diag = FALSE )
if ( p < 20 ) sizev = 15 else sizev = 2
row_plot = ifelse( is.null( bdgraph.obj2 ), 1, 2 )
op = par( mfrow = c( row_plot, 2 ), pty = "s", omi = c( 0.3, 0.3, 0.3, 0.3 ), mai = c( 0.3, 0.3, 0.3, 0.3 ) )
plot.igraph( G_igraph, layout = layout.circle, main = colnames[1], vertex.color = "white", vertex.size = sizev, vertex.label.color = 'black' )
plot.igraph( est_igraph, layout = layout.circle, main = colnames[2], vertex.color = "white", vertex.size = sizev, vertex.label.color = 'black' )
if( !is.null( bdgraph.obj2 ) )
{
est2_igraph <- graph.adjacency( as.matrix(est2), mode = "undirected", diag = FALSE )
plot.igraph( est2_igraph, layout = layout.circle, main = colnames[3], vertex.color = "white", vertex.size = sizev, vertex.label.color = 'black' )
}
if( !is.null( bdgraph.obj3 ) )
{
est3_igraph <- graph.adjacency( as.matrix(est3), mode = "undirected", diag = FALSE )
plot.igraph( est3_igraph, layout = layout.circle, main = colnames[4], vertex.color = "white", vertex.size = sizev, vertex.label.color = 'black' )
}
par( op )
}
return( round( result, 3 ) )
}
## ------------------------------------------------------------------------------------------------|
# To compare measures the performance of estimated graphs based on true graph
## ------------------------------------------------------------------------------------------------|
compute_measures = function( G, est_G )
{
upper_G = G[ upper.tri( G ) ]
upper_est_G = est_G[ upper.tri( est_G ) ]
tp = sum( ( upper_G != 0 ) * ( upper_est_G != 0 ) )
tn = sum( ( upper_G == 0 ) * ( upper_est_G == 0 ) )
fp = sum( ( upper_G == 0 ) * ( upper_est_G != 0 ) )
fn = sum( ( upper_G != 0 ) * ( upper_est_G == 0 ) )
# harmonic mean of precision and recall, called F-measure or balanced F-score
F1score = ( 2 * tp ) / ( 2 * tp + fp + fn )
specificity = tn / ( tn + fp )
sensitivity = tp / ( tp + fn )
# Matthews Correlation Coefficients (MCC)
mcc = ( ( tp * tn ) - ( fp * fn ) ) / ( sqrt( ( tp + fp ) * ( tp + fn ) ) * sqrt( ( tn + fp ) * ( tn + fn ) ) )
return( c( tp, tn, fp, fn, F1score, specificity, sensitivity, mcc ) )
}