/*********  NEGATIVE MULTINOMIAL REGRESSION  *************
Author:
    Daniel A. Powers  mailto:dpowers@mail.la.utexas.edu
    Associate Professor, Sociology
    University of Texas at Austin
Date: 16 Nov 1999

This code is provided gratis without performance guarantees.

This program does negative multinomial regression
for clustered counts as described by
Guang Guo (1996) _Sociological Methodology_ 113-132.
It assumes a gamma distributed random effect w with
variance 1/phi.

The general model is
E[y(ij)|w(i)] = T(ij) exp{x(ij)B} w(i)
where
y(ij) is the count for the jth
individual in the ith cluster,
T(ij)is the length of the exposure interval
(an offset term).
Conditional on the cluster-level random effect,
individual counts within clusters are
independent.

Individuals with the same "id" are assumed to belong to the same cluster.

To apply this to panel data, each indivdual's repreated measures would
share the same id number.


*/



new;
/* define data */

load dta[2370,16]=medic.dat;
output file=negm.out reset;
m=rows(dta);
id=dta[.,1];
 y=dta[.,7];
 t=dta[.,3];

c=unique(id,1);
n=rows(c);

/* new data matrix (ndta) for use in model building */
const=ones(m,1);
i89=dta[.,9].*dta[.,13];
i90=dta[.,9].*dta[.,14];
i91=dta[.,9].*dta[.,15];

ndta=const~dta[.,9]~dta[.,13:15]~dta[.,6]~dta[.,2]~dta[.,4:5]~
i89~i90~i91;


/* variable list */


let dlist = const ovall yr89 yr90 yr91 pci hmoden bedden rmdden
  i89 i90 i91;

/* model variables (from data list) */

let model = const ovall yr89 yr90 yr91 pci hmoden bedden rmdden
i89 i90 i91;

names=model|"LOGPHI";

ni=indcv(model,dlist);
 if ismiss(ni);
 "  MODEL VARIABLE NAME DOES NOT MATCH NAME IN DATA LIST";?;
stop;
 endif;
    k=1;
      do until k > rows(ni);
        if k==1;
          x=ndta[.,ni[k]];
           elseif k >1;
            x=x~ndta[.,ni[k]];
           endif;
         k=k+1;
      endo;

if rows(model) /= cols(x);

/* don't know how this could happen but.... */

"NUMBER OF NAMES DOES NOT CORRESPOND TO NUMBER OF VARIABLES";?;

stop;
endif;


x[.,2]=x[.,2]./1000;
x[.,7]=x[.,7]./1000;
x[.,10]=x[.,10]./1000;
x[.,11]=x[.,11]./1000;
x[.,12]=x[.,12]./100;
k=1;
"             MEANS AND STANDARD DEVIATIONS             ";?;
 do until k > cols(x);
   $model[k];;meanc(x[.,k]);;stdc(x[.,k]);
  k=k+1;
endo;

k=cols(x);
w=ones(n,1);
w0=ones(m,1);
yx=y.*x;

{b,vb,l}=loglin(y,x,dta[.,3]);

b=b|0;


sumcount=zeros(n,1);
   y_vec=zeros(n,20);
   sumyx=zeros(n,cols(x));
      jj=zeros(n,1);

/* obtain event counts by cluster
     and make count matrix */

{sumcount,y_vec,jj}=clustcomp(y,id,c,n);
{sumyx}=clcomp(yx,id,c,n);

/* maximum cluster index (number of ind. contributions per cluster) */

mc=maxc(jj);



{b,vc1,vc2}=DFP(&fct,b,1.e-2,100,&dfct,w,0);


k=1;

se=sqrt(diag(vc2));
tv=b./se;

"PARAMETER      ESTIMATE      SE(b)         T-VAL
       ( dfp )              ";?;

  do until k > rows(b);
       $names[k];;b[k];;se[k];;tv[k];
     k=k+1;
endo;

k=1;

se=sqrt(diag(vc1));
tv=b./se;

"PARAMETER      ESTIMATE      SE(b)         T-VAL
       ( bhhh )              ";?;

  do until k > rows(b);
       $names[k];;b[k];;se[k];;tv[k];
     k=k+1;
endo;


proc(1)=grad(&fct,b,ll,dh,wt);endp;


/*CLUSTCOMP:   Sum cluster-specific information over all individuals */

