I need help using the corrcoef function to calculate PRx (pressure reactivity index)

Hi everyone! I need to calculate PRx (pressure reactivity index) using ICP (intracranial pressure) and mean ABP (mean arterial blood pressure) data. PRx is known as "pressure reactivity index is calculated as the degree of statistical correlation between the slow wave components of mean arterial pressure (ABP) and intracranial pressure (ICP)". I need to do this in order to calculate a different hypothetical reactivity index. I have PRx data- so my idea was to use the given PRx data and compare it to the calculated PRx data to ensure it is correct.

This is the script to display the PRx data that was given to me:

%PRx_display

% Load PRx data from MAT file
load('PRx,data_part1of1');

PRx_values = measurement_data;

disp('PRx Values:');
disp(PRx_values);

% Plot PRx values
figure;
plot(PRx_values);
xlabel('Time');
ylabel('PRx Value');
title('PRx');
grid on;

% Calculate average PRx
average_PRx = mean(PRx_values);

% Display average PRx
fprintf('Average PRx: %.4f\\n', average_PRx);

This is the script I am trying to calculate PRx with:

%PRx_calc

% Load ICP and MAP data
load("ICP,data_part1of1");
ICP_time_vector = round(time_vector); % Assuming time_vector in ICP file
ICP_data = measurement_data;

load("ABP,mean,data_part1of1");
ABP_time_vector = round(time_vector); % Assuming time_vector in ABP file
ABP_data = measurement_data;

% Remove duplicate time points from ICP
[ICP_time_vector, uniqueIdx] = unique(ICP_time_vector);
ICP_data = ICP_data(uniqueIdx);

% Remove duplicate time points from ABP
[ABP_time_vector, uniqueIdx] = unique(ABP_time_vector);
ABP_data = ABP_data(uniqueIdx);

% Ensure consistent dimensions
if length(ICP_time_vector) \~= length(ICP_data)
    error('ICP_time_vector and ICP_data have different lengths');
end

if length(ABP_time_vector) \~= length(ABP_data)
    error('ABP_time_vector and ABP_data have different lengths');
end

% Interpolating ABP to match ICP time vector
if \~isequal(ICP_time_vector, ABP_time_vector)
    ABP_data = interp1(ABP_time_vector, ABP_data, ICP_time_vector, 'linear', 'extrap');
end

% Combine ICP data and aligned ABP data
ICP = \[ICP_time_vector(:), ICP_data(:)\]; % Ensure column vectors
ABP = \[ICP_time_vector(:), ABP_data(:)\]; % ABP_data is now interpolated to match ICP_time_vector

% Parameters
windowSize = 30 \* 10; % 30 windows of 10 seconds
stepSize = 60; % Update every 60 seconds
n = length(ICP_time_vector); % Number of data points

% Preallocate PRx array
PRx = NaN(floor((n - windowSize) / stepSize) + 1, 1);

% Compute moving correlation coefficient
for i = 1:stepSize:(n - windowSize + 1)
    % Extract current window of data
    windowICP = ICP(i:i + windowSize - 1, 2);
    windowABP = ABP(i:i + windowSize - 1, 2);
    
    % Calculate correlation coefficient for the current window
    R = corrcoef(windowICP, windowABP);
    
    % Store PRx (correlation coefficient of ICP and ABP)
    PRx(floor(i / stepSize) + 1) = R(1, 2);
    
end

% Calculate average PRx ignoring NaN values
averagePRx = sum(PRx, 'omitnan') / sum(\~isnan(PRx));

% Display average PRx in the command window
fprintf('Average PRx: %.4f\\n', averagePRx);

% Plot PRx vs Time in minutes
figure;
plot((0:length(PRx)-1)\*stepSize/60, PRx); % Convert to minutes
title('PRx vs Time');
xlabel('Time (minutes)');
ylabel('PRx Values');

