HELP: How to evenly distribution of product ID and other parameters between two groups

[vicky07]

Hello,

I’ve a challenge where I need to do even distribution of multiple parameters with in two groups. In my below data I have 23 unique patient ID groups. Each patient ID group has multiple records with varying balance, risk and prod type’s combinations.

I want to split the total unique patient id in to two groups. Group A is 80% and group B is 20%. So in this case Group A should be 18 unique patient ID and group B should be 5 patient group. This split percent could change every month. Next run  could be 70/30 split.

What I am trying to achieve is similar mean balance, similar mean risk and uniform prod split between test and control groups. I was told that SAS/OR might solve this problem. I’ve SAS 9.4  Level 1M3 version.

Desired output:

Test						Control
Mean Bal	Mean Risk	AH	KO	ZU		Mean Bal	Mean Risk	AH	KO	ZU
2700	0.676	38%	22%	40%		2650	0.68	38%	22%	40%

Any help is much appreciated. Thank you in advance for your time and attention!!

PS: I have close to 480 unique patient ID groups but I am sharing only few for sample purposes. So please let me know if you need more data.  

My Data:

Patient Group	Balance	Risk	prod
3116290	4842.86	0.2222	ZU
3116290	1423.83	0.0200	AH
3116290	10040.3	0.0857	AH
3116290	3081.93	0.1514	ZU
3116290	4393.63	0.2171	AH
3116290	12199.3	0.2829	ZU
3116290	600.71	0.3486	KO
3116290	5082.37	0.4143	AH
3116290	13443.3	0.4800	KO
3116290	8202.96	0.5457	ZU
3116290	4601	0.6114	KO
3480975	5169.2	0.8889	ZU
3480975	6011.33	0.0429	AH
3480975	1832.1	0.1086	KO
3480975	2859.38	0.1743	KO
3480975	8462.93	0.2400	ZU
3480975	8318.86	0.3057	AH
3480975	7715.96	0.3714	ZU
3480975	2522.06	0.4371	AH
3480975	1637.89	0.5029	ZU
3480975	1942.42	0.5686	KO
3480975	2478.48	0.6343	ZU
3954070	14022.8	0.6667	AH
3954070	4831.72	0.6667	ZU
3954070	1000.32	0.0543	AH
3954070	2466.47	0.1200	AH
3954070	7685.72	0.1857	ZU
3954070	2011.07	0.2514	KO
3954070	9222.99	0.3171	ZU
3954070	5234.82	0.3829	ZU
3954070	7727.08	0.4486	AH
3954070	9481.1	0.5143	KO
3954070	6005.09	0.5800	AH
3954070	1738.4	0.6457	KO
4183977	5520.07	0.2222	AH
4183977	4995.73	0.7778	KO
4183977	10158	0.0629	ZU
4183977	7960.01	0.1286	KO
4183977	4713.21	0.1943	ZU
4183977	10567.3	0.2600	KO
4183977	9962.74	0.3257	ZU
4183977	2846.64	0.3914	KO
4183977	4343.59	0.4571	AH
4183977	5405.71	0.5229	ZU
4183977	10297.6	0.5886	AH
5692230	6074.19	0.6667	KO
5692230	5847.36	0.0114	ZU
5692230	3724.57	0.0771	AH
5692230	5670	0.1429	ZU
5692230	6360.26	0.2086	KO
5692230	3020.87	0.2743	AH
5692230	773.24	0.3400	AH
5692230	3565.45	0.4057	ZU
5692230	1610.1	0.4714	KO
5692230	7046.04	0.5371	ZU
5692230	2315.76	0.6029	KO
5852387	2407.18	0.7778	ZU
5852387	16275	0.0343	AH
5852387	5783.73	0.1000	KO
5852387	7513.21	0.1657	AH
5852387	1604.84	0.2314	KO
5852387	5114.36	0.2971	AH
5852387	4473.8	0.3629	ZU
5852387	7395.58	0.4286	AH
5852387	13924	0.4943	ZU
5852387	8483.98	0.5600	KO
5852387	14856	0.6257	AH
5981283	3942.05	0.1111	AH
5981283	8109.96	0.0314	KO
5981283	2604.99	0.0971	ZU
5981283	1627.6	0.1629	KO
5981283	10124.5	0.2286	ZU
5981283	4503.83	0.2943	KO
5981283	2077.9	0.3600	AH
5981283	4858.14	0.4257	KO
5981283	15987.3	0.4914	KO
5981283	5338.41	0.5571	ZU
5981283	5411.3	0.6229	AH
6681096	2671.5	0.4444	ZU
6681096	2014.46	0.8889	AH
6681096	10174.2	0.0657	KO
6681096	7755.11	0.1314	AH
6681096	9176.04	0.1971	KO
6681096	6199.91	0.2629	AH
6681096	5959.14	0.3286	KO
6681096	4568.93	0.3943	AH
6681096	4901.56	0.4600	ZU
6681096	7179.79	0.5257	KO
6681096	4385.1	0.5914	AH
6712099	4458.79	0.8889	KO
6712099	1081.77	0.0286	ZU
6712099	1882.14	0.0943	KO
6712099	4956.5	0.1600	ZU
6712099	1923.54	0.2257	AH
6712099	5555.1	0.2914	ZU
6712099	5857.5	0.3571	KO
6712099	951.77	0.4229	ZU
6712099	4765.94	0.4886	ZU
6712099	3899.5	0.5543	AH
6712099	1566.08	0.6200	KO
6915827	10802.4	0.1111	AH
6915827	6550.25	0.3333	KO
6915827	6538.26	0.0457	ZU
6915827	9412.25	0.1114	AH
6915827	11126.8	0.1771	ZU
6915827	3565.03	0.2429	KO
6915827	5312.21	0.3086	AH
6915827	2552.77	0.3743	KO
6915827	9565.37	0.4400	ZU
6915827	3142.98	0.5057	KO
6915827	2162.13	0.5714	AH
6915827	5518.75	0.6371	KO
7081267	1990.64	0.1111	KO
7081267	2306.08	0.0229	ZU
7081267	1774.27	0.0886	AH
7081267	5419.36	0.1543	KO
7081267	5446.57	0.2200	ZU
7081267	5760.97	0.2857	KO
7081267	8158.38	0.3514	AH
7081267	3845.6	0.4171	ZU
7081267	7852.74	0.4829	AH
7081267	9937.26	0.5486	KO
7081267	9644.65	0.6143	AH
7712089	8489.78	0.4444	AH
7712089	11872.2	0.0143	KO
7712089	6018.03	0.0800	ZU
7712089	2368.36	0.1457	KO
7712089	7418.02	0.2114	ZU
7712089	1725.35	0.2771	ZU
7712089	5759.08	0.3429	AH
7712089	5251.2	0.4086	KO
7712089	5822.4	0.4743	AH
7712089	5829.46	0.5400	KO
7712089	3064.08	0.6057	AH
7823183	1314.14	0.1111	KO
7823183	10605.1	0.5556	ZU
7823183	734.27	0.0600	KO
7823183	3428.32	0.1257	ZU
7823183	8687.86	0.1914	AH
7823183	6893.45	0.2571	ZU
7823183	5771.71	0.3229	AH
7823183	7104.48	0.3886	ZU
7823183	1906.69	0.4543	KO
7823183	7374.7	0.5200	AH
7823183	2243.76	0.5857	KO
7854689	5999.57	0.0235	ZU
7854689	6328.44	0.0086	AH
7854689	2025.56	0.0743	KO
7854689	9011.13	0.1400	AH
7854689	2396.85	0.2057	ZU
7854689	1499.99	0.2714	AH
7854689	3843.56	0.3371	KO
7854689	1807.62	0.4029	AH
7854689	1381.71	0.4686	ZU
7854689	2279.52	0.5343	AH
7854689	3627.87	0.6000	ZU
7893195	4203.73	0.3333	ZU
7893195	6851.77	0.0257	KO
7893195	4976.24	0.0914	ZU
7893195	16976.4	0.1571	AH
7893195	6642.81	0.2229	KO
7893195	5181.4	0.2886	AH
7893195	8186.85	0.3543	ZU
7893195	543.88	0.4200	KO
7893195	2767.34	0.4857	AH
7893195	3503.75	0.5514	AH
7893195	5646.97	0.6171	ZU
7985432	1625.6	0.3333	AH
7985432	6007.09	0.0171	AH
7985432	9309.43	0.0829	KO
7985432	2554.18	0.1486	AH
7985432	14081	0.2143	KO
7985432	1938.99	0.2800	KO
7985432	2127.68	0.3457	ZU
7985432	6490.22	0.4114	AH
7985432	13093.2	0.4771	ZU
7985432	2370	0.5429	AH
7985432	3034.79	0.6086	ZU
8609457	2629.98	0.8889	AH
8609457	3266.02	0.0057	KO
8609457	6573.52	0.0714	ZU
8609457	5726.06	0.1371	KO
8609457	13919.8	0.2029	AH
8609457	5175.39	0.2686	KO
8609457	1563.75	0.3343	ZU
8609457	3317.71	0.4000	KO
8609457	11042.2	0.4657	AH
8609457	6191.83	0.5314	KO
8609457	3106.54	0.5971	KO
8762365	10710.1	0.4444	ZU
8762365	10390.2	0.4444	AH
8762365	3769.6	0.0486	KO
8762365	3877.63	0.1143	ZU
8762365	12776.3	0.1800	KO
8762365	7184.28	0.2457	AH
8762365	570.75	0.3114	ZU
8762365	2574.02	0.3771	AH
8762365	3385.32	0.4429	KO
8762365	2440.11	0.5086	AH
8762365	2588.24	0.5743	ZU
8762365	1598.54	0.6400	AH
8764589	6235.53	0.8889	ZU
8764589	7318.75	0.2222	KO
8764589	3556.55	0.0571	ZU
8764589	8079.75	0.1229	AH
8764589	4582.75	0.1886	KO
8764589	2720.45	0.2543	AH
8764589	2249.32	0.3200	KO
8764589	4456.55	0.3857	KO
8764589	10328.1	0.4514	ZU
8764589	14134.1	0.5171	AH
8764589	9372.82	0.5829	ZU
8908712	2846.48	0.5556	KO
8908712	6355.41	0.7778	AH
8908712	8377.28	0.0514	AH
8908712	1724.81	0.1171	KO
8908712	3780.54	0.1829	AH
8908712	3842.75	0.2486	ZU
8908712	3941.74	0.3143	KO
8908712	10141.9	0.3800	AH
8908712	6379.06	0.4457	AH
8908712	9135.39	0.5114	ZU
8908712	10832	0.5771	KO
8908712	4100.04	0.6429	ZU
8923394	16299.4	0.1111	KO
8923394	1665.49	0.0371	ZU
8923394	11165.4	0.1029	AH
8923394	6492.05	0.1686	AH
8923394	1379.13	0.2343	AH
8923394	8787.07	0.3000	ZU
8923394	11570.8	0.3657	KO
8923394	2021.29	0.4314	ZU
8923394	5785.57	0.4971	KO
8923394	3702.73	0.5629	AH
8923394	1947.92	0.6286	ZU
9129938	2434.89	0.5556	KO
9129938	929.46	0.0029	ZU
9129938	11754.9	0.0686	AH
9129938	3575.18	0.1343	ZU
9129938	6434.74	0.2000	AH
9129938	7206.08	0.2657	ZU
9129938	4188.49	0.3314	AH
9129938	5805.92	0.3971	ZU
9129938	796.84	0.4629	KO
9129938	10338.7	0.5286	ZU
9129938	4977.87	0.5943	ZU
9871098	5769.37	0.9822	AH
9871098	1052.96	0.0400	KO
9871098	5279.43	0.1057	ZU
9871098	2122.15	0.1714	ZU
9871098	3295.09	0.2371	AH
9871098	8973.24	0.3029	KO
9871098	5101.27	0.3686	AH
9871098	7284.64	0.4343	KO
9871098	1726.34	0.5000	AH
9871098	8804.15	0.5657	ZU
9871098	10659.6	0.6314	KO

