5
$\begingroup$

so here I am with a time series of data (hours (t) and corresponding measurements (a)).

a = {0, 2, 5, 6, 3, 6, 5, 8, 2, 1, 10};
t = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10};

The plot looks like this:

ListLinePlot[Partition[Riffle[timepoints, a], 2],  GridLines -> {{}, {4}}]

enter image description here

Using a trapezoidal rule function, I get an AUC of 43.

auc[data_, times_] := (
merged = Inner[List, times, data, List];
merged = DeleteCases[merged, {x_, y_} /; Not[NumericQ[y]]]; 
N@Total[Partition[merged, 2, 1] /. {{x1_, y1_}, {x2_, y2_}} -> (x2 - x1) ((y1 + y2)/2)] )

What I am trying to do now is modify the function to get the AUC for values above a certain threshold, let's say 4. This is indicated by the horizontal grid line in the graph.

First thought was replacing the values exceeding 4 by 4, calculating a second AUC and substracting. Of course this changes the shape of the curve.

a2  = a /. x_ /; x > 4 -> 4;
GraphicsGrid[{
{ListLinePlot[Partition[Riffle[timepoints, a], 2], 
  PlotRange -> {All, {0, 10}}, GridLines -> {{}, {4}}],
ListLinePlot[Partition[Riffle[timepoints, a2], 2], 
  PlotRange -> {All, {0, 10}}, GridLines -> {{}, {4}}] }
}]

enter image description here

I figure the solution is along the lines of linear interpolation of where the curve intersects the threshold and adding intermediate points and I really have no idea where to start implementing this.

Any pointers? Or are there biomedical packages that will calculate this stuff from raw data?

Thx!

$\endgroup$
4
$\begingroup$
a = {0, 2, 5, 6, 3, 6, 5, 8, 2, 1, 10};
t = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
f = Interpolation[Transpose@{t, a}, InterpolationOrder -> 1];
Plot[f[t], {t, 0, 10}, ImageSize -> 200]
f2[t_] := Min[f[t], 4];
Plot[f2[t], {t, 0, 10}, ImageSize -> 200]
NIntegrate[f[t], {t, 0, 10}]
NIntegrate[f2[t], {t, 0, 10}]

enter image description here

$\endgroup$
  • $\begingroup$ Chenminqi, you just solved a whole lot of problems at once. Thanks! $\endgroup$ – blinkenstack Jun 21 '14 at 14:32
1
$\begingroup$

Here is a solution which would work only for ListPlot. PolygonArea taken from here, inter taken from here.

Graphics`Mesh`MeshInit[];
a = {0, 2, 5, 6, 3, 6, 5, 8, 2, 1, 10};
t = Range[0, 10];
raw = Reverse /@ Thread[{a, t}];
(* find the intersection between the threshold and the ListLinePlot *)
f = Interpolation[raw, InterpolationOrder -> 1];
inter[n_] := 
  x /. FindRoot[f[x] - n, {x, #[[1, 1]], #[[2, 1]]}] & /@ 
    Select[Partition[Sort@Last@
      Last@Reap[Plot[f[x] - n, {x, 0, 10}, EvaluationMonitor :> Sow[{x, f[x] - n}]]], 2, 1], 
      #[[2, 2]] #[[1, 2]] <= 0 &]
coord[n_] := Thread[{inter[n], Array[n &, Length@inter[n]]}]
(* Plot it and calculate the area *)
data[n_] := 
  DeleteDuplicates@
    Append[Sort@Join[Select[raw, Last@# <= n &], coord[n], {{Max@raw, n}}], {Max@raw, 0}]
Manipulate[
 Show[ListLinePlot@raw, 
 Graphics[{Text["Area = " <> ToString@Abs@PolygonArea[data[n]], {Mean[First /@ raw], 2}], 
   FaceForm@Opacity@.1, EdgeForm@Thin, Polygon@data[n]}, Axes -> True]], {n, 0, Max@a}]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.