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