Intervalle de prédiction pour le modèle à effets mixtes lmer () dans R

37

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

http://glmm.wikidot.com/faq

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.merModdocumentation, 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 bootMercette tâche", mais, bon sang, je ne vois pas comment utiliser bootMerpour 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 secmoyen, 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.0dans 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 bootMerutilisations use.u=FALSEpar 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 armpaquet.

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 latticepaquet:

FishID Estimations

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()

Intervalle de prévision pour le poisson

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?

possible
la source
1
predictIntervalinclut l'erreur / incertitude pour les termes à effets fixes et aléatoires. En dotplotvous 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 fixe fishWtet 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 les predictIntervalrésultats.
jknowles

Réponses:

19

Cette question et cet excellent échange ont été le moteur de la création de la predictIntervalfonction dans le merToolspackage. bootMerC’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 predictIntervalest conçu pour utiliser les arm::simfonctions 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 valeur newdatafournie par l’utilisateur. C'est simple à utiliser - tout ce dont vous avez besoin est:

library(merTools)
preds <- predictInterval(lme1, newdata = newDat, n.sims = 999)

Vous pouvez spécifier toute une série d'autres valeurs à predictIntervalinclure, 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 thetaparamètres de l' lmerobjet n'est pas incluse, mais toutes les autres variations sont capturées par cette méthode, ce qui donne une approximation assez décente.

jknowles
la source
3
Cela a l'air génial! En parcourant la vignette maintenant. Merci!
possible
Les intervalles de prédiction ne sont pas tout à fait ce à quoi je m'attendais. Voir la mise à jour 3 ci-dessus.
possible
N'aime predictInterval()pas les effets aléatoires imbriqués? Par exemple, en utilisant l' msleepensemble de données du ggplot2package: 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
possible
Je parie que ça n'aime pas les effets imbriqués. Je ne pense pas que nous ayons eu des tests dans notre suite de tests. Je vais déposer un problème sur GitHub pour examiner cela. Je vous recommande également d'essayer d'abord la version dev de GitHub devtools::install_github("jknowles/merTools").
jknowles
2
En tant que mise à jour, la dernière version de développement de merTools autorise les effets imbriqués. Il sera poussé à CRAN sous peu.
jknowles
15

Faites ceci en faisant bootMergénérer un ensemble de prédictions pour chaque réplique d’amorçage paramétrique:

predFun <- function(fit) {
    predict(fit,newDat)
}
bb <- bootMer(lme1,nsim=200,FUN=predFun,seed=101)

La sortie de bootMerest dans un objet pas très transparent "boot", mais nous pouvons obtenir les prédictions brutes du $tcomposant.

Combien de temps Fish E a-t-il battu Fish D?

predMat <- bb$t
dim(predMat) ## 200 rows (PB reps) x 10 (predictions)

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:

mean(predMat[,5]<predMat[,4])  ## 0.57
Ben Bolker
la source
Je reçois des résultats inattendus. Si je fixe fishWt = 1 pour tous les poissons dans newDat, je m'attendrais à ce que le temps moyen / médian pour le poisson A soit ~ 10,1, le poisson B ~ 10.2, ..., le poisson J ~ 11,0 (car leur temps dans les données de formation est défini comme:) sec <- 9 + rep(1:10, rep(100,10))/10 + oz + rnorm(1000)/10. Lorsque j'utilise predict(), 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?
hossibley
Utiliser use.u=TRUEcomme dans: bb <- bootMer(lme1,nsim=200,FUN=predFun,seed=101,use.u=TRUE)semble me donner ce que je veux. Merci!
possible
OK, ça devient un peu délicat. Vous devez regarder l' use.uargument à 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 utiliser use.u=TRUE, mais je ne pense pas que cela fera nécessairement ce que vous voulez ...
Ben Bolker
Si j'utilise 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. bootMersuppose qu'il n'y a pas d'erreur dans notre estimation ponctuelle d'ER. Si je l’utilise use.u=FALSE, bootMerprend-il en compte les estimations de points RE? Il semble que les bootMerrésultats d'utilisation use.u=FALSEsoient équivalents (ou asymptotiquement équivalents) à ceux utilisés re.form=NAdans l' predict()instruction. Est-ce vrai?
possible
1
Je pense que ce n'est pas implémenté dans ATM, mais vous pouvez extraire les variances conditionnelles des modes conditionnels / BLUP via c(attr(ranef(lme1,condVar=TRUE)[[1]],"postVar"))(elles sont toutes identiques dans cet exemple), puis échantillonner ces valeurs.
Ben Bolker