*******************************************

Program: BDM
Purpose: Fit the Bivariate Dale Model in SAS.

REQUIRES SAS VERSION 8.2 (OR GREATER) AND SAS/STAT

Garnett P McMillan	gnet@unm.edu
and
Tim Hanson	hanson@math.unm.edu

Created:  10/24/02
Modified: 11/17/04


 DISCLAIMER
 -----------
 THIS MACRO IS PROVIDED "AS IS".  THERE ARE NO WARRANTIES, EXPRESSED 
 OR IMPLIED, AS TO MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE 
 REGARDING THE ACCURACY OF THE MATERIALS OR CODE CONTAINED HEREIN.
 -----------




-------------------------------------------------------------------------------------------------------
	Input
		dat=,			*** Name of input dataset. See documentation for data formatting.
		cond=,			*** Name of outcome measure for bernoulli portion of two-part model.
		condpred=,		*** Predictors used for the bernoulli piece, seperated by blanks. 
							Leaving this blank will only model the binary response 
							with intercept terms.
		rowvar=,		*** Name of the ordinal response for the row variable.
		colvar=,		*** Name of the ordinal response for the column variable.
		ct=,			*** Name of the cell count variable.
		rowpred=,		*** Predictors used for the row variable marginal model, seperated by blanks. 
							Leaving this blank will only model the Row response 
							with intercept terms.
		colpred=,		*** Predictors used for the column variable marginal model, seperated by blanks.
							Leaving this blank will only model the Column response 
							with intercept terms.
		gcrpred=,		*** Predictors used for the GCR model.
		gcrrow=,		*** 'row' Indicates the GCR model will contain row level terms. 
							Leaving this blank will omit the row level in the GCR model.
					
		gcrcol=			*** 'col' Indicates the GCR model will contain column level terms. 
							Leaving this blank will omit the column level in the GCR model.
IMPORTANT:
	-IT IS ADVISABLE TO LIMIT VARIABLE NAME LENGTHS TO 8 CHARACTERS OR LESS. 
		UNEXPECTED RESULTS MAY OCCUR OTHERWISE.
	-VARIABLE NAMES SHOULD NOT END WITH DIGITS. THIS CAUSES PROBLEM IN THE CODE.
	-FOR DRINKING PATTERN ANALYSIS (SUCH AS IN MCMILLAN ET AL.), REORDER THE QUANTITY 
		AND FREQUENCY VARIABLES SO THAT THE HIGHEST LEVEL OF DRINKING CORRESPONDS TO LEVEL 1,
		THE NEXT HIGHEST CORRESPONDS TO LEVEL 2, ETC.

Output datasets include:
BDM_BIN_EST = Covariance matrix for the Bernoulli part (if applicable).
BDM_BIN_PARMEST = Parameter estimates for the Bernoulli part (if applicable).
BDM_SPECS = Specification for the BDM.
BDM_FITS = Fit statistics for the BDM.
BDM_EST = Parameter estimates for the BDM.
BDM_COV = Covariance matrix for the BDM.
BDM_FINAL = Observed and predicted counts, with raw and Pearson residuals.

Notes: 
- The macro will automatically summarize data into categories by covariate combination 
 	 so that MLE procedure is improved. The data need not be in summarized form when submitted 
	 to the macro.
- The macro automatically assigns 0 starting values for all predictors in the GCR model
	other than row, column, or intercept effects. This is necessary to start the optimization 
	routine.

References:
McMillan, GP, Hanson, T, Bedrick, E and Lapham, SC (nd) "Application of the Bivariate Dale 
					Model to Estimating Risk Factors for Alcohol Consumption Patterns." 

Dale, JC (1985) "A Bivariate Discrete Model of Changing Colour in Blackbirds."
		in Statistics in Ornithology, Morgan, BJT and North, PM eds. Springer-Verlag, Berlin.

Dale, JC (1986) "Global Cross-Ratio Models for Bivariate, Discrete, Ordered Responses."
		Biometrics 42, 909-917.

Molenberghs, G and Lesaffre, E (1994) "Marginal Modeling of Correlated Ordinal Data Using 
					a Multivariate Plackett Distribution." Journal of the American Statistical 
					Association 89, 633-644.

*******************************************;


*OPTIONS MPRINT SYMBOLGEN;
*options nodate pageno = 1;

proc optsave;  RUN;

%macro BDM(dat,cond,condpred,rowvar,colvar,ct,rowpred,colpred,gcrpred,gcrrow,gcrcol);

proc datasets library = work;
	delete BDM_BIN_PARMEST BDM_SPECS BDM_FITS BDM_EST BDM_COV BDM_FINAL;
RUN; QUIT;


