ARL-TR-7925 • JAN 2017

US Army Research Laboratory

Quantal Response: Nonparametric

Modeling

by Joseph C Collins

Approved for public release; distribution is unlimited.
NOTICES

Disclaimers

The findings in this report are not to be construed as an official Department of the

Army position unless so designated by other authorized documents.

Citation of manufacturer’s or trade names does not constitute an official endorsement or approval of the use thereof.

Destroy this report when it is no longer needed. Do not return it to the originator.
ARL-TR-7925 • JAN 2017

US Army Research Laboratory

Quantal Response: Nonparametric

Modeling

by Joseph C Collins

Survivability/Lethality Analysis Directorate, ARL

Approved for public release; distribution is unlimited.
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704‐0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the

data needed, and completing and reviewing the collection information.  Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing

the burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704‐0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202‐

4302.  Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a

currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD‐MM‐YYYY) 2. REPORT TYPE 3. DATES COVERED (From ‐ To)

4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S) 5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION REPORT

NUMBER

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S)

11. SPONSOR/MONITOR'S REPORT NUMBER(S)

12. DISTRIBUTION/AVAILABILITY STATEMENT

13. SUPPLEMENTARY NOTES

14. ABSTRACT

15. SUBJECT TERMS

16. SECURITY CLASSIFICATION OF:   17. LIMITATION

       OF  

       ABSTRACT

18. NUMBER

       OF  

        PAGES

19a. NAME OF RESPONSIBLE PERSON

a. REPORT b. ABSTRACT c. THIS PAGE 19b. TELEPHONE NUMBER (Include area code)

Standard Form 298 (Rev. 8/98)

Prescribed by ANSI Std. Z39.18

January 2017 Technical Report

Quantal Response: Nonparametric Modeling

Joseph C Collins

ARL-TR-7925

Approved for public release; distribution is unlimited.

October 2013–December 2016

US Army Research Laboratory

ATTN: RDRL-SLB-D

Aberdeen Proving Ground, MD 21005-5068

Author’s Email: <joseph.c.collins38.civ@mail.mil>

Binary response to a continuous stimulus was originally modeled as a cumulative distribution of limit quantities interpreted as

an increasing functional relationship between stimulus and probability of response. The Generalized Linear Model approach

does not make use of the limit distribution but allows arbitrary functional relationships including both classical higher-order

polynomials and also nonparametric spline models based on penalized maximum likelihood estimation.

quantal response, polynomial, penalized likelihood, smoothing, spline

70

Joseph C Collins

Unclassified Unclassified Unclassified UU 410-278-6832

ii
Approved for public release; distribution is unlimited.

Contents

List of Figures iv

List of Tables v

1. Introduction 1

2. General QR Models 1

3. The Chang-Bodt QR Model 2

4. Parametric QR Models 5

5. Nonparametric QR Models 8

5.1 The P-spline QR Model 8

5.2 The S-spline QR Model 12

6. Model Comparison 16

7. Conclusions and Recommendations 18

8. References 19

Appendix A. The Linear Model 21

Appendix B. The Generalized Linear Model 33

Appendix C. B-splines 47

Appendix D. Matrix Differentiation Operators 53

List of Symbols, Abbreviations, and Acronyms 58

Distribution List 59

iii
Approved for public release; distribution is unlimited.

List of Figures

Fig. 1 Chang-Bodt model ..................................................................4

Fig. 2 Regression basis sets................................................................6

Fig. 3 Logistic regression ..................................................................7

Fig. 4 P-spline optimization ...............................................................9

Fig. 5 Optimal and adjusted P-spline solutions ...................................... 10

Fig. 6 Adjusted P-spline solution with 90% confidence intervals................ 11

Fig. 7 S-spline optimization ............................................................. 13

Fig. 8 Optimal S-spline solutions....................................................... 14

Fig. 9 Adjusted S-spline solution with 90% confidence intervals................ 15

Fig. 10 All solutions ........................................................................ 17

Fig. A-1 P-spline optimization, linear model............................................ 26

Fig. A-2 S-spline optimization, linear model............................................ 28

Fig. A-3 Smoothing parameter selection, linear model................................ 30

Fig. A-4 GCV-optimal fits, linear model ................................................. 31

Fig. B-1 P-spline optimization, GLM..................................................... 40

Fig. B-2 S-spline optimization, GLM..................................................... 41

Fig. B-3 Smoothing parameter selection, GLM ........................................ 42

Fig. B-4 GCV-optimal fits, GLM .......................................................... 43

Fig. C-1 B-spline basis, uniform........................................................... 49

Fig. C-2 B-spline basis, compact .......................................................... 50

Fig. C-3 B-spline basis derivatives, uniform ............................................ 51

Fig. C-4 B-spline basis derivatives, compact ............................................ 52

iv
Approved for public release; distribution is unlimited.

List of Tables

Table 1 Chang-Bodt model parameters....................................................3

Table 2 Goodness of fit comparison...................................................... 16

v
Approved for public release; distribution is unlimited.

Intentionally left blank.

vi
Approved for public release; distribution is unlimited.

1. Introduction

The quantal response (QR) model seeks to characterize the value of a binary response y ∈ {0,1} as a function of continuous stimulus x. The original formulation

postulates an unobservable stimulus limit L, which determines the response via

y =









0, x < L

1, x > L .

(1)

If L is known and constant, then y is a step function of x with jump at x = L.

Otherwise, suppose L is random with cumulative distribution function (CDF)

FL (t) = Pr [ L 6 t ] . (2)

The distribution of L determines the QR model probability of response as

P(x) = Pr [ y = 1 | x ] = Pr [ L 6 x ] = FL (x) . (3)

Assuming a specific distribution family for FL such as logistic, G(x) = (1+e

−x

)

−1

,

or normal, G(x) = (2π)

−1/2

 x

−∞

e

−t

2

/2 dt, fitting the model amounts to estimating

the limit distribution location and scale parameters m and s in

P(x) = FL (x) = G



x − m

s



. (4)

The CDF of L is interpreted as a functional model for the probability of response.

The existence of L is not necessary, and as such the QR model is a special case of

the Generalized Linear Model (GLM). See Collins1,2for details. The requirement

of increasing P can be relaxed, and then P is no longer equivalent to any CDF

(which, of course, must always be increasing). What remains is a purely functional

model with no implicit limit distribution. Several such formulations follow.

2. General QR Models

A binary response y ∈ {0,1} has an expected value depending on stimulus x,

E [ y | x] = P(x) . (5)

1
Approved for public release; distribution is unlimited.

The conditional distribution of y is Bernoulli with the indicated mean, so

Pr [ y = 1 | x] = E [ y | x] = 1 − Pr [y = 0 | x] . (6)

In the usual case, response is an increasing function of x, so the response function

P must be a CDF. A standard choice for P uses the logistic CDF G(x) = (1 +

e

−x

)

−1 with the linear parameterization P(x) = G(a + bx) or the location-scale

parameterization P(x) = G ((x − m)/s).

This is inadequate, for example, in the presence of the shattergap phenomenon,

when the probability of penetration (y) is observed to decrease in some velocity

