ANOVA fatorial de medidas repetidas

Olá prof. Marcos,

Fiz uma simulação de dados básica para nós experimentarmos rodar a análise. Os dados seguem abaixo, qualquer problema em abrir, me avisa que eu envio por outra forma.
https://www.dropbox.com/s/7uz2tmtc7kli893/Exclus%C3%A3o.txt?dl=0

Então, tenho um experimento de exclusão, onde fiz cercas para isolar herbívoros de médio porte na Caatinga (bode). Tenho várias áreas, em cada uma um bloco com e sem herbívoros e áreas com pressão de pastejo baixa, média e alta. Além disso, como a sazonalidade na Caatinga é marcante e a dinâmica anual da vegetação muito variável, faço 4 medidas ao ano para acompanhar toda essa dinâmica. Até agora tenho seis medidas, ou seja, um ano e meio de experimento. Os dados são a média de biomassa por área.

Gostaria de fazer uma análise separada, considerando apenas a presença/ausência de herbívoros e, depois, adicionar o fator de pressão de pastejo.

Uma outra dúvida, para a análise daqui eu retirei, mas nos dados reais em algumas áreas mais secas ocorre muitos zeros. Como isso pode atrapalhar a análise e o que deveria fazer para contornar isso? Será que daria problema?

Bom! É isso! Ansiosa para dar uma conferida em como está indo o experimento.

Oi, Tatiane, desculpe pela demora!

Então, o desenho experimental "pede" mesmo uma anova de medidas repetidas, o que não é difícil de fazer com os seus dados. O único inconveniente é que vamos ter que mudar o formato, mas é coisa simples de se fazer: é que a variável resposta terá que ficar toda em uma coluna, e cada repetição deve ser identificada como uma nova variável (tempo, ou como preferir chamar).

Essa coisa do formato na verdade é confusa, pois algumas funções do R permitem analisar no formato que você colocou, e outras não. :p Mas para usarmos as nossas funções mais familiares, como a aov, vamos precisar do formato novo, certo?

Vou começar a elaborar aqui o exemplo, começando com os dados no formato que eu quero:

Código:

dados<-structure(list(Área = c(1L, 1L, 2L, 2L, 3L, 3L, 4L, 4L, 5L, 5L,
6L, 6L, 7L, 7L, 8L, 8L, 9L, 9L, 10L, 10L, 11L, 11L, 12L, 12L,
1L, 1L, 2L, 2L, 3L, 3L, 4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L, 8L, 8L,
9L, 9L, 10L, 10L, 11L, 11L, 12L, 12L, 1L, 1L, 2L, 2L, 3L, 3L,
4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L, 8L, 8L, 9L, 9L, 10L, 10L, 11L,
11L, 12L, 12L, 1L, 1L, 2L, 2L, 3L, 3L, 4L, 4L, 5L, 5L, 6L, 6L,
7L, 7L, 8L, 8L, 9L, 9L, 10L, 10L, 11L, 11L, 12L, 12L, 1L, 1L,
2L, 2L, 3L, 3L, 4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L, 8L, 8L, 9L, 9L,
10L, 10L, 11L, 11L, 12L, 12L, 1L, 1L, 2L, 2L, 3L, 3L, 4L, 4L,
5L, 5L, 6L, 6L, 7L, 7L, 8L, 8L, 9L, 9L, 10L, 10L, 11L, 11L, 12L,
12L), Tratamento = structure(c(2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L), .Label = c("Com", "Sem"), class = "factor"),
    Pressão = structure(c(2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L,
    3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L,
    2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L,
    3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L,
    2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L,
    3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L,
    2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L,
    3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L,
    2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L,
    3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L
    ), .Label = c("alta", "baixa", "media"), class = "factor"),
    Tempo = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
    3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
    4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L,
    5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
    5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
    6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L), .Label = c("T0",
    "T1", "T2", "T4", "T5", "T6"), class = "factor"), Biomassa = c(26.29,
    19.64, 45.27, 20.05, 25.91, 0, 14.65, 38.63, 75.17, 97.13,
    30.21, 10.2, 19.98, 14.2, 57.87, 51.95, 13.29, 7.3, 20.58,
    36.22, 34.23, 23.6, 13.38, 11.48, 23.45, 31.57, 0, 15.43,
    49.83, 14.5, 35.6, 24.5, 55.39, 59.33, 47.1, 41.94, 57.86,
    32.5, 23.84, 32.02, 53.7, 29.03, 65.62, 9.57, 35.76, 23.49,
    6, 6.87, 30.31, 12.44, 59.18, 65.02, 41.71, 3.6, 117.63,
    87.4, 166.05, 117.63, 87.65, 64.35, 16.43, 54.33, 67.63,
    45.58, 42.55, 5.67, 45.6, 30.29, 67.89, 35.76, 35.4, 22.3,
    17.69, 9.04, 25.67, 7.8, 10.31, 5.88, 15.79, 10.23, 14.49,
    6.96, 33.86, 27.65, 22.26, 23.45, 10.54, 5.21, 12.76, 1.31,
    6.29, 10.62, 29.93, 15.75, 19.1, 18.35, 28.1, 16.45, 24.27,
    10.54, 18, 24.34, 60.11, 40.95, 30.29, 25.67, 69.81, 20,
    59.39, 47.54, 14.98, 6.04, 15.76, 10.69, 28.52, 0, 60.5,
    45.72, 37.91, 9.44, 55.65, 27.84, 62.37, 21.15, 66.27, 30.75,
    125.76, 32.29, 70.54, 29.77, 95.42, 59.96, 114.55, 78.54,
    47.22, 26.11, 60.48, 29.71, 68.76, 45.65, 75.69, 55.65, 45.65,
    20.49)), .Names = c("Área", "Tratamento", "Pressão", "Tempo",
"Biomassa"), class = "data.frame", row.names = c(NA, -144L))




Continuo abaixo, aguarde.

Certo, vamos lá!

Código:

summary(dados)

#Entendeu o novo formato? Olhe só:
head(dados)

attach(dados)

#Agora é só criar os modelos! Vamos começar com o mais simples, que você queria
#Olhando apenas o efeito da presença/ausência:

resultado1<-aov(Biomassa~Tratamento*Tempo + Error(factor(Área)))

#Vamos ver os resultados:
summary(resultado1)

Aqui estão os resultados:

Error: factor(Área)
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 11 23093 2099

Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Tratamento 1 8566 8566 21.838 7.77e-06 ***
Tempo 5 31391 6278 16.006 4.36e-12 ***
Tratamento:Tempo 5 4012 802 2.046 0.077 .
Residuals 121 47463 392
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Fica meio torto quando colamos aqui no fórum, mas acho que dá pra ver. Reproduza aí no R para ver nele, vai ficar melhor.

O tratamento tem um efeito, e o Tempo também. Não há interação entre os dois, apesar de passar beeeeeem perto de ter.

Se não há interação, o modelo está te dizendo duas coisas: primeiro, que a biomassa é afetada pelos tratamentos; segundo, que ela muda ao longo do tempo.

Como não ha interações, você pode avaliar tudo isso de maneira separada, fazendo gráficos de médias individuais, um para o tratamento e outro para o tempo, por exemplo. Mas também pode fazer tudo junto via sciplot, se preferir:

Código:

library(sciplot)
lineplot.CI(Tempo, Biomassa, group=Tratamento)

Acho que fica bem legal assim, economiza tempo fazendo tudo no mesmo gráfico.

Agora, avançando e detalhando um pouco:

- no comando da anova, a função Error() foi usada para indicar que existe a dependência. Ou seja, que as áreas se repetem, e por isso temos o modelo de medidas repetidas. Como Área estava codificada em números, eu também usei a função factor() para que ela fosse reconhecida como um fator. Se ela estivesse com classes (tipo um, dois, três, por extenso mesmo), bastaria o Error(Área), certo?

