* version 1.0.0 05dec2000 MER+JanBrogger modified from sskdlg STB
capture program drop sskapp
program define sskapp , rclass
	version 6.0

	**The lines enclosed in comments are the ones that are different
	**from the dialog box.
	**Most of what they do is set the global macros
	**The rest is done by the original script
	**Some parts of the original have been removed - that have to do with
	**the dialog box and graphing

	** Two lines are added at the bottom to make it an r-class program
	** Most of the capture clear statements have been removed and replace
	** by preserve/restore at the beginning and end
	** A drop _all statement has been added in the middle of the program

	**Some of the display lines have _n newlines eliminated from them
	**to shorten output

	syntax , p1(real) p2(real) diff(real) kapp(real) [n(integer 0) level(integer 95)]

	preserve

	global D_kappa=`kapp'
	global D_p1=`p1'
	global D_p2=`p2'
	if `n'>0 {
		global D_cb8=1
		global D_ssize=`n'
	}
	global D_level=`level'

	**End of command-line additions (two more lines at bottom)
	
	global z = invnorm(1-(1-$D_level/100)/2)

*  Exit situations + Warning messages 

	if $D_kappa < -1 {
		di _n in re "  "  _n
		di    in re "  The specified value for the minimum kappa worth"  _n
   		di    in re "  finding is less than -1, which is meaningless."  _n
   		di    in re "  Sample size will not be calculated." _n 
		di _n in re "  "  _n
		*capture clear
	   	exit
 	}

 	if $D_kappa > 1 {
		di _n in re "  "  _n
		di    in re "  The specified value for the minimum kappa"  _n
   		di    in re "  worth finding is above 1, which is meaningless."  _n
   		di    in re "  Sample size will not be calculated." _n 
		di _n in re "  "  _n
		*capture clear
	   	exit
 	}

 	if $D_p1 <= 0 {
		di _n in re "  "  _n
		di    in re "  The specified value for p1 is less or equal to 0,"  _n
   		di    in re "  which is meaningless. Sample size will not be"  _n
   		di    in re "  calculated." _n 
		di _n in re "  "  _n
		*capture clear
	   	exit
 	}

 	if $D_p1 > 1 {
		di _n in re "  "  _n
		di    in re "  The specified value for p1 is above 1,"  _n
   		di    in re "  which is meaningless. Sample size will not be"  _n
   		di    in re "  calculated." _n 
		di _n in re "  "  _n
		*capture clear
	   	exit
 	}

 	if $D_p2 <= 0 {
		di _n in re "  "  _n
		di    in re "  The specified value for p1 is less or equal to 0,"  _n
   		di    in re "  which is meaningless. Sample size will not be"  _n
   		di    in re "  calculated." _n 
		di _n in re "  "  _n
		*capture clear
   	exit
 	}

 	if $D_p2 > 1 {
		di _n in re "  "  _n
		di    in re "  The specified value for p1 is above 1,"  _n
   		di    in re "  which is meaningless. Sample size will not be"  _n
   		di    in re "  calculated." _n 
		di _n in re "  "  _n
		*capture clear
	   	exit
 	}

*  Evaluating maximum kappa given a set o specified marginals 

	if "$D_p2" == "" {$D_p1 == p} 

  	if $D_p1 >= $D_p2 { 
     		global kmax = (($D_p2+(1-$D_p1))-(($D_p1*$D_p2)+                  /*
*/           		     ((1-$D_p1)*(1-$D_p2))))/(1-(($D_p1*$D_p2)+           /*
*/           		     ((1-$D_p1)*(1-$D_p2)))) } 
  	if $D_p1 < $D_p2 {
     		global kmax = (($D_p1+(1 -$D_p2))-(($D_p1*$D_p2)+                 /*
*/                		((1-$D_p1)*(1-$D_p2))))/(1-(($D_p1*$D_p2)+        /*
*/                		((1-$D_p1)*(1-$D_p2)))) }

* Evaluating minimum kappa given a set o specified marginals 

  	if $D_p1 >= $D_p2 { 
     		global kmin = (((1-$D_p2)-$D_p1)-(($D_p1*$D_p2)+                 /*
*/                		((1-$D_p1)*(1-$D_p2))))/(1-(($D_p1*$D_p2)+       /*
*/                		((1-$D_p1)*(1-$D_p2)))) }
  	if $D_p1 < $D_p2 {
     		global kmin = (((1-$D_p1)-$D_p2)-(($D_p1*$D_p2)+                 /*
*/                		((1-$D_p1)*(1-$D_p2))))/(1-(($D_p1*$D_p2)+       /*
*/                		((1-$D_p1)*(1-$D_p2)))) }

