# Determining values of correction factor based on x bins in observed vs. actual data

I am trying to automate a problem I usually solve by hand.

I have a sensor that collects data from the field. Every 6 months or so, I have to do a calibration on that sensor by collecting x data points. Using a calibrated device, I measure the same data measured by the field sensor and I have the following table:

    Field  | Actual
#   spd val val   err
1    89 25  28  -10.71%
2   100 32  39  -17.95%
3    54 64  58   10.34%
4    65 32  34   -5.88%
5    70 12  11    9.09%
6   100 32  37  -13.51%
7   110 48  50   -4.00%
8    78 32  32    0.00%
9    95 18  22  -18.18%
10   92 82  81    1.23%
11   78 32  30    6.67%
12   98 38  40   -5.00%
13  103 42  45   -6.67%
14   99 54  57   -5.26%
15   88 61  57    7.02%
16   79 65  62    4.84%
17   83 58  54    7.41%
18  111 39  44  -11.36%
19  112 33  38  -13.16%
20   68 28  26    7.69%


The following is the Error versus Speed plot:

How I currently solve this problem:

• look at the data visually and determine how the data is "trending".
• I can split the data in x bins (depending on how new the sensor is, newer sensors can divide in to multiple bins)
• For the example above, I notice 3 bins:
• 50 to 70 - around 10% with one outlier
• 70 to 95 - around 3-4% with one outlier
• 95 to 120 - around -10% with one outlier.
• I then create the three bins for the sensor
• based on these bins, I determine a fudge or correction factor for each bin
• I then retest the sensor for more counts and redraw the graph. and repeat the process
• this process is time consuming especially when I have 10 sensors to look after, and that process can take up to 3 hours to complete

Question:

How can I create a statistical or mathematical solution to aid me in solving my problem?

Edit 1:

Some clarification:

• Speed: speed as the object hits the sensor
• Value: read by the sensor
• Final result is to obtain an average error of less than 6% on all reads
• Sometimes I collect more than 20 reads
• I would like to have uniformly calibrated across all ranges of speed as in the difference between actual vs measured to be relatively consistent so that every measured value could be used without significant error

Edit 2:

• tools I have access to are Excel, Python, and R (in that order of usability for myself)

Edit 3:

• Each object that's hitting the sensor is different in shape, size, and speed, so therefore the object hitting at any speed is not related to the size of the object. The sensor as it gets older tends to reads values differently at different speed bins.

Edit 4:

• Based on @Glen_b 's answer, I created the following MS Excel file for non binned results: http://97.107.136.148/results_stats_se.xls . Disregard the steer columns in the excel sheet. as it does not impact this problem.

Edit 5:

• I collected more data as shown below.
• Sub is a partial measurement of the object, so sub is part of the full_measurement.
• There is no ratio relationship between sub and full other than it's sub

Data:

Speed    Non Calibrated sensor        Calibrated sensor
Speed    Sub_1 Measurement_full    Sub_1    Measurement_full
sensor 1
88     5500      51800            5610    49010
95     5400      66000            5630    62700
sensor 2
99     4700      34600           43100    35850
103     5400      38200            5160    35360
84     4900      47700            4550    48330
101     5400      44800            5190    43050
106     5200      28400            5290    24390
96     4900      44900            4730    42900
104     5300      35000            5350    35290
102     5400      47700            5290    46810
99     5500      35500            5270    34410
89     5400      35600            5290    35490
98     5300      39300            5240    37580
sensor 3
93     5000      35600            4860    31050
79     5500      41000            5130    37810
92     4600      30700            4750    29720
99     4900      47100            4820    46230
94     5100      48300            4990    48230
93     5300      51600            5460    54850
99     5100      41700            4910    38960
85     5100      42200            4960    42810
81     5500      44100            5120    37960
98     5600      64300            5550    62510
102     5300      33400            5380    31220
sensor 4
106     4800      35000            4430    30440
99     5900      56300            5590    52740
78     5900      56300            5730    54380
96     6300      49100            5620    47790
98     5500      42100            5340    38340
70     5300      45700            4750    43670
77     5100      37900            5100    34470
106     4300      12100            3930     9750
107     5800      35600            5270    28930
90     3800      32300            3310    29390
75     5800      62200            5490    61470
88     4900      39200            5010    35720
sensor 5
104     4500      44200            5000    45120
95     5460      49400            5440    48520
104     4700      33900            5210    33420
86     5200      36400            5280    35140
88     5000      49500            5760    52250
94     5300      44400            5200    43310
95     5200      54000            5370    54460
99     5100      37000            5410    35600
91     5600      62400            5540    61860
95     5400      38000            5310    32510
94     5100      49300            4990    46210
101     4800      44300            5380    45130
100     6600      23700            7500    24630
105     5100      29200            5400    26740
94     5500      42400            5730    40760
89     5200      49600            5270    49570
99     5300      47100            5700    47970
100     4900      36600            5220    36790
90     7200      25500            7340    23480
100     5200      36500            5410    36210

