1
$\begingroup$

I have data that looks like this: there is a group of 27 subjects with one dichotomous variable $y_1$ at 3 times points. The propbability of $y_1$ is different between the 3 time points (100%, 85%, and 40% respectively). I want to prove that their is difference a significant difference between each of the time point (between $t_1-t_2$, $t_1-t_3$ and $t_2-t_3$). For this I would like to use a logistic regression (ideally a mixed model). I know there is a statistical difference between the following pairs $t_1-t_3$ and $t_2-t_3$ as showed by the result of the post-hoc chi-square test. However I can't find a logistic model that proves this. Instead, the logistic models seem to tell that there are no difference in the pair $t_1-t_3$

Here is first 1. the simulation of the data, 2. the post-hoc chi-square test, 3. fitting of a logistic mixed model, 4. fitting logistic regression 5. fitting logistic regression with a score test as pointed in this answer

1. Simulation of the data

library(dplyr)
library(multcomp)
library(lme4)
library(fifer)

# N.B.: if you don't have the package "fifer", it 
# has to be downloaded there: 
# "https://cran.r-project.org/src/contrib/Archive/fifer/"
# and install following the instructions there:
# "http://www.ryantmoore.org/files/ht/htrtargz.pdf"
# the whole operation takes 1 minute (really)

set.seed(12345)

funda = 25
y_sub1 = rbinom(funda, 1, 0.85)
y_sub2 = rbinom(funda, 1, 0.4)
y1 = c(rep(1,funda), y_sub1, y_sub2)
y1[2*funda + 3] = NA
time = c(rep("t1", funda), rep("t2", funda), rep("t3", funda))
ID_init = c()
for (i in 1:funda){
  ID_init = c(ID_init, paste0("ID", i))
}
ID = rep(ID_init, 3)
df = data.frame(y1, time, ID)

y_sub1 = rbinom(funda, 1, 0.87)
y_s1 = as.data.frame(table(y_sub1))

y_sub2 = rbinom(funda, 1, 0.37)
y_s2 = as.data.frame(table(y_sub2))

M <- as.table(rbind(c(0, funda), c(y_s1[1,2], y_s1[2,2]), c(y_s2[1,2], y_s2[2,2])))
dimnames(M) <- list(time=c("t1","t2","t3"),proc=c("No PV isol", "PV isol"))
M

Output:

proc
time No PV isol PV isol
  t1          0      25
  t2          5      20
  t3         15      10

2. Post hoc chi square test

chisq.post.hoc(M)


Adjusted p-values used the fdr method.

  comparison  raw.p  adj.p
1  t1 vs. t2 0.0502 0.0502
2  t1 vs. t3 0.0000 0.0000
3  t2 vs. t3 0.0086 0.0129

3. Logistic mixed model

mod_glmer = glmer(y1 ~ time + (1|ID), data = df, family = binomial)
mod_glmer = glht(mod_glmer, linfct = mcp(time = "Tukey"))
mod_glmer = summary(mod_glmer)
mod_glmer

Output:

unable to evaluate scaled gradientModel failed to converge: degenerate  Hessian with 1 negative eigenvaluesvariance-covariance matrix computed from finite-difference Hessian is
not positive definite or contains NA values: falling back to var-cov estimated from RX

 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: glmer(formula = y1 ~ time + (1 | ID), data = df, family = binomial)

Linear Hypotheses:
              Estimate Std. Error z value Pr(>|z|)   
t2 - t1 == 0  -16.3995  1573.5380  -0.010  0.99993   
t3 - t1 == 0  -18.5930  1573.5379  -0.012  0.99991   
t3 - t2 == 0   -2.1935     0.7119  -3.081  0.00412 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

4. logistic regression:

mod_glm = glm(y1 ~ time, data = df, family = "binomial")
mod_glm <- glht(mod_glm, linfct = mcp(time = "Tukey")) %>% summary()
mod_glm

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: glm(formula = y1 ~ time, family = "binomial", data = df)

Linear Hypotheses:
          Estimate Std. Error z value Pr(>|z|)    
t2 - t1 == 0  -18.1798  2150.8027  -0.008   1.0000  
t3 - t1 == 0  -19.9025  2150.8026  -0.009   0.9999  
t3 - t2 == 0   -1.7228     0.6492  -2.654   0.0159 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

5. logistic regression:

mod_glm = glm(y1 ~ time, data = df, family = "binomial")
mod_glm = anova(mod_glm, test = "Rao")
mod_glm
mod_glm = glht(mod_glm, linfct = mcp(time = "Tukey")) %>% summary()
mod_glm

However i get this error message:

