Regresi efek campuran non-linear dalam R

14

Anehnya, saya tidak dapat menemukan jawaban untuk pertanyaan berikut menggunakan Google:

Saya memiliki beberapa data biologis dari beberapa individu yang menunjukkan perilaku pertumbuhan sigmoid secara kasar pada waktunya. Jadi, saya ingin memodelkannya menggunakan pertumbuhan logistik standar

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

dengan p0 menjadi nilai awal pada t = 0, k menjadi batas asimptotik pada t-> infinity dan r menjadi kecepatan pertumbuhan. Sejauh yang saya bisa lihat, saya dapat dengan mudah memodelkan ini menggunakan nls (kurangnya pemahaman pada bagian saya: mengapa saya tidak bisa memodelkan sesuatu yang serupa menggunakan regresi logit standar dengan penskalaan waktu dan data? EDIT: Terima kasih Nick, tampaknya orang melakukannya misalnya untuk proporsi, tetapi jarang http://www.stata-journal.com/article.html?article=st0147 Pertanyaan berikutnya tentang garis singgung ini adalah apakah model tersebut dapat menangani outlier> 1).

Sekarang saya ingin mengizinkan beberapa efek tetap (terutama kategorikal) dan beberapa acak (ID individual dan mungkin juga ID studi) pada tiga parameter k, p0 dan r. Apakah saya cara terbaik untuk melakukan ini? Model SSlogis tampaknya masuk akal untuk apa yang saya coba lakukan, apakah itu benar? Apakah salah satu dari model berikut masuk akal untuk memulai? Saya sepertinya tidak bisa mendapatkan nilai awal yang benar dan memperbarui () hanya berfungsi untuk efek acak, bukan yang tetap - ada petunjuk?

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

Karena saya baru mengenal model campuran non-linear pada khususnya dan model non-linear pada umumnya, saya akan sangat menghargai beberapa rekomendasi membaca atau tautan ke tutorial / FAQ dengan pertanyaan pemula.

Rob Hall
sumber
2
Jika Anda menganggap k sebagai diketahui, Anda dapat mengukur populasi Anda dengan P / k. Jika k adalah sesuatu yang harus diperkirakan, itu saja berarti bahwa masalah Anda bukan regresi logit standar.
Nick Cox
1
Nick terima kasih Ya, pada akhirnya saya percaya bahwa k perlu diperkirakan dan dimasukkan dalam regresi. Ketertarikan saya dalam menggunakan regresi logit adalah murni akademis. Saya pikir ini mungkin model yang bagus untuk memulai sebelum pergi ke pemodelan non-linear, tetapi saya tidak dapat menemukan contoh untuk regresi logit untuk data non-biner menggunakan Google. Saya bertanya-tanya apakah ada beberapa alasan (misalnya asumsi distribusi tentang kesalahan) yang membuatnya menjadi ide yang buruk untuk digunakan misalnya glmer dengan tautan logit dengan data kontinu.
Rob Hall
3
Pemodelan logit untuk respons yang proporsinya terus menerus telah ada selama beberapa waktu, tetapi tampaknya tidak dikenal. Lihat misalnya Baum di stata-journal.com/sjpdf.html?articlenum=st0147 Namun demikian itu bukan situasi Anda. Saya tidak bisa mengomentari implementasi R.
Nick Cox
Terima kasih, Nick untuk tautan menarik ini - yang menjelaskan beberapa hal untuk saya. Sedihnya, sepertinya saya belum bisa mengubah tanggapan Anda. (Kalau-kalau ada orang yang kesulitan mengakses tautan langsung, yang berikut ini berhasil bagi saya: stata-journal.com/article.html?article=st0147 )
Rob Hall
1
Pertumbuhan logistik menyiratkan kurva naik monoton. Jika data tidak cocok, Anda akan mendapatkan kecocokan atau perangkat lunak tidak akan diputar, tergantung pada apa yang Anda lakukan.
Nick Cox

Jawaban:

12

Saya ingin membagikan beberapa hal yang saya pelajari sejak mengajukan pertanyaan ini. nlme tampaknya merupakan cara yang wajar untuk memodelkan efek campuran non-linear dalam R. Mulai dengan model basis sederhana:

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
)

Kemudian gunakan pembaruan untuk meningkatkan kompleksitas model. Parameter awal agak sulit untuk dikerjakan, mungkin perlu sedikit waktu untuk mencari tahu urutannya. Perhatikan bagaimana efek tetap baru untuk var1 pada Asym mengikuti efek tetap reguler untuk 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 tampak lebih kuat terhadap outlier dalam dataset saya dan tampaknya menawarkan konvergensi yang lebih andal untuk model yang lebih kompleks. Namun, tampaknya downside adalah bahwa fungsi kemungkinan yang relevan perlu ditentukan secara manual. Berikut ini adalah model pertumbuhan logistik dengan efek tetap dari var1 (biner) pada Asym. Anda dapat menambahkan efek tetap pada xmid dan scal dengan cara yang serupa. Perhatikan cara aneh menentukan model menggunakan rumus ganda sebagai hasil ~ efek tetap ~ efek acak.

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
sumber
1
Terima kasih Rob untuk posting Anda, ini sebenarnya yang saya coba lakukan dengan data saya. Saya tidak mengerti apa itu var1 di nestedModel di Asym dan bagaimana Anda menghitungnya?
Ini hanyalah sebuah contoh tentang bagaimana memasukkan efek dari beberapa variabel pada Asym: "Berikut ini adalah model pertumbuhan logistik dengan efek tetap dari var1 (biner) pada Asym." Misalnya, Anda memiliki variabel "Diperlakukan" yang memiliki dua nilai 0 dan 1, jadi gantilah "Diperlakukan" untuk "var1".
PA6OTA