Caracterização do problema

Como foi dito anteriormente, em razão de melhorarmos a quantidade de óleo temos que melhorar a sua qualidade, ou seja, a quantidade extraída por amendoim. Queremos entender estatísticamente a quantidade de óleo em função de: pressão do CO2, temperatura do CO2, o fluxo de CO2, a umidade do amendoim e o tamanho do amendoim.

Análise Exploratória dos Dados

Veremos na tabela seguinte que temos para cada varíavel regressora duas opções ou dois estados. O experimento foi conduzido de tal maneira, sendo isto comum, pois estamos pensando de maneira industrial a exemplo de uma máquina que foi programada para aumentar de 125 bar em 125 bar a pressão do CO2. Não há opções verdadeiramente contínuas.

Resultado dos Experimentos

x1	x2	x3	x4	x5	y
415	95	5	60	1.28	80
550	25	5	40	4.05	21
550	95	15	60	1.28	96
550	25	15	40	1.28	66
550	25	15	60	4.05	21
550	95	5	40	1.28	99
415	25	15	60	1.28	63
550	95	5	60	4.05	33
550	25	5	60	1.28	74
415	25	5	40	1.28	63

Visualmente pela tabela não podemos afirmar nada, além do que já foi dito, então prosseguiremos com nossa exploração dos dados, para assim formularmos nossas hipóteses sobre como seria a relação entre y e os x’s.

Sendo assim, iremos pensar em cada setup das variáveis regressoras (x1, x2,…) como 0 ou 1, pois queremos apenas medir as diferenças de se mudar o estado de cada variável em relação a outra.

Resultado dos Experimentos

x1	x2	x3	x4	x5	y
0	1	0	1	0	80
0	1	1	1	1	44
0	0	1	0	1	24
0	1	0	0	1	36
1	1	0	0	0	99
1	0	1	1	1	21
0	0	1	1	0	63
1	1	0	1	1	33
0	0	0	1	1	23
1	0	1	0	0	66

Primeirmente, vemos aqui que embora a mediana do rendimento de óleo extraído tenha sido maior quando se aumenta a pressão, também temos um aumento da variância. Mesmo que a maioria dos dados tenha ficado concentrado em rendimentos maiores do óleo, é preciso notar esse aumento da variabilidade. Esperamos que o aumento de pressão leve a um aumento do óleo extraído.

De fato, o aumento da temperatura do CO2, de uma temperatura ambiente de 25 graus Celsius, para 95 graus Celsius, a qual é muitíssimo quente, aumenta o rendimento de òleo. Não somente a mediana indica isto, dizendo que 50% dos dados estão acima de 60% do rendimento de óleo de amendoim, porém a variância é muito mais bem comportada, tendo um comportamento bastante gaussiano. Esperamos que o aumento da temperatura, aumente, significativamente, de óleo extraído.

Embora haja um aumento de 50% no fluxo, há uma diminuição na mediana do rendimento de óleo extraido. Além disto, há um aumento da variabilidade. A expressividade disso, contudo, não é muitíssimo grande entre as duas opções de fluxo do CO2. Esperamos que o aumento de fluxo, diminua muito pouco a quantidade de óleo extraído.

Naturalmente, esperamos que a umidade do amendoim seja um influenciador na quantidade de óleo extraída. Vemos que quando se tem um amendoim mais úmida, a mediana fica mais acima, ou seja, pelo menos 50% dos dados estão tem um rendimento maior do que quando a umidade é menor. Além diso, a variância é claramente menor, ou seja, o aumento da umidade produz tem mais chance de não só produzir mais óleo, porém é mais certo de produzir mais óleo. Mesmo assim, a discrepância não é enorme o que traz uma leve surpresa aos leigos.

Aqui é preciso uma certa interpretação, não só estatística, mas do próprio procedimento industrial. Nós estamos medindo o rendimento de óleo por comboio de amendoins. Esperávamos que quanto maior o tamanho do amendoim, maior a sua quantidade de óleo, porém, não só os dados dizem o contrário, como podemos ver pela diferença gritante entre as medianas, mas como também essa intuitividade mostra que se desconhece o procedimento como um todo. O tamanho do amendoim não necessariamente determina sua densidade, a quantidade verdadeira de matéria orgânica está determinado pelo tamanho e densidade, não somente um ou outro. Concluímos visivelmente que a redução do tamanho do amendoim faça com que o óleo extraído seja mais abudante.

Modelo Proposto

Faremos um modelo de regressão linear múltiplo. Ou seja, reescreveremos \(y\) como: \[ y = \begin{bmatrix} y_1\\ ...\\ y_{16} \end{bmatrix}, \] \[ X = \begin{bmatrix} x_{1,1} &... &... &x_{5,1} \\ ... & ... & ... &... \\ ... & ... & ... &... \\ x_{1,16} &... &... &x_{5,16} \end{bmatrix}, \],

\[ \beta = \begin{bmatrix} \beta_1\\ ...\\ \beta_{5} \end{bmatrix}, \]

