Calculating the area under a curve, but above a certain threshold value

Question

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}}]

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}}] }
}]

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!

Answer 1

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}]

Answer 2

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@[email protected], EdgeForm@Thin, Polygon@data[n]}, Axes -> True]], {n, 0, Max@a}]