title1 ' ';
title2 ' ';
footnote1 ' ';
footnote2 ' ';
footnote3 ' ';

***** seperate bernoulli from bivariate ordinal data *******;
data BDM_bivdat BDM_bindat;
	set &dat;
		%if %length(&cond) %then %do;
			output BDM_bindat;
			if &cond then output BDM_bivdat;
		%end;
		%else %do; output BDM_bivdat; %end;
run;
*************************************************************;

******************************************************************************************;
****** Linear logistic model of the probability of having a binary ordinal response ******;
%if %length(&cond) %then %do;

ODS TRACE ON;
ODS OUTPUT PARAMETERESTIMATES = BDM_BIN_PARMEST;
	title1 "Logistic regression model of the probability of &cond. = 1.";
	proc logistic data = BDM_bindat descending outest = BDM_bin_est covout;
		model &cond = &condpred;
		WEIGHT &CT;
		output out = BDM_binpred p = pred;
	run;
ODS TRACE OFF;
%end;
******************************************************************************************;

***** begin by summarizing the data so that likelihood is easier to compute ***;
proc summary nway sum data = BDM_bivdat;
	var &ct;
	class &rowvar &colvar &rowpred &colpred &gcrpred &condpred;
	output out = BDM_data_out sum = &ct;
run;

proc sort data = BDM_data_out out = BDM_d1; by &rowvar. &colvar.; run;

***** determine levels of row variable and levels of column variable ****;
***** then create new variables with number of free parameters       ****;
data _null_;
	set BDM_d1 end = eof;
		if eof;
		rxc = &rowvar.*&colvar.;

		call symput("rxc",trim(left(rxc)));
		call symput("r",trim(left(&rowvar.)));
		call symput("c",trim(left(&colvar.)));
run;

%let rmin1 = %eval(&r. -1 );	*** The number of row intercepts ****;
%let cmin1 = %eval(&c. -1 );	*** The number of column intercepts ***;
%let rmin2 = %eval(&r. -2 );	*** The number of row terms in GCR model ***;
%let cmin2 = %eval(&c. -2 );	*** The number of column terms in GCR model ***;

*************************************************************************;

title1 ' ';
title2 ' ';


***** Macro %NVAR determines the number of variables in the list of predictors ****;
%macro nvar(varlist);
   %local vars count;
   %let count = 1;
   %let vars = %nrbquote(%scan(&varlist, &count));
   %do %until(&vars= );
      %let count = %eval(&count + 1); 
      %let vars = %nrbquote(%scan(&varlist, &count));
   %end;
   %eval(&count - 1)
%mend nvar;
*************************************************************************;


***** Macro %LINEAR writes out the linear predictor for each outcome measure and GCR model ****;
%macro linear(predlist,resp);
%if %length(&predlist) %then %do;
	%let nvars = %nvar(&predlist);


				%if &nvars =  1 %then %do;
					%let last = %scan( &predlist, 1 );
					%let list =  + &&resp._&last. * &last.;
				%end;

		%else %DO k = 1 %to %eval(&nvars-1);
			%let var&k. = %scan( &predlist, &k );
			%put &&var&k;	
			%let pred&k. = &&resp._&&var&k. ; 
			%put &&pred&k;  

			%if &k. = 1 %then %do; 
				%let list = &&resp._&&var&k. * &&var&k. +;
			%end;

			%else  %do;
				%let list = %cmpres(&list.)  &&resp._&&var&k. * &&var&k. +;
			%end;
		
		%end;

		%if &nvars ne 1 %then %do;
				%let last = %scan( &predlist, &nvars );
				%let list = + %cmpres(&list.)  &&resp._&last. * &last. ;
				%put &list.;
		%end;

		&list
%end;			
%mend linear;

%let rowpredict =  %linear(&rowpred.,&rowvar.); 	***** The linear predictor for the row outcome ****;
%let colpredict =  %linear(&colpred.,&colvar.);  ***** The linear predictor for the column outcome ****;
%let gcrpredict =  %linear(&GCRpred.,GCR); 		***** The linear predictor for the GCR model *****;
*************************************************************************;


****** starting values and parameters table for marginal models *********;
ods select parameterestimates;
ods output parameterestimates = BDM_&rowvar.init;
proc logistic data = BDM_d1 outest= BDM_&rowvar.fit;
	model &rowvar. = &rowpred.;
	weight &ct;
run;
ods output close;

ods select parameterestimates;
ods output parameterestimates = BDM_&colvar.init;
proc logistic data = BDM_d1 outest = BDM_&colvar.fit;
	model &colvar. = &colpred.;
	weight &ct;
run;
ods output close;

