The metafor Package

A Meta-Analysis Package for R

User Tools

Site Tools


tips:clogit_paired_binary_data

Differences

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

Link to this comparison view

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:
 +<code rsplus>
 +### 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)
 +</code>
 +
 +Now we can fit the conditional logistic regression model with ''clogit()'' from the ''survival'' package:
 +<code rsplus>
 +library(survival)
 +summary(clogit(event ~ group + strata(id)))
 +</code>
 +Results:
 +<code output>
 +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
 +</code>
 +
 +For paired binary data, the same results can be obtained by fitting a mixed-effects logistic regression model:
 +<code rsplus>
 +library(lme4)
 +summary(glmer(event ~ group + (1 | id), family=binomial, nAGQ=17))
 +</code>
 +Note that it is necessary to increase the number of quadrature points quite a bit to get sufficient accuracy here. Results:
 +<code output>
 +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
 +</code>
 +
 +We can also use ''glmmML()'' for this:
 +<code rsplus>
 +library(glmmML)
 +glmmML(event ~ group, cluster=id, family=binomial, method="ghq", n.points=17)
 +</code>
 +Results:
 +<code output>
 +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
 +</code>
 +
 +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:
 +<code rsplus>
 +tab <- table(ai,ci)
 +</code>
 +Then we can obtain the MLE of the log odds ratio with:
 +<code rsplus>
 +round(log(tab[2,1]/tab[1,2]),4)
 +</code>
 +<code output>
 +[1] 0.8109
 +</code>
 +And the corresponding SE:
 +<code rsplus>
 +round(sqrt(1/tab[2,1] + 1/tab[1,2]),4)
 +</code>
 +<code output>
 +[1] 0.4249
 +</code>
 +
 +Finally, let's use ''rma.glmm()'' for this:
 +<code rsplus>
 +library(metafor)
 +rma.glmm(measure="OR", ai=ai, bi=bi, ci=ci, di=di, model="CM.EL", method="FE")
 +</code>
 +<code output>
 +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
 +</code>
 +Results again match what we obtained above.
tips/clogit_paired_binary_data.txt · Last modified: 2019/07/03 08:12 (external edit)