/*
**                         FILE BPM19.ARC
**
** This file contains the bivariate probit model procedures.
** Part 19.
*/


/*CASE INFLUENCE*/



/* Declaration of globals. */

declare _rpar ?=0;
declare _cases ?=0;
declare _file ?=0;


/* Selection of appropriate libraries. */

library bpm,pgraph;
graphset;
hardcopy("5");

/*
** INFL
**
** FORMAT
**                           infl(dataset,indvar,depvar);
**
** INPUT      dataset  : Either a string, the name of the dataset that
**                       contains the variables or the null string.
**            indvar   : Either an Lx1 or (L-1)x1 character vector containing
**                       the names of the independent variables or,
**                       in case dataset is the null string, an
**                       NxL or Nx(L-1) matrix of independent variables.
**                       (L-1 if the first column isn't a unit column; in that
**                       case a unit column is appended and L-1 is incremented).
**            depvar   : Either a 2x1 character vector containing the names
**                       of the dependent variables or, in case dataset is the
**                       null string, an Nx2 matrix of dependent 0-1 variables.
**
**
** GLOBALS    _rpar:  0: The correlation coefficient is a constant to be
**                       estimated (default).
**                    1: The correlation depends on the independent variables
**                       in the same way as the variables do.
*/



proc (0)= infl(dataset,indvar,depvar);

local x,y,teta,ll,teta1,ll1,teta2,ll2,tetra,llt,se,i,j,k,g1,g2,cnt,s,yhulp,
      indexi,indexd,nvec,handle,hulp1,hulp2,cnt,test,t,hulp,afg1,teta1,var,
      covar,sortv,afg2,mult,dist,teller,teta11,xas,yas,string,l,yas1,cooka,
      yas2,yas3,yas4,cookb,cookr,hata,hatb,hatr,xa,xb,xr,jj,v,boven,nummer,
      p,een,twee,drie,vier,obscol,q,maximum,grond,plafond;


disable;

/*
** If the input is on a dataset, it is processed in the subsequent, if not
** it is read from the matrices passed in indvar and depvar.
*/

if dataset$=="";

   x=indvar;
   y=depvar;
   let depvar="y1" "y2";
   let indvar="constant" "x1" "x2" "x3" "x4" "x5" "x6" "x7" "x8" "x9";
   if x[.,1]==ones(rows(x),1);
      indvar=indvar[1:cols(x)];
   else;
      indvar=indvar[1:cols(x)+1];
   endif;

else;

   {indvar,indexi}=indices(dataset,indvar);
   {depvar,indexd}=indices(dataset,depvar);
   nvec=setvars(dataset);
   open handle=^dataset for read;
   do until eof(handle);
      hulp1=readr(handle,100);
      i=1;
      do while i<=rows(indvar);
         hulp2=varget(indvar[i]);
         if rows(hulp2)>1;
            hulp2=hulp2|hulp1[.,indexi[i]];
         else;
            hulp2=hulp1[.,indexi[i]];
         endif;
         cnt=varput(hulp2,indvar[i]);
         i=i+1;
      endo;

      i=1;
      do while i<=rows(depvar);
         hulp2=varget(depvar[i]);
         if rows(hulp2)>1;
            hulp2=hulp2|hulp1[.,indexd[i]];
         else;
            hulp2=hulp1[.,indexd[i]];
         endif;
         cnt=varput(hulp2,depvar[i]);
         i=i+1;
      endo;

   endo;
   handle=close(handle);
   y=varget(depvar[1])~varget(depvar[2]);
   x=varget(indvar[1]);
   i=2;
   do while i<=rows(indvar);
      x=x~varget(indvar[i]);
      i=i+1;
   endo;

   if not(x[.,1]==ones(rows(x),1));
      indvar="constant"|indvar;
   endif;

   nvec=setvars(dataset);

endif;



/* Checking for faulty data and deletion of rows with missing values. */

if (rows(x)/=rows(y));
   errorlog "x en y must have the same number of rows.";
   end;
elseif (cols(y)/=2);
   errorlog "y does not stand for a table of 0-1 variables.";
   end;
endif;