************* SUM LOG LIKELIHOOD FOR EACH MARGINAL ONLY MODEL ******************;
******* used to evaluate the drop in deviance when including the GCR portion ***;
data BDM_MARG_FIT;
	SET BDM_&rowvar.fit(keep = _lnlike_ rename = (_LNLIKE_ = ROWL)); 
	SET BDM_&colvar.fit(keep = _lnlike_ rename = (_LNLIKE_ = COLL));
		sumfit = -2*(ROWL+COLL);
		keep sumfit;
	label sumfit = "Sum -2 LogL for Marginal ONLY Models";
run;

****************** INITIAL VALUES FOR MLE ROUTINE **************;
data BDM_inits(keep = parameter estimate);
	length parameter $100.;

	set BDM_&colvar.init(keep = estimate variable classval0 in = &colvar.)
		BDM_&rowvar.init(keep = estimate variable classval0 in = &rowvar.);

		if variable = "Intercept" then do;
			if &rowvar. then parameter = "&rowvar." || "_int" || trim(left(classval0));
			if &colvar. then parameter = "&colvar." || "_int" || trim(left(classval0));
		end;
		else do;
			if &rowvar. then parameter = "&rowvar._" || trim(left(variable));
			if &colvar. then parameter = "&colvar._" || trim(left(variable));
		end;
run;
************************************************************;

****** starting values and parameters for GCR model ***********************;
**** Note that 1 is added to all cell counts so that the OR is estimable. *;
**** This is only done for the purposes of getting starting values. *******;

%macro gcrinits(gcrrow,gcrcol);

proc delete data = BDM_modeldata; run;

data BDM_modeldata;
	retain value . rowlevel . collevel . level .;
run;

%do i = 1 %to %eval(&r.-1);
	%do j = 1 %to %eval(&c.-1);
		data BDM_gcrdat;
			set BDM_d1;
				&ct = &ct + 1;
				if &rowvar. <= &i. then row = 1;
				if &rowvar. > &i. then row = 2;
				if &colvar. <= &j. then col = 1;
				if &colvar. > &j. then col = 2;
		run;

		ods trace on;
		ods select relativerisks;
		ods output relativerisks = BDM_orest;
			proc freq data = BDM_gcrdat ;
				tables row*col / plcorr;
				weight &ct;
			run;
		ods output close;
		ods trace off;

		data BDM_indat;	
			set BDM_orest;
			retain level 0;
				if index(studytype,"Odds") > 0;
				rowlevel = &i.;
				collevel = &j.;

				keep value rowlevel collevel level;
		run;

		proc append base = BDM_modeldata data = BDM_indat force; run;
	%end;

ods trace on;
ods output parameterestimates = BDM_pe;
proc genmod data = BDM_modeldata;
	class rowlevel collevel;
		model value = &gcrrow.level &gcrcol.level / link = log;
run;
ods output close;
ods trace off;

data BDM_gcrinits(where = (parameter ne ' ') keep =  estimate parameter);
length parameter $100.;
	set BDM_pe(rename = (parameter = parm));

		if upcase(parm) =  "ROWLEVEL" and input(level1,8.) <= %eval(&r.-2) then
			parameter = "GCR" || "&rowvar." || trim(left(level1)) ;
		if upcase(parm) = "COLLEVEL" and input(level1,8.) <= %eval(&c.-2) then
			parameter = "GCR" || "&colvar." || trim(left(level1)) ;
		if upcase(parm) = "INTERCEPT" then parameter = trim(left("GCR" || "DELTA"));
run;

%end;
	
%mend;
%gcrinits(&gcrrow.,&gcrcol.);

*****************starting values for other GCR parms ***********;

%macro otherparms(predlist,resp);
%local nvars k;
%if %length(&predlist) %then %do;
	%let nvars = %nvar(&predlist);
		%DO k = 1 %to %eval(&nvars);
			%let var&k. = %scan( &predlist, &k );
			%let pred&k. = &&resp._&&var&k. ; 
			%put &&pred&k;  
			parameter =  "&&pred&k";
			estimate = 0;
			output;
		%end; 
%end;
%mend otherparms;

data BDM_oinits;
length parameter $100;
	%otherparms(&gcrpred. , GCR); 
run;

************************************************************;

***** Concatenate table of starting values and parameters ****;
data BDM_parminits;
	set BDM_inits BDM_gcrinits %if %length(&gcrpred) %then %do; BDM_oinits; %end;
	;
run;
************************************************************;


**** Define an array for the row / column effects in the GCR model *****;
%macro arraylistrow;
%local list1 ;
	%let list1 =;	
	%if %length(&gcrrow) %then %do;			
		%let list1 = %str(array row_GCR{&r} GCR&rowvar.1 - GCR&rowvar.&rmin2. 0 0 ;) ;
	%end;
	&list1
