32
Universidade de Brasília IE-Departamento de Estatística Planejamento e Pesquisa 1 Professor George Von Borries Grupo: André Ramos 09/89533 Bruno Santos 09/0107969 João Renato Falcão 10/0107575 Lucas Cusinato 09/122267 Lista Extra PP1 21/05/2012

ANOVA com R

Embed Size (px)

Citation preview

Page 1: ANOVA com R

Universidade de Brasília IE-Departamento de Estatística Planejamento e Pesquisa 1 Professor George Von Borries Grupo: André Ramos 09/89533

Bruno Santos 09/0107969 João Renato Falcão 10/0107575 Lucas Cusinato 09/122267

Lista Extra PP1

21/05/2012

Page 2: ANOVA com R

2

Contents 1)Entrada de Dados ................................................................................................................................... 3

2) Exemplos .............................................................................................................................................. 4

2.1) Exemplo 1- Supressão do crescimento de bactérias em carne armazenada ................................... 4

2.2) Exemplo 2 – Detecção de inflamação em terapia com Amidarone ............................................... 7

2.3) Exemplo 3 - Tamanho de amostra por tratamento ........................................................................ 9

2.4) Exemplo 4 – Tamanho da amostra em função do poder do teste (Letra H) ................................ 11

2.5) Exemplo 5 - ‘Flow Rate Through Filters’ .................................................................................. 13

2.6) Exemplo 6 – ‘Tissue Growth of Tomatoes’ ................................................................................. 16

2.7) Exemplo 7 - ‘Pituitary Function of Hens’ ................................................................................... 18

2.8) Exemplo 8 - Carangueijos .......................................................................................................... 24

3)ANEXO-Programação Completa ........................................................................................................ 26

Page 3: ANOVA com R

3

1)Entrada de Dados Para que o R possa ler os dados é recomendado fazer com o ele leia esses dados de um

arquivo txt ou Excel, já que para colocá-los diretamente na programação ficaria muito confuso e não é muito agradável visualmente. Dessa forma fica mais fácil também fazer alterações nos dados. Para os exemplos seguintes os dados foram colocados em formato txt. Mais detalhes sobre a leitura serão comentados ao longo do trabalho.

Para evitar dificuldades é necessário que o diretório do R esteja com destino na pasta aonde os dados se encontrem. Isso pode ser feito de duas maneiras:

• A primeira é clicando com o botão direito no ícone do R e selecionando “propriedades” e colar o endereço da pasta dos dados em “Start in”.

• A outra a opção seria abrir o programa R e seguir o caminho File->change Dir..

Page 4: ANOVA com R

4

• E então selecionar as pastas onde estão salvos os dados.

2) Exemplos

2.1) Exemplo 1- Supressão do crescimento de bactérias em carne armazenada Primeiro o dados foram colocados em um arquivo txt para que a leitura fosse feita no R. Os dados foram colocados da seguinte maneira:

Ao contrário do SAS, para que as análises possam ser feitas, o R necessita que cada valor observado seja acompanhado do seu respectivo tratamento, como visto acima. É importante que os dados estejam separados por tabulação (TAB) e não por espaço simples! Também se deve verificar o arquivo txt para que não exista espaços ou linhas em branco não desejados e assim evitar problemas futuros.

Abaixo segue a programação para que os dados sejam importados para o R e depois lidos:

Page 5: ANOVA com R

5

A função read.delim lê dados separados por tabulação, existem outras derivações como por exemplo ler dados direto do excel, que seria a função read.csv, mas para os exemplos vamos utilizar sempre a função read.delim. Os argumentos dentro da função são o arquivo que o programa deve ler no caso “ex1.txt’, seguido de hearder=TRUE(ou = T) que indica que a primeira linha dos dados é o cabeçalho, e por último stringsAsFactors=TRUE , diz para o R que os caracteres ou strings são fatores. A variável ‘ex1’ foi onde foram armazenados os dados. Sendo assim quando executamos ex1 e assim podemos ver como R leu os dados. A saída do R ficou da seguinte forma:

Agora que já foi feita a leitura de dados pode-se passar as análises feitas no SAS em sala de aula.

Em seguida vamos fazer uma ANOVA:

Para fazer a ANOVA simplesmente utilizamos o comando aov. Os argumentos da função devem ser a formula, ‘neste caso queremos comparar o número de bactérias por tratamento, em seguida colocamos a base de dados que queremos utilizar, aqui ‘ex1’.

Nota-se que temos apenas a soma quadradas e os graus de liberdade. Para obter a tabela da ANOVA completa basta usar o seguinte comando: summary. Assim temos a seguinte saída:

A partir da anova podemos utilizar o gráfico de resíduos e o QQ-plotpara vermos como

se comportam os dados. Para isso utilizamos a função:

Page 6: ANOVA com R

6

Para que todos os quatro gráficos na mesma página utiliza-se a função par(mfrow=c(2,2)). Que basicamente organiza 2 gráficos por linha e 2 gráficos por coluna.

Para que possamos ter os intervalos para as médias usaremos a função lm, que gera uma regressão com os interceptos sendo nossas próprias médias. Em seguida usaremos a função confint para gerar os intervalos.

Pelo teste F conclui-se que as médias do número de bactérias de cada tratamento são

diferentes.

Page 7: ANOVA com R

7

2.2) Exemplo 2 – Detecção de inflamação em terapia com Amidarone O exemplo 2 segue a mesma ideia do exemplo 1, por isso usaremos basicamente o mesmo código. Primeiro começamos criando o arquivo txt com os dados (ex2.txt) e em seguida o importamos para o R através do comando read.delim.

Lembrando que os argumentos dentro da função são o arquivo que o programa deve ler no caso “ex2.txt’, seguido de hearder=TRUE(ou = T) que indica que a primeira linha dos dados é o cabeçalho, e por último stringsAsFactors=TRUE , diz para o R que os caracteres ou strings são fatores. A variável ‘ex1’ foi onde foram armazenados os dados.

Em seguida fazemos novamente a ANOVA para testar a igualdade de médias através do aov e para termos a tabela completa usamos o comando summary.

Page 8: ANOVA com R

8

Para que possamos ter os intervalos para as médias usaremos a função lm e confint para gerar os intervalos.

Pelo teste F concluímos que as temperaturas médias dos tratamentos são diferentes.

Page 9: ANOVA com R

9

2.3) Exemplo 3 - Tamanho de amostra por tratamento Para este exemplo teremos que baixar um novo pacote para o R, devido ao mesmo não possuir em seu sistema original a função para o calculo do tamanho de amostra para certo nível de poder. Primeiro, vamos a Pacotes>Instalar Pacote(s)... Será aberto uma nova janela onde aparece os diversos servidores de pacotes R, escolhemos um.

Em seguida outra janela aparece mostrando todos os pacotes presentes no servidor escolhido. No nosso caso precisaremos apenas do pacote {pwr}. Em seguida o R abrirá uma janela onde aparecerá o progresso download do pacote. Por ultimo o R imprimira em seu logo se o pacote foi instalado com sucesso e local onde ele se encontra no computador.

Para que possamos usar o pacote em nossa programação será necessário carrega-lo através do comando library. O exemplo já nos fornece alguns dados que são necessários para o calculo do tamanho de amostra. É desejado um poder equivalente a 0.90 para o exemplo então através de erro e tentativa

Page 10: ANOVA com R

10

mudamos o valor de r até acharmos uma quantidade de repetições que geram o poder desejado. Em nosso exemplo vemos que 9 é um tamanho de amostra que gera um resultado próximo o suficiente ao poder desejado.

Page 11: ANOVA com R

11

2.4) Exemplo 4 – Tamanho da amostra em função do poder do teste (Letra H) Importamos os dados presentes no txt para o R, como nos exemplos passados.

Em seguida usamos o comando oav e summary para gerar a ANOVA e conseguirmos o

valor do MSE.

Page 12: ANOVA com R

12

Logo, usando as informações fornecidas e os valores encontramos montamos as funções para tentar encontrar o número de cruzamentos necessários para a rejeição da hipótese nula. Ao final encontramos que 27 cruzamentos são suficientes para rejeitarmos a hipótese nula ao nível de significância de 0.01.

Page 13: ANOVA com R

13

2.5) Exemplo 5 - ‘Flow Rate Through Filters’ Neste exemplo já rejeitamos a hipótese de igualdade de médias do teste F, logo desejamos comparar essas médias a uma média controle e verificar qual o tratamento é mais eficaz. Em nosso exemplo o tratamento F será o nosso controle. Primeiro teremos que baixar e instalar o pacote {multcomp}, como explicado no exemplo 3, para que possamos aplicar o teste de comparação de médias de Dunnett. Em seguida carregamos o pacote com o comando library e importamos os dados para o R através do comando read.delim.

Em seguida calcularemos as medias de cada tratamento com o comando aggregate e pomo nome no cabeçalho com o mando names.

Page 14: ANOVA com R

14

Rodamos agora a ANOVA, oav, para que possamos obter os valores para o calculo do

teste de Dunnett e criamos a matriz de contrates com o comando contrMat.

Note que na matriz de contraste os tratamento aparecem como números invés de letras, logo temos A=1, B=2, C=3, D=4, E=5 e F=6.

Por último, rodamos o teste de Dunnett pelo comando glht e representamos graficamente os intervalos de confiança pelos comandos confint e plot.

Page 15: ANOVA com R

15

Através desses resultados vemos que os tratamentos B e C são praticamente iguais ao tratamento F e que os tratamentos A, D e E são efetivos sendo o E o mais eficaz de todos eles.

Page 16: ANOVA com R

16

2.6) Exemplo 6 – ‘Tissue Growth of Tomatoes’ Este exemplo segue a mesma ideia do exemplo anterior. Queremos fazer comparações de médias de tratamentos com um tratamento controle e ver a eficácia. Começamos carregando o pacote {multcomp} com o comando library e importando os dados para o R com o comando read.delim.

Em seguida usamos o comando aov para gerar a ANOVA e o comando glht para fazer o teste de Dunnett.

Page 17: ANOVA com R

17

Por último montamos os intervalos de confiança para as médias através do comando confint e os representamos graficamente através do comando plot.

Estes resultados nos mostram que os três tratamentos são efetivos, sendo o C o tratamento mais eficaz.

Page 18: ANOVA com R

18

2.7) Exemplo 7 - ‘Pituitary Function of Hens’ 2.7 Glândulas Pituitárias

Para este exemplo é necessário instalar os pacotes {agricolae} e {multicomp} . Como nos exemplos anteriores é preciso ler o os dados de um arquivo texto. O banco de dados fica da seguinte forma:

Primeiro vamos fazer o boxplot desses dados. Para isso basta utilizar o comando plot() que o boxplt será feito, para colocar título gráfico fazemos plot(gal, main=’Nível de T3 por Tratamento’), o gráfico segue na página seguinte:

Page 19: ANOVA com R

19

Agora temos que fazer a ANOVA para ver se existe diferença significativa entre os tratamentos. Como anteriormente utilizamos o comando aov(), e depois para ver a tabela completa da ANOVA temos que utilizar o comando summary(anova). Obtemos a saída seguinte:

Vemos que o p-valor do teste F foi muito pequeno, logo rejeitamos a hipótese nula, ou seja, não existe igualdade entre as médias.

Para encontrar as médias de quadrados mínimos temos que utilizar a função lm() que faz uma regressão, aonde os coeficiente de cada tratamento será a média de quadrados mínimos para aquele tratamento , temos que dizer que não queremos intercepto , caso contrário o R entende o que o primeiro tratamento é o intercepto, para isso dentro da formula da função colocamos ’-1’.

