messages d'avertissement dans lme4 pour l'analyse de survie qui ne s'est pas produite il y a 3 ans

Question

J'essaye d'adapter un modèle linéaire d'effets mixtes généralisé à mes données, en utilisant le paquet lme4.

Les données peuvent être décrites comme suit (voir l'exemple ci-dessous): Données de survie du poisson sur 28 jours. Les variables explicatives dans l'exemple de jeu de données sont:

• Region C'est la région géographique d'où proviennent les larves.
• treatment Les températures auxquelles les sous-échantillons de poissons de chaque région ont été relevés.
• replicate Une des trois réplications de l'expérience entière
• tub Variable aléatoire. 15 cuves (utilisées pour maintenir les températures expérimentales dans les aquariums) au total (3 replicate s pour chacune des 5 températures treatment s). Chaque bac contenait 1 aquarium pour chaque Region (4 aquariums au total) et était situé au hasard dans le laboratoire.
• Day Explicite, nombre de jours depuis le début de l'expérience.
• stage n'est pas utilisé dans l'analyse. Peut être ignoré.

variable de réponse

• csns survie cumulative. c'est-à-dire remaining fish/initial fish at day 0.
• start poids utilisés pour indiquer au modèle que la probabilité de survie est relative à ce nombre de poissons au début de l'expérience.
• aquarium Deuxième variable aléatoire. C'est l'identifiant unique de chaque aquarium contenant la valeur de chaque facteur auquel il appartient. par exemple. N-14-1 signifie Region N, Treatment 14, replicate 1.

Mon problème est inhabituel dans la mesure où j'installé le modèle suivant avant:

dat.asr3<-glmer(csns~treatment+Day+Region+ 
     treatment*Region+Day*Region+Day*treatment*Region+ 
     (1|tub)+(1|aquarium),weights=start, 
     family=binomial, data=data2) 

Cependant, maintenant que je tente de relancer le modèle, pour générer des analyses pour la publication, je suis obtenir les erreurs suivantes avec la même structure de modèle et le même package. La sortie est listé ci-dessous:

> Warning messages: 

1: Dans eval (expr, envir, Enclos): #successes non entières dans un GLM binomial!
2: Dans checkConv (attr (opt, "derivs"), optez $ par, ctrl = contrôle $ checkConv,:
modèle n'a pas réussi à converger avec max | grad | = 1,59882 (tol = 0,001, composante> 1)
3: Dans checkConv (attr (opt, "derivs"), optez $ par, ctrl = contrôle $ checkConv,:
modèle est presque non identifiable: très grande valeur propre
- Rescale des variables, le modèle est presque non identifiable: grande valeur propre ratio
- Variables de remise à zéro?

Ma compréhension est la suivante:

message d'avertissement 1.

non-integer #success in a binomial glm fait référence au format de la proportion de la variable csns. J'ai consulté plusieurs sources, ici incluses, github, r-help, etc, et tous suggéré cela. Le chercheur qui m'a aidé dans cette analyse il y a trois ans est inaccessible. Cela peut-il avoir à voir avec les changements dans le paquet lme4 au cours des 3 dernières années?

message d'avertissement 2.

Je comprends que c'est un problème, car il y a des points de données insuffisantes pour adapter le modèle, en particulier à
L-30-1, L-30-2 et L-30-3,
où seulement deux observations sont faits:
Day 0 csns=1.00 et Day 1 csns=0.00
pour les trois aquariums. Par conséquent, il n'y a pas de variabilité ou de données suffisantes pour adapter le modèle.

Néanmoins, ce modèle dans lme4 a fonctionné avant, mais ne fonctionne pas sans ces avertissements maintenant.

Message d'avertissement 3

Celui-ci est tout à fait inconnu pour moi. Je ne l'ai jamais vu auparavant.

données Exemple:

Region treatment replicate tub Day stage csns start aquarium 
N  14   1 13 0 1 1.00 107 N-14-1 
N  14   1 13 1 1 1.00 107 N-14-1 
N  14   1 13 2 1 0.99 107 N-14-1 
N  14   1 13 3 1 0.99 107 N-14-1 
N  14   1 13 4 1 0.99 107 N-14-1 
N  14   1 13 5 1 0.99 107 N-14-1 

Les données en question 1005cs.csv est ici par nous transférons: http://we.tl/ObRKH0owZb

Toute aide à déchiffrer ce problème, serait grandement apprécié. En outre, toutes les suggestions alternatives pour des paquets ou des méthodes appropriés pour analyser ces données seraient également utiles.

Answer 1

tl; dr l'avertissement "succès non-entiers" est précis; C'est à vous de décider si l'ajustement d'un modèle binomial à ces données a vraiment du sens. Les autres avertissements suggèrent que l'ajustement est un peu instable, mais la mise à l'échelle et le centrage de certaines variables d'entrée peuvent faire disparaître les avertissements. Il est à vous, encore une fois, de décider si les résultats de ces différentes formulations sont assez différentes pour vous inquiéter ...

data2 <- read.csv("1005cs.csv") 
library(lme4) 

modèle Fit (formulation du modèle légèrement plus compact)

dat.asr3<-glmer(
    csns~Day*Region*treatment+ 
     (1|tub)+(1|aquarium), 
      weights=start, family=binomial, data=data2) 

Je reçois les avertissements que vous signalez.

Jetons un coup d'oeil à les données:

library(ggplot2); theme_set(theme_bw()) 
ggplot(data2,aes(Day,csns,colour=factor(treatment)))+ 
     geom_point(aes(size=start),alpha=0.5)+facet_wrap(~Region) 

Rien évidemment problématique ici, même si il montre clairement que les données sont très proches de 1 pour certaines combinaisons de traitement, et que la les valeurs de traitement sont loin de zéro. Essayons d'échelle & centrage certaines des variables d'entrée:

data2sc <- transform(data2, 
        Day=scale(Day), 
        treatment=scale(treatment)) 
dat.asr3sc <- update(dat.asr3,data=data2sc) 

Maintenant, l'avertissement « très grande valeur propre » est parti, mais nous avons encore l'avertissement « # non entiers succès », et un max | grad | = 0,082. Essayons un autre optimiseur:

dat.asr3scbobyqa <- update(dat.asr3sc, 
        control=glmerControl(optimizer="bobyqa")) 

Maintenant, seul l'avertissement "non-integer #successes" subsiste.

d1 <- deviance(dat.asr3) 
d2 <- deviance(dat.asr3sc) 
d3 <- deviance(dat.asr3scbobyqa) 
c(d1,d2,d3) 
## [1] 12597.12 12597.31 12597.56 

Ces déviances ne diffèrent pas de beaucoup (0,44 sur l'échelle de la déviance est plus que pourrait être pris en compte par erreur d'arrondi, mais pas beaucoup de différence dans la qualité d'ajustement); En fait, le premier modèle donne le meilleur (le plus bas) deviance, ce qui suggère que les mises en garde sont des faux positifs ...

resp <- with(data2,csns*start) 
plot(table(resp-floor(resp))) 

Cela montre clairement qu'il ya vraiment des réponses non entières, de sorte que le l'avertissement est correct.