(x) range. An example data set is taken from Chang and Bodt,3 described therein

as “Results of 69 Ballistic Shots on Phase II Al2O3/Kevlar Armor Plates”, although

there are only 68 data points in the report.

Chang and Bodt also develop a specific parametric model that must be analyzed

from first principles (not GLM). This is presented in Section 3.

The remaining approaches are all based on GLM. Standard Bernoulli GLM can

estimate an arbitrarily complex response, not necessarily increasing, as in Section 4.

Nonparametric estimation of the response is accomplished with penalized B-spline

models, as in Section 5.1, or smoothing spline models, as in Section 5.2.

3. The Chang-Bodt QR Model

The Chang-Bodt3 model is

P(x) =



1 − Pz(x,mz,sz)



· P1(x,m1,s1) + Pz(x,mz,sz) · P2(x,m2,s2) , (7)

where the Pi are location-scale CDFs. P1 is the (monotonic) probability of penetration for an unshattered threat, and P2 is the probability of penetration for a shattered

threat. Pzis the probability of shatter. For low x, shatter is unlikely, Pz ∼ 0, and

P ∼ P1. For high x, shatter is likely, Pz ∼ 1, and P ∼ P2. In the intermediate x

range, the mixture is weighted according to the increasing shatter probability Pz.

For convenience, let R = 1 − P.

2
Approved for public release; distribution is unlimited.

The likelihood for a single observation (xi, yi) is

Li = Rz(xi) · P1(xi)

yi R1(xi)1−yi + Pz(xi) · P2(xi)yi R2(xi)1−yi

, (8)

where the parameters of Pk are (mk ,sk ).

Li =









Rz(xi) · R1(xi) + Pz(xi) · R2(xi) , yi = 0

Rz(xi) · P1(xi) + Pz(xi) · P2(xi) , yi = 1

. (9)

For the logistic CDF, we have G(x) =



1 + exp(−x)

−1

. The usual location-scale

parameterization is P(x) = G



(x − m)/s



.

Parameter estimates can be obtained by numerical optimization of the negative loglikelihood Λ = −

P

log Li. The parameters reported by Chang and Bodt are not

optimal estimates. Both sets of parameter estimates are given in Table 1.

Table 1 Chang-Bodt model parameters

Model m1 s1 m2 s2 mz sz Λ

Reported 1650 100 2550 100 2050 100 22.58

Estimated 1646.97 90.49 2529.34 75.75 2020.15 87.80 21.91

See Fig. 1 for the response curves.

3
Approved for public release; distribution is unlimited.

1400 1600 1800 2000 2200 2400 2600 2800

0.0

0.2

0.4

0.6

0.8

1.0

x

P

report

estimate

Fig. 1 Chang-Bodt model

4
Approved for public release; distribution is unlimited.

4. Parametric QR Models

GLM with the Bernoulli response provides a way to characterize a possibly nonincreasing function as P(x) = G



f (x)



, where G is a monotone function such as the

standard logistic, normal, or Cauchy CDF. Finite-dimensional parametric models

are obtained by using some basis ( f1, f2,. . .), and then f (x) = X

t β for fixed k

where β = ( β1,. . . , βk ) and X =



f1(x),. . . , f k (x)



. So, we get

P(x) = G(X

t

β) . (10)

This accounts for the (monotonic) basic linear or location-scale model

P(x) = G(b0 + b1 x) = G



x − m

s



(11)

and polynomials of arbitrary degree

P(x) = G *

,

X

k

i=0

bi x

i+

-

. (12)

The canonical polynomial basis is given by fi(x) = x

i

. In practice, we use an orthogonal polynomial basis, with degree( fi) = i and 

fif j = 0 if i , j. This

eliminates numerical problems and provides the same solution as the canonical basis for each k. Other popular choices for basis sets include the natural spline and

B-spline. See Fig. 2 for basis set examples with dimension k = 5 where various

colors distinguish the basis elements. Figure 3 shows logistic response estimates

for these 3 basis sets with dimension d = 3,4,5,6.

5
Approved for public release; distribution is unlimited.

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.2

0.4

0.6

0.8

1.0 canonical polynomial

x

y

0.0 0.2 0.4 0.6 0.8 1.0

−0.3

−0.2

−0.1

0.0

0.1

0.2

0.3 orthogonal polynomial

x

y

0.0 0.2 0.4 0.6 0.8 1.0

−0.4

−0.2

0.0

0.2

0.4

0.6

natural spline

x

y

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.2

0.4

0.6

0.8

1.0 B−spline

x

y

Fig. 2 Regression basis sets

6
Approved for public release; distribution is unlimited.

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

d = 3

polynomial

B−spline

N−spline

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

d = 4

polynomial

B−spline

N−spline

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

d = 5

polynomial

B−spline

N−spline

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

d = 6

polynomial

B−spline

N−spline Fig. 3 Logistic regression

7
Approved for public release; distribution is unlimited.

5. Nonparametric QR Models

Nonparametric linear models are presented in Appendix A. These arise from a process that can be described as penalized maximum likelihood estimation, penalized

least squares, or smoothing.

In Appendix B, the iteratively reweighted least squares (IRLS) GLM estimation

procedure for QR models is adapted to use nonparametric penalized linear models.

This gives rise to nonparametric penalized QR models.

5.1 The P-spline QR Model

Eilers and Marx4 propose using an overfitted B-spline model and then penalizing

to produce a smooth fit. This is called a P-spline model. In this application, we fit

a penalized GLM with a degree 3 B-spline basis of size p = 32 penalized with the

second derivative D. The smoothing operator is S = e

λD.

Selection of the smoothing parameter is usually accomplished by optimizing some

information or cross-validation quantity such as the Akaike information criterion

(AIC), ordinary cross-validation (OVC), or generalized cross-validation (GCV) as

described in Appendix A. There are other such quantities and other methods, and no

single procedure is known to give the best solution. In fact, no single procedure even

works for all data sets. So the choice of smoothing parameter selection procedure

is itself somewhat subjective. To compensate for the perception that the procedures

tend to oversmooth the response, an adjusted smoothing parameter is computed as

the minimum of the 3 values obtained less their range. For this data, λaic = −6.2,

λgcv = −6.0, λocv = −5.7, and λadj = −6.7.

See Fig. 4 for fits with various values of the smoothing parameter. Optimal solutions

for the 3 smoothing selection methods along with the adjusted solution are shown

in Fig. 5. The adjusted solution and confidence intervals are shown in Fig. 6.

8
Approved for public release; distribution is unlimited.

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −9.2

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −8.6

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −7.9

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −7.3

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −6.7

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −6.1

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −5.4

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −4.8

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −4.2

Fig. 4 P-spline optimization

9
Approved for public release; distribution is unlimited.

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

AIC opt

GCV opt

OCV opt

adjusted

Fig. 5 Optimal and adjusted P-spline solutions

10
Approved for public release; distribution is unlimited.

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

Fig. 6 Adjusted P-spline solution with 90% confidence intervals

11
Approved for public release; distribution is unlimited.

5.2 The S-spline QR Model

