### presentation

```CONSENSUS MULTIVARIATE
CALIBRATION OR MAINTENANCE
WITHOUT REFERENCE SAMPLES
USING TIKHONOV TYPE
REGULARIZATION APPROACHES
John Kalivas, Josh Ottaway, Jeremy Farrell, Parviz Shahbazikah
Department of Chemistry
Idaho State University
Pocatello, Idaho 83209 USA
Outline
• Multivariate calibration
• Tikhonov regularization (TR)
• TR calibration maintenance with reference samples to
form full wavelength or sparse models
– Selecting “a” model
– Selecting a collection of models
– Comparison to PLS
• TR calibration or maintenance without reference
samples
– Examples with comparison to PLS
• Summary TR variant equations
2
Spectral Multivariate Calibration
• y = Xb
y = m x 1 vector of analyte reference values for m calibration
samples
X = m x n matrix of spectra for n wavelengths
b = n x 1 regression (model) vector
•
ˆ   XT X 
b
1
T
X y
MLR solution; requires m ≥ p (wavelength selection)
• bˆ  X  y
Biased regression solutions such as TR, RR (a TR variant),
PLS, and PCR
• Requires meta-parameter (tuning parameter) selection
T
• yˆ u n k  x u n k bˆ
3
Quantitation by Tikhonov Regularization (TR)
• General TR in 2-norm
m in
 Xb  y
2
2
η
2
2
Lb
2
y X 
 
b
 0   ηL 
t
2 t
bˆ   X X  η L L 
1
= Euclidian vector 2-norm
(vector magnitude or
length)

2
Depending on the
calibration goal, L can
have different forms
t
X y
• Ridge regression (RR) when L = I
m in

Xb  y
2
2
η
2
b
2
2
y X 
    b
 0   ηI 
t
2
bˆ   X X  η I 
1
X y
t

RR is regularized by using I
and selecting η to minimize
prediction errors (low bias)
simultaneously shrinking the
model vector (low variance)
4
Selecting η
• Cross-validation
• L-curve graphic (can use with RMSEC)
• Bias/Variance can be assessed
• Useful for putting RR, PLS, etc. on one plot for objective
comparison
– C.L. Lawson, et.al.,
Solving Least-Squares
Problems. Prentice-Hall,
(1974)
– P. C. Hansen,
L bˆ
Rank-Deficient and
Discrete Ill-Posed
Problems: Numerical
Aspects of Linear
Inversion, SIAM Press
(1998)
overfitting
best model
2
underfitting
yˆ  y
2
5
Calibration Maintenance
• Need primary model to function over time
and/or under new secondary conditions
1. Prepare calibration samples to span all potential
spectral variances
•
Not possible with a seasonal or geographical effects in
some data sets
2. Preprocess primary and secondary data to be
robust to new conditions
3. Adjust spectra measured under new conditions to
fit the primary model
4. Update the primary model to predict in the new
conditions
6
Calibration Maintenance with TR2
• Model updating a RR model requires a new penalty term
– Minimize prediction errors for a few samples from new secondary
conditions

m in X b  y
2
2
η b
2
2
2
 λ Mb  yM
2
2
2

M = spectra from secondary conditions
yM = analyte reference values
• Avoid measuring many samples by tuning with λ
• Local centering
– Respectively mean center X, y, M, and yM
– Validation spectra centered to M
7
Pharmaceutical Example
• M. Dyrby, et. al., Appl. Spectrosc. 56 (2002) 579-585
– http://www.models.life.ku.dk/datasets; Dept. of Food Sciences, Univ. of
Copenhagen
• 310 Escitolopram tablets measured in NIR from 7,400-10,507 cm-1 at
resolution 6 cm-1 for 404 wavelengths
• Four tablet types based on nominal weight: type 1, type 2, type 3, and
type 4
• Three tablet batches (production scale): laboratory, pilot, and full
• 30 tablets for each batch tablet type combination
0.15
 Lab, type1
 Lab, type 2
 Lab, type 3
 Lab, type 4
 Full, type 1
 Full, type 2
 Full, type 3
 Full, type 4
0.1
PC2
0.05
0
-0.05
-0.1
-0.1
0
PC1
0.1
8
Objective
• Using laboratory produced tablets as the primary
calibration set
– Determine active pharmaceutical ingredient (API)
concentration in new tablets produced in full production
(secondary condition)
Primary Calibration Space: 30 random lab batch
samples with 15 from types 1 and 2 each
Secondary Calibration Space: 30 random full batch
samples with 15 from types 1 and 2 each
Standardization Set M: 4 random full batch samples with
2 from types 1 and 2 each
Validation Space: Remaining 30 full batch types 1 and 2
• Other batch type combinations studied
9
Example Model Merit Landscapes

m in X b  y
2
2
η b
2
2
2
 λ M b  yM
2
2
2

η
RMSEC
η
λ
λ
η
RMSEM
η
λ
λ
10
Model Merit Landscapes

m in X b  y
2
η b
2
2
2
 λ M b  yM
2
2
2
η
η
RMSEC
2
λ
RMSEM
λ
Convergence at small λ
Prediction of primary
prediction of
secondary improves

• Secondary conditions are not
included in new model
• Amounts to using primary RR with
local centering where secondary
validation samples are centered to
the mean of M
Best local
centered models
11
Model Merit Landscapes

m in X b  y
η b
2
2
2
2
 λ M b  yM
2
2
2

η
η
RMSEC
2
bˆ
λ
λ
bˆ
η
2
RMSEM
2
Prediction of primary
prediction of
secondary improves
too large
and RMSEC and RMSEM
• Can use an L-curve at a
fixed λ value
λ
12
2
Model Merit Evaluations
• Multiple merits can be used to assess tradeoff
– Respective RMSEC and RMSEM landscapes for R2, slope,
and intercept
– L-curves at selected η and λ values
2
ˆb
– H i   RM SEM i RM SEM m ax  2  bˆ i
m ax

2
2

1
H
0.8
λ = 54.29
η
RMSEV
0.6
0.4
0.2
0
0
0.004 0.257 16.037 1000
η
λ
λ = 54.29
13
Model Updating Results
Updating Primary Lab Batch Types 1 and 2 to Predict
Secondary Full Batch Types 1 and 2
bˆ
Method
2
TR2
2.10
RR local
centering
RR no
update
η
λ
RMSEC
RMSEM
RMSEV
R2
0.731
0.014
0.264
0.966
2.70
0.468
0.245
0.487
0.935
16.40
0.096
-
0.653
0.925
1.701
54.29
0.588
0
0.0359
-
Lab and Full Batches Types 1 and 2 Self
Predicting Using RR
Batch
bˆ
Lab
Full
RMSEC
RMSEV
R2
η
16.40
0.096
0.276
0.972
0.0359
1.34
0.205
0.239
0.968
0.215
2
• Updated primary
models predicts
equivalently to the
secondary model
predicting the
secondary
validation samples
14
Model Vectors
RR Lab Batch
3
2
bˆ i
TR2
1
0
-1
-2
0.2
8000
9000
Wavelength,
10000
cm-1
RR Full Batch
0
-0.2
0.1
bˆ i
bˆ i
-0.4
0
8000
9000
10000
Wavelength, cm-1
bˆ
-0.1
2
-0.2
8000
9000
10000
Wavelength, cm-1
15
Using PLS
• PLS (and other methods) can also be used
• With PLS, the PLS latent vectors (PLS factors) replace
the η values
TR2
 y

0

 λy
M

  X 
 

 ηI b
 

  λM 
 

1
t
2
2
t
t
ˆ
b  X X  η I  λ M M  X y
PLS
 y   X 


b
 λyM   λM 
y '  X 'b

bˆ   X '  y '
16
PLS Model Merit Landscapes
RMSEM
RMSEC
1
0.8
0.5
5
5
Factors
1
0.6
0.4
15
0.2
20
0
0.4
10
10
0.002
λ
1.0
0.3
15
0.2
20
0
37
0.1
0.002
bˆ
1.0
37
RMSEV
2
1
1
10
60
10
15
40
15
20
20
20
Factors
5
80
0
λ
0.002
λ
1.0
37
0.6
5
0
• Similar
landscape trends
• The discrete
factor aspect of
PLS can make it
difficult to
capture the
underlying
continuity of the
landscapes
0.5
0.4
0.3
0.002
λ
1.0
37
17
PLS and TR2 Model Updating Results
Updating Primary Lab Batch Types 1 and 2 to
Predict Secondary Full Batch Types 1 and 2
Method
bˆ
2
RMSEC RMSEM RMSEV
R2
TR2
2.10
0.731
0.014
0.264
0.966
PLS
3.29
0.658
0.024
0.266
0.964
η
λ
1.70
54.29
3
factors
19.31
• PLS prediction equivalent to TR2
• The discrete factor aspect of PLS can make it difficult to
capture the underlying continuity of the landscapes
18
Sparse TR Calibration Maintenance
• TR2:
m in
 Xb  y
2
2
η
2
b
2
2
 λ M b  yM
2
2
2

• TR2b (sparse model):

m in
Xb  y
2
2
η
2
Lb
2
2
2
 λ M b  yM
2
2

• L = diagonal matrix with lii = 1/│bi│
Gorodnitsky IF, Rao BD. IEEE Transactions on Signal Processing 1997; 45:
600-616.
• TR2-1 (sparse model):
m in
 Xb  y
2
2
η b
 λ M b  yM
2
1
2
2

19
TR2-1 Sparse Model Merit Landscapes

m in X b  y
Models with increasing η
RMSEC
bˆ
2
λ
2
2
η b
 λ M b  yM
2
1
RMSEM
RMSEV
λ
• Similar landscape
trends
• For small λ values,
the η values are the
same across λ
• At greater λ values,
the η values vary
across λ
20
2
2

TR2 and TR2-1 Model Updating Results
Updating Primary Lab Batch Types 1 and 2 to Predict
Secondary Full Batch Types 1 and 2
Method
bˆ
R2
RMSEC RMSEM RMSEV
2
TR2-1
10.65
0.666
0.029
0.227
0.970
TR2
2.10
0.731
0.014
0.264
0.966
PLS
3.29
0.658
0.024
0.266
0.964
η
λ
7650
2442
1.70
54.29
3
factors
19.31
TR2-1 prediction
results improve
over TR2 and PLS
0.5
TR2-1
5
bˆ i
TR2
0.2
0
PLS
0
0
-0.2
-5
-0.5
-0.4
8000
9000
cm-1
10000
8000
9000
cm-1
10000
8000
9000
cm-1
10000
21
Other TR Sparse Maintenance Methods
• TR2-1b (sparse models when L = I or L ≠ I):
m in

2
Xb  y
2
 η Lb
2
 λ M b  yM
2
1
2

• TR1-2b (full wavelength when L = I):
m in

Xb  y
2
2
η
2
Lb
2
2
 λ M b  yM
1

• TR1 (sparse models when L = I or L ≠ I):
m in

Xb  y
2
2
 η Lb
1
 λ M b  yM
1

22
Other TR Applications
• Updating a primary model:
– for extra virgin olive oil adulterant quantitation to a
new geographical region (applicable to new
seasons)
– to a new temperature
– formed on one instrument to work on another
23
Summary
• Only a few samples needed for M with appropriate
weighting
• Same samples measured in primary and
secondary conditions are not needed
– Avoids long term stability issue
• PLS and other methods can be used
– Discrete nature (PLS factors) can limit landscapes
• Need to select a pair of tuning parameters for “a”
model
• Requires reference values for yM
24
Consensus (Ensemble) Modeling
• Samples predicted with a collection of models
– Composite (fused) prediction is formed
– Simple mean prediction used here
• Typically form models by random sampling
across calibration samples and/or variables
• From collection, filter for model quality
• Ideal models:
1. High degree of prediction accuracy
2. Small but noteworthy difference between selected
models (model diversity)
25
Consensus TR and PLS Modeling
• Models formed from varying tuning
parameter values
• Plot predicted values against reference
values for X,y and M,yM
• Use respective R2, slope, and intercept
model merit values
• Natural target values:
–
→1
– Slope → 1
– Intercept → 0
R2
Merit
Min
Max
0.85
0.95
Slope X,y 0.85
0.95
R2 X,y
Intercept
0.30
1.00
0.98
0.99
0.01
0.20
0.95
0.99
X,y
R2 M,yM
Intercept
M,yM
Slope
M,yM
26
TR and PLS Consensus Models (RMSEV)
1 PLS Model
348 TR2 Models
1
0.6
0.5
Factors
η
5
10
0.4
15
20
λ
Model with increasi8ng η
628 TR2-1 Models
λ
0
0.3
0.002
1.0
37
• Fewer PLS models selected
due to sharpness of
landscapes from the discrete
factor nature of PLS
• Number of “good” models can
reducing the increment sizes
of η and λ
27
Consensus Mean Model Updating Results
Updating Primary Lab Batch Types 1 and 2 to Predict
Secondary Full Batch Types 1 and 2
No.
Method
Models
bˆ
2
RMSEC RMSEM RMSEV
R2
TR2
348
4.32
0.552
0.016
0.284
0.955
TR2-1
628
16.32
0.580
0.007
0.274
0.958
PLS
1
3.29
0.658
0.024
0.266
0.964
η
Λ
0.591
207
0.379
1619
3 factors
19.31
• The one PLS model predicts best
• PLS limited to discrete factors where TR
allows 0 ≤ η < ∞ to more fully resolve the
landscape
28
Consensus Models and Correlations
2.5
TR2
348 models
1.5
i
0.5
0
4
2
Frequency
bˆ 1
x 10
1.5
1
0.5
-0.5
8000
9000
0
10000
0.6
x 10
TR2-1
628 models
bˆ i
10
0
-10
8000
9000
cm-1
10000
0.9
1
4
6
Frequency
20
0.7
0.8
Correlation
4
2
0
0
0.5
Correlation
1
29
Summary
• Only a few samples needed for M with appropriate
weighting
• Same samples measured in primary and secondary
conditions are not needed
– Avoids long term stability issue
• Can select “a” model or a collection of models
– Natural target values (thresholds) with model merits R2, slope, and
intercept for primary and secondary standardization sets
– Work in progress
• Requires reference values for yM
30
Without Reference Samples
Beer’s law: x = yaka + yiki + m +
+n
ka = pure component (PC) analyte spectrum
ki = PC spectrum of ith interferent
(drift, background, etc.)
m = rest of the sample matrix
n = spectral noise
T
yˆ  x bˆ





T
T
T
T
yˆ  y a k a bˆ  y i k i bˆ  m b  ...  n bˆ
• Ideal situation:
WHEN: 1 .  k Ta bˆ   1


T
T
2. k i bˆ an d m bˆ  0
3. bˆ
T
su ch th at n bˆ  0
2
THEN: 4. yˆ  y a
• Cannot simultaneously satisfy 1, 2, and 3 to obtain 4
31
Compromise PCTR2 Model

m in k b  1
T
a
2
2
η b
2
2
2
 λ Nb  0
2
2
2

N = spectra without analyte, e.g., ki
• Minimizing the sum requires a tradeoff between
the three conditions
– The closer the three conditions are met, the more
likely yˆ  y a
• Updating the non-matrix effected PC ka to predict
in current conditions (spanned by N)
32
Sources of N
• PC interferent spectra
– Reference values are 0
• Matrix effected samples without the analyte
– Reference values are 0
• Constant analyte samples
– Reference values are 0 after spectra are mean centered
• Estimate using samples with reference values
T

yy 
N  I T X
y y

Goicoechea et al. Chemom. Intell. Lab. Syst.
56 (2001) 73-81
• Samples for N need to be measured at current
conditions
33
• EVOO samples: Crete, Peloponnese, and Zakynthos
–
–
–
–
RR calibration: 56 samples spiked 5, 10, and 15% (wt/wt) sunflower oil
Primary: PC sunflower oil, 1 sample
Secondary: EVOO, 25 samples
Validation: 22 spiked samples
• Synchronous fluorescence
spectra 270 to 340 nm
at Δλ=20 nm
3
x 10
6
Zakynthos
Intensity
Sunflower
EVOO
2
1
0
280
300
320
Excitation Wavelength (nm)
34
Model Merit Landscapes
RMSEN
RMSEPC
λ
η
bˆ
η
λ
RMSEV
2
η
λ
η
λ
35
H Values
PCTR2 at η = 9.1e3
RR with PCTR2 cal samples
1
1
0.8
0.8
H
0.6
0.6
0.4
0.4
0.2
0.2
1
3.6e-6
3.6e-3
λ
3.6
1.8e3
0
95
9.7e4
1.0e8
η
RR full cal
1
H
0.8
H i   RM SEN i RM SEN m ax  
2
0.6

H i   R M SEC i R M SEC m ax  
2
0.4
bˆ i

bˆ m ax
2
bˆ i
2
2
bˆ m ax

2
2

2
0.2
95
9.7e4
η
1.0e8
36
Model Updating From PC Sunflower
Method
(No.
Samples)
bˆ
PCTR2 (26)
RMSEV
R2
2.6e-7
0.031
0.882
RR (56)
4.0e-7
0.028
0.649
RR with
PCTR2
samples (26)
2.3e-7
0.077
0.787
0.2
0.15
yˆ i
0.1
2
RR
PCTR
LS line
LS line
Equality
η
λ
9.1e3
0.0036
1.9e5
-
• Updated PC predicts better
than a full calibration
1.6e6
-
yi = 0.422xi + 0.048
yi = 0.807xi - 0.0074
bˆ
2
0.05
0
0.05
0.1
yi
0.15
37
Temperature Data Set
• Wülfert, et al., Anal. Chem. 70 (1998) 1761-1767
• hhttp://www.models.life.ku.dk/datasets; Dept. of Food Sciences,
Univ. of Copenhagen
•
•
•
•
•
•
Water, 2-propanol, ethanol (analyte)
850 to 1049 nm at 1 nm intervals at 30, 40, 50, 60, and 70°C
Calibration: 13 mixtures from 0% to 67% at 30°C
Validation: 6 mixtures from 16% to 66% at 70°C
Primary: PC ethanol at 30°C
Non-analyte matrix (standardization set) N at 70°C
PC interferents water and 2-propoanol (2 samples)
Blanks (3 samples)
Constant analyte (CA, 5 samples)
38
PCTR2 Model Merit Landscapes
RMSEPC
RMSEN
100
100
0.8
0.37
η
1
0.37
0.8
0.6
0.6
1e-3
1e-3
0.4
2.73e-6
0
0.2
8.7e-7 9.5e-5 1.1e-2
1.2
0.4
2.73e-6
0.2
0
100
bˆ
8.7e-7 9.5e-5 1.1e-2
1.2
100
RMSEV
2
100
14
100
0.6
12
0.37
10
η
8
1e-3
0.5
0.37
0.4
1e-3
0.3
6
4
2.73e-6
0.2
2.73e-6
2
0
8.7e-7 9.5e-5 1.1e-2
λ
1.2
100
0.1
0
8.7e-7 9.5e-5 1.1e-2
λ
1.2
100
39
Model Updating From PC 30°C to 70°C
Method
(No. Samples)
N
PCTR2 (6)
Blanks
Int PC
PCTR2 (9)
bˆ
RMSEV
R2
13.41
0.037
0.97
Blanks
CA
17.23
0.050
0.99
PCTR2 (8)
CA
Int PC
15.42
0.093
0.98
PCTR2 (11)
Blanks
CA
Int PC
23.32
RR at 30°C
no update (13)
-
RR at 70°C (13)
-
2
4.52
4.93
0.069
0.99
0.258
0.66
0.115
0.86
Updated PC predicts as well as
secondary model predicting the
secondary validation samples
η
λ
5.2e-5
0.021
5.2e-5
0.054
5.2e-5
0.013
5.2e-5
0.054
0.043
0.034
-
40
PCTR2 and PLS Modeling Temperature
Updating analyte PC at 30°C to 70°C using interferent PC and blanks
Method
No.
Models
PCTR2
1454
PLS
84
RMSEPC
RMSEN
RMSEV
R2
13.70
0.0003
0.004
0.036
0.973
12.33
0.0003
0.013
0.054
0.963
η
Λ
0.00044
4.94
5 factors
0.17
PLS and PCTR2 predict similarly
PCTR2 RMSEV
PLS RMSEV
0.37
η
0.6
1
0.6
0.5
2
0.5
3
0.4
4
0.3
0.2
5
0.2
0.1
6
0.1
0.4
1e-3
0.3
bˆ
2.73e-6
2
0
8.7e-7 9.5e-5 1.1e-2
λ
1.2
100
Factors
100
0
2.8e-6
1.0e-3
λ
9.1e-1
100
41
PCTR2 Consensus Modeling Temperature
• On-going work
1. Cannot use R2, slope, and intercept for respective
predicted values of the PC and N
– Set thresholds for RMSEN, RMSPC, and bˆ based on
2
preliminary inspection of landscapes
bˆ
, RMSEN and RMSEPC
2
– Can further filter based on predicted values
• Majority vote
• Remove outliers
2. Combine predicted value of analyte pure component
sample with predicated non-analyte samples to obtain
R2, slope, andˆ intercept
b
2
42
PCTR Variants (Calibration or Maintenance)
1.No reference values

2
m in k b  1
T
a
2
η b
2
2
2
 λ Nb  0
2
2
2

2.With current condition sample reference values

2
m in k b  1  η b
T
a
2
2
2
2
 λ M b  yM
2
2
2

3.A combination of N and M
4.Replace b 2 with b 1 or Lb 2to obtain sparse models
43
Summary
• PCTR2 calibrates (updates) to current conditions without
reference samples
• Only a few new samples needed
• Can predict better than a full calibration
– More focused to orthogonalize to the sample matrix
T
T
T
T
yˆ  y k bˆ  y k bˆ  m b  ...  n bˆ
a

a

i

i

bias
variance
• Requires PC analyte spectrum
– Does not have to be matrix effected
• Requires non-analyte samples
– Can be estimated with reference samples
44
Other TR Variants
Expression
m in

m in
 Xb  y
m in
 Xb  y
m in



2
2
2
2
Xb  y
2
2
T
Xb  y
 Xb  y
m in  X b  y
m in
2
2
m in k a b  1
m in
2
Xb  y
2
2
2
2
2
2
2
η
η
2
2
η
η
2
2
2
2
2
Lb
 η Lb
η
2
Lb
2

RR when L = I; includes variable selection
when L = diag and approximates an 1-norm
2
2
1
2
2
2
Lb
2
 λ M b  yM
 λ Nb
2
L b  b * 
η b
1
2
2


2
2

2
2
2
1

Model updating; includes variable selection

Model updating with variable selection;
approximates 0-norm when L = diag

Model updating with robustness to the
standardization set M
Calibration or updating without reference
samples
Calibration to target model b*
Adaptive LASSO and LASSO when L = I
λ b
2
1
2
 λ M b  yM
 λ M b  yM
Lb
 η Lb
2
2

Claerbout JF, Muir F. Geophysics 1973; 38: 826-844
Elastic net
Other On-Going Consensus Modeling
• In addition to combining a set of models, can combine
TR2, PLS, PCTR2, … sets of model predictions
Consensus TR2 models
Consensus PLS models
Consensus PCTR2 models
Final prediction
Consensus TR2-1 models
46
```