



(*:Mathematica:: V2.1 *)

               

(*:Context: "pseplot`" *)          



(*Title: periodic screeening *)



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



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



(*Requirements:  matrixlx *)



(*History: 1998 Stuart G. Baker *)



 <<psedata.m

 

BeginPackage["pseplot`","statplot`","matrixlx`"]





PSEPlot::usage="PSEPlot[detectdata,detectvec,a,z,type]"





Clear[PSEPlot]





Begin["Private`"] 





PSEPlot[detectdata_,detectvec_,a_,z_,type_:"RAN"]:=

  Module[{font,

           resf,resi,ress,resc,resr,pi,pix,

           qf,qf1,qI,qS,qfdel,qa,

           varI,varS,varFdel,va,

           qc,qr,qb,

           vc,vr,vb,

           svec,

           namefx,nameix,namesx,

           sf,sf1,si,ss,sr,

           ff,ff1,fi,fs,fr,

           DIF,vDIF,seDIF,

           plotF,plotI,plotS,plotDIF,plotN,plotNrct,

           ticksx,ticksy},

    {resf,resi,ress,resr,resc,pi,pix}=detectvec;

    display="None";

    font={"Courier",10};

    f=Exp[z];

    qF=PolyMean[resf,a,z,f];

    qI=PolyMean[resi,a,z,f];

    qS=PolyMean[ress,a,z,f];

    qFdel=PolyMeanDel[resf,a,z,f];

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

    qN=qI+qS-qFdel;

  

    varF=PolyVar[resf,a,z,f];

    varI=PolyVar[resi,a,z,f];

    varS=PolyVar[ress,a,z,f];

    varFdel=PolyVarDel[resf,a,z,f];

    varN=varI+varS+varFdel;

      

     

   age=Range[25]-1;

   ticksx={{.5,"40"},{10,"50"},{20,"60"}};

   size=2;

   basic={age,size,ticksx,font,display};

   {namefx,nameix,namesx,namerx,namecx}=

         {"q(F|a) ","q(I|a) ", "q(S|a)", "q(R|a)", "q(C|a)"};

  {{qFx,nfx},{qIx,nix},{qSx,nsx}}=FractionFIS[detectdata];

   plotF=PointPlot[basic,{nfx,qFx,qF,varF,namefx},a];

   plotI=PointPlot[basic,{nix,qIx,qI,varI,nameix},a];

   plotS=PointPlot[basic,{nsx,qSx,qS,varS,namesx},a];

 

   seN=Sqrt[varN];

   ticksy={{.002,".002"},{.004,".004"},{.006,".006"}};

   plotN=PlotCI[qN,seN,{a,0,24},

         PlotName->"    q(N|a)    ",

         TicksX->ticksx,

         TicksY->ticksy,

         PlotRangeY->{0,.006},

         AxesLabelX->"a",

         Font->font,

        ShowPlot->False];

 

  

  If[type=="RAN",  

    qC=PolyMean[resc,a,z,f];

    qR=PolyMean[resr,a,z,f];

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

    varC=PolyVar[resc,a,z,f];

    varR=PolyVar[resr,a,z,f];

    varNrct=(varC + (1-pi)^2 varR)/pi^2;

    seNrct=Sqrt[varNrct];

    plotNrct=PlotCI[qNrct,seNrct,{a,0,24},

         PlotName->"    q*(N|a)   ",

         TicksX->ticksx,

         TicksY->ticksy,

         PlotRangeY->{0,.006},

         AxesLabelX->"a",

         Font->font,

         ShowPlot->False];

    DIF=qN-qNrct;

    vDIF=varN+varNrct;

    seDIF=Sqrt[vDIF];

    ticksx={{.5,"40"},{20,"60"}};

    ticksyDIF={{-.002,"-.002"},{0,"0"},{.002,".002"}};

    plotDIF=PlotCI[DIF,seDIF,{a,0,24},PlotName->"q(N|a)-q*(N|a)",

         TicksX->ticksx,TicksY->ticksyDIF,AxesLabelX->"a",Font->font,

         PlotRangeY->{-.006,.006},ShowPlot->False];

     Show[GraphicsArray[{{plotF,plotN},{plotI,plotNrct},{plotS,plotDIF}}],

                      DisplayFunction->$DisplayFunction],

 (*else nonrandomized*)

    Show[GraphicsArray[{{plotF,plotI},{plotS,plotN}}],    

                      DisplayFunction->$DisplayFunction]          

                      ];

                      

                      



   

Return[Null]]  

 

  

    

PointPlot[{age_,size_,ticksx_,font_,display_},{nx_,qx_,q_,v_,namex_},a_]:=

  Module[{mx,plot1,plot2,plot,ticksy},

   mx=Flatten[(nx^.25) / (Min[nx]^.25)];

      plot1=ListPlotSize[Transpose[{age,qx,mx size}],Identity];

   ticksy={{.002,".002"},{.004,".004"},{.006,".006"}};

   plot2=PlotCI[q,Sqrt[v],{a,0,24},PlotName->namex,

      TicksX->ticksx,TicksY->ticksy,AxesLabelX->"a",PlotRangeY->{0,.008},

      Font->font,ShowPlot->False,Prolog->PointSize[.1]];

   plot=Show[plot2,plot1,DisplayFunction->Identity];    

 Return[plot]]

 







FractionFIS[{dataF_,dataI_,dataS_,datanumFIS_,datacensFIS_,

               dataR_,dataC_,datanumR_,datanumC_}]:=

  Module[{resf,resi,ress,resc,resr},

   resf={dataF/datanumFIS//N,datanumFIS};

   resi=RateCens[dataI,datanumFIS,datacensFIS];

   ress=RateCens[dataS,datanumFIS,datacensFIS];

   resdata={resf,resi,ress};

 Return[resdata]]

  

  

Rate[x_,m_]:= 

Module[{r,c,xadd,n,nadd,qf},

  {r,c}=Dimensions[x];

  xadd=ShiftAdd[x];

  n=Transpose[{m}] ~Dir~ J[1,c];

  nadd=ShiftAdd[n];

  qadd=Flatten[xadd/nadd]//N;

 Return[{qadd,nadd}]]



RateCens[x_,datanumFIS_,datacensFIS_]:= 

  Module[{nc,n,xadd,nadd,q},

   n0=datanumFIS;

   cr=datacensFIS;    

  {cr1,cr2,cr3}=Transpose[cr];

  {r,c}=Dimensions[x];

   n=Switch[c,

    1, n0,

    2, {n0,n0-cr1},

    3, {n0,n0-cr1,n0-cr1-cr2},

    4, {n0,n0-cr1,n0-cr1-cr2,n0-cr1-cr2-cr3}];

   n1=Transpose[n];

   xadd=ShiftAdd[x];

   nadd=ShiftAdd[n1];

   qadd=Flatten[xadd/nadd]//N;

Return[{qadd,nadd}]]  



ShiftAdd[matrix_]:=

Module[{a,b,mat,vec},

{a,b}=Dimensions[matrix];

 mat=Table[shiftadd[matrix,i],{i,1,b}];

 vec=Apply[Plus,mat];

Return[vec]]              



shiftadd[matrix_,i_]:=

Module[{col,col1,col2},

 col=SelectColumns[matrix,{i}];

 If[i>1,

 col1=RemoveRows[col,Range[i-1]] ~Vcat~ J[i-1,1,0],

 col1=col];

Return[Flatten[col1]]]





(*----------------Mean and Variance from Poisson Regression Output-----*)





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

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



    

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

   Module[{p,derivf,deriv},

     p=Length[beta];

     powervec=Range[p]-1;

     derivf=D[f,a] /. a->(beta .(t^powervec));

     deriv=derivf (t^powervec);   

     varl=deriv . var. deriv//Simplify;

    Return[varl]]





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

  Module[{eqn,eqn1,del},

    eqn=PolyMean[{beta,var},t,a,f];

    eqn1=eqn/.(t->(t+1));

    del=eqn1-eqn//Simplify;

  Return[del]]



  

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

   Module[{p,derivf,deriv0,deriv1,deriv},

     p=Length[beta];

     powervec=Range[p]-1;

     derivf=(t^powervec) (D[f,a] /. a->(beta . (t^powervec)));

     derivf1=derivf/.(t->t+1);

     deriv=derivf -derivf1;

     varl=deriv . var. deriv//Simplify;

    Return[varl]]   

   









End[]

EndPackage[]





