CAUSALDISCOVERY Procedure

This example shows how to use the CAUSALDISCOVERY procedure to study the performance of the TOP algorithm and focuses on DAGs that contain 100 variables. The number of 100-variable DAGs is huge, and far beyond the maximum number that a double-precision number can express. For a data set that contains 100 variables, finding the DAG that is used in the data generating process (DGP) is an almost impossible challenge.

This example checks the accuracy of the TOP algorithm in three cases of different number of edges:

For each case, 30 trials are examined. In each trial, first the SIMULATE statement simulates a DAG and then generates a 10,000-observation data set based on the simulated DAG; then the LEARN statement estimates the DAG according to the simulated data set; and finally, the EVALUATE statement evaluates the accuracy by comparing the estimated DAG from the LEARN statement to the true simulated DAG from the SIMULATE statement. In the data set simulation, a parameter value that is associated with an edge is uniformly drawn from . This specific range is chosen to avoid the numerical issues for the simulated data set, because when the number of edges increases, especially in the full DAG case, the chance of getting a numerically unstable data set also increases. Here, the numerically unstable data set means that the data set that is expressed as a matrix is not numerically full rank, although in theory it should be full rank. In the accuracy evaluation, the counts of three possible mistakes in DAG estimation—namely, missing edges, extra edges, and reverse edges—are reported; the sum of these three counts, the so-called structural Hamming distance (SHD), is also reported; and the final reported criterion is the true positive rate (TPR) for checking how many truly existing directed edges are found.

The SAS macro for a trial is as follows:

%macro dagTrial(simId,nNodes,nEdges,nObs,tblEval);
   * simulate DAG and data;
   proc causaldiscovery;
      var col1 - col&nNodes.;
      simulate outdata=mylib.sdata out=mylib.sdag / nEdges=&nEdges. nObs=&nObs.
         parm(dist=uniform(lb=-0.75 ub=0.75 excllb=-0.4 exclub=0.4))
	         seed=%eval(12345678+&simId.);
   run;

   * learn the DAG structure;
   proc causaldiscovery;
      var col1 - col&nNodes.;
      learn data=mylib.sdata outdag=mylib.outdag / maxiter=%eval(5*&nNodes)
         alpha=(0.1 0.01 0.0001 0.000001 0.00000001)
	         out=mylib.oest ;
   run;

   * evaluate the accuracy of structure learning;
   proc causaldiscovery;
      var col1 - col&nNodes.;
      evaluate estdag=mylib.outdag truedag=mylib.sdag / out=mylib.oeval;
   run;

   data mylib.oeval;
      set mylib.oeval; simID=&simId.;
   run;

   %if &simId.=1 %then %do;
      %let dataread=mylib.oeval;
   %end;
   %else %do;
      %let dataread= mylib.&tblEval. mylib.oeval;
   %end;

   data mylib.&tblEval.;
      set &dataread;
   run;
%mend;

Alternatively, you can specify all three statements (SIMULATE, LEARN, and EVALUATE) in a single procedure call as follows:

%macro dagTrial(simId,nNodes,nEdges,nObs,tblEval);

   proc causaldiscovery;
      var col1 - col&nNodes.;
      * simulate DAG and data;
      simulate out=mylib.sdag outdata=mylib.sdata /
         nEdges=&nEdges. nObs=&nObs. seed=%eval(12345678+&simId.)
         parm(dist=uniform(lb=-0.75 ub=0.75 excllb=-0.4 exclub=0.4));

      * learn the DAG structure;
      learn outdag=mylib.outdag / maxiter=%eval(5*&nNodes)
         alpha=(0.1 0.01 0.0001 0.000001 0.00000001)
         out=mylib.oest ;

      * evaluate the accuracy of structure learning;
      evaluate / out=mylib.oeval;
   run;

   data mylib.oeval;
      set mylib.oeval; simID=&simId.;
   run;

   %if &simId.=1 %then %do;
      %let dataread=mylib.oeval;
   %end;
   %else %do;
      %let dataread= mylib.&tblEval. mylib.oeval;
   %end;

   data mylib.&tblEval.;
      set &dataread;
   run;
