Nieliniowa regresja efektów mieszanych w R

14

Co zaskakujące, nie byłem w stanie znaleźć odpowiedzi na następujące pytanie za pomocą Google:

Mam pewne dane biologiczne od kilku osób, które pokazują z grubsza esicy zachowanie wzrostu w czasie. Dlatego chcę go modelować przy użyciu standardowego wzrostu logistycznego

P(t) = k*p0*exp(r*t) / (k+p0*(exp(r*t)-1))

gdzie p0 jest wartością początkową przy t = 0, k oznacza asymptotyczną granicę przy t-> nieskończoności, zaś r jest prędkością wzrostu. O ile widzę, mogę łatwo modelować to za pomocą nls (brak zrozumienia z mojej strony: dlaczego nie mogę modelować czegoś podobnego za pomocą standardowej regresji logit poprzez skalowanie czasu i danych? EDYCJA: Dzięki Nick, najwyraźniej ludzie robią to np. Dla proporcje, ale rzadko http://www.stata-journal.com/article.html?article=st0147 . Następnym pytaniem dotyczącym tej stycznej byłoby pytanie, czy model jest w stanie obsłużyć wartości odstające> 1).

Teraz chcę pozwolić na pewne stałe (głównie kategoryczne) i niektóre losowe (indywidualny identyfikator i ewentualnie także identyfikator badania) wpływ na trzy parametry k, p0 i r. Czy nlme jest najlepszym sposobem na zrobienie tego? Model SSlogis wydaje się rozsądny z punktu widzenia tego, co próbuję zrobić, czy to prawda? Czy któryś z poniższych jest rozsądnym modelem na początek? Nie mogę poprawnie wyregulować wartości początkowych, a update () wydaje się działać tylko w przypadku efektów losowych, a nie stałych - jakieś wskazówki?

nlme(y ~ k*p0*exp(r*t) / (k+p0*(exp(r*t)-1)), ## not working at all (bad numerical properties?)
            data = data,
            fixed = k + p0 + r ~ var1 + var2,
            random = k + p0 + r ~ 1|UID,
            start = c(p0=1, k=100, r=1))

nlme(y ~ SSlogis(t, Asym, xmid, scal), ## not working, as start= is inappropriate
            data = data,
            fixed = Asym + xmid + scal ~ var1 + var2, ## works fine with ~ 1
            random = Asym + xmid + scal ~ 1|UID,
            start = getInitial(y ~ SSlogis(Dauer, Asym, xmid, scal), data = data))

Ponieważ jestem nowy w nieliniowych modelach mieszanych, aw szczególności w modelach nieliniowych, doceniłbym niektóre zalecenia dotyczące czytania lub linki do samouczków / często zadawanych pytań z pytaniami dla początkujących.

Rob Hall
źródło
2
Jeśli uważasz, że k jest znane, możesz skalować populacje o P / k. Jeśli k jest czymś do oszacowania, to samo oznacza, że ​​twój problem nie jest standardową regresją logit.
Nick Cox,
1
Dziękuję Nick. Tak, ostatecznie uważam, że k należy oszacować i uwzględnić w regresji. Moje zainteresowanie wykorzystaniem regresji logitów było czysto akademickie. Pomyślałem, że może to być dobry model na początek przed przejściem do modelowania nieliniowego, ale nie byłem w stanie znaleźć żadnych przykładów regresji logit dla danych nie-binarnych za pomocą Google. Zastanawiałem się, czy istnieje jakiś powód (np. Założenia dystrybucyjne dotyczące błędów), które sprawiają, że używanie np. Brokera z łączem logit z ciągłymi danymi jest złym pomysłem.
Rob Hall
3
Modelowanie logitów dla odpowiedzi, które mają ciągłe proporcje, istnieje już od pewnego czasu, ale najwyraźniej nie jest dobrze znane. Zobacz np. Baum w stata-journal.com/sjpdf.html?articlenum=st0147 Niemniej jednak to nie twoja sytuacja. Nie mogę komentować implementacji R.
Nick Cox,
Dziękuję Nickowi za ten interesujący link - który wyjaśnia mi kilka rzeczy. Niestety wydaje się, że nie mogę jeszcze głosować za twoją odpowiedzią. (W przypadku, gdy ludzie mają problemy z dostępem do bezpośredniego linku, działało dla mnie: stata-journal.com/article.html?article=st0147 )
Rob Hall
1
Wzrost logistyczny implikuje monotoniczną krzywą wzrostu. Jeśli dane się nie zgadzają, źle się dopasujesz lub oprogramowanie nie będzie odtwarzane, w zależności od tego, co robisz.
Nick Cox,

Odpowiedzi:

12

Chciałem podzielić się niektórymi rzeczami, których nauczyłem się od czasu zadawania tego pytania. nlme wydaje się rozsądnym sposobem modelowania nieliniowych efektów mieszanych w R. Zacznij od prostego modelu podstawowego:

library(nlme)
data <- groupedData(y ~ t | UID, data=data) ## not strictly necessary
initVals <- getInitial(y ~ SSlogis(t, Asym, xmid, scal), data = data)
baseModel<- nlme(y ~ SSlogis(t, Asym, xmid, scal),
    data = data,
    fixed = list(Asym ~ 1, xmid ~ 1, scal ~ 1),
    random = Asym + xmid + scal ~ 1|UID,
    start = initVals
)

Następnie użyj aktualizacji, aby zwiększyć złożoność modelu. Parametr początkowy jest nieco trudny w obsłudze, może zająć trochę majsterkowania, aby ustalić kolejność. Zwróć uwagę, że nowy stały efekt dla var1 w Asymie jest zgodny ze zwykłym stałym efektem dla Asym.

 nestedModel <- update(baseModel, fixed=list(Asym ~ var1, xmid ~ 1, scal ~ 1), start = c(fixef(baseModel)[1], 0, fixef(baseModel)[2], fixef(baseModel)[3]))

lme4 wydawał się bardziej odporny na wartości odstające w moim zbiorze danych i wydawał się oferować bardziej niezawodną zbieżność dla bardziej złożonych modeli. Jednak wydaje się, że wadą jest to, że odpowiednie funkcje prawdopodobieństwa muszą być określone ręcznie. Poniżej znajduje się logistyczny model wzrostu ze stałym wpływem var1 (binarnie) na Asym. Możesz dodawać stałe efekty do Xmid i skalować w podobny sposób. Zwróć uwagę na dziwny sposób określania modelu przy użyciu podwójnej formuły jako wyniku ~ efekty ustalone ~ efekty losowe.

library(lme4) ## careful loading nlme and lme4 concurrently
customLogitModel <- function(t, Asym, AsymVar1, xmid, scal) {
    (Asym+AsymVar1*var1)/(1+exp((xmid-t)/scal))
}

customLogitModelGradient <- deriv(
    body(customLogitModel)[[2]], 
    namevec = c("Asym", "AsymVar1", "xmid", "scal"), 
    function.arg=customLogitModel
)

## find starting parameters
initVals <- getInitial(y ~ SSlogis(t, Asym, xmid, scal), data = data)

# Fit the model
model <- nlmer(
    y ~ customLogitModelGradient(t=t, Asym, AsymVar1, xmid, scal, var1=var) ~ 
    # Random effects with a second ~
    (Asym | UID) + (xmid | UID) + (scal | UID), 
    data = data, 
    start = c(Asym=initVals[1], AsymVar1=0, xmid=initVals[2], scal=initVals[3])
)
Rob Hall
źródło
1
Dziękuję Rob za twój post, to jest dokładnie to, co próbuję zrobić z moimi danymi. Nie rozumiem, co to jest var1 w nestedModel w Asym i jak go obliczyłeś?
To tylko przykład, jak uwzględnić wpływ jakiejś zmiennej na Asym: „Poniżej znajduje się logistyczny model wzrostu z ustalonym efektem var1 (binarnie) na Asym”. Na przykład masz zmienną „Traktowane”, która ma dwie wartości 0 i 1, więc zamiast „var1” zamień „Leczone”.
PA6OTA