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] (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)