\[ \epsilon = \begin{bmatrix} \epsilon_1\\ ...\\ \epsilon_{16} \end{bmatrix} \],

\[ y = X\beta + \epsilon \]

tal que a nossa predição será dada por: \(\hat{y} = X\hat{\beta}\), o erro será descrito como e terá a seguinte distribuição: \(\epsilon \sim N(0, \sigma^{2})\) e Sendo estes, y, x1, … , x5, os mesmo daquela tabela inicial.

Desejamos achar os betas, os nossos coeficiente lineares para assim acharmos as nossas predições sobre o impacto marginal de cada varíavel em y. Além disto, diremos que cada variável regressora, ou seja, x1, …, x5, é lineramente independente da outra.

Modelo Ajustado

Minimizando o \(\sum_{i = 1}^{i = 16}{\epsilon^2}\) teremos que \(\hat{\beta} = (X'X)^{-1}Xy\). Fazendo esta conta, teremos que o modelo ajustado será: \[ \hat{y} = 66.625 + 19.750*x2 -44.500*x5 \] com as demais variáveis não sendo estatísticamente significantes para o modelo. Estes valores ligados a x2 e a x5 representam o impacto marginal de cada diferente população gerada pelas nossas variáveis explicativas no rendimento da extração de óleo de amendoim, ou seja, no valor médio esperado incrementado no valor de y por estarmos em outra população. Analisamos estes resultados através do programa de software \(R\).

Começamos colocando todas as variáveis no modelo, pois entendemos que elas tem a sua contribuição mesmo que algumas delas suspeitamos que sejam pequenas. Eis o modelo inicial:

## 
## Call:
## lm(formula = y ~ x2 + x5 + x1 + x3 + x4, data = .)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.250  -4.438   0.125   5.250   9.500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   62.250      4.939  12.604 1.84e-07 ***
## x2            19.750      4.033   4.897 0.000625 ***
## x5           -44.500      4.033 -11.035 6.40e-07 ***
## x1             7.500      4.033   1.860 0.092544 .  
## x3             1.250      4.033   0.310 0.762949    
## x4             0.000      4.033   0.000 1.000000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.065 on 10 degrees of freedom
## Multiple R-squared:  0.9372, Adjusted R-squared:  0.9058 
## F-statistic: 29.86 on 5 and 10 DF,  p-value: 1.055e-05

Na primeira parte vemos que o modelo proposto é aceitável pela estatística F, que diz que o p-valor do modelo como um todo ser desprezível está com 1.055*e-7, ou seja, a probabilidade de este modelo estar errada é baixíssimo. Podemos analisar também que pelos testes marginais da contribuição de cada variável, tomando como nível de significância de 5%, veríamos que apenas as variáveis x2 e x5 são estatísticamente significantes para o modelo.

Além disto, podemos ver pelo ANOVA, na qual vê se a contribuição sequencial marginal, ou seja, vê-se como x2 contribui para o modelo, depois, como x5, contribui em relação a x2, depois, x1 em relação a x2 com x5, etc. Vemos que a Soma de Quadrados de Regressão corrigidos, da forma que mencionamos sequencial, ao realizarmos os testes F de cada uma, novamente temos que somente x2 e x5 são significantes ao nível de significância de 5%. As contribuições delas foram muito maiores que a dos demais. Sendo assim, reescreveremos o nosso modelo somente com x2 e x5 e jogaremos novamente no \(R\):

## 
## Call:
## lm(formula = y ~ x2 + x5, data = .)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.375  -4.188  -0.875   3.438  12.625 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   66.625      3.566  18.683 8.96e-11 ***
## x2            19.750      4.118   4.796 0.000349 ***
## x5           -44.500      4.118 -10.807 7.26e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.236 on 13 degrees of freedom
## Multiple R-squared:  0.9149, Adjusted R-squared:  0.9018 
## F-statistic: 69.89 on 2 and 13 DF,  p-value: 1.107e-07

Vemos pelos t testes marginais que as variáveis juntamente com o intercepto são todas estatísticamente significantes para o modelo proposto, além da estatística F, na qual mediu se a adequeação do modelo como um todo. Além disto, conseguimos estimar aquele \(\sigma^2\) do nosso modelo, sendo ele: 8.236 aproximadamente.

Iremos construir intervalos de confiança de nível de significância de 5% para os nossos parâmetros.

##                 2.5 %    97.5 %
## (Intercept)  58.92076  74.32924
## x2           10.85391  28.64609
## x5          -53.39609 -35.60391

De fato, nossos parâmetros são significativos, pois nenhum dos intervalos contém o 0.

Novamente, nos faremos valer do \(R\) para fazer o intervalo de confiança para uma nova observação. Só podemos fazer predições para o espaço amostral ou seja o espaço dos experimentos.

##      fit      lwr      upr
## 1 86.375 66.98643 105.7636

Vejamos os outros experimentos que também tem x2 = 95, e x5 = 1.28.

x2	x5	y
1	0	99
1	0	71
1	0	80
1	0	96

Vemos, então que o intervalo de confiança está condizente com os experimentos observados.

Análise de Resíduo e Diagnóstico

Desejamos agora analisar os resíduos que são descritos como \(e = y - \hat{y}\), podemos ver estes como realizações daquele nosso erro \(\epsilon\). Esperamos que estes sejam normais, com média 0 e variância fixa também, e ao studentizarmos eles que eles sigam uma distribuição t-student. Faremos isto pois ao studentizarmos os resíduos teremos ele com variância fixa unitária, com média 0, e também, esperaremos os seus valores estarem concentrados entre 2 e -2.

Inicialmente, desejamos ver se eles tem média zero, pois é o que esperamos segundo as hipóteses assumidas e teoria construída no modelo proposto.

## [1] -1.110223e-16

De fato, o valor é desprezível e praticamente zero, como esperávamos teóricamente. Plotaremos o resíduos studentizados para vermos se é o que estamos esperando.

Temos poucas observações, então, não esperávamos uma distribuição perfeita, porém, o seu padrão lembra muitíssimo a t-student, para somente 16 observações. Novamente, poucas observações, porém, está seguindo o padrão que desejávamos.

Prosseguiremos na análise.

Nesta análise, estamos atrás de checar a aleatoriedade dos resíduos, identificar outliers e também vermos se os resíduos estão concentrados em torno do 0, para que a média seja 0. Podemos ver que há aleatoriedade, e também que há um outliers. Ele se configura outlier pela probabilidade dele acontecer ser muitíssimo pequena, pois esperamos que 95% dos dados estejam entre 2 e -2, ele está fora. O seu índice é o 7. Vejamos como ele é:

x2	x5	y
1	0	71

Podemos ver que o seu x2 = 95 e x5 = 1.28, com y = 71. Todos os outros experimentos (índices 4, 16 e 11) que tiveram x2 = 95 e x5 = 1.28 tiveram resultados acima de 80, sendo dois deles acima de 90. De fato, ele é uma observação ‘estranha’. Tentaremos classificar este outlier.

## 
## Call:
## lm(formula = y ~ x2 + x5, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.0577  -4.8077   0.0577   4.5673  10.9423 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   67.808      3.077  22.037 4.48e-11 ***
## x2            22.115      3.641   6.074 5.55e-05 ***
## x5           -46.865      3.641 -12.872 2.20e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.017 on 12 degrees of freedom
## Multiple R-squared:  0.9413, Adjusted R-squared:  0.9315 
## F-statistic:  96.2 on 2 and 12 DF,  p-value: 4.094e-08

Vemos que os coeficientes mudaram pouco, do intercepto e do x2 porém, o x5 mudou em até uma unidade, além do R2 ter aumentado também. Suspeitamos de que seja um outlier influente.

Embora um desses pontos foi um que tivesse mais influência entres os demais, não o classificaremos como um outlier influente, pois não houve uma mudança expressiva nos coeficientes, nem houve mudanças consideráveis em relação aos desvios padrões do erros e também ao R2.

Utilizando do \(R\) veremos se o experimento 7 está muito afastado do espaço gerado pelos experimentos.

X <- model.matrix(fit)
n <- nrow(X)
p <- ncol(X)
H <- X%*%solve(t(X)%*%X)%*%t(X)
h <- diag(H)
s <- summary(fit)$sigma
r <- resid(summary(fit))
ts <- r/(s*sqrt(1-h))
di <- (1/p)*(h/(1-h))*(ts^2)

diag_media = 2*(p/n)
diag_outlier7 = h[7]

diag_outlier7 > diag_media

##     7 
## FALSE

Vemos que de fato, ele ainda permanece no espaço. Sendo assim, nosso outlier é apenas o que é um outlier, o qual não pode ser capturado pelo nosso modelo, mas não o influencia nem o pertuba.

Continuando nossa análise dos resíduos.

Vemos um leve comportamento heterocedástico, porém, pela falta de observações não podemos concluir afirmativamente, por esse gráfico. Novamente, vemos o experimento 7 como um outlier, pois prevíamos que seu valor seria muito mais do que 71, pelas suas propriedades do x2 e x5.

Nosso modelo conseguiu captar melhor as observações quando x2 = 25, do que quando x2 = 95, onde ocorre o outlier novamente. Vemos isto pelas medianas, além da diferença entre as variâncias.

Segue-se a mesma interpretação para x5. Há algo diferente neste plot, pois embora o modelo tenha descrito melhor quando x5 = 1.28, vemos ainda o outlier nessa faixa. Sendo este mesmo, aquele que suspeitávamos ser influente, embora não o fosse, o experimento 7.

## Gaussian model (lm object)

Vemos aqui o ‘envelopamento’ dos resíduos juntamente com o QQ Plot dele com as distribuição teórica t-student. Vemos que todos os pontos estão dentro da banda de confiança de 95%, através do algoritmo de simulação dada por Atkison (1981). Há uma certa tentativa expressiva de fugir da distribuição, da linha de apoio, porém, ainda estaria dentro dos parâmetros da distribuição studentizada esperada.