%mend;

%macro arraylistcol;
%local list1 ;
	%let list1 =;	
	%if %length(&gcrcol) %then %do;
		%let list1 = %str(array col_GCR{&c} GCR&colvar.1 - GCR&colvar.&cmin2. 0 0;);
	%end;
	&list1
%mend;
***************************************************************************;

***************************************************************************;
**** Specify Row and Column terms in the GCR model ************************;
%macro gcrlistr;
	%local list;
		%if %length(&gcrrow) %then %do;
			%let list = + row_GCR[i];
		%end;
	&list
%mend;

%macro gcrlistc;
	%local list;
		%if %length(&gcrcol) %then %do;
			%let list = + col_GCR[j];
		%end;
	&list
%mend;	
******************************************************************;

******************************************************************;
***********PROC NLMIXED MAXIMIZES THE LIKELIHOOD FUNCTION ********;
***** See Dale (1986) for full development of the likelihood *****;
******************************************************************;

ods trace on;
ods output parameterestimates = BDM_EST
		   FITSTATISTICS = BDM_GOFIT
		   CONVERGENCESTATUS = BDM_CONV
		   specifications = BDM_SPECS
		   dimensions = BDM_DIM
		   CovMatParmEst=BDM_COV;

proc nlmixed data = BDM_D1 TECH = TRUREG OPTCHECK cov;

	parms  / data = BDM_parminits;

array cellprobs{&r.,&c.} 	cp1 - cp&rxc.;
array cumprobs{&r.,&c.}  	F1  - F&rxc.;
array S_terms{&r.,&c.} 		S1  - S&rxc.;
array log_psi{&r.,&c.} 		log_psi1 - log_psi&rxc.;
array psi{&r.,&c.} 			psi1 - psi&rxc.;
array row_margin{&r.} 		r_1 - r_&r.;
array col_margin{&c.} 		c_1 - c_&c.;

array row_ints{&rmin1.} &rowvar._int1 - &rowvar._int&rmin1.;
array col_ints{&cmin1.} &colvar._int1 - &colvar._int&cmin1.;

array eta{&rmin1.} eta1 - eta&rmin1.;
array ksi{&cmin1.} ksi1 - ksi&cmin1.;

%arraylistrow; 
%arraylistcol; 

do i = 1 to &rmin1.;
	eta[i] = row_ints[i]  &rowpredict.;
	row_margin[i] = exp(eta[i]) / (1+ exp(eta[i])); 
end;

do j = 1 to &cmin1.;
	ksi[j] = col_ints[j]  &colpredict.;
 	col_margin[j] = exp(ksi[j]) / (1+ exp(ksi[j])); 
end;

do i = 1 to &r.;
	do j = 1 to &c.;
	log_psi[i,j] = GCRdelta  &GCRpredict %gcrlistr %gcrlistc;
	end;
end;

row_margin[&r.] = 1;
col_margin[&c.] = 1;
cumprobs[&r.,&c.] = 1;

psi[1,1]       = exp(log_psi[1,1]); 
S_terms[1,1]   = sqrt(  (1+(r_1+c_1)*(psi[1,1]-1))**2 + 4*psi[1,1]*(1-psi[1,1])*r_1*c_1);
cumprobs[1,1]  = ( 1 +     (r_1+c_1)*(psi[1,1]-1) - S_terms[1,1] ) / (2*(psi[1,1]-1));
cellprobs[1,1] = cumprobs[1,1];
if &rowvar. = 1 and &colvar. = 1 then z = cellprobs[1,1];

do i = 2 to &r.;
	psi[i,1] = exp(log_psi[i,1]); 
	S_terms[i,1] = sqrt(  (1+(row_margin[i]+col_margin[1])*(psi[i,1]-1))**2 + 4*psi[i,1]*(1-psi[i,1])
							*row_margin[i]*col_margin[1]);	
	cumprobs[i,1] = ( 1 +    (row_margin[i]+col_margin[1])*(psi[i,1]-1) - S_terms[i,1] ) / (2*(psi[i,1]-1));
	cellprobs[i,1] = cumprobs[i,1] - cumprobs[i-1,1];
	if &rowvar. = i and &colvar. = 1 then z = cellprobs[i,1];
end;

