fonction qqmath appliquée sur un objet de classe ranef.lmer

Questions sur les fonctions statistiques de R

Modérateur : Groupe des modérateurs

Tillard
Messages : 87
Enregistré le : 17 Déc 2004, 10:32

fonction qqmath appliquée sur un objet de classe ranef.lmer

Messagepar Tillard » 22 Aoû 2006, 21:34

Merci Renaud pour ta reponse detaillée a propos de la fonction qqmath cachée dans Matrix, il fallait la trouver celle la;

j'ai refait les calculs a la main et ceux ci sont toujours differents de ceux obtenus avec qqmath

un exemple

Code : Tout sélectionner

library(lme4)
clinic <- data.frame(
    clinic = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,,
    1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,
    2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
    2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,
    4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,
    5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,8,8,8),
    trt = c("drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug",
    "drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug",
    "drug","drug","drug","drug","drug","drug","drug","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl",
    "cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl",
    "cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","drug",
    "drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug",
    "drug","drug","drug","drug","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl",
    "cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl",
    "cntl","cntl","cntl","cntl","cntl","cntl","drug","drug","drug","drug","drug","drug","drug","drug","drug",
    "drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","cntl","cntl","cntl","cntl","cntl",
    "cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","drug",
    "drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug",
    "cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl",
    "cntl","cntl","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug",
    "drug","drug","drug","drug","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl","cntl",
    "cntl","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","drug","cntl","cntl","cntl",
    "cntl","cntl","cntl","cntl","cntl","cntl","cntl","drug","drug","drug","drug","drug","cntl","cntl","cntl",
    "cntl","cntl","cntl","cntl","cntl","cntl","drug","drug","drug","drug","drug","drug","cntl","cntl","cntl",
    "cntl","cntl","cntl","cntl"),
    resp =  c(1,1,1,1,1,1,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,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,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,1,1,1,1,1,1,1,1,
    1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,
    0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,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,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,
    1,1,1,1,0,0,1,1,1,1,1,1,0))                                                                                                         

fm1 <- lmer(resp ~ trt + (1|clinic), data=clinic, family=binomial, method = "Laplace")
summary(fm1)
random <- ranef(fm1, postVar=T)
qqmath(random)

#la meme chose mais "à la main"
mat.random <- data.frame(mean = unlist(random[[1]]),
                low = unlist(random[[1]]) - (1.96 * sqrt(c(attr(random[[1]], "postVar")))),
                high = unlist(random[[1]]) + (1.96 * sqrt(c(attr(random[[1]], "postVar")))))
mat.random <- mat.random[order(mat.random$mean),]

plot(1:8, mat.random$mean, ylim=c(-4,4))
for(i in 1:8){
    segments(i, mat.random$low[i], i, mat.random$high[i])}
abline(h = 0)


les intervalles de confiance diffèrent

ne faudrait-il pas légèrement modifier la fonction qqmath de la façon suivante:
voir ligne #a modifier et #modifie

Code : Tout sélectionner

setMethod("qqmath", signature(x = "ranef.lmer"),
          function(x, data, ...) {
              prepanel.ci <- function(x, y, se, subscripts, ...) {
                  y <- as.numeric(y)
                  se <- as.numeric(se[subscripts])
                  hw <- 1.96 * se
                  list(ylim = range(y - hw, y + hw, finite = TRUE))
              }
              panel.ci <- function(x, y, se, subscripts, pch = 16, ...)  {
                  panel.grid(h = -1,v = -1)
                  panel.abline(h = 0)
                  x <- as.numeric(x)
                  y <- as.numeric(y)
                  se <- as.numeric(se[subscripts])
                  ly <- y - 1.96 * se
                  uy <- y + 1.96 * se
                  panel.segments(x, y - 1.96*se, x, y + 1.96 * se,
                                 col = 'black')
                  panel.xyplot(x, y, pch = pch, ...)
              }
              f <- function(x) {
                  if (!is.null(attr(x, "postVar"))) {
               #       require("lattice", quietly = TRUE)
                      pv <- attr(x, "postVar")
                      cols <- 1:(dim(pv)[1])
                      se <- unlist(lapply(cols, function(i) sqrt(pv[i, i, ])))
                      nr <- nrow(x)
                      nc <- ncol(x)
                      ord <- unlist(lapply(x, order)) +
                          rep((0:(nc - 1)) * nr, each = nr)
                      rr <- 1:nr
                      ind <- gl(ncol(x), nrow(x), labels = names(x))
                      xyplot(unlist(x)[ord] ~
                             rep(qnorm((rr - 0.5)/nr), ncol(x)) | ind[ord],
#                     se = se, prepanel = prepanel.ci, panel = panel.ci,            #à modifier
                      se = se[ord], prepanel = prepanel.ci, panel = panel.ci,       #modifié
                             scales = list(y = list(relation = "free")),
                             xlab = "Standard normal quantiles",
                             ylab = NULL, aspect = 1, ...)
                  } else {
                      qqmath(~values|ind, stack(x),
                             scales = list(y = list(relation = "free")),
                             xlab = "Standard normal quantiles",
                             ylab = NULL, ...)
                  }
              }
              lapply(x, f)
          })