[RobPratt · Accepted solution]

Does the following PROC OPTMODEL code do what you want?

proc optmodel;
   /* declare parameters and read data */
   set GROUPS = /Test Control/;
   num groupSize {GROUPS} = [18, 5];

   set OBS;
   str PatientGroup {OBS};
   num Balance {OBS};
   num Risk {OBS};
   str prod {OBS};
   read data indata into OBS=[_N_] PatientGroup Balance Risk prod;

   set PATIENT_GROUPS = setof {i in OBS} PatientGroup[i];
   set PRODS = setof {i in OBS} prod[i];
   set OBS_p {PATIENT_GROUPS} init {};
   num patientCount {PATIENT_GROUPS} init 0;
   num sumBalance {PATIENT_GROUPS} init 0;
   num sumRisk {PATIENT_GROUPS} init 0;
   num sumProd {PATIENT_GROUPS, PRODS} init 0;
   str p_this, pr_this;
   for {i in OBS} do;
      p_this = PatientGroup[i];
      OBS_p[p_this] = OBS_p[p_this] union {i};
      patientCount[p_this] = patientCount[p_this] + 1;
      sumBalance[p_this] = sumBalance[p_this] + Balance[i];
      sumRisk[p_this] = sumRisk[p_this] + Risk[i];
      pr_this = prod[i];
      sumProd[p_this, pr_this] = sumProd[p_this, pr_this] + 1;
   end;
   num totalBalance = sum {p in PATIENT_GROUPS} sumBalance[p];
   num totalRisk = sum {p in PATIENT_GROUPS} sumRisk[p];
   num totalProd {pr in PRODS} = sum {p in PATIENT_GROUPS} sumProd[p,pr];

   /* Assign[p,g] = 1 if patient group p is assigned to group g, 0 otherwise */
   var Assign {PATIENT_GROUPS, GROUPS} binary;

   /* assign each patient group to exactly one group */
   con AssignOnce {p in PATIENT_GROUPS}:
      sum {g in GROUPS} Assign[p,g] = 1;

   /* assign groupSize[g] patient groups to group g */
   con GroupSizeCon {g in GROUPS}:
      sum {p in PATIENT_GROUPS} Assign[p,g] = groupSize[g];

   /* penalize absolute difference from fair share for Balance */
   var BalanceErrorPlus {GROUPS} >= 0;
   var BalanceErrorMinus {GROUPS} >= 0;
   con BalanceErrorCon {g in GROUPS}:
      (sum {p in PATIENT_GROUPS} sumBalance[p] * Assign[p,g]) / totalBalance - groupSize[g] / card(PATIENT_GROUPS) 
    = BalanceErrorPlus[g] - BalanceErrorMinus[g];
   impvar BalanceError = sum {g in GROUPS} (BalanceErrorPlus[g] + BalanceErrorMinus[g]);

   /* penalize absolute difference from fair share for Risk */
   var RiskErrorPlus {GROUPS} >= 0;
   var RiskErrorMinus {GROUPS} >= 0;
   con RiskErrorCon {g in GROUPS}:
      (sum {p in PATIENT_GROUPS} sumRisk[p] * Assign[p,g]) / totalRisk - groupSize[g] / card(PATIENT_GROUPS) 
    = RiskErrorPlus[g] - RiskErrorMinus[g];
   impvar RiskError = sum {g in GROUPS} (RiskErrorPlus[g] + RiskErrorMinus[g]);

   /* penalize absolute difference from fair share for Prod */
   var ProdErrorPlus {GROUPS, PRODS} >= 0;
   var ProdErrorMinus {GROUPS, PRODS} >= 0;
   con ProdErrorCon {g in GROUPS, pr in PRODS}:
      (sum {p in PATIENT_GROUPS} sumProd[p,pr] * Assign[p,g]) / totalProd[pr] - groupSize[g] / card(PATIENT_GROUPS) 
    = ProdErrorPlus[g,pr] - ProdErrorMinus[g,pr];
   impvar ProdError = sum {g in GROUPS, pr in PRODS} (ProdErrorPlus[g,pr] + ProdErrorMinus[g,pr]);

   /* minimize total error */
   min TotalError = BalanceError + RiskError + ProdError;

   /* call MILP solver */
   solve;

   /* create output data */
   impvar MeanBalance {g in GROUPS} = 
      (sum {p in PATIENT_GROUPS} sumBalance[p] * Assign[p,g]) / (sum {p in PATIENT_GROUPS} patientCount[p] * Assign[p,g]);
   impvar MeanRisk {g in GROUPS} = 
      (sum {p in PATIENT_GROUPS} sumRisk[p] * Assign[p,g]) / (sum {p in PATIENT_GROUPS} patientCount[p] * Assign[p,g]);
   impvar MeanProd {g in GROUPS, pr in PRODS} = 
      (sum {p in PATIENT_GROUPS} sumProd[p,pr] * Assign[p,g]) / (sum {p in PATIENT_GROUPS} patientCount[p] * Assign[p,g]);
   create data GroupOutdata from [Group] groupSize MeanBalance MeanRisk
      {pr in PRODS} <col(pr)=MeanProd[Group,pr]>;

   str assignedGroup {PATIENT_GROUPS};
   for {p in PATIENT_GROUPS} do;
      for {g in GROUPS: Assign[p,g].sol > 0.5} do;
         assignedGroup[p] = g;
         leave;
      end;
   end;
   create data PatientGroupOutdata from [PatientGroup=p] assignedGroup patientCount sumBalance sumRisk
      {pr in PRODS} <col(pr)=sumProd[p,pr]>;
