\documentclass[a4paper,11pt]{article}
\voffset=-2.54cm%Zera margemsuperior
\hoffset=-2.54cm%zera margem esquerda
\topmargin=0.5cm%margem superior
\textheight=25cm%%Altura do texto
\textwidth=17cm%largura do texto
\usepackage[brazilian]{babel}
\usepackage[ansinew]{inputenc}
\pagenumbering {arabic}
\usepackage{amsthm,amsfonts,bm}
\usepackage{graphicx}
\usepackage[T1]{fontenc}
\usepackage{ae}
\usepackage{Sweave1}






             
\title{PLANEJAMENTO E ANÁLISE DE EXPERIMENTOS}
\author{Éder David Borges da Silva, Renato de Oliveira Gonçalves}

\begin{document}
\maketitle

% Definindo o Sumário
\newpage
\thispagestyle{empty}
\begin{center}
\renewcommand{\contentsname}{SUMÁRIO}%\\
%\begin{flushright} {\normalsize pg} \end{flushright}}
\tableofcontents
\end{center}

\newpage
\section{Introdução ao R}
\begin{Schunk}
\begin{Sinput}
 ########## Introdução ao ambiente de programação no software R ################
 # O material deste arquivo foi baseado nas apostilas sobre o software, R básico
 # de autoria do professor Adílson dos Anjos e do professor Paulo Justiniano 
 
 # o "R" trabalha com programação orientada à objetos
 # para entender melhor podemos executar o comando abaixo
 a <- 2
 # o comando significa que o objeto "a" recebe o valor "2"
 a
 # digitando o nome do objeto é mostrado o valor que foi armazenado
 # o "R" é sensível a letras maiúsculas e minúsculas
 A <- 3
 a
 A
 # os objetos podem receber nomes maiores mas que não podem começar
 # com números ou com "."
 
 # comandos podem ser separados em uma mesma linha por ";"
 x <- 2;y <- 3;x;y
 # se um comando não for completado o sinal de "+" aparecerá
 dados <- c(2,3
 #se isso ocorrer, digite o comando que está faltando ou "ESC" para cancelar
 )
 dados
 #o uso do método recursivo possibilita a execução de expressões dentro
 # de expressões
 a <- 6
 c <- a/(b <- 2);c
 # para listar os objetos armazenados dentro do diretório de trabalho
 # utilaza-se o comando abaixo
 ls()
 #para remover os objetos armazenados utiliza-se o comando abaixo
 #rm(list=ls())
 #operações matemáticas simples podem ser executadas na  linha de comando
 2+2
 2*3
 4/2
 2^2
 #algumas funções matemáticas estão prontas para o uso
 sqrt(4)
 sin(x)
 #algumas aplicações com vetores (vetores serão discutidos com mais detalhes)
 a <- c(2,4,6)
 b <- c(1,2,3)
 a+b
 a/b
 x <- 2
 y <- c(4,8,10)
 y*x
 y/x 
 #outros comandos sobre vetores
 y
 sum(y)#executa a sima do vetor
 length(y)
 media <- sum(y)/length(y)
 media
 mean(y)
 var(y)
 sd(y) 
 # pode-se trabalhar com operadores lógicos definidos de acordo com seguintes
 # símbolos
 # <: menor
 # <=: menor igual
 # >: maior
 # >=: maior igual
 # ==: igual
 # !=: diferente
 # &: "e"
 # ||: "ou"
 
 a <- c(1,2,3,4,5)
 a>3
 a!=2
 b <- c(1,2,3,5,5)
 (b>2) & (b==5)
 b
 which(b==3)
 # existem várias formas de se obter ajuda no "R", por exemplo:
 #help(anova)
 #?anova
 # ainda é possível buscar ajuda das seguintes formas
 #help.search("anova")
 #help.search('anova')
 # para a visualização de exemplos pode-se usar:
 #example(t.test)
 
 # existem de pacotes demonstração que apresentam as funcionalidades do software
 #demo(graphics)
 # digitando o comando abaixo pareceráu uma lista das demos disponíveis
 #demo()
 
 # os tipos de objeto mais importantes para o nosso trabalho são:
 # Vetores;
 # Matrizes;
 # Fatores;
 # Listas;
 # Data frames;
 # Funções;
 x <- c(1,2,3,4,5)
 x
 is.numeric(x)
 class(x)
 # objetos podem ser convertidos com comando:
 x <- as.factor(x)
 class(x)
 is.numeric(x)
 # o "R" possui várias formas para gerar sequencias de números
 a <- 1:20
 a
 b <- 1:20
 b
 # uma sequencia mais elaborada pode ser construída usando a função "seq"
 seq(0,10, by=0.2)
 seq(length=12, from=0, by=3)
 seq(0,20,length=12)
 # usando a função "rep"
 rep(c(1,2,3),2)
 rep(c(1,2,3), each=2)
 rep(1:3, each=3, times=2)
 rep("r", 5)
 # a função "rep" gera um vetor numérico ou alfanumérico
 # para ser utlizado como um fator em um experimento é necessário convertê-lo
 # em um fator
 # pode-se utilizar o comando "gl" que cria um vetor da classe "factor"
 
 A <- gl(4,2)
 A
 class(A)
 levels(A)
 gl(2,4)
 gl(3,3,length=18)
 # para cada elemento de um vetor é possível atribuir um nome
 notas <- c(7,5,8,5,9)
 names(notas) <- c("Joao","Pedro","Paulo","Jose","Tiago")
 notas
 aprovados <- notas>=7
 aprovados
 ###--------------------------------------------------------------------------###
 ############################## GRAFICOS ########################################
 #demo(graphics)
 #demo(image)
 #demo(persp)
 
 # plotando um gráfico simples
 x <- 1:10
 y <- sqrt(x)
 plot(x,y)
 # pode-se acrescentar ao gráfico textos e outros elementos
 plot(x,y,main="Titulo")
 plot(x,y,main="Titulo",sub="sub titulo\n")
 plot(x,y,main="Titulo",xlab="Eixo x",ylab="Eixo y")
 plot(x,y,main="Titulo",xlab="Eixo x",ylab="Eixo y",type="l")
 plot(x,y,main="Titulo",xlab="Eixo x",ylab="Eixo y",type="l",col="red")
 plot(x,y,main="Titulo",xlab="Eixo x",ylab="Eixo y",type="l",col=2, axes=FALSE)
 plot(x,y,main="Titulo",xlab="Eixo x",ylab="Eixo y",type="b",col=2)
 # pode-se acrescentar pontos e linhas em um gráfico
 plot(x,y)
 points(rev(x),y)
 lines(x,3-y)
 lines(x,4-y)
 lines(x,5-y)
 # pode-se adicionar linhas com comando abline
 plot(x,y)
 abline(h=c(3,2))
 abline(v=4)
 # o argumento "pch" define o símbolo que será utilizado no gráfico
 x<-1:20
 y<-1:20
 plot(x,y,pch=1:20)
 # o comando "cex=" aumenta ou diminui o tamanho dos caracteres gráficos
 plot(x,y,pch='E',cex=1.5)
 plot(x,y,pch=1:20,cex=0.5)
 # identificar um ponto no gráfico
 plot(x,y)
 #locator(n=1)
 
 # identificar vários pontos no gráfico
 #identify(x,y)
 
 # inserir uma legenda no gráfico
 plot(x,y)
 legend(x=8,y=2,legend="pontos",pch=1,cex=1.5)
 plot(x,y,type="l")
 lines(x,y+1,col=2)
 legend(8,2,c("Grã 1","linha 2"),pch=19 ,col=1:2,cex=1.5)
 # inserir um texto no meio do gráfico
 text(locator(2),"Ponto")
 pdf("Meugrafico2.pdf")
 plot(x,y,type="l")
 lines(x,y+1,col=2)
 legend(8,2,c("Grã 1","linha 2"),pch=19 ,col=1:2,cex=1.5)
 dev.off()
 # para dividir a janela gráfica usa-se o comando "par()
 par(mfrow=c(2,2))
 plot(x,y)
 plot(x,y)
 # para retornar a janela gráfica à configuração normal basta fechá-la
 
 ###--------------------------------------------------------------------------###
 ############################### VETORES #######################################
 #o objeto mais simples do "R" é um vetor
 x.n<-c(89,65,76.5)
 x.l<-c('a','b','c')
 x.l<-letters[1:3]
 x.l
 x.L<-LETTERS[4:5]
 x.L
 # outras formas de secriar um vetor
 assign ("x.n",c(89,65,76.5))
 c('a','b','c')->x.l
 # pode-se criar um vetor utilizando o comando "scan"
 dados<-scan()
 dados
 # pode-se trabalhar com os elementos dentro do vetor
 x<-1:10
 x
 y<-x[5:10]
 x; y
 #retira a terceira observação de y
 z<-y[-3];z
 z
 z[1]<-NA
 z 
 mean(na.omit(z)[-1])
  a<-1:5
 c<-a[a>2]
 c
 ###--------------------------------------------------------------------------###
 ########################## LISTAS E VETORES ####################################
 
 # uma lista é uma coleção de objetos no "R"
 dados<-list(nome="Adilson",civil="casado", cachorros=2)
 dados
 # os objetos de uma lista não precisam ser da mesma classe ou tamanho
 dados[3]
 dados[[3]]
 dados$cachorro
 # a função "length" fornece o número de objetos na lista ou o tamanho da lista
 length(dados)
 # pode-se digitar apenas um número mínimo de letras para identificar o objeto
 dados$c
 dados$ca
 # para alterar os valores da lista
 dados[[3]]<-1
 dados$cachorros<-1
 dados$profissao<-"professor"
 dados
 ###--------------------------------------------------------------------------###
 ############################### MATRIZES #######################################
 # criando uma matriz
 m1 <- matrix(1:12, ncol = 3)
 m1
 matrix(1:12, ncol = 3, byrow = T)
 # obtendo alguns parâmetros da matriz
 length(m1)
 dim(m1)
 nrow(m1)
 ncol(m1)
 # Extraindo alguns valores da matriz
 m1[1, 2]
 m1[2, 2]
 m1[, 2]
 m1[3, ]
 m1[1:2, 2:3]
 #atribuindo nomes as dimensões da matriz
 dimnames(m1)
 dimnames(m1) <- list(c("L1", "L2", "L3", "L4"), c("C1", "C2", "C3"))
 dimnames(m1)
 m1[c("L1", "L3"), ]
 m1[c(1, 3), ]
 m2 <- cbind(1:5, 6:10)
 m2
 m3 <- cbind(1:5, 6)
 m3
 # operações com matrizes
 m4 <- matrix(1:6, nc = 3)
 m5 <- matrix(10 * (1:6), nc = 3)
 m4
 m5
 m4*m5
 m4+m5
 m5-m4
 m4/m5
 # a multiplicação de matrizes é feita usando o operador %*%. A funçãoo t() 
 #faz transposição e a inversão é obtida com solve().
 t(m4) %*% m5
 #a função solve() fornece a solução de sistemas de equações lineares.
 # x + 3y - z = 10
 # 5x - 2y +z = 15
 # 2x + y - z = 7
 
 mat <- matrix(c(1, 5, 2, 3, -2, 1, -1, 1, -1), nc = 3)
 vec <- c(10, 15, 7)
 solve(mat, vec)
 ###--------------------------------------------------------------------------###
 ################################# DATA FRAMES ##################################
 
 d1 <- data.frame(X = 1:10, Y = c(51, 54, 61, 67, 68, 75, 77, 75, 80, 82))
 d1
 names(d1)
 d1$X
 d1$Y
 plot(d1)
 plot(d1$X, d1$Y)
 d2 <- data.frame(Y = c(10 + rnorm(5, sd = 2), 16 + rnorm(5, sd = 2), 14 + rnorm(5, sd = 2)))
 d2
 d2$lev <- gl(3, 5)
 d2
 by(d2$Y, d2$lev, summary)
 d3 <- expand.grid(1:3, 4:5)
 d3
 expand.grid(1:3, 1:2)
 #Listas
 
 lis1 <- list(nomes = c("Pedro", "Joao", "Maria"), mat = matrix(1:6, nc = 2))
 lis1
