ARMY RESEARCH LABORATORY F 

Modeling Nonmonotonic 

Dose-Response Curves 

Quinton J. Nottingham, Jeffrey B. Birch, and Barry A. Bodt 

ARL-TR- 197.2 January 200 1 

Approved for public release: distribution unlimited. 
. 

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. 

DESTRUCTION NOTICE-Destroy by any method that 

will prevent disclosure of contents or reconstruction of this 

document. 
Army Research Laboratory 

Aberdeen Proving Ground, MD 21005-5425 

ARL-TR-1972 January 2001 

Modeling Nonmonotonic 

Dose-Response Curves 

Quinton J. Nottingham, Jeffrey B. Birch 

Virgina Polytechnic Institute & State University 

Barry A. Bodt 

Information Science and Technology Directorate, ARL 

Approved for public release; distribution unlimited. 
Abstract 

A number of procedures have been used to analyze nonmonotonic binary 

data to predict the probability of response. Some classical procedures are 

the Up and Down strategy, the Robbins-Monro procedure, and other sequential optimization designs. Recently, nonparametric procedures such as 

kernel regression and local linear regression have been applied to this type 

of data. It is well known that kernel regression has problems fitting the data 

near the boundaries, and a drawback with local linear regression is that it 

may be too linear when fitting data from a curvilinear function. This report 

introduces a procedure called local logistic regression, which fits a logistic 

regression function at each of the data points. United States Army projectile 

data are used in an example that supports the use of local logistic regression for analyzing nonmonotonic binary data for certain response curves. 

Properties of local logistic regression are presented along with simulation 

results that indicate some of the strengths of the procedure. 

ii 
Contents 

1 Introduction 

2 Possible Models 

2.1 Parametric Models . . . . . . . . . . . . . . . . . . . . . . . . 

2.2 Nonparametric Models . . . . . . . . . . . . . . . . . . . . . . 

3 Some Simulation Results 

4 Discussion and Summary 

Bibliography 

Appendices 

A Numerical Considerations 

B Program Documentation & Comments 

Distribution 

Report Documentation Page 

Figures 

1 

3 

3 

5 

11 

13 

14 

17 

17 

19 

33 

35 

1 Army macrocomposite data, illustrating shatter gap phenomenon 2 

2 Linear logistic regression fit and local linear regression fit to 

macrocomposite data . . . . . . . . . . . . . . . . . . . . . . . . . 4 

3 Comparison of fits between llogrs and cdf mixture model . . . . 5 

4 Local logistic regression fit to macrocomposite data with 95% 

confidence bands . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 

5 Local logistic regression fit to macrocomposite data using linear 

predictor for confidence interval . . . . . . . . . . . . . . . . . . 9 

6 Curves 1 to 3 from table 1 generated for Monte Carlo simulation 11 

..* 

111 
Tables 

1 Parameter values for simulated curves ............... 11 

2 Lower V& estimation summary ................... 12 

3 Upper Vf0 estimation summary ................... 12 

iv 
1. Introduction 

Consider modeling the dose-response curve for the ungrouped binary response variable y. This dose-response curve represents the probability that 

y is 1 for each given value of y dose,” the single regressor 21. In many applications in such areas as biology, industry, and business, this response 

curve may be modeled successfully by a monotonic parametric function 

such as the normal or logistic cumulative distribution function (cdf). However, in this report, we present an application where the response curve is 

nonmonotonic, rendering traditional methods such as logistic regression 

inappropriate. We suggest the use of local logistic regression (llogr), a nonparametric method that simultaneously enables modeling a nonmonotonic 

dose-response curve while maintaining the restriction that all estimated 

responses take values between 0 and 1. The following real-life example 

demonstrates the flexibility and applicability of this new procedure. 

An experiment of importance to the development of a kinetic energy penetrator (a a projectile”) and the armor to resist it is to model the probability 

that the projectile will penetrate the plate of armor as a function of the velocity of the projectile. Normally, one expects the probability of penetration 

to increase as velocity increases. Routinely, testers are asked to estimate the 

V&, the velocity where the probability is 0.50 that the projectile will penetrate the armor. Experimenters use variants of the Up and Down strategy 

(Dixon and Mood, 1948), the Robbins-Monro (Robbins and Monro, 1951), 

or other sequential designs to gather observations, and the maximum likelihood estimate (MLE) for the mean (V&J is formed with a logistic or normal response function. An interlaboratory study of V& estimation revealed 

widely varying estimates for some penetrator-armor matches (Chang and 

Bodt, 1997). Initially, the blame was thought to lie with varying test procedures among experimental testing sites. Recent careful study of the data 

suggests a phenomenon at work here that one material scientist refers to as 

the * shatter gap.” 

The shatter gap is not precisely defined, and the exact physical mechanism 

for significant bullet shattering is unknown. Simply put, when bullet shattering occurs against certain armor materials, kinetic energy is diffused, 

thereby reducing the ability of the penetrator to defeat the armor. This reduction is reflected over a range of velocity in which a decrease in response 

probability accompanies increasing velocity Eventually, as velocity continues to increase, the kinetic energy developed is so much an overmatch that 

the probability of penetration rises again. One working definition, suitable 

for data currently available, characterizes the shatter gap in terms of an upper and lower V&, V~OU, and V&L, respectively. In those data, penetrator 

shatter does not cause the probability of penetration to decrease until after 

1 
a V~O is achieved. A second, higher, 1 $0 is achieved by partially intact penetrators with increased velocity The interval they bound is thought of as the 

shatter gap. 

These ideas are illustrated in figure 1, showing the results of firing 303 

projectiles at macrocomposite armor with velocities ranging from 1432 to 

2773 ft/s. The response variable is whether the projectile penetrated the armor (scored as a * 1”) or failed to penetrate the armor (scored as a * 0” ). 

Typical of penetrator-against-plate response curves, low velocities result in 

a low proportion of penetrations, and high velocities result in a high proportion of penetrations. Unusual is that over intermediate velocities, the 

proportion of penetrations first rises with velocity, then falls, and then rises 

once more. Note especially the response activity between 2100 and 2300 

ft/s relative to response activity for neighboring velocities. 

The shatter gap phenomenon is more clearly seen from the curve in figure 

1; this curve (obtained by the method presented in sect. 2) represents the 

estimated probability of penetration as a function of velocity. 

We can represent the Vs:30 graphically by extending a line at the probability 

of 0.5 over to the estimated probability curve and dropping lines down to 

the velocity axis at the intersection of the line and the curve. The interval 

bounded by these two 1+,0’s on the velocity axis is the shatter gap. Figure 1 

shows that the shatter gap extends from V50f, of approximately 1713 ft/s to 

V&T of approximately 2510 ft/s, a distance of 797 ft/s, far larger than the 

usual error of approximately 100 ft/s. 

In section 2 we discuss several analytical methods to deal with data of this 

type. 

Figure 1. Army l- . -- . . . . . . . ". . .-.-... I 

macrocomposite data /- 

(solid dots), illustrating 0.9. ,’ 

shatter gap 

o 8 

phenomenon. Curve is 0.7 

estimated probability of g 0.6. 

penetration. s 8' \ 

0.5.. 

‘\ $$ 

E 0.4. 

0.3. <\ 

0.2. . , 

0.1 ,j/' 

0, .' 

1300 1500 V50~ 1900 2100 2300 V5ou 2700 2900 

Velocity (fth) 

2 
2. Possible Models 

An obvious approach to fitting 0 to 1 response data is to use the standard 

logistic regression model. However, as will be demonstrated, such a model 

will fail to identify the shatter gap in the penetration data because of the inherent monotonically increasing (or decreasing) nature of the logistic curve. 

Two more promising possible models for fitting data of the type given in 

the curve in figure 1 are the cdf mixture model (CDFMM) (Bodt and Chang, 

1997), a parametric method model procedure, and local Ilogr, a nonparametric method. 

2.1 Parametric Models 

Consider the model 

yi = P(q) + &i, i = 1, . ..) n ) 

