        



(*:Mathematica:: V2.1 *)



(*:Context: "pse`" *)          



(*Title: periodic screeening *)



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

                Baker, S.G. 1998 Biometrics to appear*)



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



(*Requirements:  psederiv           derivatives

                 psefit             fitting the data

                 psepreg            Poisson regresSion

                 psedata            data

                 pseplot            plots of incidence                

                 statplot           plot for confidence intervals

                 matrixlx           matrix programming language

                 psefalse           falsepositives*)



(*History: 1998 Stuart G. Baker *)



(* To reproduce calculations in the manuscript



   PSEHIP[dataHIP, a, z]





To modify the data, change data in  psedata.m

 

  nonrandomized study:    dataOBS

  randomized study:       dataRAN                        



To run PSE for ages 40 to 50 with mortality endpoint use

 

   nonrandomized study:   PSEOBS[dataOBS,a,z,{40,50},"MORT"]

   randomized study:      PSERAN[dataRAN,a,z,{40,50},"MORT"] *)

   



<<psedata.m



BeginPackage["pse`","psederiv`","psefit`","psefalse`",

             "psepreg`","pseplot`","statplot`","matrixlx`"]







PSEHIP::usage="PSE[data,a,z] is the analysis in the manuscript"



PSEOBS::usage="PSEOBS[dataOBS,a,z,model,ageRange,outcome] for the PSE analysis

               of data from a nonrandomized study"



PSERAN::usage="PSERAN[dataRAN,a,z,model,ageRange,outcome] for the PSE analysis

               of data from a randomized study"



PSEMort1::usage="PSEMort1"

PSENode1::usage="PSENode1"







Clear[PSEHIP,PSEOBS,PSERAN]





  

     

Begin["Private`"] 



PSEHIP[data_,a_,z_,model_:"Linear"]:= 

  Module[{resvec},

  resvec=FitAll[data,model,"RAN"];

  PSEPlot[data[[1]],resvec[[1]],a,z,"RAN"];

  PSECore[resvec,a,z];

 Return[Null]]

  

PSEOBS[data_,a_,z_,ageRange_:{40,50}]:=

  PSEOBS[data,a,z,"Linear",ageRange,"MORT"]



PSEOBS[data_,a_,z_,model_,ageRange_,outcome_:"MORT"]:=

  Module[{resvec},

  resvec=FitAll[data,model,"OBS"];

  PSEPlot[data[[1]],resvec[[1]],a,z,"OBS"];

  res={PSE[resvec,a,z,ageRange,outcome,"OBS"]};

  Print@ TableForm[res,

         TableHeadings->{{"NONRAN"},{"est","low","up","pvalue"}}];

 Return[Null]]



PSERAN[data_,a_,z_,ageRange_:{40,50}]:=

  PSERAN[data,a,z,"Linear",ageRange,"MORT"]



PSERAN[data_,a_,z_,model_,ageRange_,outcome_:"MORT"]:=

  Module[{resvec},

  resvec=FitAll[data,model,"RAN"];

  PSEPlot[data[[1]],resvec[[1]],a,z,"RAN"];

  PSESet[resvec,a,z,ageRange,outcome];

 Return[Null]]



  



(*-------------------HIP---------------------------*)



PSECore[resvec_,a_,z_]:=

 Module[{},

  Print[" "];            Print["Ages 40 versus 50"];

  Print["Node"];         PSENode1[resvec,a,z];

  Print["Mortality"];    PSEMort1[resvec,a,z];             

  Print[" "];            Print["Ages 50 versus 60"];

  Print["Node"];         PSENode2[resvec,a,z];

  Print["Mortality"];    PSEMort2[resvec,a,z];

Return[Null]]  

 



PSENode1[resvec_,a_,z_]:=   PSESet[resvec,a,z,{40,50},"NODE"]

PSENode2[resvec_,a_,z_]:=   PSESet[resvec,a,z,{50,60},"NODE"]

PSEMort1[resvec_,a_,z_]:=   PSESet[resvec,a,z,{40,50},"MORT"]

PSEMort2[resvec_,a_,z_]:=   PSESet[resvec,a,z,{50,60},"MORT"]

 



 

PSESet[resvec_,a_,z_,ageRange_,outcome_]:=

 Module[{res},

   res={PSE[resvec,a,z,ageRange,outcome,"OBS"],

        PSE[resvec,a,z,ageRange,outcome,"RAN1"],

        PSE[resvec,a,z,ageRange,outcome,"RAN2"]};

   Print@ TableForm[res,

          TableHeadings->{{"NONRAN","RANI","RANII"},{"est","low","up","pvalue"}}];

  Return[Null]]  

  

  

  