\end{Sinput}
\end{Schunk}
\newpage
\section{Revisão de estatistica basica}
\begin{Schunk}
\begin{Sinput}
 ###--------------------------------------------------------------------------###
 library(fBasics)
 require(agricolae)
 require(MASS)
 ###--------------------------------------------------------------------------###
 # Estatistica descitiva
 set.seed(12)
 x <- rnorm(15,m=900,sd=50)
 x      
 sum(x)
 min(x)
 max(x)
 length(x)
 mean(x)
 median(x)
 range(x)
 var(x)
 sd(x)
 summary(x)
 quantile(x)
 quantile(x,probs = seq(0, 1, 0.05))
 ###--------------------------------------------------------------------------###
 # Distribuições de probabilidades
 # Gerando uma amostra aleatorio normal
 set.seed(123)
 x <- rnorm(100)
 par(mfrow=c(2,1))
 hist(x,prob=TRUE)
 lines(density(x))
 curve(dnorm(x),from=min(x),to=max(x),col='red',add=TRUE,lwd=2)
 boxplot(x, horizontal =TRUE)
 # Grafico QQPlot
 qqnorm(x)
 qqline(x)
 # Exemplo
 # Uma variedade de tomate possui produtividade media de 7,9 kg/planta  e variância de
 # 0,97 (kg/planta). Supor que a distribuição seja normal e calcular a probabilidade
 # produtividade (X) de uma planta sorteada dessa variedade estar de acordo com os
 # seguintes eventos:
 
 # a) X > 9,0 kg
 1-pnorm(9,m=7.9,sd=0.97^0.5)
 # b) 8,0 < X < 9,5 kg
 pnorm(9.5,m=7.9,sd=0.97^0.5)-pnorm(8,m=7.9,sd=0.97^0.5)
 # c) X < 7 kg
 pnorm(7,m=7.9,sd=0.97^0.5)
 # d) 6,5 < X < 8,5 kg
 pnorm(8.5,m=7.9,sd=0.97^0.5)-pnorm(6.5,m=7.9,sd=0.97^0.5)
 # Intervalo com 95%
 qnorm(c(0.025,0.975),m=7.9,sd=0.97^0.5)
 # Gerando uma amostra aleatorio Poisson
 set.seed(123)
 x <- rpois(50,l=2)
 par(mfrow=c(2,1))
 hist(x,prob=TRUE)
 lines(density(x))
 lines(seq(0,6,by=1),dpois(seq(0,6,by=1),lam=2),col='red',lwd=2)
 boxplot(x, horizontal =TRUE)
 # Exemplo
 # Uma raça bovina apresenta uma doença virotica em 2% dos animais;
 # a) Qual é a probabilidade de uma amostra de 100 animais não apresentar animal doentes?
 ppois(0,lam=2)
 # b) Qual é a probabilidade de se obterem, no maximo, 1 animais doentes nessa amostra?
 ppois(1,lam=2)
 # b) Qual é a probabilidade de se obterem, no maximo, 2 animais doentes nessa amostra?
 ppois(2,lam=2)
 # b) Qual é a probabilidade de se obterem, 5 animais doentes nessa amostra?
 ppois(5,lam=2)-ppois(4,lam=2)
 # Gerando uma amostra aleatorio binomial
 set.seed(123)
 x <- rbinom(100,50,0.5)
 par(mfrow=c(2,1))
 hist(x,prob=TRUE)
 lines(density(x))
 lines(seq(0,50,by=1),dbinom(seq(0,50,by=1),50,0.5),col='red',lwd=2)
 boxplot(x, horizontal =TRUE)
 # Exemplo
 # A probabilidade do sucesso de uma semente germinar é p=0.90. Determine a probabilidade de
 # germinação, em uma amostra de tamanho n=4, de todas as possibilidades de realizações de uma
 # variavel aleatoria X, defina o numero de sementes germinadas na amostra.
 
 # X = 0
 pbinom(0,4,0.90)
 # X = 1
 pbinom(1,4,0.90)-pbinom(0,4,0.90)
 # X = 2
 pbinom(2,4,0.90)-pbinom(1,4,0.90)
 # X = 3
 pbinom(3,4,0.90)-pbinom(2,4,0.90)
 # X = 4
 pbinom(4,4,0.90)-pbinom(3,4,0.90)
 ###--------------------------------------------------------------------------### 
 # Explorando as distribuições de probabilidade
 set.seed(123)
 x <- rnorm(100)
 x
 FX <- graph.freq(x)
 table.freq(FX)
 M <- mean(x)
 M
 Sd <- sd(x)
 Sd
 hist(x,prob=TRUE,ylim=c(0,1))
 xx <- seq(min(x),max(x),l=200)
 lines(xx,dnorm(xx,M,Sd),lwd=2,col="red")
 plot.ecdf(x,add=TRUE)
 ###--------------------------------------------------------------------------###  
 ### Intervalos de confiança
 
 IC <- function(M=M,Sd=Sd,SIC=0.95)
 {
 hist(x,prob=TRUE,main=SIC)
 xx <- seq(min(x),max(x),l=200)
 lines(xx,dnorm(xx,M,Sd),lwd=3)
 xI <- seq(qnorm((1-SIC)/2,M,Sd),qnorm((1-(1-SIC)/2),M,Sd),l=1000)
 xY <- dnorm(xI,M,Sd)
 xY[c(1,1000)] <- 0
 
 xI1 <- seq(qnorm(0.00000001,M,Sd),qnorm((1-SIC)/2,M,Sd),l=1000)
 xY1 <- dnorm(xI1,M,Sd)
 xY1[c(1,1000)] <- 0
 
 xI2 <- seq(qnorm((1-(1-SIC)/2),M,Sd),qnorm(0.99999999,M,Sd),l=1000)
 xY2 <- dnorm(xI2,M,Sd)
 xY2[c(1,1000)] <- 0
 
 polygon(xI,xY,col='gray')
 polygon(xI1,xY1,col='lightblue')
 polygon(xI2,xY2,col='lightblue')
 }
 ### Intervalo de confiança
 IC(M=M,Sd=Sd,SIC=0.95)
 IC(M=M,Sd=Sd,SIC=0.85)
 for ( i in seq(0.001,0.999999,by=0.01))
  {
      IC(M=M,Sd=Sd,SIC=i)
       #Sys.sleep(0.2)
  }
 # Ajustando paramentros via maxima verossimilhança
 # Teoria do metodo
 
 x <- rnorm(1500,456.25,125)
 fitdistr(x,'normal')
 mean(x)
 sd(x)
 ###--------------------------------------------------------------------------### 
 ## Gerando serie de dados ficticios
 n <- 150
 set.seed(123)
 dados <- data.frame(Cultivar=factor(rep(c("A","B"),each=n)),
 Prod=c(rnorm(n,m=9500,sd=300),rnorm(n,m=8000,sd=500)))
 summary(dados)
 with(dados,by(Prod,Cultivar,summary))
 boxplot(dados$Prod~dados$Cultivar)
 ### Histograma dos dados
 basicStats(dados$Prod, ci = 0.95)
 histPlot(as.timeSeries(dados$Prod))
 ### Plotando a distribuição das duas cultivares
 hist(dados$Prod[dados$Cultivar=="A"],prob=TRUE,col='lightblue',
 xlim=c(min(dados$Prod),max(dados$Prod)),main='Cultivar A e B',xlab='Produtividade')
 lines(density(dados$Prod[dados$Cultivar=="A"]),lwd=3)
 hist(dados$Prod[dados$Cultivar=="B"],prob=TRUE,col='blue',add=TRUE)
 lines(density(dados$Prod[dados$Cultivar=="B"]),lwd=3)
 ## Igualdade de Variancias
 var.test(dados$Prod~dados$Cultivar,conf.level=0.95)
 ## Igualdade de medias
 t.test(dados$Prod~dados$Cultivar,var.equal = TRUE,conf.level=0.95)
 model.av <- aov(Prod~Cultivar,data=dados)
 summary(model.av)
 ###--------------------------------------------------------------------------### 
 # Amostragem
 # GERA UMA POPULAÇÃO DE APARTIR DESTA FAZ AMOSTRAGENS NESTAS DE TAMANHO 1,2,3,...
 # ATE O TAMANHO DA AMOSTRA, DEPOIS CALCULA o COMPORTAMENTO DA MÉDIA E VARIANCIA DESTAS
 # AMOSTRAS E COMPARA GRAFICAMENTE COM MEDIA E VARIANCIA DA POPULAÇÃO
 rm(list=ls())
 set.seed(123)##Para Gerar sempre a mesma população, podendo ser alterado
 ###Parametros da Simulação
 tpop<-100 #tamanho da Poupulação a ser amostrada
 mpop<-125.5  #Media da população
 vpop<-12  #Variância da População
 ###Grafico com amostras
 #amostra1<-5  #Grafico com XXX amostras
 #amostra2<-30 #Grafico com XXX amostras
 ###Calculos
 populacao<-rnorm(tpop*20,mpop,sqrt(vpop))
 medpop<-mean(populacao)
 varpop<-var(populacao)
 tabela<-cbind(namostra=1:tpop,mean=c(rep(NA,tpop)),var=c(rep(NA,tpop)))
 for (i in 2:tpop)
 	{
 		amostra<-sample(populacao,i)
 		tabela[i,2]=mean(amostra)
 		tabela[i,3]=var(amostra)
                 par(mfrow=c(3,1))
 		plot(function(x) dnorm(x,medpop,sqrt(varpop)),xlim=c(mpop-3*sqrt(vpop),
     mpop+3*sqrt(vpop)),ylim=c(0,0.35),lwd=3,
     main="Distribuição População(Preto) x Amostra(Vermelho)",
     xlab=c("Tamanho de amostra",tabela[i,1]),ylab="Frequencia")
 		plot(function(x) dnorm(x,tabela[i,2],sqrt(tabela[i,3])),
     xlim=c(mpop-3*sqrt(vpop),
     mpop+3*sqrt(vpop)),col="red",lwd=2,add=T)
 		hist(amostra,prob=T,nclass=80,add=T)
 		rug(amostra)
     plot(tabela[,1],tabela[,2],main="Media",ylab=mean(populacao),
     xlab="Tamanho da Amostra")
     abline(h=medpop,lwd=2,col='red')
     plot(tabela[,1],tabela[,3],main="Variancia",ylab=var(populacao),
     xlab="Tamanho da Amostra")
     abline(h=varpop,lwd=2,col='red')
 	}
 par(mfrow=c(2,2))
 plot(tabela[,1],tabela[,2],main="Media",ylab=mean(populacao),
 xlab="Tamanho da Amostra")
 abline(h=medpop,lwd=2,col='red')
 plot(tabela[,1],tabela[,3],main="Variancia",ylab=var(populacao),
 xlab="Tamanho da Amostra")
 abline(h=varpop,lwd=2,col='red')
 hist(tabela[,2],main="Media",xlab=c("Média da População",mean(populacao)),
 nclass=100,prob=TRUE)
 lines(density(tabela[,2],na.rm=TRUE),col="blue",lwd=2)
 rug(tabela[,2])
 hist(tabela[,3],main="Variância",xlab=c("Variância da População",var(populacao)),
 nclass=100,prob=TRUE)
 lines(density(tabela[,3],na.rm=TRUE),col="blue",lwd=2)
 rug(tabela[,3])
