Stat Computing > SAS > Textbook Examples > Intro to Generalized Linear Models by Dobson

SAS Textbook Examples
An Introduction to Generalized Linear Models by Annette J. Dobson
Chapter 1:  Introduction

1.6.5 Example: Tropical cyclones Table 1.2 on page 13.

data table1_2;
  input years $1-7	season	number;
datalines;
1956/7	1	6
1957/8	2	5
1958/9	3	4
1959/60	4	6
1960/1	5	6
1961/2	6	3
1962/3	7	12
1963/4	8	7
1964/5	9	4
1965/6	10	2
1966/7	11	6
1967/8	12	7
1968/9	13	4
;
run;

Figure 1.1 on page 14.

%let start = 3;
%let end = 8;
proc iml;
  use table1_2;
  read all var {number} into y ;  
  n = nrow(y) ;  
  c = j(&end-&start+1, 2, 0);  
  do i = 1 to &end - &start + 1;
     c[i, 1]= sum(y)*log(i+&start - 1) - n*(i+&start - 1);
	 c[i, 2] = i+&start - 1;
  end;
  create fig1_1 from c[colname={'lstar' 'theta'}];
  append from c;
quit;
symbol i = spline v=none;
axis1 order = (35 to 55 by 5) minor = none label=none;
proc gplot data = fig1_1;
  plot lstar*theta /vaxis = axis1 hminor=0;
run;
quit;

Table 1.3 on page 15. We use proc iml to implement the bisection calculation.

%let criterion = .00001;
%let left = 5;
%let right = 6;
%let itnum = 20;
proc iml;
  use table1_2;
  read all var {number} into y ; /*read the number of count into a matrix*/
  n = nrow(y) ; /*number of rows*/
  c = j(&itnum, 3, 0); /*create a matrix containing informatin for Table 1.3*/
  left = &left;
  right = &right;

  c[1,1] = 1;
  c[1, 2] = left;
  c[1, 3]= sum(y)*log(left) - n*left;

  c[2,1] = 2;
  c[2, 2] = right;
  c[2, 3]= sum(y)*log(right) - n*right;
   lb = c[1,3];
   ub = c[2,3];
  diff = ub - lb;
   	do i = 1 to &itnum while(abs(diff)>&criterion);
    	    mid = (left + right)/2;
			lstar = sum(y)*log(mid) - n*mid;
		c[i+2, 1] = i+2;
		c[i+2, 2] = mid;
		c[i+2, 3] = lstar;
    		if ub > mid & mid > lb then do;
        		lb = lstar; 
				left = mid;
			end;
	    	else if lb > mid & mid > ub then do ;
				right = mid;
				ub = lstar;
			end;
			else if lb < ub then do;
				left = mid;
				lb = lstar;
			end;
			else do;
				right = mid;
				ub = lstar;
			end;
			diff = ub - lb;
  end;
	 create c  var {k theta_k lstar};
	 append from c;
quit;
proc print data = c noobs;
  where k ~=0;
  format lstar 8.5;
run;

 K    THETA_K       LSTAR
 1    5.00000    50.87953
 2    6.00000    51.00668
 3    5.50000    51.24186
 4    5.75000    51.19239
 5    5.62500    51.23491
 6    5.56250    51.24293
 7    5.53125    51.24355
 8    5.54688    51.24352
 9    5.53906    51.24361
10    5.53516    51.24359
11    5.53711    51.24360

UCLA Researchers are invited to our Statistical Consulting Services
We recommend others to our list of Other Resources for Statistical Computing Help
These pages are Copyrighted (c) by UCLA Academic Technology Services