The smoothing spline, or S-spline, QR model has the form P(x) = G



f (x)



, where

f is an S-spline linear model with knots at the stimulus data points, and G is an

arbitrary link function. Wahba5is the standard reference for smoothing splines.

Smoothing splines are obtained through optimization in certain function spaces,

and the GLM implementation described in Section B.3 accomplishes this using

the standard IRLS GLM algorithm. In this application, we fit a cubic smoothing

spline GLM by penalizing with the second derivative D. The smoothing operator is

S = e

λD.

We allow for multiple observations at a single stimulus by averaging the response

and multiplying the weight by the observation multiplicity at that level. However,

the sum of squared errors (SSE) and Λ are computed from the original data, so they

are comparable with the other models.

Smoothing parameter selection methodology is the same as for the P-spline, Section 5.1. Typically, different λ values are obtained.

See Fig. 7 for fits with various values of the smoothing parameter. Optimal solutions

for all 3 smoothing selection methods and the adjusted solution are shown in Fig. 8.

The adjusted solution and confidence intervals are shown in Fig. 9.

12
Approved for public release; distribution is unlimited.

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −9.4

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −8.8

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −8.2

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −7.5

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −6.9

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −6.3

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −5.7

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −5.0

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −4.4

Fig. 7 S-spline optimization

13
Approved for public release; distribution is unlimited.

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

AIC opt

GCV opt

OCV opt

adjusted

Fig. 8 Optimal S-spline solutions

14
Approved for public release; distribution is unlimited.

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

Fig. 9 Adjusted S-spline solution with 90% confidence intervals

15
Approved for public release; distribution is unlimited.

6. Model Comparison

Table 2 displays fit statistics and Fig. 10 shows the superimposed solutions for selected methods: degree 6 parametric models, reported and estimated Chang-Bodt

(CB), and adjusted P-spline and S-spline. The adjusted smoothing spline solution

is the best fit in terms of squared error and maximum likelihood measures.

Table 2 Goodness of fit comparison

Model SSE Λ

Polynomial (6) 7.237 22.49

B-spline (6) 7.212 22.34

Natural spline (6) 7.201 22.07

CB reported 7.310 22.58

CB estimated 7.186 21.91

P-spline adjusted 7.045 21.98

S-spline adjusted 6.847 21.37

16
Approved for public release; distribution is unlimited.

1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

poly(6)

B−spline(6)

N−spline(6)

CB rpt

CB est

P−spline adj

S−spline adj

Fig. 10 All solutions

17
Approved for public release; distribution is unlimited.

7. Conclusions and Recommendations

The simple location-scale quantal response model may not adequately capture the

behavior of observed phenomena.

Higher-order polynomial and finite-dimensional spline basis models allow for more

complicated responses as the polynomial degree or spline basis dimension increases.

Penalized B-spline (P-spline) and smoothing spline (S-spline) models offer the most

flexibility as these are nonparametric (not constrained to any particular functional

form). These should be useful in identifying nonstandard behavior via statistical

goodness-of-fit tests.

18
Approved for public release; distribution is unlimited.

8. References

1. Collins J. Quantal response: practical sensitivity testing. Aberdeen Proving

Ground (MD): Army Research Laboratory (US); 2012 Jun. Report No.: ARLTR-6022. 1

2. Collins J. Quantal response: estimation and inference. Aberdeen Proving

Ground (MD): Army Research Laboratory (US); 2014 Sep. Report No.: ARLTR-7088. 1

3. Chang AL, Bodt BA. JTCG/AS interlaboratory ballistic test program—final report. Aberdeen Proving Ground (MD): Army Research Laboratory (US); 1997

Dec. Report No.: ARL-TR-1577. 2

4. Eilers PHC, Marx BD. Flexible smoothing with B-splines and penalties. Statistical Science. 1996;11(2):89–121. 8

5. Wahba G. Spline models for observational data. Philadelphia (PA): Society for

Industrial and Applied Mathematics; 1990. 12

19
Approved for public release; distribution is unlimited.

Intentionally left blank.

20
Approved for public release; distribution is unlimited.

Appendix A. The Linear Model

21
Approved for public release; distribution is unlimited.

A.1 The Basic Linear Model

The usual linear model is

Y = X β + ε , (A-1)

where Y is an n×1 response, X is an n×p independent matrix, β is a p×1 parameter,

and the n × 1 error ε ∼ N(0,σ2In) is normally distributed.

For any positive-definite symmetric matrix W, the corresponding weighted inner

product and norm are, respectively, hx, yiW = x

tW y and kxkW = hx, xi

1/2

W

. So,

in general, the density of a normal vector with mean M and variance-covariance

matrix V can be written as

f (x) = (2π)

−n/2

|V |

−1/2

expf

−

1

2

kx − Mk

2

V−1

g

. (A-2)

For the model of Eq. A-1, we have EY = M = X β and VarY = V = σ

2

In. Each

column of X is a linear predictor, and the model is

yi =

X

p

j=1

Xij βj. (A-3)

We can work with a p-parameter model for a single predictor v by choosing a set of

fixed basis functions { f1,. . . , f p} and setting Xij = f j(vi).

yi =

X

p

j=1

βjf j(vi) . (A-4)

For example, the choice of f j(v) = v

j−1 gives the polynomial model

y = β0 + β1v + β2v

2 + · · · βp−1vp−1

. (A-5)

See Eq. A-2. Solution by least squares is equivalent to maximum likelihood estimation for normal error, and the criterion is to choose u that minimizes Q = ε

t

ε =

kεk

2 = kY − X βk2

since ε = Y − X β. This is

Q = kY k

2 − 2βtXtY + kX βk2

, (A-6)

22
Approved for public release; distribution is unlimited.

and by setting the derivative to 0,

dQ

dβ

= −2X

tY + 2XtX β = 0 , (A-7)

we obtain the normal equations

X

tX β = XtY (A-8)

with solution

βˆ = (X

tX)−1XtY (A-9)

and response estimate

Yˆ = HY , (A-10)

where the so-called hat matrix is

H = X(X

tX)−1Xt

. (A-11)

Note that E βˆ = β and Var βˆ = σ

2

(X

tX)−1

.

A.2 The Weighted Model

When the error is N(0,Σ), the correct inner product is weighted by the symmetric

W = Σ

−1

, and so Q = ε

tWε = kεk2

W

= kY − X βk

2

W

. This is

Q = kY k

2

W − 2β

tXtWY + kX βk2

W . (A-12)

Then

dQ

dβ

= −2X

tWY + 2XtW X β = 0 , (A-13)

the normal equations are

X

tW X β = XtWY , (A-14)

the solution is

β = (X

tW X)−1XtWY , (A-15)

and the response estimate is Yˆ = HY, where

H = X(X

tW X)−1XtW . (A-16)

Note that E βˆ = β and Var βˆ = (X

tW X)−1

.

23
Approved for public release; distribution is unlimited.

For these models, the likelihood function is

L( β) = (2π)

−n/2

|Σ|

−1/2

expf−

1

2

kY − X βk

2

Σ−1

g

, (A-17)

its log derivative is

d

dβ

log L( β) = X

0Σ−1

(Y − X β) , (A-18)

and the information matrix is

Mβ = X

0Σ−1X . (A-19)

A.3 Smoothing

If the system is ill-conditioned, we can maximize a penalized likelihood function

LS ( β) = L( β) · exp−

1

2

kS βk

2



. (A-20)

Equivalently, penalization can be applied to the least-squares formulation by minimizing Q = kεk

2

W

+ kS βk

2

for some smoothing operator S. This is

Q = kY k

2

W − 2β

tXtWY + kX βk2

W + k βk

2

S

t S

. (A-21)

Then

dQ

dβ

= −2X

tWY + 2XtW X β + 2St

S β = 0 , (A-22)

and the normal equations are

(X

tW X + St

S) β = X

tWY (A-23)

with solution

β = (X

tW X + St

S)

−1XtWY (A-24)

and response estimate Yˆ = HY where

H = X(X

tW X + St

S)

−1XtW . (A-25)

This is equivalent to the model Y

∗ = X∗ β + ε∗ with weight W∗ where Y∗ =

f

Y

0

g

,

24
Approved for public release; distribution is unlimited.

X

∗ =

f

X

S

g

, and W∗ =

f

W 0

0 I

g

. Then the normal equations are

f

X

t St

g













W 0

0 I

























X

S













β =

f

X

t St

g













W 0

0 I

























Y

0













f

X

t St

g













W X

S













β =

f

X

t St

g













WY

0













(X

tW X + St

S) β = X

tWY (A-26)

as required. This so-called augmented representation is useful for computation.

Modern software for linear least-squares estimation operates on Eqs. A-8, A-14,

and A-26 through the response vector Y, the weight vector W, the smoothing vector

S, and the design matrix X. The normal equations are solved efficiently without

inverting the design matrix, and we get parameter estimates and diagnostics such as

the parameter variance and hat matrix diagonal.

A.4 Nonparametric P-spline

Smoothing or penalization can be applied to an overfitted parametric model as given

by Eq. A-4. Suppose f has the particular form

f (x) =

X

p

j=1

βjf j(xi) (A-27)

for a B-spline basis ( f1,. . . , f p) and D is a linear differential operator. Computational details for the B-spline are in Appendix C, and for the differential operator

see Appendix D. The solution is the minimizer ( β1,. . . , βp) of

Xn

i=1



yi −

X

p

j=1

βjf j(xi)

2

+ λ

2

 b

a



D

X

p

j=1

βjf j



(u)

2

du (A-28)

and is again given exactly by Eq. A-21. This works because of the ordered partial partition property of B-spline bases: smoothing the coefficient vector in fact

smooths the solution. Eilers and Marx1 call this a P-spline (penalized B-spline)

model. Figure A-1 depicts the effect of smoothing parameter variation for the Pspline model.

1Eilers PHC, Marx BD. Flexible smoothing with B-splines and penalties. Statistical Science.

1996;11(2):89–121.

25
Approved for public release; distribution is unlimited.

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −10.5

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −9.5

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −8.5

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −7.5

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −6.5

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −5.5

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −4.5

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −3.5

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −2.5

Fig. A-1 P-spline optimization, linear model

26
Approved for public release; distribution is unlimited.

A.5 Nonparametric S-spline

The smoothing spline is another manifestation of penalized maximum likelihood or

least squares estimation. A smoothing spline, or S-spline, f is the minimizer of

Xn

i=1



yi − f (xi)

2 + λ

2

 b

a



(D f )(u)

2

du , (A-29)

where f is in some suitable function space, and D is a linear differential operator.

The solution is a spline with knots at the data points. If D has order p, the solution

is piecewise polynomial of degree 2p − 1 with continuous derivatives of orders

0,1,2,. . . ,2p − 2. See Wahba2.

The discrete representation of this problem is Eq. A-21, with X = I and S taken to

be a scalar multiple of a differential operator D, so S = λD and S

tS = λ2DtD. The

representation D of D is given in Appendix D. The discrete formulation provides

the exact solution of Eq. A-29.

When D is the second derivative, the solution is a cubic spline. The solution is given

at the data points, and spline interpolation can be used to evaluate the response at

other values.

Figure A-2 depicts the effect of smoothing parameter variation for the S-spline

model.

2Wahba G. Spline models for observational data. Philadelphia (PA): Society for Industrial and

Applied Mathematics; 1990.

27
Approved for public release; distribution is unlimited.

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −9.7

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −9.0

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −8.2

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −7.5

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −6.7

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −6.0

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −5.2

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −4.5

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

λ = −3.7

Fig. A-2 S-spline optimization, linear model

28
Approved for public release; distribution is unlimited.

A.6 Selecting the Smoothing Parameter

Choice of smoothing parameter λ is usually based on minimization of one of the

following quantities. The ordinary cross-validation (OCV) is

OCV =

1

n

Xn

i=1

 

yi − yˆi

1 − Hii!2

. (A-30)

The generalized cross-validation (GCV) is

GCV =

1

n

Pn

i=1



yi − yˆi

2

f

1

n

Tr(I − H)

g2

. (A-31)

The Akaike information criterion (AIC) is

AIC = ∆ + 2 · Tr H , (A-32)

where the deviance ∆ = −2 log(Lreduced/Lfull) is defined in terms of the likelihood

function L.

For normal error, Lfull = 1, and based on Eq. A-2, we have

log Lreduced = −

n

2

log(2π) −

1

2

log |Σ| −

1

2

kεk

2

Σ−1

, (A-33)

where the residuals are ε = Y − X β and the error variance is Σ.

For independent and identically distributed (IID) errors, Σ = σ

2

In and |Σ| = (σ

2

)

n

.

And so

log Lreduced = −

n

2

log(2π) −

n

2

log(σ

2

) −

1

2σ2

ε

t

ε . (A-34)

We replace σ

2 by the estimator εt

ε/n and get

log Lreduced = −

n

2

log(2π) −

n

2

log ε

t

ε

n



−

n

2

= −

n

2

f

log2π

ε

t

ε

n



+ 1

g

(A-35)

so

∆ = n

f

log2π

ε

t

ε

n



+ 1

g

. (A-36)

See Fig. A-3 for graphs of OCV, GCV, and AIC as functions of λ. Note that AIC

optimization fails for the S-spline, as no minimum value is obtained. See Fig. A-4

for GVC-optimal P-spline and S-spline solutions. With a maximum y-difference of

less than 0.01, the solutions are practically indistinguishable.

29
Approved for public release; distribution is unlimited.

−10 −8 −6 −4 −2 0

−20

−10

0

10

20

30

40

λ

aic, P−spline

−10 −8 −6 −4 −2 0

0.0008

0.0010

0.0012

0.0014

0.0016

0.0018

λ

gcv, P−spline

−10 −8 −6 −4 −2 0

3

4

5

6

7

λ

ocv, P−spline

−10 −8 −6 −4 −2 0

−150

−100

−50

0

50

λ

aic, S−spline

−10 −8 −6 −4 −2 0

0.0008

0.0010

0.0012

0.0014

0.0016

0.0018