\end{Sinput}
\end{Schunk}
\newpage
\section{Delineamento inteiramente ao acaso}
\begin{Schunk}
\begin{Sinput}
 ################## Delineamento Completamente Casualizado ######################
 # Os dados foram retirados do livro Estatística aplicada à pesquisa agrícola
 # do autor Francisco Zimmermann, pg.:56.
 # São referentes ao experimento do pesquisador Albert Baêta Santos sobre adubação
 # nitrogenada em arroz irrigado.
 # Tratamentos:
 # 1) 80 kg/ha no plantio;
 # 2) 40 kg/ha no plantio e 40 kg/ha aos 40 DAE;
 # 3) 13,2 kg/ha no plantio e 66,8 kg/ha aos 40 DAE;
 # 4) 13,2 kg/ha no plantio e 33,4 kg/ha aos 40 e aos 60 DAE.
 ###---------   Análise do Experimento   ----------###
 
 # Relizando a limpeza de todas os dados armazenados na memoria do programa
 rm(list=ls())
 # Extraindo os dados da planilha excel
 require(RODBC)
 xlscon <- odbcConnectExcel("exemplo dic.xls")
 dados <- sqlFetch(xlscon, "dados")
 odbcClose(xlscon)
 head(dados)  
 # inspecionando o classe de dados
 is.data.frame(dados)
 #Verificando as variaveis contida no conjunto de dados
 names(dados)
 # Verificando o conteudo das variaveis
 dados$resp
 dados$trat
 # Verificando a classe de cada variaveis
 is.factor(dados$trat)
 is.numeric(dados$resp)
 # Lembramos que resp deve ter conteudo numerico e trat deve ser fator
 dados$trat<-as.factor(dados$trat)
 dados$resp<-as.numeric(dados$resp)
 # Com o comando "?summary" obtemos uma analise descritivva dos dados
 summary(dados)
 # Com a função "?by", podemos examinar os dados de uma forma mais eficiente
 sumarygeral<-by(dados$resp,dados$trat, summary)
 sumarygeral
 # Podemos também extrair algumas de forma isolada
 by(dados$resp,dados$trat,mean)
 by(dados$resp,dados$trat, var)
 by(dados$resp,dados$trat, sd)
 by(dados$resp,dados$trat, sum)
 by(dados$resp,dados$trat, median)
 by(dados$resp,dados$trat, range)
 by(dados$resp,dados$trat, min)
 by(dados$resp,dados$trat, max)
 # Para facilitar a manipulação dos dados podemos attachar o objeto "dados" que contem
 # nossa variaveis de interesse, após este comando as variaveis ficaram na memoria do R
 attach(dados)
 resp
 # Para realizar o inverso do attach utilizamos a função detach que tira o objeto da memoria
 # detach(dados)
 # A inspeção das variveis contidas na memoria do R pode ser feita com o comando "search"
 search()
 # Outra forma similar ao comando "by" e utilizar o c comando "?tapply"
 dados.m <-tapply(resp, trat, mean)
 dados.m
 dados.v <-tapply(resp, trat, var)
 dados.v
 dados.sd <-tapply(resp, trat, sd)
 dados.sd
 # Como inspeção grafica dos dados podemos utilizar o comando "?plot" ou "?boxplot"
 plot(dados,main="Producao de arroz em funcao da adubacao",
 xlab= expression('Tratamentos'),ylab="Producao (kg/ha)")
 # Com a função points podemos colocar a media sobre o boxplot, outra função similar
 # e interessante e a função "?lines"
 points(dados.m, pch="x", col="blue", cex=1.5)
 # Da mesma maneira como foi colocado a medias podemos acresentar o desvio padrão
 points(dados.m+dados.sd, pch="-", col=3, cex=2.5)
 points(dados.m-dados.sd, pch="-", col=3, cex=2.5)
 # A ANOVA e relizada com o comando "?aov" detalhes sobre a especificação de modelos
 # serão feitos durante o curso conforme estes aparecerão
 dados.av <- aov(resp ~ trat, data = dados)
 dados.av
 ## Extraindo CV
 require(agricolae)
 cv.model(dados.av)
 # A obtenção do quadro tipico da ANOVA a qual conheçemos pode ser feito com o
 # comando "?summary" ou "?anova"
 summary(dados.av)
 anova(dados.av)
 # Análise dos pressupostos da ANOVA
  #análise gráfica
 par(mfrow=c(2,2))
 plot(dados.av)
 # Teste de Bartllet para homocedasticidade
 bartlett.test(dados.av$res, trat)
 # Teste de Shapiro-Wilk para Normalidade
 shapiro.test(dados.av$res)
 # Teste de Tukey para comparações múltiplas
 # No R, o teste de Tukey é apresentado através de intervalos de confiança.
 # A interpretação é: se o intervalo de confiança para a diferença entre duas
 # médias não incluir o valor zero, significa que se rejeita a hipótese nula,
 # caso contrário, não se rejeita. O resultado pode ser visto através de uma 
 # tabela e/ou graficamente:
 dados.tu <- TukeyHSD(dados.av)
 dados.tu
 plot(dados.tu)
 # HSD.test
 # O pacote agricolae tem a função de comparação atraves do Tukey implementada com
 # uma saida mais interessante, o "HSD.test", nesta função fornecemos os tratamento,
 # os graus de liberdade e soma de quadrados do erro e significancia.
 require(agricolae)
 dados.anova<-anova(dados.av)
 names(dados.anova)
 tuk<-HSD.test(resp,trat,dados.anova$Df[2],dados.anova$Mean[2],alpha=0.05)
 tuk
 # Grafico implementados no pacote agricolae
 bar.group(tuk,border="blue",col="green",ylim=c(0,8000),main="Meu grafico",
 xlab="tratamento",y="media")
 bar.err(tuk,border="blue",col="green",ylim=c(0,8000),main="Meu grafico",
 xlab="tratamento",y="media",add=T)
 # Compração via Scott-Knott
 require(ScottKnott)
 ?SK
 dados.SK<- SK(trat,resp, model='resp ~ trat', which='trat',error='Within',
 sig.level=0.05)
 summary(dados.SK)
 # Comparando resultados com o Tukey
 tuk
 # Comparação via contrastes
 # OBS: Explicar teoria e tecnica dos Balões
 cont.dados<-matrix(c(3, -1, -1, -1,#1° Contraste
                      0, 2, -1, -1,#2° Contraste
                      0, 0, 1, -1) #3° Contraste
                      ,nrow=4,ncol=3,byrow=F)
 cont.dados
 # Definindo os contrastes
 contrasts(dados$trat)<-cont.dados
 contrasts(dados$trat)
 dados$trat
 # Executando a ANOVA atraves de contrastes
 dados.avC <- aov(resp ~ trat, data = dados)
 effects(dados.avC)
 summary(dados.avC)
 summary(dados.av)
 # Contrastes estabelecidos
 dados.avC$con
 dados.av$con
 summary(dados.avC,split=list(trat=list(C1=1,C2=2,C3=3)))
