Jewel or Junk?

Data balancing and outlier detection with DataModeler

Mark Kotanchek

Evolved Analytics LLC

Abstract

Lots of variables. Not much in the way of data. The data variables are correlated and coupled. Truth is unknown. What does it all mean?

This is a common situation for many real-world data sets. . In this talk, Mathematica's CountryData[ ] is used as a surrogate for typical social science data to demonstrate how DataModeler’s unique capabilities can be used to efficiently and effectively attack such data. Although basic SymbolicRegression and data exploration features will be briefly reviewed, the emphasis will be recently added capabilities:

Data Balancing — which records and attributes contain the most information? The SMITS algorithm lets us explore the information content of data sets and identify key data records. This lets us quickly build balanced data sets from real-world data which leads to improved and empirical models as well as faster model development.

Variable Selection — which variables count from a nonlinear perspective? Symbolic regression is attractive because it naturally handles correlated inputs and pull out the driving variables and present them without artificial constraints and assumptions. In other words, the data tells us which inputs and input combinations are important.

Niching — what is the trade-off on number of variables or variable combinations? The ability to explore the trade-off between the number of variables or possible variable combination leads to increased insight into the underlying mechanisms in the modeled system.

Outlier Detection — outliers are either the most interesting or nuisance nuggets of data. However, identifying outliers is difficult in multivariate data sets — especially when the outliers are inliers. Although the data balancing processing can identify some outliers, model-based outlier detection lets us quickly identify the data records which warrant special attention.

In the beginning was the data

We are using CountryData[ ] to illustrate social science data modeling

