Régression logistique

[Rollot Fabien]

Bonjour,
j'ai quelques interrogations pour faire une régression logistique.
J'ai un jeu de données et je veux expliquer la probabilité de décès par rapport à des concentrations de produits. J'ai 4 colonnes, une avec les proba de décès,et trois autres représentants les concentration de différents produit.
Je souhaite établir une régression logistique mais je ne sais pas quel fonction prendre entre lrm du package Design et entre glm.
Je ne vois pas trop la différence entre les deux fonctions.
Si quelqu'un pourrait m'aider à ce sujet
Merci

[Nicolas Péru]

Pour ma part, je conseillerais la fonction glm tout simplement parce que de nombreux outil de diagnostic existe via les fonctions assiociées à glm. Je ne sais pas ce qu'il en est pour lrm mais tous les stateux que j'ai pu croiser utilise glm.

Nicolas

[Rollot Fabien]

ok merci,
J'aurais une autre question à ce sujet, en cours nous avons vu la régression logistique binaire mais là je ne pense pas qu'elle le soit ( on veut expliquer la probabilité de décès), donc family=binomial dans gjms ne marche pas?
Peut-etre utilisé l'option weights?

[Nicolas Péru]

Si la réponse est bien binaire :
-0 : survie
-1 : décès

proba de décès = nb de décès/nb de patient

Dans le cas où on a la proba directement, il faut effectivement indiquer dans l'arguments weights le n correspondant à chaque valeur.

[Renaud Lancelot]

ou utiliser la syntaxe

Code : Tout sélectionner

cbind(succes, echec) ~ x1 + x2...

pour la formule, ou dans cet exemple

Code : Tout sélectionner

cbind(deces, n - deces) ~ x1 + x2...

Renaud

[Rollot Fabien]

Renaud je ne comprend pas ta réponse. Où dois je mettre ce code? J'utilise quand même la fonction glm avec l'option weights?

[Renaud Lancelot]

Exemple:

Code : Tout sélectionner

> set.seed(12321)
> n <- rep(20, 10)
> y <- rbinom(10, 20, .2)
> x <- factor(sample(c("A", "B"), 10, replace = TRUE))
>
> glm(cbind(y, n - y) ~ x, family = binomial)

Call:  glm(formula = cbind(y, n - y) ~ x, family = binomial)

Coefficients:
(Intercept)           xB 
    -1.1371      -0.4724 

Degrees of Freedom: 9 Total (i.e. Null);  8 Residual
Null Deviance:      16.24
Residual Deviance: 14.76        AIC: 48.05
> glm(y/n ~ x, weights = n, family = binomial)

Call:  glm(formula = y/n ~ x, family = binomial, weights = n)

Coefficients:
(Intercept)           xB 
    -1.1371      -0.4724 

Degrees of Freedom: 9 Total (i.e. Null);  8 Residual
Null Deviance:      16.24
Residual Deviance: 14.76        AIC: 48.05


Renaud

[Rollot Fabien]

Ok, donc si je comprend bien la commande cbind remplace en quelque sorte l'option weights?
Je pensais faire comme suis pour expliquer la probabilité de décès,avec p le nombre de décés, n le nombre total d'insecte pour l'essai et x1,x2,x3 des concentrations de produits.
Dois-je utiliser cbind alors? j'espere que vous arriverez à lire!

Code : Tout sélectionner

> p=c(10,33,20,9,6,56,39,24)
> n=c(59,92,57,32,13,95,61,35)
> y=p/n
> x1=c(1,3,4,4,4,5,5,6)
> x2=c(5,9,9,9,11,13,14,14)
> x3=c(21,21,21,23,23,23,23,23)
> y=factor(y)
> X1=log(x1)
> X2=log(x2)
> X3=log(x3)
> reglog=data.frame(y,X1,X2,X3)
> model=glm(y~(X1+X2+X3)^3,reglog,family=binomial,weights=n)
> model