proc(3) = clustcomp(y,id,c,n);
 local i, ii, j, jj, clust_sum, sum_cl, y_vec;
  m=rows(y);
 jj=zeros(n,1);
 sum_cl=zeros(n,1);
  y_vec=-ones(n,20);
 clust_sum=0;
 i=1;
 ii=0;
 j=0;
 jj=zeros(n,1);
do until i > m;
      if i eq 1; clust_sum=y[i];
        ii = 1;
         j = 1;
   else;
        if id[i] ne id[i-1];
         clust_sum=y[i];
        ii = ii+1;
         j = 1;
      elseif id[i] eq id[i-1];
        clust_sum=clust_sum+y[i];
        j=j+1;
     ii=ii;
   endif;
 endif;
jj[ii]=j;
   sum_cl[ii]   = clust_sum;
    y_vec[ii,j] = y[i];
i=i+1;
endo;
retp(sum_cl, y_vec ,jj);
endp;

proc(1) = clcomp(x,id,c,n);
 local i, ii, j, jj, clust_sum, sum_cl;
  m=rows(x);
  k=cols(x);
 jj=zeros(n,1);
 sum_cl=zeros(n,k);
 clust_sum=0;
 i=1;
 ii=0;
 j=0;
 jj=zeros(n,1);
do until i > m;
      if i eq 1; clust_sum=x[i,.];
        ii = 1;
         j = 1;
   else;
        if id[i] ne id[i-1];
         clust_sum=x[i,.];
        ii = ii+1;
         j = 1;
      elseif id[i] eq id[i-1];
        clust_sum=clust_sum + x[i,.];
        j=j+1;
     ii=ii;
   endif;
 endif;
jj[ii]=j;
   sum_cl[ii,.] = clust_sum;
i=i+1;
endo;
retp(sum_cl);
endp;




/* proc to do ML */


proc(1)=fct(b);
local a, u, u_vec, sumu, dlogu, dlogsum, dloggam1, dloggam2,
cc, d, logl, sum, i, j, np;
external proc clustcomp;

np=rows(b);
sumu=zeros(n,1);
u_vec=zeros(n,mc);
sum=zeros(n,1);
dlogsum=zeros(n,1);
dlogu=zeros(n,1);
cc = zeros(n,1);
logl=zeros(n,1);

/* reparam */

a=exp(b[np]);

/* reparam */

u=t.*exp(x*b[1:np-1]);

{sumu,u_vec,jj}=clustcomp(u,id,c,n);

/* i loop */
i=1;
 do until i > n;

  j=1;
    do until j > jj[i];
       dlogu[j] = ln(u_vec[i,j]);
  j=j+1;
 endo;

  dlogsum[i]= ln(sumu[i] + a);

  j=1;
    do until j > jj[i];
          cc[j] = y_vec[i,j] * (dlogu[j]-dlogsum[i]);
    sum[i]=sum[i]+cc[j];
  j=j+1;
  endo;

dloggam1 = lnfact((sumcount[i] + a) - 1 );
dloggam2 = lnfact(a-1);

       d = a * (ln(a) - ln(sumu[i]+a));

logl[i]  =  dloggam1 - dloggam2 + d + sum[i];
i=i+1;
endo;
retp(logl);
/* retp(-sumc(logl)); */
endp;


proc(1) = dfct(&fct,b,ll,dh,w);
/* proc(1) = dfct(b); */
local a, u, u_vec, sumu, dlogu, dlogsum, dloggam1, dloggam2, ux,
cc, d, logl, sum, i, j, np, fpd, sumux, digama1, digama2, nx,
ux_vec;
external proc clustcomp;
external proc clcomp;
external proc digamma;

np=rows(b);
nx=cols(x);
sumu=zeros(n,1);
u_vec=zeros(n,mc);
sumux=zeros(n,nx);
fpd=zeros(n,np);
a=exp(b[np]);
u=t.*exp(x*b[1:np-1]);
yx=y.*x;
ux=u.*x;


{sumu,u_vec,jj}=clustcomp(u,id,c,n);
{sumux}=clcomp(ux,id,c,n);


/* i loop */


i=1;
 do until i > n;
  digama1 = digamma(sumcount[i]+a);
  digama2 = digamma(a);
  fpd[i,np] =    digama1 - digama2 + 1 + ln(a)
          - ln(a + sumu[i] )
          - (a + sumcount[i]) / (a + sumu[i]);
   j=1;
     do until j > np-1;
            fpd[i,j] =   -(sumux[i,j]/(sumu[i]+a)) * (a+sumcount[i])
      + sumyx[i,j];
      j=j+1;
