On va utiliser la librairie survival qui normalement sert à faire les modèles de survie, mais qui a aussi la fonction clogit (conditional logit)1.
library(package = "survival")## Loading required package: splinesOn importe les données et on garde que les paires où les deux individus ont eu un transfert.
tab <- read.csv("~/Documents/bdr_tim/CJPOK.csv")
# on enlève la première colonne
tab <- tab[-1]
names(tab) <- c("pair",
"nTrans_0", "nSacs_0", "trans_0", "pImpl_0",
"nTrans_1", "nSacs_1", "trans_1", "pImpl_1")
tab <- subset(tab, trans_0 & trans_1)On remet les données en forme pour avoir une ligne par transfert.
pair_0 <- with(tab, rep(pair, nTrans_0))
nTrans_0 <- with(tab, rep(nTrans_0, nTrans_0))
impl_0 <- unlist(with(tab,
mapply(function(x, y) c(rep(0, x), rep(1, y)),
x = nTrans_0 - nSacs_0,
y = nSacs_0)))
ttt_0 <- with(tab, rep(rep("frais", nrow(tab)), nTrans_0))
pair_1 <- with(tab, rep(pair, nTrans_1))
nTrans_1 <- with(tab, rep(nTrans_1, nTrans_1))
impl_1 <- unlist(with(tab,
mapply(function(x, y) c(rep(0, x), rep(1, y)),
x = nTrans_1 - nSacs_1,
y = nSacs_1)))
ttt_1 <- with(tab, rep(rep("tdiff", nrow(tab)), nTrans_1))
tabReg <- data.frame(pair = c(pair_0, pair_1),
nTrans = as.factor(c(nTrans_0, nTrans_1)),
impl = c(impl_0, impl_1),
ttt = c(ttt_0, ttt_1))
head(tabReg, 10)## pair nTrans impl ttt
## 1 1 2 0 frais
## 2 1 2 0 frais
## 3 2 2 0 frais
## 4 2 2 0 frais
## 5 4 2 0 frais
## 6 4 2 0 frais
## 7 6 2 0 frais
## 8 6 2 0 frais
## 9 7 2 0 frais
## 10 7 2 0 fraisProportion d’échec des transferts selon le traitement et le nombre de transferts.
round(with(tabReg, tapply(impl, list(ttt, nTrans),
function(x) prop.table(table(x))[1])), 2)## 1 2 3
## frais 0.56 0.76 1
## tdiff 0.59 0.84 1Pas vraiment de différences entre les groupes de traitement, par contre plus t’en transfères plus ça rate. Enfin je suppose que moins t’as de chances que ça s’implante plus t’en transfères.
Comme la probabilité d’échec est de 100% pour les individus avec 3 transferts, on le groupe avec le groupe à 2 transferts, sinon le modèle va pas converger.
tabReg <- within(tabReg, {
nTrans_2 <- as.factor(ifelse(nTrans == 3, 2, nTrans))
})On fait la regression logistique conditionnelle. En gros avec cette regression tu as un risque de base par strate (ici les paires). Mais ce risque de base n’est pas calculé, il est juste pris en compte par le modèle. C’est pour ça que dans les résultats tu n’as pas d’intercept, et pas de test de l’effet paire.
mod1 <- clogit(impl ~ ttt * nTrans_2 + strata(pair), data = tabReg)
summary(mod1)## Call:
## coxph(formula = Surv(rep(1, 622L), impl) ~ ttt * nTrans_2 + strata(pair),
## data = tabReg, method = "exact")
##
## n= 622, number of events= 160
##
## coef exp(coef) se(coef) z Pr(>|z|)
## ttttdiff -0.375 0.687 0.479 -0.78 0.433
## nTrans_22 -1.069 0.343 0.512 -2.09 0.037 *
## ttttdiff:nTrans_22 0.298 1.347 0.644 0.46 0.644
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## ttttdiff 0.687 1.455 0.269 1.757
## nTrans_22 0.343 2.912 0.126 0.937
## ttttdiff:nTrans_22 1.347 0.742 0.381 4.759
##
## Rsquare= 0.013 (max possible= 0.312 )
## Likelihood ratio test= 8.14 on 3 df, p=0.0431
## Wald test = 7.88 on 3 df, p=0.0486
## Score (logrank) test = 8.23 on 3 df, p=0.0414Pas d’interaction entre le traitement et le nombre de transferts.
mod2 <- clogit(impl ~ ttt + nTrans_2 + strata(pair), data = tabReg)
summary(mod2)## Call:
## coxph(formula = Surv(rep(1, 622L), impl) ~ ttt + nTrans_2 + strata(pair),
## data = tabReg, method = "exact")
##
## n= 622, number of events= 160
##
## coef exp(coef) se(coef) z Pr(>|z|)
## ttttdiff -0.186 0.831 0.243 -0.76 0.4460
## nTrans_22 -0.890 0.411 0.334 -2.67 0.0077 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## ttttdiff 0.831 1.20 0.515 1.34
## nTrans_22 0.411 2.44 0.213 0.79
##
## Rsquare= 0.013 (max possible= 0.312 )
## Likelihood ratio test= 7.93 on 2 df, p=0.019
## Wald test = 7.7 on 2 df, p=0.0213
## Score (logrank) test = 8.02 on 2 df, p=0.0182Pas d’effet traitement.
On peut aussi tenter avec un modèle linéraire généralisé avec une intercept aléatoire par paire.
library(package = "lme4")## Loading required package: Matrix
## Loading required package: RcppmodMix_1 <- glmer(impl ~ ttt * nTrans_2 + (1 | pair),
data = tabReg,
family = binomial(link = logit))
summary(modMix_1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: impl ~ ttt * nTrans_2 + (1 | pair)
## Data: tabReg
##
## AIC BIC logLik deviance df.resid
## 683.5 705.7 -336.7 673.5 617
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.022 -0.530 -0.422 0.793 2.455
##
## Random effects:
## Groups Name Variance Std.Dev.
## pair (Intercept) 0.595 0.771
## Number of obs: 622, groups: pair, 194
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.276 0.356 -0.78 0.4386
## ttttdiff -0.207 0.409 -0.51 0.6129
## nTrans_22 -1.088 0.392 -2.77 0.0056 **
## ttttdiff:nTrans_22 -0.274 0.506 -0.54 0.5885
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ttttdf nTr_22
## ttttdiff -0.845
## nTrans_22 -0.897 0.798
## ttttdf:T_22 0.695 -0.834 -0.766Pas d’interaction non plus.
modMix_2 <- glmer(impl ~ ttt + nTrans_2 + (1 | pair),
data = tabReg,
family = binomial(link = logit))
summary(modMix_2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: impl ~ ttt + nTrans_2 + (1 | pair)
## Data: tabReg
##
## AIC BIC logLik deviance df.resid
## 681.8 699.5 -336.9 673.8 618
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.008 -0.525 -0.419 0.801 2.383
##
## Random effects:
## Groups Name Variance Std.Dev.
## pair (Intercept) 0.613 0.783
## Number of obs: 622, groups: pair, 194
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.143 0.256 -0.56 0.578
## ttttdiff -0.392 0.224 -1.75 0.081 .
## nTrans_22 -1.252 0.252 -4.97 6.6e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ttttdf
## ttttdiff -0.663
## nTrans_22 -0.786 0.439Pas d’effet traitement non plus. Attention en comprarant les coefficients, clogit utilise une fonction de lien exponentielle comme le modèle de Cox, alors que glmer utilise une fonction de lien logit.
Dans les exemples la regression logistique conditionelle estime donc des risques relatifs, alors que le modèle linéaire généralisé estime des odds-ratios.
En fait c’est parce que une régression logistique conditionnelle est équivalente à un modèle de survie stratifié sans censure à temps unique… Peu importe.↩