\end{Sinput}
\end{Schunk}
\newpage
\section{Delineamento Blocos ao Acaso}
\begin{Schunk}
\begin{Sinput}
 ###-------------------------------------------------------------------------####
 ################################## DBA #########################################
 # Delineamento Blocos ao acaso
 # Os dados foram retirados do livro Estatística aplicada à pesquisa agrícola
 # do autor Francisco Zimmermann, pg.: 76.
 # São referentes ao ensaio regional de feijão preto
 # Tratamentos:
 #  1) 12 cultivares de feijão preto
 # limpando os dados armazenados
 rm(list=ls())
 #Extraindo os dados da planilha excel
 require(RODBC)
 xlscon <- odbcConnectExcel("exemplo dba.xls")
 dados <- sqlFetch(xlscon, "dados")
 odbcClose(xlscon)
 head(dados)
 #explorando os dados
 summary(dados)
 dados$bloco<-as.factor(dados$bloco)
 dados$trat<-as.factor(dados$trat)
 summary(dados)
 summary(dados$resp)
 require(fBasics)
 basicStats(dados$resp, ci = 0.95)
 histPlot(as.timeSeries(dados$resp))
 densityPlot(as.timeSeries(dados$resp))
 qqnormPlot(as.timeSeries(dados$resp))
 par(mfrow = c(1,2))
 plot.design(dados,xlab="fatores", ylab="medias dos genótipos",)
 attach(dados)
 search()
 dados.m <-tapply(resp, trat, mean)
 dados.m
 dados.v <-tapply(resp, trat, var)
 dados.v
 boxplot(resp~trat)
 points(dados.m, pch="x", col=2, cex=1.5)
 # Executando a analise de variância
 dados.av<-aov(resp~bloco+trat)
 summary(dados.av)
 # ANOVA não considerando os blocos - blocos não significativos
 dados.avSB<-aov(resp~trat)
 summary(dados.avSB)
 # Extraindo CV
 require(agricolae)
 cv.model(dados.av)
 # Verificando os residuos
 par(mfrow=c(2,2))
 plot(dados.av)
 # Homocedaciade
 bartlett.test(dados.av$res~trat)
 # normalidade 
 shapiro.test(dados.av$res)
 # Teste Tukey
 dados.anova<-anova(dados.av)
 dados.anova
 names(dados.anova)
 tuk<-HSD.test(resp,trat,dados.anova$Df[3],dados.anova$Mean[3],alpha=0.05)
  bar.group(tuk,border="blue",col="green",ylim=c(0,3000),main="Meu grafico",
  xlab="tratamento",y="Produtividade media")
 #teste Scott-Knott
 require(ScottKnott)
 dados.SK<- SK(trat,resp, model='resp ~ trat + bloco', which='trat',error='Within',
 sig.level=0.05)
 summary(dados.SK)
 #Comparando resultados com o Tukey
 tuk