* More global macros 

  	global rkmax = round($kmax, .01)
  	global kplusd = ($D_kappa+$D_d)
  	global kmind = ($D_kappa-$D_d)
  	global rkmin = round($kmin, .01)

* Calculating Q, conditional on a permissible maximum value of kappa  
*     (given specified marginals) 

 	if $kmax < $D_kappa {
		di _n in re "  "  _n
		di    in re "  Marginals p1 = $D_p1 and p2 = $D_p2 are incompatible"  _n
   		di    in re "  with the specified kappa ($D_kappa), since the "  _n
   		di    in re "  maximum possible value is $rkmax."  _n
   		di    in re "  Sample size will not be calculated." _n 
		di _n in re "  "  _n
		*capture clear
	   	exit
 	}

 	if $D_kappa < $kmin {
		di _n in re "  "  _n
		di    in re "  Marginals p1 = $D_p1 and p2 = $D_p2 are incompatible"  _n
   		di    in re "  with the specified kappa ($D_kappa), since the "  _n
   		di    in re "  minimum possible value is $rkmin. " _n
   		di    in re "  Sample size will not be calculated." _n 
		di _n in re "  "  _n
		*capture clear
	   	exit
 	}

  	if $kmax < $kplusd & $kmin > $kmind {
		di _n in re "  "  _n
		di    in re "  The specified value of d ($D_d) is incompatible"  _n  
   		di    in re "  with the kappa you selected ($D_kappa), since"  _n 
   		di    in re "  both confidence limits ($kmind and $kplusd)"  _n 
                di    in re "  exceed the possible boundary values for kappa "  _n 
   		di    in re "  ($rkmin and $rkmax, respectively), given the  "   _n 
   		di    in re "  specified marginals (p1 = $D_p1 and p2 = $D_p2). "  _n 
                di    in re "  Sample size will not be calculated." _n 
		di _n in re "  "  _n
		*capture clear
	   	exit
 	}

 		global Q = ( ( ( (-1 + $D_kappa) * ( (-2*$D_kappa*$D_p1) +                /*    
*/		  	    (($D_kappa^2)*$D_p1) + (4*$D_kappa*($D_p1^2)) -               /*
*/		  	    (2*($D_kappa^2)*($D_p1^2)) - (2*$D_kappa*$D_p2) +             /*
*/		  	    (($D_kappa^2)*$D_p2) - (4*$D_p1*$D_p2) +                      /*
*/		  	    (8*$D_kappa*$D_p1*$D_p2) - (6*($D_kappa^2)*$D_p1*$D_p2) +     /*
*/		  	    (4*($D_p1^2)*$D_p2) - (12*$D_kappa*($D_p1^2)*$D_p2) +         /*
*/		  	    (8*($D_kappa^2)*($D_p1^2)*$D_p2) + (4*$D_kappa*($D_p2^2)) -   /*
*/		  	    (2*($D_kappa^2)*($D_p2^2)) + (4*$D_p1*($D_p2^2)) -            /*
*/		  	    (12*$D_kappa*$D_p1*($D_p2^2)) +                               /*
*/		  	    (8*($D_kappa^2)*$D_p1*($D_p2^2)) -                            /*
*/		  	    (4*($D_p1^2)*($D_p2^2)) +                                     /*
*/		  	    (12*$D_kappa*($D_p1^2)*($D_p2^2)) -                           /*
*/		  	    (8*($D_kappa^2)*($D_p1^2)*($D_p2^2)))                         /*
*/		  	    ) /  ( ($D_p1 + $D_p2 - (2*$D_p1*$D_p2) )^2 )   )  )

  	global rQ = round($Q,.001)

* Calculating sample size 

  	global S = ($Q*($z/$D_d)^2)
  	global rS = round($S,1)

* Calculating dout from sample(S) 

	if "$D_cb8" == "1" { 
		global dout = sqrt(($Q * $z) / $D_ssize)
		global rdout = round($dout, .001)
	}

  	if $kmax < $kplusd {
   		di _n in bl "  Note: the specified value of d ($D_d) is "  
   		di    in bl "  partially incompatible with the kappa you" 
   		di    in bl "  selected ($D_kappa), since the upper confidence "
                di    in bl "  limit ($kplusd) would be beyond the maximum " 
   		di    in bl "  possible value for kappa ($rkmax) for the  "  
   		di    in bl "  specified marginals (p1 = $D_p1 and p2 = $D_p2). "
                di    in bl "  Yet, sample size is calculated, since d is "
		di    in bl "  compatible with kappa regarding the lower"
                di    in bl "  limit." _n 
 	}

 	if $kmin > $kmind {
   		di _n in bl "  Note: the specified value of d ($D_d) is "  
   		di    in bl "  partially incompatible with the kappa you" 
   		di    in bl "  selected ($D_kappa), since the lower confidence "
                di    in bl "  limit ($kmind) would be below the minimum  " 
   		di    in bl "  possible value for kappa ($rkmin) for the  "  
   		di    in bl "  specified marginals (p1 = $D_p1 and p2 = $D_p2). "
                di    in bl "  Yet, sample size is calculated, since d is "
		di    in bl "  compatible with kappa regarding the upper"
                di    in bl "  limit." _n 
 	}

  	if $D_kappa < 0 & $D_kappa > -1 {
   		di _n in bl "  Atention! The specified values for the"  
   		di    in bl "  minimum kappa worth finding is less than 0." 
   		di    in bl "  Although negative values down to -1 are" 
   		di    in bl "  possible, it is very unusual to specify a"
                di    in bl "  value less than 0." _n 
 	}


