Je veux obtenir un intervalle de prédiction autour d'une prédiction à partir d'un modèle lmer (). J'ai trouvé des discussions à ce sujet:
http://rstudio-pubs-static.s3.amazonaws.com/24365_2803ab8299934e888a60e7b16113f619.html
mais ils semblent ne pas tenir compte de l'incertitude des effets aléatoires.
Voici un exemple spécifique. Je cours des poissons d'or. J'ai des données sur les 100 dernières courses. Je veux prédire la 101ème, en tenant compte de l'incertitude de mes estimations RE et des estimations FE. J'inclus une interception aléatoire pour les poissons (il y a 10 poissons différents) et un effet fixe pour le poids (les poissons moins lourds sont plus rapides).
library("lme4")
fish <- as.factor(rep(letters[1:10], each=100))
race <- as.factor(rep(900:999, 10))
oz <- round(1 + rnorm(1000)/10, 3)
sec <- 9 + rep(1:10, rep(100,10))/10 + oz + rnorm(1000)/10
fishDat <- data.frame(fishID = fish,
raceID = race, fishWt = oz, time = sec)
head(fishDat)
plot(fishDat$fishID, fishDat$time)
lme1 <- lmer(time ~ fishWt + (1 | fishID), data=fishDat)
summary(lme1)
Maintenant, pour prédire la 101e course. Les poissons ont été pesés et sont prêts à partir:
newDat <- data.frame(fishID = letters[1:10],
raceID = rep(1000, 10),
fishWt = 1 + round(rnorm(10)/10, 3))
newDat$pred <- predict(lme1, newDat)
newDat
fishID raceID fishWt pred
1 a 1000 1.073 10.15348
2 b 1000 1.001 10.20107
3 c 1000 0.945 10.25978
4 d 1000 1.110 10.51753
5 e 1000 0.910 10.41511
6 f 1000 0.848 10.44547
7 g 1000 0.991 10.68678
8 h 1000 0.737 10.56929
9 i 1000 0.993 10.89564
10 j 1000 0.649 10.65480
Fish D s’est vraiment laissé aller (1,11 oz) et il est en fait prévu de perdre face à Fish E et Fish F, qui sont tous deux meilleurs que par le passé. Cependant, je veux maintenant pouvoir dire: "Le poisson E (pesant 0,91 oz) battra le poisson D (pesant 1,11 oz) avec une probabilité p." Existe-t-il un moyen de faire une telle déclaration en utilisant lme4? Je veux que ma probabilité p prenne en compte mon incertitude à la fois pour l'effet fixe et pour l'effet aléatoire.
Merci!
En regardant la predict.merMod
documentation, PS suggère que "Il n’ya pas d’option pour calculer les erreurs types des prévisions car il est difficile de définir une méthode efficace qui intègre l’incertitude dans les paramètres de variance; nous recommandons bootMer
cette tâche", mais, bon sang, je ne vois pas comment utiliser bootMer
pour faire cela. Cela semble bootMer
être utilisé pour obtenir des intervalles de confiance initialisés pour les estimations de paramètres, mais je peux me tromper.
MISE À JOUR Q:
OK, je pense que je posais la mauvaise question. Je veux pouvoir dire: "Le poisson A, pesant 1 kg, aura un temps de course de (lcl, ucl) 90% du temps."
Dans l'exemple que j'ai présenté, le poisson A, pesant 1,0 oz, aura un temps de course 9 + 0.1 + 1 = 10.1 sec
moyen, avec un écart type de 0,1. Ainsi, son temps de course observé sera entre
x <- rnorm(mean = 10.1, sd = 0.1, n=10000)
quantile(x, c(0.05,0.50,0.95))
5% 50% 95%
9.938541 10.100032 10.261243
90% du temps. Je veux une fonction de prédiction qui tente de me donner cette réponse. Régler tous fishWt = 1.0
dans newDat
, réexécuter la carte SIM, et en utilisant (comme suggéré par Ben Bolker ci - dessous)
predFun <- function(fit) {
predict(fit,newDat)
}
bb <- bootMer(lme1,nsim=1000,FUN=predFun, use.u = FALSE)
predMat <- bb$t
donne
> quantile(predMat[,1], c(0.05,0.50,0.95))
5% 50% 95%
10.01362 10.55646 11.05462
Cela semble être réellement centré sur la moyenne de la population? Comme si cela ne tenait pas compte de l'effet FishID? Je pensais que c'était peut-être un problème de taille d'échantillon, mais quand j'ai augmenté le nombre de courses observées de 100 à 10 000, j'ai toujours des résultats similaires.
Je noterai les bootMer
utilisations use.u=FALSE
par défaut. D'un autre côté, en utilisant
bb <- bootMer(lme1,nsim=1000,FUN=predFun, use.u = TRUE)
donne
> quantile(predMat[,1], c(0.05,0.50,0.95))
5% 50% 95%
10.09970 10.10128 10.10270
Cet intervalle est trop étroit et semblerait être un intervalle de confiance pour le temps moyen du poisson A. Je veux un intervalle de confiance pour le temps de course observé de Fish A, et non pour son temps de course moyen. Comment puis-je l'obtenir?
MISE À JOUR 2, PRESQUE:
Je pensais avoir trouvé ce que je cherchais dans Gelman and Hill (2007) , page 273. Besoin d'utiliser le arm
paquet.
library("arm")
Pour le poisson A:
x.tilde <- 1 #observed fishWt for new race
sigma.y.hat <- sigma.hat(lme1)$sigma$data #get uncertainty estimate of our model
coef.hat <- as.matrix(coef(lme1)$fishID)[1,] #get intercept (random) and fishWt (fixed) parameter estimates
y.tilde <- rnorm(1000, coef.hat %*% c(1, x.tilde), sigma.y.hat) #simulate
quantile (y.tilde, c(.05, .5, .95))
5% 50% 95%
9.930695 10.100209 10.263551
Pour tous les poissons:
x.tilde <- rep(1,10) #assume all fish weight 1 oz
#x.tilde <- 1 + rnorm(10)/10 #alternatively, draw random weights as in original example
sigma.y.hat <- sigma.hat(lme1)$sigma$data
coef.hat <- as.matrix(coef(lme1)$fishID)
y.tilde <- matrix(rnorm(1000, coef.hat %*% matrix(c(rep(1,10), x.tilde), nrow = 2 , byrow = TRUE), sigma.y.hat), ncol = 10, byrow = TRUE)
quantile (y.tilde[,1], c(.05, .5, .95))
5% 50% 95%
9.937138 10.102627 10.234616
En fait, ce n'est probablement pas exactement ce que je veux. Je ne prends en compte que l'incertitude globale du modèle. Dans une situation où, par exemple, j'ai observé 5 races observées pour le poisson K et 1 000 races observées pour le poisson L, l'incertitude associée à ma prédiction pour Fish K devrait être beaucoup plus grande que celle associée à ma prédiction pour Fish L.
Nous examinerons plus en détail Gelman et Hill 2007. Je pense que je pourrais finir par devoir passer à BUGS (ou à Stan).
METTRE À JOUR LE 3:
Peut-être que je conceptualise mal les choses. Utiliser la predictInterval()
fonction donnée par Jared Knowles dans une réponse ci-dessous donne des intervalles qui ne sont pas tout à fait ce à quoi je m'attendais ...
library("lattice")
library("lme4")
library("ggplot2")
fish <- c(rep(letters[1:10], each = 100), rep("k", 995), rep("l", 5))
oz <- round(1 + rnorm(2000)/10, 3)
sec <- 9 + c(rep(1:10, each = 100)/10,rep(1.1, 995), rep(1.2, 5)) + oz + rnorm(2000)
fishDat <- data.frame(fishID = fish, fishWt = oz, time = sec)
dim(fishDat)
head(fishDat)
plot(fishDat$fishID, fishDat$time)
lme1 <- lmer(time ~ fishWt + (1 | fishID), data=fishDat)
summary(lme1)
dotplot(ranef(lme1, condVar = TRUE))
J'ai ajouté deux nouveaux poissons. Fish K, pour qui nous avons observé 995 races, et Fish L, pour qui nous avons observé 5 races. Nous avons observé 100 courses pour Fish AJ. Je correspond le même lmer()
qu'avant. En regardant dotplot()
le lattice
paquet:
Par défaut, dotplot()
réordonne les effets aléatoires en fonction de leur estimation ponctuelle. L'estimation pour Fish L est sur la ligne supérieure et présente un intervalle de confiance très large. Le poisson K est sur la troisième ligne et a un intervalle de confiance très étroit. Cela a du sens pour moi. Nous avons beaucoup de données sur Fish K, mais pas beaucoup de données sur Fish L, nous sommes donc plus confiants dans nos prévisions concernant la vitesse de nage réelle de Fish K. Maintenant, je pense que cela conduirait à un intervalle de prédiction étroit pour Fish K et à un intervalle de prédiction large pour Fish L lors de l'utilisation predictInterval()
. Howeva:
newDat <- data.frame(fishID = letters[1:12],
fishWt = 1)
preds <- predictInterval(lme1, newdata = newDat, n.sims = 999)
preds
ggplot(aes(x=letters[1:12], y=fit, ymin=lwr, ymax=upr), data=preds) +
geom_point() +
geom_linerange() +
labs(x="Index", y="Prediction w/ 95% PI") + theme_bw()
Tous ces intervalles de prédiction semblent être identiques en largeur. Pourquoi notre prédiction pour Fish K n'est-elle pas plus étroite que les autres? Pourquoi notre prédiction pour Fish L n'est-elle pas plus large que les autres?
predictInterval
inclut l'erreur / incertitude pour les termes à effets fixes et aléatoires. Endotplot
vous ne voyez l'incertitude due à la partie aléatoire de la prédiction, essentiellement l'incertitude autour de l'estimation des intersections spécifiques poissons. Si votre modèle a beaucoup d'incertitude dans le paramètre fixefishWt
et que ce paramètre détermine la majeure partie de la valeur prévue, l'incertitude entourant une interception de poisson spécifique est triviale et vous ne verrez pas une grande différence dans la largeur des intervalles. Nous devrions rendre cela plus clair dans lespredictInterval
résultats.Réponses:
Cette question et cet excellent échange ont été le moteur de la création de la
predictInterval
fonction dans lemerTools
package.bootMer
C’est la voie à suivre, mais pour certains problèmes, il n’est pas faisable sur le plan informatique de générer des remises en état amputées du modèle entier (dans les cas où le modèle est volumineux).Dans ces cas, il
predictInterval
est conçu pour utiliser lesarm::sim
fonctions afin de générer des distributions de paramètres dans le modèle, puis d’utiliser ces distributions pour générer des valeurs simulées de la réponse compte tenu de la valeurnewdata
fournie par l’utilisateur. C'est simple à utiliser - tout ce dont vous avez besoin est:Vous pouvez spécifier toute une série d'autres valeurs à
predictInterval
inclure, notamment en définissant l'intervalle pour les intervalles de prédiction, en choisissant si vous souhaitez indiquer la moyenne ou la médiane de la distribution et en choisissant d'inclure ou non la variance résiduelle du modèle.Ce n'est pas un intervalle de prédiction complet car la variabilité des
theta
paramètres de l'lmer
objet n'est pas incluse, mais toutes les autres variations sont capturées par cette méthode, ce qui donne une approximation assez décente.la source
predictInterval()
pas les effets aléatoires imbriqués? Par exemple, en utilisant l'msleep
ensemble de données duggplot2
package:mod <- lmer(sleep_total ~ bodywt + (1|vore/order), data=msleep); predInt <- predictInterval(merMod=mod, newdata=msleep)
Renvoie une erreur:Error in '[.data.frame'(newdata, , j) : undefined columns selected
devtools::install_github("jknowles/merTools")
.Faites ceci en faisant
bootMer
générer un ensemble de prédictions pour chaque réplique d’amorçage paramétrique:La sortie de
bootMer
est dans un objet pas très transparent"boot"
, mais nous pouvons obtenir les prédictions brutes du$t
composant.Combien de temps Fish E a-t-il battu Fish D?
Les temps du poisson E sont dans la colonne 5, ceux du poisson D dans la colonne 4. Il suffit donc de savoir dans quelle proportion la colonne 5 est inférieure à la colonne 4:
la source
sec <- 9 + rep(1:10, rep(100,10))/10 + oz + rnorm(1000)/10
. Lorsque j'utilisepredict()
, les temps de prédiction pour les poissons A, E et J sont 10.09, 10.49 et 10.99, comme prévu. Cependant, les temps médians pour la méthode bootMer que vous avez décrite sont: 10.52, 10.59 et 10.50. J'aurais attendu plus d'accord?use.u=TRUE
comme dans:bb <- bootMer(lme1,nsim=200,FUN=predFun,seed=101,use.u=TRUE)
semble me donner ce que je veux. Merci!use.u
argument àbootMer
. La question qui se pose est la suivante: lorsque vous parlez d'incertitude quant à l'effet fixe et à l'effet aléatoire, qu'entendez-vous par «effet aléatoire»? Voulez-vous dire incertitude dans la variance à effets aléatoires ou dans les modes conditionnels (c.-à-d. Les effets spécifiques au poisson)? Vous pouvez utiliseruse.u=TRUE
, mais je ne pense pas que cela fera nécessairement ce que vous voulez ...use.u=TRUE
, alors les "valeurs de u [rester] fixées à leurs valeurs estimées". J'interprète cela comme une signification, quelle que soit notre estimation ponctuelle d'effet aléatoire pour le poisson A, elle est considérée comme la vérité honnête de Dieu, si vous voulez.bootMer
suppose qu'il n'y a pas d'erreur dans notre estimation ponctuelle d'ER. Si je l’utiliseuse.u=FALSE
,bootMer
prend-il en compte les estimations de points RE? Il semble que lesbootMer
résultats d'utilisationuse.u=FALSE
soient équivalents (ou asymptotiquement équivalents) à ceux utilisésre.form=NA
dans l'predict()
instruction. Est-ce vrai?c(attr(ranef(lme1,condVar=TRUE)[[1]],"postVar"))
(elles sont toutes identiques dans cet exemple), puis échantillonner ces valeurs.