I noticed the other day that I can mimic a stratified logistic regression using a meta-analytic approach at least when there is only one stratification variable with two different values. This is the aggregated data I use with metabin from the package meta.

> dmeta
  verum_event verum_N plac_event plac_N sex
1          25      50         30    100   0
2          30      80         25     70   1

I get the following (leaving out the heterogeneity stuff)

> meta::metabin(verum_event,verum_N,plac_event,plac_N,data=dmeta,sm="OR",method="Inverse",comb.random=F)
      OR           95%-CI %W(fixed)
1 2.3333 [1.1585; 4.6996]      47.5
2 1.0800 [0.5547; 2.1027]      52.5

Number of studies combined: k = 2

                       OR           95%-CI    z p-value
Fixed effect model 1.5574 [0.9611; 2.5236] 1.80  0.0720

Using an expanded data set d

> car::some(d)
    y trt sex
55  0   0   0
62  0   0   0
65  0   0   0
83  1   0   0
131 1   1   0
137 1   1   0
168 0   0   1
178 0   0   1
210 1   0   1
218 1   0   1

I get with clogit:

> summary(clogit(y~trt+strata(sex),data=d,method="exact"));
    exp(coef) exp(-coef) lower .95 upper .95
trt     1.557     0.6421    0.9612     2.524

exactly the same. Well, almost. But the SAS proc freq results are exactly the same as the meta-analysis. I expected as much since stratified analysis is computing the effect in each stratum separately and combine them afterwards. Which is also what a meta-analysis does. However the same approach does not work with relative risk (for risk differences the jury is still out). As I don't know of any way to compute stratified relative risks with R (is there one?) I rely on the SAS output which says

RR     lower  upper
1.1595 0.9585 1.4027

But the meta analysis looks different:

> meta::metabin(verum_event,verum_N,plac_event,plac_N,data=dmeta,sm="RR",method="Inverse",comb.random=F)
      RR           95%-CI %W(fixed)
1 1.6667 [1.1083; 2.5063]      51.8
2 1.0500 [0.6879; 1.6026]      48.2

Number of studies combined: k = 2

                       RR           95%-CI    z p-value
Fixed effect model 1.3338 [0.9945; 1.7890] 1.92  0.0545

Why do these two methods agree when using the logit link but don't when using log link?

  • $\begingroup$ Can you give more detail about the SAS approach? I would expect the approaches to agree for both links. $\endgroup$ – Thomas Lumley Jul 15 '20 at 21:26
  • $\begingroup$ @ThomasLumley As I don't have much experience with SAS I asked a colleague for the syntax. He gave me this: proc freq data=d; tables sex*trt*y / CMH cl riskdiff(common); I thought that one would need proc glm for this but as the odds ratio results where the same as those with clogit I think it is correct. SAS actually gives two results for RR but I guess that's due to the fact that contrary to OR the estimate is not symmetric when switching event and non-event. But now that I look at the separate result I'm confused, they are different from what I get. sex=0 sas RR=1.4 (or 0.6) R RR=1.6 $\endgroup$ – diffset Jul 16 '20 at 7:10

It turns out the problem was in my inexperience in coding SAS and understanding its output. After changing the input of the data such that one of the lines 1 1 1 comes first and adding order=data in the proc freq call I indeed also get the same result for relative risks (column 1) as in the meta-analysis. As for risk difference I get the same estimate (when using Mantel-Haenszel) but a somewhat larger standard error in SAS which might be due to an appropriate option I am still missing regarding SAS.

In summary, I am very satisfied to see that these two methods are indeed equivalent. This allows one to calculate stratified risk ratios in R quite easily, at least in such simple cases.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.