data abcc;
  input A S fd1-fd6;
  array f{6} fd1-fd6;
  array dose{6} (0 5 30 75 150 300);
  do d=0 to 5;
    freq=f{d+1};
    dlvl=dose{d+1};
    output; drop fd1-fd6;
  end;

cards;
0 0    0     7    3   1   4  11
0 1 5015 10752 2989 694 418 387
1 0    5     4    6   1   3   6
1 1 5973 11811 2620 711 792 820
2 0    2     8    3   1   3   7
2 1 5669 10828 2798 797 596 624
3 0    3    19    4   2   1  10
3 1 6158 12645 3566 972 694 608
4 0    3     7    3   2   2   6
4 1 3695  9053 2415 655 393 289
;;;;

proc print data=abcc;
run;

data catdat;
  set abcc;
  if freq=0 then freq=1.e-14;

proc catmod data=catdat;
  weight freq;
  model a*s*d=_response_ / pred=freq;
  loglin  a|d s|d;
run;

* collapsed table:

proc freq data=abcc noprint;
  weight freq;
  tables s*d*dlvl / out=sdtable(rename=(count=freq));


* Linear by linear association model;

Set up design variables;

data sdtable;
  set sdtable;
  array da{5} d1-d5;
  array sda{5} sd1-sd5;
  sd=s*dlvl;
  u=1;
  do j=1 to 5;
    da{j}=(d=j);
   sda{j}=(d=j)*(s=1);
  end;
run;

proc print data=sdtable;
run;

* Poisson Regression by IRLS;

options nocenter ls=76 ps=54;
proc iml;
  eps=.0000001;
  panic=25;
  use sdtable;
    read all var {freq} into n;
    read all
     var ("u"||"s"||("d1":"d5")||"sd")
     into x;
  p=nrow(n);
  r=ncol(x);

  beta=log(n[+]/p)//J(r-1,1,0);
  nupd=0;
  do until((abs(beta-oldb)[+] < eps) | (nupd > panic));
    nupd=nupd+1;
    mu=exp(x*beta);
    W=diag(mu);
    ystar=X*beta + (n-mu)/mu;
    oldb=beta;
    Info = X`*W*X;
    beta=solve(Info,X`*W*ystar);
  end;
  sebeta=sqrt(vecdiag(inv(info)));
  varmat= (   W  ||  W*X) //
          ( X`*W || Info);
  seres=sqrt(vecdiag(sweep(varmat,(p+1:p+r)))[1:p]);
  g2=2*(n#log((n||J(p,1,1))[,<>]/mu))[+]; df=p-r;
  resid=n-mu;
  zres=resid/seres;
  score=abs(X`*resid)[<>];
  print nupd score;
  print g2 df;
  print beta sebeta;
  print n mu resid seres zres;
quit;
run;