Call:  glm(formula = y ~ (X1 + X2 + X3)^3, family = binomial, data = reglog,      weights = n)

Coefficients:
(Intercept)           X1           X2           X3        X1:X2        X1:X3        X2:X3     X1:X2:X3 
    -291.07      -849.83       776.27        32.35       -58.82       321.05      -221.54           NA 

Degrees of Freedom: 7 Total (i.e. Null);  1 Residual
Null Deviance:      347.9
Residual Deviance: 3.063e-10    AIC: 14

[Renaud Lancelot]

Ok, donc si je comprend bien la commande cbind remplace en quelque sorte l'option weights?

Oui. Dans les cas extrêmes, rien n'interdit de lire l'aide de la fct glm ;-) Voir les arguments formula et weights, ainsi que la section "details".

Je pensais faire comme suis pour expliquer la probabilité de décès,avec p le nombre de décés, n le nombre total d'insecte pour l'essai et x1,x2,x3 des concentrations de produits.
Dois-je utiliser cbind alors?

Les deux syntaxes sont équivalentes.

j'espere que vous arriverez à lire!

Code : Tout sélectionner

> p=c(10,33,20,9,6,56,39,24)
> n=c(59,92,57,32,13,95,61,35)
> y=p/n
> x1=c(1,3,4,4,4,5,5,6)
> x2=c(5,9,9,9,11,13,14,14)
> x3=c(21,21,21,23,23,23,23,23)
> y=factor(y)
> X1=log(x1)
> X2=log(x2)
> X3=log(x3)
> reglog=data.frame(y,X1,X2,X3)
> model=glm(y~(X1+X2+X3)^3,reglog,family=binomial,weights=n)
> model

Call:  glm(formula = y ~ (X1 + X2 + X3)^3, family = binomial, data = reglog,      weights = n)

Coefficients:
(Intercept)           X1           X2           X3        X1:X2        X1:X3        X2:X3     X1:X2:X3 
    -291.07      -849.83       776.27        32.35       -58.82       321.05      -221.54           NA 

Degrees of Freedom: 7 Total (i.e. Null);  1 Residual
Null Deviance:      347.9
Residual Deviance: 3.063e-10    AIC: 14

Apparemment, votre modèle est trop complexe pour les données dont vous disposez. Peut-être commencer par un modèle additif et ajouter des termes d'interaction si nécessaire, selon exploration graphique ou votre connaissance de la biologie des insectes ou du mécanisme d'action de l'insecticide.

Renaud

[Rollot Fabien]

Oui oui, justement ma mission est de trouver le meilleur modèle qui soit, de plus j'ai rétréci mes données car j'en ai beaucoup plus!
Mais je ne comprend pas pourquoi ça ne me donne pas les mêmes résultats lorsque j'essaye avec cbind:

Code : Tout sélectionner

> p=c(10,33,20,9,6,56,39,24)
> n=c(59,92,57,32,13,95,61,35)
> y=p/n
> x1=c(1,3,4,4,4,5,5,6)
> x2=c(5,9,9,9,11,13,14,14)
> x3=c(21,21,21,23,23,23,23,23)
> y=factor(y)
> X1=log(x1)
> X2=log(x2)
> X3=log(x3)
> reglog=data.frame(y,X1,X2,X3)
> model=glm(cbind(p,n-p)~(X1+X2+X3)^3,reglog,family=binomial)
> model

Call:  glm(formula = cbind(p, n - p) ~ (X1 + X2 + X3)^3, family = binomial,      data = reglog)

Coefficients:
(Intercept)           X1           X2           X3        X1:X2        X1:X3        X2:X3     X1:X2:X3 
    -22.349      -14.896       21.958        5.795        1.209        3.981       -6.576           NA 

Degrees of Freedom: 7 Total (i.e. Null);  1 Residual
Null Deviance:      54.27
Residual Deviance: 0.09681      AIC: 47.5