- para a outra análise, é tudo igual, basta adicionar a nova variável com mais um '*'. Só que aí teremos uma interação tripla, e isso pode dar trabalho pra interpretar se ela for significativa. Mas aí é outro caso, por enquanto basta saber montar o modelo.

- se você tiver muitos zeros, certamente teremos problemas com a normalidade dos resíduos, o que pode comprometer a análise. Mas o bacana é que esta mesma lógica de modelo pode ser replicada em um GLM, no qual você pode trabalhar com outras distribuições.

Qualquer coisa vamos nos falando aqui no tópico para elaborar mais, certo?

Olá Pessoal,
Trabalho com ecologia vegetal, mais especificamente comunidade arbóreas de um fragmento de Mata Atlântica. Preciso fazer um Anova de Medidas Repetidas no seguinte desenho experimental: Uma área de estudo com três tratamentos (Interior – Controle; Borda 1; Borda 2). Dentro de cada tratamento tenho 4 parcelas, em cada uma das parcelas eu faço medidas de DAP de cada árvore anualmente (série temporal com 10 medidas para cada árvore). Tenho que avaliar a Dinâmica da Comunidade, para isso preciso saber o se existe diferença entre os anos e entre tratamentos.

Código:



###### PRIMEIRO EXEMPLO EU UTILIZO OS DADOS DOS INDIVÍDUOS:
##### DADOS UTILIZADOS:
##### ID - INDIVIDUO (ÁRVORE)
##### ÁREA - TRATAMENTOS ( IN- INTERIOR / RE – B1 / GA – B2 )
##### PARCELA - AMOSTRA DENTRO DOS TRATAMENTOS (INTERIOR - 1,2,3,4 ; B1 – 1,2,3,4 ); B2 (1,2,3,4)
##### TEMP - ANO DE AMOSTRAGEM (2001/ 2012 )
##### DAP - (TAMANHO DO DIÂMETRO DA ÁRVORE) - VARIÁVEL CONTÍNUA

biom.dados<-read.table("dados_final3.txt",header=T,dec=",")
biom.dados

######### MODELO DOS DADOS #############
id area parcela temp dap
1 IN EC 1 18,30
1 IN EC 12 18,80
2 IN EC 1 17,20
2 IN EC 12 NA
3 IN EC 1 28,80
3 IN EC 12 NA
4 IN EC 1 16,70
4 IN EC 12 18,30
5 IN EC 1 14,00
5 IN EC 12 16,60
6 IN EC 1 13,20
6 IN EC 12 13,50
7 IN EC 1 35,00
7 IN EC 12 36,50
8 IN EC 1 12,00
8 IN EC 12 12,10
9 IN EC 1 10,50
9 IN EC 12 10,60
10 IN EC 1 11,00
10 IN EC 12 14,00
11 IN EC 1 79,90
11 IN EC 12 81,50
12 IN EC 1 31,30
12 IN EC 12 33,80
13 IN EC 1 24,95
13 IN EC 12 25,40
14 IN EC 1 10,20
14 IN EC 12 NA
15 IN EC 1 10,00
15 IN EC 12 NA

################## ANOVA ANINHADA - 1 VARIÁVEL - DADOS DAP ######
anos<-as.factor(biom.dados$temp)         # Transformando anos em fator #
area<-as.factor(biom.dados$area)           # Transformando área em fator #
parcela<-as.factor(biom.dados$parcela) # Transformando parcela em fator #
ind<-as.numeric(biom.dados$id)               # Transformando id em número #
dap<-as.numeric(biom.dados$dap)           # Transformando dap em número #

################ MODELOS DE ANOVA REPETIDAS – INDIVÍDUOS #######
resultado1<-aov(dap~ind*anos*area + Error(ind/area))
resultado2<-aov(dap~ind*anos*area + Error(ind/anos))
resultado3<-aov(dap~ind*(anos*area) + Error(ind/(anos*area)))

#######################################################

