————————————————————————————————————————————————————————————————————————-

————————————————————————————————————————————————————————————————————————-

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)

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

————————————————————————————————————————————————————————————————————————-