Page 20: ANOVA com R

20

Em seguida para obter a tabela com a estimativa, o valor de t observado e o p-valor, ultilizamos a função summary():

Para fazer um gráfico com essas médias estimadas, temos que retirar o valor da média da tabela gerada pelo output do R. Essa tabela gerada pelo summary(reg7) é um data frame, então precisamos informar ao R qual a coluna do data frame que queremos, para isso basta colocar o nome do data frame seguido de um ‘$’ e seguida o nome da coluna. Neste caso então o data frame se chama summary(reg7) e a coluna que queremos é a parte ‘coefficients’,que está assinalada acima, logo fazemos summary(reg7)$coefficients obtemos o seguinte:

Agora queremos retirar apenas a coluna ‘Estimate’, porém agora o R considera a parte assinalada como uma matriz, então temos que informar qual a coluna que queremos, para isso utilizamos tab[,1], o que indica que queremos todas as linhas da primeira coluna:

Page 21: ANOVA com R

21

Agora finalmente para fazer o gráfico com as médias utilizamos o comando plot(), mas agora temos que dar nomes ao eixo x, senão o R chamara cada tratamento de 1,2,3,4,5. Então vamos utilizar o argumento xaxt=”n” dentro da função plot()indicando que não queremos que o R nomeie o x. E depois para nomear o eixo x utilizamos a função axis(1,at=1:5,lab=c('60g','80g','Fasting','Mash','Premolt')), sendo o ‘1’ para informar que vamos alterar o eixo x, at=1:5 é para dizer que do primeiro até o quinto ponto vamos dar nomes, e lab=c('60g','80g','Fasting','Mash','Premolt'), são os nomes de cada um dos pontos. Finalmente obtemos o gráfico:

Pela a ANOVA vimos que a hipóteses de igualdade de médias foi rejeitada. Por isso temos que fazer testes de comparações múltiplas. Veremos como são feitos os testes LSD de Fisher, SNK e Tukey. O pacote {agricolae} fornece os dois primeiros testes, o teste de Tukey já é padrão do R.

Infelizmente as funções para o teste LDS e SNK somos obrigado a informar o MSE da tabela ANOVA, isso pode ser feito a mão, ou seja, nós mesmos escrevemos o MSE dentro da função ou podemos fazer com q o R retire o valor da tabela e assim pode-se usar sempre esse programa para diferentes dados. Para fazer com o que o R retire o valor da tabela temos que informa-lo onde está esse valor, lembrando que a tabela ANOVA e dado por summary(anova7)

Page 22: ANOVA com R

22

e o MSE se encontra na terceira coluna e na segunda linha fazemos: mse7=summary(anova7)[[1]][[3]][[2]]. Também precisamos informar os graus de liberdade, pelo mesmo pensamento temos gl7= summary(anova7)[[1]][[1]][[2]].

Vamos fazer os testes LSD e SNK, LSD.test() e SNK.test() , nessas duas funções do pacote {agricolae} temos que informar primeiro a variável resposta, neste caso gal$T3, depois os tratamentos, gal$Tratamento. Em seguida temos que informar o MSE e os graus de liberdade, mse7 e gl7 respectivamente. Obtemos:

Por último vamos fazer o teste de Tukey, temos apenas que informar a anova que estamos analisando e o tratamento:

Page 23: ANOVA com R

23

A seguintes funções são do pacote {multcomp}. Temos que pegar a média de cada grupo então vamos utilizar a função aggregate(), informando de qual variável queremos tirar a média:

Agora com a função contrMat{} monta a matriz de contrastes, colocando pelo método de tukey para que seja possível comparar todas duas a duas:

Agora podemos fazer o teste de comparação múltipla utilizando a matriz de contrastes com a função ghlt(), aonde temos que informar a tabela ANOVA que estamos utilizando, neste caso anova7 e a variável que vamos comparar, aqui o Tratamento:

Page 24: ANOVA com R

24

2.8) Exemplo 8 - Carangueijos Neste exemplo utilizaremos os pacotes {geoR} e {car} para que possamos utilizar o método de Box-Cox. Primeiro baixamos e instalamos estes pacotes e importamos os dados para o R pelo comando read.delim. Note que fazemos uma transformação, somamos um ruído, nos dados para que não tenhamos zeros e assim podemos fazer a transformação sem problemas.

Em seguida carregaremos a função {geoR}e a utilizaremos para calcular o melhor lambda

pra realizar a transformação.

Page 25: ANOVA com R

25

Agora carregando a função {car} fazemos a transformação dos dados.

Page 26: ANOVA com R

26

3)ANEXO-Programação Completa

################################# Exemplo1 ################################# #Lendo os dados ex1=read.delim('ex1.txt', header=T,stringsAsFactors=T ) ex1 #Boxplot plot(Bacterias~Tratamento,data=ex1) #Fazendo a ANOVA anova=aov(Bacterias~Tratamento,ex1) anova summary(anova) #Gráficos de diánostico par(mfrow=c(2,2)) plot(anova) #regressão para obter as médas e intervalo de confiança reg1=lm(Bacterias~Tratamento-1,ex1) inter.ex1=confint(reg1,level=0.95) inter.ex1 #Teste de Tukey tukey=TukeyHSD(anova,'Tratamento') tukey plot(tukey) ################################## Exemplo2 ################################## #Lendo os dados ex2=read.delim('ex2.txt', header=T,stringsAsFactors=T ) ex2 #Boxplot plot(Temp~Tratamento,data=ex2)

Page 27: ANOVA com R

27

#Fazendo a ANOVA anova2=aov(Temp~Tratamento,ex2) anova2 summary(anova2) #Gráficos de diagnótico da ANOVA par(mfrow=c(2,2)) plot(anova2) ################################## Exemplo3 ############################################################# # Lendo o pacote library(pwr) #Informando os dados e montando a função phi mse=0.22 trat=3 r=9 N=trat*r ui=c(0.5,-0.2,-0.3) u=mean(ui) ti=ui-u t=sum(ti^2) a=0.05 fsq= r*t/(N*mse) fsq f=sqrt(fsq) f #Função para obter o poder poder=pwr.anova.test(k=trat,n=r,f=f,sig.level=a) poder #Utilizando o pacote library(pwr) #Lendo os dados ex4=read.delim('ex4.txt', header=T,stringsAsFactors=T ) ex4

Page 28: ANOVA com R

28

#Boxplot plot(Tempo~Tratamento,data=ex4) #Fazendo a ANOVA anova4=aov(Tempo~Tratamento,ex4) anova4 summary(anova4) #Montando uma matriz para deixar todos os gráficos numa saída só par(mfrow=c(2,2)) # gráficos de diagnóstico da ANOVA plot(anova4) # Retirando o MSE da ANOVA, Informando médias e Montando a funçao phi mse4=summary(anova4)[[1]][[3]][[2]] trat4=3 mi=c(20,16,18) m=mean(mi) taui=mi-m tau=sum(taui^2) a4=0.01 #Fazendo o Looping para ir testando quantas repetições até o poder 0.9 rep=0 for(rep in 1:30 ){ rep=rep+1 tr=trat4*rep fisq= rep*tau/(tr*mse4) fi=sqrt(fisq) pod=pwr.anova.test(k=trat4,n=rep,f=fi,sig.level=a4) if(pod$power>=0.9) {break} } pod$n

Page 29: ANOVA com R

29

