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 probability 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
the simulation of the data,
the post-hoc chi-square test,
fitting of a logistic mixed model,
fitting logistic regression
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!
set.seed(12345), so that everyone gets the same results? $\endgroup$seedand adapted the ouput. Any idea why there is a difference between output of post hoc chi square and logistic regression? $\endgroup$glmnet? $\endgroup$