t=x~y;
_cases=rows(t);
g1=cols(x);
g2=cols(t);
nummer=seqa(1,1,rows(x));
t=packr(t~nummer);
nummer=t[.,cols(t)];
_cases=_cases|rows(t);
y=t[.,g1+1:g2];
x=t[.,1:g1];

if sumc(sumc((y[.,1 2] .ne 1) .and (y[.,1 2] .ne 0)))/=0;
   errorlog "The dependent variables are not of the desired form.";
   end;
endif;


/* Transforming the dependent variables into contingency tables. */

hulp=zeros(rows(y),4);
hulp[.,1]=((y[.,1] .== 0) .and (y[.,2] .== 0));
hulp[.,2]=((y[.,1] .== 1) .and (y[.,2] .== 0));
hulp[.,3]=((y[.,1] .== 0) .and (y[.,2] .== 1));
hulp[.,4]=((y[.,1] .== 1) .and (y[.,2] .== 1));

y=hulp;
clear hulp;
g2=g2+2;


/* Multiple columns sort on the independent variable(s). */

t=x~y~nummer;
sortv=seqa(1,1,g1);
t=sortmc(t,sortv);
nummer=t[.,cols(t)];
t=t[.,1:cols(t)-1];

obscol=0~0~0;
teller=0;


/* Summing over equal independent variable vectors. */

i=1;
k=1;
s=zeros(rows(t),g2);
do while i<=rows(t);
   cnt=t[i,1:g1];
   j=i;
   test=1;
   do while (j<=rows(t)) and test;
      j=j+1;
      if j<=rows(t);
         test=t[i,1:g1]==t[j,1:g1];
      else;
         test=0;
      endif;
   endo;
   j=j-1;
   if j<=rows(t);
      /* Observatienummer bijhouden */
      een=0;
      twee=0;
      drie=0;
      vier=0;
      p=i;
      do while p<=j;
         if t[p,g1+1]==1;
            een=een|nummer[p];
         elseif t[p,g1+2]==1;
            twee=twee|nummer[p];
         elseif t[p,g1+3]==1;
            drie=drie|nummer[p];
         elseif t[p,g1+4]==1;
            vier=vier|nummer[p];
         endif;
         p=p+1;
      endo;
      if rows(een)>1;
         teller=teller+1;
         obscol=obscol|(een[2:rows(een)]~((teller~1).*ones(rows(een)-1,2)));
      endif;
      if rows(twee)>1;
         teller=teller+1;
         obscol=obscol|(twee[2:rows(twee)]~((teller~2).*ones(rows(twee)-1,2)));
      endif;
      if rows(drie)>1;
         teller=teller+1;
         obscol=obscol|(drie[2:rows(drie)]~((teller~3).*ones(rows(drie)-1,2)));
      endif;
      if rows(vier)>1;
         teller=teller+1;
         obscol=obscol|(vier[2:rows(vier)]~((teller~4).*ones(rows(vier)-1,2)));
      endif;
      /* Einde observatienummer bijhouden */
      s[k,g1+1:g2]=sumc(t[i:j,g1+1:g2])';
      s[k,1:g1]=cnt;
      k=k+1;
   endif;
   i=j+1;
endo;

s=s[1:k-1,.];
obscol=obscol[2:rows(obscol),.];


/*
** Separating dependent and independent variables and adding a column
** of ones to the matrix of independent variables if necessary.
*/

y=s[.,cols(s)-3:cols(s)];
x=s[.,1:cols(s)-4];
if not(x[.,1]==1);
   x=ones(rows(x),1)~x;
endif;
l=cols(x);


/* Calling BIVPROB. */

{teta,ll,tetra,llt,teta1,ll1,teta2,ll2}=bivprob(x,y,indvar,depvar);


/* Case deletion. */

se=teta[.,2];
teta=teta[.,1];
teta1=zeros(rows(teta),1);
covar=zeros(cols(x),1);
var=zeros(2,1);
mult=0;
teller=0;
dist=zeros(rows(teta)+1,1);
cooka=zeros(2,1);
cookb=zeros(2,1);
cookr=zeros(2,1);

cls;