where y~i represents the response variable, P(Q) is the penetration probability that is a function of the velocity vi at the ith experimental run, and Q 

is the random error term. The response is recorded as a * 1” if armor penetration has occurred and a y 0” otherwise. The errors are assumed to have 

expectation of zero. 

In the logistic regression model, the form of P (t)i) is based on the logistic 

cumulative distribution function, where 

P(Vi) = (1 + exp [- (PO + p1Vi)])-I = (1 + exp (-Q))-’ = F (2)#), 

where gi = (1 vi) and p = 

0 

PO where underlining indicates a vector - Pl ’ 

and a prime indicates a transpose. 

The method of maximum likelihood is used to obtain estimates (MLEs), 

of the unknown parameters ,& and pi. To obtain the MLEs, the method 

requires an iterative procedure, for example, the Gauss-Newton method 

using iterated reweighted least squares (IRLSs). Details concerning logistic 

regression may be found in McCullagh and Nelder (1983), Myers (1990) 

and Ryan (1997), among many others. 

While the logistic regression fit guarantees an estimated response between 

0 and 1, the fit, as seen in figure 2 (dashed line), is entirely inadequate to 

the armor penetration data- resulting in a nearly flat line fit, completely 

missing the up and down nature of the response seen in the llogr fit (solid 

line). 
Figure 2. Linear logistic 1 .I 

regression fit (dashed . -.--..- *.. . 

line) and local linear 0.9 

regression fit (bold solid 

line) to macrocomposite 

data. 0.7. 

% 

s 

a 0.5. 

: 

u 

03 

0.1 

/ 

-0.1 

1300 1500 1700 1900 2100 2300 2500 2700 2900 

Velocity (ftk) 

The CDFMM approach, using mixtures of three cdfs, models the penetration probability P(,u) as a function of velocity Y as 

where FT. FI, and F, represent appropriately chosen cdfs for the transition, 

intact, and shattered mechanisms. 

One interpretation of the above model is to consider P(~J) as the probability of response due to the mixture of two penetration mechanisms. The 

first mechanism is penetration from an intact (I) projectile. The second 

mechanism is penetration from a shattered (S) projectile. The penetrator 

being intact or shattered is also a function of velocity in this formulation, 

and so the proportions of the two penetration mechanisms in the mixture 

will vary also with velocity. It is assumed in the model that the probability 

of projectile shatter will increase monotonically. This proportion change is 

modeled as the transition (T) from the intact mechanism to the shattered 

mechanism. 

The CDFMM model fit to the Army penetration data is illustrated in figure 

3, which also shows for comparison the llogr fit. 

The dashed line represents the CDFMM fit using logistic cdfs with estimated means of 1972, 1664, and 2523 and estimated standard deviations 

of 186, 56, and 46 for the transition, intact, and shattered mechanisms, respectively The curves are obviously quite similar, with the CDFMM fit 

resulting in slightly larger (at around 1800 ft/s) and slightly smaller (at 

around 2400 ft/s) penetration probabilities. The shatter gap obtained with 

the CDFMM fit is only slightly wider than that obtained with the llogr fit. 

The estimates are the MLEs for this nonlinear model with binary response, 

an example of a generalized nonlinear model. 

One problem with the CDFMM approach is that six parameter estimates 

are required to fit the model. This requires considerable effort to achieve 

convergence. We used the SAS Proc NLIN (SAS Version 7) with an appropriately iterated weight matrix to obtain the fit for the Army penetration 

4 
Figure 3. Comparison of 1 

fits between llogrs (solid 

line) and cdf mixture 

O 9 

’ 

model (dashed line). 0.8 

0.7 

g 0.6 

z CL 0.5 

z 

cc 0.4 

0.3 

0.2 

0.1 

0 .I1300 1500 1700 1900 '2100 2300 2500 2700 2900 

Velocity (ft/s) 

data. Our experience suggests that this procedure is extremely sensitive to 

the starting values, often failing to converge if the starting values are only 

slightly different from the values at convergence. Additionally, the resulting fit, as can be seen in figure 3, is somewhat inflexible, resulting in fits that 

look like a mixture of logistic cdfs, even when the data do not exhibit such 

a pattern. That is, the curve may appear to be too smooth, as illustrated in 

figure 3, where the llogr method seems to capture the wiggle resulting from 

the cluster of data around 2250 ft/s. 

2.2 Nonparametric Models 

Another approach to modeling these data is to use nonparametric regression techniques. Such methods as kernel regression and local linear regression can trace out smooth curves with a multitude of peaks and valleys. The 

shatter gap effect that must be captured in the armor penetration example 

could then be identified by such y local” fi tting techniques. 

Consider again the model 

yi = P(q) +&i, i = l)...) n, 

where P(vi) is no longer specified by a parametric function such as the 

CDFMM function given in section 2.1. It is assumed only that P(Q) is some 

arbitrary smooth function. If E(Q) = 0, then E(yi) = P(vi) = Pi, the armor penetration probability at velocity vi. It follows that $i = P(vi) = & 

is the estimated penetration probability at velocity vi. Kernel regression is 

a method of approximating P(Q) (where ~0 is any arbitrary velocity including Q, the velocity for the ith observation) using a weighting sequence 

on the response variable, where the weights are functions of the distances 

between the values of the regressor variable and ~0. One form of the kernel weighting sequence, proposed by Nadaraya (1964) and Watson (1964), 

assigns weights of h& to yj when estimating the response at ZIO by 

5 
where I<() represents some appropriate kernel function and b is bandwidth. 

The kernel estimate of response at ?I,, p:, is given by a weighted average of 

the observed yi’s as 

where I,$ (I+ replaces ~~~ in the /I;~ expression above) are the kernel weights 

on the observations yj for estimating the mean response at velocity 11; using 

the kernel function I<(). j = 1.2. . . . . 71. 

All n estimates of mean response may be expressed as 8” = H”;lg, where 

6” is the n x 1 vector of kernel predictions and H” is the 71 x II matrix of - 

weights with rows bik’ = (11:~ 11:~ . . . /$r,), for i = l,..., 71,. The matrix 

H’, called the kernel hat matrix, plays a role in kernel regression similar to 

that of the ordinary least square (OLS) ” hat” matrix H in linear regression. 

The kernel estimate of P(u~) depends on the choice of the kernel function K() and the bandwidth 0. For example, one popular choice for kernel 

function is the simple Gaussian kernel defined as I<(U) = cxp (-I?). The 

method for selecting the bandwidth is extremely critical. The magnitude of 

the bandwidth determines the smoothness, or lack thereof, of the regression function. One of the more commonly used methods, referred to as the 

method of cross-validation, selects the bandwidth to minimize the PRESS 

(Allen, 1974) statistic. A related method finds b by minimizing a penalized 

PRESS statistic, called PRESS*(b), given by 

? (Yi - p;-i (b)) 2 

PRESS*(b) = ‘=l 

n - t,r [H”] ’ 

where p2-i (b) is the a minus-i” predicted penetration probability based on 

kernel regression ZJ~ for the current value of b with the ith observation removed, and tr [H”] is the trace of the n x n kernel weight matrix H”. Since 

tr [Hk] reflects the kernel fits’ ” model degrees of freedom” (Cleveland, 

1978), it is seen that the denominator of PRESS*(b) penalizes the PRESS 

statistic for choosing b too small. Empirical studies by the authors and others (Einsporn and Birch, 1993; Mays, 1995) have demonstrated that using 

PRESS*(b) is often superior to using PRESS as a bandwidth selector. 

The kernel regression fit to the armor penetration data results in a smooth 

curve that follows the up and down pattern exhibited in figures 1 and 3. 

While this is a desirable characteristic, the resulting curve is not appealing, 

