/* Program for plotting the power curve of an F-test
   of the null hypothesis of equal treatment means
   against the sample size for a one-way ANOVA.
   Here  n  represents the number of observations
   for each treatment.                             */


proc iml;

    start powerf;

    alpha = .05;         /* Enter the Type I error level  */

    k = 3;               /* Enter the number of treatments  */

    n1 = 2;              /* Enter the lower bound on the sample size */
    n2 = 50;             /* Enter the upper bound on the sample size */
    ninc = 2;            /* Enter the sample size increment */

    nstep = int((n2-n1)/ninc) + 1;
    power = j(nstep,1,0);
    nobs = j(nstep,1,0);

    do  i= 1 to nstep;
      n = n1 + ninc*(i-1);   /*  compute the sample size */
      a=k-1; b=k*(n-1);      /*  Set the degrees of freedom  */

      d=n/2;                 /*  Evaluate the noncentrality parameter */
                        
                             /*  Evaluate power  */
      f = finv(1-alpha,a,b);
      power[i,1] =  1 - probf(f,a,b,d);
      nobs[i,1] = n;
      end;

  power = power||nobs;
  create out1 from power;
  append from power;

  finish;

run powerf;



goptions  cback=white colors=black device=WIN target=WINPRTC;

  axis1 label = (h=3 r=0 a=90 f=swiss 'Power')
        length = 5 in
        color = black width=8.0 style=1;

  axis2 label = (h=2.5 f=swiss 'Sample size')
        length = 5.5 in
        color=black width=8.0 style=1;

  symbol1 v=none i=spline l=1 h=2 c=black width=8;

proc gplot data=out1;
  plot col1*col2 / vaxis=axis1 haxis=axis2;
  title h=3 f=swiss  'Power versus Sample Size';
  title2 h=2 f=swiss  'F-test for equal means';
  run;


/* Program for computing the smallest sample size
   that produces a desired power value of an F-test
   of a null hypothesis of the form CB=0.  */
  
proc iml;

    start powerfun;

    alpha = .05;         /* Enter the Type I error level  */
	power = .90;         /* Enter desired power */

                         /* Enter C matrix for a null 
	                        hypothesis of the form CB=0 */

     C = { 0 1 0 -1,
	       0 0 1 -1};
						 /* Enter deviations as a multiple of
							the error standard deviation  */
	 d= { .5, 1};

	                     /* Enter the rows of the model matrix that
	                        be replicated as the sample size increases */

    W={1 1 0 0,
	   1 0 1 0,
	   1 0 0 1};
	   					 /* Compute the rank of W */
	PW = W*GINV(T(w)*W)*T(w);
    call svd(U,Q,V,W);
	rankw = sum(q);
     
 
    n1 = 2;              /* Enter the lower bound on initial sample size*/
    npower=n1;

	X = W;               /* Set up initial model matrix  */
    do i = 2 to n1;
	  X = X//W;
	  end;

    n2 = 100;            /* Enter the upper bound on the sample size */
    ninc = 1;            /* Enter the sample size increment */

    nstep = int((n2-n1)/ninc) + 1;
   

    do  i= 1 to nstep;
      n = n1 + ninc*(i);   /*  compute the sample size */

	  do j = 1 to ninc;      /*  Expand the model matrix */
	  X=X//W;
	  end;

      a=nrow(C);            /*  Set the degrees of freedom  */
      b=nrow(X)-rankw;

                             /*  Evaluate the noncentrality parameter */
	  nc=t(d)*inv(C*ginv(t(X)*X)*t(C))*d;

                             /*  Evaluate power  */

      f = finv(1-alpha,a,b);
      if(1 - probf(f,a,b,nc)< power) then npower=n;
      end;

  npower=npower+1;
  if(npower< n1+ninc*nstep+1) then do;
  print "sample size; " npower;
  stop;
  end;

  print "Equired sample size is larger than " npower; 

  finish;

run powerfun;