PCMCI+ with RegressionCI

[DavidQiwi]

Hello everyone,

I'm currently using PCMCI+ as part of my bachelor's thesis on process optimization in injection molding using causal machine learning. I encountered a limitation when applying the Regression CI method with discrete variables.

Some of my variables have exactly two categories, and others have three. The issue arises when I set three discrete variables (regardless of their specific values): the method throws an error and fails to run (ValueError: Unknown label type: continuous. Maybe you are trying to fit a classifier, which expects discrete classes on a regression target with continuous values). However, if I only specify two variables as discrete—regardless of whether they have two or three categories—Regression CI runs without problems and yields reasonable results.

Is there a known limitation in Regression CI regarding the number of discrete variables that can be used simultaneously? Could this be related to how the data are formatted or encoded?

I’d appreciate any clarification or suggestions on how to handle this case.

Thanks in advance!

[jakobrunge]

Can you check on the latest master? I fixed a bug last week...

[DavidQiwi]

Hello,

Thank you for your recommendation. I have upgraded tigramite version and tried to set 3 discrete values in my data. The method throws an error (ValueError: This solver needs samples of at least 2 classes in the data, but the data contains only one class: 17). The number of samples mentioned in the error differ, depending on variables, which were set to discrete. These discrete values have even more than 2 classes (more than 2 states). That is why, I could not understand the problem and why it occurs.

It would be great if you could give me adivce!

[jakobrunge]

Can you post a minimal code example with synthetic data?

[DavidQiwi]

Hello!
It has required considerable time to thoroughly explore the method and gain a deeper understanding of its underlying principles. Example code is provided alongside representative data, which is intended to closely reflect real-world data relevant to the specific use case.
import numpy as np
import pandas as pd

Set random seed for reproducibility

np.random.seed(42)

Number of time points (the actual number of samples in the real world data is 11000)

n = 1000

def add_outliers(data, outlier_values, num_outliers):
indices = np.random.choice(len(data), size=num_outliers, replace=False)
outliers = np.random.choice(outlier_values, size=num_outliers)
data[indices] = outliers
return data

Variable 1: 710–730 with outliers at 690

var1 = np.random.uniform(710, 730, n)
var1 = add_outliers(var1, [690], 10)

Variable 2: 2680–2730 with outliers at 2600

var2 = np.random.uniform(2680, 2730, n)
var2 = add_outliers(var2, [2600], 10)

Variable 3: 160.15–160.50 with outliers in 161.0–161.25

var3 = np.random.uniform(160.15, 160.50, n)
var3[:10] = np.random.uniform(161.0, 161.25, 10) # inject outliers at the start

Variable 4: 450–470 (no outliers)

var4 = np.random.uniform(450, 470, n)

Variable 5: 8–8.7 with outliers in 9.0–9.5

var5 = np.random.uniform(8.0, 8.7, n)
var5[:10] = np.random.uniform(9.0, 9.5, 10)

Variable 6: 178.45–178.67

var6 = np.random.uniform(178.45, 178.67, n)

Variable 7: around 12.6 with little noise

var7 = np.random.normal(12.6, 0.01, n)

Variable 8: 2422–2428

var8 = np.random.uniform(2422, 2428, n)

Variable 9: discrete {0.280, 0.290, 0.300}

var9 = np.random.choice([0.280, 0.290, 0.300], size=n)

Variable 10: discrete {16.01, 16.02}, with more 16.01

var10 = np.random.choice([16.01, 16.02], size=n, p=[0.8, 0.2])

Variable 11: 85–95 with some outliers down to 60

var11 = np.random.uniform(85, 95, n)
var11 = add_outliers(var11, [60], 10)

Variable 12: binary {0, 1}, mostly 0

var12 = np.random.choice([0, 1], size=n, p=[0.85, 0.15])

Time index (hourly) (for clarity)

time_index = pd.date_range(start='2023-01-01', periods=n, freq='H')

df = pd.DataFrame({
'var1': var1,
'var2': var2,
'var3': var3,
'var4': var4,
'var5': var5,
'var6': var6,
'var7': var7,
'var8': var8,
'var9': var9,
'var10': var10,
'var11': var11,
'var12': var12,
}, index=time_index)

[jakobrunge]

There might still be a bug, but a remark: if all your variables are categorical, then rather use the Gsquared or CMIsymb CI tests.

[DavidQiwi]

Hello,

There are only three categorical variables in the time series. When two of them are specified as discrete (i.e., data_type[:,i]=1), the method performs as expected. However, when all three variables are set as discrete, the method does not behave correctly.

Additionally, I encountered an issue when using the regression-based conditional independence (CI) test to obtain the causal graph: the resulting graph contains cycles and is therefore not a CPDAG. As I understand from the corresponding paper, the method should produce a CPDAG, which is acyclic and partially directed.

Could this behavior be due to violated assumptions of the method? Or might it be related to the treatment of discrete variables?

[jakobrunge]

Hi, I am on holidays this week and can have a look next week. It must be related to the sklearn logistic regression function and not inherently to tigramite.

Regarding the cycles. In the asymptotic case, PCMCI+ results in CPDAGs, but for finite data not necessarily. Does the true process underlying your synthetic data contain feedbacks?

[jakobrunge]

Regarding your example: I don't see how you used PCMCI+ and regressionCI, please provide the full code. For example, did you correctly specify the data_type in the Data frame?

[DavidQiwi]

Hello,