-
Why are you binning the data rather than say using a regression to come up with a calibration line/curve? –  Glen_b Jun 18 '13 at 2:47
@Glen_b I reflected on some of our discussion in the clarification. Perhaps, we should start deleting some of our comments –  dassouki Jun 18 '13 at 3:23
Where you say 'final result' -- you can't guarantee this; if your sensor has a lot of noise in its measurement, you may not even get close. You can't manufacture information that's not present. Presumably what you mean is that you'd like to have that, and if you have it you're satisfied and don't need to keep trying to improve the situation. –  Glen_b Jun 18 '13 at 4:36
Re Edit 4: what do the "Steer" columns mean? –  whuber Jun 18 '13 at 15:03
@whuber I forgot to delete the steer columns. In some standards we are asked to calibrate by the steer value, in others we do it by the reading. For my sensors, calibration is done based on the reading variable. I will reflect that in an edit –  dassouki Jun 18 '13 at 15:05

It looks to me like it might potentially make sense to account for the way the error changes with the second (calibrated) val as well as with speed. If you're only interested in the relationship with speed, it looks like linear regression would be a good first suggestion.

First, let's look at your data; this is every variable against every other variable, usually called either a scatterplot matrix or a pairs plot:

Features in the pairs plot:

The plot of err vs spd, bottom left, outlined in reddish brown. This is the same plot you show and indicates a clear downward trend. This is what prompts me to suggest linear regression of some kind as a 'first attempt'.

Consider the plot of val1 against val2 (val2 is the calibrated-instrument measurement, the one you want your field measured-val to be close to), the plot near the middle, outlined in red. It has a nice straight-looking relationship, as you'd hope. There's a 'flat' part where a whole bunch of different calibrated values near 32 were recorded as '32'; this potentially suggests some kind of issue with the device (though it could also be mere chance combined with the rounding that's there).

Consider the err vs val2, the plot in blue. Values to the left of $\text{val2}\leq 32$ (roughly) look very noisy there and then there's less noise but a distinct upward trend after that. You can see a similar effect in the plot of err vs val1 (the purple-outlined plot in the last row; this is potentially useful, because val2 isn't available normally).

Edit: there's also a suggestion of nonlinearity in the relationship of the two 'val' variables with spd (unsurprisingly).

---

Linear regression (first model):

This fits a straight line. You could just fit a straight line of err against spd, but this misses a bunch of potential value in the information you have in the other variables.

The usual form of linear regression is via least squares, and this is probably easiest to do, but given the possibility of outliers you mention, you may want to consider a more robust alternative.

Here's a linear regression done in R. Note that my err variable isn't multiplied by 100, so if you did it with your numbers, your coefficient would be 100 times as big:

 mdl1 <- lm(err~spd,calib)
summary(mdl1)

Call:
lm(formula = err ~ spd, data = calib)

Residuals:
Min       1Q   Median       3Q      Max
-0.12761 -0.02540  0.01683  0.04742  0.09650

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.3242267  0.0865953   3.744 0.001485 **
spd         -0.0039833  0.0009619  -4.141 0.000613 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06867 on 18 degrees of freedom
Multiple R-squared:  0.4879,    Adjusted R-squared:  0.4594
F-statistic: 17.15 on 1 and 18 DF,  p-value: 0.0006134


The residual standard error, 0.0687 suggests that the typical error (let alone the maximum error!), even after adjustment is still a good deal larger than you'd like.

Looking at the diagnostics, the model doesn't look too bad. There are low outliers, and the variation looks bigger near the middle, with a decrease at larger fitted values.

---

The relationship between val2 and val1

Before we plow on with that, though, we should investigate the relationship between the field measured val and the calibrated one. The pairs plot above doesn't make it obvious, but I could discern a slight bend... and given the pattern in the err vs val2 plot, I decided to take a closer look:

Clearly there's a 'bend' here! You shouldn't ignore that; it's a source of a fair amount of your error.

It's clearer still if you look at the raw error:

Edit: This effect remains if you remove the clump of values at 32 --

The problem is, of course that in the field you don't know val2 and you need val1 to estimate the 'gold standard' value -- yet val1 is the 'noisy' measurement. This is sometimes called an 'inverse regression' problem.

It's probably worth your time to google on this; or rather on one of these:

inverse regression -sliced --- (because this avoids some not-so-relevant hits)