PSE[{detectvec_,mortvec_,nodevec_},a_,z_,{startage_,endage_},outcome_,type_]:=

   Module[{resF,resI,resS,resR,resC,

           pi,gamma,

           hF,hI,hS,hR,hC,hF1,hNRCT,

           nF,nI,nS,nR,nC,

           pF,pI,pS,pR,pC,

           wF,wI,wS,wR,wC,

           ageRange,f,

           qF,qF1,qI,qS,qR,qC,qRCT,           

           qFlist,qF1list,qIlist,qSlist,qNlist,qNRCTlist,

           dqF,dqI,dqS,   

           mF,mF1,mI,mS,mN,mNRCT,

           dhF,dhI,dhS,

           vmF,vmI,vmS,

           vF,vI,vS,

           vqF,vqI,vqS,vqR,vqC,

           g,g0,g1,seg,pvalue},

          

   {resF,resI,resS,resR,resC,pi,gamma}=detectvec;

  {{hF,nF},{hI,nI},{hS,nS},{hR,nR},{hC,nC}}=mortvec; 

  {{pF,wF},{pI,wI},{pS,wS},{pR,wR},{pC,wC}}=nodevec;         

  

 

 (*-------------AGE SPECIFIC INCIDENCE------------------------*)

   ageRange={startage,endage}; 

   f=Exp[z]; 

   qF=polymean[resF,a,z,f];

   qF1=qF/.a->(a+1);

   qI=polymean[resI,a,z,f];

   qS=polymean[resS,a,z,f];

   qN=qF+qI+qS-qF1;

   qFlist=FunctionToList[{qF,a},ageRange];

   qIlist=FunctionToList[{qI,a},ageRange];

   qSlist=FunctionToList[{qS,a},ageRange];

   qF1list=FunctionToList[{qF1,a},ageRange];

   vF=resF[[2]];

   vI=resI[[2]];

   vS=resS[[2]];

 

 (*-----------CANCER MORTATLITY AFTER DETECTION*) 

  If[outcome=="MORT",

    (*mortality probabilities*)

      hF1=Join[Drop[hF,-1],{0}];

      mF= CreateMList[hF,ageRange];

      mF1= Drop[CreateMList[hF,{startage,endage+1}],1];

      mI= CreateMListAvg[hI,ageRange];

      mS= CreateMList[hS,ageRange];

    (*derivatives*)

      dhF=DerivH[qFlist,hF,ageRange]-DerivH[qF1list,hF1,ageRange];

      dhI=DerivHAvg[qIlist,hI,ageRange];

      dhS=DerivH[qSlist,hS,ageRange];

    (*variances*)

      vmF=Apply[Plus,dhF^2 hF (1-hF)/nF];

      vmI=Apply[Plus,dhI^2 hI (1-hI)/nI];

      vmS=Apply[Plus,dhS^2 hS (1-hS)/nS]];

  If[outcome=="NODE",

    (*probabilities*)

      {mF,mF1,mI,mS,mR,mC}={pF,pF,pI,pS,pR,pC};

    (*derivatives*)                                                      

      dhF=qFlist-qF1list;

      dhI=qIlist;

      dhS=qSlist;

    (*variances*)

      vmF=Apply[Plus,dhF^2 pF (1-pF)/wF];

      vmI=Apply[Plus,dhI^2 pI (1-pI)/wI];

      vmS=Apply[Plus,dhS^2 pS (1-pS)/wS]];

       

  

       

   (*ESTIMATES: INCIDENCE*)

     If[type=="OBS" || type=="RAN1",

        qNlist=qFlist+qIlist+qSlist-qF1list]; 

     If[type=="RAN2",

        qC=polymean[resC,a,z,f];

        qR=polymean[resR,a,z,f];

        qNRCT=(qC - (1-pi) qR) / pi;

        vC=resC[[2]];

        vR=resR[[2]];

        qNRCTlist=FunctionToList[{qNRCT,a},ageRange]];

        

     

    (*ESTIMATES: MORTALITY*) 

     If[type=="OBS",  

        If[outcome=="MORT", mN=CreateMListAvg[hR,ageRange]];

        If[outcome=="NODE", mN=mR]];

     If[type=="RAN1" || type=="RAN2", 

        If[outcome=="MORT",

           hNRCT=ReplaceLowerBound[(hC - (1-gamma) hR) / gamma,0,0];

           mNRCT=CreateMListAvg[hNRCT,ageRange]];

        If[outcome=="NODE",    

           mNRCT=(mC - mR (1-gamma))/gamma]]; 

  

     

      (*VARIANCE: AGE-SPECIFIC INCIDENCE*)

       If[type=="OBS",

        dqF=DerivQ[{qF,a},mF-mN,ageRange]-DerivQ[{qF1,a},mF1-mN,ageRange];

        dqI=DerivQ[{qI,a},mI-mN,ageRange];

        dqS=DerivQ[{qS,a},mS-mN,ageRange];

        vqF=dqF . vF . dqF;

        vqI=dqI . vI . dqI;

        vqS=dqS . vS . dqS];

      If[type=="RAN1",

        dqF=DerivQ[{qF,a},mF-mNRCT,ageRange]-DerivQ[{qF1,a},mF1-mNRCT,ageRange];

        dqI=DerivQ[{qI,a},mI-mNRCT,ageRange];

        dqS=DerivQ[{qS,a},mS-mNRCT,ageRange];

        vqF=dqF . vF . dqF;

        vqI=dqI . vI . dqI;

        vqS=dqS . vS . dqS];   

      If[type=="RAN2",

        dqF=DerivQ[{qF,a},mF,ageRange]-DerivQ[{qF1,a},mF1,ageRange];

        dqI=DerivQ[{qI,a},mI,ageRange];

        dqS=DerivQ[{qS,a},mS,ageRange];

        dqC=DerivQ[{qC,a},mNRCT,ageRange]/pi;

        dqR=DerivQ[{qR,a},mNRCT,ageRange](-1)(1-pi)/pi;

        vqF=dqF . vF . dqF;

        vqI=dqI . vI . dqI;

        vqS=dqS . vS . dqS;

        vqC=dqC . vC . dqC;

        vqR=dqR . vR . dqR];

 

      (*VARIANCE CANCER MORTALITY*)

      If[type=="OBS",  

         If[outcome=="MORT", 

           dhR=DerivHAvg[qNlist,hR,ageRange] (-1);          

           vmR=Apply[Plus,dhR^2 hR (1-hR)/nR]];

         If[outcome=="NODE",

           dhR=qNlist (-1);

           vmR=Apply[Plus,dhR^2 pR (1-pR)/wR]]];            

     If[type=="RAN1",     

        If[outcome=="MORT",

           dhC=DerivHAvg[qNlist,hC,ageRange]/gamma;

           dhR=DerivHAvg[qNlist,hR,ageRange](-1)(1-gamma)/gamma;       

           vmC=Apply[Plus,dhC^2 hC (1-hC)/nC];

           vmR=Apply[Plus,dhR^2 hR (1-hR)/nR]];

        If[outcome=="NODE",

           dhC=qNlist /gamma;

           dhR=qNlist (-1)(1-gamma)/gamma;

           vmC=Apply[Plus,dhC^2 pC (1-pC)/wC];

           vmR=Apply[Plus,dhR^2 pR (1-pR)/wR]]];             

     If[type=="RAN2",     

        If[outcome=="MORT",

           dhC=DerivHAvg[qNRCTlist,hC,ageRange]/gamma;

           dhR=DerivHAvg[qNRCTlist,hR,ageRange](-1) (1-gamma)/gamma;       

           vmC=Apply[Plus,dhC^2 hC (1-hC)/nC];

           vmR=Apply[Plus,dhR^2 hR (1-hR)/nR]];

        If[outcome=="NODE",

           dhC=qNRCTlist /gamma;

           dhR=-qNRCTlist (1-gamma)/gamma;

           vmC=Apply[Plus,dhC^2 pC (1-pC)/wC];

           vmR=Apply[Plus,dhR^2 pR (1-pR)/wR]]];     

  

  (*GAIN*) 

   gFlist=qFlist mF;

   gF1list=qF1list mF1;

   gIlist=qIlist mI;

   gSlist=qSlist mS;

   If[type=="OBS", gNlist=qNlist mN];

   If[type=="RAN1",gNlist=qNlist mNRCT];

   If[type=="RAN2",gNlist=qNRCTlist mNRCT];

   g0=gFlist[[1]]+Apply[Plus,gIlist]+Apply[Plus,gSlist];

   g1=Apply[Plus,gNlist]+Take[gF1list,-1][[1]];

   g=g1-g0;

 

   (*check*)

   (*   g0c=Apply[Plus,gFlist]+Apply[Plus,gIlist]+Apply[Plus,gSlist];

        g1c=Apply[Plus,gNlist]+Apply[Plus,gF1list];

        gc=g1c-g0c;

        Print[{g,gc}]; *)

        

    If[type=="OBS", varg=vqF+vqI+vqS         +vmF+vmI+vmS+vmR];

    If[type=="RAN1",varg=vqF+vqI+vqS         +vmF+vmI+vmS+vmR+vmC];

    If[type=="RAN2",varg=vqF+vqI+vqS+vqC+vqR +vmF+vmI+vmS+vmR+vmC];

   

 

   res=report[g,varg];

   seg=Sqrt[varg]//N;

   gL=g-1.96 seg;

   gU=g+1.96 seg;  

   pvalue=PValue[g/seg]; 

   pval=Round[pvalue 100]/100//N;             

   vec=Flatten@Round[{g,gL,gU}10000];

   res=Join[vec,{pval}];

Return[res]]