qqmath(random)


la, on retrouve les mêmes résultats que ceux obtenus à la main !!!

J'ai également fait tourner la procedure mcmc (fonction mcmcsamp)

Code : Tout sélectionner

samp.fm1 <- mcmcsamp(fm1, n=1000, saveb=T)

library(coda)
summary(samp.fm1)

> summary(samp.fm1)

Iterations = 1:1000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 1000

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                   Mean     SD Naive SE Time-series SE
(Intercept)     -1.1222 0.6689  0.02115        0.03933
trtdrug          0.7414 0.2963  0.00937        0.01509
log(clnc.(In))   0.9186 0.6471  0.02046        0.03541
deviance       294.0704 5.7647  0.18230        0.23163
b.1             -0.1824 0.6764  0.02139        0.03479
b.2              1.7997 0.7154  0.02262        0.04275
b.3              0.8812 0.7237  0.02289        0.04266
b.4             -1.4377 0.7981  0.02524        0.04489
b.5             -0.7042 0.7408  0.02343        0.04511
b.6             -1.8918 1.0352  0.03273        0.08299
b.7             -0.8608 0.8957  0.02832        0.04217
b.8              1.7889 0.8019  0.02536        0.04248

2. Quantiles for each variable:

                   2.5%      25%      50%      75%    97.5%
(Intercept)     -2.4772  -1.5161  -1.1138  -0.7254   0.2769
trtdrug          0.1466   0.5443   0.7282   0.9512   1.3115
log(clnc.(In))  -0.2637   0.4725   0.8965   1.3059   2.2845
deviance       284.4250 289.8622 293.4111 297.5706 306.6304
b.1             -1.5281  -0.5971  -0.1843   0.2404   1.1706
b.2              0.4760   1.3376   1.7706   2.2621   3.2568
b.3             -0.6914   0.4477   0.8840   1.3527   2.3560
b.4             -2.9606  -2.0078  -1.4226  -0.8916   0.1073
b.5             -2.3730  -1.1950  -0.6561  -0.2059   0.7170
b.6             -4.6164  -2.4431  -1.8137  -1.2115  -0.1152
b.7             -2.7983  -1.4526  -0.7782  -0.2633   0.6969
b.8              0.1587   1.2791   1.8146   2.3033   3.3381


est ce que la suite est valide (plot des effets aléatoires et de leur intervalle de confiance obtenu par mcmcsamp)?

Code : Tout sélectionner

res.samp.random <- summary(samp.fm1)$statistics[5:12,1:2]

#plot des résidus standardisés façon "MLWin"
qqnorm(res.samp.random[,1] / res.samp.random[,2])

ci.mcmc <- data.frame(id=1:8,mean = res.samp.random[,1],
            low = res.samp.random[,1] - (1.96 * res.samp.random[,2]),
            high = res.samp.random[,1] + (1.96 * res.samp.random[,2]))
ci.mcmc <- ci.mcmc[order(ci.mcmc$mean),]

#graphe "qqmath" prenant en compte la variance des effets aléatoires obtenues par mcmc
plot(1:8, ci.mcmc$mean, ylim=c(-4,4))
for(i in ci.mcmc$id){
    segments(i, ci.mcmc$low[i], i, ci.mcmc$high[i])}
abline(h = 0)


Désolé pour ce long post !
amicalement
Emmanuel Tillard
UMR ERRC (Elevage des Ruminants en Regions Chaudes)
CIRAD - St PIERRE (La Réunion)
tel: 02 62 49 92 54

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

Messagepar Renaud Lancelot » 23 Aoû 2006, 07:59

Salut Manu,

1. Tu dois avoir raison mais pas le temps de regarder en détail. Peux-tu faire un mail à D Bates pour lui signaler le pb ?

2. Idem pour mcmcsamp, mais qques rques:

* voir les différents posts de D Bates et autres sur la question des effets aléatoires, IC et tests à partir des sorties de mcmcsamp, par exemple (non exhaustif):

http://finzi.psych.upenn.edu/R/Rhelp02a/archive/81159.html
http://finzi.psych.upenn.edu/R/Rhelp02a/archive/60072.html
http://finzi.psych.upenn.edu/R/Rhelp02a/archive/76742.html
http://finzi.psych.upenn.edu/R/Rhelp02a/archive/82046.html

* une chaîne de longueur 1000 est très (trop !) courte. Il y a différentes fonctions de diagnostic dans coda qui permettent d'estimer la longueur minimal pour avoir une précision donnée (pas le temps de rentrer dans les détails). J'ai un bouquin sur MCMC (Gilks, Richardson et Spiegelhalter): passe le prendre à l'occasion.

* la fct HPDinterval de coda (écrite par Bates) simplifie le boulot de calcul des IC.

Amicalement,

Renaud


Retourner vers « Archives : Fonctions statistiques »

Qui est en ligne

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

cron