Я хочу получить интервал прогнозирования вокруг прогноза из модели lmer (). Я нашел некоторое обсуждение по этому поводу:
http://rstudio-pubs-static.s3.amazonaws.com/24365_2803ab8299934e888a60e7b16113f619.html
но они, похоже, не учитывают неопределенность случайных эффектов.
Вот конкретный пример. Я гоняю золотую рыбку. У меня есть данные о последних 100 гонках. Я хочу предсказать 101-й, принимая во внимание неопределенность моих оценок RE и оценок FE. Я включил случайный перехват для рыбы (есть 10 различных рыб) и фиксированный эффект для веса (менее тяжелая рыба быстрее).
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)
Теперь прогнозируем 101-ю гонку. Рыба была взвешена и готова к отправке:
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
Рыба D действительно отпустила себя (1,11 унции) и, по прогнозам, проиграет Рыбе E и Рыбе F, с которыми он был лучше, чем в прошлом. Однако теперь я хочу сказать: «Рыба E (весом 0,91 унции) с вероятностью p победит Рыбу D (весом 1,11 унции)». Есть ли способ сделать такое заявление, используя lme4? Я хочу, чтобы моя вероятность p учитывала мою неопределенность как фиксированного эффекта, так и случайного эффекта.
Благодарность!
PS, глядя на predict.merMod
документацию, он предлагает «Нет никакой возможности для вычисления стандартных ошибок прогнозов, потому что трудно определить эффективный метод, который включает неопределенность в параметрах отклонения; мы рекомендуем bootMer
для этой задачи», но, черт возьми, я не вижу как использовать, bootMer
чтобы сделать это. Похоже bootMer
, будет использоваться для получения доверительных интервалов для оценки параметров, но я могу ошибаться.
ОБНОВЛЕНО Q:
ОК, я думаю, что задавал неправильный вопрос. Я хочу иметь возможность сказать: «У рыбы А весом в унцию будет время гонки, которое (lcl, ucl) составляет 90% времени».
В примере, который я выложил, у рыбы А весом 1,0 унции будет время гонки в 9 + 0.1 + 1 = 10.1 sec
среднем со стандартным отклонением 0,1. Таким образом, его наблюдаемое время гонки будет между
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% времени Я хочу функцию прогнозирования, которая пытается дать мне этот ответ. Установка олл- fishWt = 1.0
ин newDat
, повторный запуск сима и использование (как предложено Беном Болкером ниже)
predFun <- function(fit) {
predict(fit,newDat)
}
bb <- bootMer(lme1,nsim=1000,FUN=predFun, use.u = FALSE)
predMat <- bb$t
дает
> quantile(predMat[,1], c(0.05,0.50,0.95))
5% 50% 95%
10.01362 10.55646 11.05462
Это, кажется, на самом деле сосредоточено вокруг среднего населения? Как будто он не учитывает эффект FishID? Я подумал, что, возможно, это проблема размера выборки, но когда я увеличил количество наблюдаемых гонок со 100 до 10000, я все равно получил похожие результаты.
Я отмечу bootMer
использование use.u=FALSE
по умолчанию. С другой стороны, используя
bb <- bootMer(lme1,nsim=1000,FUN=predFun, use.u = TRUE)
дает
> quantile(predMat[,1], c(0.05,0.50,0.95))
5% 50% 95%
10.09970 10.10128 10.10270
Этот интервал слишком узок и может показаться доверительным интервалом для среднего времени рыбы А. Мне нужен доверительный интервал для наблюдаемого времени гонки Рыбы А, а не его среднего времени гонки. Как я могу получить это?
ОБНОВЛЕНИЕ 2, ПОЧТИ:
Я думал, что нашел то, что искал, в Gelman and Hill (2007) , стр. 273. Нужно использовать arm
пакет.
library("arm")
Для рыбы А:
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
Для всех рыб:
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
На самом деле, это, вероятно, не совсем то, что я хочу. Я только принимаю во внимание общую неопределенность модели. В ситуации, когда у меня, скажем, 5 наблюдаемых гонок для рыбы K и 1000 наблюдаемых гонок для рыбы L, я думаю, что неопределенность, связанная с моим прогнозом для рыбы K, должна быть намного больше, чем неопределенность, связанная с моим прогнозом для рыбы L.
Посмотрим дальше на Gelman and Hill 2007. Я чувствую, что в итоге мне может понадобиться переключиться на BUGS (или Stan).
ОБНОВЛЕНИЕ 3-го:
Возможно, я плохо концептуализирую вещи. Использование predictInterval()
функции, данной Джаредом Ноулзом в ответе ниже, дает интервалы, которые не совсем то, что я ожидаю ...
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))
Я добавил две новые рыбы. Fish K, для которого мы наблюдали 995 гонок, и Fish L, для которого мы наблюдали 5 гонок. Мы наблюдали 100 гонок за рыбой AJ. Я подхожу так же, lmer()
как и раньше. Глядя на dotplot()
из lattice
пакета:
По умолчанию dotplot()
случайные эффекты переупорядочиваются по их точечной оценке. Оценка для Рыбы L находится в верхней строке и имеет очень широкий доверительный интервал. Рыба K находится на третьей строчке и имеет очень узкий доверительный интервал. Это имеет смысл для меня. У нас есть много данных о Fish K, но не много данных о Fish L, поэтому мы более уверены в наших предположениях об истинной скорости плавания Fish K. Теперь, я думаю, это приведет к узкому интервалу предсказания для Рыбы K и широкому интервалу предсказания для Рыбы L при использовании 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()
Все эти интервалы прогнозирования кажутся одинаковыми по ширине. Почему наш прогноз для Fish K уже, чем другие? Почему наш прогноз для Fish L не шире, чем у других?
источник
predictInterval
включает в себя ошибку / неопределенность для фиксированных и случайных членов эффекта. Вdotplot
вас вы видите только неопределенность из-за случайной части прогноза, по существу неопределенность вокруг оценки конкретных перехватов рыбы. Если ваша модель имеет большую неопределенность в фиксированном параметре,fishWt
и этот параметр определяет большую часть прогнозируемого значения, тогда неопределенность вокруг любого конкретного перехвата рыбы является тривиальной, и вы не увидите большой разницы в ширине интервалов. Мы должны сделать это более ясным вpredictInterval
результатах.Ответы:
Этот вопрос и отличный обмен послужили толчком для создания
predictInterval
функции вmerTools
пакете.bootMer
это путь, но для некоторых проблем не представляется возможным в вычислительном отношении сгенерировать начальные модификации всей модели (в случаях, когда модель велика).В этих случаях
predictInterval
предназначен для использованияarm::sim
функций для генерации распределений параметров в модели, а затем для использования этих распределений для генерации смоделированных значений отклика с учетомnewdata
предоставленных пользователем. Он прост в использовании - все, что вам нужно сделать, это:Можно указать целый ряд других значений, в
predictInterval
том числе установить интервал для интервалов прогнозирования, выбрать, сообщать ли среднее значение или медиану распределения, и выбрать, включать ли остаточную дисперсию из модели.Это не полный интервал прогнозирования, потому что изменчивость
theta
параметров вlmer
объекте не учитывается, но все остальные вариации фиксируются с помощью этого метода, давая довольно приличное приближение.источник
predictInterval()
любит вложенные случайные эффекты? Например, используяmsleep
набор данных изggplot2
пакета:mod <- lmer(sleep_total ~ bodywt + (1|vore/order), data=msleep); predInt <- predictInterval(merMod=mod, newdata=msleep)
Возвращает ошибку:Error in '[.data.frame'(newdata, , j) : undefined columns selected
devtools::install_github("jknowles/merTools")
сначала попробовать версию dev из GitHub .Сделайте это,
bootMer
сгенерировав набор прогнозов для каждой параметрической копии начальной загрузки:Вывод
bootMer
находится в не очень прозрачном"boot"
объекте, но мы можем получить необработанные прогнозы из$t
компонента.Сколько времени Fish E побеждает Fish D?
Время рыбы Е в столбце 5, время рыбы D в столбце 4, поэтому нам просто нужно знать, какая доля в столбце 5 меньше столбца 4:
источник
sec <- 9 + rep(1:10, rep(100,10))/10 + oz + rnorm(1000)/10
. Когда я используюpredict()
, время прогнозирования для рыб A, E и J составляет 10,09, 10,49 и 10,99, как и ожидалось. Тем не менее, среднее время для метода bootMer, который вы описали, составляет: 10,52, 10,59 и 10,50. Я бы ожидал большего соглашения?use.u=TRUE
как в:bb <- bootMer(lme1,nsim=200,FUN=predFun,seed=101,use.u=TRUE)
кажется, дает мне то, что я хочу. Благодарность!use.u
аргументbootMer
. Вопрос в том, что когда вы говорите «неопределенность фиксированного эффекта и случайного эффекта», что вы подразумеваете под «случайным эффектом»? Вы имеете в виду неопределенность в дисперсии случайных эффектов или в условных режимах (то есть специфических для рыб эффектов)? Вы можете использоватьuse.u=TRUE
, но я не думаю, что он обязательно сделает то, что вы хотите ...use.u=TRUE
, то «значения u [остаются] фиксированными на их оценочные значения». Я интерпретирую это как значение, какой бы ни была наша оценка точки случайного воздействия для Рыбы А, она принимается за Честную Истину Бога, если хотите.bootMer
Предполагается, что в нашей оценке RE нет ошибок. Если я используюuse.u=FALSE
,bootMer
учитывает ли вообще точечные оценки RE? Кажется, чтоbootMer
результаты при использованииuse.u=FALSE
эквивалентны (или, асимптотически эквивалентны) использованиюre.form=NA
вpredict()
выражении. Это правда?c(attr(ranef(lme1,condVar=TRUE)[[1]],"postVar"))
(все они идентичны в этом примере), а затем сэмплировать эти значения.