###### SEGUNDO EXEMPLO EU UTILIZO OS DADOS DAS PARCELAS:
##### DADOS UTILIZADOS:
##### PARCELA - AMOSTRA DENTRO DOS TRATAMENTOS (INTERIOR - 1,2,3,4 ; B1 – 1,2,3,4 ); B2 (1,2,3,4)
##### ÁREA - TRATAMENTOS ( IN - INTERIOR / RE - REDE ELETRICA / GA - GASODUTO)
##### DAP (TAMANHO DO DIÂMETRO DA ÁRVORE) - VARIÁVEL CONTÍNUA
##### DENSIDADE (N° DE INDIVIDUOS/AREA ) - VARIÁVEL CONTÍNUA
##### TEMP - ANO DE AMOSTRAGEM ( 2001/ 2012)

biom.anos<-read.table("DAP_NUPEM_TEST.txt",header= T,dec="," , fill=TRUE)
biom.anos

######### MODELO DOS DADOS #############
PARCELA AREA DAP DENS ANO
1 IN 24.16 157.5 2001
2 IN 23.7 150 2001
3 IN 21.8 167.5 2001
4 IN 24.4 135 2001
1 IN 26.88 130 2012
2 IN 25.44 145 2012
3 IN 22.37 167.5 2012
4 IN 25.69 122.5 2012
5 B1 18.41 167.5 2001
6 B1 16.89 252.5 2001
7 B1 18.57 227.5 2001
8 B1 24.78 167.5 2001
5 B1 19.98 252.5 2012
6 B1 19.42 235 2012
7 B1 17.64 225 2012
8 B1 24.54 162.5 2012
9 B2 15.5 190 2001
10 B2 19.8 152.5 2001
11 B2 19.4 252.5 2001
12 B2 17.9 212.5 2001
9 B2 18.8 160 2012
10 B2 22.2 122.5 2012
11 B2 19.2 217.5 2012
12 B2 18.4 202.5 2012

############### ANOVA ANINHADA - 2 VARIÁVEL - DADOS DAP X DENSIDADE ############

anos<-as.factor(biom.anos$ANO) # Transformando anos em fator #
area<-as.factor(biom.anos$AREA) # Transformando area em fator #
dap<-as.numeric(biom.anos$DAP) # Transformando dap em numero #
dens<-as.numeric(biom.anos$DENS) # Transformando densidade em numero #
parcela<-as.factor(biom.anos$PARCELA) # Transformando parcela em fator #

################ MODELOS DE ANOVA REPETIDAS - PARCELA #######################
bio.anova.1<-aov(dap~parcela*anos+ Error(parcela/anos))
bio.anova.2<-aov(dap~parcela*anos + Error(area))
bio.anova.3<-aov(dap~parcela*anos + Error(anos/parcela))



Minha dúvida é saber sobre o fator Error. Você pode explicar como definir esse fator?

Professor Marcos, apresentei para dois professores que trabalham na área de estatística do programa de pós que estudo a possibilidade de trabalhar os meus dados de medidas repetidas no tempo com dois tratamentos com ANOVA de medidas repetidas. Eles disseram que não conhecem referência na literatura de ANOVA de medidas repetidas e que fere os pressusposos do teste de ANOVA.
Já vi em algumas publicações da área de pesquisa que estou desenvolvendo em que usam ANOVA para medidas repetidas no tempo, mas não vi citação de nenhuma referência que me dê um suporte para que eu possa defender essa idéia.
Pode falar um pouco a respeito disso. Como posso mostrar que não fere os pressupostos da ANOVA e ter respaldo da literatura que é possível fazer isso?

robertamfs escreveu:Professor Marcos, apresentei para dois professores que trabalham na área de estatística do programa de pós que estudo a possibilidade de trabalhar os meus dados de medidas repetidas no tempo com dois tratamentos com ANOVA de medidas repetidas. Eles disseram que não conhecem referência na literatura de ANOVA de medidas repetidas e que fere os pressusposos do teste de ANOVA.
Já vi em algumas publicações da área de pesquisa que estou desenvolvendo em que usam ANOVA para medidas repetidas no tempo, mas não vi citação de nenhuma referência que me dê um suporte para que eu possa defender essa idéia.
Pode falar um pouco a respeito disso. Como posso mostrar que não fere os pressupostos da ANOVA e ter respaldo da literatura que é possível fazer isso?

