################################################### ### chunk number 1: ################################################### options(prompt=" ", continue=" ") ################################################### ### chunk number 2: ################################################### ########## 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 ################################################### ### chunk number 3: ################################################### ###--------------------------------------------------------------------------### 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]) ###--------------------------------------------------------------------------### ################################################### ### chunk number 4: ################################################### ################## 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))) ################################################### ### chunk number 5: ################################################### ###-------------------------------------------------------------------------#### ################################## 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 ################################################### ### chunk number 6: ################################################### ########################## 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) ################################################### ### chunk number 7: ################################################### ################################## 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) ################################################### ### chunk number 8: ################################################### ###-------------------------------------------------------------------------#### ############################# 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()) ################################################### ### chunk number 9: ################################################### ###--------------------------------------------------------------------------### ######################### 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()) ################################################### ### chunk number 10: ################################################### ####------------------------------------------------------------------------#### ######################### 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) ###--------------------------------------------------------------------------### ################################################### ### chunk number 11: ################################################### options(prompt=">", continue="+")