données binaires avec glmer

[margot julien]

Bonjour,

j'ai des données non normales alors j'ai transformé ces données en données binaires (0 = pas de sucre et 1 = du sucre), j'aimerais les analyser, sauf que j'ai des messages d'erreur que j'utilise glm (sans tenir compte des répétitions) ou bien glmer avec ou sans interaction et en tenant compte des répétitions.

ma variable mesurée est : sug (binary data: 0 pour pas de sucre mesuré et 1 si du sucre à été mesuré)
mes effets fixes: temp (2 Temperatures différentes) et var (6 variétés différentes testées) + l'interaction des deux
mes réplicats = rep (le sucre à été mesuré sur 5 réplicats biologiques)

J'ai donc un split plot design et j'aimerais voir les différences variétales, entre températures et pour les interactions au niveau du sucres.

Pouvez-vous m'aider svp? Voici mon script ci-dessous et le csv en format code ci-dessous.
Normalement je dois obtenir un effet de l'interaction temp * var car je sais que j'ai 3 variétés qui ne sucrent pas pour temp = 8 mais qui sucrent pour temp = 4 et pour les 3 autres variétés elles ne sucrent pas pour les 2 températures.

Merci par avance.
Margot

Code : Tout sélectionner

file<-read.csv2("file.csv", na.strings = "NA",header = TRUE)
head(file)
str(file)
file$sug<-as.numeric(file$sug)
file$temp<-as.factor(file$temp)
library(lme4)
library(car)
##Measured variable : sug (binary data: 0 means no sugar and 1 means sugar)
##fixed effect to observed: temp, var and interaction temp, var
##replicates = rep (5 biological replicates per sug analysis)

#split plot design?
   # 2 temperatures (temp)
   # 6 varieties (var)
   #observation of the sugar content (sug) on 5 biological replicates

#model without random factor (rep) / glm function
model_global = glm(sug ~ temp * var , family="binomial", data =file)
Anova(model_global)
    #warning:
        # Warning messages:
               #   1: glm.fit: fitted probabilities numerically 0 or 1 occurred
                # 2: glm.fit: fitted probabilities numerically 0 or 1 occurred
               # 3: glm.fit: fitted probabilities numerically 0 or 1 occurred

#model with random factor (rep) / glmer fuction
model_global_int = glmer(sug ~ temp * var + (1 | rep), family=binomial, data =file)
Anova(model_global_int)

     #Warning messages:
       #   1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
          #                     unable to evaluate scaled gradient
            #                   2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
                                   # Hessian is numerically singular: parameters are not uniquely determined
model_global_plus = glmer(sug ~ temp + var + (1 | rep), family=binomial, data =file)
Anova(model_global_plus)

      #Error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,  :
                     #  (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

VOICI mon fichier au format code:

Code : Tout sélectionner

> dput(file)
structure(list(temp = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L), .Label = c("4", "8"), class = "factor"),
    var = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L,
    3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
    5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L,
    6L, 6L, 6L, 6L, 6L), .Label = c("var1", "var2", "var3", "var4",
    "var5", "var6"), class = "factor"), rep = structure(c(1L,
    2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L,
    2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L,
    2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L,
    2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L), .Label = c("rep1",
    "rep2", "rep3", "rep4", "rep5"), class = "factor"), sug = c(1,
    1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0,
    0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0)), row.names = c(NA, -60L), class = "data.frame")

[Eric Wajnberg]

Votre question relève d'un problème de statistique, pas de l'usage de R.

Si vous tapez :

Code : Tout sélectionner

summary(model_global)