Block[
{numericHeaders, numericData},
{numericHeaders, numericData} = ImportDataMatrix[
Join[{augmentedDataFields},augmentedRawData],
DeleteNonNumericColumns -> 0.95,
DeleteNonNumericRows-> False
 ];
numericData =Map[If[MatchQ[#,_?NumericQ],1,0,0]&,numericData,{2}];

rawPlot = ArrayPlot[Transpose@numericData, 
FrameLabel ->LabelForm@{"Attributes","Countries"},
PlotLabel ->LabelForm["CountryData[ ] Numeric Data Distribution",FontSize -> 18,FontColor ->Red],
AspectRatio ->Automatic,
Mesh -> False, 
ImageSize ->600]
];

Block[
{rejectionThresholds, possibleDataSizes,fractionOfDataSize,countryTrajectory,attributeTrajectory},
rejectionThresholds = Range[0.0,0.95,0.05];

possibleDataSizes =Dimensions@Last[ImportDataMatrix[
Join[{augmentedDataFields},augmentedRawData],
DeleteNonNumericColumns -> (# /. 0.0 ->True),
DeleteNonNumericRows-> True
 ]]&/@rejectionThresholds;
possibleDataSizes = possibleDataSizes /. {0} ->{0,0};

fractionOfDataSize= ((#)/({Length@countries,Length@dataFields}))&/@possibleDataSizes ;

{countryTrajectory,attributeTrajectory}=Transpose[
{rejectionThresholds,#}
]&/@Transpose@fractionOfDataSize;

countryTrajectory= Tooltip[MapThread[Tooltip[#1,#2]&,{countryTrajectory,possibleDataSizes[[All,1]]}],"number of countries"];
attributeTrajectory= Tooltip[MapThread[Tooltip[#1,#2]&,{attributeTrajectory,possibleDataSizes[[All,2]]}],"number of attributes"];

tradeoffPlot = ListPlot[
{countryTrajectory,attributeTrajectory},
Frame->True,
FrameLabel->{"fraction non-numeric allowed initially allowed in columns","fraction of possible countries with complete data"},
PlotLabel -> LabelForm@"Effect of varying DeleteNonNumericColumns threshold",
PlotStyle ->{{Thick,PointSize[Large]}},
GridLines ->Automatic,
Joined -> True,
PlotMarkers->{Automatic,Medium},
ImageSize ->600
]
];

Block[
{numericHeaders, numericData, numberCountriesWithThatAttribute},
{numericHeaders, numericData} = ImportDataMatrix[
Join[{augmentedDataFields},augmentedRawData],
DeleteNonNumericColumns -> 0.5,
DeleteNonNumericRows-> False
 ];

numericData =Map[If[MatchQ[#,_?NumericQ],1,0,0]&,numericData,{2}];
numberCountriesWithThatAttribute = MapThread[Tooltip[#1,#2]&,{Plus@@@Transpose@numericData,numericHeaders}];

attributeCompletenessPlot = ListPlot[Reverse@Sort[numberCountriesWithThatAttribute], 
PlotLabel ->LabelForm["Data Completeness vs. Attribute",FontSize -> 14,FontColor ->Red],
Frame->True,
FrameLabel ->LabelForm@{"sorted by # countries with that attribute","number of countries"},
ImageSize ->600
]
];

Block[
{numericHeaders, numericData,numberOfAttributesPerCountry},
{numericHeaders, numericData} = ImportDataMatrix[
Join[{augmentedDataFields},augmentedRawData],
DeleteNonNumericColumns -> 0.5,
DeleteNonNumericRows-> False
 ];

numericData =Map[If[MatchQ[#,_?NumericQ],1,0,0]&,numericData,{2}];
numberOfAttributesPerCountry=MapThread[Tooltip[#1,#2]&,{Plus@@@numericData,countries}];

countryCompletenessPlot = ListPlot[Reverse@Sort[numberOfAttributesPerCountry], 
PlotLabel ->LabelForm["Data Completeness vs. Country",FontSize -> 14,FontColor ->Red],
Frame->True,
FrameLabel ->LabelForm@{"sorted countries","# attributes"},
ImageSize ->600
]
];

TabView[{
"Raw Data" -> rawPlot,
"Attribute Completeness" -> attributeCompletenessPlot,
"Country Completeness" -> countryCompletenessPlot,
"The Trade-off"-> tradeoffPlot
} ]

Today: Two Views of the data — information content & model-based

The data at the sweet spot — We have a FAT data array

From the trade-off analysis we did while exploring the raw data, it appears that the best bang-for-the-buck in terms of maximizing our data set size is if we set DeleteNonNumericColumns -> 0.15. As we can see, we still have a substantial number of attributes to use in our modeling. Here we have opted to try to model the GDPPerCapita. We also exclude fields containing “GDP”, “ValueAdded” and NationalIncome from the input data set to mitigate the risk of  generating trivial solutions.

Block[
{selectedHeaders, selectedData,fieldsContainingLife,fieldsContainingAge,  fieldsContainingMortality, fieldsContainingDeath, fieldsContainingGDP,fieldsContainingValueAdded,targetAttributePosn, suppressedAttributePosn },

{selectedHeaders, selectedData}=ImportDataMatrix[
augmentedDataFields,augmentedRawData,
DeleteNonNumericColumns ->0.15,
DeleteNonNumericRows-> True
 ];

fieldsContainingGDP =Flatten[StringCases[#,___~~"GDP"~~___]&/@selectedHeaders];
fieldsContainingValueAdded = Flatten[StringCases[#,___~~"ValueAdded"~~___]&/@selectedHeaders];

fieldsContainingLife =Flatten[StringCases[#,___~~"Life"~~___]&/@selectedHeaders];
fieldsContainingAge =Flatten[StringCases[#,___~~"Age"~~___]&/@selectedHeaders];
fieldsContainingMortality =Flatten[StringCases[#,___~~"Mortality"~~___]&/@selectedHeaders];
fieldsContainingDeath =Flatten[StringCases[#,___~~"Death"~~___]&/@selectedHeaders];

(* 
suppressedAttributePosn = Flatten@Position[
selectedHeaders,Alternatives@@Flatten[
{fieldsContainingLife,fieldsContainingAge,fieldsContainingMortality,fieldsContainingDeath}
]
];
targetAttributePosn = First@Flatten@Position[selectedHeaders,"LifeExpectancy"];
*)

suppressedAttributePosn = Flatten@Position[
selectedHeaders,Alternatives@@Flatten[
{fieldsContainingGDP,"NationalIncome",fieldsContainingValueAdded}
]
];
targetAttributePosn = First@Flatten@Position[selectedHeaders,"GDPPerCapita"];


targetedResponse = (Subscript[selectedData, [[All, targetAttributePosn]]]);
targetLabel = (Subscript[selectedHeaders, [[targetAttributePosn]]]);

keepers = Complement[Range@Length@selectedHeaders,{targetAttributePosn},suppressedAttributePosn];
inputLabels = (Subscript[selectedHeaders, [[keepers]]]);
inputMatrix=(Subscript[selectedData, [[All, keepers]]]);

attributeTable = GridTable[
LabelForm[
Partition[Join[inputLabels,{targetLabel}],4,4,{1,1},{}],
FontSize ->9,FontFamily->"Helvetica"],
Label-> LabelForm[{"The selected data consists of",Length@inputLabels+1,"attributes from",Length@inputMatrix,"countries"},Joined -> True,FontSize -> 14, FontColor ->Red]
]
];

correlationPlot =CorrelationMatrixPlot[
MergeInputResponseData[inputMatrix,targetedResponse],
DataVariableLabels-> None,
PlotLabel->LabelForm[{"Correlations of", Length@inputLabels,"input attributes and target response"},Joined->True,
FontSize -> 20]
];

selectedDataPlot = ArrayPlot[
Transpose@RescaleData[
MergeInputResponseData[inputMatrix,targetedResponse],
SelectionFunction->({Min@#,Max@#}&)],
ColorFunction->(Hue[(2/3)#]&),
FrameLabel -> LabelForm[{"Attributes","Countries"},FontSize->Large],
PlotLabel->LabelForm["Which countries are outliers? Which attributes are most important",FontSize->20],
ImageSize ->600
];


TabView[{
"the attributes" -> attributeTable,
"selected data correlation matrix" -> correlationPlot,
"the selected data" -> selectedDataPlot
}]

Data Summary

A FAT data array — nearly as many inputs as data records

Highly correlated inputs

Which variables are important?

Which variables are sufficient?

Which countries are most typical?

Which countries are most atypical?

A linear modeling approach is not viable

Information Content — An Introduction to Data Balancing

The default InformationMetric for BalanceData is NeighborhoodAsymmetry using a Neighborhood of the five nearest data points.  Below we show the information content of a synthetic data set which has been balanced (ordered) using the SMITS (Simple Multi-dimensional Iterative Technique for Sub-sampling) algorithm. Moving the slider will change the number of data points retrieved and displayed in the plot on the right.

Manipulate[
Module[
{numZones,maxRadius,minRadius,numSamples,zoneWidth,lumpyData,uniformData,randomData,wts},

(* balancedRandomData,indexedWts should also be localized; however, Manipulate cannot handle objects which are not on the global context when nested *) 

numZones = 5;
maxRadius =5;
minRadius =3;
numSamples = 700;

zoneWidth = ((2 π)/(2* numZones));

lumpyData=(#1*{Cos@#2,Sin@#2})&@@@({#1,#2+#3}&@@@Table[{RandomReal[{minRadius,maxRadius}],2*zoneWidth*RandomInteger[{1,numZones}],RandomReal[(zoneWidth/4){-1,1}]},{numSamples}]);
uniformData = (#1*{Cos@#2,Sin@#2})&@@@Table[{RandomReal[{minRadius,maxRadius}],RandomReal[{0,2π}]},{Round[(numSamples/2)]}];

randomData = Join[lumpyData,uniformData];

{wts,balancedRandomData} = BalanceData[
randomData,
InformationMetric -> metric,
Neighborhood->neighborhoodSize
];
indexedWts = Transpose@{Range@Length@wts,wts};
];

Manipulate[
GraphicsRow[{
ListPlot[
{
Take[indexedWts,-(Length@indexedWts-n) /. 0->-1],
Take[indexedWts,n]
},
AspectRatio ->1,
Frame->True,
FrameLabel -> {"balanced sample number","Information Content"},
PlotStyle -> {{Blue,PointSize@Small},{Red,PointSize@Medium}},
GridLines ->Automatic,
PlotLabel ->LabelForm[{"info content =",indexedWts[[n,-1]]},Joined ->True],
PlotRange -> {All,{0,1}}
],
ListPlot[
Take[balancedRandomData,n],
PlotStyle -> {Red,PointSize@Medium},
AspectRatio ->1,
PlotRange ->{{-5.05,5.05},{-5.05,5.05}},
PlotLabel ->LabelForm[{"number of pts =",n},Joined ->True]
]
}]
,
{{n,Length@balancedRandomData},1,Length@balancedRandomData,1},

Deployed->True,
TrackedSymbols:>None
]
,
{{metric,NeighborhoodAsymmetry},{NeighborhoodAsymmetry,LocalNonlinearity,AverageNeighborDistance}},
{{neighborhoodSize,5},Range[2,10],ControlType -> SetterBar},
Deployed->True,
TrackedSymbols-> None
]

Information Content — the country perspective

The information content distribution

countryInformation = GridTable[{{
Module[
{indexWts},
indexWts =BalanceDataIndices[inputMatrix,targetedResponse];
ListPlot[
MapThread[Tooltip[#2,countries[[inputMatrix[[#1,1]]]]]&,Transpose@indexWts],
Frame -> True,
FrameLabel ->{"Ranking","Incremental Information"},
PlotLabel ->LabelForm[{"Cumulative information of",Length@inputMatrix,"countries","based upon",Length@inputLabels,"attributes"},Joined->True],
GridLines->Automatic,
PlotRange ->{All,{0,1}},
PlotStyle ->{Red,PointSize@Medium},
ImageSize->400
]
],
Module[
{indexWts},
countryWts = indexWts =BalanceDataIndices[inputMatrix,targetedResponse,Normalize ->False];
ListLogPlot[
MapThread[Tooltip[#2,countries[[inputMatrix[[#1,1]]]]]&,Transpose@indexWts],
Frame -> True,
FrameLabel ->{"Ranking","Incremental Information"},
PlotLabel ->LabelForm[{"Information Metric from",Length@inputLabels,"attributes from",Length@inputMatrix,"countries"},Joined->True],
GridLines->Automatic,
PlotRange ->All,
PlotStyle ->{Red,PointSize@Medium},
ImageSize ->400
]
]
}},Label-> "The country (record) Information Content"];

topAndBottomCountries = GridTable[{{
GridTable[
{Tooltip[#,LabelForm[#,FontSize->18]]&@countries[[inputMatrix[[#1,1]] ]],#2}&@@@Take[countryWts,20],
TableHeadings->{Automatic,{"Rank","Attribute","Information Score"}},
Label ->"Top Countries"
],
GridTable[
{Tooltip[#,LabelForm[#,FontSize->18]]&@countries[[inputMatrix[[#1,1]] ]],#2}&@@@Take[countryWts,-20],
TableHeadings->{Automatic,{"Rank","Attribute","Information Score"}},
Label ->"Bottom Countries"
]
}},Label->LabelForm["The top & bottom ranked countries (data records)",FontColor -> Red, FontSize->18]];

topCountryBivariatePlot = BivariatePlot[
Transpose[RescaleData@MergeInputResponseData[inputMatrix,targetedResponse]][[All,#]],
DataVariableLabels -> countries[[inputMatrix[[#,1]]]],
PlotLabel->LabelForm["Country Comparison (after application of RescaleData)", FontSize ->14, FontColor->Blue]
]&@Flatten@Take[countryWts,6,1];

colorFunc[country_]:= Switch[country,
"Lesotho",Red,
_,Gray];
lesothoPlot=  WorldPlot[{Africa,colorFunc}];


TabView[{
"Record (country) Information" ->countryInformation,
"Top & Bottom Countries" -> topAndBottomCountries,
"BivariatePlot"-> topCountryBivariatePlot,
"Lesotho?"->lesothoPlot
}]

Notes

Introduction

BalanceData lets us rank data records — in this case, countries — in terms of their generality. This is shown in the plot below where the incremetal information added by considering additional countries is displayed. The countries are ranked in terms of their incremental information contribution; mousing over the dots displays the country name. As we might expect, given the number of attributes, we have a fairly long tail of information gain.

Information Content Profiles

The above left curve displays the incremental normalized cumulative information benefit from adding additional countries. If we look at the unnormalized evaluation of the data contribution, shown to the right, we get a slightly different view of the data. This view is useful for outlier detection; in this case, we don't see extreme deviations so, as we should expect from curated data, there doesn't appear to be erroneous data. As we shall eventually see, however, there are countries which are not easily modeled and, therefore are outliers from the point of view of typical behavior.

The Top & Bottom Countries

The SMITS algorithm used for the data balancing is an iterative scheme wherein the most typical (least incremental information content) is iteratively removed. This multi-pass scheme means that if we select the top K data records we will always have a balanced data set. From the first graph, it appears that by the time we get to the top 20 countries, we have captured over 60% of the information content of the 135 countries being considered. Looking at this selection, we would conclude that we have a reasonable cross-section of the world —with a country like Leosotho representing a whole segment of countries!

BivariatePlot

We compare the top selected countries. We need to preprocess the data using DataModeler‘s RescaleData function since the individual attributes have wildly different natural ranges. The key in looking at this is that the select countries are (pretty much) uncorrelated in their mutual behavior..

Information Content — the attribute perspective

attributeInfoProfile=GridTable[{{Module[
{indexWts},
attributeWts = indexWts =BalanceDataIndices[Transpose@RescaleData@MergeInputResponseData[inputMatrix,targetedResponse],
Normalize -> False,
SelectionFunction -> None];
ListLogPlot[
MapThread[Tooltip[#2,MergeInputResponseData[inputLabels,targetLabel][[#1]]]&,Transpose@indexWts],
Frame -> True,
FrameLabel ->{"Ranking","Information Score"},
PlotLabel ->LabelForm[{"Information scores of",Length@inputLabels+1,"attributes","based upon",Length@inputMatrix,"countries"},Joined->True],
GridLines->Automatic,
PlotRange ->All,
PlotStyle ->{Red,PointSize@Medium},
ImageSize ->400
]
],Module[
{indexWts}, indexWts =BalanceDataIndices[Transpose@RescaleData@MergeInputResponseData[inputMatrix,targetedResponse],
SelectionFunction -> None];
ListPlot[
MapThread[Tooltip[#2,MergeInputResponseData[inputLabels,targetLabel][[#1]]]&,Transpose@indexWts],
Frame -> True,
FrameLabel ->{"Ranking","Cumulative Information"},
PlotLabel ->LabelForm[{"Cumulative information of",Length@inputLabels+1,"attributes","based upon",Length@inputMatrix,"countries"},Joined->True],
GridLines->Automatic,
PlotRange ->{0,1},
PlotStyle ->{Red,PointSize@Medium},
ImageSize ->400
]
]}},
Label-> "The Attribute Information Content"];

topAndBottomAttributes = GridTable[{{
GridTable[
{
Tooltip[#,LabelForm[#,FontSize->18]]&@MergeInputResponseData[inputLabels,"GDPPerCapita"][[#1]],
MergeInputResponseData[AutoSymbolList["x",Length@inputLabels],"y"][[#1]],
#2}&@@@Take[attributeWts,20],
TableHeadings->{Automatic,{"Rank","Attribute","symbol","Information Score"}},
Label ->"Top Attributes for Country Discrimination"
],
GridTable[
{
Tooltip[#,LabelForm[#,FontSize->18]]&@MergeInputResponseData[inputLabels,"GDPPerCapita"][[#1]],
MergeInputResponseData[AutoSymbolList["x",Length@inputLabels],"y"][[#1]],
#2}&@@@Take[attributeWts,-#],
TableHeadings->{Take[Range[Length@inputLabels+1],-#],{"Rank","Attribute","symbol","Information Score"}},
Label ->"Bottom Attributes for Country Discrimination"
]&@20
}},
Label->LabelForm["The top & bottom ranked attributes",FontColor -> Red, FontSize->18]];



TabView[{
"Variable (attribute) Information" -> attributeInfoProfile,
"Top & Bottom Attributes"-> topAndBottomAttributes
}]

Notes

Attribute information

Now let us turn the information on its side and look at the data from the perspective of attributes rather than countries. The key here is that we need to RescaleData before we do the data balancing to the various attributes are on a common scale — otherwise, the attribute with the largest natural range would dominate the data balancing; the corrolary is that we need to turn off the automatic column rescaling since the data is already normalized prior to being transposed. Again, mousing over the attributes will show a tooltip label; however, the key observation is that we do not have an extreme dynamic range — which implies that we are not seeing a lot of redundant attributes.

Top & Bottom Attributes

We can also look at the cumulative informaton content of the various attributes. The absence of attribute dominance is confirmed by the slow cumulative information behavior in the plot below. The attributes, indeed, appear to be measuring different things.

Information Content & Data Balancing — Summary

We have a means to create a balanced data set from unbalanced data

Balanced data avoids inappropriate weighting of data regions ⟹ better models AND faster modeling

We can also detect (some) outliers during the data balancing

The SMITS algorithm does not assume a data model (just an information metric)

Identifying the dominant records can be insightful — these are the prototypes

NKM — A New Kind of Modeling

Also known as Symbolic Regression via Genetic Programming

Start with random expressions

Mutate & breed expressions based upon how well they fit the data

Two goals: accuracy AND simplicity (avoid over or under fitting)

Evolving BOTH model structure AND coefficients — this is NOT conventional regression

Diverse model forms will be discovered

Models are human interpretable

Trustable models — detect operating in new regimes or fundamental system changes

A quick SymbolicRegression

A quick modeling run (30 seconds) — looking for GDPPerCapita

TabView[{
"A Quick Modeling Run"-> ParetoFrontLogPlot[
quickModels = SymbolicRegression[inputMatrix,targetedResponse,
RunsPerCascade ->{4,1,1},
TimeConstraint-> 30
],
DataVariableLabels-> inputLabels,
ImageSize -> 600
],
"The Dominant Variables" -> GridTable[{{
VariablePresenceMap[
quickModels,
BoxRegion -> {50,0.5},
ImageSize ->400
],
VariablePresenceMap[
quickModels,
BoxRegion -> {50,0.5},
ModelingVariables->Automatic,
DataVariableLabels->inputLabels,
ImageSize ->400,
AspectRatio->1
]
}},Label->"Driving variables identified during the quick modeling"]
}]

Lots of 10 minute runs — Looking for GDPPerCapita

allDataView=GridTable[{{
ParetoFrontLogPlot[#,
PlotLabel->LabelForm[{Length@#,"models retrieved from disk"},Joined->True],
DataVariableLabels-> inputLabels
]&[
archivedModels = UniqueAndFitModels@MergeModelSets@RetrieveModelSets[ProjectName-> "GDPPerCapita"]
],(
interestingModels = SelectModels[
archivedModels,
BoxRegion->{#1,#2}
];
ParetoFrontPlot[interestingModels,
PlotLabel->LabelForm[{Length@interestingModels,"interesting models selected from BoxRegion ->",ToString@{#1,#2}},Joined->True],
DataVariableLabels-> inputLabels,
ToolTipLimit->All
]
)&[120,0.04]
}}];

allDataVariablePresence =  GridTable[{{
VariablePresenceTable[
interestingModels,
DataVariableLabels->inputLabels,
Label->LabelForm@"Variables present in >5% of models",
SignificanceLevel -> 0.05
],
VariablePresenceMap[
interestingModels,
Mesh -> None,
AspectRatio -> 1,
PlotLabel->LabelForm@"Variable presence in the interesting models"
]
}},
Label-> "VariablePresence in interesting models (Target is GDPPerCapita)"
];

nichedModelDimensionality = NicheModels[
interestingModels,
Split->ModelDimensionality
];

dimensionalNiche1=GraphicsGrid[{{BoxWhiskerPlot[Sequence@@(
Last[Transpose@#]&/@ModelFitness/@nichedModelDimensionality
),
AspectRatio -> 1,
BoxOrientation ->Vertical,
FrameLabel -> {"Dimensionality","1-(R^2)"},
PlotLabel ->LabelForm@"Error vs. Dimensionality",
BoxLabels-> ModelDimensionality/@First/@nichedModelDimensionality
],
BoxWhiskerPlot[Sequence@@(
First[Transpose@#]&/@ModelFitness/@nichedModelDimensionality
),
AspectRatio -> 1,
BoxOrientation ->Horizontal,
FrameLabel -> {"Complexity","Dimensionality"},
PlotLabel ->LabelForm@"Complexity vs. Dimensionality",
BoxLabels-> ModelDimensionality/@First/@nichedModelDimensionality
]
}},
ImageSize ->800,
PlotLabel->LabelForm["Niching by Dimensionality",FontSize->18,FontColor->Red]
];

lowDModels = Block[
{presenceRequired = 0.05,
varCountLimit = 5, 
limitedDimensionalityModels, targetedDrivers, excludedInputs, focusedModels},

limitedDimensionalityModels = SelectModels[interestingModels, 
NumberOfVariables->varCountLimit,
SelectionStrategy->AllModels];
targetedDrivers = VariablePresence[interestingModels,
PresenceMetric -> "Drivers" ,
SignificanceLevel ->presenceRequired];
excludedInputs = Complement[AutoSymbolList["x",Length@inputLabels],targetedDrivers];

limitedDimensionalityAndDriverModels =focusedModels = SelectModels[limitedDimensionalityModels, 
ExcludedVariables->excludedInputs,
SelectionStrategy->AllModels];

GridTable[{{
VariablePresenceMap[
focusedModels,
ModelingVariables ->Intersection[AutoSymbolList["x",Length@inputLabels],targetedDrivers] ,
DataVariableLabels->inputLabels,
PlotLabel -> LabelForm[{"Variables in >",Round[100*presenceRequired],"% of models having at most",varCountLimit,"variables"},Joined ->True],
AspectRatio->1
],
ParetoFrontPlot[
limitedDimensionalityAndDriverModels,
DataVariableLabels->inputLabels,
PlotLabel -> LabelForm[{Length@limitedDimensionalityAndDriverModels,"Models with limited dimensionality & dominant variables"},Joined->True]
]
}}]
];

varLabelMapping =MapThread[
Rule[#1,#2]&,
{AutoSymbolList["x",Length@inputLabels],inputLabels}
];

popularSubspaceNiching = GridTable[{{
ListLogPlot[
Tooltip[Length@#,ModelSubspace@First@# /. varLabelMapping]&/@#,
PlotStyle->{PointSize@Large,Red},
Frame ->True,
FrameLabel->{"subspace index","# of models in subspace"},
PlotLabel->LabelForm[{Length@#,"different subspaces are present"},Joined->True],
GridLines->Automatic,
ImageSize->400,
PlotRange ->All
]&[
nichedModelSubspace=Reverse@SortBy[
NicheModels[
interestingModels,
Split->ModelSubspace
],Length@#&
]
],
ListLogPlot[
Tooltip[{ModelDimensionality@First@#,Length@#},ModelSubspace@First@# /. varLabelMapping]&/@#,
PlotStyle->{PointSize@Large,Red},
Frame ->True,
FrameLabel->{"subspace dimensionality","# of models in subspace"},
PlotLabel->LabelForm[{Length@#,"different subspaces are present"},Joined->True],
GridLines->Automatic,
PlotRange ->All,
ImageSize->400
]&@nichedModelSubspace
}},
Label->"Niching by Subspace"
];

dominantVarNiche=GridTable[{{
Block[
{presenceThreshold = 0.5},
dominantVariables = VariablePresence[
interestingModels,
PresenceMetric ->"Drivers",
SignificanceLevel ->presenceThreshold
];

ParetoFrontPlot[#,
PlotLabel->LabelForm[{Length@#,"models containing only the variables with at least",Round[100*presenceThreshold],"% presence"},Joined->True],
DataVariableLabels-> inputLabels
]&[
dominantVariableModels =SelectModels[interestingModels,
SelectionStrategy -> AllModels,
ExcludedVariables-> Complement[
AutoSymbolList["x",Length@inputLabels],
dominantVariables
]
]
]
],
ParetoFrontPlot[#2,
PlotLabel->LabelForm[
{#1,"seconds to produce",Length@#2,"unique models after optimization"},Joined->True,
NumberFormatting->{5,1}
],
DataVariableLabels-> inputLabels
]&@@Timing[
optimizedDominantVarModels=
Quiet@UniqueAndFitModels@RemoveIntrons@OptimizeModel[
dominantVariableModels,
inputMatrix,
targetedResponse,
SimplificationFunction -> (FullSimplify@#&)
]
]
}}];


TabView[{
"Many 10 Minute Runs" -> allDataView,
"Variable Presence" ->allDataVariablePresence,
"Dimensionality Niching" -> dimensionalNiche1,
"Low-D Models" ->  lowDModels,
"Popular Subspaces" -> popularSubspaceNiching,
"Dominant Vars" -> dominantVarNiche
}]

Model-Based Variable Selection & Niching Summary

Despite the fact that we have a fat array and correlated inputs, we managed to identify and isolate driving variables

Diverse model structures were uncovered that fit the data

The discovered models featured different numbers of models (dimensionality) and different variable combinations (subspaces)

We can select models for use which have limited dimensionality or targeted variable combinations

Model-Based Outlier Detection — The selected models

Candidate Models …

Low-dimensional models

Popular subspace

Identify dominant variables & focus on those models from candidates

Build an ensemble of diverse models from focused subset

ensembleVarDisplay = Module[
{varCountLimit = 4, 
subspacePopularityThreshold = 20,
dominantModelThreshold = 0.1,
lowDimensionalModels,popularSubspaceModels, compositeModelSet, dominantVariables,excludedVariables
},
lowDimensionalModels = SelectModels[interestingModels, 
NumberOfVariables->varCountLimit,
SelectionStrategy->AllModels];
popularSubspaceModels = MergeModelSets@Select[
NicheModels[
interestingModels,
Split->ModelSubspace
],
Length@#≥subspacePopularityThreshold&
];
compositeModelSet = UniqueModels@MergeModelSets[
lowDimensionalModels,
popularSubspaceModels
];
dominantVariables = VariablePresence[
compositeModelSet,
PresenceMetric ->"Drivers",
SignificanceLevel ->dominantModelThreshold
];
excludedVariables= Complement[
AutoSymbolList["x",Length@inputLabels],
dominantVariables
];
candidateModels= SelectModels[compositeModelSet, 
ExcludedVariables -> excludedVariables,
SelectionStrategy->AllModels];

ensemble=CreateModelEnsemble[
candidateModels,
inputMatrix,
targetedResponse
];

GridTable[{{
VariablePresenceTable[candidateModels,
DataVariableLabels->inputLabels,
Label -> LabelForm[{"Candidate model variable set"},Joined ->True]],
VariablePresenceMap[
ensemble,
ModelingVariables -> Automatic,
DataVariableLabels-> inputLabels,
Mesh ->True,
AspectRatio -> 1,
PlotLabel ->LabelForm@"Ensemble Variables"
]
}}]
];

ensembleModelReport = ModelSelectionReport[
ensemble,
DataVariableLabels->LabelForm@inputLabels,
Label->LabelForm["Models in the Selected Ensemble",FontSize -> 18]
];

outlierIndices =First/@(
outlierDetails = DataOutlierAnalysis[ensemble,inputMatrix,targetedResponse,SignificanceLevel->{Automatic,3}]
);
outlierRecords = inputMatrix[[outlierIndices]];
outlierResponse = targetedResponse[[outlierIndices]];

overPerformers = Flatten@Position[
EvaluateEnsemblePrediction[ensemble,outlierRecords] - outlierResponse,
a_/;(a< 0)];
underPerformers = Complement[Range@Length@outlierResponse,overPerformers];

ensemblePredictions = GridTable[{{
Show[
EnsemblePredictionPlot[ensemble,inputMatrix,targetedResponse,
ToolTipFunction->(countries[[ inputMatrix[[#,1 ]] ]]&),
PlotStyle->{{Darker@Green,PointSize[0.015]}}],
EnsemblePredictionPlot[
ensemble,outlierRecords[[overPerformers]],outlierResponse[[overPerformers]],
ToolTipFunction->(countries[[ outlierRecords[[overPerformers]][[#,1 ]] ]]&),
PlotStyle->{{Blue,PointSize[0.02]}}],
EnsemblePredictionPlot[
ensemble,outlierRecords[[underPerformers]],outlierResponse[[underPerformers]],
ToolTipFunction->(countries[[ outlierRecords[[underPerformers]][[#,1 ]] ]]&),
PlotStyle->{{Red,PointSize[0.02]}}],
PlotLabel ->LabelForm["\!\(\*
StyleBox[\"Underpredicted\",\nFontColor->RGBColor[0, 0, 1]]\) and \!\(\*
StyleBox[\"Overpredicted\",\nFontColor->RGBColor[1, 0, 0]]\) countries",FontSize ->14],
AspectRatio->1
],
GridTable[
LabelForm[{
countries[[First@#4]],
#3,
If[#3<0,"Underpredicted","Overpredicted"]
},
FontColor -> If[#3<0,Blue,Red]
]&@@@outlierDetails,
Label->"GDP per Capita Model-based Outliers",
TableHeadings -> {Automatic, {"oulier rank","country","oulier distance","classification"}}
]
}}];

TabView[{
"Ensemble Variables" -> ensembleVarDisplay,
"Ensemble Models" -> ensembleModelReport,
"Ensemble Predictions" -> ensemblePredictions
}]

Outlier Detection Summary

Outlier detection inherently HAS to be model-based

ArrayPlot[
test =Transpose@RescaleData[
MergeInputResponseData[inputMatrix,targetedResponse],
SelectionFunction->({Min@#,Max@#}&)],
ColorFunction->(Hue[(2/3)#]&),
FrameLabel -> LabelForm@{"Attributes","Countries"},
PlotLabel->LabelForm["Which countries are outliers?",FontSize->18],
ImageSize ->600
]

The outlier is either VERY important or VERY distracting

This is a chicken & egg problem

Good models make outlier detection easy

Outliers make good model development hard

The answer is not to develop THE model but to use many diverse models

The more models that flag a data record as being an outlier, the more likely it IS an outlier

Ultimately, deciding a data records is an outlier is an external decision

Jewel or Junk — Summary

Two views: information content & model-based

Each view is powerful & insightful

Iteratively integrating these views is also a good idea

Initialization

Load DataModeler

<<DataModeler`

$VersionOfDataModeler

Load WorldPlot

<<WorldPlot`;

Load the raw CountryData[ ] column & row headings

The countries of the world — countries

We have 237 countries

Length[
countries =CountryData["Countries"]
]

The (possible) properties of those countries — dataFields

We have 225 attributes associated with the countries — not all of them numeric scalars and some of them missing.

Length[
dataFields=CountryData["Properties"]
]

Retrieve all of the available information — rawData

Note that the initial retrieval of the data from WRI can take MUCH longer. Fortunately, once it is on your machine it shouldn't be too onerous to actually extract the data.

Timing@Dimensions[
rawData = Outer[CountryData[#1, #2]&,countries,dataFields]
]

Retrieve the data for the USA

Also, out of an expectation that the US will have the most complete data set, let us pull out the US data and check for missing information. The result shows some surprising omissions

usaData=CountryData["UnitedStates",#]&/@dataFields;
TableForm@Cases[MapThread[List,{dataFields,usaData}],{_,Missing[__]}]

Augment the raw data

The SymbolicRegression[ ] wants to see a matrix of scalars. However, there is some information in the available data which is embedded as rules within specific data fields.

Look at the available extra fields

TableForm@Cases[MapThread[List,{dataFields,usaData}],{_,{__Rule}}]

Extract the interesting extra fields and define labels

{additionalDataFields,additionalRawData}=Block[
{interestingFields,ruleBasedData, ruleFields, ruleList,hiddenDataMatrix, hiddenDataLabels},
interestingFields = {"ExpenditureFractions","GDPSectorFractions","SectorLaborFractions","PavedAirportLengths","Pipelines","ReligionsFractions"};
ruleBasedData=Outer[CountryData[#1, #2]&,countries,interestingFields];
ruleFields =DeleteCases[
Union/@Flatten/@ Transpose[ruleBasedData /. ruleList:{__Rule} :> ruleList[[All,1]]],
Missing[_String],∞];
hiddenDataMatrix = Join[
Sequence@@(
MapThread[
#1 /. #2&,
{ruleFields, Transpose@ruleBasedData /. a:Missing[_]:>{2->2 }}] /. _String -> Missing["NotAvailable"
]
),
2];
hiddenDataLabels = Flatten@MapThread[#2 /. a_String:> StringJoin[#1,"•",a]&,{interestingFields,ruleFields}];
{hiddenDataLabels, hiddenDataMatrix}
] ;
additionalDataFields

Form an augmentedRawData set & associated labels

Below we merge the additional data with the original raw set and also add an index into the records so we can later associate records with their source country.

augmentedDataFields = Join[{"CountryIndex"},dataFields,additionalDataFields];
Dimensions[
augmentedRawData = Join[Transpose@{Range@Length@rawData},rawData,additionalRawData,2]
]