————————————————————————————————————————————————————————————————————————-
————————————————————————————————————————————————————————————————————————-
Configura base de dados…
Apaga variáveis do ambiente R…
Carregando bibliotecas R…
# Pacotes que serao utilizados
library(agricolae)
library(ggplot2)
library(asbio) # Teste Tukey de Aditividade
library(doBy)
library(car) # Teste de Fligner Killeen
library(effects)
library(pwr)
————————————————————————————————————————————————————————————————————————-
QUESTÕES:
Desenvolvimento/rascunho:
#Para melhorar a visualização via notação científica dos dados:
options(scipen=10)
# Leitura do arquivo:
tecnicas<-read.table("tecnicas.txt",h=T)
tecnicas
## mix vlr
## 1 M1 3129
## 2 M1 3000
## 3 M1 2865
## 4 M1 2890
## 5 M2 3200
## 6 M2 3300
## 7 M2 2975
## 8 M2 3150
## 9 M3 2800
## 10 M3 2900
## 11 M3 2985
## 12 M3 3050
## 13 M4 2600
## 14 M4 2700
## 15 M4 2600
## 16 M4 2765
str(tecnicas)
## 'data.frame': 16 obs. of 2 variables:
## $ mix: Factor w/ 4 levels "M1","M2","M3",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ vlr: int 3129 3000 2865 2890 3200 3300 2975 3150 2800 2900 ...
summary(tecnicas)
## mix vlr
## M1:4 Min. :2600
## M2:4 1st Qu.:2791
## M3:4 Median :2938
## M4:4 Mean :2932
## 3rd Qu.:3070
## Max. :3300
# Análise Gráfica dos Dados:
# Analise via BOXPLOT das medianas de resistencia para cada caixa:
#boxplot(vlr~mix,data=tecnicas,col="lightgray")
# Boxplot : y x misturas (mix)
#ggplot(tecnicas,aes(x=mix,y= vlr,fill=mix)) + geom_boxplot()
# função aov: para estimação do modelo linear:
mod1<-aov(vlr~mix,data=tecnicas)
#mod1
# Exibe a Tabela da ANAVA:
summary(mod1)
## Df Sum Sq Mean Sq F value Pr(>F)
## mix 3 489740 163247 12.73 0.000489 ***
## Residuals 12 153908 12826
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Estimativas dos Efeitos Principais (Media das misturas) e do Intercepto:
coef(mod1)
## (Intercept) mixM2 mixM3 mixM4
## 2971.00 185.25 -37.25 -304.75
# dados da coluna:
summary(tecnicas$vlr)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2600 2791 2938 2932 3070 3300
# Comando SUMMARY Via biblioteca doBy:
# Com médias:
summaryBy(vlr~mix,data=tecnicas,FUN=c(length,mean))
## mix vlr.length vlr.mean
## 1 M1 4 2971.00
## 2 M2 4 3156.25
## 3 M3 4 2933.75
## 4 M4 4 2666.25
# Com medianas:
summaryBy(vlr~mix,data=tecnicas,FUN=c(length,median))
## mix vlr.length vlr.median
## 1 M1 4 2945.0
## 2 M2 4 3175.0
## 3 M3 4 2942.5
## 4 M4 4 2650.0
# Previsão do valor para cada elemento da amostra:
#predict(mod1)
# Valores reais da amostra:
#tecnicas$vlr
# Diferença (resíduo) entre o valor previsto e o valor real:
#residuals(mod1)
# CV do Pacote Agricolae:
cv.model(mod1)
## [1] 3.862817
#names(mod1)
#mod1
#mod1$qr
# Configura tela para todos os gráficos (caso contrário, vai de um em um):
#par(mfrow = c(2,2))
# Diagnósticos Gráficos:
#plot(mod1)
# Teste de Shapiro-Wilks: (aproximadamente normal)
# Teste de Shapiro quando Ho: tem distribuição aprox. normal. (grandes) amostras:
#shapiro.test(mod1$residuals)
#### Teste de Bartlett : H0 : Variancia dos residuos é homogênea:
#bartlett.test(vlr~mix,data=tecnicas)
### Teste de Fisher com correção de BONFERRONI (para o campo "mix"):
#LSD.test(mod1,"mix",p.adj="bon",console=TRUE)
#Visualização gráfica:
#modeloLSD=LSD.test(mod1,"mix",p.adj="bon",console=TRUE)
#modeloLSD$groups
#bar.group(modeloLSD$groups,ylim=c(0,4000))
#### Teste de Tukey:
#test2 = TukeyHSD(mod1)
#test2
#plot(test2)
#### Teste de Tukey do Pacote Agricolae
#test3 = HSD.test(mod1,"mix",console=TRUE)
#bar.group(test3$groups,ylim=c(0,4000))
Respostas (Questão 01):
Letra (a):
# Boxplot : y x misturas (mix)
ggplot(tecnicas,aes(x=mix,y= vlr,fill=mix)) + geom_boxplot()
# Exibe a Tabela da ANAVA:
summary(mod1)
## Df Sum Sq Mean Sq F value Pr(>F)
## mix 3 489740 163247 12.73 0.000489 ***
## Residuals 12 153908 12826
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Letra (b):
# Teste de Fisher com correção de BONFERRONI (para o campo "mix") via variável modLSD1:
modLSD1=LSD.test(mod1,"mix",p.adj="bon",console=TRUE)
##
## Study: mod1 ~ "mix"
##
## LSD t Test for vlr
## P value adjustment method: bonferroni
##
## Mean Square Error: 12825.69
##
## mix, means and individual ( 95 %) CI
##
## vlr std r LCL UCL Min Max
## M1 2971.00 120.55704 4 2847.624 3094.376 2865 3129
## M2 3156.25 135.97641 4 3032.874 3279.626 2975 3300
## M3 2933.75 108.27242 4 2810.374 3057.126 2800 3050
## M4 2666.25 80.97067 4 2542.874 2789.626 2600 2765
##
## alpha: 0.05 ; Df Error: 12
## Critical Value of t: 3.152681
##
## Least Significant Difference 252.4675
## Means with the same letter are not significantly different.
##
## Groups, Treatments and means
## a M2 3156
## a M1 2971
## a M3 2934
## b M4 2666
bar.group(modLSD1$groups,ylim=c(0,4000))
# Teste de Tukey do Pacote Agricolae
#modTUKEYagr=HSD.test(mod1,"mix",console=TRUE)
#bar.group(modTUKEYagr$groups,ylim=c(0,4000))
# Teste de Tukey (HSD):
TukeyHSD(mod1,ordered=TRUE)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
## factor levels have been ordered
##
## Fit: aov(formula = vlr ~ mix, data = tecnicas)
##
## $mix
## diff lwr upr p adj
## M3-M4 267.50 29.74971 505.2503 0.0261838
## M1-M4 304.75 66.99971 542.5003 0.0115923
## M2-M4 490.00 252.24971 727.7503 0.0002622
## M1-M3 37.25 -200.50029 275.0003 0.9652776
## M2-M3 222.50 -15.25029 460.2503 0.0693027
## M2-M1 185.25 -52.50029 423.0003 0.1493561
modTUKEY1= TukeyHSD(mod1,ordered=TRUE)
plot(modTUKEY1)
Letra (c):
# Configura tela para todos os gráficos (caso contrário, vai de um em um):
par(mfrow = c(2,2))
# Diagnósticos Gráficos:
plot(mod1)
Pela análise do gráfico quantil-quantil (Normal Q-Q), observamos uma distribuição normal dos resíduos das amostras.
Afirmativa lastreada pelos resultados do teste de Shapiro sobre os resíduos do modelo (vide abaixo):
# Teste de Shapiro-Wilks: (análise dos resíduos com distribuição aproximadamente normal)
shapiro.test(mod1$residuals)
##
## Shapiro-Wilk normality test
##
## data: mod1$residuals
## W = 0.97046, p-value = 0.846
Letra (d):
# Configura tela para todos os gráficos (caso contrário, vai de um em um):
par(mfrow = c(2,2))
# Valores reais das amostras:
plot(tecnicas$vlr,main="Valores reais das amostras")
abline(h=mean(tecnicas$vlr),col="blue",lwd=1)
# Previsão do valor para cada elemento das amostras:
plot(predict(mod1),main="Previsão (resíduos no modelo)")
abline(h=mean(predict(mod1)),col="blue",lwd=1)
# Diferença (resíduo) entre o valor previsto e o valor real:
plot(residuals(mod1),main="Diferença (real x previsto)")
abline(h=mean(residuals(mod1)),col="blue",lwd=1)
# Gráfico da distribuição normal:
qqnorm(residuals(mod1),main="Distribuição Normal")
qqline(residuals(mod1))
Letra (e):
————————————————————————————————————————————————————————————————————————-
Desenvolvimento/rascunho:
# Leitura do arquivo:
refino<-read.table("refino.txt",h=T)
refino
## vel frn vlr
## 1 V05 F1 8
## 2 V05 F2 4
## 3 V05 F3 5
## 4 V05 F4 6
## 5 V10 F1 14
## 6 V10 F2 5
## 7 V10 F3 6
## 8 V10 F4 9
## 9 V15 F1 14
## 10 V15 F2 6
## 11 V15 F3 9
## 12 V15 F4 2
## 13 V20 F1 17
## 14 V20 F2 9
## 15 V20 F3 3
## 16 V20 F4 6
str(refino)
## 'data.frame': 16 obs. of 3 variables:
## $ vel: Factor w/ 4 levels "V05","V10","V15",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ frn: Factor w/ 4 levels "F1","F2","F3",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ vlr: int 8 4 5 6 14 5 6 9 14 6 ...
summary(refino)
## vel frn vlr
## V05:4 F1:4 Min. : 2.000
## V10:4 F2:4 1st Qu.: 5.000
## V15:4 F3:4 Median : 6.000
## V20:4 F4:4 Mean : 7.688
## 3rd Qu.: 9.000
## Max. :17.000
# Comando SUMMARY Via biblioteca doBy:
# Com médias:
summaryBy(vlr~vel,data=refino,FUN=c(length,mean))
## vel vlr.length vlr.mean
## 1 V05 4 5.75
## 2 V10 4 8.50
## 3 V15 4 7.75
## 4 V20 4 8.75
# Com medianas:
summaryBy(vlr~vel,data=refino,FUN=c(length,median))
## vel vlr.length vlr.median
## 1 V05 4 5.5
## 2 V10 4 7.5
## 3 V15 4 7.5
## 4 V20 4 7.5
# Boxplot : y x fornos
#ggplot(refino,aes(x=frn,y=vlr,fill=frn)) + geom_boxplot()
#Transformando as colunas vel e frn em fatores:
#vel<-factor(refino$vel)
#frn<-factor(refino$frn)
# função aov: para estimação do modelo linear com dois fatores:
#mod2<-aov(vlr~vel+frn+vel:frn,data=refino)
#mod2<-aov(vlr~vel+frn,data=refino)
#mod2
# Exibe a Tabela da ANOVA:
#summary(mod2)
# Verificando Interacao entre Blocos e Niveis do Fator (velocidade nos fornos x graos)
#interaction.plot(refino$vel,refino$frn,refino$vlr,fixed=TRUE)
# Verificando Interacao entre Blocos e Niveis do Fator (fornos x velocidade x graos)
#interaction.plot(refino$frn,refino$vel,refino$vlr,fixed=TRUE)
# Teste de Aditividade de Tukey
# Testa H0 : Nao ha efeito interacao
#tukey.add.test(refino$vlr,refino$vel,refino$frn)
# Interpretação do TESTE TUKEY ACIMA: como o p-value é alto (maior que o nivel de significancia de 5%, pois está dando 15...%), a chance de encontrar um F ALTO é baixa, sendo assim, aceita-se a hipotese nula (pois se em outra situação, o F é alto em um teste de hipoteses, você rejeita a nula)
# Estimacao do Modelo como se fosse ANAVA 1 Fator (velocidade):
#modelo<-aov(vlr~vel,data=refino)
# Exibe a Tabela da ANAVA:
#summary(modelo)
# Estimacao do Modelo com blocagem dos fornos:
#mdbca<-aov(vlr~vel+frn,data=refino)
#summary(mdbca)
# CV do Pacote Agricolae
#cv.model(mdbca)
# Configura tela para todos os gráficos (caso contrário, vai de um em um):
#par(mfrow = c(2,2))
# Analise Grafica da Normalidade e Homogeneidade das Variâncias
#plot(mdbca,which=1) # Resíduos x Valores previstos
#plot(mdbca,which=2) # gráfico qqplot
#### Teste de Shapiro-Wilks
#shapiro.test(mdbca$residuals)
## Qual solucao deveria ser recomendada? Use o Teste de Tukey.
### Te s t e de Tukey Ag r i c o l a e
#HSD.test(mdbca,"vel",group=TRUE,console=TRUE)
Respostas (Questão 02):
Letra (a):
# Analise Grafica dos Dados
# Boxplot : yx velocidades
ggplot(refino,aes(x=vel,y=vlr,fill=vel))+geom_boxplot()
# Boxplot : yx fornos
ggplot(refino,aes(x=frn,y=vlr,fill=frn))+geom_boxplot()
# Configura tela para todos os gráficos (caso contrário, vai de um em um):
#par(mfrow = c(2,2))
# Verificando Interacao entre Blocos e Niveis do Fator (velocidade nos fornos x graos)
interaction.plot(refino$vel,refino$frn,refino$vlr,fixed=TRUE)
interaction.plot(refino$frn,refino$vel,refino$vlr,fixed=TRUE)
# Teste Tukey de Aditividade
# Hipótese Ho: Não há efeito de interação
tukey.add.test(refino$vlr,refino$vel,refino$frn)
##
## Tukey's one df test for additivity
## F = 2.9804177 Denom df = 8 p-value = 0.1225519
Letra (b):
# Estimação do Modelo como se fosse ANAVA 1 Fator (não utilizado)
#modANAVA1F<-aov(vlr~vel,data=refino)
# Tabela ANAVA:
#summary(modANAVA1F)
#Estimação do modelo com blocagem dos fornos:
modBLC<-aov(vlr~vel+frn,data=refino)
summary(modBLC)
## Df Sum Sq Mean Sq F value Pr(>F)
## vel 3 22.19 7.40 0.853 0.4995
## frn 3 165.19 55.06 6.348 0.0133 *
## Residuals 9 78.06 8.67
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# CV do modelo:
cv.model(modBLC)
## [1] 38.31024
plot(modBLC,which=1,main="RESÍDUOS X VALORES PREVISTOS")
plot(modBLC,which=2,main="Gráfico de normalidade")
shapiro.test(modBLC$residuals)
##
## Shapiro-Wilk normality test
##
## data: modBLC$residuals
## W = 0.9301, p-value = 0.2447
# Teste de Bartlett (não pode ser usado neste modelo - não há repetição em cada bloco:
#bartlett.test(vlr~interaction(vel,frn),data=refino)
Letra (c):
### Teste de Tukey Agricolae
HSD.test(modBLC,"vel",group=TRUE,console=TRUE)
##
## Study: modBLC ~ "vel"
##
## HSD Test for vlr
##
## Mean Square Error: 8.673611
##
## vel, means
##
## vlr std r Min Max
## V05 5.75 1.707825 4 4 8
## V10 8.50 4.041452 4 5 14
## V15 7.75 5.057997 4 2 14
## V20 8.75 6.020797 4 3 17
##
## alpha: 0.05 ; Df Error: 9
## Critical Value of Studentized Range: 4.41489
##
## Honestly Significant Difference: 6.501145
##
## Means with the same letter are not significantly different.
##
## Groups, Treatments and means
## a V20 8.75
## a V10 8.5
## a V15 7.75
## a V05 5.75
————————————————————————————————————————————————————————————————————————-
Desenvolvimento/rascunho:
# Leitura do arquivo:
processos<-read.table("processos.txt",h=T)
#processos
#str(processos)
#summary(processos)
# Comando SUMMARY Via biblioteca doBy:
# Com médias:
#summaryBy(vlr~tmp,data=processos,FUN=c(length,mean))
# Com medianas:
#summaryBy(vlr~tmp,data=processos,FUN=c(length,median))
# Convertendo Temperatura e Pressão para FATORES:
#processos$temp <-as.factor(processos$temp)
#processos$pres<-as.factor(processos$pres)
#str(processos)
# Analise Grafica dos Dados:
#interaction.plot(processos$prs,processos$tmp,processos$vlr,type="b",pch=19,fixed=T,xlab="Pressão (psig)",ylab="Temperatura (ºC )" )
# As notações (duas formas) abaixo retornam o mesmo resultado:
# Estimacao do Modelo : Forma 1
#modelo.aov<-aov(vlr~temp+pres+temp:pres,data=processos)
# Estimacao do Modelo : Forma 2
#modelo.aov<-aov(vlr~temp*pres,data=processos)
# Tabela da ANAVA
#summary(modelo.aov)
# Verifica Normalidade e Variancias Homogeneas
#plot(modelo.aov,which=1) # Residudos x Valores previstos
#plot(modelo.aov,which=2) # grafico qqnorm
# Os gráficos reforçam a evidência de que os residuos tem uma distribuição linear normal.
#### Teste de Shapiro-Wilks : H0 : os dados tem dist. aprox. normal
#shapiro.test(modelo.aov$residuals)
#### Teste de Fligner-Killeen : H0 : Variancia dos resíduos é homogenea
#fligner.test(vlr~interaction(temp,pres),data=processos)
# Graficos dos Efeitos R base
#plot.design(vlr~temp*pres,data=processos)
# Grafico dos Efeitos Pacote effects
#plot(allEffects(modelo.aov))
# Necessario para os testes com o agricolae
#inter<-with(processos,interaction(temp,pres))
#inter
#amod<-aov(vlr~inter,data=processos)
#amod
#HSD.test(amod,"inter",group=TRUE,console=TRUE)
# cálculo PWR:
# Cálculo do poder para o teste em questão com: k=3 -> niveis de fatores e não fatores (que são dois: temperatura e pressão) - n=18 valores
#pwr.anova.test(k=3,f=.25,sig.level=.05,n=18)
# Se quisermos um teste com força de 80%:
#pwr.anova.test(k=3,power=.80,f=.25,sig.level=.05)
# Teremos que usar um n=53 !
Respostas (Questão 03):
Letra (a):
# Analise Grafica dos Dados:
interaction.plot(processos$prs,processos$tmp,processos$vlr,type="b",pch=19,fixed=T,xlab="Pressão (psig)",ylab="Temperatura (ºC )" )
# Teste Tukey de Aditividade
# Hipótese Ho: Não há efeito de interação
tukey.add.test(processos$vlr,processos$tmp,processos$prs)
##
## Tukey's one df test for additivity
## F = 0.0055883 Denom df = 12 p-value = 0.9416415
O teste Tukey de Aditividade rejeita a hipótese de não haver efeito de interação com p-value em 94% (vide acima).
Análise estatística via Estimação do modelo:
# Estimacao do Modelo (forma 1):
modEST<-aov(vlr~tmp+prs+tmp:prs,data=processos)
# Tabela da ANAVA
summary(modEST)
## Df Sum Sq Mean Sq F value Pr(>F)
## tmp 2 0.3011 0.1506 8.469 0.008539 **
## prs 2 0.7678 0.3839 21.594 0.000367 ***
## tmp:prs 4 0.0689 0.0172 0.969 0.470006
## Residuals 9 0.1600 0.0178
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Letra (b):
# Verifica Normalidade e Variancias Homogeneas
plot(modEST,which=1,main="Resíduos x Valores Previstos") # Residudos x Valores previstos
plot(modEST,which=2,main="Normalidade") # grafico qqnorm
#Teste de Shapiro-Wilks : H0 : os dados tem dist. aprox. normal
shapiro.test(modEST$residuals)
##
## Shapiro-Wilk normality test
##
## data: modEST$residuals
## W = 0.87366, p-value = 0.02046
# Teste de Fligner-Killeen : H0 : Variancia dos resíduos é homogenea
fligner.test(vlr~interaction(tmp,prs),data=processos)
##
## Fligner-Killeen test of homogeneity of variances
##
## data: vlr by interaction(tmp, prs)
## Fligner-Killeen:med chi-squared = 14.869, df = 8, p-value =
## 0.06174
Letra (c):
#Graficos dos Efeitos R base
plot.design(vlr~tmp*prs,data=processos,ylim=c(90,91))
abline(h=max(processos$vlr),col="blue",lwd=1)
# Grafico dos Efeitos (Pacote effects)
plot(allEffects(modEST),ylim=c(90,91.2))
# Pacote agricolae:
inter<-with(processos,interaction(tmp,prs))
amod<-aov(vlr~inter,data=processos)
HSD.test(amod,"inter",group=TRUE,console=TRUE)
##
## Study: amod ~ "inter"
##
## HSD Test for vlr
##
## Mean Square Error: 0.01777778
##
## inter, means
##
## vlr std r Min Max
## T150.P200 90.30 0.14142136 2 90.2 90.4
## T150.P215 90.65 0.07071068 2 90.6 90.7
## T150.P230 90.30 0.14142136 2 90.2 90.4
## T160.P200 90.20 0.14142136 2 90.1 90.3
## T160.P215 90.55 0.07071068 2 90.5 90.6
## T160.P230 90.00 0.14142136 2 89.9 90.1
## T170.P200 90.60 0.14142136 2 90.5 90.7
## T170.P215 90.85 0.07071068 2 90.8 90.9
## T170.P230 90.25 0.21213203 2 90.1 90.4
##
## alpha: 0.05 ; Df Error: 9
## Critical Value of Studentized Range: 5.594712
##
## Honestly Significant Difference: 0.5274745
##
## Means with the same letter are not significantly different.
##
## Groups, Treatments and means
## a T170.P215 90.85
## ab T150.P215 90.65
## ab T170.P200 90.6
## ab T160.P215 90.55
## bc T150.P200 90.3
## bc T150.P230 90.3
## bc T170.P230 90.25
## bc T160.P200 90.2
## c T160.P230 90
————————————————————————————————————————————————————————————————————————-