Regressão de efeitos mistos não lineares em R

14

Surpreendentemente, não consegui encontrar uma resposta para a seguinte pergunta usando o Google:

Eu tenho alguns dados biológicos de vários indivíduos que mostram um comportamento de crescimento aproximadamente sigmóide com o tempo. Assim, desejo modelá-lo usando um crescimento logístico padrão

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

com p0 sendo o valor inicial em t = 0, k sendo o limite assintótico em t-> infinito er sendo a velocidade de crescimento. Tanto quanto eu posso ver, posso modelar isso facilmente usando nls (falta de entendimento da minha parte: por que não posso modelar algo semelhante usando a regressão de logit padrão ao escalonar tempo e dados? EDIT: Obrigado Nick, aparentemente as pessoas fazem, por proporções, mas raramente http://www.stata-journal.com/article.html?article=st0147 A próxima pergunta sobre essa tangente seria se o modelo pode lidar com valores extremos> 1).

Agora, desejo permitir alguns efeitos fixos (principalmente categóricos) e aleatórios (um ID individual e possivelmente também um ID de estudo) nos três parâmetros k, p0 e r. O nlme é a melhor maneira de fazer isso? O modelo SSlogis parece sensato para o que estou tentando fazer, está correto? Um dos seguintes é um modelo sensato para começar? Parece que não consigo acertar os valores iniciais e update () parece funcionar apenas para efeitos aleatórios, não fixos - alguma dica?

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 sou novo nos modelos mistos não lineares em particular e não lineares em geral, gostaria de receber algumas recomendações de leitura ou links para tutoriais / FAQs com perguntas para iniciantes.

Rob Hall
fonte
2
Se você considerar k como conhecido, poderá escalar suas populações em P / k. Se k é algo a ser estimado, isso significa que seu problema não é uma regressão logit padrão.
Nick Cox
1
Obrigado Nick. Sim, no final, acredito que k precisa ser estimado e incluído na regressão. Meu interesse em usar a regressão logit era puramente acadêmico. Eu pensei que esse poderia ser um bom modelo para começar antes de ir para a modelagem não linear, mas não consegui encontrar exemplos de regressão de logit para dados não binários usando o Google. Fiquei me perguntando se havia algum motivo (por exemplo, suposições de distribuição sobre os erros) que tornasse uma má idéia usar, por exemplo, glmer com um link de logit com dados contínuos.
Rob Hall
3
A modelagem do Logit para respostas de proporções contínuas já existe há algum tempo, mas aparentemente não é bem conhecida. Veja, por exemplo, Baum em stata-journal.com/sjpdf.html?articlenum=st0147 No entanto, essa não é a sua situação. Não posso comentar sobre implementações de R.
Nick Cox
Obrigado Nick por esse link interessante - que esclarece algumas coisas para mim. Infelizmente, parece que ainda não posso aprovar sua resposta. (No caso de pessoas têm problemas para acessar o link direto, o seguinte trabalhou para mim: stata-journal.com/article.html?article=st0147 )
Rob Hall
1
O crescimento logístico implica uma curva ascendente monotônica. Se os dados não corresponderem, você terá um ajuste inadequado ou o software não será reproduzido, dependendo exatamente do que você está fazendo.
Nick Cox

Respostas:

12

Eu queria compartilhar algumas das coisas que aprendi desde que fiz essa pergunta. O nlme parece uma maneira razoável de modelar efeitos mistos não lineares em R. Comece com um modelo básico simples:

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
)

Em seguida, use a atualização para aumentar a complexidade do modelo. O parâmetro de início é um pouco complicado de trabalhar, pode levar alguns ajustes para descobrir a ordem. Observe como o novo efeito fixo para var1 no Asym segue o efeito fixo regular do Asym.

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

O lme4 parecia mais robusto contra os outliers no meu conjunto de dados e parecia oferecer uma convergência mais confiável para os modelos mais complexos. No entanto, parece que o lado negativo é que as funções de probabilidade relevantes precisam ser especificadas manualmente. A seguir, é apresentado o modelo de crescimento logístico com um efeito fixo de var1 (binário) no Asym. Você pode adicionar efeitos fixos no xmid e scal de maneira semelhante. Observe a maneira estranha de especificar o modelo usando uma fórmula dupla como resultado - efeitos fixos - efeitos aleatórios.

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
fonte
1
Obrigado Rob pela sua postagem, na verdade é exatamente o que estou tentando fazer com meus dados. Eu não entendo o que é var1 no nestedModel no Asym e como você calculou?
Este é apenas um exemplo de como incluir o efeito de alguma variável no Asym: "A seguir está o modelo de crescimento logístico com um efeito fixo de var1 (binário) no Asym". Por exemplo, você tem a variável "Tratado" que possui dois valores 0 e 1, portanto, substitua "Tratado" por "var1".
PA6OTA 16/03/19