* ===========  Stuff for qmax, graph(q) and graph(s) ============ 

  	*capture clear
	drop _all
  	qui set obs 1000
  	global k  (((_n/1000)/.5)-1)
  	gen k = $k

* Putting globals into vars

  	gen p1 = $D_p1
  	gen p2 = $D_p2
  	gen d = $D_d
  	gen z = $z
  	gen kmax = $kmax
  	gen kmin = $kmin
  	gen kplusd = $kplusd
  	gen kmind = $kmind

* Generating Q. When incompatible situations occur, Q is set to missing
*    and plots for those will be absent. If all are missing, there will
*    be no plot at all (only an empty frame)

      	gen Q = ( ( ( (-1 + k) * ( (-2*k*p1) +                   /*    
*/		    ((k^2)*p1) + (4*k*(p1^2)) -                  /*
*/		    (2*(k^2)*(p1^2)) - (2*k*p2) +                /*
*/		    ((k^2)*p2) - (4*p1*p2) +                     /*
*/		    (8*k*p1*p2) - (6*(k^2)*p1*p2) +              /*
*/		    (4*(p1^2)*p2) - (12*k*(p1^2)*p2) +           /*
*/		    (8*(k^2)*(p1^2)*p2) + (4*k*(p2^2)) -         /*
*/		    (2*(k^2)*(p2^2)) + (4*p1*(p2^2)) -           /*
*/		    (12*k*p1*(p2^2)) +                           /*
*/		    (8*(k^2)*p1*(p2^2)) -                        /*
*/		    (4*(p1^2)*(p2^2)) +                          /*
*/		    (12*k*(p1^2)*(p2^2)) -                       /*
*/		    (8*(k^2)*(p1^2)*(p2^2)))                     /*
*/		    ) /  ( (p1 + p2 - (2*p1*p2) )^2 )   )  )

   	qui replace Q=. if kmax < k
   	qui replace Q=. if k < kmin
   	qui replace Q=. if kmax < kplusd & kmin > kmind 

   	qui egen qmax = max(Q)
   	qui gen rqmax = round(qmax, .001)
	global rqmax rqmax

   	qui gen Smax = (rqmax*(z/d)^2)
   	qui gen rSmax = round(Smax,1)
	global rSmax rSmax
 
   	qui gen S = (Q*(z/d)^2)
   	qui gen rS = round(S,1)

   	qui gen kqmax=.
   	qui replace kqmax=k if Q==qmax & Q~=.
   	qui gen rkqmax=round(kqmax,.001)
	global rkqmax rkqmax

   	qui sort kqmax
   	qui egen Qbar=mean(Q)


* ============ OUTPUTS (sample size, Q, Qmax and dout) ======================

	set more 1

	di _n in gr "             " _n
	di  in gr " Sample size for kappa" _n

   	di  in wh "   Results for kappa=$D_kappa, p1=$D_p1, p2=$D_p2," 
   	di  in wh "   d=$D_d [$D_level% Conf. Interval]:" 
   	di  in gr "    *  Sample size = " in ye $rS 

	if "$D_cb1" == "1" {
   		di _n in gr "    *  Value of Q = " in ye $rQ _n 
	}

 	if "$D_cb2" == "1" {
   		if Qbar~=. {
   			di  in wh "  " _n
  			di  in wh "   Given the values specified above, the" 
    			di  in wh "   maximum sample size is " in ye $rSmax in wh " for a" 
    			di  in wh "   kappa of " in ye $rkqmax _n

  		}
  		if Qbar==. {
 
   			di _n in re "  The specified values preclude calculating any Q." 
   			di    in re "  " _n
   			di _n in re "  Sample size cannot be calculated." _n
  		exit
  		}
	}

	if "$D_cb8" == "1" { 
   		di  in gr "  " 
   		di  in wh "   If the sample size is fixed at " in gr $D_ssize in wh ", " 
   		di  in wh "   given the kappa, p1, p2 and CI stated " 
   		di  in wh "   above," in gr " d = " in ye $rdout  in wh "." _n
	}
   	*di _n in gr " --- End -------------------------------------------------" _n
	*di _n in gr "             " _n

	sort k