\end{Sinput}
\end{Schunk}
\newpage
\section{Fatorial 2x3}
\begin{Schunk}
\begin{Sinput}
 ########################## Esquema fatorial 2x3 ################################
 # O experimento a seguir, refere-se a um estudo de mudas de Eucalipto onde foi 
 # estudado o crescimento (cm) de duas espécies (E1 e E2) em três tipos de recipiente
 # (R1, R2 e R3). O delineamento foi completamente casualizado com 4 repetições. 
 # extraido "Experimentação Agricola" Kronka & Banzatto pg 106.
 rm(list=ls())
 fat2x3 <- read.table("fatorial2x3.txt", header=T)
 names(fat2x3)
 summary(fat2x3)
 attach(fat2x3)
 is.factor(rec)
 is.factor(esp)
 is.factor(resp)
 is.numeric(resp)
 # Explorando os dados
 fat2x3.m<-tapply(resp, list(rec,esp), mean)
 fat2x3.m
 fat2x3.mr<-tapply(resp, rec, mean)
 fat2x3.mr
 fat2x3.me<-tapply(resp, esp, mean)
 fat2x3.me
 # Verificando interações
 par(mfrow=c(1,2))
 interaction.plot(rec, esp, resp)
 interaction.plot(esp, rec, resp)
 # Análise de variância
 fat2x3.av<-aov(resp ~ rec + esp + rec * esp)
 # Resumindo o comando:
 fat2x3.av<-aov(resp ~ rec * esp)
 summary(fat2x3.av)
 # A função a seguir ?model.tables? produz tabelas das médias definidas pelo modelo.
 fat2x3.mt <- model.tables(fat2x3.av, ty="means")
 fat2x3.mt
 ### Análise de resíduos
 par(mfrow=c(2,2))
 plot(fat2x3.av)
 ### NOrmalidade
 hist(fat2x3.av$res)
 shapiro.test(fat2x3.av$res)
 # Box-plot para os tratamentos vs resíduos, se existe homocedastidade, espera-se
 # que os Box-plot sejam semelhantes, ou seja, a variabilidade é a mesma em todas as caixas.
 par(mfrow=c(2,2))
 boxplot(fat2x3.av$res~esp)
 boxplot(fat2x3.av$res~rec)
 boxplot(fat2x3.av$res~esp:rec)
 bartlett.test(fat2x3.av$res, esp)
 bartlett.test(fat2x3.av$res, rec)
 ##----------------------------------------------------------------------------##
 ### Desdobrando interações
 ### Desdobrando os graus de liberdade de um fator dentro de cada nível do outro, 
 # utilizando a notação / que indica efeitos aninhados.
 par(mfrow=c(1,2))
 interaction.plot(rec, esp, resp)
 interaction.plot(esp, rec, resp)
 ###Construindo o quadro auxiliar para compreender o desdobramento
 QuadA<-tapply(resp,list(esp,rec),sum)
 QuadA
 addmargins(QuadA)
 fat2x3.avr <- aov(resp ~ rec/esp)
 #Para que o comando "summary" possa mostrar nossa analise adequadamente e necessario
 # definir os efeitos de interesse.
 effects(fat2x3.avr)
 summary(fat2x3.avr, split=list("rec:esp"=list(esp.dentro.r1=1,
                                               esp.dentro.r2=2,
                                               esp.dentro.r3=3)))
 fat2x3.ave <- aov(resp ~ esp/rec)
 effects(fat2x3.ave)
 summary(fat2x3.ave, split=list("esp:rec"=list(rec.dentro.e1=c(1,3),
                                               rec.dentro.e2=c(2,4))))
 require(agricolae)
 ### Tukey
 fat2x3.anova<-anova(fat2x3.av)
 fat2x3.anova
 ##----------------------------------------------------------------------------##
 #Recipientes dentro de esp1
 rec1<-subset(fat2x3,fat2x3$esp=="e1")
 HSD.test(rec1$resp,rec1$rec,fat2x3.anova$Df[4],fat2x3.anova$Mean[4],alpha=0.05)
 #Recipientes dentro de esp2
 rec2<-subset(fat2x3,fat2x3$esp=="e2")
 HSD.test(rec2$resp,rec2$rec,fat2x3.anova$Df[4],fat2x3.anova$Mean[4],alpha=0.05)
 ##----------------------------------------------------------------------------##
 #Especie dentro de recipiente1
 esp1<-subset(fat2x3,fat2x3$rec=="r1")
 HSD.test(esp1$resp,esp1$esp,fat2x3.anova$Df[4],fat2x3.anova$Mean[4],alpha=0.05)
 #Especie dentro de recipiente2
 esp2<-subset(fat2x3,fat2x3$rec=="r2")
 HSD.test(esp2$resp,esp2$esp,fat2x3.anova$Df[4],fat2x3.anova$Mean[4],alpha=0.05)
 #Especie dentro de recipiente3
 esp3<-subset(fat2x3,fat2x3$rec=="r3")
 HSD.test(esp3$resp,esp3$esp,fat2x3.anova$Df[4],fat2x3.anova$Mean[4],alpha=0.05)
