66
Pàg 1 de 66 REGRESSÃO LINEAR SIMPLES – Noções preliminares

REGRESSÃO LINEAR SIMPLES – Noções …palhares2.flu.angelfire.com/uniceub/ia/iat010.pdfPàg 2 de 66 REGRESSÃO LINEAR SIMPLES – Noções preliminares Quando temos duas variáveis

  • Upload
    votram

  • View
    217

  • Download
    0

Embed Size (px)

Citation preview

Pàg 1 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

Pàg 2 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Quando temos duas variáveis numéricas podemos ter:

• Relação funcional entre as duas variáveis

Y = f(X)

• Relação estatística entre as duas variáveis

A relação entre X e Y não é “perfeita”

� Exemplo de relação funcional e relação estatística

Pàg 3 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

Pàg 4 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

Pàg 5 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

Pàg 6 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

Pàg 7 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

Pàg 8 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

Pàg 9 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

Pàg 10 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

Pàg 11 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

Pàg 12 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

Após executar o script acima obtemos, dentre outros, o seguinteresultado:

� A equação de regressão acima nos diz que:

• Se não aplicarmos dose alguma de fertilizante (isto éX = 0) então esperamos que o talhão de 1 ha produza 69.867kg da cultura plantada.

• A cada aumento de uma unidade de dose de fertilizante, isto éa cada kg de fertilizante aplicado, espera-se que, além daprodução de 69.867 kg, seja produzido adicionalmente 8.442kg da cultura plantada.

Pàg 13 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� E este outro resultado apresentado abaixo, o que significa?

� Bem, estamos lidando com uma amostra (no caso n = 28talhões).

� Foi com base nesta amostra que estimamos β0 e β1 através doalgoritmo de mínimos quadrados. Este procedimentomatemático foi feito via comando lm do R.

� Ok, mas por se tratar de uma amostra, devemos realizar umteste estatístico de hipótese, ou simplesmente um teste dehipótese.

� Aliás, o procedimento lm já realiza dois testes de hipótese. Umteste para β0 e outro para β1.

Pàg 14 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� O teste de hipótese é para nós uma caixa preta.

� Mas vamos tentar entender o resultado que ele nos indica.Isto é importante para considerarmos se a reta de regressãoobtida pode ser, de fato, representativa para o nossoproblema.

� O primeiro teste de hipótese testa a seguinte hipótesedenominada H0 (hipótese nula) contra a hipótese alternativaHa:

H0 : β1 = 0 X Ha : β1 ≠ 0

� Em essência, queremos saber se existe uma reta de regressãoou, em outras palavras, se a variável preditora X (em nossoexemplo, DOSE de fertilizante) “tem relação estatística” coma produção da cultura (Y), aumentando ou diminuindoconforme o valor de X.

� Se X “tem relação estatística” Y, então vai existir um valornão nulo para β1. Então o teste tende a apontar a hipótesealternativa Ha. Lembre-se que β1 é o coeficiente angular dareta !!!!

� Mas se X não “tem relação estatística” com Y, então β1énulo, ou seja, β1= 0.

Pàg 15 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

Como interpretar o resultado do teste de hipótese

para β1 sem entender a lógica do mesmo?

Passo 1 – Com base no valor estimado de β1, que é o coeficienteda variável preditora X (em nosso exemplo, DOSE) é calculadoum t valor (t value).

Passo 2 – O t valor (t value) é uma espécie de “representante” dovalor estimado β1. Então se calcula a probabilidade de se obtereste t valor (t value) ou um valor maior (Pr > |t|), supondo queH0 é a hipótese correta !!!

Passo 3 – Normalmente se decide assim: Se Pr > |t| for menor ouigual a 5%, então aceitamos a hipótese Ha. Caso contrárioaceitamos H0.

Pàg 16 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Veja na figura abaixo o resultado que obtivemos:

Passo 1 – O valor estimado de β1 é 8.4415. O t valorcorrespondente é t = 10.613

Passo 2 – Supondo H0 correta (isto é supondo que β1 = 0, o quesignifica que não há relação linear) a probabilidade de se obter tvalue maior ou igual a 10.613 é de 6.06e-11 = 6.06 x 10-11 =0.0000000000606 = 0.000000000000606 % .