I thought the calculation script would work, but it doesn't. For instance, for a certain set of data, PRx_calc.m gives me an average PRx value of 0.32, while PRx_display.m gives me 0.52. The trend of the plots look the same, but the PRx_calc.m's for some reason looks shifted down and distorted a little bit. I am looking for help with the way I am using the corrcoef function and general debugging help.

Here are the first 20 values of each type of data: First 20 values of

ICP: 10 10 10 10 11 11 11 11 10 11 10 11 11 11 11 11 10 11 11 11

First 20 values of

ABP: 89 89 89 89 89 89 89 89 89 89 89 89 89 89 89 89 89 89 89 89

First 20 calculated

PRx values: NaN NaN NaN NaN NaN NaN -0.0838 -0.2111 -0.0586 0.1704 0.4276 0.3193 0.3417 0.4128 0.5488 0.4790 0.5207 0.6822 0.7107 0.5235
Average PRx: 0.3209

The PRx_calc script, which produces the above values, should produce the given PRx values(below).

First 20 PRx Values: 0.6150 0.5996 0.6010 0.6276 0.7499 0.4495 0.4488 0.4492 0.5078 0.7263 0.7290 0.7343 0.7490 0.7561 0.6997 0.4609 0.2567 0.5935 0.7317 0.6666
Average PRx: 0.5247

Please let me know if there is anything else I should add, more data points, graphs, etc. Thank you!

Attached are the ABP, ICP, and PRx files.

[dpb]

Attach the .mat file would most help.

Is there a reference that would illuminate the calculation that created the other? Is it known those are actually from the identical data?

There are a bunch of "\" characters in your code that are illegal MATLAB syntax...please correct those so the code could be run as you intend once we have the data.

ICP=[10 10 10 10 11 11 11 11 10 11 10 11 11 11 11 11 10 11 11 11];
ABP=[89 89 89 89 89 89 89 89 89 89 89 89 89 89 89 89 89 89 89 89];
corrcoef(ICP,ABP)

ans = 2x2

1 NaN NaN NaN

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Since there is no change in the ABP variable, there is no defined correlation coefficient for the above dataset--from all the NaN in the beginning of the shown output, we must infer that is probably the case for quite a long stretch of the input data.

One presumes probably the other calculation removed all those data from consideration and may have done other preprocessing as well; we don't know and have no way to tell what might be different.

[William Rose]

Excellent comment from @dpb. As usual.

The reactivity index definition is from this paper.

I have worked with cerebral pressure and arterial blood pressure signals - see here and here.

As @dpb said, post some data.

[Caleb]

@dpb @William Rose

Thank you both for your comments! I have attached the datafiles foer PRx, ICP, and ABP

[William Rose]

@Caleb, Are the time values in seconds?

[Caleb]

@William Rose Yes they are!

Answer 1

[William Rose]

@Caleb,

The pressure reactivity index is defined as "a moving correlation coefficient from 30 consecutive 10-s averages of ICP and ABP waveforms. We called this the PRx index (pressure reactivity index)" here. The authors clarify that they use MAP, not instantaneous BP, for the ABP waveform, which make sense, since they compute the average over 10 seconds. Thus they use data from 5 minutes MAP and ICP for each estimate of PRx. They update the index every minute. Therefor consecutive esitmates of PRx use 4 minutes of overlapping data. PRx is designed to reveal how well cerebral autoregulation is working. They assume, based on other work, that only frequencies of 0.05 Hz and slower carry information about cerebral autoregulation. That is why they sample the data every 10 seconds.

Suppose you have MAP and ICP. I will make two simulated signals, with duration 10 minutes.

% Simulated data
fs=2;           % sampling rate (Hz)
dt=1/fs;        % sampling interval (s)
T=600;          % signal duration (s)
N=T/dt;         % number of samples
t=(0:N-1)*dt;   % time vector (s)
RR=12;          % respiratory rate (breaths per minute)
fRR=RR/60;      % resp rate (Hz)
% Simulated mean arterial pressure: mean=90, +-3 mmHg due to respiration.
map=90+3*sin(2*pi*fRR*t);
% Simulated intracranial pressure: mean=15, +-2 mmHg due to respiration.
icp=15+2*cos(2*pi*fRR*t);