\end{Sinput}
\end{Schunk}
\newpage
\section{Parcelas Subdivididas}
\begin{Schunk}
\begin{Sinput}
 ################################## Parcela sub-dividida ########################
 # O experimento a seguir refere-se a um estudo obtido por meio de um delineamento 
 # em blocos completos casualizados no arranjo de parcelas subdivididas. Neste 
 # experimento estudou-se o efeito de diferentes variedades de aveia (A) e 
 # diferentes tratamentos de sementes (B) sobre o rendimento em kg por parcela. 
 # Portanto, considere o experimento em parcelas subdivididas de dados de produção 
 # de aveia. Lembrando que foi conduzido em um delineamento em blocos completos casualizados.
 subd<- read.table("subdividida.txt", header=T)
 subd
 summary(subd)
 subd$aveia <- as.factor(subd$aveia)
 subd$trat <- as.factor(subd$trat)
 subd$bloco <- as.factor(subd$bloco)
 summary(subd)
 attach(subd)
 # Análise de variância
 subd.av <- aov(resp ~ bloco + aveia*trat + Error(bloco:aveia))
 summary(subd.av)
 # Visualização gráfica do comportamento de variedades dentro de tratamentos de
 # sementes e vice-versa.
 par(mfrow=c(1,2))
 interaction.plot(aveia, trat, resp, xlab="Variedades")
 interaction.plot(trat, aveia, resp,xlab="Trat. de Sementes")
 # Em experimento de parcela subddividida, a soma de quadrado e graus de liberdade 
 # dos erros parciais podem ser obtidos, interpletando simplemente isso como uma 
 # interação, no caso do exemplo acima o "Error(bloco:a)" ou "bloco*a" reproduzem 
 # os mesmos valores, sendo a diferença que a primeira sintese já organiza a tabela
 # da ANOVA de forma correta aplicando o teste F de forma certa, isso não é obtido
 # na segunda sintese. Mas para analisar os pressupostos, é comun trabalhar com o 
 # erro principal onde a secunda sintese facilita a extração deste, exemplificado 
 # no comando a seguir:
 subd.avb <- aov(resp ~ bloco + aveia*trat + bloco*aveia)
 summary(subd.avb)
 ### As médias de cada fator podem ser obtidas com o comando model.tables
 subd.m <- model.tables(subd.av, ty="means")
 subd.m
 ### Análise de resíduos
 par(mfrow=c(2,2))
 plot(subd.avb)
 ### Normalidade
 shapiro.test(subd.avb$res)
 ### Homocedacicidade
 par(mfrow=c(2,2))
 boxplot(subd.avb$res~trat)
 boxplot(subd.avb$res~aveia)
 boxplot(subd.avb$res~aveia:trat)
 # O bartlet pode ser executado, mas sempre lembramos de suas restrições pois ele 
 # não e capaz de analisar as variações oriundas da interação
 bartlett.test(subd.avb$res, trat)
 bartlett.test(subd.avb$res, aveia)
 ##----------------------------------------------------------------------------##
 #Testes de comparação multipla:
 require(agricolae)
 # O Correto seria apenas realizar a analise da interação, mas vamos realizar também
 # para todas os fatores e interação como forma de ilustrar os comandos, lembrando 
 # que para analisar desta formas devemos considerar que a interação não foi 
 # significativa (Lembrando que a nossa foi).
 
 # Fator principal aveia
 subd.anovaRA<-anova(subd.av$'bloco:a')
 subd.anovaRA
 # extraindo os QM e DF de "a"
 QMEa<-subd.anovaRA$Mean[2]
 QMEa
 Dfea<-subd.anovaRA$Df[2]
 Dfea
 HSD.test(resp,aveia,Dfea,QMEa,alpha=0.05)
 ##----------------------------------------------------------------------------##
 # Fator secundario trat
 subd.anovaRB<-anova(subd.av$'Within')
 subd.anovaRB
 # extraindo os QM e DF de "a"
 QMEb<-subd.anovaRB$Mean[2]
 QMEb
 Dfeb<-subd.anovaRB$Df[2]
 Dfeb
 HSD.test(resp,trat,Dfeb,QMEb,alpha=0.05)
 ##----------------------------------------------------------------------------##
 # Desdobrando e testanto a interação "aveia:trat"
 # Primeiro vamos estudar o fatores secundarios dentro dos principais
 interaction.plot(aveia,trat, resp, xlab="Variedades")
 subd.avISP<-aov(resp ~ bloco +aveia/trat + aveia*trat -trat  + Error(bloco:aveia))
 summary(subd.avISP)
 effects(subd.avISP$'Within')
 summary(subd.avISP,split=list("aveia:trat"=list(b.dentro.a1=c(1,5,9),
                                                 b.dentro.a2=c(2,6,10),
                                                 b.dentro.a3=c(3,7,11),
                                                 b.dentro.a4=c(4,8,12))))   
 ## Aplicação do teste Tukey
 require(agricolae)
 QMRb<-20.311
 DFRb<-36
 tapply(subd$resp,list(Aveia=subd$aveia,Tratamento=subd$trat),mean)
 #               Tratamento
 #Aveia      1      2      3     4
 #    1 36.050 50.625 45.850 37.30
 #    2 50.850 55.375 53.100 54.30
 #    3 53.925 51.375 55.875 56.05
 #    4 61.925 63.425 57.675 61.25
 
 ## b dentro de A1
 A1<-subset(subd,subd$aveia=="1")
 A1
 HSD.test(A1$resp,A1$trat,DFRb,QMRb)
 ## b dentro de A2
 A2<-subset(subd,subd$aveia=="2")
 HSD.test(A2$resp,A2$trat,DFRb,QMRb)
 ## b dentro de A3
 A3<-subset(subd,subd$aveia=="3")
 HSD.test(A3$resp,A3$trat,DFRb,QMRb)
 ## b dentro de A4
 A4<-subset(subd,subd$aveia=="4")
 HSD.test(A4$resp,A4$trat,DFRb,QMRb)
 ##----------------------------------------------------------------------------##
 ## Agora vamos estudar a interação dos efeitos principais dentro dos secundarios
 interaction.plot(trat,aveia, resp, xlab="Tratamentos de Sementes")
 #               Tratamento
 #Aveia      1      2      3     4
 #    1 36.050 50.625 45.850 37.30
 #    2 50.850 55.375 53.100 54.30
 #    3 53.925 51.375 55.875 56.05
 #    4 61.925 63.425 57.675 61.25
 
 subd.avIPS<-aov(resp ~ bloco + trat/aveia + trat*aveia -aveia)
 summary(subd.avIPS)
 effects(subd.avIPS)
 summary(subd.avIPS,split=list("trat:aveia"=list(a.dentro.b1=c(1,5,9),
                                                 a.dentro.b2=c(2,6,10),
                                                 a.dentro.b3=c(3,7,11),
                                                 a.dentro.b4=c(4,8,12))))
 # Ponderações devem ser feitas, sobre está ANOVA, a soma de quadrados do desdodramentos
 # estão corretas, porem o teste F e p-valor, não!!! visto que o error não está 
 # calculado adequadamente, neste tipo de desdobramento deve calcular um erro medio,
 #  dado por: QMex = (QMa+(b-1)*QMb)/b , sendo b número de tratamentos secundarios.
 #  Os graus de liberdade do erros também devem ser corrigidos pela formula de SATTERTHWAITE:
 #  [QMa+(b-1)*QMb]^2
 #--------------------
 #(QMa)^2+[(b-1)*QMb]^2
 #------  ------------
 # gla        glb
 
 summary(subd.av)
 ##Erro A
 QMa<-68.70
 Gla<-9
 ##Erro B
 QMb<-20.311
 Glb<-36
 #N° de tratamentos secundarios
 NTb<-4
 ##Erro médio
 QMex<-(QMa+(NTb-1)*QMb)/NTb
 QMex
 ##Graus de liberdade
 Glx<-((QMa+(NTb-1)*QMb)^2) / ( ((QMa^2)/Gla)  + ((((NTb-1)*QMb)^2)/Glb) )
 Glx
 ##Com o novo erro e Gl podemos agora aplicar o teste F adequadamente e realizar 
 # também as comparações multiplas
 interaction.plot(trat,aveia, resp, xlab="Tratamentos de Sementes")
 tapply(subd$resp,list(Aveia=subd$aveia,Tratamento=subd$trat),mean)
 #               Tratamento
 #Aveia      1      2      3     4
 #    1 36.050 50.625 45.850 37.30
 #    2 50.850 55.375 53.100 54.30
 #    3 53.925 51.375 55.875 56.05
 #    4 61.925 63.425 57.675 61.25
 
 ## aveia dentro de trat1
 B1<-subset(subd,subd$trat=="1")
 B1
 HSD.test(B1$resp,B1$aveia,Glx,QMex)
 ## aveia dentro de trat2
 B2<-subset(subd,subd$trat=="2")
 HSD.test(B2$resp,B2$aveia,Glx,QMex)
 ## aveia dentro de trat3
 B3<-subset(subd,subd$trat=="3")
 HSD.test(B3$resp,B3$aveia,Glx,QMex)
 ## aveia dentro de trat4
 B4<-subset(subd,subd$trat=="4")
 HSD.test(B4$resp,B4$aveia,Glx,QMex)
