# 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] (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) + ​ + + [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.