setwd("C:\\Users\\Pessoal\\Desktop\\ESTATÍSTICA\\UFPB\\8º PERÍODO\\PLANEJAMENTO DE EXPERIMENTOS I\\BOOK")Planejamento de Experimentos
O objetivo é fazer um exemplo de cada tipo de experimento descrito no sumário deste trabalho, tipos de experimentos estudados na disciplina de Planejamentos de Experimentos I, no período 2022.2 na Universidade Federal da Paraíba. Disciplina que foi ministrada pelo docente Dr. João Agnaldo do Nascimento.
Aluno: Paulo Manoel da Silva Junior
Matrícula: 20190041314
1 Anova - One Way
Exemplo retirado do livro do Douglas C. Montgomery, que tem por título Design and Analysis of Experiments, a versão que foi utilizada foi a oitava versão.
- Segue a descrição do problema utilizado para a aplicação da Anova one way.
– Um desenvolvedor de produto está investigando a resistência à tração de uma nova fibra sintética que será usada para fazer tecidos para camisas masculinas. A força é geralmente afetada pela porcentagem de algodão utilizado na mistura de materiais para a fibra. O engenheiro conduz um experimento completamente aleatório com cinco níveis de conteúdo de algodão e replica o experimento cinco vezes. Os dados são mostrados na tabela a seguir.
Esse modelo se encaixa em um experimento de efeito fixo, pois as porcentagens são discretas.
| Cotton Weight Percent | Observations | ||||
|---|---|---|---|---|---|
| 15 | 7 | 7 | 15 | 11 | 9 |
| 20 | 12 | 17 | 12 | 18 | 18 |
| 25 | 14 | 19 | 19 | 18 | 18 |
| 30 | 19 | 25 | 22 | 19 | 23 |
| 35 | 7 | 10 | 11 | 15 | 11 |
Definindo o diretório e carregando o banco de dados
dados <- read.csv2("dados.csv", header = T, col.names = c("Percentual", "Observação"))Carregando a biblioteca dplyr para visualização dos dados
dados %>%
gt::gt()| Percentual | Observação |
|---|---|
| 15 | 7 |
| 15 | 7 |
| 15 | 15 |
| 15 | 11 |
| 15 | 9 |
| 20 | 12 |
| 20 | 17 |
| 20 | 12 |
| 20 | 18 |
| 20 | 19 |
| 25 | 14 |
| 25 | 19 |
| 25 | 19 |
| 25 | 18 |
| 25 | 18 |
| 30 | 19 |
| 30 | 25 |
| 30 | 22 |
| 30 | 19 |
| 30 | 23 |
| 35 | 7 |
| 35 | 10 |
| 35 | 11 |
| 35 | 15 |
| 35 | 11 |
1.1 Estatística Descritiva do banco
DescTools::Desc(dados$Observação, plotit = F, digits = 4)------------------------------------------------------------------------------
dados$Observação (integer)
length n NAs unique 0s mean meanCI'
25 25 0 13 0 15.0800 12.9420
100.0% 0.0% 0.0% 17.2180
.05 .10 .25 median .75 .90 .95
7.0000 7.8000 11.0000 15.0000 19.0000 20.8000 22.8000
range sd vcoef mad IQR skew kurt
18.0000 5.1794 0.3435 5.9304 8.0000 -0.0042 -1.1572
lowest : 7 (3), 9, 10, 11 (3), 12 (2)
highest: 18 (3), 19 (5), 22, 23, 25
heap(?): remarkable frequency (20.0%) for the mode(s) (= 19)
' 95%-CI (classic)Visualização gráfica
ggplot2::ggplot(dados, ggplot2::aes(y = Observação, x = factor(Percentual), fill = factor(Percentual)), legend.tile = "Percentual") + ggplot2::scale_fill_brewer(palette="Dark2") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F, alpha = 0.8)+
ggplot2::ggtitle("Boxplot da força de acordo com o percentual de algodão")+
ggplot2::xlab("Percentual") + ggplot2::ylab("Observação")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Percentual")
- De acordo com a visualização gráfica podemos observar que existe indícios de que existe uma diferença na resistência à tração de acordo com o percentual de algodão. Existindo evidências de que o melhor percentual é 30.
tapply(dados$Observação, dados$Percentual, mean) 15 20 25 30 35
9.8 15.6 17.6 21.6 10.8 Conforme podemos ver em média a resistência do percentual 30 é de 21.6.
1.2 Aplicação do modelo One Way
- O objetivo é verificar se há diferença significativa entre os grupos, ou seja, se existe diferença significativa devido aos tipos diferentes de percentual. E a função também ajusta um modelo de regressão para possíveis novos valores de resistência com base em outros percentuais de algodão.
attach(dados)
ExpDes::crd(Percentual, Observação, quali = F,
mcomp = "tukey", nl = FALSE, hvar = "levene",
sigT = 0.05, sigF = 0.05)------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Treatament 4 476.64 119.16 14.254 1.1722e-05
Residuals 20 167.20 8.36
Total 24 643.84
------------------------------------------------------------------------
CV = 19.17 %
------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value: 0.1444829
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
------------------------------------------------------------------------
Homogeneity of variances test
p-value: 0.5749962
According to the test of levene at 5% of significance, residuals can be considered homocedastic.
------------------------------------------------------------------------
Adjustment of polynomial models of regression
------------------------------------------------------------------------
Linear Model
=========================================
Estimate Standard.Error tc p.value
-----------------------------------------
b0 11.0800 2.1247 5.2148 0.00004
b1 0.1600 0.0818 1.9565 0.0645
-----------------------------------------
R2 of linear model
--------
0.067137
--------
Analysis of Variance of linear model
================================================
DF SS MS Fc p.value
------------------------------------------------
Linear Effect 1 32 32 3.83 0.06452
Lack of fit 3 444.6400 148.2133 17.73 1e-05
Residuals 20 167.2000 8.3600
------------------------------------------------
------------------------------------------------------------------------
Quadratic Model
==========================================
Estimate Standard.Error tc p.value
------------------------------------------
b0 -40.1771 8.2275 -4.8833 0.0001
b1 4.6171 0.6960 6.6339 0
b2 -0.0891 0.0138 -6.4487 0
------------------------------------------
R2 of quadratic model
--------
0.796528
--------
Analysis of Variance of quadratic model
===================================================
DF SS MS Fc p.value
---------------------------------------------------
Linear Effect 1 32 32 3.83 0.06452
Quadratic Effect 1 347.6571 347.6571 41.59 0
Lack of fit 2 96.9829 48.4914 5.8 0.01031
Residuals 20 167.2000 8.3600
---------------------------------------------------
------------------------------------------------------------------------
Cubic Model
==========================================
Estimate Standard.Error tc p.value
------------------------------------------
b0 58.8229 37.7096 1.5599 0.1345
b1 -8.5095 4.9289 -1.7264 0.0997
b2 0.4609 0.2049 2.2490 0.0359
b3 -0.0073 0.0027 -2.6901 0.0141
------------------------------------------
R2 of cubic model
--------
0.923458
--------
Analysis of Variance of cubic model
===================================================
DF SS MS Fc p.value
---------------------------------------------------
Linear Effect 1 32 32 3.83 0.06452
Quadratic Effect 1 347.6571 347.6571 41.59 0
Cubic Effect 1 60.5000 60.5000 7.24 0.01408
Lack of fit 1 36.4829 36.4829 4.36 0.0497
Residuals 20 167.2000 8.3600
---------------------------------------------------
------------------------------------------------------------------------Resposta: Depois de procedido o teste, podemos observar que ao nível de confiança de 95% e com base na amostra rejeitamos a hipótese nula de que não há diferença significativa devido ao percentual diferente de algodão, concluímos com isso que o percentual de algodão interfere diretamente na resistência da camisa. O resultado do p-valor foi de 1.1722 \times 10^{-5}. O modelo de regressão ajustado foi melhor sendo o cúbico, tendo um R^2, coeficiente de determinação que quanto mais próximo de um 1 melhor, igual a 0.9235. Todavia, no modelo cúbico o Lack of fit não foi aceito, sendo assim não deveria ser utilizado o modelo de regressão nessa questão.
Outros resultados importantes foram o teste de normalidade, que ao nível de 5% de significância, sendo aplicado o teste de Shapiro-Wilk, foi verificado que os resíduos estão distribuidos normalmente. E o outro resultado importante ao nível de significância de 5%, e aplicado o teste de Levene foi que os resíduos são homocedastiscos.
1.3 Estimação do intervalo de confiança para os percentuais
knitr::kable(rcompanion::groupwiseMean(Observação~Percentual,data = dados, conf = 0.95, traditional = T,digits =4))| Percentual | n | Mean | Conf.level | Trad.lower | Trad.upper |
|---|---|---|---|---|---|
| 15 | 5 | 9.8 | 0.95 | 5.645 | 13.96 |
| 20 | 5 | 15.6 | 0.95 | 11.430 | 19.77 |
| 25 | 5 | 17.6 | 0.95 | 15.030 | 20.17 |
| 30 | 5 | 21.6 | 0.95 | 18.360 | 24.84 |
| 35 | 5 | 10.8 | 0.95 | 7.244 | 14.36 |
Comentário: Com a estimação do intervalo de confiança para a resistência da camisa de acordo com o percentual de algodão, podemos já de início observar que existe grupos que são formados, todavia isso será melhor observado graficamente e analiticamente através das comparações múltiplas.
1.4 Suposições do modelo
- Para verificar as suposições do modelo é preciso que padronizemos os resíduos, para que valores, para que não seja observado valores por conta da escala na qual os dados estão inseridos.
Padronizando os resíduos, é necessário principalmente para verificação dos outliers.
anova <- aov(Observação~factor(Percentual), data = dados)
res <- anova$residuals
SSE <- sum(res^2)
MSE <- SSE/anova$df.residual
resid.padronizado <- res/sqrt(MSE)
resid.padronizado 1 2 3 4 5 6
-0.96840025 -0.96840025 1.79845761 0.41502868 -0.27668579 -1.24508603
7 8 9 10 11 12
0.48420012 -1.24508603 0.83005736 1.17591459 -1.24508603 0.48420012
13 14 15 16 17 18
0.48420012 0.13834289 0.13834289 -0.89922880 1.17591459 0.13834289
19 20 21 22 23 24
-0.89922880 0.48420012 -1.31425748 -0.27668579 0.06917145 1.45260037
25
0.06917145 Temos acima os resíduos padronizados
1.4.1 Outliers
Utilizaremos a função dixon.test do pacote outliers, testamos se existe outlier nos resíduos. As hipóteses são as seguintes.
H_0: Não \hspace{0.1cm} existe \hspace{0.1cm} outliers
H_1: Existe \hspace{0.1cm} outliers
Teste utilizando os resíduos sem ser padronizados
outliers::dixon.test(res)
Dixon test for outliers
data: res
Q.3 = 0.20455, p-value = 0.9427
alternative hypothesis: highest value 5.2 is an outlierTeste utilizando os resíduos padronizados
outliers::dixon.test(resid.padronizado)
Dixon test for outliers
data: resid.padronizado
Q.3 = 0.20455, p-value = 0.9427
alternative hypothesis: highest value 1.79845760605179 is an outlierResultado: Não rejeitamos a hipótese nula, ou seja, com base na amostra e com 95% de confiança, não temos outliers nos resíduos, e com a verificação do teste, podemos obervar que tanto utilizando o resíduo padronizado, como o mesmo sem estar padronizado, obtivemos a mesma estatística de teste, que foi igual a 0.2045, e o p-valor foi igual a 0.9427.
1.4.2 Normalidade
- Testando a normalidade dos resíduos, utilizaremos a função
shapiro.test, que aplica o teste de Shapiro-Wilk, para verificar a normalidade dos dados, temos então as seguintes hipóteses:
H_0: Os \hspace{0.1cm} resíduos \hspace{0.1cm} estão \hspace{0.1cm} distribuidos \hspace{0.1cm}normalmente H_1: Os \hspace{0.1cm} resíduos \hspace{0.1cm} não \hspace{0.1cm} estão \hspace{0.1cm} distribuidos \hspace{0.1cm}normalmente
Aplicando o teste nos resíduos padronizados
shapiro.test(resid.padronizado)
Shapiro-Wilk normality test
data: resid.padronizado
W = 0.93954, p-value = 0.1445Resposta: Não rejeitamos a hipótese nula, ou seja, com 95% de confiança e com base na amostra, concluímos que os resíduos estão distribuídos normalmente, o resultado do p-valor foi de 0.1445.
1.4.3 Homogeneidade de Variância
- Vamos começar utilizando o teste de Levene, que é utilizado quando não se tem normalidade nos dados, que não é o caso no qual está sendo considerado. O teste de Levene, encontra-se na biblioteca
lawstate na funçãolevene.test.
Teste de Levene
No teste de Levene, as hipóteses são as seguintes:
H_0: Possui \hspace{0.1cm} Homocedasticidade H_1: Não \hspace{0.1cm} Possui \hspace{0.1cm} Homocedasticidade
lawstat::levene.test(Observação, Percentual)
Modified robust Brown-Forsythe Levene-type test based on the absolute
deviations from the median
data: Observação
Test Statistic = 0.39474, p-value = 0.81Resposta: Não rejeitamos a hipótese nula, ou seja, com 95% de confiança e com base na amostra, não temos problema com falta de homocedasticidade, o resultado do p-valor foi de 0.81.
Teste de Fligner
No teste de Fligner, o teste encontra-se na biblioteca stats, função fligner.test as hipóteses são as seguintes:
H_0: Possui \hspace{0.1cm} Homocedasticidade H_1: Não \hspace{0.1cm} Possui \hspace{0.1cm} Homocedasticidade
fligner.test(Observação, Percentual)
Fligner-Killeen test of homogeneity of variances
data: Observação and Percentual
Fligner-Killeen:med chi-squared = 1.963, df = 4, p-value = 0.7426Resposta: Não rejeitamos a hipótese nula, ou seja, com 95% de confiança e com base na amostra, não temos problema com falta de homocedasticidade, o resultado do p-valor foi de 0.7426.
Teste de Bartlett
No teste de Bartlett, o teste encontra-se na biblioteca stats, função bartlett.test as hipóteses são as seguintes:
H_0: Possui \hspace{0.1cm} Homocedasticidade H_1: Não \hspace{0.1cm} Possui \hspace{0.1cm} Homocedasticidade
bartlett.test(Observação, Percentual)
Bartlett test of homogeneity of variances
data: Observação and Percentual
Bartlett's K-squared = 1.0797, df = 4, p-value = 0.8975Resposta: Não rejeitamos a hipótese nula, ou seja, com 95% de confiança e com base na amostra, não temos problema com falta de homocedasticidade, o resultado do p-valor foi de 0.8975.
Foi verificado que não temos problemas com a heterocedasticidade da variância, o que foi comprovado através dos três testes acima.
1.5 Anova Robusta
- Supondo a violação da homocedasticidade, precisamos utilizar a anova Robusta OneWay para resultados mais confiáveis, sendo assim vamos supor que foi violada essa importante característica.
1.5.1 Teste de Welch
No teste de Welch, temos as seguintes hipóteses:
H_0: \tau_{1} = \tau_2 = \tau_3 = .... \tau_i = 0 H_1: \tau_{i} \neq \tau_{j} \hspace{0.1cm}, \quad para \hspace{0.1cm} pelo \hspace{0.1cm} menos \hspace{0.1cm} um \hspace{0.1cm} i \hspace{0.1cm} \neq j
O teste de Welch, se encontra no pacote onewaytests, com a função welch.test.
onewaytests::welch.test(Observação~Percentual, data = dados)
Welch's Heteroscedastic F Test (alpha = 0.05)
-------------------------------------------------------------
data : Observação and Percentual
statistic : 12.42802
num df : 4
denom df : 9.905549
p.value : 0.0007070279
Result : Difference is statistically significant.
------------------------------------------------------------- Resultado: Rejeitamos a hipótese nula, ou seja, com 95% de confiança e com base na amostra, supondo a violação da homocedasticidade da variância, concluímos que existe diferença significativa devido ao percentual de algodão na resistência a tração da camisa.
- Ainda supondo a violação da homocedasticidade das variâncias, aplicaremos o método de comparações múltiplas de Welch.
knitr::kable(onewaytests::paircomp(onewaytests::welch.test(Observação~Percentual, data = dados)))
Welch's Heteroscedastic F Test (alpha = 0.05)
-------------------------------------------------------------
data : Observação and Percentual
statistic : 12.42802
num df : 4
denom df : 9.905549
p.value : 0.0007070279
Result : Difference is statistically significant.
-------------------------------------------------------------
Bonferroni Correction (alpha = 0.05)
-----------------------------------------------------
Level (a) Level (b) p.value No difference
1 15 20 0.256810809 Not reject
2 15 25 0.034148385 Reject
3 15 30 0.003223292 Reject
4 15 35 1.000000000 Not reject
5 20 25 1.000000000 Not reject
6 20 30 0.146065454 Not reject
7 20 35 0.419052058 Not reject
8 25 30 0.289912826 Not reject
9 25 35 0.032427474 Reject
10 30 35 0.002587586 Reject
----------------------------------------------------- | Level (a) | Level (b) | p.value | No difference |
|---|---|---|---|
| 15 | 20 | 0.2568108 | Not reject |
| 15 | 25 | 0.0341484 | Reject |
| 15 | 30 | 0.0032233 | Reject |
| 15 | 35 | 1.0000000 | Not reject |
| 20 | 25 | 1.0000000 | Not reject |
| 20 | 30 | 0.1460655 | Not reject |
| 20 | 35 | 0.4190521 | Not reject |
| 25 | 30 | 0.2899128 | Not reject |
| 25 | 35 | 0.0324275 | Reject |
| 30 | 35 | 0.0025876 | Reject |
1.5.2 Teste de Bryan Forsythe
No teste de Bryan Forsythe, temos as seguintes hipóteses:
H_0: \tau_{1} = \tau_2 = \tau_3 = .... \tau_i = 0 H_1: \tau_{i} \neq \tau_{j} \hspace{0.1cm}, \quad para \hspace{0.1cm} pelo \hspace{0.1cm} menos \hspace{0.1cm} um \hspace{0.1cm} i \hspace{0.1cm} \neq j
O teste de Bryan Forsythe, se encontra no pacote onewaytests, com a função bf.test.
onewaytests::bf.test(Observação~factor(Percentual), data = dados)
Brown-Forsythe Test (alpha = 0.05)
-------------------------------------------------------------
data : Observação and factor(Percentual)
statistic : 14.25359
num df : 4
denom df : 18.14843
p.value : 1.978908e-05
Result : Difference is statistically significant.
------------------------------------------------------------- Resultado: Rejeitamos a hipótese nula, ou seja, com 95% de confiança e com base na amostra, supondo a violação da homocedasticidade da variância, concluímos que existe diferença significativa devido ao percentual de algodão na resistência a tração da camisa.
1.6 Comparações Múltiplas
Aplicaremos alguns métodos de comparação múltiplas, dois a dois, para saber quais são os iguais e quais os diferentes, um dos métodos mais utilizados é o método de Tukey.
1.6.1 Teste de Tukey
TukeyHSD(anova) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Observação ~ factor(Percentual), data = dados)
$`factor(Percentual)`
diff lwr upr p adj
20-15 5.8 0.3279622 11.2720378 0.0344652
25-15 7.8 2.3279622 13.2720378 0.0030995
30-15 11.8 6.3279622 17.2720378 0.0000244
35-15 1.0 -4.4720378 6.4720378 0.9810840
25-20 2.0 -3.4720378 7.4720378 0.8076796
30-20 6.0 0.5279622 11.4720378 0.0273435
35-20 -4.8 -10.2720378 0.6720378 0.1033135
30-25 4.0 -1.4720378 9.4720378 0.2246058
35-25 -6.8 -12.2720378 -1.3279622 0.0105536
35-30 -10.8 -16.2720378 -5.3279622 0.0000791Resposta: De acordo com o teste de comparações múltiplas de Tukey, observamos que os percentuais 15 e 35 formam um grupo entre si, 25 e 20 formam outro grupo entre si, 35 e 20 e 30 e 25 formam outro grupo entre si, respectivamente.
Verificando com base em uma análise gráfica
plot(agricolae::HSD.test(anova, "factor(Percentual)"), main = "Teste de Tukey - Comparações Múltiplas", ylab = "Observações", xlab = "Percentuais")
De acordo com a visualização gráfica é reforçado o que foi analisado acima.
1.6.2 Teste de Fisher
- O teste de Fisher encontra-se na biblioteca
agricolae, e pode ser utilizado através da funçãoLSD.test.
out1 <- agricolae::LSD.test(anova, "factor(Percentual)", p.adj = "bonferroni")
out1$statistics
MSerror Df Mean CV t.value MSD
8.36 20 15.08 19.17352 3.153401 5.7665
$parameters
test p.ajusted name.t ntr alpha
Fisher-LSD bonferroni factor(Percentual) 5 0.05
$means
Observação std r LCL UCL Min Max Q25 Q50 Q75
15 9.8 3.346640 5 7.102727 12.49727 7 15 7 9 11
20 15.6 3.361547 5 12.902727 18.29727 12 19 12 17 18
25 17.6 2.073644 5 14.902727 20.29727 14 19 18 18 19
30 21.6 2.607681 5 18.902727 24.29727 19 25 19 22 23
35 10.8 2.863564 5 8.102727 13.49727 7 15 10 11 11
$comparison
NULL
$groups
Observação groups
30 21.6 a
25 17.6 ab
20 15.6 bc
35 10.8 cd
15 9.8 d
attr(,"class")
[1] "group"Resposta: O teste de Fisher, encontrou os mesmos grupos que o método de Tukey de comparações múltiplas.
Visualização gráfica
plot(out1, main = "Teste de Fisher - Comparações Múltiplas", ylab = "Observações", xlab = "Percentuais")
1.6.3 Teste de Scheffé
- O teste de Scheffé, pode ser encontrado na biblioteca
agricolae, através da funçãoscheffe.test.
out2 <- agricolae::scheffe.test(anova, "factor(Percentual)")
out2$statistics
MSerror Df F Mean CV Scheffe CriticalDifference
8.36 20 2.866081 15.08 19.17352 3.385901 6.191664
$parameters
test name.t ntr alpha
Scheffe factor(Percentual) 5 0.05
$means
Observação std r Min Max Q25 Q50 Q75
15 9.8 3.346640 5 7 15 7 9 11
20 15.6 3.361547 5 12 19 12 17 18
25 17.6 2.073644 5 14 19 18 18 19
30 21.6 2.607681 5 19 25 19 22 23
35 10.8 2.863564 5 7 15 10 11 11
$comparison
NULL
$groups
Observação groups
30 21.6 a
25 17.6 a
20 15.6 ab
35 10.8 b
15 9.8 b
attr(,"class")
[1] "group"Resposta: O teste de Scheffé, encontrou apenas dois grupos, sendo o grupo de percentuais 30,25 e 20, junto com o outro grupo 20,35 e 15, ou seja, o percentual 20 pertence aos dois grupos.
Visualização gráfica
plot(out2, main = "Teste de Scheffé - Comparações Múltiplas", xlab = "Percentuais", ylab= "Observações")
Para os intervalos de confiança temos:
DescTools::ScheffeTest(anova)$`factor(Percentual)` %>%
knitr::kable()| diff | lwr.ci | upr.ci | pval | |
|---|---|---|---|---|
| 20-15 | 5.8 | -0.3916641 | 11.9916641 | 0.0738903 |
| 25-15 | 7.8 | 1.6083359 | 13.9916641 | 0.0089344 |
| 30-15 | 11.8 | 5.6083359 | 17.9916641 | 0.0001003 |
| 35-15 | 1.0 | -5.1916641 | 7.1916641 | 0.9890732 |
| 25-20 | 2.0 | -4.1916641 | 8.1916641 | 0.8750947 |
| 30-20 | 6.0 | -0.1916641 | 12.1916641 | 0.0606444 |
| 35-20 | -4.8 | -10.9916641 | 1.3916641 | 0.1845557 |
| 30-25 | 4.0 | -2.1916641 | 10.1916641 | 0.3430893 |
| 35-25 | -6.8 | -12.9916641 | -0.6083359 | 0.0265515 |
| 35-30 | -10.8 | -16.9916641 | -4.6083359 | 0.0003039 |
1.6.4 Teste de Duncan
- Para o teste de Duncan, vamos utilizar a função
duncan.testque se encontra no pacoteagricolae
out3 <- agricolae::duncan.test(anova, "factor(Percentual)")
out3$statistics
MSerror Df Mean CV
8.36 20 15.08 19.17352
$parameters
test name.t ntr alpha
Duncan factor(Percentual) 5 0.05
$duncan
Table CriticalRange
2 2.949998 3.814519
3 3.096506 4.003964
4 3.189616 4.124360
5 3.254648 4.208450
$means
Observação std r Min Max Q25 Q50 Q75
15 9.8 3.346640 5 7 15 7 9 11
20 15.6 3.361547 5 12 19 12 17 18
25 17.6 2.073644 5 14 19 18 18 19
30 21.6 2.607681 5 19 25 19 22 23
35 10.8 2.863564 5 7 15 10 11 11
$comparison
NULL
$groups
Observação groups
30 21.6 a
25 17.6 b
20 15.6 b
35 10.8 c
15 9.8 c
attr(,"class")
[1] "group"Resultado: De acordo com o teste de Duncan, temos os seguintes grupos formados: 25 e 20 formando um grupo e 35 e 15 formando outro grupo. E o percentual 30 ficando como um único grupo.
Visualização gráfica do teste de Duncan
plot(out3, main = "Teste de Duncan - Comparações Múltiplas", xlab = "Percentuais", ylab = "Observações")
1.6.5 Teste SNK
- O teste SNK pode ser utilizado com a função
SNK.testdo pacoteagricolae, como segue abaixo.
out4 <- agricolae::SNK.test(anova, "factor(Percentual)")
out4$statistics
MSerror Df Mean CV
8.36 20 15.08 19.17352
$parameters
test name.t ntr alpha
SNK factor(Percentual) 5 0.05
$snk
Table CriticalRange
2 2.949998 3.814519
3 3.577935 4.626478
4 3.958293 5.118305
5 4.231857 5.472038
$means
Observação std r Min Max Q25 Q50 Q75
15 9.8 3.346640 5 7 15 7 9 11
20 15.6 3.361547 5 12 19 12 17 18
25 17.6 2.073644 5 14 19 18 18 19
30 21.6 2.607681 5 19 25 19 22 23
35 10.8 2.863564 5 7 15 10 11 11
$comparison
NULL
$groups
Observação groups
30 21.6 a
25 17.6 b
20 15.6 b
35 10.8 c
15 9.8 c
attr(,"class")
[1] "group"Resultado: De acordo com o teste SNK, temos os seguintes grupos formados: 25 e 20 formando um grupo e 35 e 15 formando outro grupo. E o percentual 30 ficando como um único grupo.
Visualização gráfica do teste SNK
plot(out4, main = "Teste SNK - Comparações Múltiplas", xlab = "Percentuais", ylab = "Observações")
1.6.6 Teste de Waller
- O teste de Waller pode ser realizado através da função
waller.testdo pacoteagricolae, como segue abaixo:
out5 <- agricolae::waller.test(anova, "factor(Percentual)")
out5$statistics
Mean Df CV MSerror F.Value Waller CriticalDifference
15.08 20 19.17352 8.36 14.25359 1.949 3.56406
$parameters
test name.t ntr K
Waller-Duncan factor(Percentual) 5 100
$means
Observação std r Min Max Q25 Q50 Q75
15 9.8 3.346640 5 7 15 7 9 11
20 15.6 3.361547 5 12 19 12 17 18
25 17.6 2.073644 5 14 19 18 18 19
30 21.6 2.607681 5 19 25 19 22 23
35 10.8 2.863564 5 7 15 10 11 11
$comparison
NULL
$groups
Observação groups
30 21.6 a
25 17.6 b
20 15.6 b
35 10.8 c
15 9.8 c
attr(,"class")
[1] "group"Resultado: De acordo com o teste de Waller, temos os seguintes grupos formados: 25 e 20 formando um grupo e 35 e 15 formando outro grupo. E o percentual 30 ficando como um único grupo.
Visualização gráfica do teste de Waller
plot(out5, main = "Teste de Waller - Comparações Múltiplas", xlab = "Percentuais", ylab = "Observações")
Observação: Esses testes de comparações múltiplas foram realizado por causa da aceitação de homocedasticidade das variâncias, se isso não ocorresse deveríamos utilizar o método de comparação de Welch.
1.6.7 Teste de Dunnet (Utilizado quando tem grupo controle)
Esse teste geralmente se é utilizado quando temos um grupo de controle, e como estamos analisando a resistência a tração da fibra em que será utilizada para ser feita a camisa masculina, verificamos que quanto maior o valor de interesse melhor, já foi observado que o grupo que tem a maior média é o grupo do percentual 30, logo, este valor será o valor de referência utilizado como controle.
DescTools::DunnettTest(Observação, Percentual, control = "30",conf.level = 0.95)
Dunnett's test for comparing several treatments with a control :
95% family-wise confidence level
$`30`
diff lwr.ci upr.ci pval
15-30 -11.8 -16.650493 -6.949507 4.6e-06 ***
20-30 -6.0 -10.850493 -1.149507 0.0130 *
25-30 -4.0 -8.850493 0.850493 0.1248
35-30 -10.8 -15.650493 -5.949507 4.0e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1É observado que como o objetivo é ter uma alta tração a resistência, não existe evidências estatísticas a um nível de significância de 5% para se dizer que os percentuais 30 e 25 são diferentes entre si, sendo preferível o uso destes, para que seja obtido um material de alta qualidade.
1.7 Contrastes ortogonais
Como no exemplo que está sendo considerado temos 5 grupos, temos 4 contrastes ortogonais, que serão verificados, que são eles.
\begin{cases} & C_1: \mu_{15} = \mu_{20} \\ & C_2: \mu_{25} = \mu_{30} \\ & C_3: \mu_{15} + \mu_{20} = \mu_{25} + \mu_{30} \\ & C_4: \mu_{15} + \mu_{20} + \mu_{25} + \mu_{30} = \mu_{35} \end{cases}
Criando os contrastes, depois prosseguindo para saber se são ortogonais.
C1 <- c(1,-1,0,0,0)
C2 <- c(0,0,1,-1,0)
C3 <- c(1,1,-1,-1,0)
C4 <- c(-1,-1,-1,-1,4)ibd::check.orthogonality(rbind(C1,C2,C3,C4))[1] 1Como o resultado do teste, foi observado que os contrastes são ortogonais.
Prosseguindo com a análise dos contrastes ortogonais
Percentual <-factor(Percentual)
c2 = rbind("C1: 15 = 20" = C1,
"C2: 25 = 30" = C2,
"C3: 15 + 20 = 25 + 30" = C3,
"C4: 15 + 20 + 25 + 30 = 35 " = C4)
contraste2 = gmodels::make.contrasts(c2)Registered S3 method overwritten by 'gdata':
method from
reorder.factor DescToolsAnovac2 = aov(Observação ~ Percentual,
contrasts = list(Percentual = contraste2))
summary(Anovac2, split = list(Percentual=list("C1: 15 = 20"=1, "C2: 25 = 30"=2, "C3: 15 + 20 = 25 + 30" =3, "C4: 15 + 20 + 25 + 30 = 35" = 4))) Df Sum Sq Mean Sq F value Pr(>F)
Percentual 4 476.6 119.16 14.254 1.17e-05 ***
Percentual: C1: 15 = 20 1 84.1 84.10 10.060 0.00480 **
Percentual: C2: 25 = 30 1 40.0 40.00 4.785 0.04076 *
Percentual: C3: 15 + 20 = 25 + 30 1 238.0 238.05 28.475 3.19e-05 ***
Percentual: C4: 15 + 20 + 25 + 30 = 35 1 114.5 114.49 13.695 0.00142 **
Residuals 20 167.2 8.36
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A um nível de significância de 5 % rejeitamos todos os contrastes, ou seja, nenhum dos teste de contrastes foi aceito.
2 Ancova
- A análise de covariância, pode ser considerado como uma extensão do modelo de Anova One-Way.
Exemplo retirado do livro do Douglas C. Montgomery, que tem por título Design and Analysis of Experiments, a versão que foi utilizada foi a oitava versão.
- Segue a descrição do problema utilizado para a aplicação da Ancova.
– Quatro formulações diferentes de uma cola industrial são sendo testado. A resistência à tração da cola quando aplicada unir peças também está relacionado com a espessura da aplicação. Cinco observações sobre resistência (y) em libras e espessura (x) em 0,01 polegada são obtidos para cada formulação. os dados são mostrado na tabela a seguir. Analise esses dados e desenhe conclusões apropriadas.
Glue Formulation
| 1 | 1 | 2 | 2 | 3 | 3 | 4 | 4 |
|---|---|---|---|---|---|---|---|
| y | x | y | x | y | x | y | x |
| 46.5 | 13 | 48.7 | 12 | 46.3 | 15 | 44.7 | 16 |
| 45.9 | 14 | 49.0 | 10 | 47.1 | 14 | 43.0 | 15 |
| 49.8 | 12 | 50.1 | 11 | 48.9 | 11 | 51.0 | 10 |
| 46.1 | 12 | 48.5 | 12 | 48.2 | 11 | 48.1 | 12 |
| 44.3 | 14 | 45.2 | 14 | 50.3 | 10 | 48.6 | 11 |
Tratamento <- c(rep(1,5),rep(2,5),rep(3,5),rep(4,5))
y <- c(46.5,45.9,49.8,46.1,44.3,48.7,49,50.1,48.5,45.2,46.3,47.1,48.9,48.2,50.3,44.7,43,51,48.1,48.6)
x <- c(13,14,12,12,14,12,10,11,12,14,15,14,11,11,10,16,15,10,12,11)
Tratamento <- as.factor(Tratamento)
banco <- data.frame(Tratamento, y,x)2.1 Visualização e Descritiva do banco
banco %>%
gt::gt()| Tratamento | y | x |
|---|---|---|
| 1 | 46.5 | 13 |
| 1 | 45.9 | 14 |
| 1 | 49.8 | 12 |
| 1 | 46.1 | 12 |
| 1 | 44.3 | 14 |
| 2 | 48.7 | 12 |
| 2 | 49.0 | 10 |
| 2 | 50.1 | 11 |
| 2 | 48.5 | 12 |
| 2 | 45.2 | 14 |
| 3 | 46.3 | 15 |
| 3 | 47.1 | 14 |
| 3 | 48.9 | 11 |
| 3 | 48.2 | 11 |
| 3 | 50.3 | 10 |
| 4 | 44.7 | 16 |
| 4 | 43.0 | 15 |
| 4 | 51.0 | 10 |
| 4 | 48.1 | 12 |
| 4 | 48.6 | 11 |
Estatística Descritiva do banco
DescTools::Desc(banco[2:3], plotit = F, digits = 3)------------------------------------------------------------------------------
Describe banco[2:3] (data.frame):
data frame: 20 obs. of 2 variables
20 complete cases (100.0%)
Nr ColName Class NAs Levels
1 y numeric .
2 x numeric .
------------------------------------------------------------------------------
1 - y (numeric)
length n NAs unique 0s mean meanCI'
20 20 0 = n 0 47.515 46.487
100.0% 0.0% 0.0% 48.543
.05 .10 .25 median .75 .90 .95
44.235 44.660 46.050 48.150 48.925 50.120 50.335
range sd vcoef mad IQR skew kurt
8.000 2.196 0.046 2.595 2.875 -0.314 -1.019
lowest : 43.000, 44.300, 44.700, 45.200, 45.900
highest: 49.000, 49.800, 50.100, 50.300, 51.000
' 95%-CI (classic)
------------------------------------------------------------------------------
2 - x (numeric)
length n NAs unique 0s mean meanCI'
20 20 0 7 0 12.450 11.598
100.0% 0.0% 0.0% 13.302
.05 .10 .25 median .75 .90 .95
10.000 10.000 11.000 12.000 14.000 15.000 15.050
range sd vcoef mad IQR skew kurt
6.000 1.820 0.146 2.224 3.000 0.300 -1.212
value freq perc cumfreq cumperc
1 10 3 15.0% 3 15.0%
2 11 4 20.0% 7 35.0%
3 12 5 25.0% 12 60.0%
4 13 1 5.0% 13 65.0%
5 14 4 20.0% 17 85.0%
6 15 2 10.0% 19 95.0%
7 16 1 5.0% 20 100.0%
' 95%-CI (classic)Visualização gráfica
ggplot2::ggplot(banco, ggplot2::aes(y = y, x = factor(Tratamento), fill = factor(Tratamento)), legend.tile = "Percentual") + ggplot2::scale_fill_brewer(palette="Dark2") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F)+
ggplot2::ggtitle("Boxplot da resistência a cola industrial quanto a quatro formulações diferentes")+
ggplot2::xlab("Formulações") + ggplot2::ylab("Resistência")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Formulações")
2.2 Adequação da covariável
- Neste caso testaremos a hipótese de que a covariável x é adequada no modelo OneWay, ou seja, ela é significativa no ajuste. As hipóteses são as seguintes:
H_0: \beta_i = 0
H_1: \beta_i \neq 0
ancova <- aov(y~factor(Tratamento)+x, data = banco)
summary(ancova) Df Sum Sq Mean Sq F value Pr(>F)
factor(Tratamento) 3 11.06 3.69 2.637 0.0876 .
x 1 59.57 59.57 42.624 9.54e-06 ***
Residuals 15 20.96 1.40
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Resposta: Rejeita H_0, ou seja, como podemos analisar a um nível de significância de 5% e com base na amostra, a covariável x é significativa para o auxílio da resistência à tração da cola industrial.
2.3 Significância das formulações
- Verificando se há efeito devido as quatro formulações diferentes, definindo as hipóteses:
H_0: \tau_1= \tau_2 = \tau_3 = \tau_4 = 0 H_1: \tau_i \neq \tau_j, \hspace{0.1cm} para \hspace{0.1cm} pelo \hspace{0.1cm} menos \hspace{0.1cm} um \hspace{0.1cm} i \hspace{0.1cm} \neq j
ancova2 <- aov(y~x+factor(Tratamento), data = banco)
summary(ancova2) Df Sum Sq Mean Sq F value Pr(>F)
x 1 68.85 68.85 49.269 4.15e-06 ***
factor(Tratamento) 3 1.77 0.59 0.422 0.74
Residuals 15 20.96 1.40
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Resposta: Não rejeitamos H_0, ou seja, com base na amostra e com 95% de confiança não existe efeito significativo devido ao tipo diferente de formulação da cola.
- Concluimos então a utilização da Ancova.
3 Bloco Aleatório
Experimento em bloco aleatório.
Exemplo retirado do livro do Douglas C. Montgomery, que tem por título Design and Analysis of Experiments, a versão que foi utilizada foi a oitava versão.
Segue a descrição do problema utilizado para a aplicação do modelo de blocos aleatórios.
Três soluções de lavagem diferentes estão sendo comparadas para estudar sua eficácia em retardar o crescimento de bactérias em Recipientes de leite de 5 galões. A análise é feita em laboratório, e apenas três testes podem ser executados em qualquer dia. Porque os dias poderiam representam uma fonte potencial de variabilidade, o experimentador decide usar um delineamento de blocos aleatórios. As observações são tomadas por quatro dias, e os dados são mostrados aqui.
Days
| Solution | 1 | 2 | 3 | 4 |
|---|---|---|---|---|
| 1 | 13 | 22 | 18 | 39 |
| 2 | 16 | 24 | 17 | 44 |
| 3 | 5 | 4 | 1 | 22 |
Segue então, que tomamos os dias como bloco e as soluções como o tratamento, segue abaixo a inserção dos dados.
dia <- c(rep(1:4,3))
solution <- c(rep(1,4),rep(2,4),rep(3,4))
y <- c(13,22,18,39,16,24,17,44,5,4,1,22)
dados <- data.frame(dia,solution,y)Visualizando o banco de dados
dados %>%
gt::gt()| dia | solution | y |
|---|---|---|
| 1 | 1 | 13 |
| 2 | 1 | 22 |
| 3 | 1 | 18 |
| 4 | 1 | 39 |
| 1 | 2 | 16 |
| 2 | 2 | 24 |
| 3 | 2 | 17 |
| 4 | 2 | 44 |
| 1 | 3 | 5 |
| 2 | 3 | 4 |
| 3 | 3 | 1 |
| 4 | 3 | 22 |
3.1 Estatística Descritiva do banco
- Média do crescimento de bactérias de acordo com a solução aplicada
tapply(dados$y, dados$solution, mean) 1 2 3
23.00 25.25 8.00 De acordo com o resultado acima, sem levar em consideração os dias como blocos, temos que existe uma diferença grande entre a solução 3 com as demais.
Visualização gráfica
ggplot2::ggplot(dados, ggplot2::aes(y = y, x = factor(solution), fill = factor(solution)), legend.tile = "Solução") + ggplot2::scale_fill_brewer(palette="Dark2") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F)+
ggplot2::labs(title = "Boxplot do crescimento das bactérias em galões de leite", subtitle = "Com relação a três tipos de lavagens diferentes")+
ggplot2::xlab("Solução") + ggplot2::ylab("Crescimento")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Soluções")
Média do crescimento de bactérias de acordo com os dias.
tapply(dados$y, dados$dia, mean) 1 2 3 4
11.33333 16.66667 12.00000 35.00000 Visualização gráfica
ggplot2::ggplot(dados, ggplot2::aes(y = y, x = factor(dia), fill = factor(dia)), legend.tile = "Dia") + ggplot2::scale_fill_brewer(palette="Dark2") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F)+
ggplot2::labs(title = "Boxplot do crescimento das bactérias em galões de leite", subtitle = "Com relação aos quatro dias diferentes")+
ggplot2::xlab("Dia") + ggplot2::ylab("Crescimento")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Dia")
Comentário: Podemos observar através da análise gráfica e da estimação pontual que foi feita acima, como vemos a média da resistência do crescimento, sem levar em consideração diferença de solução é observado que o valor do dia 1 foi 11.33.
DescTools::Desc(dados[3], plotit = F, digits = 3)------------------------------------------------------------------------------
Describe dados[3] (data.frame):
data frame: 12 obs. of 1 variables
12 complete cases (100.0%)
Nr ColName Class NAs Levels
1 y numeric .
------------------------------------------------------------------------------
1 - y (numeric)
length n NAs unique 0s mean meanCI'
12 12 0 11 0 18.750 10.483
100.0% 0.0% 0.0% 27.017
.05 .10 .25 median .75 .90 .95
2.650 4.100 11.000 17.500 22.500 37.500 41.250
range sd vcoef mad IQR skew kurt
43.000 13.011 0.694 8.154 11.500 0.492 -0.793
value freq perc cumfreq cumperc
1 1 1 8.3% 1 8.3%
2 4 1 8.3% 2 16.7%
3 5 1 8.3% 3 25.0%
4 13 1 8.3% 4 33.3%
5 16 1 8.3% 5 41.7%
6 17 1 8.3% 6 50.0%
7 18 1 8.3% 7 58.3%
8 22 2 16.7% 9 75.0%
9 24 1 8.3% 10 83.3%
10 39 1 8.3% 11 91.7%
11 44 1 8.3% 12 100.0%
' 95%-CI (classic)Visualização gráfica
ggplot2::ggplot(dados, ggplot2::aes(y))+
ggplot2::geom_histogram(breaks = hist(dados$y, plot = F)$breaks, fill = "white", color = "black")+
ggplot2::labs(title = "Histograma do crescimento de bactérias em galões de leite")+
ggplot2::xlab("Crescimento")+ggplot2::ylab("Frequência") + ggplot2::theme_test()
3.2 Verificação de Significância - Sem interação
- Verificando se existe diferença significativa nos tipos diferentes de soluções, alocando os dias como blocos para reduzir a variabilidade concernente aos dias como bloco.
ExpDes::rbd(solution, dia, y, quali = TRUE, mcomp = "tukey", nl = FALSE, sigT = 0.05, sigF = 0.05)------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Treatament 2 703.50 351.75 40.717 0.00032316
Block 3 1106.92 368.97 42.711 0.00019248
Residuals 6 51.83 8.64
Total 11 1862.25
------------------------------------------------------------------------
CV = 15.68 %
------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value: 0.4026655
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
------------------------------------------------------------------------
Homogeneity of variances test
p-value: 0.6965755
According to the test of oneillmathews at 5% of significance, the variances can be considered homocedastic.
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a 2 25.25
a 1 23
b 3 8
------------------------------------------------------------------------Resposta: Depois de procedido o teste, verificamos que existe uma diferença significativa com relação as soluções, ou seja, as soluções produzem um efeito diferente, ao nível de significância de 5%. É observado o que já existia indícios através da visualização gráfica separando a solução 1 e 2 como estatisticamente equivalentes, enquanto a solução 3 é estatisticamente diferente das demais.
Observando o resultado, vemos que é significativa a diferença entre os blocos, o que também já existia evidência na visualização gráfica, pois o dia 4 em média era diferente dos outros dias.
Outro resultado importante é que os resíduos estão distribuidos normalmente com 95% de confiança, verificado através do teste de Shapiro-Wilk, e que não tem problema de heterocedasticidade com a variância.
3.3 Verificação de Significância - Com interação
- Verificando se existe diferença significativa nos tipos diferentes de soluções, alocando os dias como blocos para reduzir a variabilidade concernente aos dias como bloco, e definindo interação dos dias com os tipos diferentes de soluções.
anova <- aov(y~dia+solution+dia*solution, data = dados)
summary(anova) Df Sum Sq Mean Sq F value Pr(>F)
dia 1 660.0 660.0 7.181 0.0279 *
solution 1 450.0 450.0 4.896 0.0578 .
dia:solution 1 16.9 16.9 0.184 0.6794
Residuals 8 735.3 91.9
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Resposta: Conforme o resultado acima a interação entre dia e o tipo diferente de solução não é estatisticamente significativa ao nível de confiança de 95%.
4 Quadrado Latino
Exemplo retirado do livro do Douglas C. Montgomery, que tem por título Design and Analysis of Experiments, a versão que foi utilizada foi a oitava versão.
Segue a descrição do problema utilizado para a aplicação do modelo de quadrado latino.
O efeito de cinco ingredientes diferentes (A, B, C, D, E) sobre o tempo de reação de um processo químico está sendo estudado. Cada lote de novo material é grande o suficiente para permitir cinco corridas a serem feitas. Além disso, cada execução requer aproximadamente 1 hora e meia, portanto, apenas cinco corridas podem ser feitas em um dia. O pesquisador decide executar o experimento como um quadrado latino para que os efeitos do dia e do lote podem ser sistematicamente controlados. Ele obtém os dados a seguir.
Day
| Batch | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|
| 1 | A = 8 | B = 7 | D = 1 | C = 7 | E = 3 |
| 2 | C = 11 | E = 2 | A = 7 | D = 3 | B = 8 |
| 3 | B = 4 | A = 9 | C = 10 | E = 1 | D = 5 |
| 4 | D = 6 | C = 8 | E = 6 | B = 6 | A = 10 |
| 5 | E = 4 | D = 2 | B = 3 | A = 8 | C = 8 |
Inserção dos dados
Lotes <- c(rep(1,5), rep(2,5), rep(3,5), rep(4,5),rep(5,5))
ingredientes <-c("A","B","D","C","E","C","E","A","D","B","B","A","C","E","D","D","C","E","B","A","E","D","B","A","C")
ingredientes <- as.factor(ingredientes)
experimento <- c(rep(1:5,5))
y <- c(8,7,1,7,3,11,2,7,3,8,4,9,10,1,5,6,8,6,6,10,4,2,3,8,8)
dados <- data.frame(experimento, Lotes, ingredientes, y)Visuzalizando o banco de dados
As primeiras 10 linhas do banco de dados.
dados %>%
slice_head(n=10) %>%
gt::gt()| experimento | Lotes | ingredientes | y |
|---|---|---|---|
| 1 | 1 | A | 8 |
| 2 | 1 | B | 7 |
| 3 | 1 | D | 1 |
| 4 | 1 | C | 7 |
| 5 | 1 | E | 3 |
| 1 | 2 | C | 11 |
| 2 | 2 | E | 2 |
| 3 | 2 | A | 7 |
| 4 | 2 | D | 3 |
| 5 | 2 | B | 8 |
4.1 Estatística Descritiva do banco
DescTools::Desc(dados[4], plotit = F)------------------------------------------------------------------------------
Describe dados[4] (data.frame):
data frame: 25 obs. of 1 variables
25 complete cases (100.0%)
Nr ColName Class NAs Levels
1 y numeric .
------------------------------------------------------------------------------
1 - y (numeric)
length n NAs unique 0s mean meanCI'
25 25 0 11 0 5.88 4.67
100.0% 0.0% 0.0% 7.09
.05 .10 .25 median .75 .90 .95
1.20 2.00 3.00 6.00 8.00 9.60 10.00
range sd vcoef mad IQR skew kurt
10.00 2.93 0.50 2.97 5.00 -0.12 -1.23
value freq perc cumfreq cumperc
1 1 2 8.0% 2 8.0%
2 2 2 8.0% 4 16.0%
3 3 3 12.0% 7 28.0%
4 4 2 8.0% 9 36.0%
5 5 1 4.0% 10 40.0%
6 6 3 12.0% 13 52.0%
7 7 3 12.0% 16 64.0%
8 8 5 20.0% 21 84.0%
9 9 1 4.0% 22 88.0%
10 10 2 8.0% 24 96.0%
11 11 1 4.0% 25 100.0%
' 95%-CI (classic)- Média do tempo de Reação de acordo com o tipo de ingrediente
tapply(dados$y, dados$ingredientes, mean) %>%
knitr::kable(col.names = c("Média"))| Média | |
|---|---|
| A | 8.4 |
| B | 5.6 |
| C | 8.8 |
| D | 3.4 |
| E | 3.2 |
- Variância do tempo de Reação de acordo com o tipo de ingrediente
tapply(dados$y, dados$ingredientes, var) %>%
knitr::kable(col.names = c("Variância"))| Variância | |
|---|---|
| A | 1.3 |
| B | 4.3 |
| C | 2.7 |
| D | 4.3 |
| E | 3.7 |
- Boxplot do tempo de reação químico com diferente tipos de ingredientes.
ggplot2::ggplot(dados, ggplot2::aes(y = y, x = factor(ingredientes), fill = factor(ingredientes)), legend.tile = "Ingredientes") + ggplot2::scale_fill_brewer(palette="Dark2") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F)+
ggplot2::labs(title = "Boxplot do tempo de reação químico", subtitle = "Com relação aos cinco diferentes ingredientes")+
ggplot2::xlab("Ingredientes") + ggplot2::ylab("Tempo de Reação")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Ingredientes")
4.2 Aplicação
ExpDes::latsd(ingredientes, Lotes, experimento, y, quali = TRUE, mcomp = "tukey", sigT = 0.05, sigF = 0.05, unfold = 1)------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Treatament 4 141.44 35.360 11.3092 0.00049
Row 4 15.44 3.860 1.2345 0.34762
Column 4 12.24 3.060 0.9787 0.45501
Residuals 12 37.52 3.127
Total 24 206.64
------------------------------------------------------------------------
CV = 30.07 %
------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value: 0.5476372
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a C 8.8
a A 8.4
ab B 5.6
b D 3.4
b E 3.2
------------------------------------------------------------------------Resposta: De acordo com o experimento conduzido, tem-se como conclusão que existe diferença significativa de acordo com o tipo diferente de Ingrediente, o teste de comparação múltipla de Tukey, já forneceu que estatisticamente os ingredientes C, A e B são iguais entre si no tempo de reação, e os ingredientes B, D e E também são iguais entre si.
Outro resultado importante é que ao nível de confiança de 95%, os resíduos estão distribuidos normalmente.
Foi obtido um coeficiente de variação elevado, sendo o mesmo observado de 30.07%.
5 Planejamento Fatorial 2x2
Exemplo retirado do livro do Douglas C. Montgomery, que tem por título Design and Analysis of Experiments, a versão que foi utilizada foi a oitava versão.
Segue a descrição do problema utilizado para a aplicação do modelo de quadrado latino.
Um artigo no AT&T Technical Journal (março/abril 1986, vol. 65, pp. 39–50) descreve a aplicação de dois níveis projetos fatoriais para fabricação de circuitos integrados. Uma etapa básica de processamento é crescer uma camada epitaxial em silício polido bolachas. Os wafers montados em um susceptor são posicionados dentro de uma redoma de vidro, e vapores químicos são introduzidos. O receptor é girado e o calor é aplicado até que a camada epitaxial seja grosso o suficiente. Um experimento foi executado usando dois fatores: arsênico vazão (A) e tempo de deposição (B). Quatro réplicas foram executadas, e a espessura da camada epitaxial foi medida (m). Os dados são mostrados na Tabela abaixo.
Replicate
| A | B | I | II | III | IV |
|---|---|---|---|---|---|
| - | - | 14.037 | 16.165 | 13.972 | 13.907 |
| + | - | 13.880 | 13.860 | 14.032 | 13.914 |
| - | + | 14.821 | 14.757 | 14.843 | 14.878 |
| + | + | 14.888 | 14.921 | 14.415 | 14.932 |
Inserindo os dados
A <- c(rep(0,4),rep(1,4),rep(0,4),rep(1,4))
B <- c(rep(0,4),rep(0,4),rep(1,4),rep(1,4))
replica <- c(rep(1:4,4))
y <- c(14.037, 16.165, 13.972, 13.907, 13.880, 13.860, 14.032, 13.914, 14.821, 14.757, 14.843, 14.878, 14.888, 14.921, 14.415, 14.932)
dados <- data.frame(A,B, replica, y)Visualizando as primeiras 10 linhas do banco de dados
dados %>%
slice_head(n=10) %>%
gt::gt()| A | B | replica | y |
|---|---|---|---|
| 0 | 0 | 1 | 14.037 |
| 0 | 0 | 2 | 16.165 |
| 0 | 0 | 3 | 13.972 |
| 0 | 0 | 4 | 13.907 |
| 1 | 0 | 1 | 13.880 |
| 1 | 0 | 2 | 13.860 |
| 1 | 0 | 3 | 14.032 |
| 1 | 0 | 4 | 13.914 |
| 0 | 1 | 1 | 14.821 |
| 0 | 1 | 2 | 14.757 |
5.1 Estatística descritiva do banco
- Média das observações de acordo com a presença ou ausência do fator A. Onde 0 = 55% e 1 = 59%
tapply(dados$y, factor(dados$A, levels = c(0,1), labels = c("55%","59%")), mean) %>%
knitr::kable(col.names = "Média")| Média | |
|---|---|
| 55% | 14.67250 |
| 59% | 14.35525 |
Visualização gráfica
ggplot2::ggplot(dados, ggplot2::aes(y = y, x = factor(A, levels = c(0,1), labels = c("55%", "59%")), fill = factor(A, levels = c(0,1), labels = c("55%", "59%"))), legend.tile = "Arsênico") + ggplot2::scale_fill_brewer(palette="Dark2") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F)+
ggplot2::labs(title = "Boxplot da espessura da camada epitaxial (m) ", subtitle = "Com relação a porcentagem de Arsênico")+
ggplot2::xlab("Arsênico") + ggplot2::ylab("Espessura da camada")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Arsênico")
- Média das observações de acordo com a presença ou ausência do fator B. Onde 0 = Tempo curto e 1 = Presença
tapply(dados$y, factor(dados$B, levels = c(0,1), labels = c("Tempo curto/10 min", "Tempo longo/15 min")), mean) %>%
knitr::kable(col.names = "Média")| Média | |
|---|---|
| Tempo curto/10 min | 14.22087 |
| Tempo longo/15 min | 14.80687 |
Visualização gráfica
ggplot2::ggplot(dados, ggplot2::aes(y = y, x = factor(B, levels = c(0,1), labels = c("Curto/10 min", "Longo/15 min")), fill = factor(B, levels = c(0,1), labels = c("Curto/10 min", "Longo/15 min"))), legend.tile = "Tempo") + ggplot2::scale_fill_brewer(palette="Dark2") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F)+
ggplot2::labs(title = "Boxplot da espessura da camada epitaxial (m) ", subtitle = "Com relação ao tempo de decomposição")+
ggplot2::xlab("Tempo de decomposição") + ggplot2::ylab("Espessura da camada")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Tempo")
5.2 Aplicação do Planejamento Fatorial 2x2
- Vamos utilizar a função
fat2.crddo pacoteExpDespara aplicação do método.
ExpDes::fat2.crd(A,B,y, quali = c(TRUE, TRUE), mcomp = "tukey", fac.names = c("Arsênico", "Tempo"), sigT = 0.05, sigF = 0.05)------------------------------------------------------------------------
Legend:
FACTOR 1: Arsênico
FACTOR 2: Tempo
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Arsênico 1 0.4026 4 1.2619 0.28327
Tempo 1 1.3736 5 4.3054 0.06016
Arsênico*Tempo 1 0.3170 2 0.9935 0.33856
Residuals 12 3.8285 3
Total 15 5.9216 1
------------------------------------------------------------------------
CV = 3.89 %
------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value: 0.0002643226
WARNING: at 5% of significance, residuals can not be considered normal!
------------------------------------------------------------------------
No significant interaction: analyzing the simple effect
------------------------------------------------------------------------
Arsênico
According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
Levels Means
1 0 14.67250
2 1 14.35525
------------------------------------------------------------------------
Tempo
According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
Levels Means
1 0 14.22087
2 1 14.80687
------------------------------------------------------------------------Resposta: De acordo com a aplicação do método de Planejamento Fatorial 2x2, é observado que não existe diferença estatísticamente significante a nenhum dos dois fatores tomados individualmente e também a interação deles, sendo assim temos que com a presença de cada um de maneira isolada, a presença de ambos e a ausência de ambos estatisticamente é igual a espessura da camada epitaxial.
Outro resultado que é ressaltado é que os resíduos não estão distribuidos normalmente com o nível de significância de 5%.
6 Planejamento Fatorial 2^3 ou 2x2x2
Exemplo retirado do livro do Douglas C. Montgomery, que tem por título Design and Analysis of Experiments, a versão que foi utilizada foi a oitava versão.
Segue a descrição do problema utilizado para a aplicação do modelo de planejamneto fatorial 2^3.
Uma empresa comercializa seus produtos por mala direta. Um experimento foi conduzido para estudar os efeitos de três fatores na taxa de resposta do cliente para um determinado produto. Os três fatores são: A um tipo de correio usado (3ª classe, 1ª classe), tipo B de brochura descritiva (cor, preto e branco) e C preço oferecido ($19.95, $24.95). As correspondências são feitas para dois grupos de 8.000 clientes selecionados aleatoriamente, com 1.000 clientes tomers em cada grupo recebendo cada combinação de tratamento. Cada grupo de clientes é considerado uma réplica. O variável de resposta é o número de pedidos feitos. A experiência os dados mentais são mostrados na Tabela abaixo.
Dados do Exemplo
| Coded Factors | Number of Orders | Factor Levels | ||||||
|---|---|---|---|---|---|---|---|---|
| Run | A | B | C | Replicate I | Replicate II | Low (-1) | High (+1) | |
| 1 | - | - | - | 50 | 54 | A (Class) | 3rd | 1st |
| 2 | + | - | - | 44 | 42 | B (Type) | Black-White | Color |
| 3 | - | + | - | 46 | 48 | C($) | $19.95 | $24.95 |
| 4 | + | + | - | 42 | 43 | |||
| 5 | - | - | + | 49 | 46 | |||
| 6 | + | - | + | 48 | 45 | |||
| 7 | - | + | + | 47 | 48 | |||
| 8 | + | + | + | 56 | 54 |
Inserindo os dados
A <- c(rep(-1,2),rep(1,2),rep(-1,2),rep(1,2),rep(-1,2),rep(1,2),rep(-1,2),rep(1,2))
B <- c(rep(-1,2),rep(-1,2),rep(1,2),rep(1,2),rep(-1,2),rep(-1,2),rep(1,2),rep(1,2))
C <- c(rep(-1,2),rep(-1,2),rep(-1,2),rep(-1,2),rep(1,2),rep(1,2),rep(1,2),rep(1,2))
replica <- c(rep(1:2,8))
y <- c(50,54,44,42,46,48,42,43,49,46,48,45,47,48,56,54)
dados <- data.frame(A,B,C,replica,y)Visualizando as primeiras 10 linhas do banco de dados
dados %>%
slice_head(n=10) %>%
gt::gt()| A | B | C | replica | y |
|---|---|---|---|---|
| -1 | -1 | -1 | 1 | 50 |
| -1 | -1 | -1 | 2 | 54 |
| 1 | -1 | -1 | 1 | 44 |
| 1 | -1 | -1 | 2 | 42 |
| -1 | 1 | -1 | 1 | 46 |
| -1 | 1 | -1 | 2 | 48 |
| 1 | 1 | -1 | 1 | 42 |
| 1 | 1 | -1 | 2 | 43 |
| -1 | -1 | 1 | 1 | 49 |
| -1 | -1 | 1 | 2 | 46 |
6.1 Estatística descritiva do banco
- Média das observações de acordo com o fator A (tipo de correio utilizado). Onde -1 = 3ª classe e 1 = 1ª classe.
tapply(dados$y, factor(dados$A, levels = c(-1,1), labels = c("3ª classe","1ª classe")), mean) %>%
knitr::kable(col.names = "Média")| Média | |
|---|---|
| 3ª classe | 48.50 |
| 1ª classe | 46.75 |
Visualização gráfica
ggplot2::ggplot(dados, ggplot2::aes(y = y, x = factor(A, levels = c(-1,1), labels = c("3ª classe", "1ª classe")), fill = factor(A, levels = c(-1,1), labels = c("3ª classe", "1ª classe"))), legend.tile = "Tipo de correio usado") + ggplot2::scale_fill_brewer(palette="Dark2") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F, alpha = 0.5)+
ggplot2::labs(title = "Boxplot da Quantidade de pedidos feitos ", subtitle = "Com relação ao tipo de correio utilizado")+
ggplot2::xlab("Tipo de correio utilizado") + ggplot2::ylab("Quantidade de pedidos feitos")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Tipo de correio")
- Média das observações de acordo com o fator B (tipo de brochura descritiva). Onde -1 = Preto e Branco e 1 = Colorida.
tapply(dados$y, factor(dados$B, levels = c(-1,1), labels = c("Preto e Branco","Colorida")), mean) %>%
knitr::kable(col.names = "Média")| Média | |
|---|---|
| Preto e Branco | 47.25 |
| Colorida | 48.00 |
Visualização gráfica
ggplot2::ggplot(dados, ggplot2::aes(y = y, x = factor(B, levels = c(-1,1), labels = c("Preto e Branco", "Colorida")), fill = factor(B, levels = c(-1,1), labels = c("Preto e Branco", "Colorida"))), legend.tile = "Tipo de brochura") + ggplot2::scale_fill_brewer(palette="Dark2") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F, alpha = 0.5)+
ggplot2::labs(title = "Boxplot da Quantidade de pedidos feitos ", subtitle = "Com relação ao tipo de brochura descritiva")+
ggplot2::xlab("Tipo de brochura descritiva") + ggplot2::ylab("Quantidade de pedidos feitos")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Tipo de brochura")
- Média das observações de acordo com o fator C (preço oferecido, como só tem dois valores poderia ser categorizado em mais barato e mais caro). Onde -1 = $ 19.95 e 1 = 24.95.
tapply(dados$y, factor(dados$C, levels = c(-1,1), labels = c("$19.95","$24.95")), mean) %>%
knitr::kable(col.names = "Média")| Média | |
|---|---|
| $19.95 | 46.125 |
| $24.95 | 49.125 |
Visualização gráfica
ggplot2::ggplot(dados, ggplot2::aes(y = y, x = factor(C, levels = c(-1,1), labels = c("$19.95", "$24.95")), fill = factor(C, levels = c(-1,1), labels = c("$19.95", "$24.95"))), legend.tile = "Preço oferecido") + ggplot2::scale_fill_brewer(palette="Dark2") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F, alpha = 0.5)+
ggplot2::labs(title = "Boxplot da Quantidade de pedidos feitos ", subtitle = "Com relação ao preço oferecido")+
ggplot2::xlab("Preço oferecido") + ggplot2::ylab("Quantidade de pedidos feitos")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Preço oferecido")
6.2 Aplicação do Planejamento Fatorial 2^3
- Vamos utilizar a função
fat3.rbddo pacoteExpDespara aplicação do método, pois em nosso exemplo tem a presença de dois blocos.
ExpDes::fat3.rbd(A,B,C, replica, y, quali = c(TRUE, TRUE, TRUE), mcomp = "tukey", fac.names = c("Tipo do correio", "Tipo de brochura", "Preço oferecido"), sigT = 0.05, sigF = 0.05)------------------------------------------------------------------------
Legend:
FACTOR 1: Tipo do correio
FACTOR 2: Tipo de brochura
FACTOR 3: Preço oferecido
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc
Block 1 0.25 0.25 0.0737
Tipo do correio 1 12.25 12.25 3.6105
Tipo de brochura 1 2.25 2.25 0.6632
Preço oferecido 1 36.00 36 10.6105
Tipo do correio*Tipo de brochura 1 42.25 42.25 12.4526
Tipo do correio*Preço oferecido 1 100.00 100 29.4737
Tipo de brochura*Preço oferecido 1 49.00 49 14.4421
Tipo do correio*Tipo de brochura*Preço oferecido 1 4.00 4 1.1789
Residuals 7 23.75 3.39286
Total 14 269.75
Pr>Fc
Block 0.7939
Tipo do correio 0.0992
Tipo de brochura 0.4423
Preço oferecido 0.0139
Tipo do correio*Tipo de brochura 0.0096
Tipo do correio*Preço oferecido 0.001
Tipo de brochura*Preço oferecido 0.0067
Tipo do correio*Tipo de brochura*Preço oferecido 0.3135
Residuals
Total
------------------------------------------------------------------------
CV = 3.87 %
------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value: 0.2862962
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
Significant Tipo do correio*Tipo de brochura interaction: analyzing the interaction
------------------------------------------------------------------------
Analyzing Tipo do correio inside of each level of Tipo de brochura
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Tipo do correio:Tipo de brochura -1 1 50.00 50.00000 14.7368 0.0064
Tipo do correio:Tipo de brochura 1 1 4.50 4.50000 1.3263 0.2873
Residuals 7 23.75 3.39286
------------------------------------------------------------------------
Tipo do correio inside of the level -1 of Tipo de brochura
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a -1 49.75
b 1 44.75
------------------------------------------------------------------------
Tipo do correio inside of the level 1 of Tipo de brochura
According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
Levels Means
1 -1 47.25
2 1 48.75
------------------------------------------------------------------------
Analyzing Tipo de brochura inside of each level of Tipo do correio
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Tipo de brochura:Tipo do correio -1 1 12.50 12.50000 3.6842 0.0964
Tipo de brochura:Tipo do correio 1 1 32.00 32.00000 9.4316 0.018
Residuals 7 23.75 3.39286
------------------------------------------------------------------------
Tipo de brochura inside of the level -1 of Tipo do correio
According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
Levels Means
1 -1 49.75
2 1 47.25
------------------------------------------------------------------------
Tipo de brochura inside of the level 1 of Tipo do correio
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a 1 48.75
b -1 44.75
------------------------------------------------------------------------
Significant Tipo do correio*Preço oferecido interaction: analyzing the interaction
------------------------------------------------------------------------
Analyzing Tipo do correio inside of each level of Preço oferecido
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Tipo do correio:Preço oferecido -1 1 91.125 91.12500 26.8579 0.0013
Tipo do correio:Preço oferecido 1 1 21.125 21.12500 6.2263 0.0413
Residuals 7 23.750 3.39286
------------------------------------------------------------------------
Tipo do correio inside of the level -1 of Preço oferecido
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a -1 49.5
b 1 42.75
------------------------------------------------------------------------
Tipo do correio inside of the level 1 of Preço oferecido
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a 1 50.75
b -1 47.5
------------------------------------------------------------------------
Analyzing Preço oferecido inside of each level of Tipo do correio
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Preço oferecido:Tipo do correio -1 1 8.00 8.00000 2.3579 0.1685
Preço oferecido:Tipo do correio 1 1 128.00 128.00000 37.7263 5e-04
Residuals 7 23.75 3.39286
------------------------------------------------------------------------
Preço oferecido inside of the level -1 of Tipo do correio
According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
Levels Means
1 -1 49.5
2 1 47.5
------------------------------------------------------------------------
Preço oferecido inside of the level 1 of Tipo do correio
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a 1 50.75
b -1 42.75
------------------------------------------------------------------------
Significant Tipo de brochura*Preço oferecido interaction: analyzing the interaction
------------------------------------------------------------------------
Analyzing Tipo de brochura inside of each level of Preço oferecido
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Tipo de brochura:Preço oferecido -1 1 15.125 15.12500 4.4579 0.0726
Tipo de brochura:Preço oferecido 1 1 36.125 36.12500 10.6474 0.0138
Residuals 7 23.750 3.39286
------------------------------------------------------------------------
Tipo de brochura inside of the level -1 of Preço oferecido
According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
Levels Means
1 -1 47.50
2 1 44.75
------------------------------------------------------------------------
Tipo de brochura inside of the level 1 of Preço oferecido
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a 1 51.25
b -1 47
------------------------------------------------------------------------
Analyzing Preço oferecido inside of each level of Tipo de brochura
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Preço oferecido:Tipo de brochura -1 1 0.50 0.50000 0.1474 0.7125
Preço oferecido:Tipo de brochura 1 1 84.50 84.50000 24.9053 0.0016
Residuals 7 23.75 3.39286
------------------------------------------------------------------------
Preço oferecido inside of the leve -1 of Tipo de brochura
According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
Levels Means
1 -1 47.5
2 1 47.0
------------------------------------------------------------------------
Preço oferecido inside of the leve 1 of Tipo de brochura
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a 1 51.25
b -1 44.75
------------------------------------------------------------------------Resposta: Conforme o procedimento do planejamento fatorial 2x2x2, temos como resultado que o fator C, sendo Preço oferecido, ao nível de confiança de 95% foi estatisticamente significante, ou seja, o preço oferecido afeta a quantidade de pedidos, tivemos também como resultado importante que as interações dois a dois foram significativas, ou seja analisando as interações dois a dois obtemos os seguintes resultados:
Podemos observar que em média o número de pedidos é maior dos individuos que preferem utilizar a 3ª classe de correio e também por um tipo de brochura preto e branco. E que em média o número de pedidos é maior dos individuos que preferem o tipo de brochura colorida e também utilizar a 1ª classe dos correios. Observa-se que em média o número de pedidos é maior dos individuos que escolhem a 1ª classe dos correios e juntamente com isso escolhem o preço oferecido de $24.95
Observa-se que o número médio de pedidos é maior dos individuos que selecionam o tipo de brochura colorida e o preço sugerido de $24.95.
Podemos observar ainda do resultado na análise residual que os reíduos estão distribuidos normalmente, ao nível de significância de 5%.
Observação: Os resultados de comparação das interações dois a dois devem ser tidos como referência com o grupo complementar do respectivo fator.
7 Planejamento Fatorial 3^3 ou 3x3x3
Exemplo retirado do livro do Douglas C. Montgomery, que tem por título Design and Analysis of Experiments, a versão que foi utilizada foi a oitava versão.
Segue a descrição do problema utilizado para a aplicação do modelo de planejamneto fatorial 2^3.
Um pesquisador médico está estudando o efeito de lidocaína no nível de enzima no músculo cardíaco de cães beagle. Três marcas comerciais diferentes de lidocaína (A), três níveis de dosagem (B), e três cães (C) são usados no experimento, e duas réplicas de um planejamento fatorial 3^3 são executadas.
Os dados estão dispostos na tabela abaixo:
Replicate I
| Dog | Dog | Dog | ||
|---|---|---|---|---|
| Lidocaine Brand | Dosage Strength | 1 | 2 | 3 |
| 1 | 1 | 96 | 84 | 85 |
| 1 | 2 | 94 | 99 | 98 |
| 1 | 3 | 101 | 106 | 98 |
| 2 | 1 | 85 | 84 | 86 |
| 2 | 2 | 95 | 98 | 97 |
| 2 | 3 | 108 | 114 | 109 |
| 3 | 1 | 84 | 83 | 81 |
| 3 | 2 | 95 | 97 | 93 |
| 3 | 3 | 105 | 100 | 106 |
Replicate II
| Dog | Dog | Dog | ||
|---|---|---|---|---|
| Lidocaine Brand | Dosage Strength | 1 | 2 | 3 |
| 1 | 1 | 84 | 85 | 86 |
| 1 | 2 | 95 | 97 | 90 |
| 1 | 3 | 105 | 104 | 103 |
| 2 | 1 | 80 | 82 | 84 |
| 2 | 2 | 93 | 99 | 95 |
| 2 | 3 | 110 | 102 | 100 |
| 3 | 1 | 83 | 80 | 79 |
| 3 | 2 | 92 | 96 | 93 |
| 3 | 3 | 102 | 111 | 108 |
Inserindo os dados
replica <- c(rep(1,27),rep(2,27))
A <- c(rep(1,9),rep(2,9),rep(3,9),rep(1,9),rep(2,9),rep(3,9))
B <- c(rep(1,3),rep(2,3),rep(3,3),rep(1,3),rep(2,3),rep(3,3),rep(1,3),rep(2,3),rep(3,3),rep(1,3),rep(2,3),rep(3,3),rep(1,3),rep(2,3),rep(3,3),rep(1,3),rep(2,3),rep(3,3))
C <- c(rep(1:3,18))
y <- c(96,84,85,94,99,98,101,106,98,85,84,86,95,98,97,108,114,109,84,83,81,95,97,93,105,100,106,84,85,86,95,97,90,105,104,103,80,82,84,93,99,95,110,102,100,83,80,79,92,96,93,102,111,108)
dados <- data.frame(replica, A, B, C, y)Visualizando as primeiras 10 linhas do banco de dados
dados %>%
slice_head(n=10) %>%
gt::gt()| replica | A | B | C | y |
|---|---|---|---|---|
| 1 | 1 | 1 | 1 | 96 |
| 1 | 1 | 1 | 2 | 84 |
| 1 | 1 | 1 | 3 | 85 |
| 1 | 1 | 2 | 1 | 94 |
| 1 | 1 | 2 | 2 | 99 |
| 1 | 1 | 2 | 3 | 98 |
| 1 | 1 | 3 | 1 | 101 |
| 1 | 1 | 3 | 2 | 106 |
| 1 | 1 | 3 | 3 | 98 |
| 1 | 2 | 1 | 1 | 85 |
7.1 Estatística Descritiva do banco
- Média das observações de acordo com o fator A (marcas diferentes de lidocaína).
Observação: Nessa breve estatística Descritiva não estamos considerando os blocos, ou seja, estamos colocando todos na mesma análise sem separar por blocos.
tapply(dados$y, factor(dados$A, levels = c(1,2,3), labels = c("Marca A","Marca B", "Marca C")), mean) %>%
knitr::kable(col.names = "Média")| Média | |
|---|---|
| Marca A | 95.00000 |
| Marca B | 95.61111 |
| Marca C | 93.77778 |
Visualização gráfica
ggplot2::ggplot(dados, ggplot2::aes(y = y, x = factor(A, levels = c(1,2,3), labels = c("Marca A", "Marca B", "Marca C")), fill = factor(A, levels = c(1,2,3), labels = c("Marca A", "Marca B", "Marca C"))), legend.tile = "Marcas de Lidocaína") + ggplot2::scale_fill_brewer(palette="Accent") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F, alpha = 0.8)+
ggplot2::labs(title = "Boxplot do nível de Enzimas no músculo cardíaco de cães Beagle", subtitle = "Com relação as marcas de Lidocaína")+
ggplot2::xlab("Marcas de Lidocaína") + ggplot2::ylab("Nível de Enzima")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Marcas de Lidocaína")
- Média das observações de acordo com o fator B (Níveis de Dosagem).
tapply(dados$y, factor(dados$B, levels = c(1,2,3), labels = c("Nível de Dosagem I","Nível de Dosagem II", "Nível de Dosagem III")), mean) %>%
knitr::kable(col.names = "Média")| Média | |
|---|---|
| Nível de Dosagem I | 83.94444 |
| Nível de Dosagem II | 95.33333 |
| Nível de Dosagem III | 105.11111 |
Visualização gráfica
ggplot2::ggplot(dados, ggplot2::aes(y = y, x = factor(B, levels = c(1,2,3), labels = c("Dosagem I", "Dosagem II", "Dosagem III")), fill = factor(B, levels = c(1,2,3), labels = c("Dosagem I", "Dosagem II", "Dosagem III"))), legend.tile = "Níveis de Dosagem") + ggplot2::scale_fill_brewer(palette="Set1") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F, alpha = 0.8)+
ggplot2::labs(title = "Boxplot do nível de Enzimas no músculo cardíaco de cães Beagle", subtitle = "Com relação aos níveis de dosagem")+
ggplot2::xlab("Níveis de Dosagem") + ggplot2::ylab("Nível de Enzima")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Níveis de Dosagem")
- Média das observações de acordo com o fator C (Cachorros).
tapply(dados$y, factor(dados$C, levels = c(1,2,3), labels = c("Cão I","Cão II", "Cão III")), mean) %>%
knitr::kable(col.names = "Média")| Média | |
|---|---|
| Cão I | 94.83333 |
| Cão II | 95.61111 |
| Cão III | 93.94444 |
Visualização gráfica
ggplot2::ggplot(dados, ggplot2::aes(y = y, x = factor(C, levels = c(1,2,3), labels = c("Cão I", "Cão II", "Cão III")), fill = factor(C, levels = c(1,2,3), labels = c("Cão I", "Cão II", "Cão III"))), legend.tile = "Cachorros diferentes") + ggplot2::scale_fill_brewer(palette="Set2") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F, alpha = 0.8)+
ggplot2::labs(title = "Boxplot do nível de Enzimas no músculo cardíaco de cães Beagle", subtitle = "Com relação aos três diferentes cães")+
ggplot2::xlab("Cachorros") + ggplot2::ylab("Nível de Enzima")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Cachorro")
- Verificando a média das observações de acordo com os blocos, ou réplicas para observar se há diferença em média
tapply(dados$y, factor(dados$replica,levels = c(1,2),labels = c("Bloco I", "Bloco II")), mean) %>%
knitr::kable(col.names = "Média")| Média | |
|---|---|
| Bloco I | 95.59259 |
| Bloco II | 94.00000 |
Visualização Gráfica
ggplot2::ggplot(dados, ggplot2::aes(y = y, x = factor(replica, levels = c(1,2), labels = c("Bloco I", "Bloco II")), fill = factor(replica, levels = c(1,2), labels = c("Bloco I", "Bloco II"))), legend.tile = "Blocos") + ggplot2::scale_fill_brewer(palette="Set3") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F, alpha = 0.8)+
ggplot2::labs(title = "Boxplot do nível de Enzimas no músculo cardíaco de cães Beagle", subtitle = "Com relação aos blocos ou replicas")+
ggplot2::xlab("Blocos") + ggplot2::ylab("Nível de Enzima")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Blocos")
7.2 Aplicação do planejamento fatorial 3^3
Neste caso será utilizada a função fat3.rbd, pois o exemplo possui blocos, da biblioteca ExpDes.
ExpDes::fat3.rbd(A,B,C,replica,y, quali = c(TRUE, TRUE, TRUE), mcomp = "tukey", fac.names = c("Marcas de Lidocaína", "Níveis de Dosagem", "Cães"), sigT = 0.05, sigF = 0.05)------------------------------------------------------------------------
Legend:
FACTOR 1: Marcas de Lidocaína
FACTOR 2: Níveis de Dosagem
FACTOR 3: Cães
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc
Block 1 34.24074 34.24074 2.8694
Marcas de Lidocaína 2 31.37037 15.68519 1.3144
Níveis de Dosagem 2 4040.03704 2020.01852 169.2793
Cães 2 25.03704 12.51852 1.0491
Marcas de Lidocaína*Níveis de Dosagem 4 112.51852 28.12963 2.3573
Marcas de Lidocaína*Cães 4 11.85185 2.96296 0.2483
Níveis de Dosagem*Cães 4 56.51852 14.12963 1.1841
Marcas de Lidocaína*Níveis de Dosagem*Cães 8 68.92593 8.61574 0.722
Residuals 26 310.25926 11.93305
Total 52 4690.75926
Pr>Fc
Block 0.1022
Marcas de Lidocaína 0.2859
Níveis de Dosagem 0
Cães 0.3646
Marcas de Lidocaína*Níveis de Dosagem 0.0799
Marcas de Lidocaína*Cães 0.9081
Níveis de Dosagem*Cães 0.341
Marcas de Lidocaína*Níveis de Dosagem*Cães 0.6707
Residuals
Total
------------------------------------------------------------------------
CV = 3.64 %
------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value: 0.03848143
WARNING: at 5% of significance, residuals can not be considered normal!
------------------------------------------------------------------------
No significant interaction: analyzing the simple effect
------------------------------------------------------------------------
Marcas de Lidocaína
According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
Levels Means
1 1 95.00000
2 2 95.61111
3 3 93.77778
------------------------------------------------------------------------
Níveis de Dosagem
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a 3 105.1111
b 2 95.33333
c 1 83.94444
------------------------------------------------------------------------
Cães
According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
Levels Means
1 1 94.83333
2 2 95.61111
3 3 93.94444
------------------------------------------------------------------------Resposta: Ao nível de significância de 5% obtivemos que apenas os níveis de dosagem são significativos, ou seja, apenas eles produzem diferença dado aos diferentes níveis de aplicação da lidocaína na enzima muscular dos cães beagle, as demais não foram significativas, nem a interação de ambas.
Outro resultado que pode ser destacado é que os resíduos não estão distribuidos normalmente.
8 Planejamento Gage R&R
Exemplo retirado do livro do Douglas C. Montgomery, que tem por título Introdução ao controle estatístico de qualidade, a versão que foi utilizada foi a sétima versão.
Segue a descrição do problema utilizado para a aplicação do modelo de planejamneto Gage R&R.
Em um estudo para o isolamento da repetibilidade e da reprodutibilidade do medidor, dois operadores usam o mesmo medidor para medir dez peças, três vezes cada. Os dados são mostrados na tabela a seguir
Operador I
Operador II
Número da peça
Medidas
Medidas
| 1 | 2 | 3 | 1 | 2 | 3 | |
|---|---|---|---|---|---|---|
| 1 | 50 | 49 | 50 | 50 | 48 | 51 |
| 2 | 52 | 52 | 51 | 51 | 51 | 51 |
| 3 | 53 | 50 | 50 | 54 | 52 | 51 |
| 4 | 49 | 51 | 50 | 48 | 50 | 51 |
| 5 | 48 | 49 | 48 | 48 | 49 | 48 |
| 6 | 52 | 50 | 50 | 52 | 50 | 50 |
| 7 | 51 | 51 | 51 | 51 | 50 | 50 |
| 8 | 52 | 50 | 49 | 53 | 48 | 50 |
| 9 | 50 | 51 | 50 | 51 | 48 | 49 |
| 10 | 47 | 46 | 49 | 46 | 47 | 48 |
Inserindo os valores
operador <- c(rep(1,3),rep(2,3),rep(1,3),rep(2,3),rep(1,3),rep(2,3),rep(1,3),rep(2,3),rep(1,3),rep(2,3),rep(1,3),rep(2,3),rep(1,3),rep(2,3),rep(1,3),rep(2,3),rep(1,3),rep(2,3),rep(1,3),rep(2,3))
peça <- c(rep(1,6),rep(2,6),rep(3,6),rep(4,6),rep(5,6),rep(6,6),rep(7,6),rep(8,6),rep(9,6),rep(10,6))
repetição <- c(rep(1:3,20))
y <- c(50,49,50,50,48,51,52,52,51,51,51,51,53,50,50,54,52,51,49,51,50,48,50,51,48,49,48,48,49,48,52,50,50,52,50,50,51,51,51,51,50,50,52,50,49,53,48,50,50,51,50,51,48,49,47,46,49,46,47,48)Criando um data frame
dados <- data.frame(operador, repetição, peça, y)Visualizando as primeiras 10 linhas do banco de dados
dados %>%
slice_head(n=10) %>%
gt::gt()| operador | repetição | peça | y |
|---|---|---|---|
| 1 | 1 | 1 | 50 |
| 1 | 2 | 1 | 49 |
| 1 | 3 | 1 | 50 |
| 2 | 1 | 1 | 50 |
| 2 | 2 | 1 | 48 |
| 2 | 3 | 1 | 51 |
| 1 | 1 | 2 | 52 |
| 1 | 2 | 2 | 52 |
| 1 | 3 | 2 | 51 |
| 2 | 1 | 2 | 51 |
8.1 Estatística Descritiva do banco de dados
- Média da medida da peça de acordo com o operador.
tapply(dados$y, factor(dados$operador, levels = c(1,2), labels = c("Operador 1", "Operador 2")),mean) %>%
round(2) %>%
knitr::kable(col.names = "Média da medida")| Média da medida | |
|---|---|
| Operador 1 | 50.03 |
| Operador 2 | 49.87 |
Visualização gráfica
ggplot2::ggplot(dados, ggplot2::aes(y = y, x = factor(operador, levels = c(1,2), labels = c("Operador I", "Operador II")), fill = factor(operador, levels = c(1,2), labels = c("Operador I", "Operador II"))), legend.tile = "Operador") + ggplot2::scale_fill_brewer(palette="Paired") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F, alpha = 0.8)+
ggplot2::labs(title = "Boxplot da medida da peça", subtitle = "Com relação aos operadores")+
ggplot2::xlab("Operadores") + ggplot2::ylab("Medida da Peça")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Operador")
Comentário: O operador 2 tem uma variabilidade maior do que o operador 1.
- Média da medida da peça de acordo com a medida.
tapply(dados$y, factor(dados$repetição, levels = c(1,2,3), labels = c("Medida I", "Medida II", "Medida III")),mean) %>%
round(2) %>%
knitr::kable(col.names = "Média da Medida")| Média da Medida | |
|---|---|
| Medida I | 50.40 |
| Medida II | 49.60 |
| Medida III | 49.85 |
Visualização gráfica
ggplot2::ggplot(dados, ggplot2::aes(y = y, x = factor(repetição, levels = c(1,2,3), labels = c("Medida I", "Medida II", "Medida III")), fill = factor(repetição, levels = c(1,2,3), labels = c("Medida I", "Medida II","Medida III"))), legend.tile = "Medida") + ggplot2::scale_fill_brewer(palette="Set2") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F, alpha = 0.8)+
ggplot2::labs(title = "Boxplot da medida da peça", subtitle = "Com relação as três medidas diferentes")+
ggplot2::xlab("Medidas") + ggplot2::ylab("Medida da Peça")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Medida")
Comentário: É observado que a medida 1 os operadores tenha uma maior dificuldade, pois a variabilidade em comparação com as demais medidas é maior.
- Operadores e a medida
ggplot2::ggplot(dados, ggplot2::aes(y = y, x = factor(repetição, levels = c(1,2,3), labels = c("Medida I", "Medida II", "Medida III")), fill = factor(operador, levels = c(1,2), labels = c("Operador I", "Operador II"))), legend.tile = "Operador") + ggplot2::scale_fill_brewer(palette="Set2") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F, alpha = 0.8)+
ggplot2::labs(title = "Boxplot da medida da peça", subtitle = "Com relação as três medidas diferentes e os operadores")+
ggplot2::xlab("Medidas") + ggplot2::ylab("Medida da Peça")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Operador")
8.2 Aplicação do planejamento Gage R&R
- Para a aplicação do planejamento gage R&R, vamos utilizar a função
ss.rrdo pacoteSixSigma.
SixSigma::ss.rr(var = y,
part = peça,
appr = operador,
data = dados,
main = "Gage R&R - Exemplo do Montgomery",
sub = "Book - Planejamento de Experimentos I")Complete model (with interaction):
Df Sum Sq Mean Sq F value Pr(>F)
peça 9 99.02 11.002 18.280 9.38e-05
operador 1 0.42 0.417 0.692 0.427
peça:operador 9 5.42 0.602 0.401 0.927
Repeatability 40 60.00 1.500
Total 59 164.85
alpha for removing interaction: 0.05
Reduced model (without interaction):
Df Sum Sq Mean Sq F value Pr(>F)
peça 9 99.02 11.002 8.241 2.46e-07
operador 1 0.42 0.417 0.312 0.579
Repeatability 49 65.42 1.335
Total 59 164.85
Gage R&R
VarComp %Contrib
Total Gage R&R 1.335034 45.31
Repeatability 1.335034 45.31
Reproducibility 0.000000 0.00
operador 0.000000 0.00
Part-To-Part 1.611136 54.69
Total Variation 2.946170 100.00
StdDev StudyVar %StudyVar
Total Gage R&R 1.155437 6.932620 67.32
Repeatability 1.155437 6.932620 67.32
Reproducibility 0.000000 0.000000 0.00
operador 0.000000 0.000000 0.00
Part-To-Part 1.269305 7.615833 73.95
Total Variation 1.716441 10.298647 100.00
Number of Distinct Categories = 1 Warning in widths.x[pos.widths[[nm]]] <- widths.settings[[nm]] *
widths.defaults[[nm]]$x: número de itens para para substituir não é um múltiplo
do comprimento do substituto
Resposta: Podemos ver que existe uma variação devido as peças, e que também existe devido a repetibilidade, o que já havia indícios acima, quando foi realizado a estatística descritiva do banco de dados. Sendo assim, é importante, verificar a dificuldade dos operadores na primeira medição, pois ela é a que tem maior variabilidade se comparada com as demais.
9 Medidas Repetidas
Exemplo retirado de um arquivo na internet, que pode ser encontrado no link: http://jaguar.fcav.unesp.br/euclides/AL_2012/CURSO_SAS_MRT/Apostila_MRT_p1.pdf
- Cinco variedades de uma cultura (tratamentos) avaliadas ao longo do tempo (0, 30, 60, 90, 120 e 150 dias), em um experimento inteiramente casualizado, com 3 repetições. Os dados observados são apresentados a seguir:
Tempo em dias
| Variedades | Repetição | 0 | 30 | 60 | 90 | 120 | 150 |
|---|---|---|---|---|---|---|---|
| V1 | 1 | 11,82 | 14,86 | 13,84 | 15,53 | 15,49 | 15,82 |
| V1 | 2 | 12,07 | 14,44 | 13,92 | 15,47 | 16,34 | 18,64 |
| V1 | 3 | 12,45 | 14,18 | 13,76 | 14,35 | 15,93 | 16,52 |
| V2 | 1 | 12,47 | 15,19 | 15,02 | 15,54 | 18,53 | 15,76 |
| V2 | 2 | 11,07 | 13,38 | 14,61 | 14,07 | 17,84 | 16,91 |
| V2 | 3 | 10,66 | 14,22 | 13,54 | 15,93 | 15,94 | 16,81 |
| V3 | 1 | 12,92 | 14,49 | 13,40 | 13,68 | 16,26 | 14,78 |
| V3 | 2 | 10,29 | 14,42 | 14,62 | 15,84 | 16,29 | 15,62 |
| V3 | 3 | 12,83 | 13,92 | 15,69 | 15,12 | 14,91 | 17,22 |
| V4 | 1 | 11,96 | 14,71 | 14,98 | 15,25 | 16,21 | 15,53 |
| V4 | 2 | 13,38 | 15,07 | 13,62 | 15,39 | 15,77 | 16,51 |
| V4 | 3 | 10,37 | 15,78 | 13,33 | 14,50 | 16,66 | 16,34 |
| V5 | 1 | 11,05 | 13,18 | 14,61 | 14,88 | 16,51 | 16,36 |
| V5 | 2 | 10,63 | 13,14 | 14,53 | 14,21 | 16,57 | 15,24 |
| V5 | 3 | 13,43 | 14,08 | 14,23 | 14,11 | 15,86 | 17,50 |
Inserindo os dados
valores<-c(11.82,14.86,13.84,15.53,15.49,15.82,
12.07,14.44,13.92,15.47,16.34,18.64,
12.45,14.18,13.76,14.35,15.93,16.52,
12.47,15.19,15.02,15.54,18.53,15.76,
11.07,13.38,14.61,14.07,17.84,16.91,
10.66,14.22,13.54,15.93,15.94,16.81,
12.92,14.49,13.40,13.68,16.26,14.78,
10.29,14.42,14.62,15.84,16.29,15.62,
12.83,13.92,15.69,15.12,14.91,17.22,
11.96,14.71,14.98,15.25,16.21,15.53,
13.38,15.07,13.62,15.39,15.77,16.51,
10.37,15.78,13.33,14.50,16.66,16.34,
11.05,13.18,14.61,14.88,16.51,16.36,
10.63,13.14,14.53,14.21,16.57,15.24,
13.43,14.08,14.23,14.11,15.86,17.50)
tempo <- c(rep(seq(0,150,by=30),15))
tipo <- c(rep("V1",18),rep("V2",18),rep("V3",18),rep("V4",18),rep("V5",18))
repetição <- c(rep(c(rep(1,6),rep(2,6),rep(3,6)),5))Estruturando o banco de dados, e visualizando algumas linhas
banco <- data.frame(tipo, repetição, tempo, valores)
banco %>%
slice_head(n=10) %>%
gt::gt()| tipo | repetição | tempo | valores |
|---|---|---|---|
| V1 | 1 | 0 | 11.82 |
| V1 | 1 | 30 | 14.86 |
| V1 | 1 | 60 | 13.84 |
| V1 | 1 | 90 | 15.53 |
| V1 | 1 | 120 | 15.49 |
| V1 | 1 | 150 | 15.82 |
| V1 | 2 | 0 | 12.07 |
| V1 | 2 | 30 | 14.44 |
| V1 | 2 | 60 | 13.92 |
| V1 | 2 | 90 | 15.47 |
9.1 Estatística Descritiva do Exemplo
Realizando a estatística descritiva de modo mais amplo, utilizando uma função do pacote DescTools, chamada Desc.
DescTools::Desc(banco, digits = 3, plotit = F)------------------------------------------------------------------------------
Describe banco (data.frame):
data frame: 90 obs. of 4 variables
90 complete cases (100.0%)
Nr ColName Class NAs Levels
1 tipo character .
2 repetição numeric .
3 tempo numeric .
4 valores numeric .
------------------------------------------------------------------------------
1 - tipo (character)
length n NAs unique levels dupes
90 90 0 5 5 y
100.0% 0.0%
level freq perc cumfreq cumperc
1 V1 18 20.000% 18 20.000%
2 V2 18 20.000% 36 40.000%
3 V3 18 20.000% 54 60.000%
4 V4 18 20.000% 72 80.000%
5 V5 18 20.000% 90 100.000%
------------------------------------------------------------------------------
2 - repetição (numeric)
length n NAs unique 0s mean meanCI'
90 90 0 3 0 2.000 1.828
100.0% 0.0% 0.0% 2.172
.05 .10 .25 median .75 .90 .95
1.000 1.000 1.000 2.000 3.000 3.000 3.000
range sd vcoef mad IQR skew kurt
2.000 0.821 0.411 1.483 2.000 0.000 -1.533
value freq perc cumfreq cumperc
1 1 30 33.3% 30 33.3%
2 2 30 33.3% 60 66.7%
3 3 30 33.3% 90 100.0%
' 95%-CI (classic)
------------------------------------------------------------------------------
3 - tempo (numeric)
length n NAs unique 0s mean meanCI'
90 90 0 6 15 75.000 64.209
100.0% 0.0% 16.7% 85.791
.05 .10 .25 median .75 .90 .95
0.000 0.000 30.000 75.000 120.000 150.000 150.000
range sd vcoef mad IQR skew kurt
150.000 51.522 0.687 66.717 90.000 0.000 -1.307
value freq perc cumfreq cumperc
1 0 15 16.7% 15 16.7%
2 30 15 16.7% 30 33.3%
3 60 15 16.7% 45 50.0%
4 90 15 16.7% 60 66.7%
5 120 15 16.7% 75 83.3%
6 150 15 16.7% 90 100.0%
' 95%-CI (classic)
------------------------------------------------------------------------------
4 - valores (numeric)
length n NAs unique 0s mean meanCI'
90 90 0 83 0 14.674 14.307
100.0% 0.0% 0.0% 15.042
.05 .10 .25 median .75 .90 .95
11.059 12.412 13.700 14.820 15.835 16.525 17.081
range sd vcoef mad IQR skew kurt
8.350 1.756 0.120 1.557 2.135 -0.446 0.152
lowest : 10.290, 10.370, 10.630, 10.660, 11.050
highest: 17.220, 17.500, 17.840, 18.530, 18.640
' 95%-CI (classic)- Observando os valores da porcentagem de cana de açucar de acordo com os tipos variados de cana de açucar.
tapply(banco$valores, factor(banco$tipo, levels = c("V1","V2","V3","V4","V5"), labels = c("Tipo 1", "Tipo 2", "Tipo 3", "Tipo 4", "Tipo 5")), mean) %>%
knitr::kable(col.names = "Valor médio", caption = "Porcentagem de Açucar na cana de acordo com o tipo")| Valor médio | |
|---|---|
| Tipo 1 | 14.74611 |
| Tipo 2 | 14.86056 |
| Tipo 3 | 14.57222 |
| Tipo 4 | 14.74222 |
| Tipo 5 | 14.45111 |
ggplot2::ggplot(banco, ggplot2::aes(y = valores, x = factor(tipo,levels = c("V1","V2","V3","V4","V5"), labels = c("Tipo 1", "Tipo 2", "Tipo 3", "Tipo 4", "Tipo 5")), fill = factor(tipo, levels = c("V1","V2","V3","V4","V5"), labels = c("Tipo 1", "Tipo 2", "Tipo 3", "Tipo 4", "Tipo 5"))), legend.tile = "Tipos de Cana de Açucar") + ggplot2::scale_fill_brewer(palette="Set1") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F, alpha = 0.8)+
ggplot2::labs(title = "Boxplot da Porcentagem de açucar na cana de Açucar", subtitle = "Com relação aos tipos diferentes de cana de açucar")+
ggplot2::xlab("Tipos de Cana de Açucar") + ggplot2::ylab("Porcentagem de açucar")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Tipos de Cana de Açucar")
- Observando os valores da porcentagem de cana de açucar de acordo com o tempo de desenvolvimento da cultura.
tapply(banco$valores, banco$tempo, mean) %>%
knitr::kable(col.names = "Valor médio", caption = "Porcentagem de Açucar na cana de acordo com o tempo de desenvolvimento da cultura")| Valor médio | |
|---|---|
| 0 | 11.82667 |
| 30 | 14.33733 |
| 60 | 14.24667 |
| 90 | 14.92467 |
| 120 | 16.34067 |
| 150 | 16.37067 |
ggplot2::ggplot(banco, ggplot2::aes(y = valores, x = factor(tempo), fill = factor(tempo)), legend.tile = "Tempo de desenvolvimento da cultura") + ggplot2::scale_fill_brewer(palette="PuRd") +
ggplot2::geom_boxplot(outlier.colour = "red", notch = F, alpha = 0.8)+
ggplot2::labs(title = "Boxplot da Porcentagem de açucar na cana de Açucar", subtitle = "Com relação ao tempo de desenvolvimento da cultura")+
ggplot2::xlab("Tempo de desenvolvimento da cultura") + ggplot2::ylab("Porcentagem de açucar")+
ggplot2::theme_minimal() + ggplot2::labs(fill = "Tempo")
Visualização gráfica
ggplot2::ggplot(banco) +
ggplot2::aes(x = tempo, y = valores, color = factor(tipo,levels = c("V1","V2","V3","V4","V5"), labels = c("Tipo 1", "Tipo 2", "Tipo 3", "Tipo 4", "Tipo 5"))) + ggplot2::geom_line() +
ggplot2::scale_color_brewer(palette="Set1")+
ggplot2::labs(title = "Porcentagem de açucar na cana de açucar",subtitle ="Com relação ao tempo de desenvolvimento e tipo diferente de cana",color = "Tipo de cana") +
ggplot2::xlab("Tempo") + ggplot2::ylab("Porcentagem de Açucar") 
9.2 Aplicação das medidas repetidas
modelo <- nlme::lme(valores ~ tempo + repetição, random = ~ 1 | tipo, data = banco)
summary(modelo)Linear mixed-effects model fit by REML
Data: banco
AIC BIC logLik
279.8511 292.1806 -134.9255
Random effects:
Formula: ~1 | tipo
(Intercept) Residual
StdDev: 3.157454e-05 1.011496
Fixed effects: valores ~ tempo + repetição
Value Std.Error DF t-value p-value
(Intercept) 12.589206 0.3223919 83 39.04939 0.0000
tempo 0.028008 0.0020810 83 13.45853 0.0000
repetição -0.007667 0.1305836 83 -0.05871 0.9533
Correlation: Warning in abbreviate(colnames(x), minlength = rdig + 3): abbreaviate usado com
caracteres não-ASCII (Intr) tempo
tempo -0.484
repetição -0.810 0.000
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.25791500 -0.62395922 0.02273231 0.65643037 2.55813653
Number of Observations: 90
Number of Groups: 5 10 Superfície de Resposta
Exemplo retirado do livro do Douglas C. Montgomery, que tem por título Design and Analysis of Experiments, a versão que foi utilizada foi a oitava versão.
Segue a descrição do problema utilizado para a aplicação do modelo de superfície de resposta.
Os dados mostrados na Tabela abaixo foram coletados em um experimento para otimizar o crescimento de cristais em função de três variáveis x1, x2 e x3. Grandes valores de y (rendimento em gramas) são desejáveis. Ajustar um modelo de segunda ordem e analisar o ajustado superfície.
| x_1 | x_2 | x_3 | y |
|---|---|---|---|
| -1 | -1 | -1 | 66 |
| -1 | -1 | 1 | 70 |
| -1 | 1 | -1 | 78 |
| -1 | 1 | 1 | 60 |
| 1 | -1 | -1 | 80 |
| 1 | -1 | 1 | 70 |
| 1 | 1 | -1 | 100 |
| 1 | 1 | 1 | 75 |
| -1.682 | 0 | 0 | 100 |
| 1.682 | 0 | 0 | 80 |
| 0 | -1.682 | 0 | 68 |
| 0 | 1.682 | 0 | 63 |
| 0 | 0 | -1.682 | 65 |
| 0 | 0 | 1.682 | 82 |
| 0 | 0 | 0 | 113 |
| 0 | 0 | 0 | 100 |
| 0 | 0 | 0 | 118 |
| 0 | 0 | 0 | 88 |
| 0 | 0 | 0 | 100 |
| 0 | 0 | 0 | 85 |
Inserindo os dados
x1 <- c(rep(-1,4),rep(1,4),-1.682,1.682,rep(0,10))
x2 <- c(rep(-1,2),rep(1,2),rep(-1,2),rep(1,2),0,0,-1.682,1.682,rep(0,8))
x3 <- c(-1,1,-1,1,-1,1,-1,1,0,0,0,0,-1.682,1.682,rep(0,6))
y <- c(66,70,78,60,80,70,100,75,100,80,68,63,65,82,113,100,118,88,100,85)
banco <- data.frame(x1,x2,x3,y)Visualizando algumas linhas do banco de dados
gt::gt(slice_head(banco, n=10))| x1 | x2 | x3 | y |
|---|---|---|---|
| -1.000 | -1 | -1 | 66 |
| -1.000 | -1 | 1 | 70 |
| -1.000 | 1 | -1 | 78 |
| -1.000 | 1 | 1 | 60 |
| 1.000 | -1 | -1 | 80 |
| 1.000 | -1 | 1 | 70 |
| 1.000 | 1 | -1 | 100 |
| 1.000 | 1 | 1 | 75 |
| -1.682 | 0 | 0 | 100 |
| 1.682 | 0 | 0 | 80 |
10.1 Realização do Experimento de Superfície de Resposta
- Como a questão sugeriu, uma superfície de segunda ordem.
resposta <- rsm::rsm(y~SO(x1,x2,x3), data = banco)
summary(resposta)
Call:
rsm(formula = y ~ SO(x1, x2, x3), data = banco)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 100.6663 5.5638 18.0930 5.703e-09 ***
x1 1.2710 3.6912 0.3443 0.737726
x2 1.3611 3.6912 0.3687 0.720014
x3 -1.4940 3.6912 -0.4048 0.694182
x1:x2 2.8750 4.8231 0.5961 0.564363
x1:x3 -2.6250 4.8231 -0.5443 0.598190
x2:x3 -4.6250 4.8231 -0.9589 0.360205
x1^2 -3.7679 3.5929 -1.0487 0.318990
x2^2 -12.4278 3.5929 -3.4590 0.006133 **
x3^2 -9.6001 3.5929 -2.6720 0.023412 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Multiple R-squared: 0.663, Adjusted R-squared: 0.3598
F-statistic: 2.186 on 9 and 10 DF, p-value: 0.1194
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
FO(x1, x2, x3) 3 77.9 25.95 0.1395 0.9341
TWI(x1, x2, x3) 3 292.4 97.46 0.5237 0.6757
PQ(x1, x2, x3) 3 3291.7 1097.25 5.8961 0.0139
Residuals 10 1861.0 186.10
Lack of fit 5 1001.6 200.33 1.1656 0.4353
Pure error 5 859.3 171.87
Stationary point of response surface:
x1 x2 x3
0.2597353 0.1108581 -0.1400280
Eigenanalysis:
eigen() decomposition
$values
[1] -3.079142 -8.952298 -13.764404
$vectors
[,1] [,2] [,3]
x1 0.9413839 0.3309866 -0.06514712
x2 0.2100461 -0.4240105 0.88096295
x3 -0.2639639 0.8430083 0.46867915Visualização Gráfica
contour(resposta, ~x1+x2+x3, image=TRUE, img.col=terrain.colors(50))
persp(resposta, ~ x1 + x2 + x3, col=terrain.colors(50), contours="colors")
