/*                                                                */
/* Read in data file. Assign the name "pop" to county population. */
/* Assign the name "d_emg" to distance to emergency care.         */
/* Assign the name "CODE" to the subject ID.         */
/*                                                                */
/* The values for guess can be changed from {0, 0}.               */
/*                                                                */
/* The values of An^2 for various {lambda_1, lambda_2} are output */
/* to a file named power.out.                                     */
/*                                                                */
/* Rurality index scores for all subjects and for the selected    */
/* value of {lambda_1, lambda_2} are output to a file named       */
/* ri.txt.                                                        */
/*                                                                */


DM 'LOG; CLEAR; OUT; CLEAR;';

data distance;
   infile 'diabetespop.dat';
   input case d_emg pop;
   CODE=case;
   use=1;
   if d_emg = 0 then d_emg=.5;
   if d_emg =. then use=0;
   if pop = . then use=0;
   keep d_emg pop CODE;
   if use = 1 then output;

 proc means;
 proc freq;
     tables d_emg pop;
 
 proc iml;
    use distance; reset noname; read all var {d_emg pop} into YY;
    read all var {CODE} into IDnum;
    step=.05; mino=0; maxo=4.;
    k=2; c={.5, -1}; n=nrow(YY);
    do i=1 to k; YY[,i]=YY[,i]/max(YY[,i]); end;
    j=(1:n)`; p=(j-.5)/n; pi=arcos(-1);
    X=J(n,k,0);  X1=J(n,k,0);  X2=J(n,k,0);
    eps=1/2**26;


    start deriv;
    tau=lam-1; omega=(ssq(tau)); romega=sqrt(omega); IMP=I(k);
    if romega >= eps then do;
    tau=tau/romega; IMP=IMP-tau*tau`; end;
    do j=1 to k;
      y=YY[,j];
      if lam[j] =0 then do;
        t=log(y);
        t1=t#t/2;
        t2=2*t1#t/3;  end;
      else do;
        t=sign(lam[j])*y##lam[j];
        t1=t#log(y);  t2=t1#log(y); end;
      t=t-J(n,1,sum(t)/n); v0=ssq(t);  xj=t/sqrt(v0); X[,j]=xj;
      t1=t1-J(n,1,sum(t1)/n);  t2=t2-J(n,1,sum(t2)/n);
      v01=xj`*t1;  v02=xj`*t2;
      x1j=(t1-xj*v01)/sqrt(v0);  X1[,j]=x1j; v11=x1j`*t1;
      x2j=(t2-xj*v02-2*x1j*v01-xj*v11)/sqrt(v0); X2[,j]=x2j;
      end;
      z=X*c;  s2=z`*z/(n-1); s=sqrt(s2);  z=z/s;
      ii=rank(z);  B=z; z[ii]=B;
      B=X;  X[ii,]=B;  B=X1; X1[ii,]=B;  B=X2; X2[ii,]=B;
      G=sqrt(omega/s2)*diag(c)*(X1`-(X1`*z)*z`/(n-1));
      Phi=probnorm(z);
      Phi = Phi#(eps <= Phi)#(Phi <= 1-eps) + J(n,1,eps)#(Phi < eps)
          + J(n,1,1-eps)#(1-eps < Phi);
      logphi=log(Phi); logphi1=log(1-Phi);
      ph=exp(-z#z/2)/sqrt(2*pi);  logph=-(z#z+J(n,1,log(2*pi)))/2;
      A2n=2*p`*logphi1-2*p`*logphi-2*sum(logphi1)-n;
      a1=-2*exp(logph-logphi-logphi1)#(p-Phi);
      a2=2*exp(2*logph-2*logphi-2*logphi1)#((p-Phi)#(p-Phi)+p#(1-p));
      a2=a2+2*exp(logph-logphi-logphi1)#z#(p-Phi);
      d1=G*a1; f=IMP*d1;
      d2=(omega/s)*diag(c)*(X2`*a1-(X2`*z)*(z`*a1)/(n-1));
      d3=sqrt(omega/s2)*IMP*diag(c)*X2`*z;
      H=-(f*tau`+tau*f`);  H1=diag(d2)-(tau`*d1)*I(k)-(d1*d3`+d3*d1`)/(n-1);
      H1=H1-G*G`*(z`*a1)/(n-1)+(G#(J(k,1,1)*a2`))*G`;
      H=H+IMP*H1*IMP;
      finish;

     start optimize;
     delta=1;  A2n0=n;
     do while (delta > 1.e-5);
     run deriv;
     call gsorth(T,u,lindep,I(2)-tau*tau`);  T=T[,1];
     b=-inv(T`*H*T)*(T`*f);
     tau=tau+T*b;
     tau=tau/sqrt(ssq(tau));
     lam=1+tau*romega;
     delta=n*abs(A2n0-A2n);
     A2n0=A2n;
     end;
     finish;


    start search;
    filename out1 'power.out';
    file out1;
    omega0=ssq(guess-1); romega0=sqrt(omega0); tau=(guess-1)/romega0;
    lam=guess; run deriv;
    put @2 romega0 7.4 +3 A2n 9.4+3 (lam[1]) 7.4+3 (lam[2]) 7.4+3;
    lam=guess;
    maxo=(maxo-mino)+romega0; mino=romega0;
    do romega=mino to maxo by step;
    omega=romega**2;
    tau=lam-1;  tau=tau/sqrt(ssq(tau));  lam=1+romega*tau;
    run optimize;
    put @2 romega 7.4 +3 A2n 9.4+3 (lam[1]) 7.4+3 (lam[2]) 7.4+3;
    end;
    closefile out1;
    finish;


    start compute_ri;
    filename out2 'ri.txt';
    file out2;
    lam1=lam[1]; lam2=lam[2];
    if lam1=0 then do;
        tmiles=log(YY[,1]);
    end;
    else do;
        tmiles=(YY[,1]##lam1-1)/lam1;
    end;
    if lam2=0 then do;
        tpop=log(YY[,2]);
    end;
    else do;
        tpop=(YY[,2]##lam2-1)/lam2;
    end;
    m_tmiles=sum(tmiles)/n;
    m_tpop=sum(tpop)/n;
    std_miles=sqrt(ssq(tmiles-m_tmiles)/(n-1));
    std_pop=sqrt(ssq(tpop-m_tpop)/(n-1));
    tmiles=(tmiles-m_tmiles)/std_miles;
    tpop=(tpop-m_tpop)/std_pop;
    RI=.5*tmiles-tpop;
    m_RI=sum(RI)/n;
    std_RI=sqrt(ssq(RI-m_RI)/(n-1));
    RI=(RI-m_RI)/std_RI;
    do II = 1 to n;
      T1=IDnum[II]; T2=RI[II];
      put @2 T1 6.0 +3 T2 9.5+3;
    end;
    print IDnum RI [colname={"RI" "Code"}];
    closefile out2;
    finish;

     

/*  Try various values of {lambda_1, lambda_2} to choose */
/*  an initial guess.                                    */
/*                                                       */
    lam={0.2,-.1}; run deriv; print lam A2n;
    lam={0,-.15}; run deriv; print lam A2n;
    lam={0,-.2}; run deriv; print lam A2n;
    lam={1.,1.}; run deriv; print lam A2n;
    lam={0.,0.}; run deriv; print lam A2n;
    lam={-1.0,-1.5}; run deriv; print lam A2n;
    guess={0, 0};
    run search;
/*                                                       */
/*  Compute the rurality index scores for a selected     */
/*  pair {lambda_1, lambda_2}.                           */
/*                                                       */
    lam={0.0,0.0}; run compute_ri;