intervalle de confiance

[Rollot Fabien]

Bonjour,
j'ai un jeu de données composé de quatre colonnes ( probabilités de décès,et les trois autres colonnes sont des concentrations de produits). Mon problème est le suivant:
j'ai trouvé le meilleur modèle de ce jeu de données et je voudrais calculer l'intervalle de confiance de niveau 95 % pour la valeur de la probabilité de décès au point moyen des concentrations.
J'aurais voulu savoir quelle fonction prendre pour faire cela, j'ai essayé avec predict et l'option interval="confidence" mais je ne comprend pas ce que cela m'affiche.
si quelqu'un pouvait m'aider??

[Renaud Lancelot]

Faire un data.frame avec la ou les valeurs pour lesquelles vous voulez faire la prédiction, puis l'utiliser dans predict:

Code : Tout sélectionner

> set.seed(121)
> dfr <- data.frame(
+     n = rep(10, 20),
+     y = rbinom(n = 20, size = 10, prob = seq(0.05, .95, length = 20)),
+     x = seq(0.05, .95, length = 20))
>
> fm <- glm(cbind(y, n - y) ~ x, family = binomial, data = dfr)
>
> New <- data.frame(x = mean(dfr$x))
> predict(fm, newdata = New, type = "response", se = TRUE)
$fit
        1
0.4729007

$se.fit
         1
0.04111898

$residual.scale
[1] 1


Renaud

[Rollot Fabien]

Oups, là je suis perdu!!!!!
Mon meilleur modèle est le suivant:

logit(p)=beta0 + beta1*log(a) +beta2*log(b) + beta3*log(c) + beta4*log(a)*log(b)
où p est la probabilité de décès et a, b, c des concentrations de produits

Mon souci est de déterminer un intervalle de confiance pour la valeur p au point (mean(a),mean(b),mean(c))

A ce que j'ai compris je dois utiliser la fonction predict mais ensuite.....
je pensais calculer logit(p) en remplaçant par les valeurs (mean(a),mean(b),mean(c)), ça me donnerait ensuite la valeur de p que je mettrais dans predict????
:(

Code : Tout sélectionner

p=c(10,33,20,9,6,56,39,24,28,31,9,62,13,5,21,37,63,72,19,31,43,36,58,67,29,79,37,69,4,22,23,44,31,61,43,7,43,83,81,42,23,61, 37,47,16,60,14,21,49,49,88,94,81,23)
n=c(59,92,57,32,13,95,61,35,41,46,13,91,18,6,31,49,94,91,23,40,52,49,71,86,40,92,44,82,5,25,27,53,41,70,48,8,48,90,91,43,25,67,40,51,17,61,15,23,52,51,95,99,83,24)
a=c(1,3,4,4,4,5,5,6,6,7,7,7,7,7,7,7,7,8,8,8,8,8,9,9,9,9,10,10,10, 11,11,11,11,11,12,12,12,12,12, 13,13,13, 14,15,15,15,15,15,15,16,17,18,18,21)
b=c(5,9,9,9,11,13,14,14,14, 15,15,15,15,16,16,16,16,16,17,18,18,18,19,19,19,20,20,20,20,21,21,21, 22,22,22,22, 22,23,24,24,24,24,24,25,26,26,26,26,29,29,29,30,32,37)
c=c(21,21,21,23,23,23,23,23,23,24,24,24,24,25,25,25,25,25,26,26,27,28,29,29,30,30,30,30,31,31,31,32,32,32,32,32,33,33,34, 35,35,35,35, 36,36,37,38,38, 39,40,42,42,45,52)

reglog <- data.frame(p/n,log(a),log(b),log(c))
model0 <- glm( (p/n) ~ log(a) + log(b) + log(c) + log(a):log(c),reglog,family=binomial,weights=n)

[Renaud Lancelot]

Exemple:

Code : Tout sélectionner

> p <- c(10, 33, 20, 9, 6, 56, 39, 24, 28, 31, 9, 62, 13, 5, 21, 37, 63, 72, 19,
+        31, 43, 36, 58, 67, 29, 79, 37, 69, 4, 22, 23, 44, 31, 61, 43, 7, 43, 83,
+        81, 42, 23, 61, 37, 47, 16, 60, 14, 21, 49, 49, 88, 94, 81, 23)
>
> n <- c(59, 92, 57, 32, 13, 95, 61, 35, 41, 46, 13, 91, 18, 6, 31, 49, 94, 91,
+        23, 40, 52, 49, 71, 86, 40, 92, 44, 82, 5, 25, 27, 53, 41, 70, 48, 8, 48,
+        90, 91, 43, 25, 67, 40, 51, 17, 61, 15, 23, 52, 51, 95, 99, 83, 24)
>
> a <- c(1, 3, 4, 4, 4, 5, 5, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 9, 9,
+        9, 9, 10, 10, 10, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, 13, 14,
+        15, 15, 15, 15, 15, 15, 16, 17, 18, 18, 21)
>       
> b <- c(5, 9, 9, 9, 11, 13, 14, 14, 14, 15, 15, 15, 15, 16, 16, 16, 16, 16, 17,
+        18, 18, 18, 19, 19, 19, 20, 20, 20, 20, 21, 21, 21, 22, 22, 22, 22, 22,
+        23, 24, 24, 24, 24, 24, 25, 26, 26, 26, 26, 29, 29, 29, 30, 32, 37)
>       
> c <- c(21, 21, 21, 23, 23, 23, 23, 23, 23, 24, 24, 24, 24, 25, 25, 25, 25, 25,
+        26, 26, 27, 28, 29, 29, 30, 30, 30, 30, 31, 31, 31, 32, 32, 32, 32, 32,
+        33, 33, 34, 35, 35, 35, 35, 36, 36, 37, 38, 38, 39, 40, 42, 42, 45, 52)
>       
> dfr <- data.frame(p, n, a, b, c)
> rm(p, n, a, b, c)
>
> ## Modèle
> fm <- glm(cbind(p, n - p) ~ log(a) + log(b) + log(c) + log(a):log(c),
+           data = dfr, family = binomial)
>
> ## Valeurs pour lesquelles on veut une prédiction
> New <- data.frame(a = mean(dfr$a),
+                   b = mean(dfr$b),
+                   c = mean(dfr$c))
>
> ## valeur prédite et écart-type
> pred <- predict(fm, newdata = New, type = "response", se = TRUE)
>
> ## IC 95% (approximation de la loi binomiale par la loi normale)
> IC95 <- pred$fit + 1.96 * c(-1, 1) * pred$se.fit
> names(IC95) <- c("b.inf", "b.sup")
> round(IC95, 2)
b.inf b.sup
 0.81  0.86

Renaud

[Rollot Fabien]

Merci, j'ai enfin compris comment on fabrique un intervalle de confiance avec R, tout est plus clair. Je savais bien que la fonction predict servirait.
:)
Fabien R