since it is entirely possible that some of the curve’s fitted values may lie outside the 0 to 1 range, a natural restriction resulting from estimating the penetration probabilities. Indeed, the kernel regression fit to the Army macrocomposite data results in fits at the lower left that are less than zero and fits 

to the upper right that exceed one. This problem is not avoided by the use of 

more complicated nonparametric smoothers, such as local linear regression 

(Hardle, 1990), where each fit at v = 110 is obtained locally by a weighted 

6 
simple linear regression model with the h&‘s serving as the weights, or local polynomial regression (Fan and Gijbels, 1996). These methods also do 

not restrict the resulting fits to between 0 and 1, a requirement when modeling probabilities. The local linear regression fit to the macrocomposite data 

can be seen in figure 2. 

However, the logistic regression parametric method can be combined with 

the nonparametric concept of local fitting to produce a smooth curve, flexible enough to capture the up and down patterns exhibited in the armor 

penetration data and giving fits that are between zero and one. We propose 

local (linear) logistic regression where at each velocity 21 = ~0, a weighted 

linear logistic fit is obtained with h&‘s serving at the weights, in exactly the 

same manner as in local linear regression. That is, the fit at ‘u = ~0 would 

be obtained as 

p (vo) = (1 + exp [- (boo + 1;110~0)])-~ = (1 + exp [- (L@~)])-~. 

The estimated coefficients jo, which change values for each ‘u = ~0, are 

obtained via the same IRLS algorithm referred to earlier for logistic regression, with a slight adjustment to the weight matrix. Thus, one step of the 

algorithm would compute the updated value of B. as 

jol = joo + p*ox> -l x’J@oyT, , 

where the n x n diagonal matrix I/i/o has elements 

We see that llogr follows a weighting scheme that combines the logistic 

weights with the kernel weights in a local manner. The n x 1 vector $, has 

elements 

the locally adjusted response to velocity wj. Upon convergence of the IRLS 

algorithm, & is set equal to loi, and the estimated response is obtained as 

F (w,,) = F (i@,). 

Approximate inferential information can be obtained by straightforward 

application of the variance operator, the multivariate delta method, and the 

asymptotic distribution of the MLEs. For example, it can be shown that the 