############################## PARTE 2############################################ ################################################################################## ################################################################################# ####exemplo5 #Baixano o pacote library(multcomp) #Lendo os dados ex5=read.delim('ex5.txt', header=T,stringsAsFactors=T) ex5 #Boxplot plot(ex5) #Pegando a médas de cada tratamento medias= aggregate(ex5[,c('Flowrate')],ex5['Filtro'],mean) names(medias)=c('Filtro','Flowrate') medias #Fazendo a ANOVA anova5=aov(Flowrate~Filtro,data=ex5) summary(anova5) #Contrastes contrat=contrMat(medias$Flowrate,type='Dunnett',base=6) contrat #Teste de Dunnett dun <- glht(anova5, linfct = mcp(Filtro=contrat)) dun #Intravalo de Confiança inter=confint(dun,level=0.95) inter #Gráfico dos intervalos de confiança plot(inter)

Page 30: ANOVA com R

30

####exemplo6### #Lendo os Dados tomates=read.delim('tomate.txt', header=T,stringsAsFactors=T) tomates #Fazendo a ANOVA anova6=aov(resp~trat,data=tomates) summary(anova6) #Contrastes dun.to <- glht(anova6, linfct = mcp(trat='Dunnett')) dun.to #Intervalo de confiança inter.to=confint(dun.to,level=0.95) inter.to plot(inter.to,xlim=c(-20,10)) #########Exemplo7### #Mandando abrir o pacote agricolae library(agricolae) library(multcomp) #Lendo os dados gal=read.delim('gal.txt', header=T,stringsAsFactors=T) gal #Fazendo a ANOVA anova7=aov(T3~Tratamento,gal) summary(anova7) #Fazendo o Boxplot plot(gal,main='Nível de T3 por tipo de Tratamento') #Regressão para MQ reg7=lm(T3~Tratamento-1,gal) summary(reg7) tab=(summary(reg7)$coefficients) tab tab[,1]

Page 31: ANOVA com R

31

#Gráfico dos MQ plot(tab[,1],main='Médias de Mínimos Quadradros por Tratamento',xaxt='n') axis(1,at=1:5,lab=c('60g','80g','Fasting','Mash','Premolt')) #Retirando o MSE e GL mse7=summary(anova7)[[1]][[3]][[2]] mse7 gl7=summary(anova7)[[1]][[1]][[2]] gl7 #Testes de comparção múltipla LSD.test(gal$T3,gal$Tratamento,DFerro=gl7,MSerror=mse7) SNK.test(gal$T3,gal$Tratamento,DFerro=gl7,MSerror=mse7) Tukey7=TukeyHSD(anova7,'Tratamento') Tukey7 # medias7= aggregate(gal[,c('T3')],gal['Tratamento'],mean) names(medias7)=c('Tratamento','T3') medias7 #Contrastes contrat7=contrMat(medias7$T3,type='Tukey') contrat7 #Teste de comparação multipla comp7 <- glht(anova7, linfct = mcp(Tratamento=contrat7)) comp7 ##Exemplo 8 setwd('F:/PP1/trabalho') #exportei o dataset do sas para um "tab delimited" (acho que é isso) file chamado 'crabs boxcox.txt' crab = read.delim('crabs boxcox.txt') crab2 = transform(crab, crabs = crabs+1e-6) #soma um milionesimo aos dados pra nao ter zeros

Page 32: ANOVA com R

32

### o pacote mais simples e mais chato de usar (parece um pouco com o sas) library(MASS) bc.mass = MASS::boxcox(crab2$crabs ~ crab2$habitat, lambda = seq(.1, .2, length=1000)) (bc.mass$x)[which.max(bc.mass$y)] # o lambda maximo do.call(cbind, bc.mass) str(bc.mass) ### funções mais rápidas para determinar o melhor lambda e/ou fazer a transformação library(car) ?car::box.cox(crab2$crabs ~ crab2$habitat) bc.car1 = car::bcPower(crab2$crabs, 0.1416926) #com essa é preciso setar o lambda pra transformar bc.car2 = car::boxCoxVariable(crab2$crabs) #essa calcula o lambda e transforma ja library(geoR) geoR::boxcoxfit(crab2$crabs) #calcula o lambda