(* A Mathematica program that relates the PWER and FWER*)
(*  of a confidence band for net health benefits.*)

(* Run CELL 1 once followed by, in any order and as many times as desired *)
(* CELLs 2, 3, 4 or 5.*)
(* The number of programs must be 1 or 2. *)

(*Abstract*)
(*There is considerable interest in the use of net health benefit (NHB) for
decision making in cost-effectiveness analysis. A simultaneous confidence band
for NHB(Lambda) may be obtained by forming a (1-alpha) confidence interval based
on student's t statistic and letting the willingness-to-pay value, Lambda, run
over all values. The pointwise errror rate (PWER) at fixed Lambda of this
confidence interval is alpha. The familywise error rate (FWER) is the
probability that for at least one value of Lambda the simultaneous confidence
band does not cover the true NHB(Lambda). The FWER equals P(T^2> t^2), where
T^2 follows Hotelling's central distribution and t is the (1-Alpha/2) quantile
of the Student t distribution. *)



(* CELL 1 : *)
<< Statistics`MultinormalDistribution`
<< Statistics`NormalDistribution`
Off[General::spell1];
Off[General::spell];

(* CELL 2 : *)
(* This cell computes the FWER given (PWER, total n,and number of programs)*)

   (* Input : enter specific request by replacing sample values*)
      pwer = .085;      (* PWER    *)
      kp = 1;           (* number of programs*)
      totaln = 50;      (* total n *)

   (* Program : *)
      If[kp < 1 || kp > 2, {Print["kp must be 1 or 2, entered kp=", kp], 
        Goto[err]}];
      fwer = 1 - CDF[HotellingTSquareDistribution[2, totaln - kp],
        Quantile[StudentTDistribution[totaln - kp], 1 - pwer/2]^2];
      Print["Input:   PWER=", pwer, "  number of programs=", kp, "  total n=", 
          totaln];
      Print["Output:  FWER=", fwer];
      Label[err];


(* CELL 3 : *)
(* This cell computes PWER given (FWER, total n, and number of programs) *)

   (* Input : enter specific request by replacing sample values*)
      fwer = .230419;        (* FWER    *)
      kp = 1;                (* number of programs*)
      totaln = 50;           (* total n *)

   (* program *)
      If[kp < 1 || kp > 2, {Print["kp must be 1 or 2, entered kp=", kp], 
        Goto[err]}];
      pwer = 2*(1 - CDF[StudentTDistribution[totaln - kp],
            ((2*(totaln - kp)/(totaln - kp - 1))*
            Quantile[FRatioDistribution[2,totaln - kp - 1],
            (1 - fwer)])^.5]);
      Print["Input:   FWER=", fwer, "  number of programs=", kp, "  total n=", 
           totaln];
      Print["Output:  PWER=", pwer];
      Label[err];


(* CELL 4 : *)
(* This cell plots FWER as a function of PWER given (total n,and number of programs)*)

  (* Input : enter specific request by replacing sample values*)
      kp = 2;         (* number of programs *)
      totaln = 200;   (* total n vector *)

  (* Program : *)
      If[kp < 1 || kp > 2, {Print["kp must be 1 or 2, entered kp=", kp], 
        Goto[err]}];
      Print["Input:     number of programs=", kp, "  total n=", totaln];
      Plot[1 - CDF[HotellingTSquareDistribution[2, totaln - kp], 
          Quantile[StudentTDistribution[totaln - kp], 1 - pwer/2]^2],
          {pwer,0, .1},AxesLabel -> {"PWER", "FWER"}];
      Label[err];




(* Cell 5 :    *)
(* This cell computes a table of FWER for each PWER in a range specified by *)
(* user for a vector of total n and a specified number of programs*)

   (* Input : enter specific request by replacing sample values*)
      
      kp = 2;                              (*number of programs *)
      totaln = {10, 50, 100, 200, 1000};   (* total n vector *)
      (* supply lowest value, largest value,and the increment*) 
      (* for the range of the pwer below *)

  (* program *)
     If[kp < 1 || kp > 2, {Print["kp must be 1 or 2, entered kp=", kp], 
        Goto[err]}];
     fwertable = Table[{ToString[PaddedForm[pwer, {4, 3}]],
                Table[{totaln[[n]],ToString[PaddedForm[
                1 - CDF[HotellingTSquareDistribution[2, totaln[[n]] - kp],
                Quantile[StudentTDistribution[totaln[[n]] - kp], 
                1 - pwer/2]^2], {5, 4}]]}, {n, 1, Length[totaln]}]},
     (* --enter below {, lowest value, largest value,increment for pwer--}*)
     {pwer, .005, .10, .005}];
     Print["Input:  number of programs=", kp, "  total n=", totaln];
     Print["Output: Each line is a PWER followed by the pair (n,FWER) in the 
           vector total n"];
     Print[fwertable];  
     Label[err];