[Renaud Lancelot]

Je vous conseille:

1. De nettoyer votre espace de travail: ne garder que le data.frame et pas les variables isolées plus le data.frame
2. De mettre donc toutes les variables dans le data.frame, y compris p ce qui n'est pas le cas pour l'instant
3. D'utiliser "<-" et non "=" comme opérateur d'affectation
4. D'aérer votre code (espaces, indentations...) pour faciliter sa lecture et le débuggage
5. de nommer les arguments dans l'appel de la fct car quand vous ne les nommez pas, la fct les interprète comme venant dans l'ordre dans lequel ils sont programmés:

glm(formula, family = gaussian, data, weights, subset,
na.action, start = NULL, etastart, mustart,
offset, control = glm.control(...), model = TRUE,
method = "glm.fit", x = FALSE, y = TRUE, contrasts = NULL,
...)

Par exemple, data est le 3e argument alors que dans votre code, vous le donnez (sans le nommer) en second argument. Apparemment ça se passe bien ici, mais c'est une source de souci potentiel.

Renaud

[Rollot Fabien]

ok donc si je comprend bien ce que vous me dites, j'ai refait mon programme et en utilisant cbind et weights j'obtiens la même chose:

Code : Tout sélectionner

> p = c(10,33,20,9,6,56,39,24)
> n = c(59,92,57,32,13,95,61,35)
> x1 = c(1,3,4,4,4,5,5,6)
> x2 = c(5,9,9,9,11,13,14,14)
> x3 = c(21,21,21,23,23,23,23,23)
> X1 = log(x1)
> X2 = log(x2)
> X3 = log(x3)
> reglog <- data.frame(p,n,X1,X2,X3)
> model <- glm((p/n)~(X1+X2+X3)^3,data=reglog,family=binomial,weights=n)
> model1 <- glm(cbind(p,n-p)~(X1+X2+X3)^3,data=reglog,family=binomial)
> model

Call:  glm(formula = (p/n) ~ (X1 + X2 + X3)^3, family = binomial, data = reglog,      weights = n)

Coefficients:
(Intercept)           X1           X2           X3        X1:X2        X1:X3        X2:X3     X1:X2:X3 
    -22.349      -14.896       21.958        5.795        1.209        3.981       -6.576           NA 

Degrees of Freedom: 7 Total (i.e. Null);  1 Residual
Null Deviance:      54.27
Residual Deviance: 0.09681      AIC: 47.5
> model1

Call:  glm(formula = cbind(p, n - p) ~ (X1 + X2 + X3)^3, family = binomial,      data = reglog)

Coefficients:
(Intercept)           X1           X2           X3        X1:X2        X1:X3        X2:X3     X1:X2:X3 
    -22.349      -14.896       21.958        5.795        1.209        3.981       -6.576           NA 

Degrees of Freedom: 7 Total (i.e. Null);  1 Residual
Null Deviance:      54.27
Residual Deviance: 0.09681      AIC: 47.5


Mais je ne vois pas où ma régression est binaire dans le cas model (premier cas)??

[Nicolas Péru]

La proba fournie est bel et bien obtenue à partir d'une La variable aléatoire de Bernoulli (on déduit la proba d'une variable qui peut prendre comme valeur 0 ou 1)
Donc dans les deux cas on a la même chose...c'est ce que j'indiquais dans mon second post. Et c'est entre autre pour cela qu'on indique les effectifs dans l'argument weight : afin de pouvoir déduire le nb de 0 ou le
nb de 1.
Avec la commande cbind on crée une matrice dont la première colonne est le nombre de succès et la seconde le nombre d'échec (ou inversement...j'ai un doute :p).

Nicolas

[Rollot Fabien]

ok je vous remercie, les deux méthodes donnent bien les mêmes résultats, c'est rassurant.
A moi de trouver le meilleur modèle via la commande step.
Merci