\end{Sinput}
\end{Schunk}
\newpage
\section{Regressão Linear}
\begin{Schunk}
\begin{Sinput}
 ###-------------------------------------------------------------------------####
 ############################# Regressão linear #################################
 # Os dados foram retirados do livro "Curso de Estatística Experimental 15ed" do
 # Frederico Pimentel-Gomes pagina 238. Os dados são referentes ao teor de K trocável
 # no solo (X) em miliequivalentes, a variavel mensurada foi produção de Cana em t/ha (Y),
 # a avaliação foi realizada depois da aplicação de 150 kg/ha de K2O.
 
 x<-c(0.109,0.138,0.206,0.153,0.156,0.240,0.270,0.126,0.199,0.112,0.138,0.103,0.127,0.063)
 y<-c(31.8,24.5,11.8,18.8,17.3,11.0,12.2,20.6,5.3,29.3,8.0,35.8,19.6,21.4)
 ### Explorando visualmente os dados
 plot(x,y,xlab="Teor de K",ylab="Aumento de produção t/ha",main="Cana versus K")
 # Correlação linear de pearson entre as variaveis
 cor.test(x,y) 
 # Começamos relizando o ajuste de um reta y=a+b*x, sendo este o ajuste o modelo
 # mais simples, com a função "lm" podemos realizar os ajustes lineares
 ajuste1<-lm(y~x)
 summary(ajuste1)
 names(ajuste1)
 # Dado o ajuste realizado, executaremos a ANOVA pertinente a regressão, onde sob 
 # H0 todos os coeficientes são iguais a 0, e sob H1 pelo menos um deles é diferente de
  
 anova(ajuste1)
 # Com o modelo ajustado devemos começar a realizar as analises dos residuos para 
 # verificar se os pressupostos foram atendidos. Neste primeiro monento faremos isso de 
 # maneira grafica.
 
 par(mfrow=c(2,2))
 plot(ajuste1)
 # Adicionalmente a analise grafica podemos realizar o teste de Shapiro-Wilk, onde sob
 # H0 os residuos tem distribuição normal, e sob H1 não tem distribuição normal
 
 shapiro.test(ajuste1$resi)
 # Verificado o ajuste e os pressupostos podemos plotar os dados e a equação estimada.
 
 plot(x,y,xlab="Teor de K",ylab="Aumento de produção t/ha",main="Cana versus K")
 abline(ajuste1)
 # Como segunda opção podemos ajustar y=a+b*x+c*x^2, utilizamos a função "I()" para
 # especificar modelos lineares polinomiais como neste caso
 
 ajuste2<-lm(y~I(x)+I(x^2))
 summary(ajuste2)
 anova(ajuste2)
 par(mfrow=c(2,2))
 plot(ajuste2)
 shapiro.test(ajuste2$resi)
 # Gerando o grafico
 xx<-seq(min(x),max(x),by=0.01)
 yy<-predict(ajuste2,list(x=xx))
 plot(x,y,xlab="Teor de K",ylab="Aumento de produção t/ha",main="Cana versus K")
 lines(xx,yy)
 abline(ajuste1)
 ### limpando a memoria
 rm(list=ls())
\end{Sinput}
\end{Schunk}
\section{Regressão Linear Multipla}
\begin{Schunk}
\begin{Sinput}
 ###--------------------------------------------------------------------------###
 ######################### Regressão linear Multipla ############################
 # Os dados são parte de um experimento conduzido por Fernando Bortolozzo, onde foi
 # estudado a perda de nutrientes via escoamento superficial, neste para nosso objetivo
 # utilizaremos os dados de concentração de P total expressos em mg/l.
 
 dados <- read.table('Escoamento.txt',head=TRUE)
 summary(dados)
 # Explorando Visualmente os dados
 with(dados,plot(faixa,C_Ptotal))
 with(dados,plot(tempo,C_Ptotal))
 # Verificando a estrutura de correlação entre as variaveis
 require(agricolae)
 correlation(dados)
 # Na sequencia ajustaremos um modelo para cada variavel de forma independente
 # ajuste faixa
 ajuste.f<-lm(C_Ptotal~faixa,data=dados)
 summary(ajuste.f)
 anova(ajuste.f)
 par(mfrow=c(2,2))
 plot(ajuste.f)
 shapiro.test(ajuste.f$res)
 ###--------------------------------------------------------------------------###
 # Na sequencia ajustaremos um modelo para cada variavel de forma independente
 # ajuste tempo
 ajuste.t<-lm(C_Ptotal~tempo,data=dados)
 summary(ajuste.t)
 anova(ajuste.t)
 par(mfrow=c(2,2))
 plot(ajuste.t)
 shapiro.test(ajuste.t$res)
 ###--------------------------------------------------------------------------###
 # Regressão linear Multipla
 # Para ajuste multiplos simplismente adicionamos as variaveis "X1 + X2 ... Xn-1 + Xn"
 # no modelos "lm(Y~X1+X2...+Xn-1+Xn)"
 
 ajuste.m<-lm(C_Ptotal~tempo+faixa,data=dados)
 summary(ajuste.m)
 anova(ajuste.m) 
 # Analisando os residuos do ajuste
 par(mfrow=c(2,2))
 plot(ajuste.m)
 shapiro.test(ajuste.m$res)
 # Retirando o intercepto
 ajuste.m<-lm(C_Ptotal~tempo+faixa-1,data=dados)
 summary(ajuste.m)
 anova(ajuste.m)
 # Retirando o intercepto e adicionando o termo quadratico  e cubicos
 ajuste.mI<-lm(C_Ptotal~I(tempo)+I(faixa)+I(tempo^2)+I(faixa^2)+I(tempo^3)+I(faixa^3)-1,data=dados)
 summary(ajuste.mI)
 anova(ajuste.mI)
 # Analisando os residuos do ajuste
 par(mfrow=c(2,2))
 plot(ajuste.mI)
 shapiro.test(ajuste.mI$res)
 # plotando em 3d, para isso utilizaremos o pacote "scatterplot3d", caso não tenha
 # ele instalado: install.packages("scatterplot3d",dep=T)
 require(scatterplot3d)
 scatterplot3d(dados$faixa,dados$tempo,dados$C_Ptotal)
 ###--------------------------------------------------------------------------###
 #mGerando series para plotar a função ajustada
 dadosG <- expand.grid(tempoS = seq(10,120,length=50),faixaS= seq(5,30,length=50))
 plot(dadosG)
 # realizando a predição do modelo
 dadosG$Ptotal <-predict(ajuste.mI,list(faixa=dadosG$faixaS,tempo=dadosG$tempoS))
 # Plotando o modelo em 3d
 require(lattice)
 wireframe(Ptotal~tempoS*faixaS,
 data=dadosG,
 xlab=list("Tempo (min)", rot=32),
 ylab=list("Largura faixa vegetada (m)", rot=-40),
 zlab=list("Ptotal", rot=95),
 cex.lab=0.15,
 colorkey = TRUE,
 drape = TRUE,
 col.regions=rgb(colorRamp(c("green","yellow","red"))(seq(0, 1, length =100)), max = 255),
 col=0,
 scales=list(arrows = FALSE,x=list(at=seq(10,120,by=20)),y=list(at=seq(5,30,by=5))))
 # Plotando as curvas de nivel do modelo    
 contourplot(Ptotal~tempoS*faixaS,
 data=dadosG,
 col.regions=colorRampPalette(c("green","yellow","red")),
 cuts=15,
 region =TRUE,
 xlab="Tempo (min)",
 ylab="Largura da faixa vegetada (m)")
 ###--------------------------------------------------------------------------###
 # Normalmente quando ajustamos superficies de respostas estamos interessados em 
 # Encontrar pontos criticos da função como maximos e minimos, para isso de forma 
 # numerica utilizamos as funções abaixo
 # Função para máximo
 fmax<-function(dados){
     out<-dados[which.max(dados[,3]),]
     return(out)
     }
 # Função para minimo
 fmin<-function(dados){
     out<-dados[which.min(dados[,3]),]
     return(out)
     }
 # Encontrando os pontos Máximo e Minimo
 Zmax<-fmax(dadosG)
 Zmax
 Zmin<-fmin(dadosG)
 Zmin
 ###--------------------------------------------------------------------------###
 ########################### Seleção de modelos #################################
 # Na pesquisa a mensuração de inumeras variaveis acontece com frequencia, então 
 # encontrar de maneira direta quais as variaveis que estão correlacionadas com a 
 # variavel resposta torna-se trabalhosa a operação, baseado nesta necessidade
 # podemos utilizar procedimentos baseados no criteiro de AIC (Akaike Information Criterion)
 # para definir as variaveis de interesse.
 library(MASS)
 ##?stepAIC
 # Utilizando o metodo forward
 modelo.f <- lm(C_Ptotal~1,data=dados)# Modelo somente com intercpto
 result.f <- stepAIC(modelo.f,
 direction="forward",# partimos do intercepto e vamos adicionando uma variavel por vez
 scope=list(
 upper = ~I(tempo)+I(faixa)+I(tempo^2)+I(faixa^2),
 lower = ~1),trace=1)
 summary(result.f)
 # Utilizando o metodo both
 modelo.b <- lm(C_Ptotal~1,data=dados)# Modelo somente com intercpto
 result.b<- stepAIC(modelo.b,
 direction="both",# partimos do intercepto e vamos adicionando uma variavel por vez
 scope=list(
 upper = ~I(tempo)+I(faixa)+I(tempo^2)+I(faixa^2),
 lower = ~1),trace=1)
 summary(result.b)
 # Utilizando o metodo both
 modelo.ba <- lm(C_Ptotal~I(tempo)+I(faixa)+I(tempo^2)+I(faixa^2),data=dados)
 result.ba <- stepAIC(modelo.ba,
 direction="backward",# partimos do intercepto e vamos adicionando uma variavel por vez
 scope=list(
 upper = ~I(tempo)+I(faixa)+I(tempo^2)+I(faixa^2),
 lower = ~1),trace=1)
 summary(result.ba)
 rm(list=ls())
