# The metafor Package

A Meta-Analysis Package for R

### Site Tools

tips:clogit_paired_binary_data

# Differences

This shows you the differences between two versions of the page.

 tips:clogit_paired_binary_data [2019/07/03 08:12] tips:clogit_paired_binary_data [2019/07/03 08:12] (current) Line 1: Line 1: + ===== 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) + + +  0.8109 + + And the corresponding SE: + + round(sqrt(1/tab[2,1] + 1/tab[1,2]),4) + + +  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. 