Stratified logistic regression with log link

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?

• Can you give more detail about the SAS approach? I would expect the approaches to agree for both links. – Thomas Lumley Jul 15 '20 at 21:26
• @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 – 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.