[Renaud Lancelot]

Attention, ce n'est pas la seule façon de procéder: il y en a beaucoup d'autres ==> revenir à la théorie au besoin.

Renaud

[Logez Maxime]

Bonjour,

Une seconde façon de faire est calculer les bornes inférieures et supérieures dans l'espace du lien puis de revenir dans l'espace de la réponse. C'est une solution que l'on peut retrouver dans le livre : Modelling Binary Data, Second Edition, pages 99-100 mais aussi dans Applied logistic regression pages 40-42.

Code : Tout sélectionner

require(MASS)
pred <- predict(model0,type="link",se.fit=TRUE)
X <- model.matrix(model0)
W <- diag(model0$weights)
variance <- diag(X%*%ginv(t(X)%*%W%*%X)%*%t(X))
ecarttype <- sqrt(variance)

head(ecarttype)
         1          2          3          4          5          6
0.32598405 0.12724915 0.17939431 0.23235984 0.09565500 0.08167962

head(pred$se.fit)
         1          2          3          4          5          6
0.32598405 0.12724915 0.17939431 0.23235984 0.09565500 0.08167962

bornes_lien <- sweep(matrix(pred$se.fit,ncol=1)%*%matrix(c(-1,1),ncol=2)*qnorm(0.975),1,pred$fit,"+")

bornes_rep <- 1/(1+exp(-bornes_lien))

head(cbind(predict(model0,type="response"),bornes_rep))
       [,1]       [,2]      [,3]
1 0.1206475 0.06753225 0.2062965
2 0.3843358 0.32726334 0.4447808
3 0.3691831 0.29166129 0.4541005
4 0.3407518 0.24687244 0.4490453
5 0.4853694 0.43880347 0.5321907
6 0.6034310 0.56456071 0.6410377


Maxime

[Renaud Lancelot]

Oui, mais on peut faire bcp plus simple question code en utilisant l'argument type = "link" (avec se = TRUE) dans la fonction predict.

Renaud

[Logez Maxime]

Je me suis bien servi de se.fit et de la fonction predict

pred <- predict(model0,type="link",se.fit=TRUE)
weep(matrix(pred$se.fit,ncol=1)%*%matrix(c(-1,1),ncol=2)*qnorm(0.975),1,pred$fit,"+")

Maxime

[Renaud Lancelot]

Oui désolé, j'ai lu trop vite.

Renaud