if _rpar==0;
   xa=eye(l)|zeros(l+1,l);
   xb=zeros(l,l)|eye(l)|zeros(1,l);
   xr=zeros(2*l,1)|1;
   v=(-hessbp(teta,x,y));
   jj=v*xa;
   hata=jj'*inv(v)*jj;
   jj=v*xb;
   hatb=jj'*inv(v)*jj;
   jj=v*xr;
   hatr=jj'*inv(v)*jj;
else;
   xa=eye(l)|zeros(2*l,l);
   xb=zeros(l,l)|eye(l)|zeros(l,l);
   xr=zeros(2*l,l)|eye(l);
   v=(-hessbpr(teta,x,y));
   jj=v*xa;
   hata=jj'*inv(v)*jj;
   jj=v*xb;
   hatb=jj'*inv(v)*jj;
   jj=v*xr;
   hatr=jj'*inv(v)*jj;
endif;

p=1;
i=1;
do while i<=rows(y);
   j=1;
   do while j<=4;
      if y[i,j]>0;
         teller=teller+1;
         yhulp=y;
         yhulp[i,j]=yhulp[i,j]-1;
         if _rpar==0;
            afg1=gradbp(teta,x,yhulp);
            afg2=hessbp(teta,x,yhulp);
         elseif _rpar==1;
            afg1=gradbpr(teta,x,yhulp);
            afg2=hessbpr(teta,x,yhulp);
         endif;
         teta11=teta-inv(afg2)*afg1;
         teta1=teta1~teta11;
         dist=dist~(j|((teta-teta11)./(se)));
         cooka=cooka~(j|sqrt((teta-teta11)'*xa*hata*xa'*(teta-teta11)));
         cookb=cookb~(j|sqrt((teta-teta11)'*xb*hatb*xb'*(teta-teta11)));
         cookr=cookr~(j|sqrt((teta-teta11)'*xr*hatr*xr'*(teta-teta11)));
         covar=covar~x[i,.]';
         if j==1;
            var=var~(0|0);
         elseif j==2;
            var=var~(1|0);
         elseif j==3;
            var=var~(0|1);
         elseif j==4;
            var=var~(1|1);
         endif;
         mult=mult~y[i,j];
         if _file/=2;
            output on;
         endif;
         cls;
         format 4,0;
         print "Combination number:" teller;
         print "       Case number:";;
         q=sumc(obscol[.,2].==teller);
         print obscol[p:p+q-1,1]';
         p=p+q;
         print "      Multiplicity:" mult[cols(mult)];
         print "         Variables:";;
         print var[.,cols(var)]';
         format 10,4;
         print "        Covariates:";;
         print x[i,.];
         print "    Cook distances:";;
         print cooka[2,rows(cooka)] cookb[2,rows(cookb)] cookr[2,rows(cookr)];
         format 10,4;
         print "        Parameters:";;
         print  teta1[.,cols(teta1)]';
         print "       Differences:";;
         print dist[2:rows(dist),cols(dist)]';
         print "************************************************************"\
         "*****************";
         print " ";
         if _file/=2;
            output off;
         endif;
      endif;
      j=j+1;
   endo;
   i=i+1;
endo;

var=var[.,2:cols(var)];
covar=covar[.,2:cols(covar)];
teta1=teta1[.,2:cols(teta1)];
mult=mult[2:cols(mult)];
dist=dist[.,2:cols(dist)];
cooka=cooka[.,2:cols(cooka)];
cookb=cookb[.,2:cols(cookb)];
cookr=cookr[.,2:cols(cookr)];


/* Graphics. */

obscol=sortc(obscol,1);
graphset;
_plegctl=2~2~7.85~5.91;
_plegstr="(1,1)\000(2,1)\000(1,2)\000(2,2)";
_pframe=0|0;
/*_pgrid=2|2;*/
_pgrid=0|0;
_pbox=1;
_plctrl=-1;
_pnumht=.125;
_psilent=1;
_pstype=8~9~10~12;
let _pcolor=10 11 12 14;
_pnum=2;
_psymsiz=6;
_pcross=1;
xas=obscol[.,1];
if rows(xas)<=50;
   xtics(0,maxc(xas)+1,4,4);
