4

Oi Galera Tudo bem? Estou com alguns problemas a respeito de definir a estrutura de covariância no R. Eu preciso ajustar as seguintes matrizes de covariância:

  • AR
  • ARMA
  • compSymm
  • Uniforme.

Como vocês podem observar no código abaixo, eu consegui definir as estruturas AR e compSymm (o valor 0.3 coloquei aleatoriamente, posso fazer isso? Se não, como corrigir?), porém a estrutura ARMA não estou conseguindo. Poderiam verificar qual o problema ?

Obs: Alterei a extensão do arquivo da base de dados só para poder postar, mas venho utilizando csv separado por pontos e vírgula.

Observem as tentativas nos códigos abaixo.

rm(list=ls())
library(GGally)
library(splines)
library(nlme)
library(fields)
library(lattice)
require(ISLR)
library(grid)

dados = read.table("dadosnovo.csv", header = T, sep=";", dec=",")
dados$trat=dados$G
attach(dados)
head(dados)
model1 <- lme(Peso~factor(G)+Tempo,random=list(S=pdIdent(~1)),
              data=dados)
summary(model1)

model2 <- lme(Peso~factor(G)+Tempo,random=list(S=pdIdent(~1)),
              correlation = corAR1(form = ~ 1 | S ), data=dados)
summary(model2)

anova(model1,model2)
model3 <- lme(Peso~factor(G)+Tempo,random=list(S=pdIdent(~1)),
              correlation = corCompSymm(0.3, form= ~ 1 | S ), data=dados)
summary(model3)
model4 <- lme(Peso~factor(G)+Tempo,random=list(S=pdIdent(~1)),
              correlation = corARMA(form= ~ 1 | S ), data=dados) #ERROR

0

1 Resposta 1

3

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

inserir a descrição da imagem aqui

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.

0

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .