(*:Mathematica:: V2.1 *)

(*:Context:: "excens`" *)

(*Title:: informative censoring and double sampling *)

(*References::

             Baker, S.G., Wax, Y. and Patterson, B. 1992 
           Regression analysis of grouped survival data
           in the presence of informative censoring and
           double sampling, Biometrics*)

(*Keywords::  informative censoring *)

(*Requirements::  "clmmat`" *)

(*History::  1992 Stuart G. Baker *)


Clear[exbwp,createx1,createx2bar,createparametername];

y=Transpose[{{180,544,232,6,10,7,1,
              8,87,128,7,15,11,0,
              0,0,3,39,0,3,89,0,30,2,1,1,0,
              0,0,0,2,0,0,10,1,7,3,3,0,0}}];
   
              
mod1={"podds","infrl"};
mod2={"podds","noninf"};             
mod3={"podds","mcar"}; 
mod4={"podds","infall"}; 
mod5={"sat","infall"};
mod={mod1,mod2,mod3,mod4,mod5};



              
exbwp[y_,modelname_]:=
Module[{indexhazbar1,indexhazbar,
        indexhaz1,indexhaz,
        indexprob1,indexprob,
        q1,w1,g1,h1,z1,x1,
         g1bar,g1tilde,
         x1bar,x1tilde,
        func1,x1a,m1v,m1bar,m1,arg1,
        q2,w2,g2,h2,z2,x2,
         e2a,e2b,e2c,e2d,e2s,e2,
         g2tilde,z2tilde,
         x2bar,x2tilde,
         func2,x2a,m2,m2tilde,arg2,
        nv,ncore,n,q,w,x,z,g,h,
        components,mstepfunc,modelspec,
        parametername,rationame,
        c,model,ratio,name},
      
  indexhaz1={1,2,7,8,9,10,13,14,15,16,17,18,19,20,21,22,23,24};
  indexhazbar1={1,4,5,7,8,9,10,11,12};
  indexhazbar=Join[indexhazbar1,indexhazbar1+12];
  indexhaz=Join[indexhaz1,indexhaz1+24];
  indexprob1={1,2,3,4,6,7,8,11,12,13,14,15,16};
  indexprob=Join[indexprob1,indexprob1+16,
                 indexprob1+32,indexprob1+48];
  (*model 1-failure *)
   q1=Identity; 
   w1=J[52,1,0];
     gbar={{1,0,0,0,0,0},
        {0,1,1,0,0,0},
        {0,1,0,1,1,0},
        {0,1,0,1,0,1}}; 
     g1tilde=Dir[J[2,1],Id[2],J[4,1],gbar];
     (* transforms failure harzards for probabilities -- the same for
     each response category, different for drain/no drain, same
     for PF and FF samples*)
   g1=SelectRows[g1tilde,indexprob];
   h1= #2 #1 - Log[1+Exp[#2]]&;  
   z1=Dir[J[6,1],{{1},{0}}];
     x1bar=createx1[modelname];
   x1=x1bar ~Dir~ J[2,1]; 
   func1=NLR; 
    (*failure model arguments for logistic regresion*)
    m1v={{1,0,0,0, 0,0,0,0,0, 1,0,0,0},
         {0,1,0,0, 1,0,0,0,0, 0,1,0,0},
         {0,0,1,0, 0,1,0,1,0, 0,0,1,0},
         {0,0,0,1, 0,0,1,0,1, 0,0,0,1}};
    m1bar=Hcat[Id[2],Id[2]] ~ Dir~ m1v;
    (*collapse over removal categories,PF and FF samples, different
    for drain and no drain*) 
   x1a=x1bar; 
   m1=Dir[Id[2],Transpose[gbar]].m1bar;
      (*create risk sets for drain and no drain using transpose
     of gbar*)
   arg1={m1,x1a};
   
  (*model 2 -- removal*)
   q2=Identity;
   w2=J[52,1,0]; 
     e2a=Dir[Id[4],{{1,0,0,0}}];
     e2b=Dir[Id[4],{{0,1,0,0}}];
     e2c=Dir[Id[4],{{0,0,1,0}}];
     e2d=Dir[Id[4],{{0,0,0,1}}];
     e2s=Vcat[e2a,e2b,e2c,e2d];
     e2=Id[4] ~Dir~ e2s;
        (* transforms the order of probabilities
        (f=1,2,3,4|r=1;f=1,2,3,4|r=2..) to
        (r=1,2,3,4|f=1; r=1,2,3,4|r=2) *)
     g2bar=Dir[J[2,1],Id[8],gbar];
         (* makes probabilities from hazards separately for
         drain and no drain; the same for PF and FF samples*)
      g2tilde=e2.g2bar;   
   g2=SelectColumns[SelectRows[g2tilde,indexprob],indexhaz];
   h2= #2 #1 - Log[1+Exp[#2]]&;
     z2tilde=Dir[J[24,1],{{1},{0}}];
   z2=SelectRows[z2tilde,indexhaz];
    x2bar=createx2bar[modelname];
    x2tilde=x2bar ~Dir~ J[2,1];
   x2=SelectRows[x2tilde,indexhaz];
   func2=NLR;
    (*removal model arguments for logistic regression*)
    x2a=SelectRows[x2bar,indexhazbar];
     m2tilde=Transpose[g2bar].e2;
     (*transforms {f,r} to {r,f} then creates risk sets
     using transposed matrix*)
    m2=SelectRows[SelectColumns[m2tilde,indexprob],indexhaz];
    arg2={m2,x2a};
    
 (*combination*)
  ncore=BlockDiag[J[1,7],J[1,7],J[1,13],J[1,13]] . y;
  nv=ncore ~Dir~ J[13,1];
  n={nv,ncore};
  q={q1,q2};
  w={w1,w2}; 
  g={g1,g2};
  h={h1,h2};
  z={z1,z2};
  x={x1,x2};
  
   cpf=BlockDiag[J[1,4],J[1,3],J[1,2],Id[4]];
 c=BlockDiag[cpf,cpf,Id[26]];
   components={n,q,w,g,h,z,x};
   mstepfunc={func1,func2}; 
    modelspec={arg1,arg2};
  model={components,mstepfunc,modelspec};
  ratio="None";
    parametername=createparametername[modelname];
    rationame=Automatic;
  name={parametername,rationame};
  
Return[{c,model,ratio,name}]]/;MemberQ[mod,modelname]



