Chciałbym uzyskać 95% przedziały ufności dla prognoz nieliniowego nlme
modelu mieszanego . Ponieważ nie ma w tym celu standardu nlme
, zastanawiałem się, czy słuszne jest zastosowanie metody „przedziałów prognozowania populacji”, jak opisano w rozdziale książki Bena Bolkera w kontekście modeli pasujących z najwyższym prawdopodobieństwem , w oparciu o ideę ponownie próbkować parametry ustalonego efektu w oparciu o macierz wariancji-kowariancji dopasowanego modelu, symulując oparte na tym prognozy, a następnie biorąc 95% percentyle tych prognoz, aby uzyskać 95% przedziały ufności?
Kod do wykonania tego wygląda następująco: (Używam tutaj danych „Loblolly” z nlme
pliku pomocy)
library(effects)
library(nlme)
library(MASS)
fm1 <- nlme(height ~ SSasymp(age, Asym, R0, lrc),
data = Loblolly,
fixed = Asym + R0 + lrc ~ 1,
random = Asym ~ 1,
start = c(Asym = 103, R0 = -8.5, lrc = -3.3))
xvals=seq(min(Loblolly$age),max(Loblolly$age),length.out=100)
nresamp=1000
pars.picked = mvrnorm(nresamp, mu = fixef(fm1), Sigma = vcov(fm1)) # pick new parameter values by sampling from multivariate normal distribution based on fit
yvals = matrix(0, nrow = nresamp, ncol = length(xvals))
for (i in 1:nresamp)
{
yvals[i,] = sapply(xvals,function (x) SSasymp(x,pars.picked[i,1], pars.picked[i,2], pars.picked[i,3]))
}
quant = function(col) quantile(col, c(0.025,0.975)) # 95% percentiles
conflims = apply(yvals,2,quant) # 95% confidence intervals
Teraz, gdy mam już granice pewności, tworzę wykres:
meany = sapply(xvals,function (x) SSasymp(x,fixef(fm1)[[1]], fixef(fm1)[[2]], fixef(fm1)[[3]]))
par(cex.axis = 2.0, cex.lab=2.0)
plot(0, type='n', xlim=c(3,25), ylim=c(0,65), axes=F, xlab="age", ylab="height");
axis(1, at=c(3,1:5 * 5), labels=c(3,1:5 * 5))
axis(2, at=0:6 * 10, labels=0:6 * 10)
for(i in 1:14)
{
data = subset(Loblolly, Loblolly$Seed == unique(Loblolly$Seed)[i])
lines(data$age, data$height, col = "red", lty=3)
}
lines(xvals,meany, lwd=3)
lines(xvals,conflims[1,])
lines(xvals,conflims[2,])
Oto wykres z 95% przedziałami ufności uzyskanymi w ten sposób:
Czy to podejście jest poprawne, czy istnieją jakieś inne lub lepsze podejścia do obliczania 95% przedziałów ufności dla prognoz nieliniowego modelu mieszanego? Nie jestem całkowicie pewien, jak poradzić sobie ze strukturą efektu losowego modelu ... Czy powinna być średnia ponad poziomy efektów losowych? Czy też byłoby OK mieć przedziały ufności dla przeciętnego przedmiotu, które wydawałyby się bliższe temu, co mam teraz?
źródło
Odpowiedzi:
To, co tu zrobiłeś, wygląda rozsądnie. Krótka odpowiedź jest taka, że w przeważającej części kwestie przewidywania przedziałów ufności z modeli mieszanych i modeli nieliniowych są mniej więcej ortogonalne , to znaczy, że musisz się martwić o oba zestawy problemów, ale tak nie jest (wiem o tym ) wchodzą w interakcje w jakikolwiek dziwny sposób.
W całym założeniu zakładam, że przewidujesz na poziomie populacji i konstruujesz przedziały ufności jako poziom populacji - innymi słowy, próbujesz wykreślić przewidywane wartości typowej grupy, a nie uwzględniać zmienności między grupami w zakresie Twojej pewności interwały. Upraszcza to problemy z modelem mieszanym. Poniższe wykresy porównują trzy podejścia (zrzut kodu poniżej):
Możemy również zrobić ładowanie parametryczne ...
Oto wykresy CI wraz z danymi ...
... ale trudno nam dostrzec różnice.
Powiększanie odejmując przewidywane wartości (czerwony = bootstrap, niebieski = PPI, cyan = delta)
W tym przypadku interwały ładowania są w rzeczywistości najwęższe (np. Przypuszczalnie rozkłady próbkowania parametrów są w rzeczywistości nieco cieńsze niż normalne), podczas gdy interwały PPI i metody delta są bardzo do siebie podobne.
źródło