Minha dúvida também é essa: se há dependência entre as coletas no tempo, como é possível inserir em uma ANOVA comum???
Porque informar Error(dados$Area) e não Error(dados$Tempo)?

Nessa fiquei confuso....

Ca

Schiedeck escreveu:

Minha dúvida também é essa: se há dependência entre as coletas no tempo, como é possível inserir em uma ANOVA comum???
Porque informar Error(dados$Area) e não Error(dados$Tempo)?

Nessa fiquei confuso....

A variável tempo foi usada como um fator, pois no caso queremos entender se o efeito do tratamento muda ao longo do tempo e/ou interage com ele.

Já a variável área foi usada dentro de "error", porque as repetições ocorrem dentro de cada área - ou seja, é dentro de cada área que há dependência entre as medidas tomadas.

Pessoal,

Tentando aplicar meus recentes conhecimentos adquiridos no curso, podem me ajudar na construção de uma linha de raciocínio?

Tenho um conjunto de dados com 7 variáveis resposta categóricas (binária - presença e ausência), 1 variável resposta quantitativa (índice pluviométrico) e 1 variável explicativa (época do ano). A pergunta é: qual o melhor teste estatístico para avaliar a relação entre essas diferentes variáveis?
ANOVA FATORIAL ou QUI-QUADRADO?

A questão principal que eu quero saber é como a variável quantitativa interfere nas variáveis categóricas, levando em consideração a época do ano (variável explicativa categórica)

Oi, Carolina!

Fiquei com uma dúvida: a variável quantitativa é explicativa? Acho que você se confundiu, pois primeiro escreveu como se ela fosse resposta - mas a continuação do raciocínio me faz pensar que ela seja explicativa.

Se for isso mesmo, acho que um bom ponto de partida é trabalhar com um glm que seria apenas uma pequena "extensão" de uma regressão logística (já que a resposta é binária). Você poderia construir o modelo para cada variável resposta, usando época do ano e precipitação (e a interação entre eles) como explicativas. Neste sentido, se aproxima também de uma anova fatorial, mas a resposta não é quantitativa.

Olá Professor

Sim, me confundi no começo, sobe a variável quantitativa é explicativa e não resposta.

Vamos tentar trabalhar então com modelo GLM, qualquer coisa grito aqui de novo

Muito obrigada pelo norte.

Abraços.

Pessoal,

Como eu faço para visualizar meu lineplot, quando tenho muitos pontos no eixo X? Eu Não consigo manter todos os pontos em um mesmo gráfico devido a quantidade, é possível alterar essa formatação?
O Summary dos meus dados é:
individuo dia_hr concentracao
fantastico:48 1do12 : 6 Min. : 1.59
halley :48 1do16 : 6 1st Qu.: 17.32
imperio :48 1do20 : 6 Median : 28.59
lafuria :48 1do8 : 6 Mean : 33.43
vietna :48 1qf12 : 6 3rd Qu.: 44.08
vigor :48 1qf16 : 6 Max. :143.50
(Other):252 NA's :17

Oi, Carolina!

Às vezes trabalhar com gráficos assim é um imenso desafio mesmo...

Meus primeiros chutes:
- tentar omitir os pontos e usar apenas as linhas OU diminuir consideravelmente o tamanho dos pontos, sem omiti-los mas deixando-os bem pequenos mesmo, com as linhas ganhando mais destaque.
- não usar tudo em um gráfico só, e montar um painel com 6 gráficos (3 linhas por duas colunas deve ficar bom). Se tentar esta opção, lembre-se de manter a mesma escala de eixo y para todos, isto é essencial.

Obrigada professsorrr, vou tentar