(*:Mathematica:: V2.2.1 *) (*:Context: "clmnllr`" *) (*:CLM Version 1.1 *) (*:Title:: log-log regression*) (*:Summary:: LLR: log-log regression: log(-log(p)) *) (*:References:: Prentice, R.L. and Gloeckler, L.A. Biometrics, 34, 57-67.*) (*:Keywords:: log-log regression *) (*:Requirements:: matrixlx`: matrix funcitons *) (*:Author: 1993 Stuart G. Baker *) BeginPackage["clm`cdxnllr`","clm`matrixlx`"] NLLR::usage= "NLLR[ys,{m,x,k:_0},options] or NLLR[ys,{x,k:_0},options] returns the parameters for regression in models for log(-log(p)). The argument ys is the complete data, m is a matrix such that m . ys is a column matrix which represents successive pairs of failure and success, x is the design matrix, and k is a constant multiplying counts. The options are MaxIterationsNLLR,ConvergenceGoalNLLR,WorkingPrecisionNLLR, and ShowProgressNLLR." Options[NLLR]={MaxIterationsNLLR->10,ConvergenceGoalNLLR->.001, WorkingPrecisionNLLR->8,ShowProgressNLLR->False} Clear[NLLR] Begin["Private`"] NLLR[ys_,{x_,k_:0},options___Rule]:= Module[{iter,conv,prec,showprog}, {iter,conv,prec,showprog}= {MaxIterations,ConvergenceGoal,WorkingPrecision,ShowProgress}/. {options}/.Options[NLLR]; pNLLR[ys,{x,k},{iter,conv,prec,showprog}]] NLLR[ys_,{m_,x_,k_:0},options___Rule]:= NPR[m.ys,{x,k},options] (*core of NLLR iteration*) pNLLR[ys_,{x_,k_},{maxiter_,conv_,prec_,showprog_}]:= Module[{y0,y1,t0,arg,t2}, t0=pNLLRStart[ys,x,k]; {y0,y1}=pSplitStart[ys]; arg={y0,y0+y1,x,k}; t2=FixedPoint[ pNLLRStep[arg,#,showprog]&, t0, maxiter, SameTest->((Max[Abs[(#1-#2)/#1]] < conv )&)]; Return[t2]] (*starting value computed using transformation of observed proportion*) pNLLRStart[ys_,x_,k_]:= Module[{y0,y1, ya0,ya1,p,ya, xt,t0}, {y0,y1}=pSplitStart[ys]; ya0=y0+(1-Abs[Sign[y0]]) .01; ya1=y1+(1-Abs[Sign[y1]]) .01; p=ya1/(ya0+ya1); ya=Log[-Log[p]-k]; xt=Transpose[x]; t0=Inverse[xt.x].xt.ya; Return[t0]] (*one step in iteration*) pNLLRStep[{y_,n_,x_,k_},t1_,showprog_]:= Module[{eta,h,p,xt,part1,part2,s,i,t2,tmaxdif}, eta=x.t1; h=Exp[eta]; p=Exp[-h-k]; xt=Transpose[x]; part1=h (y/(1-p)-n); part2=y p h^2 / (1-p)^2; s=xt .part1; i=xt.Diag[part2-part1].x; t2=t1+pLinearSolve[i,s]; If[showprog, tmaxdif=Max[Abs[(t1-t2)/t1]]; Print[" tmaxdif=",tmaxdif]]; Return[t2]] (*--------pLinearSolve -----------------------------------*) pLinearSolve[mat1_,mat2_]:= Transpose[{LinearSolve[mat1,Flatten[mat2]]}] (*--------pSplitStart---------------------------------------------*) (*pSplitStart takes y, a column matrix which alternates numbers of successes and failures and returns y0 numbers of success and y1 numbers of failures. *) pSplitStart[y_]:= Module[{k,index0,index1, y0,y1, yt0,yt1,res}, k=Length[y]/2; index0=(2 Range[k])-1; index1=2 Range[k]; y0=y[[index0]] //N; y1=y[[index1]] //N; res={y0,y1}; Return[res]] End[] (*Protect[NLLR]*) EndPackage[]