Passo 3 – Como Pr(>|t|) = 6.06 x 10-11 < 5%, então aceitamos ahipótese Ha, isto é , β1 ≠ 0, o que significa dizer que há evidênciasde que ocorre relação linear entre as variáveis PROD e DOSE.

Pàg 17 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Uma vez que o modelo de regressão linear seja construído, épossível calcular uma medida de associação linear entre avariável resposta e a variável preditora/explicativa. Esta medidadenomina-se coeficiente de determinação e é representadapelo símbolo R2, onde

0 ≤ R2 ≤ 1

� A interpretação de R 2 é a seguinte: quanto mais próximo de 1for R2, significa que a variável preditora “explica” melhor avariabilidade da variável resposta. Quanto mais próximo R2

estiver de zero, significa que a variável preditora não é “boa”para explicar a variabilidade de valores da variável resposta.

� Em nosso exemplo R2= 81,24%. A variável DOSE explica“bem” a variabilidade de valores da variável resposta PROD

Pàg 18 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

Pàg 19 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

Pàg 20 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Sobre a influência dos pontos no modelo de regressão...

Pàg 21 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

Pàg 22 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

Pàg 23 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Sobre a alvancagem (leverage) dos pontos no modelo deregressão...

Pàg 24 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Sobre os pontos discrepantes (outliers) no modelo deregressão...

Pàg 25 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Para medir o “grau” de influência de cada ponto no modelo deregressão podemos usar a distância de Cook Di ...

REGRESSÃO LINEAR SIMPLES – Noções preliminares

Pàg 26 de 66

� Veja abaixo as distâncias de Cook Di calculadas para o modelode regressão com os 28 talhões. Note que o gráfico destacou ostalhões 7, 14 e 25. Contudo, apresentam baixo valores para Di

indicando baixa influência no modelo de regressão ajustadocom os 28 pontos.

� Então podemos considerar que neste grupo de 28 pontos nãoocorre ponto influente, se considerarmos, por exemplo, pontoinfluente com Di >= 0.5.

Pàg 27 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Se aos 28 talhões acrescentarmos o talhão 29 com DOSE = 75kg e PROD = 92 kg obteremos o resultado abaixo.

� Observa-se que o talhão 29 é realmente um ponto influente(D29 = 2.25) e é um ponto discrepante (outlier) com altaalavancagem (high leverage).

� Sua permanência ou retirada produz modelos de regressãoconsideravelmente diferentes entre si.

Pàg 28 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Agora se aos 28 talhões acrescentarmos o talhão 29 com DOSE= 32 kg e PROD = 610 kg obteremos o resultado abaixo.

� Observa-se que o talhão 29 não é um ponto influente (D29 =0.182) apesar de poder ser considerado um ponto discrepante(outlier) porém com baixa alavancagem (low leverage).

� Sua permanência ou retirada não altera significativamente omodelo de regressão

Pàg 29 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Para medir o “grau” de alavancagem (leverage) de cada pontono modelo de regressão podemos usar o hat-value hi ...

� hi com valores “altos” indicam alta alavancagem (highleverage) mas não necessariamente apresentam alta influênciano modelo

� Como regra prática, hi com valor que seja maior que o dobro dovalor médio dos hi da amostra deve ser analisado maisdetalhadamente. Se a amostra de dados é muito grande,recomenda-se considerar hi com valor que seja maior que otriplo do valor médio dos hi

Pàg 30 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Veja abaixo os hat-values hi calculadas para o modelo deregressão com os 28 talhões. Note que o ponto 16 apresentará“alta alavancagem” (high leverage) se considerarmos queh16 > 2hmedio (h16 = 0.1588 e hmedio = 0.0714)

Pàg 31 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Se aos 28 talhões acrescentarmos o talhão 29 com DOSE = 75kg e PROD = 92 kg obteremos o resultado abaixo.

� Note que o talhão 29 (h29 = 0.191) realmente é um ponto dealta alavancagem (e já vimos que é também influente...)

