Regresión de efectos mixtos no lineales en R

Sorprendentemente, no pude encontrar una respuesta a la siguiente pregunta usando Google:

Tengo algunos datos biológicos de varias personas que muestran un comportamiento de crecimiento más o menos sigmoide en el tiempo. Por lo tanto, deseo modelarlo utilizando un crecimiento logístico estándar

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

siendo p0 el valor inicial en t = 0, k es el límite asintótico en t-> infinito yr es la velocidad de crecimiento. Hasta donde puedo ver, puedo modelar esto fácilmente usando nls (falta de comprensión de mi parte: ¿por qué no puedo modelar algo similar usando la regresión logit estándar al escalar el tiempo y los datos? EDITAR: Gracias Nick, aparentemente la gente lo hace, por ejemplo, por proporciones, pero rara vez http://www.stata-journal.com/article.html?article=st0147 . La siguiente pregunta sobre esta tangente sería si el modelo puede manejar valores atípicos> 1).

Ahora deseo permitir algunos efectos fijos (principalmente categóricos) y algunos aleatorios (una identificación individual y posiblemente también una identificación de estudio) sobre los tres parámetros k, p0 y r. ¿Es nlme la mejor manera de hacer esto? El modelo SSlogis parece sensato para lo que estoy tratando de hacer, ¿es correcto? ¿Es cualquiera de los siguientes un modelo sensato para empezar? Parece que no puedo obtener los valores iniciales correctamente y actualizar () solo parece funcionar para efectos aleatorios, no fijos, ¿alguna sugerencia?

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

Como soy nuevo en los modelos mixtos no lineales en particular y en los modelos no lineales en general, agradecería algunas recomendaciones de lectura o enlaces a tutoriales / preguntas frecuentes con preguntas para novatos.

Quería compartir algunas de las cosas que aprendí desde que hice esta pregunta. nlme parece una forma razonable de modelar efectos mixtos no lineales en R. Comience con un modelo base simple:

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
)

Luego use la actualización para aumentar la complejidad del modelo. Es un poco difícil trabajar con el parámetro de inicio, puede tomar algunos retoques para determinar el orden. Observe cómo el nuevo efecto fijo para var1 en Asym sigue el efecto fijo regular para 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 parecía más robusto frente a los valores atípicos en mi conjunto de datos y parecía ofrecer una convergencia más confiable para los modelos más complejos. Sin embargo, parece que la desventaja es que las funciones de probabilidad relevantes deben especificarse manualmente. El siguiente es el modelo de crecimiento logístico con un efecto fijo de var1 (binario) en Asym. Puede agregar efectos fijos en xmid y scal de manera similar. Tenga en cuenta la extraña forma de especificar el modelo utilizando una fórmula doble como resultado ~ efectos fijos ~ efectos aleatorios.

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