λ

gcv, S−spline

−10 −8 −6 −4 −2 0

3

4

5

6

7

λ

ocv, S−spline

Fig. A-3 Smoothing parameter selection, linear model

30
Approved for public release; distribution is unlimited.

0 5 10 15

−1.0

−0.5

0.0

0.5

x

y

S−spline

P−spline

Fig. A-4 GCV-optimal fits, linear model

31
Approved for public release; distribution is unlimited.

Intentionally left blank.

32
Approved for public release; distribution is unlimited.

Appendix B. The Generalized Linear Model

33
Approved for public release; distribution is unlimited.

B.1 Generalized Linear Model Formulation

In the Generalized Linear Model (GLM), the response Y has an arbitrary distribution, η = X β is a linear function of the parameter β, and the mean response is

modeled as

μ = E[Y | X ] = G(η) (B-1)

for some monotone link function G with derivative g = G

0

. (Some authors call G

−1

the link.) Response distributions are taken to be from a single-parameter exponential

family, with the form

f (y, θ,ψ) = exp "

yθ − b(θ)

a(ψ)

+ c(y,ψ)

#

. (B-2)

The parameter θ is to be estimated, and ψ is a nuisance parameter.

With ` = log f , we calculate the moments of Y in terms of exponential family

components.

E

"

d

dθ

`(y, θ,ψ)

#

= 0 (B-3)

because f (y, θ,ψ) dy = 1, and under suitable regularity conditions d

dθ



f (y, θ,ψ) dy

=



d

dθ

f (y, θ,ψ) dy =



d

dθ

`(y, θ,ψ)· f (y, θ,ψ) dy = 0. Therefore, E[



y − b

0

(θ)



/a(ψ)]

= 0 and

E[Y] = μ = G(η) = b

0

(θ) . (B-4)

Also, as usual,

E













 

d

dθ

`(y, θ,ψ)

!2 











= − E

"

d

2

dθ

2

`(y, θ,ψ)

#

(B-5)

because Var f

d

dθ

`

g

= E

f

(

d

dθ

`)

2

g

=



(

d

dθ

`)

2

· f dy =



d

dθ

` ·

d

dθ

f dy

=

d

dθ

` · f





∞

−∞ −



d

2

dθ

2 ` · f dy = − E

f

d

2

dθ

2 `

g

.

Therefore, E[(y − μ)

2

/a(ψ)

2

] = − E[−b

00(θ)/a(ψ)] and

Var[Y] = v(μ)a(ψ) = b

00(θ)a(ψ) (B-6)

where v(μ) = Var[Y]/a(ψ) = b

00(θ) = g(η)

dη

dθ

.

In the case that η = θ, and hence μ = G(η) = G(θ), we say that G is the canonical

link function. Then μ = b

0

(θ) = G(θ) = G(η) and v(μ) = b

00(θ) = g(θ) = g(η).

34
Approved for public release; distribution is unlimited.

B.2 Estimation

Let [x1,. . . , xn] = X

t

, so the (column) vector xiis row i of X, ηi = x

t

i

β, and

θi = θ(ηi). Maximum likelihood estimation for the GLM is accomplished by maximizing the log-likelihood function

L =

Xn

i=1