elseif rows(xas)>=250;
   xtics(0,maxc(xas)+1,20,0);
else;
   xtics(0,maxc(xas)+1,4,0);
endif;
boven=maxc(maxc(ceil(10*cooka[2:rows(cooka),.])))|
      maxc(maxc(ceil(10*cookb[2:rows(cookb),.])))|
      maxc(maxc(ceil(10*cookr[2:rows(cookr),.])));
boven=maxc(boven)+0.1;
dist=dist';


/* Cook distance for first dependent variable. */

cooka=cooka';
ytics(0,boven/10,0.25,5);
string="CASE INFLUENCE. COOK DISTANCE.\l";
string=string$+"Variable: "$+depvar[1];
title(string);
_pmsgstr="Distance\000Case Number";
_pmsgctl={0.5 3 0.2 90 1 15,3.5 0.2 0.2 0 1 15};
yas=cooka[obscol[.,2],2];
yas1=yas.*(obscol[.,3].==1)+100.*(obscol[.,3]./=1);
yas2=yas.*(obscol[.,3].==2)+100.*(obscol[.,3]./=2);
yas3=yas.*(obscol[.,3].==3)+100.*(obscol[.,3]./=3);
yas4=yas.*(obscol[.,3].==4)+100.*(obscol[.,3]./=4);
maximum=maxc(yas);
hulp=maximum.==yas;
if rows(yas)-sumc(hulp.==0)==1;
   hulp=xas.*hulp;
   hulp=sumc(hulp);
   _pmsgstr=_pmsgstr$+"\000"$+ftos(hulp,"lf",4,0);
   _pmsgctl=_pmsgctl|zeros(1,6);
   _pmsgctl[3,1]=hulp;
   _pmsgctl[3,2]=maximum;
   _pmsgctl[3,3]=0.1;
   _pmsgctl[3,4]=0;
   _pmsgctl[3,5]=0;
   _pmsgctl[3,6]=15;
endif;
xy(xas,yas1~yas2~yas3~yas4);


/* Individual distances for first dependent variable. */

grond=floor((minc(minc(floor(dist[.,2:cols(dist)]*10)))-1)/5)/2;
plafond=ceil((maxc(maxc(ceil(dist[.,2:cols(dist)]*10)))+1)/5)/2;
ytics(grond,plafond,0.5,5);
j=1;
do while j<=cols(x);
   string="CASE INFLUENCE. Signed distance per covariate\l";
   string=string$+"Variable: "$+depvar[1]$+".  Covariate :"$+indvar[j];
   title(string);
   yas=dist[obscol[.,2],1+j];
   yas1=yas.*(obscol[.,3].==1)+100.*(obscol[.,3]./=1);
   yas2=yas.*(obscol[.,3].==2)+100.*(obscol[.,3]./=2);
   yas3=yas.*(obscol[.,3].==3)+100.*(obscol[.,3]./=3);
   yas4=yas.*(obscol[.,3].==4)+100.*(obscol[.,3]./=4);
   maximum=minc(yas)|maxc(yas);
   if abs(maximum[1])>abs(maximum[2]);
      maximum=maximum[1];
   else;
      maximum=maximum[2];
   endif;
   hulp=maximum.==yas;
   _pmsgstr="Distance\000Case Number";
   _pmsgctl={.5 3 0.2 90 1 15,3.5 0.2 0.2 0 1 15};
   if rows(yas)-sumc(hulp.==0)==1;
      hulp=xas.*hulp;
      hulp=sumc(hulp);
      _pmsgstr=_pmsgstr$+"\000"$+ftos(hulp,"lf",4,0);
      _pmsgctl=_pmsgctl|zeros(1,6);
      _pmsgctl[3,1]=hulp;
      _pmsgctl[3,2]=maximum;
      _pmsgctl[3,3]=0.1;
      _pmsgctl[3,4]=0;
      _pmsgctl[3,5]=0;
      _pmsgctl[3,6]=15;
   endif;
   xy(xas,yas1~yas2~yas3~yas4);
   j=j+1;
endo;


/* Cook distance for second dependent variable. */

