# The metafor Package

A Meta-Analysis Package for R

### Sidebar

tips:clogit_paired_binary_data

## Conditional Logistic Regression for Paired Binary Data

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="FE")
Fixed-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.