Now we compute the means of MAP and of ICP in 10-second segments:

Ts=10;                      % duration of segment for averaging (s)
M1=floor(N*dt/Ts);          % number of complete 10 second segments in data
x=reshape(map,[M1,fs*Ts]);
mapMn=mean(x,2);            % mean MAP in 10 second segments
y=reshape(icp,[M1,fs*Ts]);
icpMn=mean(y,2);            % mean ICP in 10 second segments

Now we compute PRx, using the 10-second-means of MAP and ICP:

M2=floor(M1/6)-4;           % number of PRx values that can be computed
PRx=zeros(1,M2);            % allocate vector for PRx
for i=1:M2
    PRx(i)=corr(mapMn(i*6-5:i*6+24),icpMn(i*6-5:i*6+24));
end
disp(PRx)

   1.0e-14 *

   -0.2615   -0.0775   -0.0350   -0.0257    0.0688    0.0385

There are 6 PRx values, as expected for 10 minutes of data. PRx(1) is derived from data from time 0 to 5 minutes. PRx(2) is from time 1 to 6 minutes. And so on. The PRx values are zero, as expected, since there is no correlation between the 10-second means of MAP and ICP, for this simulated data. In fact the 10-seond means are constant for this simulated data. But your data is probably different, so the correlations will not be zero at all times.

Good luck with your research.

[William Rose]

@Caleb,

In your original post, you said the results from PRx_display do not match results from PRx_calc, and you want to know why. The PRx values in PRx_display are read from a file, so we do not knowhow they have been calculated. The code and the raw data used to calculate the values in the file read by PRx_display would have to be supplied, if one wanted to do a detailed cmparison of PRx_display versus PRx_calc. And I do not want to do such a detailed comparison. But you could.

PRx_calc conmtans a significant amount of data pre-processing before it computes the correlations. For example: remove NaNs; compute a rounded time base; remove values that have the same time (on the rounded time scale); interpolate results to a common time base; etc. It is very possible that the differences in results are due to aspects of the data preprocessing.

I recommend that you generate simulated data that does not require any pre-processing or "clean-up". Comment out the pre-processing and clean-up code, which is not needed for this simulated data. See if the two algoritms give different results in this case. If they do give different results, you know where to focus. If they give the same results, then start adding back the pre-processing code, a little at a time, to see if and when the algorithms start to generate divergent results. If the algorithms never generate divergent results with the good simulated data, then make more simulated data, but this time with some of the flaws that the pre-processing is designed to fix, and repeat the process.

[William Rose]

@Caleb,

The simulated MAP and ICP data in my previous example had fluctiuation at the respiratory frequency (0.2 Hz). Since there were exactly 2 full breaths in each 10 second segment, the respiratory flucuations did not affect the 10-second mean values. Therefore the 10-sec means of ICP and MAP were constant, and therefore PRx=corr(mapMn,icpMn)=0.

Here is a different example, in which there is some correlation between MAP and ICP.

% Simulated data
fs=2;           % sampling rate (Hz)
dt=1/fs;        % sampling interval (s)
T=600;          % signal duration (s)
N=T/dt;         % number of samples
% Simulated mean arterial pressure (mmHg)=90 + noise with st.dev.=20
map=90+20*randn(1,N);
% Simulated intracranial pressure (mmHg)=15 + noise + portion correlated with MAP.
icp=15+3*randn(1,N)+0.05*(map-mean(map));

Compute the means of MAP and of ICP in 10-second segments:

Ts=10;                      % duration of segment for averaging (s)
M1=floor(N*dt/Ts);          % number of complete 10 second segments in data
x=reshape(map,[M1,fs*Ts]);
mapMn=mean(x,2);            % mean MAP in 10 second segments
y=reshape(icp,[M1,fs*Ts]);
icpMn=mean(y,2);            % mean ICP in 10 second segments

Compute PRx, using the 10-second-means of MAP and ICP:

M2=floor(M1/6)-4;           % number of PRx values that can be computed
PRx=zeros(1,M2);            % allocate vector for PRx
for i=1:M2
    PRx(i)=corr(mapMn(i*6-5:i*6+24),icpMn(i*6-5:i*6+24));
end

Display results:

disp(PRx)

    0.3149    0.2721    0.2370    0.3418    0.3870    0.5115

The six PRx values above are for 5-minute-long segments ending at 5, 6, ..., 10 minutes.

[William Rose]

@Caleb,

Thanks for posting data.

You have a data cleaning challenge. For example, consider the time values in the ABP file.

load ABP1.mat
t=time_vector;
dt=round(diff(t),3);  % vector of time intervals, rounded to nearest msec
fprintf('Tstart=%.1f, Tend=%.1f, duration=%.1f.\n',t(1),t(end),t(end)-t(1))

Tstart=827.9, Tend=250851.4, duration=250023.5.

Thus the data span almost 3 days (69.4 hours). Most of the sampling intervals are 1.024 or 1.025 (seconds, I assume). Too bad the sampling interval isn't 1.000, since PRx is defined using 30 consecutive, non-overlapping 10-second means. You could interpolate to 1 second sampling, or you could take 300 consecutive samples, which would be about 307 seconds instead of 300, or you could do something else.

There are 120 time intervals that are not 1.024 and not 1.025. This includes one interval=0, and intervals [249,320,1364,1682]. The latter intervals indicate gaps of about 20 to 30 minutes. How will you deal with these gaps? I suggest that you identify time ranges when there are at least 300 consecutive seconds with no gaps exceeding 1.5 s. The stem plot below shows the ending times of intervals that exceed 1.5 s.

[~,ind]=find(dt>1.5);   % indices of dt values >1.5
t_longdt=t(ind+1);      % end times of intervals with dt>1.5
stem(t_longdt,ones(size(t_longdt)),'r.'); hold on; 
stem([t(1),t(end)],[1,1],'g','filled'); hold off
ylim([0 1.5]); legend('dt>1.5','start,end')

Above: end times of intervals that exceed 1.5 s.

So far we have only discussed the time values. Next, consider the ABP values.

x=measurement_data;

There are no NaNs in this ABP recording, as showen by the result below:

disp(sum(isnan(x)))

     0

Plot ABP versus time:

plot(t,x,'-r'); xlabel('Time (s)'); ylabel('Pressure mmHg'); title('A.B.P.')

The plot of ABP versus time, above, shows numerous outliers. We need to remove them and do something about the missing values. Options include: 1. Determine if all ABP outliers correspond to the too-long intervals identified in the previous plot. If so, we may not be using those pressure values anyway, so their presence or absence won't be much problem. 2. If some outliers occur during otherwise good intervals, then we could (2A) interpolate across the outlier value(s), or (2B) add the times of outliers to the list of "bad" times from our analysis of the time vector, or (2C) do something else. We have not discussed how to identify the outliers. There are various possibilities.

Then you must repeat this entire process withe the ICP data file. Then you need to merge the ABP and ICP data sets to find the time ranges of 5 minutes or more where both ABP and ICP are good.

Then you are ready to compute PRx, in the good time ranges.

With all these options for data cleaning, and not knowing how the PRx values in the file were computed, it is hardly surprising to find differences betwen the PRx values in the file and the PRx values which you compute.

If you want to discuss this further, please send me a secure email by clicking on the envelope icon that apears in the window when you click onthe "WR" circle next to my posts. And if you do, please include your email, so I can respond directly.