� Já o ponto/talhão 16 foi considerado de alta alavancagem pelofato de seu valor ser mais que o dobro do valor médio (h16 =0.147 > (2* hmedio). Mas não é influente para o modelo.

Pàg 32 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Para medir o “grau” de discrepância (outlier) de cada ponto nomodelo de regressão poderíamos usar os resíduos“semiestudentizados” (semistudentized residuals) e*i ...

� Note que raiz(MSE) refere-se ao grupo de resíduos... Mas cadaresíduo ei pode apresentar variância(e consequentemente desviopadrão) diferente de alguns outros ei. Então os resíduossemiestudentizados, embora possam, não são usados paradetectar outliers...

Pàg 33 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Para medir o “grau” de discrepância (outlier) de cada ponto nomodelo de regressão podemos usar os “resíduos padronizados”(standardized residuals) ri também denominados “resíduosstudentizados internamente” (internally studentized residuals).

� Os resíduos padronizados (standardized residuals) dos 28talhões não apresentam outliers

Pàg 34 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Se aos 28 talhões acrescentarmos o talhão 29 com DOSE = 75kg e PROD = 92 kg obteremos o resultado abaixo.

� Note que o talhão 29 (r29 = -4.3657) realmente é um pontodiscrepante (outlier)

Pàg 35 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Porquê “Análise de Resíduos” ????

Pàg 36 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Mas antes.... vamos usar os resíduos paraver outras coisas em nossa massa dedados...

Pàg 37 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Vamos verificar a linearidade da função de regressão viagráfico da variável preditora X Resíduos. (No exemplo, avariável preditora é DOSE)...

� Observe que no gráfico DOSE versus Resíduos não se observaum “padrão”... os pontos estão “espalhados” de formarelativamente simétrica em relação ao “resíduo 0”... Então, apriori a linearidade é boa... Será?

� Um gráfico de resíduos consegue mostrar eventuaistendências... e queremos não encontrar tendência/padrão parapoder considerar o modelo adequado...

Pàg 38 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Tem casos que o gráfico variável preditora versus resíduo vaiapontar a não linearidade... veja abaixo...

� Note por exemplo que no gráfico acima Weight X Resíduo,poderíamos “passar uma reta com inclinação” pelos pontos...Isto é uma tendência (um padrão) que aponta que o modelo nãoestá adequado...

Pàg 39 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Olha o exemplo abaixo

(http://data.library.virginia.edu/diagnostic-plots/)

Pàg 40 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� No exemplo abaixo, o gráfico de resíduos X variável respostapredita mostra que a regressão linear simples Y ~ X não é boapara representar a relação estatística entre Y e X, pois pelográfico nota-se um padrão de distribuição dos resíduos...

Pàg 41 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Veja abaixo a ideia conceitual de homocedasticidade eheterocedasticidade...

Pàg 42 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Para verificar a homocedasticidade dos resíduos...

� Para os 28 talhões temos:

Pàg 43 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Para verificar a homocedasticidade observando-se os gráficosdevemos esperar que a linha vermelha (lowess curve) sejauniforme horizontalmente.

� No caso dos 28 talhões há indícios de heterocedasticidade.Assim, o requisito de homocedasticidade não é atendido, o quepode comprometer a validade do modelo de regressão linearconstruído...

� Veja abaixo outro exemplo:

Pàg 44 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Podemos também realizar testes de hipótese para verificar aocorrência ou não de homocedasticidade. Neste teste temos:

Hipótese nula H0: Ocorre homocedasticidade

Hipótese alternativa Ha: Não ocorre homocedasticidade(ocorre heterocedasticidade)

� Vamos realizar 3 teste de hipótese para averiguação dehomocedasticidade:

� Teste de Breusch-Pagan

� NCV teste

� Tesste de Brown-Forsythe

� Lembre-se que se o p-valor do teste de hipótese for menor que5%, adota-se o procedimento de rejeitar a hipótese H0,aceitando portanto válida a hipótese Ha. Em caso contrário,aceita-se H0.

Pàg 45 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Vejamos os resultados dos testes para os dados dos 28 talhões

� E agora ??? Em dois testes o p-valor foi maior que 5% (68,12%e 79,18%) indicando aceitação de H0 e no teste de Brown-Forsythe foi significativo (1,44%) indicando aceitação de Ha.

� Qual a sua decisão? Os resíduos apresentam ou nãohomocedasticidade?

Pàg 46 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Para verificar a normalidade dos resíduos...

� Para analisar a normalidade dos resíduos podemos usar (1) boxplot; (2) histograma e (3) diagrama ramo-e-folhas dos resíduos.Estes gráficos apresentam maior utilidade quando a amostra dedados é “grande”.

� Para os 28 talhões temos:

Pàg 47 de 66

REGRESSÃO LINEAR SIMPLES – Noções preliminares

� Para analisar graficamente a normalidade dos resíduos podemosusa-se também o gráfico Q-Q plot (Quantile-Quantile plot)

� Veja o gráfico QQ plot para o exemplo dos 28 talhões...

� E então? Você acha que podemos sustentar a normalidade dosresíduos com base na observação do QQ plot?

Pàg 48 de 66

APENDICE - SCRIPT R e Dados 28 talhõesTALHAO;DOSE;PROD

1;33.2;476

2;58.8;514

3;30;400

4;52.6;578

5;52.8;512

6;30;284

7;47.4;576

8;30;226

9;30;334

10;41;420

11;61;580

12;19.8;318

13;20.8;248

14;10;244

15;12.6;178

16;5.2;126

17;8;82

18;51.8;430

19;55.8;482

20;24.6;194

21;37.8;280

22;47.6;412

23;46.4;522

24;62.6;652

25;15.8;110

26;41.4;490

27;41.4;358

28;36.6;414

Pàg 49 de 66

SCRIPT R

#---> Obs.: para executar este script copie o comando abaixo e execute-o no ambiente do R

##### source("C:/Estatistica/Cursos/tut002_codR01.R")

#---------------------------------------------------------------------------------------------

# Pegar o diretório atual...

initial.dir<-getwd()

#---------------------------------------------------------------------------------------------

# Mudar para o diretório de trabalho....

setwd("C:/Estatistica/Cursos")

#---------------------------------------------------------------------------------------------

# Carregar eventuais bibliotecas necessárias ao processamento...

library("lmtest")

library("car")

library("HH")

#---------------------------------------------------------------------------------------------

# Definir o arquivo "log" com resultados do processamento...

#sink("result01.out", append = FALSE)

#---------------------------------------------------------------------------------------------

# Ler os dados...

vDados <- read.csv("tut002_dose_prod.csv", sep=";", dec=".", header=TRUE, stringsAsFactors = FALSE)

#############vDados <- read.csv("tut002_dose_prod_29outlier.csv", sep=";", dec=".", header=TRUE,

stringsAsFactors = FALSE)

Pàg 50 de 66

cat("\n \n")

cat("Variáveis disponíveis\n")

cat("-------------------- \n")

print(names(vDados))

cat("\n \n")

#--------------------------------------------------------------------------------------------

cat("\n \n")

cat("Dados da amostra\n")

cat("----------------\n")

print(vDados)

cat("\n \n")

###sink("result09.txt", append = FALSE)

#----------------------------------------------------------

#

#---> aqui executa a regressão linear (com o comando lm)

# ... o resultado fica armazenado na variável vRegrLin

#

vRegrLin <- lm(vDados$PROD ~ vDados$DOSE)

#--------------------------------------------------------------------------------------------

Pàg 51 de 66

cat("\n \n")

cat("Estimadores obtidos via regressão linear\n")

cat("----------------------------------------\n")

print(vRegrLin)

#--------------------------------------------------------------------------------------------

cat("\n \n")

cat("Inferências para os estimadores de Beta0 e Beta1\n")

cat("------------------------------------------------\n")

print(summary(vRegrLin))

#------------------------------------------------------------------------------

# Mostrar os resíduos...

#

vDataFrameAux = cbind(vDados,vRegrLin$fitted.values, vRegrLin$residuals)

cat("\n \n")

cat("Dados completos - mostrando resíduos\n")

cat("------------------------------------\n")

print(vDataFrameAux)

#------------------------------------------------------------------------------

# Mostrar variável preditora X resíduos...

#

vDataFrameAux = cbind(vDados$DOSE, vRegrLin$residuals)

vDataFrameAux = as.data.frame(vDataFrameAux)

Pàg 52 de 66

names(vDataFrameAux)[names(vDataFrameAux)=="V1"] <- "DOSE"

names(vDataFrameAux)[names(vDataFrameAux)=="V2"] <- "RESÍDUO"

cat("\n \n")

cat("Variável preditora X Resíduo\n")

cat("----------------------------\n")

print(vDataFrameAux)

#--------------------------------------------------------------------------------------------

cat("\n \n")

cat("Distância de Cook\n")

cat("-----------------\n")

vDistCook = cooks.distance(vRegrLin)

vDataFrameAux = cbind(vDados$TALHAO,vDistCook)

vDataFrameAux = as.data.frame(vDataFrameAux)

names(vDataFrameAux)[1] <- "TALHÃO"

names(vDataFrameAux)[2] <- "DIST_COOK"

print(vDataFrameAux)

cat("\n \n")

#--------------------------------------------------------------------------------------------

cat("\n \n")

cat("Hat-values\n")

cat("----------\n")

Pàg 53 de 66

vHatValues = hatvalues(vRegrLin)

vMediaHatValues = mean(vHatValues)

cat("\n")

cat("Valor médio: hmedio = ")

cat(vMediaHatValues)

cat("\n\n")

#

#---> Hat values que requerem atenção... maior que o dobro da média...

vStatusHatValues = c()

for (i in 1:length(vHatValues)) {

if (vHatValues[i] > (2 * vMediaHatValues)) {

vStatusHatValues[i] = "***"

} else {

vStatusHatValues[i] = ""

}

}

vDataFrameAux = cbind(vDados$TALHAO,vHatValues,vStatusHatValues)

vDataFrameAux = as.data.frame(vDataFrameAux)

names(vDataFrameAux)[1] <- "TALHÃO"

names(vDataFrameAux)[2] <- "Hat Value (h)"

names(vDataFrameAux)[3] <- "Obs"

print(vDataFrameAux)

cat("\n \n")

#--------------------------------------------------------------------------------------------

Pàg 54 de 66

#

f_calcula_MSE <- function(aVetNumerico) {

vSomaQuad = sum(aVetNumerico * aVetNumerico)

vMSE = vSomaQuad / (length(aVetNumerico) - 2) # 2 glib

#---

return(vMSE)

}

cat("\n \n")

cat("Resíduos semi-studentizados (semistudentized residuals)\n")

cat("-------------------------------------------------------\n")

cat("\n")

cat("---> Para melhor detectar outliers usa-se os \n")

cat(" resíduos padronizados (STANDARDIZED RESIDUALS) conhecidos por \n")

cat(" INTERNALLY STUDENTIZED RESIDUALS ... \n\n")

vResiduos = vRegrLin$residuals

vMSE = f_calcula_MSE(vResiduos)

vSemiStudentizeedResiduals = vResiduos / sqrt(vMSE)

vDataFrameAux = cbind(vDados$TALHAO, vResiduos, vSemiStudentizeedResiduals)

vDataFrameAux = as.data.frame(vDataFrameAux)

names(vDataFrameAux)[1] <- "TALHAO"

names(vDataFrameAux)[2] <- "Residuo"

names(vDataFrameAux)[3] <- "Residuo semistudentizado"

print(vDataFrameAux)

cat("\n \n")

Pàg 55 de 66

#--------------------------------------------------------------------------------------------

#

# Testes de hipótese para verificação de homocedasticidade

#

cat("\n \n")

cat("Testes de hipótese para verificação de homocedasticidade\n")

cat("--------------------------------------------------------\n")

cat("\n")

cat("1) Teste de Breush-Pagan\n")

cat(" ---------------------\n")

cat("\n")

cat(" H0: ocorre homocedasticidade nos resíduos \n")

cat(" Ha: não ocorre homocedasticidade nos resíduos \n")

vBPtest = bptest(vRegrLin)

print(vBPtest)

cat("\n")

cat("2) NCV Test\n")

cat(" --------\n")

cat("\n")

cat(" H0: ocorre homocedasticidade nos resíduos \n")

cat(" Ha: não ocorre homocedasticidade nos resíduos \n")

cat("\n")

Pàg 56 de 66

vNCVtest = ncvTest(vRegrLin)

print(vNCVtest)

cat("\n")

cat("3) Brown-Forsythe Test \n")

cat(" --------------------\n")

cat("\n")

cat(" H0: ocorre homocedasticidade nos resíduos \n")

cat(" Ha: não ocorre homocedasticidade nos resíduos \n")

cat("\n")

cat(" Obs.: Este teste não depende da normalidade dos erros.\n")

cat(" É um teste robusto quanto à forte não normalidade\n")

cat(" dos dados. \n")

vBrownForsythe = hov(PROD ~ DOSE, data=vDados)

print(vBrownForsythe)

cat("\n")

#----------------------------------------------------------------------------------

# Gráfico dos dados com reta de regressão

#

dev.new() # Or X11()

dev.1 <- as.integer(dev.cur())

dev.set(dev.1)

plot(vDados$DOSE, vDados$PROD, main="Dados com reta de regressão",

Pàg 57 de 66

xlab="Dose (Kg)", ylab="Produção (Kg)", pch=20)

abline(vRegrLin, col="red", lwd=4)

#----------------------------------------------------------------------------------

# Gráficos para verificar a distância de Cook

#

dev.new() # Or X11()

dev.2 <- as.integer(dev.cur())

dev.set(dev.2)

par(mfrow=c(1,2))

#---> Grafico 1...

plot(vDados$DOSE, vDados$PROD, main="Os dados dos 28 talhões", cex.main=0.8,

xlab="Dose (Kg)", ylab="Produção (Kg)", pch=20)

abline(vRegrLin, col="blue", lwd=4)

text(vDados$DOSE, vDados$PROD, labels=vDados$TALHAO, pos=1, cex=0.7, col="red")

#---> Grafico 2...

plot(vRegrLin, which = 4, xaxt="n", main="Distância de Cook \n Regressão com 28 talhões", cex.main=0.8 )

axis(1,vDados$TALHAO, cex.axis=0.6, xlab="Observação (talhão)")

#----------------------------------------------------------------------------------

# Gráficos para verificar alavancagem (hat-values)

Pàg 58 de 66

#

dev.new() # Or X11()

dev.3 <- as.integer(dev.cur())

dev.set(dev.3)

par(mfrow=c(1,2))

#---> Grafico 1...

plot(vDados$DOSE, vDados$PROD, main="Os dados dos 28 talhões", cex.main=0.8,

xlab="Dose (Kg)", ylab="Produção (Kg)", pch=20)

abline(vRegrLin, col="blue", lwd=4)

text(vDados$DOSE, vDados$PROD, labels=vDados$TALHAO, pos=1, cex=0.7, col="red")

#---> Grafico 2...

plot(vHatValues, xaxt="n", main="Alavancagem (leverage) \n Regressão com 28 talhões", cex.main=0.8,

ylab="Hat values (h)", type="h", xlab="observação (talhão)" )

axis(1,vDados$TALHAO, cex.axis=0.6, xlab="Observação (talhão)")

for (i in 1:length(vHatValues)) {

if (vHatValues[i] > (2 * vMediaHatValues)) {

#### points(vDados[i,1], vHatValues[i], col='red', pch=20)

text(vDados[i,1], vHatValues[i], labels=vDados[i,1], pos=3, cex=0.7)

}

}

Pàg 59 de 66

#----------------------------------------------------------------------------------

# Gráficos para verificar linearidade da função de regressão

#

dev.new() # Or X11()

dev.4 <- as.integer(dev.cur())

dev.set(dev.4)

par(mfrow=c(1,2))

#---> Grafico 1...

vResiduos = resid(vRegrLin)

plot(vDados$DOSE, vResiduos, main="Verificar linearidade via resíduos", cex.main=0.8,

xlab="Dose (Kg)", ylab="Resíduos", pch=20)

abline(0,0, col="blue", lwd=4, lty=2)

vX = vDados$DOSE

vY = vResiduos

vCurvaLowess = loess(vY ~ vX)

vXloess_seq = seq(min(vX),max(vX), (max(vX) - min(vX))/1000)

lines(vXloess_seq, predict(vCurvaLowess,vXloess_seq), col='red', lwd=1)

#---> Grafico 2...

plot(vDados$DOSE, vDados$PROD, main="Gráficos dos dados \n para verificar linearidade", cex.main=0.8,

xlab="Dose (Kg)", ylab="Produção (Kg)", pch=20)

abline(vRegrLin, col="red", lwd=4)

Pàg 60 de 66

#----------------------------------------------------------------------------------

# Gráficos para verificar homocedasticidade (constância da variância)

#

dev.new() # Or X11()

dev.5 <- as.integer(dev.cur())

dev.set(dev.5)

par(mfrow=c(1,3))

#---> Grafico 1...

vResiduos = resid(vRegrLin)

plot(vDados$DOSE, vResiduos, main="Serve para \n verificar homocedasticidade", cex.main=0.9,

xlab="Dose (Kg)", ylab="Resíduos", pch=20)

abline(0,0, col="blue", lwd=4, lty=2)

vX = vDados$DOSE

vY = vResiduos

vCurvaLowess = loess(vY ~ vX)

vXloess_seq = seq(min(vX),max(vX), (max(vX) - min(vX))/1000)

lines(vXloess_seq, predict(vCurvaLowess,vXloess_seq), col='red', lwd=1)

#---> Grafico 2...

vResiduos = resid(vRegrLin)

plot(vRegrLin$fitted, vResiduos, main="Serve para \n verificar homocedasticidade", cex.main=0.9,

xlab="Produção ajustada(Kg)", ylab="Resíduos", pch=20)

Pàg 61 de 66

abline(0,0, col="blue", lwd=4, lty=2)

vX = vRegrLin$fitted

vY = vResiduos

vCurvaLowess = loess(vY ~ vX)

vXloess_seq = seq(min(vX),max(vX), (max(vX) - min(vX))/1000)

lines(vXloess_seq, predict(vCurvaLowess,vXloess_seq), col='red', lwd=1)

#---> Grafico 3...

#

# Gráfico Scale-Location

#

#------> Padronização dos resíduos... para auxiliar detecção de homo/hetero...cedasticidade !!!

vResiduos = resid(vRegrLin) # pegar os resíduos...

vResiduos = vResiduos / (summary(vRegrLin)$sigma) # resíduos (semi)studentizados...

vResiduos = abs(vResiduos)

vResiduos = sqrt(vResiduos)

vDataFrameAux = cbind(vRegrLin$fitted.values, vResiduos)

vDataFrameAux = as.data.frame(vDataFrameAux)

names(vDataFrameAux)[1] <- "PROD_ESTIMADA"

names(vDataFrameAux)[2] <- "RAIZ_RESIDUO_PADRONIZADO"

plot(vDataFrameAux$PROD_ESTIMADA, vDataFrameAux$RAIZ_RESIDUO_PADRONIZADO, main="Gráfico

SCALE-LOCATION \n Para verificar homocedasticidade", cex.main=0.9,

xlab="Produção ajustada (Kg)", ylab="Raiz(resíduo-standardized)", pch=20)

vY = vDataFrameAux$RAIZ_RESIDUO_PADRONIZADO

Pàg 62 de 66

vX = vDataFrameAux$PROD_ESTIMADA

vCurvaLowess = loess(vY ~ vX)

vXloess_seq = seq(min(vX),max(vX), (max(vX) - min(vX))/1000)

lines(vXloess_seq, predict(vCurvaLowess,vXloess_seq), col='red', lwd=1)

#----------------------------------------------------------------------------------

# Gráficos para verificar normalidade dos resíduos

#

dev.new() # Or X11()

dev.6 <- as.integer(dev.cur())

dev.set(dev.6)

par(mfrow=c(2,2))

#---> Grafico 1... Box plot dos resíduos

vResiduos = resid(vRegrLin) # pegar os resíduos...

boxplot(vResiduos, main="Resíduos dos 28 talhões", ylab="Resíduo")

#---> Grafico 2... Histograma dos resíduos

vResiduos = resid(vRegrLin) # pegar os resíduos...

vNumClassesSturges = nclass.Sturges(vResiduos)

hist(vResiduos, main="Histograma dos resíduos (default do R)", cex.main=0.9, xlab="Resíduo", ylab="Núm.

talhões")

#---> Grafico 3... Histograma dos resíduos

vResiduos = resid(vRegrLin) # pegar os resíduos...

vNumClassesSturges = nclass.Sturges(vResiduos)

hist(vResiduos, main="Histograma dos resíduos com + ou - classes", cex.main=0.9, xlab="Resíduo",

ylab="Núm. talhões", breaks=3)

Pàg 63 de 66

#---> Grafico 4... Histograma dos resíduos

vResiduos = resid(vRegrLin) # pegar os resíduos...

vNumClassesSturges = nclass.Sturges(vResiduos)

hist(vResiduos, main="Histograma dos resíduos com + ou - classes", cex.main=0.9, xlab="Resíduo",

ylab="Núm. talhões", breaks=10)

#---> Grafico 5... Ramo-e-folhas(steam-and-leaf

cat("\n")

cat("Diagrama ramo-e-folhas (steam-and-leafs)\n")

cat("----------------------------------------\n")

vResiduos = resid(vRegrLin) # pegar os resíduos...

stem(vResiduos)

cat("\n")

#----------------------------------------------------------------------------------

# Gráfico QQ-plot para verificar normalidade dos resíduos

#

dev.new() # Or X11()

dev.7 <- as.integer(dev.cur())

dev.set(dev.7)

par(mfrow=c(2,2))

vResiduos = resid(vRegrLin) # pegar os resíduos...

vResiduos = sort(vResiduos) # tem que estar ordenado...

N = length(vResiduos)

Pàg 64 de 66

n.props = pnorm(vResiduos, mean(vResiduos), sd(vResiduos)) # here I calculate the probabilities

associated

# w/ these data if they came from a normal

# distribution w/ the same mean & SD

# I calculate the proportion of x we've gone through at each point

props = 1:N / (N+1)

n.quantiles = qnorm(props, mean=mean(vResiduos), sd=sd(vResiduos)) # this calculates the quantiles (ie

# z-scores) associated w/ the props

# here I bundle them together

my.data = data.frame(x=vResiduos, props=props,

normal.proportions=n.props,

normal.quantiles=n.quantiles)

vMedia = mean(my.data$normal.quantiles)

vDesvPad = sd(my.data$normal.quantiles)

vQuantisNormaisEscoreZ = (my.data$normal.quantiles - vMedia) / vDesvPad

my.data = cbind(my.data,vQuantisNormaisEscoreZ)

names(my.data)[names(my.data)=="vQuantisNormaisEscoreZ"] <- "QuantisTeoricos"

vDadosGrafico = round(my.data, digits=3)

cat("\n")

cat("Dados para composição 'manual' do QQ plot \n")

cat("----------------------------------------- \n")

cat("\n")

print(vDadosGrafico)

cat("\n")

Pàg 65 de 66

#---> Plot PP-plot

plot(my.data$props,my.data$normal.props, main="PP plot", xlab="Proporção dados %", ylab="proporção

normalidade %")

#---> Plot QQ-plot

plot(my.data$QuantisTeoricos,my.data$x, main="QQ plot", xlab="quantis teóricos", ylab="resíduos")

#---> QQ plot do R...

qqnorm(vResiduos)

qqline(vResiduos,col="red")

#---------------------------------------------------------------------------------------------

# Para finalizar...

cat("\n")

cat("-------------------- \n")

cat("FIM DE PROCESSAMENTO \n")

cat("-------------------- \n")

# fechar o arquivo de saida

#sink()

# descarregar as bibliotecas eventualemnte usadas...

detach("package:MASS")

detach("package:tseries")

# volta a apontar ao diretório inicial

setwd(initial.dir)

Pàg 66 de 66

#====================================== F I M

===============================================