***************************************************************************; ** PROGRAM: Gange_3MRCV.sas **; ** AUTHOR: Christopher R. Bilder **; ** Department of Statistics **; ** University of Nebraska-Lincoln **; ** chris@chrisbilder.com **; ** DATE: 5-2-05 **; ** PURPOSE: Finds a vector of multinomial probabilities subject to the **; ** constraints of the Ho model. The fit of the Ho model **; ** provides configuations for the iterative proportional fitting**; ** algorithm. Use this program with 3 MRCVs only. **; ** NOTES: **; ** 1) Do not distribute copies of this program without the permission of**; ** the author. **; ** 2) Copyright 2005 Christopher R. Bilder **; ** 3) See control_3MRCV.sas for more information **; ***************************************************************************; **********************************************************************; ** **; ** NAME: MACRO P_MAT **; ** PURPOSE: Initial multinomial probabilities **; ** **; **********************************************************************; %MACRO P_MAT; %let X = do item1 = 0 to 1%str(;); %let endit = end%str(;); *Generalize for number of items; %do j = 2 %to &ncol; %let X = &X.do item&j=0 to 1%str(;); %let endit = &endit end%str(;); %end; *Create initial matrix of all items and their corresponding multinomial probabilities; data set1_init; &X; p=0.5; output; &endit; run; %MEND P_MAT; **********************************************************************; ** **; ** NAME: MACRO ADJUST_3WAY **; ** PURPOSE: Adjust the multinomial probabilities using IPF **; ** NOTES: For the predicted across MRCV items **; ** **; **********************************************************************; %MACRO ADJUST_3WAY; %do i = 1 %to &nrow2; %do j = (&nrow2 + 1) %to (&nrow2 + &ncol2); %do k = (&nrow2 + &ncol2 + 1) %to (&nrow2 + &ncol2 + &nstrat2); *Find the sub-table; obs_pred_small = obs_pred[loc(obs_pred[,1] = &i & obs_pred[,2] = &j & obs_pred[,3] = &k),]; *print obs_pred_small; *find the 3-way table for current multinomial probabilities; three_p_small = repeat(0, 2**3, 7); counter = 1; do a = 0 to 1; do b = 0 to 1; do c = 0 to 1; marg_sum = sum( p_matrix[ (loc(p_matrix[,&i] = a & p_matrix[,&j] = b & p_matrix[,&k] = c)), &ncol+1]); three_p_small[counter,] = &i || &j || &k || a || b || c || marg_sum; counter = counter + 1; end; end; end; *print three_p_small; *Create IPF adjustment factor; * Notice use of n here - may be better to adjust for adding small constants to the * 0 sub-table cells; adj000 = (obs_pred_small[1,7]/&n) / three_p_small[1,7]; adj001 = (obs_pred_small[2,7]/&n) / three_p_small[2,7]; adj010 = (obs_pred_small[3,7]/&n) / three_p_small[3,7]; adj011 = (obs_pred_small[4,7]/&n) / three_p_small[4,7]; adj100 = (obs_pred_small[5,7]/&n) / three_p_small[5,7]; adj101 = (obs_pred_small[6,7]/&n) / three_p_small[6,7]; adj110 = (obs_pred_small[7,7]/&n) / three_p_small[7,7]; adj111 = (obs_pred_small[8,7]/&n) / three_p_small[8,7]; *print adj000 adj001 adj010 adj011 adj100 adj101 adj110 adj111; *Adjust the multinomial probabilities; do i = 1 to nrow(p_matrix); if p_matrix[i, &i]=0 & p_matrix[i, &j]=0 & p_matrix[i, &k]=0 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj000; if p_matrix[i, &i]=0 & p_matrix[i, &j]=0 & p_matrix[i, &k]=1 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj001; if p_matrix[i, &i]=0 & p_matrix[i, &j]=1 & p_matrix[i, &k]=0 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj010; if p_matrix[i, &i]=0 & p_matrix[i, &j]=1 & p_matrix[i, &k]=1 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj011; if p_matrix[i, &i]=1 & p_matrix[i, &j]=0 & p_matrix[i, &k]=0 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj100; if p_matrix[i, &i]=1 & p_matrix[i, &j]=0 & p_matrix[i, &k]=1 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj101; if p_matrix[i, &i]=1 & p_matrix[i, &j]=1 & p_matrix[i, &k]=0 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj110; if p_matrix[i, &i]=1 & p_matrix[i, &j]=1 & p_matrix[i, &k]=1 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj111; end; %end; %end; %end; %MEND ADJUST_3WAY; **********************************************************************; ** **; ** NAME: MACRO ADJUST_2WAY **; ** PURPOSE: Adjust the multinomial probabilities using IPF **; ** NOTES: For the observed within MRCV items **; ** **; **********************************************************************; %MACRO ADJUST_2WAY; *Adjust for W pairs; %do i = 1 %to (&nrow2-1); %do j = (&i + 1) %to &nrow2; obs_pred_small = obs_2way[loc(obs_2way[,1] = &i & obs_2way[,2] = &j),]; *print obs_pred_small; *find the 2-way table for current multinomial probabilities; two_p_small = repeat(0, 2**2, 5); counter = 1; do a = 0 to 1; do b = 0 to 1; marg_sum = sum(p_matrix[(loc(p_matrix[,&i] = a & p_matrix[,&j] = b)), &ncol+1]); two_p_small[counter,] = &i || &j || a || b || marg_sum; counter = counter + 1; end; end; *print two_p_small; *Create IPF adjustment factor; * Notice use of n here - may be better to adjust for adding small constants to the * 0 sub-table cells; adj00 = (obs_pred_small[1,5]/&n) / two_p_small[1,5]; adj01 = (obs_pred_small[2,5]/&n) / two_p_small[2,5]; adj10 = (obs_pred_small[3,5]/&n) / two_p_small[3,5]; adj11 = (obs_pred_small[4,5]/&n) / two_p_small[4,5]; *print adj00 adj01 adj10 adj11; *Adjust the multinomial probabilities; do i = 1 to nrow(p_matrix); if p_matrix[i, &i]=0 & p_matrix[i, &j]=0 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj00; if p_matrix[i, &i]=0 & p_matrix[i, &j]=1 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj01; if p_matrix[i, &i]=1 & p_matrix[i, &j]=0 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj10; if p_matrix[i, &i]=1 & p_matrix[i, &j]=1 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj11; end; %end; %end; *Adjust for Y pairs; %do i = (&nrow2+1) %to (&nrow2+&ncol2-1); %do j = (&i + 1) %to (&nrow2+&ncol2); obs_pred_small = obs_2way[loc(obs_2way[,1] = &i & obs_2way[,2] = &j),]; *print obs_pred_small; *find the 2-way table for current joint probabilities; two_p_small = repeat(0, 2**2, 5); counter = 1; do a = 0 to 1; do b = 0 to 1; marg_sum = sum(p_matrix[(loc(p_matrix[,&i] = a & p_matrix[,&j] = b)), &ncol+1]); two_p_small[counter,] = &i || &j || a || b || marg_sum; counter = counter + 1; end; end; *print two_p_small; *Create IPF adjustment factor; adj00 = (obs_pred_small[1,5]/&n) / two_p_small[1,5]; adj01 = (obs_pred_small[2,5]/&n) / two_p_small[2,5]; adj10 = (obs_pred_small[3,5]/&n) / two_p_small[3,5]; adj11 = (obs_pred_small[4,5]/&n) / two_p_small[4,5]; *print adj00 adj01 adj10 adj11; *Adjust the multinomial probabilities; do i = 1 to nrow(p_matrix); if p_matrix[i, &i]=0 & p_matrix[i, &j]=0 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj00; if p_matrix[i, &i]=0 & p_matrix[i, &j]=1 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj01; if p_matrix[i, &i]=1 & p_matrix[i, &j]=0 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj10; if p_matrix[i, &i]=1 & p_matrix[i, &j]=1 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj11; end; %end; %end; *Adjust for Z pairs; %do i = (&nrow2+&ncol2+1) %to (&nrow2+&ncol2+&nstrat2-1); %do j = (&i + 1) %to (&nrow2+&ncol2+&nstrat2); obs_pred_small = obs_2way[loc(obs_2way[,1] = &i & obs_2way[,2] = &j),]; *print obs_pred_small; *find the 2-way table for current joint probabilities; two_p_small = repeat(0, 2**2, 5); counter = 1; do a = 0 to 1; do b = 0 to 1; marg_sum = sum(p_matrix[(loc(p_matrix[,&i] = a & p_matrix[,&j] = b)), &ncol+1]); two_p_small[counter,] = &i || &j || a || b || marg_sum; counter = counter + 1; end; end; *print two_p_small; *Create IPF adjustment factor; adj00 = (obs_pred_small[1,5]/&n) / two_p_small[1,5]; adj01 = (obs_pred_small[2,5]/&n) / two_p_small[2,5]; adj10 = (obs_pred_small[3,5]/&n) / two_p_small[3,5]; adj11 = (obs_pred_small[4,5]/&n) / two_p_small[4,5]; *print adj00 adj01 adj10 adj11; *Adjust the multinomial probabilities; do i = 1 to nrow(p_matrix); if p_matrix[i, &i]=0 & p_matrix[i, &j]=0 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj00; if p_matrix[i, &i]=0 & p_matrix[i, &j]=1 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj01; if p_matrix[i, &i]=1 & p_matrix[i, &j]=0 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj10; if p_matrix[i, &i]=1 & p_matrix[i, &j]=1 then p_matrix[i,&ncol+1] = p_matrix[i,&ncol+1] * adj11; end; %end; %end; %MEND ADJUST_2WAY; **********************************************************************; ** **; ** NAME: MACRO CHECK_TABLES **; ** PURPOSE: Find multinomial probabilities **; ** **; **********************************************************************; %MACRO CHECK_TABLES; %do i = 1 %to &nrow2; %do j = (&nrow2 + 1) %to (&nrow2 + &ncol2); %do k = (&nrow2 + &ncol2 + 1) %to (&nrow2 + &ncol2 + &nstrat2); *Find the sub-table; obs_pred_small = obs_pred[loc(obs_pred[,1] = &i & obs_pred[,2] = &j & obs_pred[,3] = &k),]; *find the 3-way table for current multinomial probabilities; three_p_small = repeat(0, 2**3, 7); counter = 1; do a = 0 to 1; do b = 0 to 1; do c = 0 to 1; marg_sum = sum( p_matrix[ (loc(p_matrix[,&i] = a & p_matrix[,&j] = b & p_matrix[,&k] = c)), &ncol+1]); three_p_small[counter,] = &i || &j || &k || a || b || c || marg_sum; counter = counter + 1; end; end; end; obs_pred_small_divn = obs_pred_small[,7]/&n; *print three_p_small obs_pred_small_divn; %end; %end; %end; %MEND CHECK_TABLES; **********************************************************************; ** **; ** NAME: MACRO DO_IT_GANGE **; ** PURPOSE: Find multinomial probabilities to use for resampling **; ** **; **********************************************************************; %MACRO DO_IT_GANGE; %P_MAT; proc iml; edit set1_init; read all into p_matrix; *print p_matrix; edit pred_sort2; read all into obs_pred; *print obs_pred; edit out_2way_save; read all into obs_2way; *print obs_2way; *Set initial values; tol = &tol; save_p=1; iter = 1; *Loop to do iterative proportional fitting (IPF); do loopit=1 to &gange_iter while(max(abs(save_p-p_matrix[,&ncol+1]))>tol); *Used to determine convergence; save_p = p_matrix[,&ncol+1]; %ADJUST_3WAY; %ADJUST_2WAY; error = max(abs(save_p-p_matrix[,&ncol+1])); *print iter error; iter = iter + 1; end; %CHECK_TABLES; p = T(p_matrix[,&ncol+1]); *print 'Final joint probabilities:' p; *print 'Iter:' iter; *Find ending error; error = max(abs(save_p-T(p))); *print error; create save_p from p_matrix; append from p_matrix; *save information about the convergence; info = &sim_numb || error || iter; col = {'sim_numb' 'error' 'iter_conv'}; create info_set from info[colname=col]; append from info; quit; *Just rename of data set - could use this data step to save out to a file; data set1_tau; set save_p; run; *Save information about the convergence; data save_info_set; set save_info_set info_set; run; %MEND DO_IT_GANGE;