do j = 2 to &c.;
	psi[1,j] = exp(log_psi[1,j]); 
	S_terms[1,j] = sqrt(  (1+(row_margin[1]+col_margin[j])*(psi[1,j]-1))**2 + 4*psi[1,j]*(1-psi[1,j])
							*row_margin[1]*col_margin[j]);	
	cumprobs[1,j] = ( 1 +    (row_margin[1]+col_margin[j])*(psi[1,j]-1) - S_terms[1,j] ) / (2*(psi[1,j]-1));
	cellprobs[1,j] = cumprobs[1,j] - cumprobs[1,j-1];
	if &rowvar. = 1 and &colvar. = j then z = cellprobs[1,j];
end;

do i = 2 to &r.;
  	do j = 2 to &c.;
		psi[i,j] = exp(log_psi[i,j]); 
		S_terms[i,j] = sqrt(  (1+(row_margin[i]+col_margin[j])*(psi[i,j]-1))**2 + 4*psi[i,j]*(1-psi[i,j])
							*row_margin[i]*col_margin[j]);	
		cumprobs[i,j] = ( 1 +    (row_margin[i]+col_margin[j])*(psi[i,j]-1) - S_terms[i,j] ) / (2*(psi[i,j]-1));
  		cellprobs[i,j] = cumprobs[i,j]-cumprobs[i-1,j]-cumprobs[i,j-1]+ cumprobs[i-1,j-1];
		if &rowvar. = i and &colvar. = j then z = cellprobs[i,j];
 	end;
end;

      if (z > 1e-8) then ll = &CT.*log(z);
		else ll = -1e100;
	 model ll ~ general(ll);

	**** if running SAS V9.0 or greater, use :
	**** predict z out = BDM_predicts;
run;
ods trace off;
ODS OUTPUT CLOSE;
******************************************************************;

**************** predicted values for SAS version 8.2 ***********;
%macro arraylistrow2;
%local list1 ;
	%let list1 =;	
	%if %length(&gcrrow) %then %do;			
		%let list1 = %str(array row_GCR{&r} GCR&rowvar.1 - GCR&rowvar.&rmin2. arlx0 arlx0 ;) ;
	%end;
	&list1
%mend;

%macro arraylistcol2;
%local list1 ;
	%let list1 =;
	%if %length(&gcrcol) %then %do;
		%let list1 = %str(array col_GCR{&c} GCR&colvar.1 - GCR&colvar.&cmin2. arlx0 arlx0;);
	%end;
	&list1
%mend;

proc transpose data = bdm_est out = BDM_parms;
var estimate;
id parameter;
run;

proc sql;
	create table BDM_pred_dat as
	select BDM_d1.*,
		   BDM_parms.*
	from BDM_d1, BDM_parms;
quit;

******** Generate output dataset with predicted counts ***************;
data BDM_predicts;
	retain arlx0 0;
	set BDM_pred_dat;

array cellprobs{&r.,&c.} 	cp1 - cp&rxc.;
array cumprobs{&r.,&c.}  	F1  - F&rxc.;
array S_terms{&r.,&c.} 		S1  - S&rxc.;
array log_psi{&r.,&c.} 		log_psi1 - log_psi&rxc.;
array psi{&r.,&c.} 			psi1 - psi&rxc.;
array row_margin{&r.} 		r_1 - r_&r.;
array col_margin{&c.} 		c_1 - c_&c.;

array row_ints{&rmin1.} &rowvar._int1 - &rowvar._int&rmin1.;
array col_ints{&cmin1.} &colvar._int1 - &colvar._int&cmin1.;

array eta{&rmin1.} eta1 - eta&rmin1.;
array ksi{&cmin1.} ksi1 - ksi&cmin1.;

%arraylistrow2; 
%arraylistcol2; 

do i = 1 to &rmin1.;
	eta[i] = row_ints[i]  &rowpredict.;
	row_margin[i] = exp(eta[i]) / (1+ exp(eta[i])); 
end;

do j = 1 to &cmin1.;
	ksi[j] = col_ints[j]  &colpredict.;
 	col_margin[j] = exp(ksi[j]) / (1+ exp(ksi[j])); 
end;

do i = 1 to &r.;
	do j = 1 to &c.;
	log_psi[i,j] = GCRdelta  &GCRpredict %gcrlistr %gcrlistc;
	end;
end;

row_margin[&r.] = 1;
col_margin[&c.] = 1;
cumprobs[&r.,&c.] = 1;

psi[1,1]       = exp(log_psi[1,1]); 
S_terms[1,1]   = sqrt(  (1+(r_1+c_1)*(psi[1,1]-1))**2 + 4*psi[1,1]*(1-psi[1,1])*r_1*c_1);
cumprobs[1,1]  = ( 1 +     (r_1+c_1)*(psi[1,1]-1) - S_terms[1,1] ) / (2*(psi[1,1]-1));
cellprobs[1,1] = cumprobs[1,1];
if &rowvar. = 1 and &colvar. = 1 then z = cellprobs[1,1];

