A primeira coisa que eu faria, antes de proceder com a análise, é organizar o conjunto de dados. Não é necessário utilizar nem attach
e nem converter as variáveis sempre que for ajustar um modelo diferente:
library(nlme)
library(ggplot2)
dados = read.table("dadosnovo.csv", header = T, sep=",", dec=".")
dados$G <- factor(dados$G)
dados$S <- factor(dados$S)
dados$S1 <- factor(dados$S1)
str(dados)
##'data.frame': 471 obs. of 5 variables:
## $ G : Factor w/ 5 levels "G1","G2","G3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ S : Factor w/ 57 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ S1 : Factor w/ 13 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Peso : num 39.2 40.6 41.9 42.8 43 43.2 42.3 42.9 42.5 42.6 ...
## $ Tempo: int 1 2 3 4 5 6 7 8 9 10 ...
Perceba que o novo data frame dados
possui todas as suas cinco colunas com o tipo de dados que deveriam ter.
Meu segundo passo seria fazer uma análise exploratória gráfica dos dados. Para isso, vou plotar um gráfico de painéis, no qual cada painel é referente a um camundongo suíço e suas cores são referentes ao grupo do qual fazem parte. Como a gaiola não é importante nesta análise, eu não a coloquei em meu gráfico.
ggplot(dados, aes(x = Tempo, y = Peso)) +
geom_line(aes(colour = G)) +
facet_wrap(~ S) +
scale_x_continuous(breaks = seq(3, 12, 3))
Com isso, podemos ver como está o comportamento do peso de cada sujeito durante as 12 semanas do experimento. Inclusive, é possível perceber que nem todos chegaram até o final (embora eu desconfie que tu já saiba disso).
Agora sim podemos proceder com a análise. Como os modelos de 1 a 3 não tiveram problemas, vou direto para a análise do modelo 4. Ao escolher uma estrutura de correlação do tipo ARMA, é necessário definir os graus p e q dos polinômios autorregressivo e média móvel, respectivamente. Não vou entrar no mérito de como fazer isso, mas recomendo o livro Mixed Effects Models and Extensions in Ecology with R (Zuur et al., 2009) para uma discussão a este respeito.
Isto posto, para ajustar uma estrutura de correlação ARMA(1, 1) (ou seja, p=1 e q=1), basta rodar
model4 <- lme(Peso ~ G + Tempo, random=list(S=pdIdent(~1)),
correlation = corARMA(p=1, q=1, form= ~ 1 | S ), data=dados)
summary(model4)
## Linear mixed-effects model fit by REML
## Data: dados
## AIC BIC logLik
## 2069.962 2111.382 -1024.981
##
## Random effects:
## Formula: ~1 | S
## (Intercept) Residual
## StdDev: 1.655606 3.900222
##
## Correlation Structure: ARMA(1,1)
## Formula: ~1 | S
## Parameter estimate(s):
## Phi1 Theta1
## 0.7932553 0.2635397
## Fixed effects: Peso ~ G + Tempo
## Value Std.Error DF t-value p-value
## (Intercept) 43.62625 1.5002678 413 29.078976 0.0000
## GG2 -1.42735 1.7021576 52 -0.838551 0.4056
## GG3 -3.30517 1.6971294 52 -1.947506 0.0569
## GG4 -3.47787 1.6902994 52 -2.057548 0.0447
## GG5 -3.46626 1.6964361 52 -2.043260 0.0461
## Tempo 0.21538 0.0832649 413 2.586674 0.0100
## Correlation:
## (Intr) GG2 GG3 GG4 GG5
## GG2 -0.802
## GG3 -0.801 0.687
## GG4 -0.802 0.689 0.690
## GG5 -0.803 0.687 0.688 0.690
## Tempo -0.361 0.098 0.090 0.082 0.093
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.8640659 -0.3889211 0.1319522 0.7970502 2.8488140
##
## Number of Observations: 471
## Number of Groups: 57
Pronto. O teu modelo está ajustado e agora pode ser comparado com os outros, através de testes de razão de verossimilhança, para que se decida a melhor maneira de modelar estes dados.