Como encaixo um modelo de efeitos mistos não lineares para dados de medidas repetidas usando nlmer ()?

Estou tentando analisar dados de medidas repetidas e estou lutando para fazê-lo funcionar R. Meus dados são essencialmente os seguintes, eu tenho dois grupos de tratamento. Cada sujeito de cada grupo é testado todos os dias e recebe uma pontuação (a porcentagem correta em um teste). Os dados estão no formato longo:

Time Percent Subject   Group
   1       0    GK11 Ethanol
   2       0    GK11 Ethanol
   3       0    GK11 Ethanol
   4       0    GK11 Ethanol
   5       0    GK11 Ethanol
   6       0    GK11 Ethanol

Os dados se assemelham a uma curva logística, os indivíduos se saem muito mal por alguns dias, seguidos de uma rápida melhoria, seguida de um platô. Eu gostaria de saber se o tratamento afeta a curva de desempenho do teste. Meu pensamento era usar nlmer()no lme4pacote R. Posso ajustar linhas para cada grupo usando o seguinte:

print(nm1 <- nlmer(Percent ~ SSlogis(Time,Asym, xmid, scal) ~ Asym | Subject,
salinedata, start = c(Asym =.60,  xmid = 23, scal = 5)), corr = FALSE)

Posso comparar grupos analisando as estimativas para os diferentes parâmetros e desvios padrão das linhas estimadas, mas não tenho certeza se essa é a maneira correta de fazer isso. Qualquer ajuda seria muito apreciada.

Você pode usar testes de razão de verossimilhança normal. Aqui está um exemplo simples. Primeiro, vamos criar observações de 10 indivíduos com base em seus parâmetros:

Asym = .6
xmid = 23
scal = 5

n = 10
time = seq(1,60,5)

d = data.frame(time=rep(time,10),
               Asym, xmid, scal, group=0)
d$subj = factor(rep(1:n, each=length(time)))

Agora deixe metade deles ter diferentes assíntotas e parâmetros de ponto médio:

ind = (nrow(d)/2):nrow(d)
d$Asym[ind] = d$Asym[ind] + .1
d$xmid[ind] = d$xmid[ind] + 10
d$group[ind] = 1
d$group=factor(d$group)

Podemos simular valores de resposta para todos os indivíduos, com base no modelo:

set.seed(1)
d = transform(d, y = Asym/(1+exp((xmid-time)/scal)) +
                     rnorm(nrow(d), sd=.04))
library(lattice)
xyplot(y~time | group, group=subj,
       data=d, type=c("g","l"), col="black")

Podemos ver diferenças claras entre os dois grupos, diferenças que os modelos devem poder captar. Agora vamos primeiro tentar ajustar um modelo simples , ignorando grupos:

> fm1 = nls(y ~ SSlogis(time, Asym, xmid, scal), data=d)
> coef(fm1)
      Asym       xmid       scal 
 0.6633042 28.5219166  5.8286082

Talvez como esperado, as estimativas Asyme xmidestejam em algum lugar entre os valores reais dos parâmetros para os dois grupos. (Que este seria o caso não é óbvio, porém, já que o parâmetro de escala também é alterado, para ajustar a especificação incorreta do modelo.) Agora vamos ajustar um modelo completo , com parâmetros diferentes para os dois grupos:

> fm2 = nls(y ~ SSlogis(time, Asym[group], xmid[group], scal[group]),
          data=d,
          start=list(Asym=rep(.6,2), xmid=rep(23,2), scal=rep(5,2)))
> coef(fm2)
    Asym1     Asym2     xmid1     xmid2     scal1     scal2 
 0.602768  0.714199 22.769315 33.331976  4.629332  4.749555

Como os dois modelos estão aninhados, podemos fazer um teste de razão de verossimilhança:

> anova(fm1, fm2)
Analysis of Variance Table

Model 1: y ~ SSlogis(time, Asym, xmid, scal)
Model 2: y ~ SSlogis(time, Asym[group], xmid[group], scal[group])
  Res.Df Res.Sum Sq Df  Sum Sq F value    Pr(>F)    
1    117    0.70968                                 
2    114    0.13934  3 0.57034  155.54 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

O valor p extremamente pequeno mostra claramente que o modelo simples era muito simples; os dois grupos não diferem em seus parâmetros.

No entanto, as duas estimativas de parâmetros de escala são quase idênticas, com uma diferença de apenas 0,1. Talvez precisemos apenas de um parâmetro de escala? (É claro que sabemos que a resposta é sim, pois simulamos dados.)

(A diferença entre os dois parâmetros de assíntota também é apenas 0,1, mas é uma grande diferença quando levamos em consideração os erros padrão - veja summary(fm2).)

Então, ajustamos um novo modelo, com um scaleparâmetro comum para os dois grupos, mas diferentes Asyme xmidparâmetros, como antes:

> fm3 = nls(y ~ SSlogis(time, Asym[group], xmid[group], scal),
          data=d,
          start=list(Asym=rep(.6,2), xmid=rep(23,2), scal=5))
> coef(fm3)
     Asym1      Asym2      xmid1      xmid2       scal 
 0.6035251  0.7129002 22.7821155 33.3080264  4.6928316

E como o modelo reduzido está aninhado no modelo completo, podemos novamente fazer um teste de razão de verossimilhança:

> anova(fm3, fm2)
Analysis of Variance Table

Model 1: y ~ SSlogis(time, Asym[group], xmid[group], scal)
Model 2: y ~ SSlogis(time, Asym[group], xmid[group], scal[group])
  Res.Df Res.Sum Sq Df     Sum Sq F value Pr(>F)
1    115    0.13945                             
2    114    0.13934  1 0.00010637   0.087 0.7685

O valor- p grande indica que o modelo reduzido se encaixa bem como o modelo completo, conforme o esperado.

É claro que podemos fazer testes semelhantes para verificar se diferentes valores de parâmetros são necessários para apenas Asym, apenas xmidou ambos. Dito isto, eu não recomendaria fazer uma regressão gradual como esta para eliminar parâmetros. Em vez disso, basta testar o modelo completo ( fm2) em relação ao modelo simples ( fm1) e ficar satisfeito com os resultados. Para quantificar quaisquer diferenças, os gráficos serão úteis.