createx1[modelname_]:=
 Module[{x1bar0,x1bart,x1bard,x1bar},
    x1bar0=J[6,1];
    x1bart=Dir[J[2,1],{{0,0},{1,0},{0,1}}];
    x1time=Dir[J[2,1],Id[3]];
    x1bard=Dir[{{0},{1}},J[3,1]];
    x1bar=Switch[modelname[[1]],
    "sat",   Id[6],
    "podds", Hcat[x1time,x1bard]];
 Return[x1bar]]



createx2bar[modelname_]:=
 Module[{x2bar0,x2bart,x2bard,x2barr1,x2bar},
    x2bar0=J[24,1];
    x2bart=Dir[J[8,1],{{0,0},{1,0},{0,1}}];
    x2bard=Dir[{{0},{1}},J[12,1]];
    x2fr=J[2,1] ~Dir~ Transpose[{{1,1,1, 1,1,1, 1,1,1, 0,0,0}}];
    x2barrl=J[2,1] ~Dir~ Transpose[{{1,0,0, 0,1,0, 0,0,1, 0,0,0}}];
    x2barr2=J[2,1] ~Dir~ Transpose[{{0,1,1, 1,0,1, 1,1,0, 1,1,1}}];
    x2all=Id[2]~Dir~ {{1,  0,0, 0,0,0, 0,0,0},
                       {0,  0,0, 0,0,0, 0,0,0},
                       {0,  0,0, 0,0,0, 0,0,0},
                       
                       {0,  1,0, 0,0,0, 0,0,0},
                       {0,  0,1, 0,0,0, 0,0,0},
                       {0,  0,0, 0,0,0, 0,0,0},
                        
                       {0,  0,0, 1,0,0, 0,0,0},
                       {0,  0,0, 0,1,0, 0,0,0},
                       {0,  0,0, 0,0,1, 0,0,0},
                       
                       {0,  0,0, 0,0,0, 1,0,0},
                       {0,  0,0, 0,0,0, 0,1,0},
                       {0,  0,0, 0,0,0, 0,0,1}};
     x2bar=Switch[modelname[[2]],
             "infrl",   Hcat[x2bar0,x2barrl,x2bart,x2bard],
             "fr",      Hcat[x2bar0,x2fr,x2bart,x2bard],
             "mcar",    Hcat[x2bar0],
             "infall",  Hcat[x2all],
             "noninf",  Hcat[x2bar0,x2bart,x2bard]];
  Return[x2bar]]  

createparametername[modelname_]:=
    Switch[modelname[[2]],
       "infrl",   {"alpha0","alphat2","alphat3", "alphad1",
                   "beta0","betaRL","betat2","betat3","betad1"},

       "fr",      {"alpha0","alphat2","alphat3", "alphad1",
                   "beta0","betaFR","betat2","betat3","betad1"},
       "mcar",      {"alpha0","alphat2","alphat3", "alphad1",
                    "beta0"},      
        "infall",   {"alpha0","alphat2","alphat3", "alphad1",
                     "a11","a12","a22","a13","a23","a33","a14","a24","a34",
                     "b11","b12","b22","b13","b23","b33","b14","b24","b34"},                  
       "noninf",   {"alpha0","alphat2","alphat3", "alphad1",
                    "beta0","betat2","betat3","betad1"}]
     
     


        


