MATLAB

Ambulatory Glucose Profile – MATLAB tutorial


This post is meant to give a very short demonstration on how to generate a AGP-plot from measured CGM data. AGP stands for ambulatory glucose profile and is basically a visual condensed representation of multiple days of continuous glucose measurements. It allows to visually examine time ranges that are problematic and require adoption of the therapy. It does so by showing how much the glucose levels fluctuate over the course of a day. Therefore, data of multiple days is used to generate inter quartile range plots. Since I did not find too much information about how to generate these plots, i thought a quick and dirty tutorial might be helpful. For more information look here.

So let’s start!

At this point, we assume that the CGM data has already been preprocessed in a ways that each day is reduced to 24 hourly average values (x = hour, y=glucose value). Given data of 14 consecutive days (assuming no missing values), we have a total of 14*24 = 288 glucose readings.  The sample code below creates a random dataset with  a peak in placed around 6 AM (hyperglycemia in the morning) and a minimum at around 6 PM (hypoglycemia in the late afternoon).

x = 0:1:24';
y = 100+100*rand(14,25).*repmat(sin(x*2*pi/25)+0.7,14,1);
h = plot(x,y); hold on;
xlabel('hour of day');
ylabel('glucose concentration in mg/dL');
xlim([0,24]);

We can plot the mean glucose profile by calculusing the average for each day hour of all fourteen days. To do so, we use the mean function of MATLAB:

havg = plot(x,mean(y),'r','lineWidth',3);

Now, the 25/75%-IQR inter quartile ranges shall be added to the graph. This range covers 50% of the most common glucose values and gives a good indicator about how much the glucose values fluctuate. To show the upper and lower bounds, we have to calculate the 25 and 75 percentiles for each hour of the day.

p25 = prctile(y,25);
p75 = prctile(y,75);
h25 = plot(x,p25,'b--','lineWidth',3);
h75 = plot(x,p75,'b:','lineWidth',3);
legend([havg,h25,h75],'mean','25 percentile','75 percentile');

We want to fill the area between the lower and upper percentiles to better visualize the range. We can do this using Matlab’s fill function. The function takes a polygon in the form of x,y-coordinates and fills the area using a solid color. From both percentile curves, we create a polygon by simply appending the lower curve in reverse order to the upper curve. The x values must therefore also be reversed using the fliplr function.

h2 = fill([x,fliplr(x)],[p75,fliplr(p25)],'b','');
set(h,'facealpha',0.7);

The edgy appearance is due to the small number of intermediate points (we only have 24 values per curve). If we would generate more support points like for every minute of the day or so, the problem would not be so dramatic and eventually could be ignored. Here, the shape looks very unsatisfying and should be visually improved. Therefore, we create a smooth transition between the intermediate points using spline interpolation utilizing MATLAB’s interp1 function. It takes the original x and y values at first and then the new x positions for which the interpolated values shall be generated. The last argument indicates the method used for generating intermediate points.
interp1(x,y,xnew,'spline');
Thus, we first have to create the interpolated curves using 200 intermediate points which is enough to create a smooth appearance. Then we repeat the above procedure to fill the area between them. The figure below shows the smoothed results.

xInterp = linspace(0,24,100);
p25Interp = smooth(x, p25, xInterp);
p75Interp = smooth(x, p75, xInterp);
h2 = fill([xInterp, fliplr(xInterp)], [p25Interp, fliplr(p75Interp)],'b');
set(h2,'facealpha',0.7);

We can add another inter quartile range. Traditionally, this would be the 10/90%-IQR.

p10 = prctile(y,10);
p90 = prctile(y,90);
p10Interp = interp1(x,p10,xOS,'spline');
p90Interp = interp1(x,p90,xOS,'spline');
h3 = fill([xOS,fliplr(xOS)],[p90OS,fliplr(p10OS)],[0,0,1]);

The average value should also be smoothed and added to the plot at the very end so that it is not covered by the filled polygons. Note that we used a smaller alpha value for the second IQR to distinguish between them. The overall result and code is shown below.


 

%% create random sample data: hourly average glucose values for 14 days
x = 0:1:24';
y = 100+100*rand(14,25).*repmat(sin(x*2*pi/25)+0.7,14,1);

%draw data and average curve
h = plot(x,y); hold on;
xlabel('hour of day');
ylabel('glucose concentration in mg/dL');
xlim([0,24]);
havg = plot(x,mean(y),'r','lineWidth',3);
pause();

%% draw percentile curves
p25 = prctile(y,25);
p75 = prctile(y,75);
h25 = plot(x,p25,'b--','lineWidth',3);
h75 = plot(x,p75,'b:','lineWidth',3);
legend([havg,h25,h75],'mean','25 percentile','75 percentile');
pause();

%% fill 25/75 IQR
hIQR = fill([x,fliplr(x)],[p75,fliplr(p25)],[0,0,1]);
set(hIQR,'facealpha',0.7);
pause();

%% fill smoothed version of 25/75 IQR
%first, we remove the old plot
delete(h25);
delete(h75);
delete(hIQR);
%smooth using spline interpolation
xInterp = linspace(0,24,200);
p25Interp = interp1(x,p25,xInterp,'spline');
p75Interp = interp1(x,p75,xInterp,'spline');
%draw
hIQR = fill([xInterp,fliplr(xInterp)],[p75Interp,fliplr(p25Interp)],[0,0,1]);
set(hIQR,'facealpha',0.7);
legend([havg,hIQR],'mean','25/75-IQR');
pause();

%% add a smoothed 10/90-IQR
p10 = prctile(y,10);
p90 = prctile(y,90);
p10Interp = interp1(x,p10,xInterp,'spline');
p90Interp = interp1(x,p90,xInterp,'spline');
hIQR2 = fill([xInterp,fliplr(xInterp)],[p90Interp,fliplr(p10Interp)],[0,0,1]);
set(hIQR2,'facealpha',0.3);
legend([havg,hIQR,hIQR2],'mean','25/75-IQR','10/90-IQR');
pause();

%% add smoothed mean curve
yavgInterp = interp1(x,mean(y),xInterp,'spline');
havg = plot(xInterp,yavgInterp,'r','lineWidth',2);
legend([havg,hIQR,hIQR2],'mean','25/75-IQR','10/90-IQR');

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s