%mend;

The following SAS macros loop over all trials and print the final accuracy report:

%macro dagPrint(tblEval);
   data &tblEval.; set mylib.&tblEval.; run;
   proc sort data=&tblEval.; by alpha simId; run;
   proc univariate data=&tblEval. NOPRINT;
      var shd nme nee nre tpr; by alpha;
      output out=ostat
             mean=shdMean nmeMean neeMean nreMean tprMean
             std=shdStd nmeStd neeStd nreStd tprStd;
   run;
   data pstat;
      set ostat;
      pType = 'Mean'; shd = shdMean; nme = nmeMean;
      nee = neeMean; nre = nreMean; tpr = tprMean;
      output;
      pType = 'Std.'; shd = shdStd; nme = nmeStd;
      nee = neeStd; nre = nreStd; tpr = tprStd;
      output;
      keep alpha pType shd nme nee nre tpr;
   run;
   proc print data=pstat noobs label;
      format alpha f10.8;
      format shd nme nee nre tpr f7.2;
      label alpha='Significance Level'
         pType='Type' shd='Structural Hamming Distance'
         nme='N Missing Edges' nee='N Extra Edges'
         nre='N Reverse Edges' tpr='True Positive Rate';
   run;
%mend;

%macro dagTrials(sTrial,eTrial,nNodes,nEdges,nObs,tblEval);
   %do simId = &sTrial. %to &eTrial.;
      %put "Simulation &simId.";
      %dagTrial(&simId.,&nNodes.,&nEdges.,&nObs.,&tblEval.);
   %end;
   %dagPrint(&tblEval.);
%mend;

First, the 100-variable, 600-edge DAG case is checked using the following code:

%let nNodes = 100;
%let nEdges = 600;
%let nObs = 10000;
%let tblEval = evalSparse;
%dagTrials(1,30,&nNodes.,&nEdges.,&nObs.,&tblEval.);

The results are shown in Output 6.1.1. Because the proper significance level is chosen (in this case, or ), the SHD decreases to zero for all 30 trials, which means that the true 100-variable DAGs in all 30 trials are perfectly discovered and the TOP algorithm reaches 100% accuracy. Such results are not seen in any algorithms that are studied by Yu et al. (2021), who also perform the experiments on the 100-variable, 600-edge DAGs.

Second, the number of edges is increased to 2,400. The SAS code is as follows:

%let nNodes = 100;
%let nEdges = 2400;
%let nObs = 10000;
%let tblEval = evalMedium;
%dagTrials(1,30,&nNodes.,&nEdges.,&nObs.,&tblEval.);

The results are shown in Output 6.1.2. Surprisingly, because the significance level is set to or , the accuracy of the TOP algorithm reaches 100% again.

Finally, the largest possible number of edges, 4,950 edges, for a 100-variable DAG is tried. The SAS code is as follows:

%let nNodes = 100;
%let nEdges = 4950;
%let nObs = 10000;
%let tblEval = evalFull;
%dagTrials(1,30,&nNodes.,&nEdges.,&nObs.,&tblEval.);

The results are shown in Output 6.1.3. Because the significance level is chosen as 0.10 or 0.01, the TOP algorithm reaches 100% accuracy again. In this extreme case, even if other significance levels are chosen, the accuracy is still very high.

In all three cases, the TPR is always 100%, which means that in all 90 trials, the true directed edges are always found. When you run the code in this example on a 2,160-CPU grid, each trial takes only a few seconds. Note that the TOP algorithm is not an exact algorithm. Although in this example the TOP algorithm always reaches 100% accuracy for some chosen significance levels, it is not guaranteed that the TOP algorithm can reach 100% accuracy for other data sets. However, this example does show how fast, and especially how accurate, the TOP algorithm is.