anyway, the data are finite and contain more than 10000 entries (rows). The true process can actually contain feedbacks, because the variables are from manufacturing machine. The initial aim was to receive a CPDAG (or a DAG in the best case) and to calculate causal effects with Causal Inference methods. However, it is not possible to calculate causal effect on a cyclic graph. It was also possible to try ParrCorr CI test, when discrete variables were dropped (as they refer to machine's output or labels), which output a CPDAG. When using RegressionCI, there are no discrete variables in the cycles (loops), but only continuous variables.

[DavidQiwi]

`import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from tigramite.pcmci import PCMCI
from tigramite.independence_tests.regressionCI import RegressionCI
import tigramite.plotting as tp
import tigramite.data_processing as pp

#Set random seed for reproducibility
np.random.seed(42)

#Number of time points (the actual number of samples in the real world data is more than 10000 but less than 15000)
n = 1000

def add_outliers(data, outlier_values, num_outliers):
indices = np.random.choice(len(data), size=num_outliers, replace=False)
outliers = np.random.choice(outlier_values, size=num_outliers)
data[indices] = outliers
return data

#Variable 1: 710–730 with outliers at 690
var1 = np.random.uniform(710, 730, n)
var1 = add_outliers(var1, [690], 10)

#Variable 2: 2680–2730 with outliers at 2600
var2 = np.random.uniform(2680, 2730, n)
var2 = add_outliers(var2, [2600], 10)

#Variable 3: 160.15–160.50 with outliers in 161.0–161.25
var3 = np.random.uniform(160.15, 160.50, n)
var3[:10] = np.random.uniform(161.0, 161.25, 10) # inject outliers at the start

#Variable 4: 450–470 (no outliers)
var4 = np.random.uniform(450, 470, n)

#Variable 5: 8–8.7 with outliers in 9.0–9.5
var5 = np.random.uniform(8.0, 8.7, n)
var5[:10] = np.random.uniform(9.0, 9.5, 10)

#Variable 6: 178.45–178.67
var6 = np.random.uniform(178.45, 178.67, n)

#Variable 7: around 12.6 with little noise
var7 = np.random.normal(12.6, 0.01, n)

#Variable 8: 2422–2428
var8 = np.random.uniform(2422, 2428, n)

#Variable 9: discrete {0.280, 0.290, 0.300}
var9 = np.random.choice([0.280, 0.290, 0.300], size=n)

#Variable 10: discrete {16.01, 16.02}, with more 16.01
var10 = np.random.choice([16.01, 16.02], size=n, p=[0.8, 0.2])

#Variable 11: 85–95 with some outliers down to 60
var11 = np.random.uniform(85, 95, n)
var11 = add_outliers(var11, [60], 10)

#Variable 12: binary {0, 1}, mostly 0
var12 = np.random.choice([0, 1], size=n, p=[0.85, 0.15])

#Time index (hourly) (for clarity)
time_index = pd.date_range(start='2023-01-01', periods=n, freq='H')

df = pd.DataFrame({
'var1': var1,
'var2': var2,
'var3': var3,
'var4': var4,
'var5': var5,
'var6': var6,
'var7': var7,
'var8': var8,
'var9': var9,
'var10': var10,
'var11': var11,
'var12': var12,
}, index=time_index)

data_type = np.zeros(df.shape, dtype='int')
data_type[:,0] = 0
data_type[:,1] = 0
data_type[:,2] = 0
data_type[:,3] = 0
data_type[:,4] = 0
data_type[:,5] = 0
data_type[:,6] = 0
data_type[:,7] = 0
data_type[:,8] = 1 # Discrete
data_type[:,9] = 1 # Discrete
data_type[:,10] = 0
data_type[:,11] = 1 # Discrete

vector_vars = {i: [(i, 0)] for i in range(len(df.columns))}
data_set_param_np = df.to_numpy()
dataframe = pp.DataFrame(data=data_set_param_np,
data_type=data_type,
vector_vars=vector_vars,
analysis_mode='single',
var_names=df.columns.to_list())

regression_ci = RegressionCI(significance='analytic')
pcmci = PCMCI(dataframe=dataframe, cond_ind_test=regression_ci, verbosity=0)

tau_min = 0
tau_max = 0
pc_alpha = 0.01
pcmci.verbosity = 0
var_map = {name: idx for idx, name in enumerate(dataframe.var_names)}
link_assumptions = {j:{(i, -tau):'o?o' for i in var_map.values() for tau in range(tau_max + 1) if (i, -tau) != (j, 0)} for j in var_map.values()}

results = pcmci.run_pcmciplus(tau_min=0, tau_max=0, pc_alpha=0.01, link_assumptions=link_assumptions, contemp_collider_rule='conservative')

tp.plot_graph(val_matrix=results['val_matrix'],
graph=results['graph'],
var_names=df.columns.to_list(), figsize=(12, 12),
node_size=0.1, node_label_size=7,
arrow_linewidth=4, show_colorbar=False)
plt.show()`

[jakobrunge]

The error was that your categorical variable contained only one class and sklearn cannot deal with this. I modified the deviance score to return 0.0 in this case. Check again.

[DavidQiwi]

Hello!

Do you mean, that if the variable is constant, then the sklearn package with linear and logistic regressions will throw an error? The problematic variable in my use case has only two distinct values over the whole dataset, and actually it is not constant. So sklearn package seems to handle it as constant, maybe because these distinct values are right next to each other.

If I set this variable to continuous, then the algorithm works and outputs a causal graph. However, the final graph is not acyclic. This is the main problem for now.

[jakobrunge]

I checked with your code and the resulting target array is indeed only containing one value, check yourself. Note that the array may slightly differ from the input data because the first few samples are cut off depending on taumax.