Regression nichtlinearer Mischeffekte in R

14

Überraschenderweise konnte ich mit Google keine Antwort auf die folgende Frage finden:

Ich habe einige biologische Daten von mehreren Personen, die mit der Zeit ein grob sigmoides Wachstumsverhalten zeigen. Daher möchte ich es mit einem logistischen Standardwachstum modellieren

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

wobei p0 der Startwert bei t = 0 ist, k die asymptotische Grenze bei t → unendlich ist und r die Wachstumsgeschwindigkeit ist. Soweit ich sehen kann, kann ich dies leicht mit nls modellieren (mangelndes Verständnis meinerseits: Warum kann ich etwas Ähnliches nicht mit Standard-Logit-Regression durch Skalieren von Zeit und Daten modellieren? Proportionen, aber selten http://www.stata-journal.com/article.html?article=st0147 Die nächste Frage zu dieser Tangente wäre, ob das Modell möglicherweise Ausreißer> 1) handhaben kann.

Nun möchte ich einige feste (hauptsächlich kategoriale) und einige zufällige (eine individuelle ID und möglicherweise auch eine Studien-ID) Effekte auf die drei Parameter k, p0 und r zulassen. Ist nlme der beste Weg, dies zu tun? Das SSlogis-Modell scheint für das, was ich versuche, sinnvoll zu sein. Ist das richtig? Ist eines der folgenden Modelle von Anfang an sinnvoll? Ich kann die Startwerte scheinbar nicht richtig ermitteln und update () scheint nur für zufällige Effekte zu funktionieren, nicht für festgelegte - irgendwelche Hinweise?

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

Da ich neu in nichtlinearen gemischten Modellen im Besonderen und nichtlinearen Modellen im Allgemeinen bin, würde ich mich über einige Leseempfehlungen oder Links zu Tutorials / FAQs mit neuen Fragen freuen.

Rob Hall
quelle
2
Wenn Sie k als bekannt ansehen, können Sie Ihre Populationen nach P / k skalieren. Wenn k geschätzt werden muss, bedeutet dies allein, dass Ihr Problem keine Standard-Logit-Regression ist.
Nick Cox
1
Danke Nick. Ja, am Ende glaube ich, dass k geschätzt und in die Regression einbezogen werden muss. Mein Interesse an der Verwendung der logit-Regression war rein akademisch. Ich dachte, dies wäre vielleicht ein gutes Modell, bevor ich zur nichtlinearen Modellierung überging, aber ich konnte mit Google keine Beispiele für die Logit-Regression für nichtbinäre Daten finden. Ich habe mich gefragt, ob es einen Grund gibt (z. B. Verteilungsannahmen zu den Fehlern), der es zu einer schlechten Idee macht, z. B. glmer mit einem Logit-Link mit kontinuierlichen Daten zu verwenden.
Rob Hall
3
Logit-Modellierung für Antworten mit stetigen Proportionen gibt es schon seit einiger Zeit, aber anscheinend ist sie nicht gut bekannt. Siehe zB Baum in stata-journal.com/sjpdf.html?articlenum=st0147 Trotzdem ist es nicht Ihre Situation. Ich kann R-Implementierungen nicht kommentieren.
Nick Cox
Danke Nick für diesen interessanten Link - das klärt ein paar Dinge für mich. Leider kann ich Ihrer Antwort noch nicht zustimmen. (Falls die Leute Probleme haben, auf den direkten Link zuzugreifen, funktionierte für mich Folgendes: stata-journal.com/article.html?article=st0147 )
Rob Hall
1
Logistisches Wachstum impliziert eine monoton ansteigende Kurve. Wenn die Daten nicht übereinstimmen, ist die Übereinstimmung schlecht, oder die Software wird nicht abgespielt, je nachdem, was Sie gerade tun.
Nick Cox

Antworten:

12

Ich wollte einige der Dinge teilen, die ich gelernt habe, seit ich diese Frage gestellt habe. nlme scheint ein vernünftiger Weg zu sein, um nichtlineare Mischeffekte in R zu modellieren. Beginnen Sie mit einem einfachen Basismodell:

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
)

Verwenden Sie dann update, um die Modellkomplexität zu erhöhen. Der Start-Parameter ist etwas schwierig zu handhaben, es kann einige Bastelarbeiten erfordern, um die Reihenfolge herauszufinden. Beachten Sie, wie der neue feste Effekt für var1 in Asym dem regulären festen Effekt für Asym folgt.

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

lme4 schien gegenüber den Ausreißern in meinem Datensatz robuster zu sein und eine zuverlässigere Konvergenz für die komplexeren Modelle zu bieten. Es scheint jedoch der Nachteil zu sein, dass die relevanten Wahrscheinlichkeitsfunktionen manuell angegeben werden müssen. Das folgende ist das logistische Wachstumsmodell mit einer festen Auswirkung von var1 (binär) auf Asym. Sie können auf ähnliche Weise feste Effekte für xmid und scal hinzufügen. Beachten Sie die seltsame Art und Weise, das Modell mithilfe einer Doppelformel als Ergebnis anzugeben - feste Effekte - zufällige Effekte.

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
quelle
1
Vielen Dank, Rob, für deinen Beitrag. Genau das versuche ich mit meinen Daten zu tun. Ich verstehe nicht, was var1 im nestedModel auf Asym ist und wie Sie es berechnet haben?
Dies ist nur ein Beispiel dafür, wie die Auswirkung einer Variablen auf Asym berücksichtigt wird: "Das folgende ist das logistische Wachstumsmodell mit einer festen Auswirkung von var1 (binär) auf Asym." Zum Beispiel haben Sie die Variable "Treated", die zwei Werte 0 und 1 hat, also ersetzen Sie "var1" durch "Treated".
PA6OTA