endo;
    i=i+1;
endo;
 retp(fpd);
/* retp(sumc(fpd)); */
endp;


proc digamma(x);
  local p;
  x=x+6;
  p=1/(x.*x);
  p=(((0.004166666666667*p-0.003968253986254).*p+
      0.008333333333333).*p-0.083333333333333).*p;
  p=p+ln(x)-0.5/x-1/(x-1)-1/(x-2)-1/(x-3)-1/(x-4)-1/(x-5)-1/(x-6);
retp(p);
endp;

/*   DFP */


Proc(3) = DFP(&f,b,tol,maxit,&grad,w,m);"DFP METHOD";?;
    local f:proc,grad:proc;
    external proc step;
    local k,s0,dg,dh,db,iter,ll,n,sll,g,s,shy,shys,yshy,ys,bn,i,vc,err,
          vc1, vc2, m1, m2;
         k = rows(b);
      s0   = 0;
      dg   = ones(k,1);
      dh   = 0;
      db   = ones(k,1);
      iter = 1;
     do until abs(dg) < tol or iter > maxit;
     ll = f(b);
      n = rows(ll);
    sll = sumc(ll);
      g = grad(&f,b,ll,dh,w);
      s = g'sqrt(w);
  /* if iter == 1; */
      if m == 0;
      m = eye(k)/abs(sll);
   endif;
     dg = s-s0;shy = db-m*dg;
   shys = shy*db'; yshy = dg'shy;
     ys = dg'db;
      m = m + ((shys + shys') - (yshy*(db*db'))/ys)/ys;
     db = m*s;

    bn = step(&f,b,db,s,iter);

?;
format /ld 8,0;  "Iteration " ;; iter;;
format /ld 12,6; "Function  " ;; sll;?;
"Estimate   Score Vector  Change Vector
_______________________________________";
  i = 1;
   do until i > k;
         bn[i];;s[i];;db[i];
     i = i+1;
   endo;
  ?;
  iter = iter + 1;
  db   = b - bn;
  b    = bn;
  s0   = s;
 vc1   = m;
endo;
   m2=moment(g,0);
   vc2=invpd(m2);

         err = diag(vc1);
      i=1;
       do until i > k;
         if err[i] <= 0;
" the variance of parameter ";; format /ld 3,0;i;;" is negative - set to 1.0 ";
         vc1[i,i] = 1;
         endif;
       i=i+1;
    endo;

         err = diag(vc2);
      i=1;
       do until i > k;
         if err[i] <= 0;
" the variance of parameter ";; format /ld 3,0;i;;" is negative - set to 1.0 ";
         vc2[i,i] = 1;
         endif;
       i=i+1;
    endo;

  retp(b,vc1,vc2);

endp;


proc(1) = grad(&f,x,f0,dh,w);
   local f:proc;
   local n,k,grdd,ax,dax,xdh,ee,i,eei;
      n = rows(f0); k = rows(x);  grdd = zeros(n,k);
        if dh == 0; ax = abs(x);
          if x /= 0; dax = x./ax; else; dax=1; endif;
     dh = (1e-5)*maxc(ax|(1e-2)*ones(k,1)).*dax;
     endif;
     xdh = x + dh ; dh = xdh - x; ee = eye(k).*dh;