vous allez constater que les pentes sont estimées avec des SE énormes. La raison est que vous avez des probas observées de 0.0 ou de 1.0 (c'est ce que vous disent les warnings), et dans ce cas le logit part à plus ou moins l'infini. Impossible d'estimer la matrice Hessienne pour calculer les variances-covariances des paramètres du modèle que vous cherchez à ajuster à vos données.

D'après McCullagh & Nelder, il faut partir dans ce cas sur une régression quasi-logisitque en ré-écrivant (customisant) la fonction de lien, et c'est possible avec R.

Mais on est clairement dans un problème hors sujet sur ce forum.

HTH, Eric.

[margot julien]

Bonjour Eric,

merci pour le retour.

Je comprends maintenant, pourriez-vous à ce moment là m'indiquer comment faire une régression quasi-logisitque customisée avec R pour pouvoir analyser mes données s'il vous plait?

Merci par avance.

[Eric Wajnberg]

Il y a plusieurs années, j'avais lancé une discussion là-dessus et avais donné la clé pour faire ceci avec R, sur un forum de statistique francophone. Vous trouverez ceci ici.

Il s'agit d'une longue discussion, et la mise en œuvre sous R n'est pas aisée..

HTH, Eric.

[margot julien]

Merci beaucoup pour ces informations, je vais étudier le lien et voir si j'arrive à appliquer cela sur R.
Bonne journée.
Cordialement,
Margot.

[Maxime Hervé]

Bonjour,

j'en profite pour rebondir sur cette question, car je découvre comment le problème peut être géré via la fonction de lien maison que mentionne Éric (merci Éric !). En repartant de l'exemple donné dans le lien d'Éric, on obtient les avertissements suivants :

Code : Tout sélectionner

> elogis <- function (M=1) {
+   linkfun <- function(mu) log((mu+0.5/M)/(1-mu+0.5/M))
+   linkinv <- function(eta) ((1+0.5/M)*exp(eta)-0.5/M)/(1+exp(eta))
+   mu.eta <- function(eta) (1+1/M)/(exp(-eta)+exp(eta)+2)
+   valideta <- function(eta) TRUE
+   link <- paste0("elogis(",M, ")")
+   structure(list(linkfun=linkfun,linkinv=linkinv,mu.eta=mu.eta,
+     valideta=valideta,name=link),class="link-glm")
+ }
> tab <- data.frame(case=letters[1:3],yes=c(25,30,0),no=c(1,0,20))
> tab
  case yes no
1    a  25  1
2    b  30  0
3    c   0 20
> mod <- glm(cbind(yes,no)~case,family=binomial(link=elogis()),data=tab)
Warning messages:
1: le pas a été tronqué : hors limites
2: glm.fit: l'algorithme s'est arrêté à la valeur limite
3: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

Puisque je n'ai jamais pratiqué cette fonction de lien, je suis un peu pris au dépourvu par ces avertissements. Si tu as l'habitude Éric, comment les interprètes-tu ? Les deux premiers sont peut-être liés à l'exemple particulier. Le troisième est plus clair.

Merci,

Maxime

[Eric Wajnberg]

Hello Maxime,

Oui, j'avais vu à l'époque ces warnings, qui ne m'ont pas trop choqué. je pense que les warnings viennent de binomial et pas de la fonction de lien.

A noter que - dans cet exemple - les SE des paramètres sont toutes les mêmes et sont plus fortes que ce qu'on attendrait (sqrt(p*(i-p)/n)) (à transposer dans l'espace logit).

Il se trouve que - dans ce cas - on trouve un meilleur ajustement en supprimant l'intercept :

Code : Tout sélectionner

glm(cbind(yes,no)~case-1,family=binomial(link=elogis()),data=tab)

Si on rajoute un "trace=T" on arrive au constat que :

Il y a toujours trois warnings indiquant que l'algorithme n'a pas convergé, qu'il s'est arrêté sur une borne (effectivement, p=1, modalité 2) et que des proba numériquement égales à 0 ou 1 ont été rencontrées. Le résultat est obtenu en 25 itérations avec de nombreuses divisions du pas. Les estimations semblent ok. La matrice de corrélation des estimateurs est excellente : nulle à la précision d'affichage (summary(mod,cor=T)). Ainsi donc, la matrice X'W^-1X semble bien conditionnée.

Bref, hormis ces warnings, je pense que le modèle est assez correctement ajusté. Je ne sais pas trop quoi répondre d'autre, sinon que je continue à penser qu'il sagit là d'une discussion purement statistique et hors-sujet ici.

J'ai - il y a des années - pas mal bossé sur ce problème, car j'avais un cas de ce type à gérer. Il n'y a guère de plus sur le web sur le traitement de ce genre de cas, ce qui est étonnant car tout ceci est clairement explicité dans le McCullagh & Nelder..

HTH, Eric.

[Maxime Hervé]

Ok merci Éric, ton explication va bien au-delà de l'idée que j'aurais pu me faire par moi-même. J'en retiens, pour faire court, que cette fonction de lien modifiée est une manière propre de traiter ce genre de cas et qu'il ne faut pas se formaliser des avertissements.

Maxime

[Logez Maxime]

Bonjour,

tu as des alternatives a cette fonction de lien du côté du package brglm et brglm2.
Ca peut-être intéressant parce que, je cite :

the resultant estimates and their corresponding standard errors are always finite while the maximum likelihood estimates can be infinite (in situations where complete or quasi separation occurs).

Cordialement,
Maxime

[Eric Wajnberg]

Bonjour,

tu as des alternatives a cette fonction de lien du côté du package brglm et brglm2.
Ca peut-être intéressant parce que, je cite :

the resultant estimates and their corresponding standard errors are always finite while the maximum likelihood estimates can be infinite (in situations where complete or quasi separation occurs).

Cordialement,
Maxime

Merci Maxime,

Je ne connais pas ces packages. Ceci dit (comme ceci est expliqué en détail dans le lien que je donne au début de ce threat), le problème dont nous parlons ici n'est pas (nécessairement) lié à un problème de séparabilité, mais juste au fait qu'il y a des probas estimées à 0.0 ou 1.0 qui conduisent la fonction logit à s'envoler à plus ou moins l'infini.

Cordialement, Eric.