I need help using the corrcoef function to calculate PRx (pressure reactivity index)
29 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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.
5 comentarios
Respuestas (1)
William Rose
el 12 de Ag. de 2024 a las 17:03
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)
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.
3 comentarios
William Rose
el 12 de Ag. de 2024 a las 20:54
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)
The six PRx values above are for 5-minute-long segments ending at 5, 6, ..., 10 minutes.
William Rose
el 13 de Ag. de 2024 a las 21:17
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))
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)))
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.
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!