Error in factor_contrasts(model) : no 'model.matrix' method for 'model' found!
$\endgroup$
  • $\begingroup$ The best model is the mixed one. $\endgroup$ – mdewey Sep 23 '18 at 13:28
  • 1
    $\begingroup$ As you say in your question, the mixed model is the preferable one. When I run your code, I get very different answers for each run. Maybe you want to start your sample code with e.g. set.seed(12345), so that everyone gets the same results? $\endgroup$ – Sal Mangiafico Sep 23 '18 at 13:39
  • $\begingroup$ Thank you for your comment @SalMangiafico. I added your seed and adapted the ouput. Any idea why there is a difference between output of post hoc chi square and logistic regression? $\endgroup$ – ecjb Sep 23 '18 at 14:11
  • $\begingroup$ If your outcome variable is genuinely always 1 as you suggest in your introduction then you have separation in your mixed model which probably accounts for (a) the huge estimates on the log scale, (b) the huge standard errors, (c) the failure of the model to converge. See stats.stackexchange.com/questions/11109/… for some more information. $\endgroup$ – mdewey Sep 23 '18 at 14:52
  • $\begingroup$ Dear @mdewey. Thank you very much for your answer. Provided the simulated data, do you have an idea of the code which should be run to perform a logistic regression using glmnet? $\endgroup$ – ecjb Sep 23 '18 at 15:51
2
$\begingroup$

In general, we would expect chi-square tests to have a different outcome than the mixed model, in just the same way that we would expect a t-test to have a different result than a paired t-test: 1) They are not the same test or model, and 2) Taking into account the repeated measures nature of the data would likely give different results.

However, in this case, there is a particular problem caused with there being no values of 0 in the dependent variable for Time 1. I'll tease this out in the code below.

The comment by @mdewey notes this also.


Install packages. I'm using the package rcompanion † rather than fifer, just because it can be downloaded from CRAN as a current package.

if(!require(dplyr)){install.packages("dplyr")}
if(!require(lme4)){install.packages("lme4")}
if(!require(emmeans)){install.packages("emmeans")}
if(!require(rcompanion)){install.packages("rcompanion")}

library(dplyr)
library(lme4)
library(emmeans)
library(rcompanion)

Create data.

set.seed(12345)

funda = 25
y_sub1 = rbinom(funda, 1, 0.85)
y_sub2 = rbinom(funda, 1, 0.4)
y1 = c(rep(1,funda), y_sub1, y_sub2)
y1[2*funda + 3] = NA
time = c(rep("t1", funda), rep("t2", funda), rep("t3", funda))
ID_init = c()
for (i in 1:funda){
  ID_init = c(ID_init, paste0("ID", i))
}
ID = rep(ID_init, 3)
df = data.frame(y1, time, ID)

Model the original data. Note, the pairwise Fisher exact tests are the same as in the question. ‡

M = xtabs(~ time + y1, data=df)
M

   ###     y1
   ### time  0  1
   ###   t1  0 25
   ###   t2  5 20
   ###   t3 14 10

pairwiseNominalIndependence(M, chisq=F, gtest=F)

   ###   Comparison p.Fisher p.adj.Fisher
   ### 1    t1 : t2 5.02e-02     5.02e-02
   ### 2    t1 : t3 2.90e-06     8.70e-06
   ### 3    t2 : t3 8.58e-03     1.29e-02

mod_glmer = glmer(y1 ~ time + (1|ID), data = df, family = binomial, nAGQ=1)

   ### Errors. These could be eliminated by using nAGQ=0,
   ###   but that's probably not a great solution.

Here, I'm going to use the emmeans package rather than multcomp. Note the huge SE on Time 1. Below I'll show that this is related to there being no 0 values for Time 1.

emmeans(mod_glmer, ~ time)


   ###  time    emmean           SE  df     asymp.LCL    asymp.UCL
   ###  t1   18.154318 1573.5378762 Inf -3065.9232479 3102.2318833
   ###  t2    1.754808    0.6023621 Inf     0.5741995    2.9354155
   ###  t3   -0.438702    0.4940877 Inf    -1.4070962    0.5296922
   ###
   ### Results are given on the logit (not the response) scale. 
   ### Confidence level used: 0.95 

The constrasts including Time 1 have a huge SE.

pairs(emmeans(mod_glmer, ~ time))

   ###  contrast estimate          SE  df z.ratio p.value
   ###  t1 - t2  16.39951 1573.537953 Inf   0.010  0.9999
   ###  t1 - t3  18.59302 1573.537923 Inf   0.012  0.9999
   ###  t2 - t3   2.19351    0.711923 Inf   3.081  0.0058

Now I'm going to change the data so that there is one value of 0 in Time 1. You'll see that the rest of the results are more satisfactory.

df2 = df

df2$y1[2]=0

M2 = xtabs(~ time + y1, data=df2)
M2

   ###     y1
   ### time  0  1
   ###   t1  1 24
   ###   t2  5 20
   ###   t3 14 10

