no way to compare when less than two revisions

Differences

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

@@ 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
+.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
+.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.