inverse regression calibration

Your problem is a calibration problem (you're trying to figure out how to estimate val2, the 'accurate' measure, from val1, the noisy one), but you have the advantage of a second predictor, which seems to add useful information to the raw calibration curve.

---

A second, simple regression model

With the decrease in variability against fitted in the first model, and the fact that the actual problem is to try to 'back out' val2 from this, consider modelling the actual error rather than the percentage error.

Further, while this model will ignore that 'bend' in the data, it will consider val2 in the model for val1. Well, we will fit two models in order to see if there's value in extending it:

summary(m1)

Call:
lm(formula = I(val1 - val2) ~ spd, data = calib)

Residuals:
Min      1Q  Median      3Q     Max
-4.8978 -1.3288 -0.3747  1.5850  4.8034

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.35791    3.22609   4.141 0.000614 ***
spd         -0.16092    0.03583  -4.491 0.000283 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.558 on 18 degrees of freedom
Multiple R-squared:  0.5284,    Adjusted R-squared:  0.5022
F-statistic: 20.17 on 1 and 18 DF,  p-value: 0.0002827


We can't really compare this model with the previous one, it's not on the same scale (however, we again see that spd is needed to 'adjust'). Now lets add val2:

summary(m2)

Call:
lm(formula = I(val1 - val2) ~ spd + val2, data = calib)

Residuals:
Min      1Q  Median      3Q     Max
-4.4331 -0.8054  0.0845  1.0342  3.4050

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.55003    2.77761   3.798  0.00144 **
spd         -0.17419    0.02954  -5.897 1.76e-05 ***
val2         0.09427    0.02976   3.168  0.00562 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.088 on 17 degrees of freedom
Multiple R-squared:  0.7034,    Adjusted R-squared:  0.6686
F-statistic: 20.16 on 2 and 17 DF,  p-value: 3.258e-05


We see that the largest and smallest residuals are smaller and that spd has a very similar effect to before. The residual standard error is about 20% smaller and the val2 effect is significant. That is, if we want to try to 'guess' val2 from spd and the noisy val1, we should not ignore the linear relationship between val1 and val2. The diagnostics (not shown) are better, though the presence of low outliers does suggest that a more robust method might get a little closer for the bulk of the data.

If you do this with the relative error instead, you don't see as much value in adding val2, but I'd still include it (and the model diagnostics aren't nearly as nice, because the size of the error really isn't constant on the percentage scale; it's much nearer to constant on the scale of actual error).

It doesn't get you down to 6% error, but that's mostly because the percentage errors for the low vals are always going to be large on the percentage scale, even though they're nice looked at as raw errors), but if you want to lower it further, you'll need more data in order to estimate a calibration curve (and your job will become more complex) - there's a suggestion of curvature with both val2 and spd; I'd probably look at natural spline regression if you went that direction.

Using Excel:

The linear regressions can be done in Excel, but some of the diagnostics won't be available, and if come to trying to calibrate a curve, it becomes very unwieldy. I think it's worth putting the data into R and learning the tools.

Some issues with sample splitting:

Take a look at some of the discussion here. In the links there, it's discussing the idea of splitting samples into either three or two groups instead of fitting a line, where the actual relationship is smooth (doesn't seem to have 'jumps'). If some of the terms aren't clear in the links, please ask.

Justification for low/high or tertiary splits in ANOVA

(However, if your example data is typical, you may actually have jumps. In that case, you might want to consider models that pick up such jumps without requiring you to identify them by hand - and you would want to account for the potential overfitting as well)

-
Thanks for the comprehensive answer. Based on my first couple of readings, I have a couple of questions for you: (1) which models would you suggest I read up on that deals with jumps? (2) you're suggesting to always use one speed bin regardless of what the data is showing? –  dassouki Jun 18 '13 at 11:46
As an addition, it's important to note that each object hitting the sensor is different. Some objects are bigger or smaller than others just like they could be faster as well –  dassouki Jun 18 '13 at 12:45
Also as an FYI, I'm trying to solve the problem using the inverse regression technique you pointed out –  dassouki Jun 18 '13 at 12:50
+1. But it's not so clear there is a "bend" in the val1 - val2 regression: another interpretation is that val1 somehow gets "stuck" at a value of $32$ or so (look at the scatterplot matrix closely). That suggests an entirely different intervention: fix the sensor. –  whuber Jun 18 '13 at 13:18
Thanks to you Glen and @whuber for the valuable discussion on this question. I must say that both of your posts comments and discussion served as an eye opener to me and simplified numerous statistical facts that were previously hard to understand (for me). –  dassouki Jun 28 '13 at 14:04