The pairwise Fisher tests are little different than the original data, but the results are similar.

pairwiseNominalIndependence(M2, chisq=F, gtest=F)

   ###   Comparison p.Fisher p.adj.Fisher
   ### 1    t1 : t2 0.189000     0.189000
   ### 2    t1 : t3 0.000034     0.000102
   ### 3    t2 : t3 0.008580     0.012900

The GLMM returns no errors.

mod_glmer2 = glmer(y1 ~ time + (1|ID), data = df2, family = binomial, nAGQ=1)

   ### No errors

The estimated marginal means don't show the inflated SE.

emmeans(mod_glmer2, ~ time)

   ###  time     emmean        SE  df  asymp.LCL asymp.UCL
   ###  t1    3.6764923 1.2467438 Inf  1.2329194 6.1200652
   ###  t2    1.7274818 0.7280539 Inf  0.3005224 3.1544412
   ###  t3   -0.4284324 0.5319783 Inf -1.4710907 0.6142258
   ### 
   ### Results are given on the logit (not the response) scale. 
   ### Confidence level used: 0.95 

The pairwise contrasts have similar results to the Fisher tests.

pairs(emmeans(mod_glmer2, ~ time))

   ###  contrast estimate        SE  df z.ratio p.value
   ###  t1 - t2  1.949010 1.1906833 Inf   1.637  0.2301
   ###  t1 - t3  4.104925 1.3742629 Inf   2.987  0.0079
   ###  t2 - t3  2.155914 0.9025752 Inf   2.389  0.0446
   ###
   ### Results are given on the log odds ratio (not the response) scale. 
   ### P value adjustment: tukey method for comparing a family of 3 estimates 

ADDED:

One solution I've found to this problem is to use the blmer package. This actually uses a Bayesian approach with a weak prior. Note that the fixef.prior option establishes the prior. In this case the 3 refers to the three parameters for the fixed effects, and the 9 indicates a reasonably large variance.

I'll leave it as an exercise for the reader to read up on this solution, as it's not an approach I'm particularly familiar with.

But the results seem reasonable.

if(!require(blme)){install.packages("blme")}

library(blme)

glmb = bglmer(y1 ~ time + (1|ID), data=df, family=binomial,
              fixef.prior = normal(cov = diag(9,3)))

joint_tests(glmb)

   ### model term df1 df2 F.ratio p.value
   ### time         2 Inf    7.41  0.0006

emmeans(glmb, ~ time)

   ### time     emmean        SE  df  asymp.LCL asymp.UCL
   ### t1    3.7347886 1.0020341 Inf  1.7708379 5.6987393
   ### t2    1.8768716 0.7072334 Inf  0.4907196 3.2630236
   ### t3   -0.3798005 0.5557619 Inf -1.4690738 0.7094727

pairs(emmeans(glmb, ~ time))

   ### contrast estimate        SE  df z.ratio p.value
   ### t1 - t2  1.857917 1.0112966 Inf   1.837  0.1575
   ### t1 - t3  4.114589 1.1025738 Inf   3.732  0.0006
   ### t2 - t3  2.256672 0.8517203 Inf   2.650  0.0220

# summary(glmb)

† Caveat, I am the creator of this package.

‡ There is a slight discrepancy in the counts for Time 3, as in the original post an NA value is counted as a 0.

$\endgroup$
  • $\begingroup$ Dear @Salmangiafico. Thank you very much for your very detailed answer and the time you took. I must say that I already suspected that the problem came from the fact that there are no values of 0 in the dependent variable for Time 1, although I don't understand why it causes mathematically a problem. The simulated data that I provided reflect the real data for which i have want show the difference between t1-t2 and t1-t2. I would like to have a solution using a mixed model without changing the data (adding a 0 instead of a 1 as you demonstrated). Does that solution exist at all? $\endgroup$ – ecjb Sep 23 '18 at 15:47
  • $\begingroup$ @ecjb , I've added a potential solution to my answer. $\endgroup$ – Sal Mangiafico Sep 23 '18 at 18:49
  • $\begingroup$ Many Thanks @Salmangiafico!!!!!! That worked! May I ask you how did you come across this solution. Do you have any particular recommendation for readings/refrences? How did you pick up the value 9 for the large variance. How should I do for other tests? $\endgroup$ – ecjb Sep 23 '18 at 19:10
  • $\begingroup$ No, it was just googling things like "r glmm complete separation". I couldn't find anything non-bayesian that worked with glmm, and I wasn't happy with anything that worked on glm. $\endgroup$ – Sal Mangiafico Sep 23 '18 at 19:17
  • $\begingroup$ The 9 is just from the online example I linked to. You can change that value, and the estimates and p-values will change. But as long as it makes for a sufficiently bland prior, the results don't change much. $\endgroup$ – Sal Mangiafico Sep 23 '18 at 19:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.