do i = 2 to &r.;
	psi[i,1] = exp(log_psi[i,1]); 
	S_terms[i,1] = sqrt(  (1+(row_margin[i]+col_margin[1])*(psi[i,1]-1))**2 + 4*psi[i,1]*(1-psi[i,1])
							*row_margin[i]*col_margin[1]);	
	cumprobs[i,1] = ( 1 +    (row_margin[i]+col_margin[1])*(psi[i,1]-1) - S_terms[i,1] ) / (2*(psi[i,1]-1));
	cellprobs[i,1] = cumprobs[i,1] - cumprobs[i-1,1];
	if &rowvar. = i and &colvar. = 1 then z = cellprobs[i,1];
end;

do j = 2 to &c.;
	psi[1,j] = exp(log_psi[1,j]); 
	S_terms[1,j] = sqrt(  (1+(row_margin[1]+col_margin[j])*(psi[1,j]-1))**2 + 4*psi[1,j]*(1-psi[1,j])
							*row_margin[1]*col_margin[j]);	
	cumprobs[1,j] = ( 1 +    (row_margin[1]+col_margin[j])*(psi[1,j]-1) - S_terms[1,j] ) / (2*(psi[1,j]-1));
	cellprobs[1,j] = cumprobs[1,j] - cumprobs[1,j-1];
	if &rowvar. = 1 and &colvar. = j then z = cellprobs[1,j];
end;

do i = 2 to &r.;
  	do j = 2 to &c.;
		psi[i,j] = exp(log_psi[i,j]); 
		S_terms[i,j] = sqrt(  (1+(row_margin[i]+col_margin[j])*(psi[i,j]-1))**2 + 4*psi[i,j]*(1-psi[i,j])
							*row_margin[i]*col_margin[j]);	
		cumprobs[i,j] = ( 1 +    (row_margin[i]+col_margin[j])*(psi[i,j]-1) - S_terms[i,j] ) / (2*(psi[i,j]-1));
  		cellprobs[i,j] = cumprobs[i,j]-cumprobs[i-1,j]-cumprobs[i,j-1]+ cumprobs[i-1,j-1];
		if &rowvar. = i and &colvar. = j then z = cellprobs[i,j];
 	end;
end;
	label z = "Pred. Prob. &Rowvar <= i AND &colvar <= j";
run;


************* OBSERVED COUNTS OVER ALL OBSERVED COVARIATE COMBO LEVELS ***********;
data BDM_newdat;
	set &dat;
		%if %length(&cond) %then %do;
			abs = &ct. * (1-&cond);
			&ct = &ct * &cond;
		%end;
run;

proc means data = BDM_newdat sum nway noprint;
	var &ct. %if %length(&cond) %then %do; abs %end;
	;			*** count variable and n abstainers = &ct. * (1 - &cond )***;
	class &CONDPRED &rowpred &colpred &gcrpred;
	output out = BDM_counts sum = n_pos %if %length(&cond) %then %do; n_abst %end;
	;
run;

	

%if %length(&rowpred) or %length(&colpred) or %length(&gcrpred) or %length(&condpred) %then %do;

proc sort data = BDM_counts;   by &condpred &rowpred &colpred &gcrpred; run;
proc sort data = BDM_predicts; by &condpred &rowpred &colpred &gcrpred; run;
	%if %length(&cond) %then %do; 
proc sort data = BDM_binpred out = BDM_bout nodupkey; by &condpred &rowpred &colpred &gcrpred; run;
	%end;

data BDM_final;
	merge BDM_predicts 
		  BDM_counts 
		    %if %length(&cond) %then %do; 
		  BDM_bout(keep = &condpred &rowpred &colpred &gcrpred pred);
		    %end;
							;
		by &condpred &rowpred &colpred &gcrpred;
run;
%end;

%else %do;
data BDM_bout;
	set BDM_binpred;
		if _n_= 1 then output;
run;

				proc sql;
					create table BDM_final as
					select BDM_counts.*,
		   					BDM_predicts.*
							%if %length(&cond) %then %do;
									,BDM_bout.pred
							%end;
					from BDM_counts, BDM_predicts
							%if %length(&cond) %then %do;
									,BDM_bout
							%end;
							;
				quit;
%end;

