(*:Mathematica:: V2.2.1 *) (*:Context: "cdxnlr`" *) (*:Version 1.1 *) (*:Title:: Logistic regression *) (*:Summary:: NLR: logistic regression *) (*References:: Agresti, A. 1990. Categorical Data Analysis. New York: John Wiley and Sons, pp 115-117. *) (*Keywords:: logistic regression *) (*Input For:: cdxem`*) (*Requirements:: matrixlx`: matrix funcitons *) (*History: 1993 Stuart G. Baker *) BeginPackage["clm`cdxnlr`","clm`matrixlx`"] NLR::usage= "NLR[ys,{m,x},options__Rule] or NLR[ys,{x},options___Rule] returns the logistic regression parameters. The argument ys is the complete data, m is a matrix such that m . y is a column matrix with sucessive pairs of counts representing success and failure, and x is the design matrix. The options are MaxIterationsNLR,ConvergenceGoalNLR,WorkingPrecisionNLR, and ShowProgressNLR" Options[NLR]={MaxIterationsNLR->10,ConvergenceGoalNLR->.0001, WorkingPrecisionNLR->8,ShowProgressNLR->False} Clear[NLR] Begin["Private`"] NLR[ys_,{x_},options___Rule]:= Module[{iter,conv,prec,showprog}, {iter,conv,prec,showprog}= {MaxIterationsNLR,ConvergenceGoalNLR, WorkingPrecisionNLR,ShowProgressNLR}/. {options}/.Options[NLR]; pNLR[ys,{x},{iter,conv,prec,showprog}]] NLR[ys_,{m_,x_},options___Rule]:= NLR[m.ys,{x},options] (*core of logistic regression*) pNLR[ys_,{x_},{maxiter_,conv_,prec_,showprog_}]:= Module[{t0,y0,y1,n,arg,tlist,t2}, t0=pNLRStart[ys,x]; {y0,y1}=pSplitStart[ys]; n=y0+y1; arg={y0,n,x}; t2=FixedPoint[ pNLRStep[arg,#,showprog]&, t0, maxiter, SameTest->((Max[Abs[(#1-#2)/ReplaceZero[#1,.001]]] < conv )&)]; Return[t2]] (*starting value computed using emprical logits*) pNLRStart[ys_,x_]:= Module[{y0,y1, z0,z1,z, p,w,n, xt,t0,arg}, {y0,y1}=pSplitStart[ys]; (*add .01 to any 0 elements*) z0=ReplaceZero[y0,.01]; z1=ReplaceZero[y1,.01]; n=z0+z1; p=z0/n; w=Diag[n p (1-p)]; z=Log[z0/z1]; xt=Transpose[x]; t0=pLinearSolve[xt.w.x,xt.w.z]; Return[t0]] (*logistic regression step*) pNLRStep[{y_,n_,x_},t1_,showprog_]:= Module[{eta,xt,p,s,i,t2,tmaxdif}, eta=x.t1; xt=Transpose[x]; p=Exp[eta]/(1+Exp[eta]); s=xt .(y-n p); i=xt.Diag[n p (1-p)].x; t2=t1+pLinearSolve[i,s]; If[showprog, tmaxdif=Max[Abs[(t1-t2)/ReplaceZero[t1,.001]]]; Print[" tmaxdif=",tmaxdif]]; Return[t2]] pLinearSolve[mat1_,mat2_]:= Transpose[{LinearSolve[mat1,Flatten[mat2]]}] (*pSplitStart takes y, a column matrix which alternates numbers of successes and failures and returns {column matrix of success, column matrix of failures} *) pSplitStart[y_]:=Transpose @ Partition[y,2] End[] (*Protect[NPS,NPSLogit,NLR,NPR,NLLR]*) EndPackage[]