* =========== Stuff for graph(p) ============================= 

  	global newp1 _n/1000
  	qui gen newp1 = $newp1
        qui gen newp2=newp1	
	qui gen kap=$D_kappa

* Evaluating maximum kappa given the marginals for w

     	qui gen newkmax = ((newp2+(1-newp1))-((newp1*newp2)+                            /*
*/            	          ((1-newp1)*(1-newp2))))/(1-((newp1*newp2)+                    /*
*/           	          ((1-newp1)*(1-newp2))))  

        qui replace newkmax=. if newkmax > -1

* Evaluating minimum kappa given a set o specified marginals

     	qui gen newkmin = (((1-newp2)-newp1)-((newp1*newp2)+                            /*
*/                         ((1-newp1)*(1-newp2))))/(1-((newp1*newp2)+                   /*
*/                         ((1-newp1)*(1-newp2)))) 

        qui replace newkmin=. if newkmin < -1

  	qui gen nkplusd = kap+d
  	qui gen nkmind = kap-d

* Generating Q for w. When incompatible situations occur, Q is set to missing
* and plots for those will be absent. If all are missing, there will
* be no plot at all (only an empty frame)

      	qui gen newQ = ( ( ( (-1 + kap) * ( (-2*kap*newp1) +                           /*    
*/		       ((kap^2)*newp1) + (4*kap*(newp1^2)) -                           /*
*/		       (2*(kap^2)*(newp1^2)) - (2*kap*newp2) +                         /*
*/		       ((kap^2)*newp2) - (4*newp1*newp2) +                             /*
*/		       (8*kap*newp1*newp2) - (6*(kap^2)*newp1*newp2) +                 /*
*/		       (4*(newp1^2)*newp2) - (12*kap*(newp1^2)*newp2) +                /*
*/		       (8*(kap^2)*(newp1^2)*newp2) + (4*kap*(newp2^2)) -               /*
*/		       (2*(kap^2)*(newp2^2)) + (4*newp1*(newp2^2)) -                   /*
*/		       (12*kap*newp1*(newp2^2)) +                                      /*
*/		       (8*(kap^2)*newp1*(newp2^2)) -                                   /*
*/		       (4*(newp1^2)*(newp2^2)) +                                       /*
*/		       (12*kap*(newp1^2)*(newp2^2)) -                                  /*
*/		       (8*(kap^2)*(newp1^2)*(newp2^2)))                                /*
*/		       ) /  ( (newp1 + newp2 - (2*newp1*newp2) )^2 )   )  )

   	qui replace newQ=. if newkmax < kap
   	qui replace newQ=. if kap < newkmin
   	qui replace newQ=. if newkmax < nkplusd
   	qui replace newQ=. if newkmin > nkmind 

   	qui egen nqmax = max(newQ)
   	qui gen rnqmax = round(nqmax, .001)

   	qui gen nSmax = (rnqmax*(z/d)^2)
   	qui gen rnSmax = round(nSmax,1)
 
   	qui gen nS = (newQ*(z/d)^2)
   	qui gen rnS = round(nS,1)

   	qui gen nkqmax=.
   	qui replace nkqmax=kap if newQ==nqmax & newQ~=.
   	qui gen rnkqmax=round(nkqmax,.001)

   	qui sort nkqmax
   	qui egen nQbar=mean(newQ)

	sort newp1


* =========== Stuff for graph(d) ============================= 

	global mult = (($D_tod)/1000)
  	qui gen srange = $mult * _n
	qui gen Qspec = $Q

	qui gen drange = sqrt((Qspec * z) / srange)
	qui gen kplusdr = kap+drange
	qui gen kmindr = kap-drange
	qui replace drange=. if kmax < kplusdr & kmin > kmindr  


* =========== double save globals ============================

   	qui sum qmax
   	global qmax `r(mean)'
   	qui sum Smax
   	global Smax `r(mean)'
   	qui sum kqmax
   	global kqmax `r(mean)'
   	qui sum Qbar 
   	global qbar `r(mean)'  

   	global S_1  = $S
   	global S_2  = $Smax 
   	global S_3  = $Q 
   	global S_4  = $qmax 
   	global S_5  = $D_kappa
    	global S_6  = $kmax 
   	global S_7  = $kmin 
   	global S_8  = $D_p1
   	global S_9  = $D_p2
   	global S_10 = $D_d
   	global S_11 = $z 


	**These two lines added to form command-line wrapper
	return scalar dfixn=$rdout
	if "$D_cb8"=="1" {
		return scalar nfixd=$rS 
	}
	*Restore the data
	restore
	**End of additions
end