data BDM_final;
	set BDM_final;
		pred_count = n_pos*z %if %length(&cond) %then %do;  %end;	;
			%if %length(&cond) %then %do;  
				if _freq_ = (n_pos + n_abst) then do;
					PRED_ABSTAINER = (1-pred)* _freq_;
				end;
				else do;
					n_abst = .;
				end;
			%end;
		raw_residual = &ct. - pred_count;
		std_res = (&ct. - pred_count) / (sqrt(pred_count*(1-z)));
		label pred_count = "Predicted count"
			  raw_residual = "Raw Residual (O-E)"
			  std_res = "Pearson Residual"
			  pred = "Predicted Prob. of &cond."
			  n_abst = "Observed Abstainers"
			  pred_abstainer = "Predicted # Abstainers"
			 ;
run;


******** BEGIN GENERATE LIST OF UNIQUE PREDICTORS FOR PRINTOUT **************;
data _null_;
	x = "%UPCASE(&condpred) %UPCASE(&rowpred) %UPCASE(&colpred) %UPCASE(&gcrpred)";
	x = trim(left(x));
	length new_x $ 200 word $50;
	DO UNTIL (WORD = " ");
		 i = sum(i,1);
		 WORD = TRIM(LEFT(scan(X,i," ") ) );
		 if indexw(new_x,word) = 0 then new_x =  left ( trim(new_x) || " " || trim(word) );
	end; 	 

	call symput("PRINTLIST",new_x);
run;

******** END GENERATE LIST OF UNIQUE PREDICTORS FOR PRINTOUT **************;


**************** PRINT RESULTS ****************;
options NODATE formdlim = " " pageno = 1 pagesize = min;
DM output 'clear';   **** CLEARS ALL OUTPUT UP TO THIS POINT ****;

%if %length(&cond) %then %do;
	title1 "Logistic regression model parameter estimates of the probability of &cond..";
	proc print data = BDM_bin_PARMest noobs; run;
%end;

proc print data = bdm_specs noobs;
	title1 "Specifications for Bivariate Dale Model";
	title2 "Response variables &rowvar &colvar";
run;

proc print data = bdm_conv noobs;
	title1 "Convergence Status of Bivariate Dale Model";
	title2 " ";
	var reason;
run;

data _null_;
	set bdm_est(where = ( substr(parameter,1,3) = "GCR") ) end=eof;
		k+1;
			if eof then call symput("DF",k);	
run;


data BDM_fits;
	set bdm_gofit; 
	set BDM_marg_fit;
		dd = sumfit - value;
		DF = &DF.;
			p=1-probchi(dd,df);
			format p pvalue6.4;
run;

proc print data = BDM_fits noobs label;
	title1 "Fit Statistics for Bivariate Dale Model and Marginal Models";
	label value = "Value for BDM"
			dd = "Drop in Deviance w/GCR"
			p = "P-value";
run;


data _null_;
	set bdm_est end = eof;
		if eof then do;
			call symput("pl2",_n_);
		end;
run;

options pagesize = %eval(&pl2+10);
proc print data = bdm_est noobs label;
	var parameter estimate standarderror probt lower upper;
	label estimate = 'Log Odds Ratio'
			lower = 'Lower CL'
			upper = 'Upper CL'
			probt = 'P-value';
	title1 "Parameter Estimates for Bivariate Dale Model";
FOOTNOTE1 "Note that regression coefficients for the marginal model";
FOOTNOTE2 "  indicate the log odds ratio of being at or below a particular response level.";
FOOTNOTE3 "  Negative coefficients indicate reduced probability of being at or BELOW a particular level.";
run;


data tabledat;
	set bdm_est;
		OR = exp(estimate);
		UCL = exp(estimate + 1.96*Standarderror);
		LCL = exp(estimate - 1.96*Standarderror);
		format or lcl ucl 5.2;
run;

proc print data = tabledat noobs;
title1 "Odds Ratios and Confidence Intervals";
var parameter or lcl ucl;
run;

options nocenter pagesize = 1000;
proc print data = BDM_final noobs label;
title "Observed & Predicted Values from Bivariate Dale Model";
title2 "Dataset = work.final";
	var &rowvar &colvar &PRINTLIST pred_count &ct. RAW_RESIDUAL std_res
			%if %length(&cond) %then %do;  
			PRED_ABSTAINER n_abst
			%end;
			;
	sum &ct.;
	footnote1 ;
run;

proc optload; run;

***** cleanup *********;
PROC DATASETS LIBRARY = WORK;
	DELETE BDM_BINDAT BDM_BIVDAT BDM_COUNTS BDM_D1 BDM_DATA_OUT BDM_BOUT BDM_BINPRED
			BDM_DIM BDM_GCRDAT BDM_GCRINITS BDM_GOFIT BDM_&ROWVAR.FIT 
			BDM_&ROWVAR.INIT BDM_&COLVAR.FIT BDM_&COLVAR.INIT
			BDM_INDAT BDM_INITS BDM_MARG_FIT BDM_MODELDATA
			BDM_OINITS BDM_OREST BDM_PARMINITS
			BDM_PARMS BDM_PE BDM_PREDICTS BDM_PRED_DAT BDM_NEWDAT BDM_CONV;
RUN; QUIT;

***********************************************;
%mend BDM;


/*
***** 	Example: Dale (1986), Table 3, P.913 ***********;

Data dale1986;
	input OPERATION $ pain NEVER SELDOM OCC REG;

IF OPERATION = 'VP' THEN VP = 1; ELSE VP = 0;
IF OPERATION = 'VA' THEN VA = 1; ELSE VA = 0;
IF OPERATION = 'VH' THEN VH = 1; ELSE VH = 0;
IF OPERATION = 'RA' THEN RA = 1; ELSE RA = 0;

	MED = 1; CT = NEVER;   OUTPUT; 
	MED = 2; CT = SELDOM;  OUTPUT; 
	MED = 3; CT = OCC; 	   OUTPUT; 
	MED = 4; CT = REG;     OUTPUT; 

DROP  NEVER SELDOM OCC REG operation;
DATALINES;
VP	1	170	7	8	0
VP	2	18	5	8	3
VP	3	7	0	4	14
VA	1	170	7	5	2
VA	2	22	7	8	1
VA	3	8	1	8	9
VH	1	176	8	5	2
VH	2	26	6	5	5
VH	3	14	3	2	9
RA	1	181	6	6	2
RA	2	17	12	7	3
RA	3	10	2	3	11
;
RUN;

proc print data = dale1986 noobs; run;

%BDM(dat=dale1986,
		cond=,
		condpred=, 
		rowvar=PAIN,
		colvar=MED,
		ct=ct,
		rowpred=,
		colpred=,
		gcrpred=VH,
		gcrrow=,
		gcrcol=col);
*/


********* Example: Dale 1985, Table 2: model a, page 32*******;
/*
Data dale1985;
	input period orbit MAN0 MAN1 MAN2 MAN3 MAN4;

ORBIT = ORBIT +1;

IF PERIOD IN (1,2,3,7,10,15,17,18) THEN DELETE;

	MANDIBLE = 1; CT = MAN0;   OUTPUT; 
	MANDIBLE = 2; CT = MAN1+MAN2;   OUTPUT; 
	MANDIBLE = 3; CT = MAN3+MAN4;   OUTPUT; 

DROP MAN0 MAN1 MAN2 MAN3 MAN4 ;
datalines;
1	0	1	0	0	0	0
1	1	0	0	0	0	0
1	2	0	0	0	0	0
2	0	2	0	0	0	0
2	1	0	0	0	0	0
2	2	0	0	0	0	0
3	0	5	0	0	0	0
3	1	0	0	0	0	0
3	2	0	0	0	0	0
4	0	9	2	1	0	0
4	1	3	1	0	0	0
4	2	0	0	0	0	0
5	0	7	2	0	0	0
5	1	0	0	0	1	0
5	2	0	0	0	0	0
6	0	4	0	0	0	0
6	1	1	2	0	0	0
6	2	0	0	0	0	0
7	0	2	0	0	0	0
7	1	3	0	0	0	0
7	2	0	0	0	0	0
8	0	6	1	0	0	0
8	1	2	2	2	2	0
8	2	0	0	0	0	0
10	0	0	0	0	0	0
10	1	0	2	0	1	0
10	2	0	0	0	0	0
11	0	0	0	0	0	0
11	1	3	0	0	0	0
11	2	0	0	1	0	0
12	0	0	0	0	0	0
12	1	1	2	0	0	0
12	2	0	0	0	1	0
14	0	1	0	0	0	0
14	1	1	2	1	2	0
14	2	0	0	1	2	4
15	0	1	0	0	0	0
15	1	1	1	0	0	0
15	2	0	0	0	0	0
17	0	0	0	0	0	0
17	1	0	0	0	1	0
17	2	0	0	0	0	0
18	0	0	0	0	0	0
18	1	0	0	0	1	0
18	2	0	0	0	0	0
;
run;

PROC PRINT DATA = dale1985 NOOBS;
RUN;


PROC summary DATA = dale1985 nway sum;
	var ct;
	class ORBIT MANDIBLE PERIOD;
	OUTPUT OUT = RAWDATA SUM = CT;
RUN;

%BDM(dat=RAWDATA,
	 rowvar=ORBIT,
	 colvar=MANDIBLE,
	 ct=ct,
	 rowpred=PERIOD,
	 colpred=PERIOD,
	 gcrpred=,
	 gcrrow=,
	 gcrcol=);

*/