quit;

[vicky07]

Hi @RobPratt,

First off thank you for your resoonse. Your solution worked perfectly when I ran with the sample data but when i expanded my input data to 1,251 unique patientgroups I am getting the below following error. What is that I am doing wrong? Please advice. 

Thanks.

NOTE: Problem generation will use 2 threads.
ERROR: A linear coefficient for constraint 'RiskErrorCon[Test]' is missing or invalid.
ERROR: A linear coefficient for constraint 'RiskErrorCon[Control]' is missing or invalid.
NOTE: The problem has 2526 variables (0 free, 0 fixed).
NOTE: The problem uses 3 implicit variables.
NOTE: The problem has 2502 binary and 0 integer variables.
NOTE: The problem has 1265 linear constraints (0 LE, 1265 EQ, 0 GE, 0 range).
NOTE: The problem has 18332 linear constraint coefficients.
NOTE: The problem has 0 nonlinear constraints (0 LE, 0 EQ, 0 GE, 0 range).
NOTE: The OPTMODEL presolver is disabled for linear problems.
NOTE: Unable to create problem instance due to previous errors.

[RobPratt]

Do you have any missing values in your data?

[vicky07]

Yes. Fixed the missing values and was able to run the code however I am getting "Out of memory" erorr this time. 

I've attached my full log, not sure if that helps.

