



(*:Mathematica:: V2.1 *)



(*:Context: "psederiv`" *)          



(*Title: periodic screeening *)



(*References:   Baker, S.G. and Chu, K.C.*)



(*Keywords:  periodic screening evaluation: incidence data *)



(*Requirements:  matrixlx *)



(*History: 1998 Stuart G. Baker *)





BeginPackage["psederiv`","matrixlx`"]





CreateMList::usage="CreateMListSurv[h,ageRange]"

CreateMListAvg::usage="CreateMListAvgSurv[h,ageRange]"

DerivQ::usage="DerivQ[{q,a},mstar,ageRange]" 

DerivH::usage="DerivH[qlist,{h,n},ageRange]"

DerivHAvg::usage="DerivHAvg[qlist,{h,n}]"

FunctionToList::usage="FunctionToList[f,alist,t]"





Clear[CreateMList,CreateMListAvg,

       DerivQ,DerivH,DerivHAvg,

       FunctionToList]

        

       

       

Begin["Private`"] 

 

    

 

 

(*-------------------------CreateMList-----------------------*)



CreateMList[h_,{startage_,endage_}]:= 

  Module[{mat},

    mat=CreateMListSurv[h,#]& /@ range[{startage,endage}];

    list=Apply[Plus,mat,1];

 Return[list]]





CreateMListAvg[h_,{startage_,endage_}]:= 

 Module[{mat},

    mat=CreateMListAvgSurv[h,#]& /@ range[{startage,endage}];

    list=Apply[Plus,mat,1];

 Return[list]]





CreateMListAvgSurv[h_,a_]:=

   Module[{h1,mlist1,mlist2,mlist},

     h1=Join[Drop[h,-1],{0}];

     mlist1=CreateMListSurv[h,a];

     mlist2=CreateMListSurv[h1,a+1];

     mlist=(mlist1+mlist2)/2;

  Return[mlist]]





CreateMListSurv[h_,a_]:= 

  Module[{n},

    n=Length[h];

    mlist= FoldList[Times,1,(1-h)] Join[h,{0}] Surv[a,n];

  Return[mlist]]

                                                            

                                                            



(*------------Deriv of Q x mlist----------------------*)





DerivQ[{q_,a_},mstar_,ageRange_]:=

  Module[{dq,dqmat,mlistage,b,dmlistmat},

   dq=derivq[q,a];

   dqmat=FunctionToList[{dq,a},ageRange];      (* b x A *)

   b=Length[dqmat];                            (* b=number of parameters*)

   dmlistmat=dqmat (J[b,1] . {mstar});         (* b x A *)

   dmlistq=Apply[Plus,dmlistmat,1];            (* b     *)

 Return[dmlistq]]





derivq[q_,a_]:=

 Module[{len},

  len=Length@q[[2]];

  der=a^ (Range[len]-1) q;

 Return[der]]





(*--------------Deriv of hazard----------------------------*)





DerivH[qlist_,hlist_,ageRange_]:=

  Module[{dmat,mat,derivlist},

   dmat=DerivMList[hlist,ageRange];     (* A x T *)

   mat=Hdir[Transpose[{qlist}],dmat];   (* A x T *)

   derivhlist=Apply[Plus,mat];          (* T  *)

Return[derivhlist]]   





DerivHAvg[qlist_,hlist_,ageRange_]:=

  Module[{hlist1,dmat,mat,derivlist},   

   hlist1=Join[Drop[hlist,1],{0}];

   dmat=(DerivMList[hlist,ageRange] +

         DerivMList[hlist1,ageRange+1])/2;                    (* A x T *)

   mat=Hdir[Transpose[{qlist}],dmat];                         (* A x T *)  

   derivhlist=Apply[Plus,mat];                                (* T     *)

Return[derivhlist]]                  







DerivMList[hlist_,{startage_,endage_}]:=

  Module[{derivMmat,derivlist},  

   derivMmat=DerivMMatrix[hlist];                (*only need to compute once*)

   derivAgebyPar= DerivMListSurv[hlist,derivMmat,#]& /@ 

                   range[{startage,endage}];

Return[derivAgebyPar]]









DerivMListSurv[hlist_,derivMmat_,a_]:=

 Module[{n,surv,mat,deriv},

  n=Length[hlist];

  surv=Surv[a,n];                   (* T+1     *)

  mat=derivMmat (J[n,1] . {surv});  (* T x T+1 *)

  deriv=Apply[Plus,mat,1];          (* T       *)

 Return[deriv]]





    (*derivative without competing risks*)

   

DerivMMatrix[hlist_]:=

 Module[{mlist,n,a,b,derivMmat},

   mlist=FoldList[Times, 1, (1-hlist)] Join[hlist,{0}];

   n=Length[hlist];

   derivMmat=Table[DerivM[mlist,hlist,{a,b}],{a,1,n},{b,1,n+1}];   (*T x T+1*)

 Return[derivMmat]]



                

DerivM[mlist_,hlist_,{a_,b_}]:=

                         -mlist[[b]]/(1-hlist[[a]])/; a < b



DerivM[mlist_,hlist_,{a_,b_}]:=

                        mlist[[a]]/hlist[[a]]/; a == b && hlist[[a]]>0



DerivM[mlist_,hlist_,{a_,b_}]:= 0/; a == b && hlist[[a]]==0



DerivM[mlist_,hlist_,{a_,b_}]:= 0/; a > b





(*-------------Survive Competing Risks-----------------*)



(* a = 0 is 40, a=1 is 41 ....*)          



Surv[a0_,n_]:=

 Module[{a,s0,s,y0},

  a=a0+1;                                  (*starts with a=0 is 40*)

  s0={.998,

      .998, .998, .998, .998, .997, 

      .997, .997, .996, .996, .996, 

      .995, .995, .994, .994, .993, 

      .993, .991, .992, .991, .991, 

      .989, .988, .987, .986, .984, 

      .983, .982, .981, .979, .978, 

      .976, .974, .972, .969, .966, 

      .963, .959, .954, .949, .943};      (*yearly survival, females, 40*)

   s=Take[s0,{a,a+n}];

   surv=Drop[FoldList[Times,1,s],-1];      (*n = T+1*)

 Return[surv]]











   

                           

(*-----------------------Support Functions---------------------*)



 





FunctionToList[{f_,a_},ageRange_]:=

 Module[{alist,flist,sum,flist1,s0,s1,s2,s3},

     alist=range[ageRange];

      flist=f/.a->alist;

    Return[flist]] 

 



 

 (*start at age 0: same as  for Poisson regression*)

 

  range[{startage_,endage_}]:= Range[endage-startage]+startage-41                            

  

  

 



  

End[]

EndPackage[]







