Nieliniowa regresja efektów mieszanych w R

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.

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