glmmPQL

Questions sur les fonctions statistiques de R

Modérateur : Groupe des modérateurs

Alexandre Villers
Messages : 3
Enregistré le : 25 Oct 2006, 06:04
Contact :

glmmPQL

Messagepar Alexandre Villers » 26 Oct 2006, 05:40

Bonjour,

Je cherche à tester une variable binomiale V (succès ou échec) en fonction d'un traitement T (1 ou 2) et en tenant compte du fait que les couples C testés sont présents dans des proportions diverses;

J'ai commencé avec un modèle de type: modele=glm(Variable~Couples+Traitement, family=binomial) suivi de: anova(modele, test="Chisq")...

Toutefois, on m'a conseillé de tenir compte de la variabilité individuelle avec un modèle à effets mixtes.
On m'a orienté vers le glmmPQL du package MASS après avoir chargé MASS et nlme

> modele=glmmPQL(Variable~Traitement, random=~1|Couples, family=binomial)
iteration 1
iteration 2
iteration 3
iteration 4

> anova(modele)
Erreur dans anova.glmmPQL(modele) : 'anova' is not available for PQL fits
>

Je ne comprends car j'ai vu sur ce forum des anova réalisées sur des glmmPQL.
Ma question est double:
- comment faire pour faire tourner cette satanée ANOVA ??? le but étant d'être sûr que le nombre de tests différent par couple ne masque pas l'effet du traitement;
-si ce n'est pas possible, puis-je me contenter de mon glm tout simple ?

Ces questions pourront paraître naïves et il y a de fortes chances pour qu'elles le soient...


Par avance merci

Alexandre

Pierre Bady
Messages : 405
Enregistré le : 02 Mai 2006, 07:46

Messagepar Pierre Bady » 26 Oct 2006, 07:51

bonjour,

les auteurs de la fonction glmmPQL ont ajouté cette fonction dans les dernières versions de MASS (notamment la version 7.2-29, date : 2006-09-05)

Code : Tout sélectionner

anova.glmmPQL <- function(objet, ...)
stop("'anova' is not available for PQL fits")


donc il me semble qu'ils ont voulu bloquer l'utilisation de la fonction anova pour les modeles glmmPQL. Dans les anciennes versions, on utilisait la fonction anova.lme.

voici une citation de Renaud qui peut expliquer également le "blocage" de cette fonction.

glmmPQL est un "bricolage" qui utilise la fonction lme pour réaliser l'essentiel des calculs. Son auteur (B Ripley) n'a pas écrit de méthode particulière pour les modèles ajustés avec glmmPQL. Je serais circonspect sur les résultats obtenus avec anova. Note d'ailleurs qu'avec les objets de classe lme, il faut utiliser l'argument type = "m" pour avoir un test de Wald correct ==> à essayer dans ton exemple pour voir si ça change qque chose. Pour plus de détails sur glmmPQL, voir la v4 du bouquin MASS (Venables et Ripley).

la citation est issue de viewtopic.php?t=164 .

sinon pour construire votre modèle, vous avez aussi la possibilité d'utiliser les fonctions lmer, glmmML, etc...

en espérant avoir aider un peu :D

@+++


Pierre