asymptotic variance of lo is var ( > i. = (x’wox)-l (X’WlX) (x’wox)-l, 

where the diagonal elements of the n x n diagonal matrix IV0 are wajj = 

F (Qo) (1 - F (~;p~)) h&, and the diagonal elements of the n x n diagonal matrix VVi are wijj = wajj h$. Using this result, one can express the 

asymptotic variance of each fit at 21 = ~0 by applying the delta method as 

var (F (z&$0)) = f2 (&PO) vb var (eo) ‘uo. 

7 
It follows then from the asymptotic distribution of F that one form 

of the approximate (1 - o) 100% confidence interval for P (~0) is 

where the standard error of F ( > IJ,~,, is the square root of VHY F *II’ r> ( (4,)) 1 

with I0 replacing & in this variance expression. 

While this approach for establishing confidence intervals is a straightforward application of Wald inference, the resulting interval, (LCL, UCL), may 

result in lower confidence limit (LCL) less than zero or an upper confidence 

limit (UCL) greater than one. To avoid this problem, we can consider a second approach where we obtain the (1 - 0) 100% confidence interval for the 

linear predictor zJ,jO by once again applying the delta method. It can be 

shown that the asymptotic variance of the linear predictor Q?,, is 

- 

with resulting (1 - (1) 100% confidence interval as 

&ii, f Zl-” se I@” = (LCL*. ucL*). 

- 2 ( > 

It follows that an approximate (1 - 0) 100% confidence interval for F ( > &/!I,, 

- 

is (F (LCL*) : F (UCL*)). 

The general procedure is to obtain the fit and corresponding (1 - 0) 1OO’X 

confidence interval (by either method) for a fine grid of values of lie throughout the range of velocity values. Connecting the fits, lower confidence limits, and upper confidence limits with straight line segments results in three 

smooth curves, the curve fit to the data and the two (1 - a) 1OO’X confidence bands. Both types of interval methods are illustrated in figures 4 and 

5. The method based on the linear predictor results in unsatisfactory wide 

intervals at the boundaries. 

Our procedure is a specific example, though developed independently, of 

the related method, local polynomial kernel regression of generalized linear models, introduced by Fan, Heckman, and Wand (1995). 

Once the response curve and the confidence bands are obtained, the V~O 

as well as confidence intervals for the V~O can be computed. Because it is 

not possible to obtain a closed functional form for the estimate of P (vo), 

the V~O and confidence intervals for the Vsc, are obtained through the y inverse” process, where a horizontal line is extended from P(V) = 0.5 to 

the three curves, representing the upper confidence limit, the estimated response curve, and the lower confidence limit. At the intersection, vertical 

lines are dropped to the velocity axis, where the V,O and the 95-percent 

confidence interval for the V~O, (1% 01,. V&p), can be obtained. (We note that 

if the shatter gap effect is present, there will be two I$” values, each with 

8 
Figure 4. Local logistic 

regression fit to 

macrocomposite data 

with 95% confidence 

bands. 

0.9. 

0.8. 

0.7. ,:' x-'- . ..___. 

0.6. 

-0.1 I 

1300 1500 1700 1900 2100 23bO 2500 2700 2900 

Velocity (Ws) 

Figure 5. Local logistic 

regression fit to 

macrocomposite data 

using linear predictor 

for confidence interval. 

l0.9. 

0.8. 

0 7. 

gz 0.6. 

2 0.5. 

g 0.4. 

f300 1500 '1700 1900 2100 2300 2500 2700 2900 

Velocity (ftk) 

9 
a corresponding 9%percent confidence interval.) This process is easily accomplished numerically with a computer. The IRLS algorithm to compute 

the llogr fit is outlined in appendix A. A SAS macro for calculating the llogr 

fit to any set of 0 to 1 data may be obtained from the authors. The complete 

SAS macro is listed in appendix B. 

10 
3. Some Simulation Results 

In this section, we present a small simulation study, generating binary data 

from several nonmonotonic dose-response curves. Because the primary interest with the Army penetration data is estimation of the V& values, we 

evaluate the llogr procedure by comparing the estimated V& values with 

the true Vi0 values. The true cdf used in this evaluation was chosen as the 

cdf mixture model: 

The values of ,UT, 1-11, and PS were varied while fST, a~, and as were set 

equal to each other at a constant value of 0.05. The cdf representing FT, F1, 

and Fs was chosen to be the logistic cdf. For example, the cdf representing FT can be written as FT = (1 + exp (- (X - /LT) /uT))-~. While our full 

simulation study considered a variety of values of PT, ,UI, and ps, for the 

sake of brevity we show only the results for the three combinations given 

in table 1 and used to generate the three curves displayed in figure 6. As 

seen in figure 6, each curve has two Vso values. 

Table 1. Parameter 

values for simulated 

curves. 

Figure 6. Curves 1 to 3 

from table 1 generated 

for Monte Carlo 

simulation. 

fhrve PT PI PS 

1 0.4 0.2 0.9 

2 0.4 0.2 0.7 

3 0.4 0.2 0.6 

0.8 

0.7 

g 0.6. 

n 

2 0.5. 

l! 0.4. 

0.3. 

0.2. 

0.1' / 

V.” 

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 

Resealed velocity (Oto 1) 
For the simulations, random binary responses were generated for 100 evenly 

spaced velocities between the values of 0 and 1. As can be seen in figure 6, the velocity values have been resealed between 0 and 1 so that the 

value of the bandwidth would be more meaningful. The cdf mixture model 

was simulated with SAS interactive matrix language (IML) for each of the 

curves in table 1. The process was repeated 500 times, and the average T/:,0 

values were computed over the 500 Monte Carlo repetitions. We view these 

average V&j values as assessing the ability of the llogr procedure to estimate 

the true I’,” values. We used the PRESS* procedure here to obtain the bandwidth for the liogr procedure. 

Tables 2 and 3 contain the results from the simulations for estimating the 

Vs, values. As can be seen in the fifth column of the tables, the local logistic 

regression does an exceptional job in estimating the lower and upper V~O 

values when they are calculable. The average absolute error of each estimate is less than 0.03. The quality of these estimates indicates the goodness 

of fit of the llogr procedure and its ability to estimate the I$) values. The 

excellent estimation can also be attributed to the ability of the PRESS* procedure in finding an appropriate value of the bandwidth to fit this family 

of curves. 

Note also that in the third column of tables 2 and 3, the number of vie 

values that were not calculable is in parentheses. In order for both V~OI, 

and VF,“~J to be estimated, the fitted llogr curve must cross the P(V) = 0.5 

line three times, as illustrated in figure 1. When the llogr curve crosses the 

P(V) = 0.5 curve line only once, then only one of the V~O values can be estimated: V&, if the crossing occurs at lower values of velocity, and V,NJ if the 

crossing occurs at upper values of velocity (with the other V~O value being 

incalculable). Data generated from a curve similar to curve 1 in figure 6 will 

more frequently fail to result in a both V 501, and &arr being estimable. Low 

penetration probabilities for high velocities in curve inhibit the movement 

of the sequential design to the highest velocities where a third crossing of 

P(w) = 0.5 could occur. Note that 31 percent of the time (156 samples out 

of 500), the V~O[J could not be calculated for curve 1: clearly an undesirable 

situation. One way to resolve this problem would be to increase the sample 

size at larger velocities. 

Table 2. Lower V5, 

estimation summary (incalculable values in 

Curve True L& Estimate L&L Average error Average abs. error 

1 0.20190 0.20827 (4 1) 0.00637 0.02590 

Parentheses). 2 0.20190 0.20782 (38) 0.00592 0.02425 

3 0.20190 0.20901 (45) 0.00711 0.02363 

Table 3. Upper I& 

estimation summary (incalculable values in 

parentheses). 

Curve True &,I, Estimate L%OL Average error Avgerage abs. error 

1 0.89999 0.88810 (156) -0.01190 0.02504 

2 0.69975 0.69845 (0) -0.00130 0.01834 

3 0.59810 0.58604 (22) -0.01205 0.02672 

12 
4. Discussion and Summary 

The research presented here has shown that the llogr procedure provides 

a more flexible and reasonable approach to fitting binary data from a nonmonotonic function than the classic logistic regression procedure, cdf mixture models, and some other nonparametric procedures (kernel and local 

linear regression). Through simulating a family of nonmonotone functions, 

it has been shown that the llogr procedure provides an outstanding estimate of V&I values with very small average and absolute errors, when 

calculable. 

13 
Bibliography 

Allen, D. A. (1974). The relationship between variable selection and data 

augmentation and a method of prediction. Technometrics 16, 125. 127. 

Chang, A. L., and Bodt, B. A. (1997). JTCG/AS Interlaboratory Ballistic Test 

Program-Final Report. U.S. Army Research Laboratory, ARL-TR-577. 

Chu, C.-K., and Marron, J. S. (1991). Choosing a kernel regression estimator. 

Stat. Sci. 6, 4, 404- 436. 

Cleveland, W. S. (1979). Robust locally weighted regression and smoothing 

scatterplots. J. Amer. Stat. Assoc. 74, 368, 829-836. 

Copas, J. B. (1983). Plotting p against X. Appl. Stat. 32, 1, 25 31. 

Dixon, W. J,. and Mood, H. M. (1948). A method for obtaining and analyzing sensitivity data. J. Amer. Stat. Assoc. 23, 109. 126. 

Einsporn, R. L., and Birch, J. B. (1993). Model-Robust Regression: Using Nonparametric Regression to Improve Parametric Regression Analyses. Virginia 

Polytechnic Institute and State University, Technical Report 93-5. 

Fan, J., Heckman, N. E., and Wand, M. P. (1995). Local polynomial kernel 

regression for generalized linear models and quasi-likelihood functions. 

J. Amer. Stat. Assoc. 90, 429, 141-150. 

Fan, J., and Gijbels, I. (1992). Variable bandwidth and local linear regression 

smoothers. Ann. Stat. 20, 2008.2036. 

Hardle, W. (1990). Applied Nonparametric Regression. Cambridge University 

Press, New York. 

Mays, James (1995). Model-Robust Regression Combining Parametric, Nonparametric, and Semiparametric Methods. Unpublished Ph.D. dissertation, Virginia Polytechnic Institute and State University, Blacksburg, Virginia. 

McCullagh, P., and Nelder, J. A. (1983). Generalized Linear Models. New York: 

Chapman and Hall. 

Myers, R. H. (1990). Classical and Modern Regression with Applications (2nd 

edition). PWS-Kent, Boston, Massachusetts. 

Nadaraya, E. A. (1964). On estimating regression. Theory Probab. Its Appf. 

USSR 9. 141- 142. 

Nottingham, Q. J., Birch, J. B., and Bodt, B. A. (2000). Local logistic regression: An application to Army penetration data, J. Stat. Comput. Simu. 66, 

35-50. 

14 
Robbins, H., and Monro, S. (1951). A stochastic approximation method. Annals Math. Stat. 22, 400-407. 

Ryan, T. P. (1997). Modern Regression Methods. New York: John Wiley and 

Sons. 

Wahba, G., and Wold, S. (1975). A completely automatic French curve: fitting spline functions by cross-validation. Commun. Stat. Ser. A 4, l- 17. 

15 
Appendix A. Numerical Considerations 

The values of the fits are obtained by using the IRLS algorithm. The procedure is outlined below. 

To find the proper bandwidth: 

For each fixed value of bandwidth b, perform the following: 

1. Compute the kernel weights for each value of velocity, Y = vi, for 

i= l,... :n,as 

2. Ateachv = zli:i = l,... , n, compute initial coefficients via local linear regression (llr) as 

&= ( xwix )-h%li: 

where the n x n diagonal matrix T/f/i has elements 1-“jj = h$. 

3. Compute n x 1 vector 311 as yrj = 

Yj -q&J 

f <&J and update the elements 

of IQ as 

~~i=F(ll:$O)(l-F(v:~io))h~, forj=l,...:n. 

4. Compute the updated estimates of the coefficients as 

5. If &, z ,&, quit and let & = Iii. Otherwise, replace &i by ,& and 

return to step 3 above. 

6. After convergence, compute the fit at 21 = vi as p (Q) = F (z&). 

7. After the fits are obtained for all n observations, compute the PRESS*(b) 

statistic as 

c (Yi - p (%)-i)’ 

PRESS* (b) = i=l 

n- tr [H”] ’ 

where p (z~i)-~ is the llogr (local logistic regression) fit based on the 

current bandwidth b, obtained with the ith data point (yi, vi) removed 

from the data set. The current version of the algorithm obtains ? (ui)-i 

by using steps 1 through 6 above, based on the n - 1 point data set. 

17 
8. Determine b*, the value of b that minimizes PRESS*(b). This can be accomplished either through a search routine (the algorithm currently 

uses a binary search to find h*) or by finding h* through a candidate 

list of bandwidths. 

Upon completion of steps 1 though 8, the value of b* has been found and 

the n fits are obtained at the 71 data points based on 2)*. Often, however, the 

fitted curve is desired over a second series of v values, which may or may 

not include the original data points. For example, the predicted probability 

of response may be desired for all velocities ranging from the minimum 

observed velocity to the maximum observed velocity at every 10 ft/s. In 

this case, steps 2 through 6 of the IRLS algorithm above are changed only by 

the replacement of Tli with ~(1, one of the velocity values where predictions 

are desired. 

18 
I Appendix B. Program Documentation & Comments 

Below is the SAS code, preceded by a few notes, to fit binary data via 

the nonparametric local logistic regression procedure (llogr) presented in 

this report. The llogr procedure uses PRESS* as the method for bandwidth 

selection. 

Some notes about the program are presented below. 

1. Lines 6, 7, and 10 read the data from a file. The file can be a text (.txt), 

print (.prn), or flat data file. The program is set up to read a file containing two columns of data: the first column contains the velocity, 

and the second column is the response (0 or 1). 

The INFILE statement reads the data set from its stored location, denoted by “fi lename.filetype.” 

2. Once the data have been read into the IML procedure, the regressor 

(in this case, the velocity) is resealed such that the values are between 

0 and 1. The program then finds the bandwidth to minimize PRESS*. 

Once the bandwidth has been found, the program does the following: 

l Obtains the llogr fit at each of the data points (lines 220- 248) 

l Finds the lower and upper V50 values (lines 252- 452) 

l Obtains predictions at lo-ft/s intervals (lines 461- 502) 

3. The program outputs the following (lines 518- 521): 

l The bandwidth (BW). 

l The chi-square statistic, the model degrees of freedom, and the 

mean squared error, denoted by CHISQR, MODELDF, and MSE, 

respectively. 

l The lower V&J (V&) and the upper V& (V&u) values. 

l The predictions at IO-ft/s intervals: Velocity (VEL), fitted value 

(YHATO), and lower and upper confidence limits for the fitted 

value (LCLO and UCLO). 

4. Once the output is obtained, the user can produce a plot of the fitted 

curve using Excel or another plotting application. However, the SAS 

code below also generates a plot with upper and lower confidence 

bounds in Computer Graphic Metafile (CGM) format. The code is on 

lines 538 to 552. The user must supply the filename and the directory. 

DO NOT CHANGE THE EXTENSION. Once the CGM file has been 

created, the user can view the plot in MS Word or another application 

that converts CGM files. 

19 
SAS Code 

The following represents the directory, filename, and filetype of the data 

set. The data set should be in text form or a flat data file. 

Options 1~330 ps=54 nodate nonumber, 

Title1 ‘Local Logistic Regression Analysis’, 

Title’2 ‘usinK Press* to Obtain BandwIdth’, 

DATA MACRO ; 

INFILE “d:\research\AR~IY\HACROCOllPtlSITE ORIGINAL DATA PRII”; 

Here variables are being read from the data set. In this case, they are velocity and response. 

INPUT VELOCITY RESPONSE; 

DATA LLOGR; 

SET HACRO, 

RUN / 

PROC SORT; 

BY VELOCITY; 

RUN; 

The following sorts the data set by velocity 

PROC IHL; 

Use LLOGR, 

Read All Into X Var {VELOCITY}; 

Read All Into PI Var {RESPONSE}; 

Now the velocity is being converted/transformed to the 0 to 1 data range. 

This is so that we can obtain a better interpretation of the bandwidth. That 

is, given a finite range for the bandwidth, one can interpret its relative size. 

However, if the data are not converted, the range of the bandwidth is from 

0 to infinity, allowing a less meaningful interpretation. 

N=NROW(X); 

MINX=HIN(X); 

HAXX=flAX(X); 

LENCTH=(MAXX-!INX)/iO + 1, 

VEL=J(LENGTH,l,O), 

v=o; 

DO I=HINX TO HAXX BY 10; 

V=V+l; 

VEL[V,]=(I-HIflX)/(HAXX-)IIAX); 

END; 

20 
X=(X-MINX)/(MAXX-MINX); * Transforms X-variable to [O,l] Range; 

*** INITIALIZATIOW ***; 

H=J(N,i,i); 

YHAT=J(N,l,l); 

LCL=J(N,l,i); 

UCL=J(N,i,i); 

PRSS=J(3,1,0); 

HO=J(N,l,O); 

TYHATMI=J(N,l,i); 

*** END OF INITIALIZATIONS ***; 

Here begins the subroutine to obtain the true PRESS statistic and the subsequent bandwidth found to minimize this PRESS statistic. 

STARTRUEYMI; 

DO I=1 TO N; 

FREE XMI; 

IF I=1 THEN 

DO; 

XE!I=X[Z:N,l; 

PHI=PI[2:N,l; 

END; 

ELSE IF I=N THEN 

DO; 

XMI=X[l:N-i,]; 

PMI=PI[i:N-l,]; 

END; 

ELSE 

DO; 

XMIi=X[l:I-i,]; 

PMIi=PI[i:I-l,]; 

XMIZ=X[I+i:N,]; 

PMIZ=PI[I+l:N,]; 

XMI=XHIi//XMI2; 

PMI=PMIl//PMIZ; 

END; 

XMI=J(N-i,i,i)IlXMI; 

xo=lllx[I,1; 

DO=X[I,]-XMI[,2]; 

DSqO=DO##2; 

DINT0 = (-i)#DSqO/(BW##2); 

EXPDO = EXP(DINT0); 
ESUHO = SUH(EXPDO), * RON SUH OF NUtlERATOR VALUES, 

KO = EXPDO/ESUHO, 

BHAT=I~IV(X~II'~DIAG(KO~~XI1I~*X~II'*DIAG(KO)*PI1I, 

DO ITER=l TO 2; 

Y=XHI*BHAT, 

P=i/(l+EXP(-Y)); 

9=1-P; 

Z=EXP(-Y)/((l+EXP(-Y))##Z), 

W=DIAG(P#Q)#DIAG(KO); 

XPXI=INV(XMI'*N*XMI), 

BHAT=BHAT+XPXI*XMI'*N*((PHI-P)/Z); 

END, 

TYHATMI[I,l]=l/(l+EXP(-XO*BHAT)); 

END, 

FINISH; 

* Subroutine to begin the search for the 

bandwidth that minimizes Press*; 

START BUSEARCH; 

BWO = IO.1, 0.5, 0.7}, 

DO A = 1 TO 3; 

BW=BVO[A,]; 

D=X*J(l,N,i)-J(N,l,l)*X'; 

D=D'; 

DSq=D##Z; 

DINT = (-1)1(DSq/(BWlt2)lJ(N,N,l)); 

EXPD = EXP(DINT); 

ESUH = EXPD[,+]; * ROW SUH OF NUMERATOR VALUES; 

K q EXPD/(ESURjJ(l,N,l)); 

RUN TRUEYNI; 

PRSS[A,] = SUM((PI-TYHATMI)N#2)/(N-TRACE(K)); 

END, 

* SECOND LOOP SEARCHES FOR THE BANDMIDTH TO MINIMIZE The True PRESS*; 

ITERB=O; 

MINI=HIN(PRSS); 

MINOLD=HAX(PRSS); 

DO WHILE (ABS(nIf~I-MIAOLD)>lE-8); 

ITERB=ITERB+i; 

MINOLD=NAX(PRSS); 

IF MINI = PRSS[l,] THEN 

DO; 

AHIN = 1; 

SET1 = BWO[l,]; 

22 
SET2 =2#BWO[l,]-BWO[2,]; 

SET3 =BWO[1,1/2; 

VEC = SET3//SET2; 

BVO[l,] q MAX(VEC); 

BM0[3,1 = BWO[2,1; 

BWO[2,1 = SETl; 

END; 

IF MINI = PRSS[2,] THEN 

DO; 

AMIN = 2; 

MAXI = MAX(PRSS); 

IF MAXI = PRSSt1,1 THEN 

DO; 

AMAX = 1; 

BWO[l,] = BUO[2,1; 

BWO[2,1 = (BWO[2,1 + BW0[3,1)/2; 

END; 

IF MAXI = PRSS[3,1 THEN 

DO; 

AHAX = 3; 

BWO[3,] q BWO[2,1; 

BW0[2,1 = (BWO[2,1 + BWO[1,1)/2; 

END; 

END; 

IF MINI = PRSS[3,1 THEN 

DO; 

AMIN = 3; 

SET1 =BWO[3,3; 

BWO[3,] = 2tBW0[3,]-BW0[2,]; 

BWO[l,] = BWO[2,]; 

BWO[Z,]=SETi; 

END; 

* OBTAIN THE LOCAL LOGISTIC FIT FOR THE NEW VALUE OF THE BANDWIDTH*: 

BW = BWO[AMIN,]; 

D=X*J(l,N,i)-J(N,l,l)*X'; 

D=D'; 

DSQ=D##2; 

DINT = (-l)#(DSQ/(BW#rt2)#J(N,N,l)); 

EXPD = EXP(DINT); 

23 
ESUH = EXPD[ ,+] , * RO!: SUll OF NUIIERATOR VALUES 

K = EXPD/(ESU~l~J(l,~l,l)), 

RUN TRUEYIII, 

IIII1PRESS = SUa((PI-TYHATnI)##2)1(A-TRACE(K)); 

* VALUE OF PRESS FOR THE CURREliT BAIIDI!IDTH FOR THIS ITERATION* 

1: RESET POSITIONS OF BAIIDKIDTH AND PRSS VECTORS DEPENDIIIC ON THE VALUE OF 

PRESS FOR THE CURRENT VALUE OF THE BANDWIDTH*; 

IF AMIN = 1 THEN 

DO, 

PRSS[3,] = PRSS[Z,]; 

PRSS[2,] q PRSS[l,]; 

PRSS[ 1 ,] = MINPRESS; 

END, 

IF AMIN = 3 THEN 

DO; 

PRSS[l,] = PRSS[Z,], 

PRSS[Z,] = PRSS[3,], 

PRSS[3,] = MINPRESS, 

END; 

IF AMIN = 2 THEN 

DO, 

IF AMAX =1 THEN PRSS[l,] = PRSS[2,], 

IF AMAX =3 THEN PRSS[3,] = PRSS[Z,]; 

PRSS[2,] = HINPRESS; 

END; 

HINI=MIN(PRSS), 

END; * SEARCH FOR OPTIMAL BANDk’IDTH ENDS, 

FINISH; 

WGT=i, 

RUN BblSEARCH; 

* Haung obtained the banduidth, ue can now obtain the local logistic regression fit*; 

X2=J(N,l,l)IIX; 

H=J(N,l,i); 

Do I=1 to N, 

D=X[I,]-X; 

DS$=D##2, 

Dlnt=(-i)RDSq/(BW##2); 

Expd=Exp(Dint), 

Esum=SUM(Expd); 

24 
K=Expd/Esum; 

H[I,l = K[I,I; 

BHAT=INV(X2'*DIAG(K)*X2)*X2'IDIAG(K)*PI; 

DO ITER=i TO 2; 

Y=X2*BHAT; 

P=i/(i+EXP(-Y)); 

Q=l-P; 

Z=P#Q; 

W=DIAC(P%Q)tDIAC(K); 

XPXI=INV(X2'*W*X2); 

BHAT=bhat+XPXI*X2'*W*((PI-P)/Z); 

END; 

Y=XZ*BHAT, 

P=i/(i+EXP(-Y)); 

Q=i-P; 

Z=EXP(-Y)/((l+EXP(-Y))##2); 

W=DIAG(PtQ)#DIAG(K); 

XPXI=INV(X2'*W*X2); 

YHAT[I,]=l/(l+Exp(-X2[I,]*BHAT)); 

END; 

* ALGORITHM FOR COHPUTING THE LOWER V50 VALUE*; 

COUNT=O; 

X5O=J(3,1,0); 

DO I=1 TO N; 

IF COUNT=0HEN 

DO; 

IF YHAT[I,])0.5 THEN 

DO; 

X50[1,]=X[I,]-0.05; 

X50[2,]=X[I,]-0.025; 

x5ot3,1=x[I,1; 

COUNT=l; 

END; 

END; 

END; 

SET50 = J(3,1,0); 

t OBTAIN PHAT FOR 3 VALUES OF X (DOSE) *; 

DO I=1 TO 3; 

D=XSO[I,]-X; 

DSQ=Dt#2; 

DINT=(-i)#DSQ/(BW##2); 

25 
EXPD=EXP(DI!IT); 

ESUtI=SUtI(EXPD); 

K=EXPD/ESWI, 

BHAT=INV(X2'*DIAG(K)*X2)*X2'*DIAG(K)*PI; 

DO ITER=l TO 2, 

Y=XZ*BHAT; 

P=l/(l+EXP(-Y)), 

9=1-P, 

Z=EXP(-Y)/((l+EXP(-Y))tt2), 

W=DIAG(P#q)#DIAG(K); 

XPXI=INV(X2'*k'*X2); 

BHAT=BHAT+XPXI*X2'*W*((PI-P)/Z), 

END; 

X0=11 IX50[1,]; 

~ET~~[I,I=~/(~+EXP(-XO~BHAT)); 

END, 

PTEBP=O; 

B=O; 

DO B=i TO 200;*WHILE (ABS(PTEMP-0.5))1E-4); 

B=B+i; 

IF SET50[2,] ) 0.5 THEN 

DO; 

INTERVAL=i; 

SET1 =X50[1,1; 

SET2 =(X50[2,l+X50[1,~)/2; 

SET3 q X50[2,]; 

X50[1,3 = SETl, 

X50[2,] = SET2; 

X50[3,1 = SET3, 

EAD; 

IF SET50[2,] ( 0.5 THEN 

DO; 

INTERVAL=2; 

SET1 =X50[2,]; 

SET2 =(X5013,] + X50[2,])/2; 

SET3 =X50[3,1; 

X50[1,]=SETl, 

X50[2,1=SET2, 

X50[3,]=SET3, 

END; 

26 
* OBTAIN YHAT FOR NEN VELOCITY*; 

V5OL =X50[2,]; 

D=VSOL-X; 

DSU=D##2; 

Dint=(-i)#DSq/(BW##2); 

Expd=Exp(Dint); 

Esum=SUM(Expd); 

K=Expd/Esum; 

BHAT=INV(X2'*DIAG(K)*X2)*X2'IDIAG(K)*PI; 

DO ITER=i TO 2; 

Y=X2*BHAT; 

P=l/(l+EXP(-Y)); 

q=1-P; 

Z=EXP(-Y)/((l+EXP(-Y))?#2); 

W=DIAG(P#$)#DIAG(K); 

XPXI=INV(X2'*W*X2); 

BHAT=BHAT+XPXI*X2'*W*((PI-P)/Z); 

END; 

xo=1Ilv5oL; 

PTEMP=i/(l+Exp(-XO*BHAT)); 

IF INTERVAL=1 THEN SET50=SET50[1,1//PTEMP//SET50[2,]; 

IF INTERVAL=2 THEN SET50=SETSO[Z,]//PTEMP//SET50[3,]; 

END; * SEARCH FOR V5OL ENDS; 

IF (ABS(PTEMP-O.~))~E-~) THEN v5oL=.; 

* ALGORITHM FORCOMPUTING THE UPPER V50 VALUE *; 

COUNT=O; 

DO I=N TO 1 BY (-1); 

IF COUNT=0HEN 

DO; 

IF YHAT[I,](O.S THEN 

DO; 

XSO[i,]=X[I,]; 

X50[2,]=X[I,]+O.O25; 

X50[3,]=X[I,]+O.O5; 

COUNT=l; 

END; 

END; 

END; 

SET50 = 3(3,1,0); 

27 
* OBTAIli PHAT FOR 3VALUES OF X (DOSE) * 

DO I=1 TO 3; 

D=XSO[I,]-X, 

DSQ=D##Z; 

DIrfT=(-l)#DSQ/(Bu##2); 

EXPD=EXP(DINT); 

ESWSU!l(EXPD), 

K=EXPD/ESUtI; 

BHAT=IIIV(X2'*DIAG(K)rXZ)tX2'IDIAG(K)*PI, 

DO ITER=l TO 2; 

Y=XZ*BHAT, 

P=i/(l+EXP(-Y)); 

4=1-P; 

Z=EXP(-Y)/((l+EXP(-Y))f#2); 

h'=DIAG(P#Q)tDIAG(K); 

XPXI=INV(X2'*k'*X2); 

BHAT=BHAT+XPXI*XZ'*W*((PI-P)/Z); 

END; 

X0=11 IX50[1,]; 

SET5O[I,]=i/(l+EXP(-XO*BHAT)); 

END; 

PTEHP=O; 

B=O; 

DO B=l TO 200; *WHILE (ABS(PTEHP-O.S)>lE-4); 

B=B+l; 

IF SET50[2,] ) 0.5 THE! 

DO; 

INTERVAL=l; 

SET1 =X50[1,1; 

SET2 =(X50[2,]+X50[1,])/2; 

SET3 =X50[2,]; 

X50[1,] = SETl, 

X50[2,] = SETZ, 

X50[3,1 = SET3, 

END, 

IF SET50[2,] < 0.5 THEN 

DO; 

IIITERVAL=2; 

SET1 q X50[2,]; 

SET2 =(X50[3,] + X50[2,])/2, 

28 
SET3 = X50[3,1; 

XSO[i,]=SETi; 

XSO[Z,]=SET2; 

X50[3,]=SET3; 

END: 

x OBTAIN YHAT FOR NEW VELOCITY*; 

V5OU = X50[2,]; 

D=VSOU-X; 

DSQ=Dt#2; 

Dint=(-i)#DSQ/(BW##2); 

Expd=Exp(Dint); 

Esum=SUM(Expd); 

K=Expd/Esum; 

BHAT=INV(XZf*DIAG(K)*X2)*X2'tDIAG(K)*PI; 

DO ITER=l TO 2; 

Y=X2*BHAT; 

P=i/(i+EXP(-Y)); 

Q=i-P; 

Z=EXP(-Y)/((i+EXP(-Y))ltt2); 

W=DIAG(P#Q)#DIAG(K); 

XPXI=INV(X2'*W*X2); 

BHAT=BHAT+XPXI*X2'*W*((PI-P)/Z); 

END; 

xo=1Ilv5ou; 

Y=X2*BHAT; 

P=l/(i+EXP(-Y)); 

q=i-P; 

Z=EXP(-Y)/((i+EXP(-Y))t#2); 

W=DIAG(P#Q)tDIAG(K); 

XPXI=INV(X2'*W*X2); 

PTEMP=l/(i+Exp(-XO*BHAT)); 

IF INTERVAL=1 THEN SET50=SET50[1,]//PTEMP//SET50[2,]; 

IF INTERVAL=2 THEN SET5O=SET50[2,]//PTEHP//SET50[3,]; 

END; *SEARCH FOR V5OU ENDS; 

IF (ABS(PTEM~-0.5))1~-2) THENv5ou=.; 

MDDELDF=SUM(H); 

DFE=N-MODELDF; 

CHISQR=SUM(((PI-YHAT)##2)/(YHAT#(i-YHAT))); 

MSE=CHISqR/DFE; 
*** BAKE PREDICTIUllS AT INCREtIE1lTS OF 10FRDll Ml VELOCITY TO IIAX VELOCITY ***, 

YHATO=J(tEIIGTH.1,0), 

LCLO=J(tEIIGTH,i,O); 

UCtO=J(tEHGTH,i,O), 

DO I=1 TO LEIGTH, 

D=VEL[I,]-X, 

D=D'; 

DSQ=D##2; 

DINT = (-i)#DSQ/(BV##2), 

EXPD = EXP(DIIIT); 

ESUB = SUH(EXPD); * ROli SUtl OF NUHERATOR VALUES, 

K = EXPD/ESUll, 

BHAT=INV(X2'*DIAG(K)*X2)*X2'rDIAG(ioxPI; 

DO iter = 1 to 2, 

Y=X2*BHAT; 

P=l/(l+EXP(-Y)); 

Q=l-P; 

Z=P#Q, 

W=DIAG(PtQ)#DIAG(K); 

XPXI=INV(X2'*W*X2), 

BHAT=BHAT+XPXI*X2'*W*((PI-PI/Z); 

END; 

Y=X2*BHAT; 

P=l/(l+EXP(-Y)); 

q=1-P; 

Z=EXP(-Y)/((i+EXP(-Y))tR2); 

Y=DIAG(PtLq)#DIAG(K), 

XPXI=INV(XZ'*!+'*X2); 

XO=llIVEL[I,]; 

YO=XO*BHAT; 

YHATO[I,]=l/(l+EXP(-YO)); 

PO=YHATO[I,]; 

VO=DIAC(K##2)lDIAG(P#Q); 

VBHATO=XPXI*X2'*VO*X2*XPXI; 

ZO=EXP(-YO)/((l+EXP(-YO))M?); 

VYHATO=(ZO#~2)#XO*VBHATO*X01; 

LCLO[I,]=YHATO[I,]-1.96*SQRT(VYHATO), 

UCLO[I,]=YHATO[I,]+l.96*SQRT(VYHATO), 

IF LCtO[I,]<O THEN LCtO[I,]=O; 

IF UCtO[I,]>l THEN UCLO[I,]=l; 

END; 

30 
VEL=VEL*(MAXX-MINX)+MINX; 

CZ={"XO" "YHATO" "LCLO" "UCLO"); 

PRED~=VEL~~YHAT~~~LCLO~~UCLO; 

CREATE YHAT2 FROH PRED:![COLNAHE=C2]; 

APPEND FROH PREOZ; 

X=X*(MAXX-MINX)+HINX; 

V50L=V50L*(MAXX-MINX)+MINX; 

V5OU=V50U*(MAXX-MINX)+MINX; 

Outputs are the bandwidth (BW), Chi-square statistic, model degrees of 

freedom, mean-squared error, I&O values, and the data (VEL, YHATO, LCLO, 

and UCLO) at lo-ft/s intervals. 

Print BW; 

Print CHISQR MODELDF MSE /; 

Print V5OL V5OU I; 

PRINT VEL YHATO LCLO UCLO; 

The following creates a SAS data set that can be then used in the Proc Glot 

procedure. 

Cl={"X" "PI"}; 

PREDl=XIIPI; 

CREATE YHAT~ FROM PRED~[COINAME=C~]; 

APPEND FROM PREDi; 

DATA PHAT3; 

MERGE YHATl YHAT2; 

PROC SORT; 

BY X; 

Below is the code to create a CGM format graph. This format gives us the 

capability to insert the plot in a Word file or any other application that reads 

the CGM file type. 

FILENAME GSASFILE 'filename.cgm'; 

GOPTIONS RESET=GOPTIONS DEVICE=CGM GACCESS=GSASFILE 

GPROTOCOL=SASGPASC GSFLEN=80; 

SYMBOL1 VALUE=DOT; 

SYMBOL2 INTERPOL=JOIN LINE=1 VALUE=NONE; 

SYMBOL3 INTERPOL=JOIN VALUE=NONE LINE=4; 

SYMBOL4 INTERPOL=JOIN VALUE=NONE LINE=4; 

AXIS1 LABEL=("VELOCITY, fps"); 

AXIS2 ORDER=(O TO 1 BY 0.1) 

LABEL=(ROTATE=gO ANGLE=-90 "RESPONSE"); 

PROC GPLOT DATA=PHAT3; 

31 
PLOT PI*X YHATO*XO UCLOtXO LCLO*XO / OVERLAY HAXIS=AXISl VAXIS=AXISZ, 

TITLE1 J=C H=2 ‘Local Logistic Regression Analym’; 

TITLE2 J=C H=2 ‘All llacrocomposlte Data’; 

Rllll, 

Quit; 

32 
Admnstr 

Defns Tech1 Info Ctr 

ATTN DTIC-OCP 

8725 John J Kingman Rd Ste 0944 

FT Belvoir VA 22060-6218 

DARPA 

ATTN S Welby 

3701 N Fairfax Dr 

Arlington VA 22203-1714 

Oft of the Secy of Defns 

ATTN ODDRE (R&AT) 

The Pentagon 

Washington DC 20301-3080 

Oft of the Secy of Defns 

ATTN OUSD(A&T)/ODDR&E(R) R J Trew 

3080 Defense Pentagon 

Washington DC 20301-7100 

AMCOM MRDEC 

ATTN AMSMI-RD WC McCorkle 

Redstone Arsenal AL 35898-5240 

US Army TRADOC 

Battle Lab Integration & Tech1 Dirctrt 

ATTN ATCD-B 

FT Monroe VA 23651-5850 

Hdqtrs Dept of the Army 

ATTN DAMO-FDT 

400 Army Pentagon Rm 3C5 14 

Washington DC 20301-0460 

US Military Acdmy 

Mathematical Sci Ctr of Excellence 

ATTN MADN-MATH MAJ M Huber 

Thayer Hall 

West Point NY 10996-1786 

Dir for MANPRINT 

Oft of the Deputy Chief of Staff for Prsnnl 

ATTN J Hiller 

The Pentagon Rm 2C733 

Washington DC 20301-0300 

SMC/CZA 

2435 Vela Way Ste 1613 

El Segundo CA 902455500 

r 

Distribution 

US Army ARDEC 

ATTN AMSTA-AR-TD 

Bldg 1 

Picatinny Arsenal NJ 07806-5000 

US Army Info Sys Engrg Cmnd 

ATTN AMSEL-IE-TD F Jenia 

FT Huachuca AZ 85613-5300 

US Army Natick RDEC Acting Tech1 Dir 

ATTN SBCN-T P Brandler 

Natick MA 01760-5002 

US Army Simulation Train & Instrmntn 

Cmnd 

ATTN AMSTI-CG M Macedonia 

ATTN J Stahl 

12350 Research Parkway 

Orlando FL 32826-3726 

US Army Tank-Automtv Cmnd RDEC 

ATTN AMSTA-TR J Chapin 

Warren MI 48397-5000 

Nav Surfc Warfare Ctr 

ATTN Code B07 J Pennella 

17320 Dahlgren Rd Bldg 1470 Rm 1101 

Dahlgren VA 22448-5100 

USAFAS 

ATTN ATSF-CBL G Durham 

FT Sill OK 73503-5600 

Marine Corps Artillery Detachment 

ATTN COL L Stuart 

759 McNair Ave 

Lawton OK 73503 

Marine Corps Detachment 

ATTN CPT Thomas 

18B Sparrow Rd 

Fallbrook CA 92028 

Univ of Texas at Austin Inst for Advncd 

Tchnlgy 

PO Box 202797 

Austin TX 78720-2797 

33 
Distribution (cont’d) 

British Aerospace Defence Sys Limited 

ATTN P Wilkins 

Vicarage Lane 

Ilford Essex IGl 4Q 

United Kingdom 

Hicks & Assoc Inc 

ATTN G Singley III 

1710 Goodrich Dr Ste 1300 

McLean VA 22102 

Roke Manor Rsrch Roke Manor 

ATTN I Sandoe 

Romsey Hampshire SO5 1 OZN 

United Kingdom 

US Army Rsrch Lab 

ATTN AMSRL-CI-CD A Neiderer 

ATTN AMSRL-CI-CD B Bodt (5 copies) 

ATTN AMSRL-CI-CD C Hansen 

ATTN AMSRL-CI-CD E Heilman 

ATTN AMSRL-CI-CD G Moss 

ATTN AMSRL-CI-CD J Forester 

ATTN AMSRL-CI-CD J OMay 

ATTN AMSRL-CI-CD M Thomas 

ATTN AMSRL-CI-CD P Jones 

US Army Rsrch Lab (cont’d) 

ATTN AMSRL-CI-CD R C Kaste 

ATTN AMSRL-CI-LP (305) (4 copies) 

Aberdeen Proving Ground MD 2 1005 

US Army Rsrch Lab 

ATTN AMSRL-CS-IO-F1 M Carr 

Aberdeen Proving Ground MD 2 1005-5066 

Director 

US Army Rsrch Oft 

ATTN AMSRL-RO-D JCI Chang 

PO Box 12211 

Research Triangle Park NC 27709 

US Army Rsrch Lab 

ATTN AMSRL-CI-AI-R Mail & Records 

Mgmt 

ATTN AMSRL-CI-AP Tech1 Pub (2 copies) 

ATTN AMSRL-CI-LL Tech1 Lib (2 copies) 

ATTN AMSRL-D D R Smith 

ATTN AMSRL-DD J M Miller 

ATTN AMSRL-SE-SA J Eicke 

ATTN AMSRL-SE-SA R Frankel 

ATTN AMSRL-SE-SA S Tenney 

Adelphi Md 20783- 1197 

34 
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, mcluding the tome for reviewing instructions. searching existing data SOurCOs. 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect Of this 

collection of Information. Including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson 

Davis Highway, Suite 1204, Arlington, VA 22202.4302. and to the Offlce of Management and Budget. Papenvoti Reduction Project (0704-0188). Washington, DC 20503. 

AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED 

January 200 1 

. TITLE AND SUBTITLE Modeling Nonmonotonic Dose-Response Curves 5. FUNDING NUMBERS 

DA PR: DYlO 

PE: 622783YlO 

;.AUTHOR(S) Quinton J. Nottingham, Jeffrey B. Birch (Virgina 

Polytechnic Institute & State University), Barry A. Bodt (ARL) 

‘. PERFORMING ORGANIZATION NAME(S) AND ADDRESS 6. PERFORMING ORGANIZATION 

U.S. Army Research Laboratory REPORT NUMBER 

Attn: AMSRL- CI-CD email: babodt@arl.mil ARL-TR-1972 

2800 Powder Mill Road 

Adelphi, MD 20783-l 197 

I. SPONSORlNWMONlTORlNG AGENCY NAME(S) AND ADDRESS 

U.S. Army Research Laboratory 

10. SPONSORING/MONITORING 

AGENCY REPORT NUMBER 

2800 Powder Mill Road 

Adelphi, MD 20783-l 197 

Il. SUPPLEMENTARY NOTES 

ARL PR: 1TEDKC 

AMS code: 622783YlOll 

12a. DISTRIBUTION/AVAILABILITY STATEMENT Approved for public release; 12b. DISTRIBUTION CODE 

distribution unlimited. 

13. ABSTRACT (Maximum 200 words) 

A number of procedures have been used to analyze nonmonotonic binary data to predict the 

probability of response. Some classical procedures are the Up and Down strategy, the 

Robbins-Monro procedure, and other sequential optimization designs. Recently, nonparametric 

procedures such as kernel regression and local linear regression have been applied to this type of 

data. It is well known that kernel regression has problems fitting the data near the boundaries, 

and a drawback with local linear regression is that it may be “too linear” when fitting data from 

a curvilinear function. This report introduces a procedure called local logistic regression, which 

fits a logistic regression function at each of the data points. United States Army projectile data are 

used in an example that supports the use of local logistic regression for analyzing nonmonotonic 

binary data for certain response curves. Properties of local logistic regression are presented along 

with simulation results that indicate some of the strengths of the procedure. 

MSUBJECTTERMS Logistic regression local smoothing, nonparametric, 15. NUMBER OF PAGES 

nonmonotonic 40 

16. PRICE CODE 

17. SECURITY CLASSIFICATION 18. SECURITY CLASSIFICATION 

OF REPORT 

19. SECURITY CLASSIFICATION 

OF THIS PAGE 

20. LIMITATION OF ABSTRACT 

Unclassified 

OFABSTRACT 

Unclassified Unclassified SAR 

.^.. --._ _. _______ ^. .- ^^^ ,- ^ ^^.