Discussion :

[Margotte02]

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 fichier "fille" en pièce-jointe.
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 :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
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

[faubry]

La valeur de sug est constante pour toutes les combinaisons sauf pour la combinaison var=var5 et temp=4, d'où les messages pour glm puisque pour toutes ces combinaisons la probabilité est 0 ou 1 comme indiqué. Ce qui implique aussi les problèmes de glmer.

[Margotte02]

La valeur de sug est constante pour toutes les combinaisons sauf pour la combinaison var=var5 et temp=4, d'où les messages pour glm puisque pour toutes ces combinaisons la probabilité est 0 ou 1 comme indiqué. Ce qui implique aussi les problèmes de glmer.

merci pour le retour. En effet les résultats sont très répètables entre les 5 réplicats.

Donc par exemple, pour var 4 quand temp = 4 j'ai toujours du sucre (=1) et pour var 4 quand temp = 8 je n'ai jamais de sucre (= 0 pour les 5 replicats), le modèle devrait donc montrer un effet de temp, un effet de var et aussi une interaction var*temp, pourquoi ici ça ne fonctionne pas? et comment faire pour montrer cela?

[faubry]

Pour glmer, c'est un problème de résolution mathématique du modèle : lors de la résolution mathématique, la matrice hessienne est alors singulière ou dégénérée donc, le modèle ne peut pas être résolue. Pour glm, les logit (log(Pr(sug=1)/Pr(sug=0)) valent l'infini en valeur absolue puisque l'une des valeurs (0 ou 1) n'existe pas dans l'échantillon.

Pour vérifier :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

aggregate( sug ~ temp + var, file, mean)

[Margotte02]

Pour glmer, c'est un problème de résolution mathématique du modèle : lors de la résolution mathématique, la matrice hessienne est alors singulière ou dégénérée donc, le modèle ne peut pas être résolue. Pour glm, les logit (log(Pr(sug=1)/Pr(sug=0)) valent l'infini en valeur absolue puisque l'une des valeurs (0 ou 1) n'existe pas dans l'échantillon.

Pour vérifier :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

aggregate( sug ~ temp + var, file, mean)

en effet avec votre formule . ça me dit "objet sug not found". Mais alors comment faire? est-ce que je pourrais arbitrairement remplacer 0 par 10 et 1 par 100 par exemple tout en gardant "family = binomial" dans la formule avec glmer?

[faubry]

en effet avec votre formule . ça me dit "objet sug not found"

Bizarre car moi ça me donne :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
structure(list(temp = c(4L, 8L, 4L, 8L, 4L, 8L, 4L, 8L, 4L, 8L, 
4L, 8L), var = c("var1", "var1", "var2", "var2", "var3", "var3", 
"var4", "var4", "var5", "var5", "var6", "var6"), sug = c(1, 0, 
0, 0, 0, 0, 1, 0, 0.8, 0, 0, 0)), row.names = c(NA, -12L), class = "data.frame")

qui montre bien qu'à part temp=4 et var=var5, sug est constant pour toutes les autres combinaisons.

La variable sug est binaire donc remplacer 0 par 10 et 1 par 100 ne changera strictement rien. Il est impossible de faire une analyse avec l'échantillon. On ne peut observer que des combinaisons qui donne 0 ou 1 et une qui donne plutôt 1.

[Margotte02]

merci pour le retour. On m'a dit que je pouvais aussi faire une régression quasi-logisitque en customisant la fonction de lien, et que c'est possible avec R. Ce sujet a été discuté ici: http://statistiques.forumpro.fr/t630...-ou-p1-0#20671

Je vais voir si j'arrive à faire ça.

[faubry]

Tu auras le même problème en changeant la fonction de lien (par exemple, passer du logit au probit) ou en utilisant du quasi binomial car pour certaines combinaisons de temp et de var (et dans ton cas, toutes sauf une), sug est constant : il n'y a donc pas de variabilité de la variable à expliquer, donc pas de variance d'où les problèmes.

[faubry]

La solution donnée par Eric Wajnberg rentre dans la catégorie des estimateurs robustes. Il existe aussi d'autres solutions dans R, notamment les packages robust (pour glm) et robustlmm (pour glmer) mais aussi glmnet ou penalized, qui utilisent d'autres critères. Il faut savoir que ces estimateurs robustes sont aussi, par nature, biaisés mais, pour faire simple, on peut souvent s'arranger que le biais soit suffisamment faible pour être négligé (à condition de bien définir cette notion de 'suffisamment faible')

Pour ton problème, si tes données sont celles du fichier, je pense qu'il est illusoire de vouloir utiliser des méthodes robustes pour résoudre ton problème, le premier problème étant à mes yeux celui de la taille de l'échantillon : le nombre de 5 répétitions pour chaque combinaison (temp, var) étant plus qu'à la limite de validité des méthodes. En tout cas, dans ce cas, le modèle mixte est à exclure car il augmente le nombre de degrés de liberté du modèle à résoudre et donc son instabilité.