cookb=cookb';
ytics(0,boven/10,0.25,5);
string="CASE INFLUENCE. COOK DISTANCE.\l";
string=string$+"Variable: "$+depvar[2];
title(string);
_pmsgstr="Distance\000Case Number";
_pmsgctl={0.5 3 0.2 90 1 15,3.5 0.2 0.2 0 1 15};
yas=cookb[obscol[.,2],2];
yas1=yas.*(obscol[.,3].==1)+100.*(obscol[.,3]./=1);
yas2=yas.*(obscol[.,3].==2)+100.*(obscol[.,3]./=2);
yas3=yas.*(obscol[.,3].==3)+100.*(obscol[.,3]./=3);
yas4=yas.*(obscol[.,3].==4)+100.*(obscol[.,3]./=4);
maximum=maxc(yas);
hulp=maximum.==yas;
if rows(yas)-sumc(hulp.==0)==1;
   hulp=xas.*hulp;
   hulp=sumc(hulp);
   _pmsgstr=_pmsgstr$+"\000"$+ftos(hulp,"lf",4,0);
   _pmsgctl=_pmsgctl|zeros(1,6);
   _pmsgctl[3,1]=hulp;
   _pmsgctl[3,2]=maximum;
   _pmsgctl[3,3]=0.1;
   _pmsgctl[3,4]=0;
   _pmsgctl[3,5]=0;
   _pmsgctl[3,6]=15;
endif;
xy(xas,yas1~yas2~yas3~yas4);



/* Individual distances for second dependent variable. */

ytics(grond,plafond,0.5,5);
j=1;
do while j<=cols(x);
   string="CASE INFLUENCE. Signed distance per covariate \l";
   string=string$+"Variable: "$+depvar[2]$+".  Covariate :"$+indvar[j];
   title(string);
   yas=dist[obscol[.,2],1+l+j];
   yas1=yas.*(obscol[.,3].==1)+100.*(obscol[.,3]./=1);
   yas2=yas.*(obscol[.,3].==2)+100.*(obscol[.,3]./=2);
   yas3=yas.*(obscol[.,3].==3)+100.*(obscol[.,3]./=3);
   yas4=yas.*(obscol[.,3].==4)+100.*(obscol[.,3]./=4);
   maximum=minc(yas)|maxc(yas);
   if abs(maximum[1])>abs(maximum[2]);
      maximum=maximum[1];
   else;
      maximum=maximum[2];
   endif;
   hulp=maximum.==yas;
   _pmsgstr="Distance\000Case Number";
   _pmsgctl={.5 3 0.2 90 1 15,3.5 0.2 0.2 0 1 15};
   if rows(yas)-sumc(hulp.==0)==1;
      hulp=xas.*hulp;
      hulp=sumc(hulp);
      _pmsgstr=_pmsgstr$+"\000"$+ftos(hulp,"lf",4,0);
      _pmsgctl=_pmsgctl|zeros(1,6);
      _pmsgctl[3,1]=hulp;
      _pmsgctl[3,2]=maximum;
      _pmsgctl[3,3]=0.1;
      _pmsgctl[3,4]=0;
      _pmsgctl[3,5]=0;
      _pmsgctl[3,6]=15;
   endif;
   xy(xas,yas1~yas2~yas3~yas4);
   j=j+1;
endo;


/* Cook distance for correlation. */

cookr=cookr';
ytics(0,boven/10,0.25,5);
string="CASE INFLUENCE. COOK DISTANCE.\l";
string=string$+"Correlation.";
title(string);
_pmsgstr="Distance\000Case Number";
_pmsgctl={0.5 3 0.2 90 1 15,3.5 0.2 0.2 0 1 15};
yas=cookr[obscol[.,2],2];
yas1=yas.*(obscol[.,3].==1)+100.*(obscol[.,3]./=1);
yas2=yas.*(obscol[.,3].==2)+100.*(obscol[.,3]./=2);
yas3=yas.*(obscol[.,3].==3)+100.*(obscol[.,3]./=3);
yas4=yas.*(obscol[.,3].==4)+100.*(obscol[.,3]./=4);
maximum=maxc(yas);
hulp=maximum.==yas;
if rows(yas)-sumc(hulp.==0)==1;
   hulp=xas.*hulp;
   hulp=sumc(hulp);
   _pmsgstr=_pmsgstr$+"\000"$+ftos(hulp,"lf",4,0);
   _pmsgctl=_pmsgctl|zeros(1,6);
   _pmsgctl[3,1]=hulp;
   _pmsgctl[3,2]=maximum;
   _pmsgctl[3,3]=0.1;
   _pmsgctl[3,4]=0;
   _pmsgctl[3,5]=0;
   _pmsgctl[3,6]=15;