polymean[{beta_List,var_List?MatrixQ},t_,a_,f_:Exp[a]]:=

   f /. a-> beta .(t^(Range[Length[beta]]-1));

         

         



CDFNormal[mu_,sigma_,x_]:=N[(Erf[(x-mu)/(Sqrt[2] sigma)] + 1) / 2]

CDFNormal01[x_]:=CDFNormal[0,1,x]

PValue[x_]:=1-CDFNormal01[x]   



End[]

EndPackage[]







 

 

detectvec25={{{-6.54155, 0.0521028}, 

              {{0.0914991, -0.00559794}, 

              {-0.00559794, 0.000429387}}}, 

              {{-7.21414, -0.00650757}, 

              {{0.0840737, -0.00562852}, 

              {-0.00562852, 0.000485849}}}, 

              {{-7.47252, 0.059098}, 

              {{0.0865703, -0.00534658}, 

              {-0.00534658, 0.000397777}}}, 

              {{-6.60132, 0.0372269}, 

              {{0.0214488, -0.00130619}, 

              {-0.00130619, 0.0000970896}}}, 

              {{-6.57246, 0.0206399}, 

              {{0.0739044, -0.00458605}, 

              {-0.00458605, 0.000352773}}}, 0.668537,.668537};

                            

                      

mortvec25=FitMort[mortRAN];

nodevec25=FitNode[nodeRAN];



resvec25={detectvec25,mortvec25,nodevec25};



dataHIP=dataRAN;

resvec=

{{{{-6.54155,

        0.0521028},{{0.0914991,-0.00559794},{-0.00559794,

          0.000429387}}},{{-7.09895,-0.00298892},{{

          0.083529,-0.00558149},{-0.00558149,0.00048179}}},{{-7.38556,

        0.0607481},{{0.0864153,-0.00533506},{-0.00533506,

          0.000396919}}},{{-6.57246,

        0.0206399},{{0.0739044,-0.00458605},{-0.00458605,

          0.000352773}}},{{-6.60132,

        0.0372269},{{0.0214488,-0.00130619},{-0.00130619,0.0000970896}}},

    0.668537,0.733875},{{{0.0363636,0.0377358,0,0.0196078,0.0612245,0.0666667,

        0,0.0263158,0,0.0857143,0.0322581,0,0.0384615,0,0},{55,53,51,51,49,45,

        41,38,36,35,31,28,26,25,24}},{{0.0416667,0.152174,0.102564,0.0294118,

        0.0606061,0.1,0,0.08,0.0869565,0.05,0,0.0526316,0,0,0},{48,46,39,34,

        33,30,26,25,23,20,19,19,18,18,16}},{{0.012987,0.0394737,0,0.0136986,

        0.0416667,0.0294118,0.0153846,0.047619,0.0169492,0.0517241,0,0.037037,

        0.0192308,0.0204082,0},{77,76,73,73,72,68,65,63,59,58,54,54,52,49,

        45}},{{0.154676,0.0772532,0.0666667,0.0575916,0.0578035,0.0516129,

        0.0503597,0.0338983,0.00961538,0.045977,0.056338,0,0,0,0},{278,233,

        210,191,173,155,139,118,104,87,71,54,47,41,29}},{{0.0861702,0.0728555,

        0.084399,0.0655271,0.0660377,0.0474453,0.0584551,0.0362319,0.0367232,

        0.0337838,0.0118577,0.0334928,0.0340909,0.0148148,0.030303},{940,851,

        782,702,636,548,479,414,354,296,253,209,176,135,99}}},{{0.265306,49},

        {0.416667,48},{0.22973,74},{0.478261,46},{0.486188,181}}};