Thanks. 

[RobPratt]

Are you able to share the data?

[vicky07]

Apologize for the delay.  Attached my data. Thanks!!

[RobPratt]

Thanks for sending the additional data.  I was able to replicate the issue and have several suggestions.

1. When the solver terminates (due to running out of memory or otherwise), it is supposed to return the best solution found so far.  There was a bug in SAS/OR 14.1 (SAS 9.4M3) where this didn't happen in the out-of-memory case.  That's why you see the "Division by zero" errors in the log.  A workaround is to use a SUBMIT block to call PROC OPTMILP, as described in this post.
2. One approach to avoid running out of memory is to use the MEMSIZE option to increase the amount of memory available to SAS, as described in the documentation here.
3. Another approach to avoid running out of memory is to use the MILP solver option NODESEL=DEPTH to do a depth-first search of the branch-and-cut tree.  This choice avoids creating a large tree, but the overall solve time might increase.
4. In SAS/OR 14.1, the default log frequency sometimes writes many lines per second.  For your problem, LOGFREQ=10000 yields a more reasonable frequency.
5. The version you are using is almost five years old, and we have had several releases since then, each with bug fixes, new features, and significant performance improvements.  I recommend upgrading to SAS/OR 15.1 (SAS 9.4M6) or SAS Optimization 8.5 (SAS Viya 3.5).
6. No matter which version you use, you might consider setting some termination criteria like RELOBJGAP=0.05 or MAXTIME=600 to stop early with a good solution.  Often with MILP problems, the solver spends a long time trying to proving optimality after it has already found an optimal solution.  In that case, letting the solver run longer would not yield any better solution.

[vicky07]

Thanks a lot!!