i = 1;
   do until i>k;
     eei = ee[.,i];
     grdd[.,i] = f(x+eei)-f(x-eei);
    i = i+1;
  endo;
   grdd = grdd./(2*dh')./sqrt(w);
 retp(grdd);
endp;


    proc(1) = step(&f,b,db,g,iter);
    local f:proc;
    local delta,dpg,vr0,ub,lb,stepone,cdelta,s,vr1,gg1,siter,vr2,vrs,
          sprev,s2prev,rprev,r2prev,sprev2,s2prev2,sp2,dsprev,vv,zz,ab,
          aa,bb,err;
   delta  = 1e-4;
   dpg    = db'(-g);
   vr0    = sumc(f(b));   format /ld 12,6; "Function Evaluations";?;-vr0;
   ub     = 0.5;
   lb     = 0.1;
  stepone = 0.01; /* may need adjusting in some cases */
   cdelta = 1 - delta;
     if iter == 1; s = stepone; else; s = 1.0; endif;
    vr1 =  -sumc(f(b+s*db));  "Step-1"; -vr1;
       gg1 = (vr1 - vr0)/(s*dpg);
          if gg1 ge delta;
            goto slast;
          endif;
     siter = 1;
         s = -dpg/(2*(vr1-vr0-dpg));
          if s < lb; s = lb; endif;
        vr2 = -sumc(f(b+s*db));"Step-2";-vr2;
          gg1 = (vr2 - vr0)/(s*dpg);
           if gg1 ge delta and gg1 le cdelta;
            goto slast;
          endif;
          sprev = s; s2prev = 1; rprev = vr2; r2prev = vr1;
     siter = 2;
       do until siter > 10;
        sprev2 = sprev*sprev; s2prev2 = s2prev*s2prev; sp2 = sprev2~s2prev2;
        dsprev = sprev - s2prev;
        vv = (1~-1|-s2prev~sprev);vv = vv./sp2;
        zz = (rprev-vr0-dpg*sprev)|(r2prev-vr0-dpg*s2prev);
        ab = (1/dsprev)*vv*zz; aa=ab[1,1]; bb=ab[2,1];
          if aa == 0;
            s=-dpg/(2*bb);
           else;
                 err = (bb*bb - 3*aa*dpg);
                 if err[1,1] < 0;
" I'm Having some problems here. Please re-try with different starting values" ;
                pause(2);
                  goto slast; endif;
             s=(-bb + sqrt(bb*bb - 3*aa*dpg))/(3*aa);
          endif;
         if s > ub*sprev; s=ub*sprev;
          elseif s < lb*sprev; s=lb*sprev;
         endif;
        vrs = -sumc(f(b+s*db)); "Step-3"; vrs;
        gg1 = (vrs - vr0)/(s*dpg);
         if gg1 ge delta and gg1 le cdelta;
           goto slast;
         endif;
        s2prev = sprev; sprev = s;
        r2prev = rprev; rprev = vrs;
      siter = siter + 1;
     endo;
   SLAST:
   retp(b + s*db);
   endp;

 proc (1) = squeeze(y,&f,s_db,s_b,mode,w);
     local maxbkst,u,s_b1,s_lm,s_itr,lb1,s_b2,lb2,null,sqz;
     local f:proc;
  maxbkst = 20;
  s_b1=s_b + s_db;s_lm=1/2;s_itr=1;
        u =  (y-f(s_b1));
    if mode == 0;
      lb1 = (u.*w)'u;
    elseif mode == 1;
      lb1 = -sumc(w.*f(s_b1));
    endif;
      format /ld 12,6; "Function Evaluations";?;lb1;

 do until s_itr > maxbkst;format 14,6;
  s_b2 = s_b + s_lm * s_db;
        u = (y-f(s_b2));
    if mode == 0;
      lb2 = (u.*w)'u;
    elseif mode == 1;
       lb2 = -sumc(w.*f(s_b2));
    endif;
    lb2;
     if lb1 >= lb2;
      if mode == 0;
        s_b1 = s_b2; s_lm = s_lm/2; s_b2 = s_b + s_lm * s_db;lb1=lb2;
      elseif mode == 1;
        s_b1 = s_b2; s_lm = s_lm/2; s_b2 = s_b + s_lm * s_db;lb2=lb1;
      endif;
       s_itr = s_itr + 1;
     else; sqz = s_itr - 1;
    retp(s_b1);
     endif;
    endo;
  sqz = 10;
 retp(s_b2);
endp;



/* proc loglin:  MLE for Log-Linear/Log-Rate/Poisson Regression Models:

call as:

   {b,h,l} = loglin(d,x,t);

    Inputs:  t time to failure (1 for loglinear models)
             d number of events
             x vector of covariates
    Output:  b coeffiecent vector
             h hessian
             l log-likelihood
*/

proc(3) = loglin(d,x,t);

      local iter,db,tol,r,w,z,b,new_b,vcb,mse,vc,seb,h,logl;

b = invpd(x'x)*(x'ln(d+.5));

/* starting conditions */

iter  = 0;
db    = 1;
tol  = .1E-3;


"IRWLS Iterations on Log Rate Model:"; /* Note: set t=ones(rows(dta),1)
                    for loglinear model */
do until abs(db) < tol;

  r   = t.*exp(x*b);
  w   = r;
  z   = x*b + (d-r)./r;

  new_b = invpd((x.*w)'x)*((x.*w)'z);
       db   = new_b - b;
        b   = new_b;
        h   = (x.*r)'x;
iter  = iter + 1;?;
b';
endo;

Logl = sumc(d.*x*b - r);
retp(b,h,logl);
endp;