\end{Sinput}
\end{Schunk}
\newpage
\section{Regressão Não Linear}
\begin{Schunk}
\begin{Sinput}
 ####------------------------------------------------------------------------####
 ######################### Regressão não linear #################################
 # Simulando dados
 x<-seq(0,5,l=50)
 set.seed(123)
 y <- 2.5*0.6^x+rnorm(50,sd=0.1)
 plot(x,y)
 data <- data.frame(x=x,y=y)
 attach(data)
 ###--------------------------------------------------------------------------###
 # Algoritimo de Gauss-Newton
 model <- nls(y~a*b^x,data=data,start=list(a=0.5,b=2),trace=TRUE)
 summary(model)
 ###--------------------------------------------------------------------------###
 # Algoritimo de Golub-Pereyra
 #model.GP <- nls(y~a*b^x,data=data,start=list(a=0.5,b=2),alg='plinear',trace=TRUE)
 #summary(model.GP)
 ###--------------------------------------------------------------------------###
 # Algoritimo de port
 model.P <- nls(y~a*b^x,data=data,start=list(a=0.5,b=2),alg='port',trace=TRUE)
 summary(model.P)
 ###--------------------------------------------------------------------------###
 # Outros argumentos podem ser modificados na "nls"
 # ?nls
 # ?nls.control
 control=nls.control(maxiter = 50,### Numero de interações
 tol = 1e-05,### Tolerancia da precisão
 minFactor = 1/1024, ### Tamanho minimo do passo
 printEval = FALSE,### mostrar os passos da interação
 warnOnly = FALSE)### mostrar erro antes do termino do processo em caso de matrix singular
 ###--------------------------------------------------------------------------###
 # Plotando o modelo
 # ?predict
 xp=seq(0,5,l=300)
 plot(x,y)
 lines(xp,predict(model,list(x=xp)))
 # Analisando o residuo
 res <- y - predict(model,list(x=x))
 hist(res)
 shapiro.test(res)
 # R2
 SQE <- summary(model)$sigma^2*summary(model)$df[2]
 SQT <- var(y)*(length(y)-1)
 R2 <- 1-SQE/SQT
 R2
 # Intervalo de confiança para parametros
 confint(model,level=0.95)
 ###--------------------------------------------------------------------------###
 ### Entendendo o processo interativo
 require(nls2)
 require(lattice)
 Pa <- seq(0.1,5,l=25)
 Pb <- seq(0.0001,1,l=25)
 Erro <- expand.grid(a=Pa,b=Pb)
 Erro <- cbind(Erro,erro=numeric(length(Erro[,1])))
 head(Erro)
 plot(Erro$a,Erro$b,pch='+',xlab="Parametro a",ylab="Parametro b",
 main='Malha para os parametros')
 ### Algoritimo de Força Bruta
 model.BT <- nls2(y~a*b^x,data=data,start=list(Erro[,1:2]),alg='brute-force')
 summary(model.BT)           
 for (i in 1:length(Erro[,1]))
   {
     Erro[i,3] <- summary(nls2(y~a*b^x,data=data,start=list(a=Erro[i,1],b=Erro[i,2]),
                          alg='brute-force',trace=FALSE))$sigma
   }
 levelplot(Erro$erro~Erro$a*Erro$b,cuts=100,xlab="Parametro a",ylab="Parametro b",
 main='Função Erro',
           col.regions= heat.colors(200))
 contourplot(Erro$erro~Erro$a*Erro$b,cuts=50,xlab="Parametro a",ylab="Parametro b",
 main='Função Erro')
 x11()
 wireframe(Erro$erro~Erro$a*Erro$b,
 shade = TRUE,xlab="Parametro a",
 ylab="Parametro b",zlab="Erro",
 main='Função Erro',
 aspect = c(61/87, 0.4),
 light.source = c(10,0,100))
 wireframe(Erro$erro~Erro$a*Erro$b,
 xlab="Parametro a",
 ylab="Parametro b",
 main='Função Erro',
 scales = list(arrows = TRUE),
 drape = TRUE, colorkey = TRUE,col.regions= heat.colors(100))
 ###--------------------------------------------------------------------------###
 ### Utilizando as derivadas parciais
 ### Obtendo as derivas parciais com "?D" ou "?deriv"
 ## Derivada em função de a
 D(expression(a*b^x),'a')
 ## Ou
 deriv(expression(a*b^x),'a')
 ## Derivada em função de b
 D(expression(a*b^x),'b')
 ## Ou
 deriv(expression(a*b^x),'b')
 ### Podemos derivar qualquer outra expressão
 D(expression(a*log(b^x^sqrt(a/tan(b)))),'b')
 #### Motando a função de atribuido das derivadas parciais a função
 der.model <- function(a,b,x)
   {
     model.func <- a*b^x ## Função do modelo
     J <- cbind(b^x,a*x*b^(x-1)) ### Jacobiana
     dimnames(J) <- list(NULL,c('a','b'))
     attr(model.func,'gradient') <- J #Atribuir ao modelo J
     model.func
   }
 ## Ajustando o modelo utilizando as derivadas
 model.Dp <- nls(y~der.model(a,b,x),data=data,start=list(a=0.5,b=2),trace=TRUE)
 summary(model.Dp)
\end{Sinput}
\end{Schunk}

\end{document}