endif;
xy(xas,yas1~yas2~yas3~yas4);


/* Individual distances for correlation. */

_pmsgstr="Distance\000Case Number";
_pmsgctl={.5 3 0.2 90 1 15,3.5 0.2 0.2 0 1 15};
ytics(grond,plafond,0.5,5);
if _rpar==0;
   string="CASE INFLUENCE. Signed distance.\l Correlation";
   title(string);
   yas=dist[obscol[.,2],1+2*l+1];
   yas1=yas.*(obscol[.,3].==1)+100.*(obscol[.,3]./=1);
   yas2=yas.*(obscol[.,3].==2)+100.*(obscol[.,3]./=2);
   yas3=yas.*(obscol[.,3].==3)+100.*(obscol[.,3]./=3);
   yas4=yas.*(obscol[.,3].==4)+100.*(obscol[.,3]./=4);
   maximum=minc(yas)|maxc(yas);
   if abs(maximum[1])>abs(maximum[2]);
      maximum=maximum[1];
   else;
      maximum=maximum[2];
   endif;
   hulp=maximum.==yas;
   _pmsgstr="Distance\000Case Number";
   _pmsgctl={.5 3 0.2 90 1 15,3.5 0.2 0.2 0 1 15};
   if rows(yas)-sumc(hulp.==0)==1;
      hulp=xas.*hulp;
      hulp=sumc(hulp);
      _pmsgstr=_pmsgstr$+"\000"$+ftos(hulp,"lf",4,0);
      _pmsgctl=_pmsgctl|zeros(1,6);
      _pmsgctl[3,1]=hulp;
      _pmsgctl[3,2]=maximum;
      _pmsgctl[3,3]=0.1;
      _pmsgctl[3,4]=0;
      _pmsgctl[3,5]=0;
      _pmsgctl[3,6]=15;
   endif;
   xy(xas,yas1~yas2~yas3~yas4);
else;
   j=1;
   do while j<=cols(x);
      string="CASE INFLUENCE. Individual distance.\l";
      string=string$+"Correlation.   Covariate :"$+indvar[j];
      title(string);
      yas=dist[obscol[.,2],1+2*l+j];
      yas1=yas.*(obscol[.,3].==1)+100.*(obscol[.,3]./=1);
      yas2=yas.*(obscol[.,3].==2)+100.*(obscol[.,3]./=2);
      yas3=yas.*(obscol[.,3].==3)+100.*(obscol[.,3]./=3);
      yas4=yas.*(obscol[.,3].==4)+100.*(obscol[.,3]./=4);
      maximum=minc(yas)|maxc(yas);
      if abs(maximum[1])>abs(maximum[2]);
         maximum=maximum[1];
      else;
         maximum=maximum[2];
      endif;
      hulp=maximum.==yas;
      _pmsgstr="Distance\000Case Number";
      _pmsgctl={.5 3 0.2 90 1 15,3.5 0.2 0.2 0 1 15};
      if rows(yas)-sumc(hulp.==0)==1;
         hulp=xas.*hulp;
         hulp=sumc(hulp);
         _pmsgstr=_pmsgstr$+"\000"$+ftos(hulp,"lf",4,0);
         _pmsgctl=_pmsgctl|zeros(1,6);
         _pmsgctl[3,1]=hulp;
         _pmsgctl[3,2]=maximum;
         _pmsgctl[3,3]=0.1;
         _pmsgctl[3,4]=0;
         _pmsgctl[3,5]=0;
         _pmsgctl[3,6]=15;
      endif;
      xy(xas,yas1~yas2~yas3~yas4);
      j=j+1;
   endo;
endif;



endp;