PS: ça sent le sapin pour cette fonction (c'est juste un avis perso)
=@===--------¬-------¬------¬-----¬
liens utiles :
http://www.gnurou.org/Writing/SmartQuestionsFr
http://neogrifter.free.fr/welcomeOnInternet.jpg
]<((((*< -------------------------------

Alexandre Villers
Messages : 3
Enregistré le : 25 Oct 2006, 06:04
Contact :

glmmPQL

Messagepar Alexandre Villers » 26 Oct 2006, 10:07

Merci pour votre réponse...
Le seul problème, du moins c'est ce que j'ai cru voir en allant sur la doc d'autres glmm, en particulier glmmML, c'est que soit je ne peux pas utiliser family=binomial, soit je ne pas mettre de random...

Mais comme je suis un novice en glm et plus, peut-être que je passe complètement à côté d'une autre façon de tester ces données.

En tous les cas, cela confirme bien ce que j'avais lu..


Encore merci

Alexandre

Pierre Bady
Messages : 405
Enregistré le : 02 Mai 2006, 07:46

Messagepar Pierre Bady » 26 Oct 2006, 12:29

rebonjour,


Le seul problème, du moins c'est ce que j'ai cru voir en allant sur la doc d'autres glmm, en particulier glmmML, c'est que soit je ne peux pas utiliser family=binomial, soit je ne pas mettre de random...


je ne comprends pas le problème (?).

Code : Tout sélectionner

require(glmmML)
require(lme4)
require(MASS)
?glmmML
?lmer
?glmmPQL
id <- factor(rep(1:20, rep(5, 20)))
y <- rbinom(100, prob = rep(runif(20), rep(5, 20)), size = 1)
x <- rnorm(100)
dat <- data.frame(y = y, x = x, id = id)
glmmML(y ~ x, data = dat, cluster = id,family=binomial)
lmer(y ~ x+(1|id), data=dat,family=binomial,method="Laplace")
summary(glmmPQL(y ~ x,random=~1|id, data=dat,family=binomial))
lmer(y ~ x+(1|id), data=dat,family=binomial,method="PQL")


Mais comme je suis un novice en glm et plus, peut-être que je passe complètement à côté d'une autre façon de tester ces données.


juste une petite remarque (faut pas le prendre mal),
on teste pas des données, on teste des hypothèses.

[edit] il est vrai que les vraisemblances sont approchées. donc ma remarque ci-dessous n'est pas valable (désolé ... j'ai répondu trop vite ) :')

sinon, tu peux utiliser le "rapport des vraisemblances pour évaluer la signification stat d'un effet fixe" ( viewtopic.php?t=164 ).


en espérant avoir aider un peu :D

@+++


Pierre
=@===--------¬-------¬------¬-----¬

liens utiles :

http://www.gnurou.org/Writing/SmartQuestionsFr

http://neogrifter.free.fr/welcomeOnInternet.jpg

]<((((*< -------------------------------

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 26 Oct 2006, 23:08

Difficile de répondre sans avoir un échantillon des données.

L'idée derrière l'utilisation d'un modèle mixte est que les couples ont été tirés au sort dans une pop de couples et qu'on souhaite faire une inférence sur la variabilité liée aux couples. Sous réserve de vérifier comment se présentent les données, la traduction de votre modèle en GLMM ajusté avec la fct lmer du package lme4 serait:

Code : Tout sélectionner

library(lme4)
fm1 <- lmer(Variable ~ Traitement + (1 | Couples), family = binomial)


ou avec la fct glmmML du package éponyme:

Code : Tout sélectionner

library(glmmML)
fm2 <- glmmML(Variable ~ Traitement, family = binomial, cluster = Couples)

Avec ces fonctions, on ne peut pas, en toute rigueur, utiliser le test du rapport des vraisemblances pour tester des effets fixes car les vraisemblances ne sont qu'approchées par les méthodes employées. Il reste possible d'utiliser le test de Wald: voir par exemple la fonction wald.test dans le package aod.

A noter que le package lme4 propose la fct mcmcsamp qui permet, à partir d'un modèle ajusté avec lmer, d'estimer ses paramètres par MCMC, puis de réaliser à l'aide des valeurs simulées n'importe quel test ou intervalle de confiance des paramètres.

Renaud

Alexandre Villers
Messages : 3
Enregistré le : 25 Oct 2006, 06:04
Contact :

glmmPQL

Messagepar Alexandre Villers » 27 Oct 2006, 08:13

Merci pour votre aide et vos lumières.
Je ne sais pas si j'ai tout tout saisi mais je vais y travailler

Cordialement

Alexandre

guillaume souchay
Messages : 9
Enregistré le : 20 Avr 2009, 07:36
Contact :

Messagepar guillaume souchay » 20 Avr 2009, 07:58

Bonjour,

j'ai un petit souci avec la fonction glmmPQL pour une régression logistique avec un modèle mixte à mesures répétées.
On cherche à savoir quels sont les effets qui permettent de sexer de facon certaine un oiseau lors d'une observation. La variable Y est coherence et est de type binaire: 1 l'observation du sexe correspond au sexe génétique, 0 sinon.
Les données sont des observations d'oiseaux, chaque ligne correspond à 1 observation, le sexe génétique (sexdna), l'âge lors de l'obs (ageobs), le type d'observation (si l'oiseau était seul ou dans un groupe, typeobs) et le code de l'oiseau sont indiquées comme suit:

ring sexdna coherence typeobs ageobs
CCAF F 1 single 6
CCAF F 1 single 7
CCAJ M 1 single 2.5
CCAL F 1 single 6
CCAL F 1 single 6
CCAL F 1 single 9.5



sexdna et type obs ont été mis en tant que facteur

sex<-as.factor(sexdna)
type<-as.factor(typeobs)


dans ce type de données, j'ai mis l'age lors de l'observation en tant qu'effet aléatoire et celui-ci est hiérarchisé dans l'individu.
mon modèle est donc :
coherence ~ µ + effet fixe sex + effet fixe type + effet aléatoire ring/ageobs

on m'a conseillé au départ d'utiliser la fonction nlme mais je n'ai pas trouvé le moyen de faire une régression logistique avec cette fonction.

je suis donc parti sur:
test<- glmmPQL(coherence~1+sex+type, random=~1|ring/ageobs,family=binomial)

cependant, je n'ai pas l'impression que cela fonctionne car je n'ai pas de valeur d'AIC:

test<- glmmPQL(coherence~1+sex+type, random=~1|ring/ageobs,family=binomial)
Loading required package: nlme
iteration 1
iteration 2
iteration 3
iteration 4
iteration 5
iteration 6
iteration 7
iteration 8
iteration 9
iteration 10

> summary(test)
Linear mixed-effects model fit by maximum likelihood
Data: NULL
AIC BIC logLik
NA NA NA

Random effects:
Formula: ~1 | ring
(Intercept)
StdDev: 2.301666

Formula: ~1 | ageobs %in% ring
(Intercept) Residual
StdDev: 2.451177 0.4209615

Variance function:
Structure: fixed weights
Formula: ~invwt
Fixed effects: coherence ~ 1 + sex + type
Value Std.Error DF t-value p-value
(Intercept) 2.3822675 0.4673510 680 5.097384 0.0000
sexM -0.2002161 0.2921041 680 -0.685427 0.4933
typesingle 0.3260595 0.4451828 286 0.732417 0.4645
Correlation:
(Intr) sexM
sexM -0.279
typesingle -0.899 -0.030

Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-4.8847171 0.2048140 0.2836963 0.3462454 5.0149940

Number of Observations: 1389
Number of Groups:
ring ageobs %in% ring
682 1102




Si quelqu'un a une idée pour comprendre pourquoi l'AIC est NA.
S'il n'y a pas d'effet de mes facteurs, c'est un résultat mais en l'état, je pense qu'il y a un souci dans le modèle (modèle erroné ou calcul du maximum de vraisemblance impossible?)
"There is no true model"

Anderson & Burnham 1999


Retourner vers « Archives : Fonctions statistiques »

Qui est en ligne

Utilisateurs parcourant ce forum : Aucun utilisateur enregistré et 1 invité