calculating paule-mandel estimator

[rmo]

Hallo to all,

i am new to SAS, i am trying to program paule mandel estimator to calculate the pooled effect size of random-effects model for meta-analysis and the between study variance (tau^2). I have done it in Excel, but i am new to SAS, i am tryring to implement it through proc optmodel. In the attachments SAS code, excel file contains the calculations. the paper of paule mandel( to optimise are equations (1, 3, 6 and 7).

Thank you very much for help 

[RobPratt · Accepted solution]

Here are three approaches.  In all three, note that w is an implicit variable (rather than a number) because it depends on the variable t.

The first approach declares y as an implicit variable and minimizes dv:

/* impvar y, min dv */
proc optmodel;
   set OBS;
   number ET{OBS};
   number EST{OBS};
   read data mydata into OBS=[_n_] ET EST;

   var t >= 0;
   impvar w{i in OBS} = 1/(EST[i]+t);
   impvar y = (sum{i in OBS} w[i]*ET[i]) / (sum{i in OBS} w[i]);
   impvar t0 = sum{i in OBS} w[i]*(ET[i]-y)**2 - (card(OBS)-1);

   con t0Con: t0 >= 0;

   min dv = t0 / (sum{i in OBS}((w[i]**2)*((ET[i]-y)**2)));
   solve;
   t = t0 + dv;
   put ET[*];
   put EST[*];
   put w[*];
   print y ET w t0 t;
quit;

The solver yields an optimal objective value close to zero, but the resulting solution differs from your Excel spreadsheet:

Solution Summary
Solver	NLP
Algorithm	Interior Point Direct
Objective Function	dv
Solution Status	Optimal
Objective Value	9.0027692E-8
Optimality Error	9.0880058E-8
Infeasibility	0
Iterations	10
Presolve Time	0.00
Solution Time	0.00

y
-0.17114

[1]	ET	w
1	0.208718	10.3286
2	-0.598677	12.0514
3	-0.020166	8.1409

t0	t
1.8787	0.0000011299

The second approach declares y as a variable, imposes a constraint, and minimizes dv:

/* var y, min dv */
proc optmodel;
   set OBS;
   number ET{OBS};
   number EST{OBS};
   read data mydata into OBS=[_n_] ET EST;

   var y init 0;
   var t >= 0;
   impvar w{i in OBS} = 1/(EST[i]+t);
   impvar t0 = sum{i in OBS} w[i]*(ET[i]-y)**2 - (card(OBS)-1);
   t = t0;

   con t0Con: t0 >= 0;

   con yCon: y * sum{i in OBS} w[i] = sum{i in OBS} w[i] * ET[i];

   min dv = t0 / (sum{i in OBS}((w[i]**2)*((ET[i]-y)**2)));
   solve;
   t = t0 + dv;
   put ET[*];
   put EST[*];
   put w[*];
   print y ET w t0 t;
quit;

The solver again yields an optimal objective value close to zero, but the resulting solution differs from your Excel spreadsheet, except that the y value now matches:

Solution Summary
Solver	NLP
Algorithm	Interior Point Direct
Objective Function	dv
Solution Status	Optimal
Objective Value	9.002769E-8
Optimality Error	9.0880059E-8
Infeasibility	2.751088E-11
Iterations	14
Presolve Time	0.00
Solution Time	0.01

y
-0.15463

[1]	ET	w
1	0.208718	10.3286
2	-0.598677	12.0514
3	-0.020166	8.1409

t0	t
1.887	0.0000011299

The third approach declares y as an implicit variable, enforces equation (5) as a constraint, and uses no objective:

/* impvar y, no objective */
proc optmodel;
   set OBS;
   number ET{OBS};
   number EST{OBS};
   read data mydata into OBS=[_n_] ET EST;

   var t >= 0;
   impvar w{i in OBS} = 1/(EST[i]+t);
   impvar y = (sum{i in OBS} w[i]*ET[i]) / (sum{i in OBS} w[i]);
   con Equation5: sum{i in OBS} w[i]*(ET[i]-y)**2 / (card(OBS)-1) = 1;

   solve noobj;
   put ET[*];
   put EST[*];
   put w[*];
   print y ET w t;
quit;

The resulting solution matches your Excel spreadsheet:

Solution Summary
Solver	NLP
Algorithm	Interior Point Direct
Objective Function	(0)
Solution Status	Optimal
Objective Value	0
Optimality Error	9.0909091E-8
Infeasibility	4.440892E-16
Iterations	6
Presolve Time	0.00
Solution Time	0.00

y
-0.15463

[1]	ET	w
1	0.208718	5.5346
2	-0.598677	5.9938
3	-0.020166	4.8380

t
0.083863

[rmo]

Thank very much for your help.
The third solution is the best.
It gives the same result in R(with the package meta), if the fixed effects model used to estimate the pooled effect size and between-study variance