The metafor Package

This is just a short illustration of how to fit the conditional logistic regression model for paired binary data using various functions, including the rma.glmm() function from the metafor package.

Fist, let's create an illustrative dataset:

### number of pairs
n <- 100
### create data
ai <- c(rep(0,n/2), rep(1,n/2))
bi <- 1-ai
ci <- c(rep(0,42), rep(1,8), rep(0,18), rep(1,32))
di <- 1-ci
### change data to long format
event <- c(rbind(ai,ci))
group <- rep(c(1,0), times=n)
id    <- rep(1:n, each=2)

Now we can fit the conditional logistic regression model with clogit() from the survival package:

library(survival)
summary(clogit(event ~ group + strata(id)))

Results:

Call:
coxph(formula = Surv(rep(1, 200L), event) ~ group + strata(id),
    method = "exact")
 
  n= 200, number of events= 90
 
        coef exp(coef) se(coef)     z Pr(>|z|)
group 0.8109    2.2500   0.4249 1.908   0.0563 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
      exp(coef) exp(-coef) lower .95 upper .95
group      2.25     0.4444    0.9783     5.175
 
Rsquare= 0.02   (max possible= 0.165 )
Likelihood ratio test= 3.95  on 1 df,   p=0.04695
Wald test            = 3.64  on 1 df,   p=0.05633
Score (logrank) test = 3.85  on 1 df,   p=0.04986

For paired binary data, the same results can be obtained by fitting a mixed-effects logistic regression model:

library(lme4)
summary(glmer(event ~ group + (1 | id), family=binomial, nAGQ=17))

Note that it is necessary to increase the number of quadrature points quite a bit to get sufficient accuracy here. Results:

Generalized linear mixed model fit by maximum likelihood (Adaptive Gauss-Hermite Quadrature, nAGQ = 17) ['glmerMod']
 Family: binomial  ( logit )
Formula: event ~ group + (1 | id)
 
     AIC      BIC   logLik deviance df.resid
   253.9    263.8   -124.0    247.9      197
 
Scaled residuals:
    Min      1Q  Median      3Q     Max
-1.1703 -0.4314 -0.2876  0.5031  1.2817
 
Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 7.207    2.685
Number of obs: 200, groups:  id, 100
 
Fixed effects:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.8110     0.4255  -1.906   0.0566 .
group         0.8109     0.4249   1.908   0.0563 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Correlation of Fixed Effects:
      (Intr)
group -0.555

We can also use glmmML() for this:

library(glmmML)
glmmML(event ~ group, cluster=id, family=binomial, method="ghq", n.points=17)

Results:

Call:  glmmML(formula = event ~ group, family = binomial, cluster = id,      method = "ghq", n.points = 17)
 
 
               coef se(coef)      z Pr(>|z|)
(Intercept) -0.8110   0.4255 -1.906   0.0566
group        0.8109   0.4249  1.908   0.0563
 
Scale parameter in mixing distribution:  2.685 gaussian
Std. Error:                              0.661
 
        LR p-value for H_0: sigma = 0:  2.461e-07
 
Residual deviance: 247.9 on 197 degrees of freedom      AIC: 253.9

For binary matched pair data, the MLE and corresponding SE can actually be given in closed form (see, for example, Agresti, 2002, Categorical data analysis, chapter 10). First, we create the 2x2 table with:

tab <- table(ai,ci)

Then we can obtain the MLE of the log odds ratio with:

round(log(tab[2,1]/tab[1,2]),4)

[1] 0.8109

And the corresponding SE:

round(sqrt(1/tab[2,1] + 1/tab[1,2]),4)

[1] 0.4249

Finally, let's use rma.glmm() for this:

library(metafor)
rma.glmm(measure="OR", ai=ai, bi=bi, ci=ci, di=di, model="CM.EL", method="EE")

Equal-Effects Model (k = 26)
Model Type: Conditional Model with Exact Likelihood
 
Tests for Heterogeneity:
Wld(df = 25) =      NA, p-val = NA
LRT(df = 25) = 32.0966, p-val = 0.1552
 
Model Results:
 
estimate      se    zval    pval    ci.lb   ci.ub
  0.8109  0.4249  1.9085  0.0563  -0.0219  1.6437  .
 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Results again match what we obtained above.