`(yi, θi,ψ) =

Xn

i=1

"

yiθi − b(θi)

a(ψ)

+ c(yi,ψ)

#

. (B-7)

This is a weighted least squares problem where the design and weight depend on

the unknown parameter, and it can be solved iteratively by the Newton-Raphson

method.

B.2.1 The Newton-Raphson Method

In 1 dimension, we find a zero of F by linearizing and updating the current argument

xo to x by solving F(xo) + (x − xo)F

0

(xo) = 0 to get x = xo − F(xo)/F

0

(xo). We

optimize F by setting F

0 = 0, so the update is x = xo − F

0

(xo)/F

00(xo).

The vector version is x = xo − [

d

2

dxxt F(xo)]

−1 d

dx F(xo). If we take x = xo + δ,

the increment δ = x − xo satisfies [ d

2

dxxt F(xo)]δ = −

d

dx F(xo). So we need some

derivatives.

B.2.2 Gradient

Differentiating, we have the gradient (vector of first derivatives)

D( β) =

dL

dβ

= a(ψ)

−1

·

Xn

i=1



yi − b

0

(θi)



·

dθi

dβ

. (B-8)

Since d

dβ

b

0

(θi) = b

00(θi)

dθi

dβ

= v(μi)

dθi

dβ

and d

dβ

b

0

(θi) =

dμi

dβ

=

d

dβG(ηi) =

d

dβG(x

t

i

β) =

g(ηi)xi, we get

dθi

dβ

=

g(ηi)

v(μi)

xi. (B-9)

So the gradient is

D( β) = a(ψ)

−1

·

Xn

i=1



yi − μi



·

g(ηi)

v(μi)

xi = a(ψ)

−1

· X

tWFUF , (B-10)

35
Approved for public release; distribution is unlimited.

where, for i = 1,. . . ,n, the diagonal weight matrix WF has elements wFi = g(ηi)

2

/v(μi)

WF = diag

g(η1)

2

v(μ1)

, · · · ,

g(ηn)

2

v(μn)



(B-11)

and the centered/scaled response vector UF has elements uFi = (yi − μi)/g(ηi)

UF =



y1 − μ1

g(η1)

, · · · ,

yn − μn

g(ηn)



. (B-12)

B.2.3 Hessian

Using d

dβ

v(μi) = v

0

(μi)

d

dβG(x

t

i

β) = v

0

(μi)g(x

t

i

β)xi,

d

2

dβ βt

θi =

g

0

(ηi)v(μi) − g(ηi)

2

v

0

(μi)

v(μi)

2

xi x

t

i

, (B-13)

and the Hessian (matrix of second derivatives) is

H ( β) =

d

2

dβ βt

L = a(ψ)

−1

·

Xn

i=1

"

−b

00(θi)

dθi

dβ

dθi

dβ

t

+



yi − μi

 d

2

dβ βt

θi

#

= a(ψ)

−1

·

Xn

i=1

"

−

g(ηi)

2

v(μi)

+



yi − μi

 g

0

(ηi)v(μi) − g(ηi)

2

v

0

(μi)

v(μi)

2

#

xi x

t

i

= −a(ψ)

−1

· X

tWN X (B-14)

where

WN = WF − WD (B-15)

and WD is a diagonal matrix with diagonal elements

wDi =



yi − μi

 g

0

(ηi)v(μi) − g(ηi)

2

v

0

(μi)

v(μi)

2

. (B-16)

Since E[WD] = 0, the expected value of the Hessian is

E H ( β) = −a(ψ)

−1

· X

tWF X . (B-17)

The Fisher Information Matrix is

Mβ = E



d

dβ

L · d

dβ

L

t



= − E



d

2

dβ βt

L



= − E H ( β) = a(ψ)

−1

· X

tWF X , (B-18)

36
Approved for public release; distribution is unlimited.

and the asymptotic estimator distribution is N



β, a(ψ) · (X

0WF X)−1



.

B.2.4 Newton-Raphson

Now, apply the Newton-Raphson algorithm to iteratively solve the optimization.

For GLM, the Newton-Raphson update is β = βo + δ where

H ( βo)δ = −D( βo)

(X

tWN X)δ = XtWFUF

(X

tWN X)δ = XtWNUN (B-19)

with

UN = W−1

N WFUF . (B-20)

These are the normal equations for minimization of Q = kUN − Xδk

2

WN

. Both WN

and UN depend on βo. The normal equations can be solved iteratively with an initial

guess βo by calculating η = X βo, μ = G(η), g, g

0

, v, v

0

, WF, UF, WD, WN, and

UN. Then solve for δ. The updated solution is β = βo + δ. Now replace βo with β,

and repeat. This is iteratively reweighted least squares with the Newton-Raphson

update.

B.2.5 Fisher Scoring

For the GLM, the Fisher scoring update uses E H in place of H to get

E H ( βo)δ = −D( βo)

(X

tWF X)δ = XtWFUF , (B-21)

which are the normal equations for minimization of Q = kUF − Xδk

2

WF

. Both WF

and UF depend on βo. The normal equations can be solved iteratively with an initial

guess βo by calculating η = X βo, μ = G(η), g, v, WF, and UF. Then solve for δ,

update, and repeat as above. This is iteratively reweighted least squares with Fisher

scoring.

B.3 Smoothing

Smoothing is accomplished by maximizing a penalized objective

LS = L − 1

2

kS βk

2

(B-22)

37
Approved for public release; distribution is unlimited.

with gradient

DS ( β) =

dLS

dβ

= X

tWFUF − St

S β , (B-23)

Hessian matrix

HS ( β) =

d

2

dβ βt

LS = −X

tWN X − St

S , (B-24)

and expected Hessian

E HS ( β) = −X

tWF X − St

S . (B-25)

The Newton-Raphson update equations are

HS ( βo)δ = −DS ( βo)

(X

tWN X + St

S)( β − βo) = X

tWFUF − St

S βo

(X

tWN X + St

S) β = (X

tWN X + St

S) βo + X

tWFUF − St

S βo

· = X

tWN X βo + XtWFUF

· = X

tWN X βo + XtWNW−1

N WFUF

· = X

tWN (X βo + W−1

N WFUF)

· = X

tWN (X βo + UN )

(X

tWN X + St

S) β = X

tWN ZN , (B-26)

where

ZN = X βo + UN . (B-27)

These are the normal equations for minimization of Q = kZN − X βk

2

WN

+ kS βk

2

.

The update increment δ is now implicit. Replace βo with β and repeat.

The Fisher update equations are E HS ( βo) β = E HS ( βo) βo − DS ( βo), or

E HS ( βo)δ = −DS ( βo)

(X

tWF X + St

S)( β − βo) = X

tWFUF − St

S βo

(X

tWF X + St

S) β = (X

tWF X + St

S) βo + X

tWFUF − St

S βo

· = X

tWF X βo + XtWFUF

· = X

tWF (X βo + UF)

(X

tWF X + St

S) β = X

tWF ZF , (B-28)

38
Approved for public release; distribution is unlimited.

where

ZF = X βo + UF . (B-29)

These are the normal equations for minimization of Q = kZF − X βk

2

WF

+ kS βk

2

.

At each step, the normal equations are solved as the augmented system of Section A.3. P-spline models are obtained with a B-spline basis as in Section A.4, and

smoothing spline models are obtained with X = I as in Section A.5.

Figure B-1 depicts the effect of smoothing parameter variation for the P-spline

model, and Fig. B-2 depicts the effect of smoothing parameter variation for the

S-spline model. See Fig. B-3 for graphs of OCV, GCV, and AIC as functions of λ,

and see Fig. B-4 for GVC-optimal P-spline and S-spline solutions.

39
Approved for public release; distribution is unlimited.

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −7.7

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −7.2

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −6.7

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −6.2

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −5.7

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −5.2

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −4.7

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −4.2

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −3.7

Fig. B-1 P-spline optimization, GLM

40
Approved for public release; distribution is unlimited.

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −7.8

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −7.3

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −6.8

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −6.3

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −5.8

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −5.3

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −4.8

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −4.3

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

P

λ = −3.8

Fig. B-2 S-spline optimization, GLM

41
Approved for public release; distribution is unlimited.

−8 −7 −6 −5 −4 −3 −2

78

80

82

84

86

λ

aic, P−spline

−8 −7 −6 −5 −4 −3 −2

0.0033

0.0034

0.0035

0.0036

0.0037

0.0038

λ

gcv, P−spline

−8 −7 −6 −5 −4 −3 −2

13.0

13.5

14.0

14.5

15.0

15.5

λ

ocv, P−spline

−8 −7 −6 −5 −4 −3 −2

76

78

80

82

84

86

λ

aic, S−spline

−8 −7 −6 −5 −4 −3 −2

0.0032

0.0033

0.0034

0.0035

0.0036

0.0037

0.0038

λ

gcv, S−spline

−8 −7 −6 −5 −4 −3 −2

13.0

13.5

14.0

14.5

15.0

15.5

λ

ocv, S−spline

Fig. B-3 Smoothing parameter selection, GLM

42
Approved for public release; distribution is unlimited.

500 1000 1500

0.0

0.2

0.4

0.6

0.8

1.0

x

y

S−spline

P−spline

Fig. B-4 GCV-optimal fits, GLM

43
Approved for public release; distribution is unlimited.

B.4 Canonical Link

For the canonical link, wFi = g(ηi) = v(μi) and uFi = (yi − μi)/wFi. Also, since

dθi/dβ = xi and d

2

θi/dβ βt = 0, it follows that E H ( β) = H ( β). So NewtonRaphson and Fisher scoring are equivalent.

B.5 Confidence Intervals

Normal-approximation 100c% confidence intervals on the mean response are given

by

G



x β ± Φ1−(1−c)/2

√

x

tV x

, (B-30)

where Φ is a standard normal quantile, V is the estimated parameter variance matrix,

and x is a row of an X matrix corresponding to the desired level. For the basis

implementation, this is

x =



f1(v),. . . , f p(v)



. (B-31)

B.6 Bernoulli Response

For example, suppose the response Y ∈ {0,1} is Bernoulli with Pr[Y = 1] = μ =

1 − Pr[Y = 0]. The Bernoulli distribution is a member of the exponential family,

Eq. B-2, with the particular form

f (y) = μ

y

(1 − μ)

1−y = exp "

y log μ

1 − μ

+ log(1 − μ)

#

. (B-32)

So a = 1, c = 0, and there is no nuisance parameter. Furthermore, θ = log μ/(1 − μ)



and μ = 1/(1 + e

−θ

), and so b(θ) = − log(1 − μ) = log(1 + e

θ

). Note that

E[Y] = b

0

(θ) = e

θ

/(1 + e

θ

) = μ and Var[Y] = v(μ) = b

00(θ) = e−θ

/(1 + e

−θ

)

2 =

μ(1 − μ) as expected.

With η = θ and μ = G(θ), we see that the canonical link for Bernoulli response is

the logistic cumulative distribution function (CDF) G(η) = 1/(1 + e

−η

) = μ. Note

that g(η) = μ(1− μ) = v(μ), so wFi = μi(1− μi) and uFi = (yi − μi)/



μi(1 − μi)



.

The resulting model is logistic regression, or the logit model.

For an arbitrary link CDF G, we take η = X β, μ = G(X β), v(μ) = μ(1 − μ), and

g(η) = g(X β).

Use of the standard normal CDF G = Φ with probability density function (PDF)

44
Approved for public release; distribution is unlimited.

g = φ gives the probit model.

Because the likelihood function is L =

Q μ

yi

i

·



1 − μi

1−yi

, we have Lfull = 1 for

the Bernoulli model, and then the deviance is ∆ = −2 log L.

As an example, consider the usual 2-parameter model with predictor (x1,. . . , xn)

and response (y1,. . . , yn). The increment δ = (d0, d1) is the solution of Mδ = A,

where

M = X

tW X =













P

wi

P

wi xi

P

wi xi

P

wi x

2

i













and A = X

tWU =













P

wi ui

P

wi ui xi













. (B-33)

To do the Fisher update of Section B.2.5, calculate the linear response ηi = b0+b1 xi,

mean μi = G(ηi), derivative gi = g(ηi), variance vi = μi(1 − μi), transformed response ui = uFi = (yi − μi)/gi

, weight wi = wFi = g

2

i

/vi, and weighted transformed

response wiui = wFi uFi = (yi − μi)gi/vi.

For the Newton-Raphson update, wDi = (yi − μi)(g

0

i

vi − g

2

i

v

0

i

)/v

2

i

and wNi = wFi −

wDi and ui = uNi = wFi/wNi(yi − μi)/gi. Then wi = wNi and wiui = wNi uNi =

wFi uFi.

For the canonical link wi = wFi = gi = vi and wiui = wFi uFi = yi − μi.

45
Approved for public release; distribution is unlimited.

Intentionally left blank.

46
Approved for public release; distribution is unlimited.

Appendix C. B-splines

47
Approved for public release; distribution is unlimited.

A degree-d B-spline basis (B1,d,. . . , Bn,d ) of dimension n uses n + d + 1 knots

(t1,. . . ,tn+d+1).

There are k = n − d − 1 internal knots (internal to the data) and 2d + 2 boundary

knots, d + 1 on each side. There are no internal knots if k = 0 and n = d + 1, and in

general k > 0 so n > d. The left knots are (t1,. . . ,td+1), the internal knots, when

n > d + 2, are (td+2,. . . ,tn), and the right knots are (tn+1,. . . ,tn+d+1).

On the data range [x0, x1] internal knots are evenly spaced so ti = x0+

i−d−1

k+1

(x1−x0)

for i = d + 2,. . . ,n. This also accounts for td+1 = x0 and tn+1 = x1.

Compact knots are constructed by replicating knots at the boundary, so t1 = · · · =

td+1 = x0 and tn+1 = · · · = tn+d+1 = x1.

Uniform knots are constructed by repeating the uniform spacing for all knots, so

ti = x0 +

i−d−1

k+1

(x1 − x0) for i = 1,. . . ,n + d + 1.

This is the Cox-de Boor recursion for B-spline basis functions:

Bi,0(x) =









1 , ti 6 x < ti+1

0 , otherwise

,

i = 1,. . . ,n + d ;

Bi,j(x) =

x − ti

ti+j − ti

Bi,j−1(x) +

ti+j+1 − x

ti+j+1 − ti+1

Bi+1,j−1(x) ,

j = 1,. . . , d , i = 0,. . . n + d − j . (C-1)

This is the relation for derivatives of B-spline basis functions:

d

dx

Bi,j(x) =

j

ti+j − ti

Bi,j−1(x) −

j

ti+j+1 − ti+1

Bi+1,j−1(x) . (C-2)

Example graphs follow. All spline bases have dimension n = 8. Figures C-1 and

C-2 demonstrate uniform and compact knots, respectively, for bases of degree d =

0,1,2, and 3. Figures C-3 and C-4 demonstrate uniform and compact knots, respectively, for bases of degree d = 3 with derivative orders p = 0,1,2, and 3.

48
Approved for public release; distribution is unlimited.

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.2

0.4

0.6

0.8

1.0

x

y

n=8, d=0, p=0, u=1

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.2

0.4

0.6

0.8

1.0

x

y

n=8, d=1, p=0, u=1

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.2

0.4

0.6

x

y

n=8, d=2, p=0, u=1

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.1

0.2

0.3

0.4

0.5

0.6

x

y

n=8, d=3, p=0, u=1

Fig. C-1 B-spline basis, uniform

49
Approved for public release; distribution is unlimited.

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.2

0.4

0.6

0.8

1.0

x

y

n=8, d=0, p=0, u=0

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.2

0.4

0.6

0.8

1.0

x

y

n=8, d=1, p=0, u=0

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.2

0.4

0.6

0.8

1.0

x

y

n=8, d=2, p=0, u=0

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.2

0.4

0.6

0.8

1.0

x

y

n=8, d=3, p=0, u=0

Fig. C-2 B-spline basis, compact

50
Approved for public release; distribution is unlimited.

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.1

0.2

0.3

0.4

0.5

0.6

x

y

n=8, d=3, p=0, u=1

0.0 0.2 0.4 0.6 0.8 1.0

−3

−2

−1

0

1

2

3

x

y

n=8, d=3, p=1, u=1

0.0 0.2 0.4 0.6 0.8 1.0

−40

−20

0

20

x

y

n=8, d=3, p=2, u=1

0.0 0.2 0.4 0.6 0.8 1.0

−400

−200

0

200

400

x

y

n=8, d=3, p=3, u=1

Fig. C-3 B-spline basis derivatives, uniform

51
Approved for public release; distribution is unlimited.

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.2

0.4

0.6

0.8

1.0

x

y

n=8, d=3, p=0, u=0

0.0 0.2 0.4 0.6 0.8 1.0

−15

−10

−5

0

5

10

15

x

y

n=8, d=3, p=1, u=0

0.0 0.2 0.4 0.6 0.8 1.0

−200

−100

0

100

x

y

n=8, d=3, p=2, u=0

0.0 0.2 0.4 0.6 0.8 1.0

−1000

−500

0

500

1000

x

y

n=8, d=3, p=3, u=0

Fig. C-4 B-spline basis derivatives, compact

52
Approved for public release; distribution is unlimited.

Appendix D. Matrix Differentiation Operators

53
Approved for public release; distribution is unlimited.

Consider that y = f (x), so y is a function of x, and n values of zi = (xi, yi) for

i = 1,. . . ,n are given, where each xi < xi+1. Let f k,p be the degree-p interpolating

polynomial through the p+1 consecutive points (zk ,. . . ,zk+p), so that f k,p(xi) = yi

for k 6 i 6 k + p. Interpolating polynomials satisfy the recursion

f k,p(x) =

(x − xk ) f k+1,p−1(x) + (xk+p − x) f k,p−1(x)

xk+p − xk

. (D-1)

Define the Newton basis polynomials nk,j by

nk,0(x) = 1 , 1 6 k 6 n

nk,j(x) =

k

Y+j−1

i=k

(x − xi) , 1 6 k 6 n , 1 6 j 6 n − k , (D-2)

so nk,j(xi) = 0 for j > 1 and k 6 i 6 k + j − 1. Note that nk,j has degree j and the

leading coefficient 1, so nk,j(x) = x

j+ lower-degree terms. Write the polynomial as

f k,p(x) =

X

p

j=0

ak,jnk,j(x) (D-3)

so that

f k,0(x) = ak,0nk,0(x) , 1 6 k 6 n

f k,j(x) = ak,jnk,j(x) + f k,j−1(x) , 1 6 j , (D-4)

and the leading coefficient of f k,jis ak,j. Then, because of Eqs. D-1 and D-4, the

coefficients obey

ak,0 = yk

ak,j =

ak+1,j−1 − ak,j−1

xk+j − xk

. (D-5)

54
Approved for public release; distribution is unlimited.

The coefficient sequence is



ak,0 , ak,1 ak,2 , ak,3 , ak,4 , · · · =

 

yk ,

yk+1 − yk

xk+1 − xk

,

yk+2−yk+1

xk+2−xk+1

−

yk+1−yk

xk+1−xk

xk+2 − xk

,

yk+3−yk+2

xk+3−xk+2

−

yk+2−yk+1

xk+2−xk+1

xk+3−xk+1

−

yk+2−yk+1

xk+2−xk+1

−

yk+1−yk

xk+1−xk

xk+2−xk

xk+3 − xk

,

yk+4−yk+3

xk+4−xk+3

−

yk+3−yk+2

xk+3−xk+2

xk+4−xk+2

−

yk+3−yk+2

xk+3−xk+2

−

yk+2−yk+1

xk+2−xk+1

xk+3−xk+1

xk+4−xk+1

−

yk+3−yk+2

xk+3−xk+2

−

yk+2−yk+1

xk+2−xk+1

xk+3−xk+1

−

yk+2−yk+1

xk+2−xk+1

−

yk+1−yk

xk+1−xk

xk+2−xk

xk+3−xk

xk+4 − xk

, · · · !.

(D-6)

For evenly spaced x j with x j+1 − x j = h and xk+j − xk = j h, we have

ak,0 = yk

ak,j =

ak+1,j−1 − ak,j−1

j h

, (D-7)

so



ak,0 , ak,1 ak,2 , ak,3 , ak,4 , · · · =

 

yk ,

yk+1 − yk

h

,

yk+2 − 2yk+1 + yk

2h

2

,

yk+3 − 3yk+2 + 3yk+1 − yk

6h

3

,

yk+4 − 4yk+3 + 6yk+2 − 4yk+1 + yk

24h

4

, · · · !. (D-8)

For evenly spaced x, the coefficients are

ak,j =

1

j!h

j

X

j

i=0

(−1)

j+i

 

j

i

!

yk+i. (D-9)

In general, the numerical derivative of order p coincides with the p

th derivative f

(p)

k,p

of f k,p,

f

(p)

k,p

(x) = p! ak,p. (D-10)

There are n − p sequences (zk ,. . . ,zk+p) and therefore n − p values of f

(p)

k,p

for

1 6 k 6 n − p.

Now, express the differentiation operator as a matrix.

55
Approved for public release; distribution is unlimited.

Define ν(j), the first-order difference matrix of dimension k > 2, size (j − 1) × j,

by ν(j)i,i = −1 and ν(j)i,i+1 = 1 for 1 6 i 6 j − 1.

ν(j) =



























−1 1

−1 1

.

.

.

.

.

.

−1 1



























. (D-11)

Then

ν(j)





















u1

.

.

.

uj





















=





















u2 − u1

.

.

.

uj − uj−1





















, (D-12)

and ν serves to evaluate the numerator of Eq. D-5.

Define δ(u, j), the lag-j (square) diagonal difference matrix of dimension n − j for

u = (u1,. . . ,un) and 1 6 j 6 n − 1 by δ(u, j)i,i = uj+i − ui for 1 6 i 6 n − j

δ(u, j) =





















uj+1 − u1

.

.

.

un − un−j





















. (D-13)

So δ serves to evaluate the denominator of Eq. D-5.

Then matrix derivative operators for the n-vector x are given recursively by

D1 = δ(x,1)

−1

ν(n)

Dj = j δ(x, j)

−1

ν(n − j + 1) Dj−1, (D-14)

where the factor of j serves to evaluate the factorial in Eq. D-10.

The matrix DP evaluates the p

th derivative.

Based on points (zk ,. . . ,zk+p) for 1 6 k 6 n− p, the individual derivative value are

f

(p)

k,p

(x) = (Dp y)k . (D-15)

56
Approved for public release; distribution is unlimited.

The derivative vector is Dp y.

Dp y =





























f

(p)

1,p

(x)

f

(p)

2,p

(x)

.

.

.

f

(p)

n−k,p

(x)





























, (D-16)

and its norm is approximately





f

(p)

(x)

2

dx ' kDp yk

2 = ytDt

pDp y . (D-17)

R language code for evaluating Dp is easy.

D <- diff(diag(n)) / diff(x, lag=1)

if (p>1) for (j in 2:p)

D <- j * (diff(diag(n - j + 1)) / diff(x, lag=j)) %*% D

Here, diff(diag(n)) is ν(n) and diff(x, lag=j)) is δ(x, j).

57
Approved for public release; distribution is unlimited.

List of Symbols, Abbreviations, and Acronyms

AIC Akaike information criterion

CB Chang-Bodt

CDF cumulative distribution function

GCV generalized cross-validation

GLM Generalized Linear Model

IID independent and identically distributed

IRLS iteratively reweighted least squares

OCV ordinary cross-validation

PDF probability density function

QR quantal response

SSE sum of squared errors

58
Approved for public release; distribution is unlimited.

1

(PDF)

DEFENSE TECHNICAL

INFORMATION CTR

DTIC OCA

2

(PDF)

DIRECTOR

US ARMY RESEARCH LAB

RDRL CIO L

IMAL HRA MAIL & RECORDS MGMT

1

(PDF)

GOVT PRINTG OFC

A MALHOTRA

1

(PDF)

DIR US ARMY EVALUATION CTR

HQ

TEAE SV

P A THOMPSON

1

(PDF)

DIR USARL

RDRL SL

P BAKER

20

(PDF)

DIR USARL

RDRL DPW

T HOLDREN

RDRL SL

D BAYLOR

P DISALVO

T STADTERMAN

RDRL SLB

R BOWEN

G MANNIX

RDRL SLB D

J COLLINS

J EDWARDS

R GROTE

L MOSS

E SNYDER

RDRL SLB E

M MAHAFFEY

RDRL SLB G

N ELDREDGE

RDRL SLB S

R DIBELKA

C KENNEDY

M PERRY

R SAUCIER

G SAUERBORN

RDRL SLB W

S SNEAD

RDRL SLE

R FLORES

59
Intentionally left blank.

60