** MACROS that we will use **; 
  %macro logi; 
     %do i=1 %to 20; 
	 data next; set a; ngerm=germ; 
	   if t= &i then ngerm = " "; 
	 proc logistic noprint order=data;
	   model ngerm=soiltemp; 
	   output out=out&i predicted=p; 
	 %end; 
 %mend logi; 
 * MACRO logi removes one response from the dataset, fits a logistic 
 regression to the modified data, and computes the predicted
 probability of germination for all points, including the one with
 missing response, from their soil temperatures. In that sense, it does 
 true "jacknifing."; 


 %macro listout; 
   %do i=1 %to 20; out&i  %end; 
 %mend listout; 
 * Because a macro is just code substitution, the macro listout
 will generate the code out1 out2 out3 etc. and substitute it
 when called; 

************  end of MACROS - program begins ***********; 


data a; input soiltemp n @; 
 do i=1 to n;  input germ $ @; 
   t+1;  * counter for each observation*; 
 output; end; cards; 
 78   1    Y
 68   9    Y N N N N N N N N
 67   9    N Y Y Y Y Y Y Y Y
 50   1    N
  ; 

   %logi; 
data results; set %listout; 
  if ngerm=" " ;   cgerm="N"; if p>.48 then cgerm="Y"; 
proc sort; by p; proc print; title "True jacknifing"; 
proc logistic data=a order=data; 
 model germ = soiltemp/ctable ; 
 output out=full predicted=p predprobs=(x) xbeta=xbeta; 
 proc print; 

 proc freq data=a order=data; 
   table soiltemp*germ/chisq; 

proc catmod data=a order=data; direct soiltemp; 
  model germ=soiltemp; 

  *** graphics ***; 
data b; set a end=eof;
  Y=(germ="Y")+(.015-.03*(germ="Y"))*(i-1); 
   output; 
if eof then do soiltemp=50 to 78; 
 germ=" "; Y=.; output; end; 
 proc logistic data=b order=data;
   model germ=soiltemp/ctable pprob=0.48; 
   output out=out2 predicted=p; 
   data out2; set out2; if Y ne . then P = .; 
   proc print data=out2; 
proc gplot; plot (y p)*soiltemp/overlay vref=0.48 href=66.4; 
  symbol1 v=dot h=1 c=black; 
  symbol2 v=none i=join c=red w=5; 
  title "Logistic on full data"; 
run; 


data bubble; set a; y=(germ="Y"); 
proc means noprint nway; 
  var Y; class soiltemp germ;
  output out=out3 n=nreps; 
data all; set out2 out3;
proc print; run; 
proc gplot;
axis1 offset=(2,2); ** offset to see full circles**; 
bubble germ*soiltemp=nreps/haxis=axis1 vaxis=axis1; 
plot2 p*soiltemp=2;
run;