{VERSION 2 3 "IBM INTEL NT" "2.3" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0
1 0 0 0 0 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0
0 }{CSTYLE "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }
{CSTYLE "2D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "
" -1 256 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 257 "" 0 1 0
0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 1 0 0 0 0 0
0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE ""
-1 260 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 261 "" 0 1 0 0
0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 262 "" 0 1 0 0 0 0 1 0 0 0 0 0 0
0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE ""
-1 264 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 265 "" 0 1 0 0
0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 266 "" 0 1 0 0 0 0 1 0 0 0 0 0 0
0 0 }{CSTYLE "" -1 267 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE ""
-1 268 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 269 "" 0 1 0 0
0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 270 "" 0 1 0 0 0 0 1 0 0 0 0 0 0
0 0 }{CSTYLE "" -1 271 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE ""
-1 272 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 273 "" 1 12 0
0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 274 "" 1 14 0 0 0 0 0 0 0 0 0
0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0
0 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 1" 0
3 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 }1 0 0 0 6 6 0
0 0 0 0 0 -1 0 }{PSTYLE "Heading 2" 3 4 1 {CSTYLE "" -1 -1 "" 1 14 0
0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 4 4 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple
Output" 0 11 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3
3 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 11 12 1 {CSTYLE "" -1 -1 ""
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }
{PSTYLE "Maple Plot" 0 13 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0
0 0 0 0 }3 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Fixed Width" 0 17
1 {CSTYLE "" -1 -1 "Courier" 1 10 0 0 0 0 0 0 0 0 0 0 3 0 0 }3 0 0 -1
-1 -1 0 0 0 0 0 0 17 0 }{PSTYLE "Title" 0 18 1 {CSTYLE "" -1 -1 "" 1
18 0 0 0 0 0 1 1 0 0 0 0 0 0 }3 0 0 -1 12 12 0 0 0 0 0 0 19 0 }
{PSTYLE "Author" 0 19 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 }3 0 0 -1 8 8 0 0 0 0 0 0 -1 0 }{PSTYLE "" 0 256 1 {CSTYLE "" -1
-1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }
{PSTYLE "" 19 257 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0
}0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 0 258 1 {CSTYLE "" -1 -1
"" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }
{PSTYLE "" 17 259 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
}2 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }}
{SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT 274 59 "File: depart\\math\\maple
\\misc\\logdown.mws Date: 2-mar-1998" }}}{PARA 259 "" 0 "" {TEXT
273 16 "22 February 1998" }}{PARA 18 "" 0 "" {TEXT -1 72 "Solving Logi
stic Regression Problems Using Downhill Simplex Minimization" }}{PARA
257 "" 0 "" {TEXT -1 119 "Leslie Wright, MD\nDepartments of Psychiatry
and Education\nQueen's University\nKingston, Ontario, Canada\nlcwrigh
t@kos.net" }}{PARA 0 "" 0 "" {TEXT -1 299 "The purpose of this sheet i
s to demonstrate an application of Francis Wright's downhill simplex m
inimization procedure to the problem of logistic regression with a dic
hotomous dependent variable. The problem basically boils down to fitt
ing, by maximum likelihood estimation, a function of the form " }}
{PARA 256 "" 0 "" {XPPEDIT 18 0 "P=1/(1+exp(-(a+b*x)))" "/%\"PG*&\"\"
\"\"\"\",&\"\"\"F&-%$expG6#,$,&%\"aGF&*&%\"bGF&%\"xGF&F&!\"\"F&F2" }}
{PARA 0 "" 0 "" {TEXT -1 6 "where " }{XPPEDIT 18 0 "x" "I\"xG6\"" }
{TEXT -1 61 " is an independent variable (continuous or dichotomous), \+
and " }{XPPEDIT 18 0 "P" "I\"PG6\"" }{TEXT -1 107 " is the probability
of certain event occuring. In the remainder of this sheet, this even
t is coded for by " }{XPPEDIT 18 0 "y" "I\"yG6\"" }{TEXT -1 88 ", whic
h takes on the value 1 in the case of an occurrence and 0 if not. The
expression " }{XPPEDIT 18 0 "a+b*x" ",&%\"aG\"\"\"*&%\"bGF$%\"xGF$F$
" }{TEXT -1 39 " found in the exponent is known as the " }{TEXT 259 5
"logit" }{TEXT -1 48 ", and we can readily rearrange the above to the \+
" }{TEXT 260 21 "logit transformation " }{TEXT -1 22 "of the logistic \+
model:" }}{PARA 258 "" 0 "" {XPPEDIT 18 0 "ln(P/(1-P))=a+b*x" "/-%#lnG
6#*&%\"PG\"\"\",&\"\"\"F(F'!\"\"F+,&%\"aGF(*&%\"bGF(%\"xGF(F(" }}
{PARA 0 "" 0 "" {TEXT -1 98 "The left-hand side of this equation is kn
own as the log-odds ratio. Since a mere substitution of " }{XPPEDIT
18 0 "y" "I\"yG6\"" }{TEXT -1 79 " for P for each data point leads to \+
the log-odds ratio being undefined--either " }{XPPEDIT 18 0 "ln(0/1)"
"-%#lnG6#*&\"\"!\"\"\"\"\"\"!\"\"" }{TEXT -1 4 " or " }{XPPEDIT 18 0 "
ln(1/0" "-%#lnG6#*&\"\"\"\"\"\"\"\"!!\"\"" }{TEXT -1 120 "--it is evid
ent that a least-squares approach to this linearization is not feasibl
e. The answer to this problem is the " }{TEXT 261 39 "method of maxim
um-likelihood estimation" }{TEXT -1 225 ". This method is derived and
applied in this worksheet by deriving the likelihood function from th
e Bernoulli distribution in question and taking the logarithm to obtai
n Log Likelihood function. We then find the parameters " }{XPPEDIT
18 0 "a" "I\"aG6\"" }{TEXT -1 5 " and " }{XPPEDIT 18 0 "b" "I\"bG6\""
}{TEXT -1 67 " that maximize this Log-Likelihood function given the da
ta at hand." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 8 "restart;" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 63 "Th
eoretical Background and Derivation of Log-Likelihood Formula" }}
{PARA 0 "" 0 "" {TEXT -1 91 "The following derivation is based on the \+
lucid treatment of Kleinbaum, Kupper, and Muller, " }{TEXT 270 59 "App
lied Regression Analysis and Other Multivariable Methods" }{TEXT -1
41 ", PWS-Kent Publishing, 1988, pp. 512-514." }}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 96 "First we start with a derivation using the Bernoulli, o
r point binomial, distribution. For the " }{XPPEDIT 18 0 "i" "I\"iG6
\"" }{TEXT -1 35 "-th observation in a collection of " }{XPPEDIT 18 0
"n" "I\"nG6\"" }{TEXT -1 20 " data points, where " }{XPPEDIT 18 0 "the
ta[i]" "&%&thetaG6#%\"iG" }{TEXT -1 50 " represents the probability of
a \"success\" (i.e., " }{XPPEDIT 18 0 "y[i]=1" "/&%\"yG6#%\"iG\"\"\"
" }{TEXT -1 6 ") and " }{XPPEDIT 18 0 "1-theta[i] " ",&\"\"\"\"\"\"&%&
thetaG6#%\"iG!\"\"" }{TEXT -1 40 " is the probabiltiy of \"failure\" (
i.e., " }{XPPEDIT 18 0 "y[i]=0" "/&%\"yG6#%\"iG\"\"!" }{TEXT -1 46 "),
the associated probability, or likelihood, " }{XPPEDIT 18 0 "l[i]" "&
%\"lG6#%\"iG" }{TEXT -1 50 ", of the particular outcome may be express
ed thus:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "l[i]:=theta[i]^y[i]*(1-
theta[i])^(1-y[i]);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>&%\"lG6#%\"iG*
&)&%&thetaGF&&%\"yGF&\"\"\"),&F.F.F*!\"\",&F.F.F,F1F." }}}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 24 "The logisitic function, " }{XPPEDIT 18 0
"1/(1+exp(-(a+b*x[i])))" "*&\"\"\"\"\"\",&\"\"\"F$-%$expG6#,$,&%\"aGF$
*&%\"bGF$&%\"xG6#%\"iGF$F$!\"\"F$F3" }{TEXT -1 47 ", is especially wel
l-suited as a candidate for " }{XPPEDIT 18 0 "theta[i]" "&%&thetaG6#%
\"iG" }{TEXT -1 94 ", since we want an expression whose values lie wi
thin the closed interval [0,1] for all real " }{XPPEDIT 18 0 "x" "I\"x
G6\"" }{TEXT -1 220 ". The logistic function also possesses some othe
r helpful properties, such as the ability to calculate log-odds (and, \+
as a result, odds) values with ease. These matters are very well cons
idered in Hosmer and Lemeshow, " }{TEXT 258 27 "Applied Logistic Regre
ssion" }{TEXT -1 84 ", Wiley, 1989. Note that I have in my substituti
ons used the equivalent expression " }{XPPEDIT 18 0 "exp(a+b*x[i])/(1+
exp(a+b*x[i]))" "*&-%$expG6#,&%\"aG\"\"\"*&%\"bGF(&%\"xG6#%\"iGF(F(F(,
&\"\"\"F(-F$6#,&F'F(*&F*F(&F,6#F.F(F(F(!\"\"" }{TEXT -1 67 " (obtained
by multiplying the numerator and denominator through by " }{XPPEDIT
18 0 "exp(a+b*x)" "-%$expG6#,&%\"aG\"\"\"*&%\"bGF'%\"xGF'F'" }{TEXT
-1 101 ") since Maple seems to produce more elegant results in the sub
sequent expansions and simplifications:" }}{PARA 0 "> " 0 "" {MPLTEXT
1 0 17 "lgt[i]:=a+b*x[i]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "theta[
i]:=exp(lgt[i])/(1+exp(lgt[i]));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21
"l[i]:=simplify(l[i]);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>&%&thetaG6#
%\"iG*&-%$expG6#,&%\"aG\"\"\"*&%\"bGF.&%\"xGF&F.F.F.,&F.F.F)F.!\"\"" }
}{PARA 11 "" 1 "" {XPPMATH 20 "6#>&%\"lG6#%\"iG*&,&\"\"\"F*-%$expG6#,&
%\"aGF**&%\"bGF*&%\"xGF&F*F*F*!\"\")-F,6#,&F/F4F0F4,$&%\"yGF&F4F*" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 4 "Let " }{XPPEDIT 18 0 "ll[i]" "&%#ll
G6#%\"iG" }{TEXT -1 50 " represent the log-likelihood associated with \+
the " }{XPPEDIT 18 0 "i" "I\"iG6\"" }{TEXT -1 15 "-th data point:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "ln(l[i]):" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 20 "simplify(expand(\"));" }}{PARA 0 "> " 0 "" {MPLTEXT
1 0 23 "ll[i]:=collect(\",y[i]);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#,(
-%#lnG6#,&\"\"\"F(-%$expG6#,&%\"aGF(*&%\"bGF(&%\"xG6#%\"iGF(F(F(!\"\"*
&&%\"yGF2F(F-F(F(*(F6F(F/F(F0F(F(" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>
&%#llG6#%\"iG,&*&,&%\"aG\"\"\"*&%\"bGF,&%\"xGF&F,F,F,&%\"yGF&F,F,-%#ln
G6#,&F,F,-%$expG6#F*F,!\"\"" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 339 "I
f we assume all of our data points to be mutually independent (not an \+
unreasonable assumption in most cases), the overall likelihood of a pa
rticular set of observations is the product of the individual likeliho
ods. Since the logarithm of a product is equal to the sum of the loga
rithms of the operands, we may derive the Log-Likelihood, " }{XPPEDIT
18 0 "LL" "I#LLG6\"" }{TEXT -1 26 ", by simply adding up all " }
{XPPEDIT 18 0 "n" "I\"nG6\"" }{TEXT -1 4 " of " }{XPPEDIT 18 0 "ll[i]
" "&%#llG6#%\"iG" }{TEXT -1 1 ":" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39
"LL:=unapply(sum(ll[i],i=1..n),(x,y,n));" }}{PARA 11 "" 1 "" {XPPMATH
20 "6#>%#LLG:6%%\"xG%\"yG%\"nG6\"6$%)operatorG%&arrowGF*-%$sumG6$,&*&,
&%\"aG\"\"\"*&%\"bGF5&9$6#%\"iGF5F5F5&9%F:F5F5-%#lnG6#,&F5F5-%$expG6#F
3F5!\"\"/F;;F59&F*F*" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 54 "Loading
Data and Constructing Function to be Optimized" }}{EXCHG {PARA 0 ""
0 "" {TEXT -1 21 "Hosmer and Lemeshow (" }{TEXT 269 27 "Applied Logist
ic Regression" }{TEXT -1 424 ", Wiley, 1989) introduce the problem wit
h a simple example involving 100 subjects, where the independent varia
ble is a subjects age in years and the dependent variable is whether o
r not he or she has coronary heart disease (CHD). We wish to simply t
o find a mathematical relationship to describe how one's risk for CHD \+
increases with age. This is the list of data pairs, which I have sto
red as two columns in a text file:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0
47 "#points:=readdata(`chd.txt`,[integer,integer]);" }}{PARA 12 "" 1 "
" {XPPMATH 20 "6#>%'pointsG7`q7$\"#?\"\"!7$\"#BF(7$\"#CF(7$\"#DF(7$F.
\"\"\"7$\"#EF(F17$\"#GF(F37$\"#HF(7$\"#IF(F7F7F7F77$F8F07$\"#KF(F:7$\"
#LF(F<7$\"#MF(F>7$F?F0F>F>7$\"#NF(FA7$\"#OF(7$FDF0FC7$\"#PF(7$FGF0FF7$
\"#QF(FI7$\"#RF(7$FLF07$\"#SF(7$FOF07$\"#TF(FQ7$\"#UF(FSFS7$FTF07$\"#V
F(FV7$FWF07$\"#WF(FY7$FZF0Fen7$\"#XF(7$FgnF07$\"#YF(7$FjnF07$\"#ZF(F\\
o7$F]oF07$\"#[F(7$F`oF0Fao7$\"#\\F(Fbo7$FcoF07$\"#]F(7$FfoF07$\"#^F(7$
\"#_F(7$F[pF07$\"#`F0F]p7$\"#aF07$\"#bF(7$FbpF0Fcp7$\"#cF0FdpFdp7$\"#d
F(Ffp7$FgpF0FhpFhpFhp7$\"#eF(7$FjpF0F[q7$\"#fF0F\\q7$\"#gF(7$F_qF07$\"
#hF07$\"#iF0Fcq7$\"#jF07$\"#kF(7$FhqF07$\"#lF07$\"#pF0" }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 911 "points := [[20, 0], [23, 0], [24, \+
0], [25, 0], [25, 1], [26, 0], [26, 0], [28, 0], [28, 0], [29, 0], [30
, 0], [30, 0], [30, 0], [30, 0], [30, 0], [30, 1], [32, 0], [32, 0], [
33, 0], [33, 0], [34, 0], [34, 0], [34, 1], [34, 0], [34, 0], [35, 0],
[35, 0], [36, 0], [36, 1], [36, 0], [37, 0], [37, 1], [37, 0], [38, 0
], [38, 0], [39, 0], [39, 1], [40, 0], [40, 1], [41, 0], [41, 0], [42,
0], [42, 0], [42, 0], [42, 1], [43, 0], [43, 0], [43, 1], [44, 0], [4
4, 0], [44, 1], [44, 1], [45, 0], [45, 1], [46, 0], [46, 1], [47, 0], \+
[47, 0], [47, 1], [48, 0], [48, 1], [48, 1], [49, 0], [49, 0], [49, 1]
, [50, 0], [50, 1], [51, 0], [52, 0], [52, 1], [53, 1], [53, 1], [54, \+
1], [55, 0], [55, 1], [55, 1], [56, 1], [56, 1], [56, 1], [57, 0], [57
, 0], [57, 1], [57, 1], [57, 1], [57, 1], [58, 0], [58, 1], [58, 1], [
59, 1], [59, 1], [60, 0], [60, 1], [61, 1], [62, 1], [62, 1], [63, 1],
[64, 0], [64, 1], [65, 1], [69, 1]]:" }}}{EXCHG {PARA 0 "" 0 ""
{TEXT -1 17 "Extract lists of " }{XPPEDIT 18 0 "x" "I\"xG6\"" }{TEXT
-1 5 " and " }{XPPEDIT 18 0 "y" "I\"yG6\"" }{TEXT -1 29 " data from th
e list of lists " }{XPPEDIT 18 0 "points" "I'pointsG6\"" }{TEXT -1 1 "
:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "Lx:=[seq(points[i,1],i=1..nops
(points))];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "Ly:=[seq(points[i,2]
,i=1..nops(points))];" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%#LxG7`q\"#?
\"#B\"#C\"#DF)\"#EF*\"#GF+\"#H\"#IF-F-F-F-F-\"#KF.\"#LF/\"#MF0F0F0F0\"
#NF1\"#OF2F2\"#PF3F3\"#QF4\"#RF5\"#SF6\"#TF7\"#UF8F8F8\"#VF9F9\"#WF:F:
F:\"#XF;\"#YF<\"#ZF=F=\"#[F>F>\"#\\F?F?\"#]F@\"#^\"#_FB\"#`FC\"#a\"#bF
EFE\"#cFFFF\"#dFGFGFGFGFG\"#eFHFH\"#fFI\"#gFJ\"#h\"#iFL\"#j\"#kFN\"#l
\"#p" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%#LyG7`q\"\"!F&F&F&\"\"\"F&F&
F&F&F&F&F&F&F&F&F'F&F&F&F&F&F&F'F&F&F&F&F&F'F&F&F'F&F&F&F&F'F&F'F&F&F&
F&F&F'F&F&F'F&F&F'F'F&F'F&F'F&F&F'F&F'F'F&F&F'F&F'F&F&F'F'F'F'F&F'F'F'
F'F'F&F&F'F'F'F'F&F'F'F'F'F&F'F'F'F'F'F&F'F'F'" }}}{EXCHG {PARA 0 ""
0 "" {TEXT -1 55 "We now have enough information to construct a functi
on " }{XPPEDIT 18 0 "f(a,b" "-%\"fG6$%\"aG%\"bG" }{TEXT -1 45 " to be \+
minimized. Since we wish to maximize " }{XPPEDIT 18 0 "LL" "I#LLG6\"
" }{TEXT -1 58 " in the given circumstances, we substitute our data li
sts " }{XPPEDIT 18 0 "Lx" "I#LxG6\"" }{TEXT -1 5 " and " }{XPPEDIT 18
0 "Ly" "I#LyG6\"" }{TEXT -1 101 " and take the negative of the result \+
in order to obtain such a function, since maximizing a function " }
{TEXT 263 2 "-f" }{TEXT -1 27 " is the same as minimizing " }{TEXT
262 1 "f" }{TEXT -1 17 ", and vice versa:" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 20 "-LL(Lx,Ly,nops(Lx)):" }}{PARA 0 "> " 0 "" {MPLTEXT 1
0 20 "f:=unapply(\",(a,b));" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%\"fG:
6$%\"aG%\"bG6\"6$%)operatorG%&arrowGF),fp9$!#V9%!%0A-%#lnG6#,&\"\"\"F6
-%$expG6#,&F.F6F0\"#?F6F6-F36#,&F6F6-F86#,&F.F6F0\"#BF6F6-F36#,&F6F6-F
86#,&F.F6F0\"#CF6F6-F36#,&F6F6-F86#,&F.F6F0\"#DF6\"\"#-F36#,&F6F6-F86#
,&F.F6F0\"#EF6FQ-F36#,&F6F6-F86#,&F.F6F0\"#GF6FQ-F36#,&F6F6-F86#,&F.F6
F0\"#HF6F6-F36#,&F6F6-F86#,&F.F6F0\"#IF6\"\"'-F36#,&F6F6-F86#,&F.F6F0
\"#KF6FQ-F36#,&F6F6-F86#,&F.F6F0\"#LF6FQ-F36#,&F6F6-F86#,&F.F6F0\"#MF6
\"\"&-F36#,&F6F6-F86#,&F.F6F0\"#NF6FQ-F36#,&F6F6-F86#,&F.F6F0\"#OF6\"
\"$-F36#,&F6F6-F86#,&F.F6F0\"#PF6F]r-F36#,&F6F6-F86#,&F.F6F0\"#QF6FQ-F
36#,&F6F6-F86#,&F.F6F0\"#RF6FQ-F36#,&F6F6-F86#,&F.F6F0\"#SF6FQ-F36#,&F
6F6-F86#,&F.F6F0\"#TF6FQ-F36#,&F6F6-F86#,&F.F6F0\"#UF6\"\"%-F36#,&F6F6
-F86#,&F.F6F0\"#VF6F]r-F36#,&F6F6-F86#,&F.F6F0\"#WF6Fht-F36#,&F6F6-F86
#,&F.F6F0\"#XF6FQ-F36#,&F6F6-F86#,&F.F6F0\"#YF6FQ-F36#,&F6F6-F86#,&F.F
6F0\"#ZF6F]r-F36#,&F6F6-F86#,&F.F6F0\"#[F6F]r-F36#,&F6F6-F86#,&F.F6F0
\"#\\F6F]r-F36#,&F6F6-F86#,&F.F6F0\"#]F6FQ-F36#,&F6F6-F86#,&F.F6F0\"#^
F6F6-F36#,&F6F6-F86#,&F.F6F0\"#_F6FQ-F36#,&F6F6-F86#,&F.F6F0\"#`F6FQ-F
36#,&F6F6-F86#,&F.F6F0\"#aF6F6-F36#,&F6F6-F86#,&F.F6F0\"#bF6F]r-F36#,&
F6F6-F86#,&F.F6F0\"#cF6F]r-F36#,&F6F6-F86#,&F.F6F0\"#dF6Fho-F36#,&F6F6
-F86#,&F.F6F0\"#eF6F]r-F36#,&F6F6-F86#,&F.F6F0\"#fF6FQ-F36#,&F6F6-F86#
,&F.F6F0\"#gF6FQ-F36#,&F6F6-F86#,&F.F6F0\"#hF6F6-F36#,&F6F6-F86#,&F.F6
F0\"#iF6FQ-F36#,&F6F6-F86#,&F.F6F0\"#jF6F6-F36#,&F6F6-F86#,&F.F6F0\"#k
F6FQ-F36#,&F6F6-F86#,&F.F6F0\"#lF6F6-F36#,&F6F6-F86#,&F.F6F0\"#pF6F6F)
F)" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 38 "Implementing an Optimizat
ion Procedure" }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 16 "Loading the Code
" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 18 "Finding values of " }{XPPEDIT
18 0 "a" "I\"aG6\"" }{TEXT -1 5 " and " }{XPPEDIT 18 0 "b" "I\"bG6\""
}{TEXT -1 139 " that minimize the above function is the essence of thi
s problem. We could take the \"classical\" approach and compute first
derivatives of " }{XPPEDIT 18 0 "f" "I\"fG6\"" }{TEXT -1 17 " with re
spect to " }{XPPEDIT 18 0 "a" "I\"aG6\"" }{TEXT -1 5 " and " }
{XPPEDIT 18 0 "b" "I\"bG6\"" }{TEXT -1 269 ", set these to zero, and s
olve the system of nonlinear equations that results. Unfortunately, t
his is a gargantuan computational challenge. It is less cumbersome to
work with the single function and attempt to find minima with any one
of a number of algorithms. The " }{TEXT 264 17 "Numerical Recipes" }
{TEXT -1 28 " books (I have Press et al, " }{TEXT 265 27 "Numerical Re
cipes in Pascal" }{TEXT -1 302 ", CUP, 1989) offer a host of routines,
and very recently Francis Wright of the University of London has port
ed the downhill simplex minimization (DSM) procedure into Maple for ou
r benefit. For a discussion of DSM, refer to this sheet (via Dr. Wrig
ht's website http://www.maths.qmw.ac.uk/~fjw) and the " }{TEXT 266 17
"Numerical Recipes" }{TEXT -1 130 " book of your choice. I have simpl
y extracted the code from Dr. Wright's sheet and saved it as an `*.m` \+
file, which we load now: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "#read \+
`downhill.m`;\nread `f:\\\\depart\\\\math\\\\maple\\\\misc\\\\downhill
.m`;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 55 "You can print out the cod
e of the relevant procedures, " }{TEXT 271 6 "amoeba" }{TEXT -1 5 " an
d " }{TEXT 272 7 "amoeba1" }{TEXT -1 278 ", by executing the following
. To make the file as compact as possible, I have removed most of Dr.
Wright's original comments, so I once again refer you to his workshee
t for further explication. I have only removed the comments--the work
ing code itself is completely unaltered:" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 13 "eval(amoeba);" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#:6&'
%\"fG%*procedureG'%&startG-%%listG6#%(numericG'%&scaleGF-'%%ftolGF-6)%
%ndimG%%fargG%\"pG%\"yG%%iterG%$iloG%\"jG6\"F:C.>8$-%%nopsG6#9%@%/-F?6
#7#-%#opG6$\"\"\"-%%evalG6#9$FJ>%,amoeba~funcGFK>FP-%(unapplyG6$-.FN6#
-%$seqG6$&8%6#8*/Fhn;FJF=Ffn@$%-amoeba~debugG-%&printG6#FP>8&7$FA-FY6$
-%'subsopG6$/Fhn,&&FAFgnFJ9&FJFAFin>8'-%'vectorG6#-%$mapG6$FPFao>Fao-%
'matrixG6#Fao>8(-%&arrayG6#;FJFJ>8)-%&roundG6#-%'evalhfG6#-%(amoeba1G6
'-%$varGFgp-Fjq6#F]pF=9'-Fjq6#Fip>%,amoeba~iterG-Faq6#&Fip6#FJ@$/F_q!
\"\"-%&ERRORG6#%5too~many~iterations.G@$F\\o-F^o6&F_qFaoF]pFar7#-FY6$&
Fao6$F_qFhnFinF:6$FPFar" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "
eval(amoeba1);" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#:6'%\"pG%\"yG%%ndimG
%%ftolG%,iter_returnG6/%%mptsG%\"jG%%inhiG%$iloG%$ihiG%\"iG%%yprrG%$yp
rG%%rtolG%#prG%$prrG%%pbarG%%iterG6\"F8C*>8--%&arrayG6#;\"\"\"9&>8.F<>
8/F<>8$,&FAF@F@F@>80\"\"!?(F8F@F@F8%%trueGC.>8'F@@%2&9%6#\"\"#&FT6#F@C
$>8(F@>8&FVC$>FenFV>FgnF@?(8)F@F@FGFMC$@$2&FT6#F\\o&FT6#FP>FPF\\o@&2&F
T6#FenF`oC$>FgnFen>FenF\\o2&FT6#FgnF`o@$0F\\oFen>FgnF\\o>8,,$*&-%$absG
6#,&FgoF@Fbo!\"\"F@,&-Fgp6#FgoF@-Fgp6#FboF@Fjp$\"#?Fjp@$2Fcp9'%&breakG
@$/FJ\"$+&C$>&9(FXFJ-%'RETURNG6#Fjp>FJ,&FJF@F@F@?(8%F@F@FAFM>&FE6#FcrF
K?(F\\oF@F@FGFM@$F`p?(FcrF@F@FAFM>Fer,&FerF@&9$6$F\\oFcrF@?(FcrF@F@FAF
MC$>Fer*&FerF@FAFjp>&F;Ffr,&FerF`q&F]s6$FenFcr$!#5Fjp>8+-%,amoeba~func
G6#F;@'1F[tFboC%?(FcrF@F@FAFM>&FCFfr,&FdsF`qFerFhs>8*-F]t6#FC@%2FgtFbo
C$?(FcrF@F@FAFM>FfsFdt>FgoFgtC$?(FcrF@F@FAFM>FfsFds>FgoF[t1F]pF[tC&@$2
F[tFgoC$?(FcrF@F@FAFM>FfsFds>FgoF[t?(FcrF@F@FAFM>Fdt,&Ffs$\"\"&FjpFerF
_v>FgtFht@%2FgtFgoC$?(FcrF@F@FAFM>FfsFdt>FgoFgt?(F\\oF@F@FGFM@$0F\\oFP
C$?(FcrF@F@FAFMC$>Fds,&F\\sF_v&F]s6$FPFcrF_v>F\\sFds>F`oF\\tC$?(FcrF@F
@FAFM>FfsFds>FgoF[t>F[rFJFPF8F8" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1
16 "Running the Code" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 29 "The DSM in
terface procedure, " }{TEXT 267 6 "amoeba" }{TEXT -1 53 ", needs to be
passed the function name (in this case " }{XPPEDIT 18 0 "f" "I\"fG6\"
" }{TEXT -1 348 "), a list or vector of starting guesses of the parame
ters to be estimated, a so-called scale factor estimating the size of \+
the initial simplex (the geometric region in which we start seeking ou
r function minimum), and a small number representing the relative accu
racy to which we want the routine to converge. Be careful with this l
ast parameter--" }{TEXT 268 6 "amoeba" }{TEXT -1 199 " may not converg
e if it is too small. The procedure is very fast if it is going to wo
rk, so if you don't get any output within a few seconds, hit stop and \+
and alter starting parameters as necessary." }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 604 "We are finally ready to find an answer. I have rather
arbitrarily chosen the origin as a starting point and unity as the sc
ale value--experiment with these to your heart's content. I have foun
d that the tolerance value I have selected strikes a nice balance betw
een speed and accuracy--in this sort of work, estimating the parameter
s to more than a few decimal places is rarely necessary. The first li
ne of output represents the parameter estimates, with the second line \+
being the estimated Log Likelihood and the third the number of iterati
ons the algorithm underwent to achieve the desired accuracy:" }}{PARA
0 "> " 0 "" {MPLTEXT 1 0 13 "start:=[0,0]:" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 123 "result:=amoeba(f,start,1,1e-14):\nsol:=[a=result[1],
b=result[2]]:sol;\n`Log-Likelihood`=-`amoeba func`(result);`amoeba ite
r`;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7$/%\"aG$!1&p7Xj_%4`!#:/%\"bG$
\"18`o98@46!#;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%/Log-LikelihoodG$!
+KYln`!\")" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"#f" }}}}{SECT 1 {PARA
4 "" 0 "" {TEXT -1 10 "The Result" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1
83 "Now we produce the maximum-likelihood-estimated logistic fit to th
e given data set:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "eq:=P=subs(sol
,1/(1+exp(-(a+b*x)))):eq;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%\"PG*$,
&\"\"\"F'-%$expG6#,&$\"1&p7Xj_%4`!#:F'%\"xG$!18`o98@46!#;F'!\"\"" }}}}
}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 30 "Plot and Discussion of Results"
}}{PARA 0 "" 0 "" {TEXT -1 71 "Let's plot the curve along with the dat
a points to see how things look:" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT
1 0 109 "plot(rhs(eq),x=18..70):crvplt:=\":\nplot(points,style=point,s
ymbol=box):ptplt:=\":\nplots[display](crvplt,ptplt);" }}{PARA 13 "" 1
"" {INLPLOT "6&-%'CURVESG6$7S7$$\"#=\"\"!$\"1fI84DA8N!#<7$$\"1nmm,_M8>
!#9$\"1&[;)[AAlRF-7$$\"1LLebi'>,#F1$\"1<=@_UT.WF-7$$\"1nm;Xf(G7#F1$\"1
PdG\"o@8&\\F-7$$\"1nm;;3_MAF1$\"1D@*fXAxc&F-7$$\"1LLee]jXBF1$\"1;p.%*)
eBD'F-7$$\"1mm\"\\=_'[CF1$\"1%zPhUpl&pF-7$$\"1++D_/KbDF1$\"16G3].]ixF-
7$$\"1nm\"p#pjlEF1$\"1ol(p-c^o)F-7$$\"1++D:'*fvFF1$\"1\\C\\t/]-(*F-7$$
\"1LLLc(4())GF1$\"13=h$Gje3\"!#;7$$\"1nm;IsL))HF1$\"1CDAmpa(>\"F[o7$$
\"1+++pT\\+JF1$\"1@&Hx'H,N8F[o7$$\"1+++N;68KF1$\"1\")*p5vTi[\"F[o7$$\"
1+++J%R;K$F1$\"1I%)*)QB4X;F[o7$$\"1mm\"p1$>?MF1$\"1!y'\\]I!4!=F[o7$$\"
1LLLsGQPNF1$\"1U-*=JF[o7$$\"1LL$e
fLZ=%F1$\"1*QQc8B,R$F[o7$$\"1LLeK-1(G%F1$\"1o;i7C\"*[OF[o7$$\"1mm;MLV(
R%F1$\"1'461$z.PRF[o7$$\"1nm\"p'*z?^%F1$\"1%*fSavKWUF[o7$$\"1++v(>!)=h
%F1$\"1G&=x!ev;XF[o7$$\"1mm;\"=n'>ZF1$\"1**=`pXA9[F[o7$$\"1*****>o@5$[
F1$\"1Rm[HM!H7&F[o7$$\"1++]O.'*R\\F1$\"1l3ARu*RU&F[o7$$\"1++D'3k`/&F1$
\"1*p*)z$\\U7dF[o7$$\"1****\\zpRi^F1$\"1!4vHh?q-'F[o7$$\"1LLLRibn_F1$
\"1`LG&y6FI'F[o7$$\"1+++&>O)z`F1$\"1*)*Q_iNze'F[o7$$\"1mm\"H_y:[&F1$\"
1j*z$fz)o$oF[o7$$\"1+++h(4Gf&F1$\"1Etm>#pu4(F[o7$$\"1KLe(opup&F1$\"1pO
.T8oItF[o7$$\"1++Dn%po!eF1$\"1%\\%QIqKhvF[o7$$\"1LLL[x#Q\"fF1$\"1\")o]
=-StxF[o7$$\"1***\\i:.e-'F1$\"1&Rw@wf4)zF[o7$$\"1nmm#y[O8'F1$\"1.XVQ(o
o;)F[o7$$\"1mm;7l$RC'F1$\"1K&))f&*pHM)F[o7$$\"1LLe.5J`jF1$\"1:$yM!>'R]
)F[o7$$\"1*****>![\"QX'F1$\"1n\")>#z@.k)F[o7$$\"1KL$o\\.!plF1$\"1-R&oi
nNy)F[o7$$\"1mmmQ'H?n'F1$\"1NLh.l[+*)F[o7$$\"1++D%*p(=y'F1$\"1wfA'[\"=
9!*F[o7$$\"1,+vH\\,()oF1$\"1yS^OF28\"*F[o7$$\"#qF*$\"1*y$QyDG4#*F[o-%'
COLOURG6&%$RGBG$\"#5!\"\"F*F*-F$6&7`q7$$\"#?F*F*7$$\"#BF*F*7$$\"#CF*F*
7$$\"#DF*F*7$F^\\l$\"\"\"F*7$$\"#EF*F*Fc\\l7$$\"#GF*F*Ff\\l7$$\"#HF*F*
7$$\"#IF*F*F\\]lF\\]lF\\]lF\\]l7$F]]lFa\\l7$$\"#KF*F*F`]l7$$\"#LF*F*Fc
]l7$$\"#MF*F*Ff]l7$Fg]lFa\\lFf]lFf]l7$$\"#NF*F*Fj]l7$$\"#OF*F*7$F^^lFa
\\lF]^l7$$\"#PF*F*7$Fb^lFa\\lFa^l7$$\"#QF*F*Fe^l7$$\"#RF*F*7$Fi^lFa\\l
7$$\"#SF*F*7$F]_lFa\\l7$$\"#TF*F*F`_l7$$\"#UF*F*Fc_lFc_l7$Fd_lFa\\l7$$
\"#VF*F*Fg_l7$Fh_lFa\\l7$$\"#WF*F*F[`l7$F\\`lFa\\lF^`l7$$\"#XF*F*7$F``
lFa\\l7$$\"#YF*F*7$Fd`lFa\\l7$$\"#ZF*F*Fg`l7$Fh`lFa\\l7$$\"#[F*F*7$F\\
alFa\\lF^al7$$\"#\\F*F*F_al7$F`alFa\\l7$$\"#]F*F*7$FdalFa\\l7$$\"#^F*F
*7$$\"#_F*F*7$F[blFa\\l7$$\"#`F*Fa\\lF^bl7$$\"#aF*Fa\\l7$$\"#bF*F*7$Fe
blFa\\lFgbl7$$\"#cF*Fa\\lFhblFhbl7$$\"#dF*F*F[cl7$F\\clFa\\lF^clF^clF^
cl7$$\"#eF*F*7$F`clFa\\lFbcl7$$\"#fF*Fa\\lFccl7$$\"#gF*F*7$FgclFa\\l7$
$\"#hF*Fa\\l7$$\"#iF*Fa\\lF]dl7$$\"#jF*Fa\\l7$$\"#kF*F*7$FddlFa\\l7$$
\"#lF*Fa\\l7$$\"#pF*Fa\\lFjz-%&STYLEG6#%&POINTG-%'SYMBOLG6#%$BOXG-%+AX
ESLABELSG6$%\"xG%!G-%%VIEWG6$;F(Ffz%(DEFAULTG" 2 612 296 296 2 0 1 0
2 9 0 4 2 1.000000 45.000000 45.000000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 1 0 0 0 0 5 0 0 0 0 0 0 }}}{PARA 0 "" 0 "" {TEXT -1 892 "This gr
aphical depiction should give you an idea of what a logistic fit is an
d what it signifies. On examining the raw data, we notice a tendency \+
for the 0's to \"bunch up\" in the lower age ranges, whereas the 1's a
ppear to be more prevalent with older age. The fitted curve represent
s the estimated probability that, at a given age, one will be a 1 (i.e
., have CHD) rather than a 0. Note that this differs from the usual r
egression task of estimating the value of the dependent variable for a
particular value of independent variable. In logistic regression (as
in the loosely related discriminant function analysis), we are intere
sted instead in the likelihood of subject belonging to a particular gr
oup (in this case, diseased vs. not diseased).\nI will not get into a \+
lengthy discussion of the meaning of the parameters, leaving it to you
to experiment here. It should be evident that " }{XPPEDIT 18 0 "a" "
I\"aG6\"" }{TEXT -1 28 " determines where along the " }{XPPEDIT 18 0 "
x" "I\"xG6\"" }{TEXT -1 55 "-axis that the curve begins to swoop up (o
r down), and " }{XPPEDIT 18 0 "b" "I\"bG6\"" }{TEXT -1 160 " specifies
how steeply it does so. In more complex multiple logistic regression
analyses, where there is more than one independent variable and there
fore more " }{XPPEDIT 18 0 "b" "I\"bG6\"" }{TEXT -1 15 "-weights, each
" }{XPPEDIT 18 0 "b" "I\"bG6\"" }{TEXT -1 318 " indicates the rate at
which the log-odds of a positive outcome changes for a unit change in
the associated independent variable. This very important property al
lows for the estimation of \"risk\" associated with each independent v
ariable. See Hosmer and Lemeshow or Kleinbaum, Kupper, and Muller for
more information." }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 15 "Closing Re
marks" }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 26 "Estimating initial guess
es" }}{PARA 0 "" 0 "" {TEXT -1 544 "At present, I know very little abo
ut the convergence properties of the DSM method when applied to this p
articular problem, and I encourage the users of this sheet to experime
nt with this. I believe that a general rule of thumb is that good ini
tial guesses will converge to the desired maximum Log-Likelihood. Hos
mer and Lemeshow demonstrate that by grouping the CHD data according t
o age ranges one can perform a linear regression on the logit transfor
mation given earlier, where the independent variable is the midpoint o
f the age range and " }{XPPEDIT 18 0 "P" "I\"PG6\"" }{TEXT -1 558 " is
the estimated frequency of CHD for that age range. This approach may
not be feasible in some cases, especially if one must make the ranges
very broad to avoid proportions of 0 or 1. Also, one always loses th
e \"fine-grained\" nature of data when subjecting it to arbitrary grou
ping. Nonetheless, this approach may produce good \"ballpark\" initia
l guesses for the parameters. Even still, as Francis Wright points ou
t, the algorithm may elude convergence even when the initial guess is \+
close to exact. This question clearly warrants closer investigation.
" }}}{SECT 0 {PARA 4 "" 0 "" {TEXT -1 38 "Extension to the multidimens
ional case" }}{PARA 0 "" 0 "" {TEXT -1 87 "The above is but a simple e
xample, with only one independent variable. In general the " }{TEXT
257 5 "logit" }{TEXT -1 123 " can comprise a linear combination of mor
e than one independent variable, leading to a multidimensional problem
of the form" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "P = 1/(1+exp
(-(b[0]+sum(b[j]*x[j],j = 1 .. m))));" }}{PARA 11 "" 1 "" {XPPMATH 20
"6#/%\"PG*$,&\"\"\"F'-%$expG6#,&&%\"bG6#\"\"!!\"\"-%$sumG6$*&&F-6#%\"j
GF'&%\"xGF6F'/F7;F'%\"mGF0F'F0" }}}{PARA 0 "" 0 "" {TEXT -1 16 "where \+
there are " }{XPPEDIT 18 0 "m" "I\"mG6\"" }{TEXT -1 44 " independent v
ariables. Finding the set of " }{XPPEDIT 18 0 "m+1" ",&%\"mG\"\"\"\"
\"\"F$" }{TEXT -1 1 " " }{XPPEDIT 18 0 "b" "I\"bG6\"" }{TEXT -1 58 "-w
eights that maximizes the Log Likelihood is the task of " }{TEXT 256
8 "multiple" }{TEXT -1 389 " logistic regression, a statistical techni
que often employed in large-scale epidemiological research where the c
ontributions of a number of risk factors to the development of a parti
cular pathological state need to be delineated. Hosmer and Lemeshow (
1989) is the standard reference for this sort of work. I will endeavo
ur to provide a multidimensional example in a subsequest worksheet." }
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "0 0 0" 54 }
{VIEWOPTS 1 1 0 3 4 1802 }