###---------------------------------------------------------------------------------------------------------------------#### ########################################## introdução ao R ################################################################# ### Para está seção usaremos material do Prof Adilson, Paulo e Daniel http://www.leg.ufpr.br/~paulojus/embrapa/Rembrapa/ http://www.dex.ufla.br/~danielff/RRC0.pdf ###---------------------------------------------------------------------------------------------------------------------#### ##################################################### DIC ################################################################## #Delineamento Completamente Casualizado ## Os dados foram retirados da apostila do professor Adilson dos Anjos, são referentes a 9 linhagens de fungos onde foi avaliado a taxa de crescimento em micras por hora. #Este é um exemplo do formato em que se deve estar os dados que será inserido no bloco de notas, para depois o software realizar a chamada (leitura da base de dados). #Relizando a limpeza de todas os dados armazenados na memoria do programa rm(list=ls()) #Carregando o arquivo que contem o conjunto de dados, e bom dar uma olhada em outras funções similares com "?read.csv" e "?read.delim" ex01<-read.table("exemplo01.txt", head=T) ex01 #inspecionando o classe de dados is.data.frame(ex01) #Verificando as variaveis contida no conjunto de dados names(ex01) #Verificando o conteudo das variaveis ex01$resp ex01$trat #Verificando a classe de cada variaveis is.factor(ex01$trat) is.numeric(ex01$resp) #Lembramos que resp deve ter conteudo numerico e trat deve ser fator ex01$trat<-as.factor(ex01$trat) ex01$resp<-as.numeric(ex01$resp) #Com o comando "?summary" obtemos uma analise descritivva dos dados summary(ex01) #Com a função "?by", podemos examinar os dados de uma forma mais eficiente sumarygeral<-by(ex01$resp,ex01$trat, summary) sumarygeral #Podemos também extrair algumas de forma isolada by(ex01$resp,ex01$trat, mean) by(ex01$resp,ex01$trat, var) by(ex01$resp,ex01$trat, sd) by(ex01$resp,ex01$trat, sum) by(ex01$resp,ex01$trat, median) by(ex01$resp,ex01$trat, range) by(ex01$resp,ex01$trat, min) by(ex01$resp,ex01$trat, max) #Para facilitar a manipulação dos dados podemos attachar o objeto "ex01" que contem nossa variaveis de interesse, após este comando as variaveis ficaram na memoria do R attach(ex01) #Para realizar o inverso do attach utilizamos a função detach que tira o objeto da memoria #detach(ex01) #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" ex01.m <- tapply(resp, trat, mean) ex01.m ex01.v <- tapply(resp, trat, var) ex01.v ex01.sd <- tapply(resp, trat, sd) ex01.sd #Como inspeção grafica dos dados podemos utilizar o comando "?plot" ou "?boxplot" plot(ex01,main="Taxa de crescimento de fungos",xlab="Fungos",ylab="micras/hora") #Com a função points podemos colocar a media sobre o boxplot, outra função similar e interessante e a função "?lines" points(ex01.m, pch="x", col=2, cex=1.5) # Da mesma maneira como foi colocado a medias podemos acresentar o desvio padrão points(ex01.m+ex01.sd, pch="-", col=3, cex=2.5) points(ex01.m-ex01.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 ex01.av <- aov(resp ~ trat, data = ex01) ex01.av #A obtenção do quadro tipico da ANOVA a qual conheçemos pode ser feito com o comando "?summary" ou "?anova" summary(ex01.av) anova(ex01.av) # Alem do quadra da anova o objeto ex01.av guarda outras informações interessantes names(ex01.av) ex01.av$coef ex01.av$res residuals(ex01.av) #Análise de resíduos #Homocedasticidade: Graficamente, pode-se analisar a homocedasticidade através de um box-plot, boxplot(ex01.av$res ~ trat,ylab="Resíduos",xlab="Linhagens") #Através de um gráfico dos resíduos vs tratamentos, plot.default(trat,ex01.av$res,ylab="Resíduos",xlab="Linhagens") #Ainda, pode-se avaliar através de um teste, como por exemplo o teste de Bartlett: bartlett.test(ex01.av$res, trat) #Normalidade: Graficamente, pode-se avaliar a normalidade dos resíduos fazendo hist(ex01.av$res, main="Histograma dos Resíduos") #ou stem(ex01.av$res) #ou qqnorm(ex01.av$res,ylab="Resíduos", main=NULL) qqline(ex01.av$res) title("Grafico Normal de Probabilidade dos Resíduos") #Teste de Shapiro-Wilk para Normalidade shapiro.test(ex01.av$res) #Independência: A independência, com algumas restrições, pode ser analisada graficamente, através de plot(ex01.av$fit, ex01.av$res, xlab="valores ajustados", ylab="resíduos") title("resíduos vs Preditos") #Ainda é possível avaliar algum tipo de dependência através da ordenação dos resíduos, caso exista uma ordem de obteção dos dados conhecida: plot(ex01.av$fit, order(ex01.av$res), xlab="valores ajustados", ylab="resíduos") title("resíduos vs Preditos") # Obtenção dos residuos padronizados, com ele podemos verificar a presença de algum outline seguindo os criterios de "x" desvios padrões names(anova(ex01.av)) s2 <- anova(ex01.av)$Mean[2] res<-ex01.av$res respad<-(res/sqrt(s2)) boxplot(respad,main="Resíduos Padronizados" ) #Verificar o residuo padronizado de cada tratamento plot.default(ex01$trat,respad, xlab="Linhagens",main="Resíduos Padronizados" ) #A localização de um ponto dentro do grafico pode ser realizada com a função "?identify" identify(ex01$trat,respad) #A analise dos residuos graficamente pode ser obtida executando um comando plot no objeto que recebeu o modelo, como ele gera automaticamente 4 graficos devemos inicialmente particionar a janela grafica par(mfrow=c(2,2)) plot(ex01.av) #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: ex01.tu <- TukeyHSD(ex01.av) ex01.tu plot(ex01.tu) #Utilização de um pacote interessante: require(agricolae) # O pacote agricolar tem a função de comparação atraves do Tukey implementada com uma saida mais interessante "?HSD.test", onde nesta função fornecemos os tratamento, os graus de liberdade e soma de quadrados do erro e significancia. ex01.anova<-anova(ex01.av) ex01.anova tuk<-HSD.test(resp,trat,ex01.anova$Df[2],ex01.anova$Mean[2],alpha=0.05) ### Grafico também são implementados no pacote agricolae, é aconselhavel dar uma boa olhada na documentação deste pacote, vale a pena!!!!!!!. bar.group(tuk,border="red",col="blue",ylim=c(0,550),main="Meu grafico",xlab="tratamento",y="media") bar.err(tuk,border="red",col="blue",ylim=c(0,550),main="Meu grafico",xlab="tratamento",y="media",add=T) ### Comparação via contrastes #OBS: Explicar teoria e tecnica dos Balões cont.ex01<-matrix(c(.25,.25,.25,.25,-.2,-.2,-.2,-.2,-.2,#1° Contraste 1,1,-1,-1,0,0,0,0,0,#2° Contraste 1,-1,0,0,0,0,0,0,0,#3° Contraste 0,0,1,-1,0,0,0,0,0,#4° Contraste 0,0,0,0,1,1,1,1,-4,#5° Contraste 0,0,0,0,1,1,-1,-1,0,#6° Contraste 0,0,0,0,1,-1,0,0,0,#7° Contraste 0,0,0,0,0,0,-1,1,0),nrow=9,ncol=8,byrow=F) ### Definindo os contrastes contrasts(ex01$trat)<-cont.ex01 contrasts(ex01$trat) ex01$trat ### Executando a ANOVA atraves de contrastes ex01.avC <- aov(resp ~ trat, data = ex01) summary(ex01.avC) summary(ex01.av) ### Contrastes estabelecidos ex01.avC$con ex01.av$con summary(ex01.avC,split=list(trat=list(C1=1,C2=2,C3=3,C4=4,C5=5,C6=6,C7=7,C8=8))) summary(ex01.avC,split=list(trat=list(C1=1,C2=2,C3=3,C4=4))) ###---------------------------------------------------------------------------------------------------------------------#### # EXERCICIO, Com o conjunto de dados "exemplo01E.txt" que contem dados retirados do livro "Estatistica Aplicada a pesquisa agícola" de Zimmermann, o exemplo enontrase na pagina 62 e foi modificado, os tratamentos são? #1-Solo sem cultivo previo de arroz, com furadam #2-Solo sem cultivo previo de arroz, sem furadam #3-Solo sem cultivo previo de arroz, com herbicida #4-Solo sem cultivo previo de arroz, sem herbicida #5-Solo com cultivo previo de arroz, com herbicida #Os dados são de materia seca por planta, em gramas, o experimento foi conduzido em casa de vegetação em vasos. #Dado O exemplo, obtenha: a analise exploratoria dos dados, analise de variancia, e teste de comparação multipla caso for necessario, e obtenha a comparação por contraste: #1-Com versus sem cultivo previo #2-Furadam versus herbicida #3-Sem furadam versus sem herbida #DEIXE TODOS OS COMANDOS NA SEQUENCIA ABAIXO PARA CORREÇÃO, OBS: Comente os comandos, ajudara na compreensão ###---------------------------------------------------------------------------------------------------------------------#### ##################################################### DBA ################################################################## ### Delinemento blocos ao acaso ### dados de produtividade de feijão em kg/ha retirado do livro Estatistica aplicada a pesquisa agrícola (Zimmmerman) pg 72. rm(list=ls()) ex02<-read.table("exemplo02.txt", head=T) ex02 #explorando os dados summary(ex02) ex02$bloco<-as.factor(ex02$bloco) ex02$trat<-as.factor(ex02$trat) attach(ex02) ex02.m <-tapply(prod, trat, mean) ex02.m ex02.v <-tapply(prod, trat, var) ex02.v boxplot(prod~trat) points(ex02.m, pch="x", col=2, cex=1.5) #Executando a analise de variância ex02.av<-aov(prod~bloco+trat) summary(ex02.av) #Obtendo o CV ex02.anova<-anova(ex02.av) CV<-((ex02.anova$Mean[3]^0.5)/(model.tables(ex02.av, "means")$tables$Gran))*100,dig=2 CV #Verificando os residuos par(mfrow=c(2,2)) plot(ex02.av) ###Homocedaciade bartlett.test(ex02.av$res~trat) ###normalidade shapiro.test(ex02.av$res) ### executando um Tukey (não realizado normalmente, lembrado que tratamento é não significativo) ex02.tu <- TukeyHSD(ex02.av) ex02.tu plot(ex02.tu) #Utilização de um pacote um interessante: require(agricolae) ex02.anova<-anova(ex02.av) ex02.anova names(ex02.anova) tuk<-HSD.test(prod,trat,ex02.anova$Df[3],ex02.anova$Mean[3],alpha=0.05) bar.group(tuk,border="red",col="blue",ylim=c(0,2550),main="Meu grafico",xlab="tratamento",y="media",cex.names=0.7) bar.err(tuk,border="red",col="blue",ylim=c(0,2550),main="Meu grafico",xlab="tratamento",y="media",cex.names=0.7,add=T) ###---------------------------------------------------------------------------------------------------------------------#### #EXERCICIO: Com o conjunto de dados "exemplo02E.txt" que contem dados retirados da tese de Adriana (Solos-UFPR), de produtividade de trigo sobre diferentes doses de esterco liquido bovino, foi conduzido com 4 blocos e 4 doses (0,60,120,180). #Obtenha:Analise exploratoria, Analise de variancia, teste de comparação multipla, e recomendações. #DEIXE TODOS OS COMANDOS NA SEQUENCIA ABAIXO PARA CORREÇÃO, OBS: Comente os comandos, ajudara na compreensão ###---------------------------------------------------------------------------------------------------------------------#### ############################################## Transformação de dados ###################################################### #Objetivo é apresentar um método de transformação de dados, para um modelo que não atendeu os pressupostos exigidos. #O exemplo é número de falhas eletrônicas em diferentes sistemas de atendimento telefônico no período de um ano de funcionamento. #O delineamento é completamente casualizado. rm(list=ls()) ex03<-read.table("exemplo03.txt", head=T) ex03 #explorando os dados summary(ex03) ex03$trat<-as.factor(ex03$trat) summary(ex03) plot(ex03) hist(ex03$resp) attach(ex03) ##Análise de variância ex03.av <- aov(resp ~ trat) ex03.av anova(ex03.av) ##Homocedasticidade plot.default(ex03$trat,ex03.av$res) bartlett.test(ex03$resp, ex03$trat) ##Normalidade qqnorm(ex03.av$res,ylab="Residuos", main=NULL) qqline(ex03.av$res) title("Grafico Normal de Probabilidade dos Resíduos") shapiro.test(ex03.av$res) #Com os resultados acima resta realizar a transformação Box-Cox, pela seguinte expressão. #A função boxcox calcula a verossimilhança perfilhada do parâmetro. Deve-se escolher o valor que maximiza esta função. #Carregando o pacote MASS, aonde pode ser obtida a função boxcox. require(MASS) box.ex03<- boxcox(resp ~ trat, data=ex03, plotit=T) box.ex03<- boxcox(resp ~ trat, data=ex03, lam=seq(-1, 1, 1/10)) #O gráfico mostra que o valor que maximiza a função é aproximadamente . #Para obter o valor exato de lambda, e não ficar somente com a análise gráfica. lambda <- box.ex03$x[which(box.ex03$y == max(box.ex03$y))] lambda Transformando os dados com lambda encontrado. ex03$respt<- (ex03$resp^(lambda)-1)/lambda ex03 attach(ex03) ex03.avt <- aov(respt ~ trat) anova(ex03.avt) par(mfrow=c(2,2)) plot(ex03.avt) ##Homocedasticidade plot.default(ex03$trat,ex03.avt$res) bartlett.test(ex03$respt, ex03$trat) ##Normalidade qqnorm(ex03.avt$res,ylab="Residuos", main=NULL) qqline(ex03.avt$res) title("Grafico Normal de Probabilidade dos Resíduos") shapiro.test(ex03.avt$res) #Caso o gráfico da verossimilhança perfilhada o intervalo de confiança para o valor 0 está contido neste intervalo. #Pode ser realizado a transformação logarítmica, abaixo os comandos. ex03.avl = aov(log(resp) ~ trat) par(mfrow=c(2,2)) plot(ex03.avl) ## Comparações multiplas ex03.anova<-anova(ex03.av) ex03.anova ex03.anovaT<-anova(ex03.avt) ex03.anovaT ex03.tu<-HSD.test(resp,trat,ex03.anova$Df[2],ex03.anova$Mean[2],alpha=0.05) ex03.tuT<-HSD.test(respt,trat,ex03.anovaT$Df[2],ex03.anovaT$Mean[2],alpha=0.05) ex03.tu ex03.tuT ###---------------------------------------------------------------------------------------------------------------------#### ##################################################### DQL ################################################################## #A análise de variância de um experimento conduzido no Delineamento em Quadrado Latino. #O experimento é sobre produção de grãos de feijão (Kg/parcela) de diferentes variedades de feijão. rm(list=ls()) ex04<- read.table("exemplo04.txt", head=T) ex04 #Explorando os dados. names(ex04) summary(ex04) ex04$coluna<-as.factor(ex04$coluna) ex04$linha<-as.factor(ex04$linha) ex04$trat<-as.factor(ex04$trat) attach(ex04) par(mfrow=c(2,2)) plot(resp ~ coluna + linha + trat) ex04.mt = tapply(resp, trat, mean) ex04.mt ex04.ml = tapply(resp, linha, mean) ex04.ml ex04.mc = tapply(resp, coluna, mean) ex04.mc plot.default(trat, resp) boxplot(resp~trat) points(ex04.mt, pch="x", col=2, cex=1.5) ### Não devem ocorrer interações entre linhas e colunas no experimento. interaction.plot(linha,trat, resp) interaction.plot(coluna,trat, resp) ### ANOVA ex04.av <- aov(resp ~ coluna+linha+trat) anova(ex04.av) ### Analisando os pressuposto par(mfrow=c(2,2)) plot(ex04.av) ### Homocedasticidade, bartlett.test(ex04.av$res,trat) ### Normalidade dos Resíduos shapiro.test(ex04.av$res) ### Comparações multiplas ex04.tu = TukeyHSD(ex04.av, "trat", ord=T) ex04.tu plot(ex04.tu) require(agricolae) ex04.anova<-anova(ex04.av) ex04.anova HSD.test(resp,trat,ex04.anova$Df[4],ex04.anova$Mean[4],alpha=0.05) ###---------------------------------------------------------------------------------------------------------------------#### ############################################### 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" Kronha & Banzatto pg 106. rm(list=ls()) ex05 <- read.table("exemplo05.txt", header=T) ex05 names(ex05) summary(ex05) attach(ex05) is.factor(rec) is.factor(esp) is.factor(resp) is.numeric(resp) ### Explorando os dados ex05.m<-tapply(resp, list(rec,esp), mean) ex05.m ex05.mr<-tapply(resp, rec, mean) ex05.mr ex05.me<-tapply(resp, esp, mean) ex05.me ### Verificando interações par(mfrow=c(1,2)) interaction.plot(rec, esp, resp) interaction.plot(esp, rec, resp) ### Análise de variância ex05.av<-aov(resp ~ rec + esp + rec * esp) ### Resumindo o comando: ex05.av<-aov(resp ~ rec * esp) summary(ex05.av) # A função a seguir “model.tables” produz tabelas das médias definidas pelo modelo. ex05.mt <- model.tables(ex05.av, ty="means") ex05.mt summary(ex04.av) ### Análise de resíduos par(mfrow=c(2,2)) plot(ex05.av) ### NOrmalidade hist(ex05.av$res) shapiro.test(ex05.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(ex05.av$res~esp) boxplot(ex05.av$res~rec) boxplot(ex05.av$res~esp:rec) bartlett.test(ex05.av$res, esp) bartlett.test(ex05.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) ex05.avr <- aov(resp ~ rec/esp) #Para que o comando "summary" possa mostrar nossa analise adequadamente e necessario definir os efeitos de interesse. effects(ex05.avr) summary(ex05.avr, split=list("rec:esp"=list(esp.dentro.r1=1,esp.dentro.r2=2, esp.dentro.r3=3))) ex05.ave <- aov(resp ~ esp/rec) effects(ex05.ave) summary(ex05.ave, split=list("esp:rec"=list(rec.dentro.e1=c(1,3), rec.dentro.e2=c(2,4)))) require(agricolae) ### Tukey ex05.anova<-anova(ex05.av) ex05.anova ##-------------------------------------------------------------------------------------------## #Recipientes dentro de esp1 rec1<-subset(ex05,ex05$esp=="e1") HSD.test(rec1$resp,rec1$rec,ex05.anova$Df[4],ex05.anova$Mean[4],alpha=0.05) #Recipientes dentro de esp2 rec2<-subset(ex05,ex05$esp=="e2") HSD.test(rec2$resp,rec2$rec,ex05.anova$Df[4],ex05.anova$Mean[4],alpha=0.05) ##-------------------------------------------------------------------------------------------## #Especie dentro de recipiente1 esp1<-subset(ex05,ex05$rec=="r1") HSD.test(esp1$resp,esp1$esp,ex05.anova$Df[4],ex05.anova$Mean[4],alpha=0.05) #Especie dentro de recipiente2 esp2<-subset(ex05,ex05$rec=="r2") HSD.test(esp2$resp,esp2$esp,ex05.anova$Df[4],ex05.anova$Mean[4],alpha=0.05) #Especie dentro de recipiente3 esp3<-subset(ex05,ex05$rec=="r3") HSD.test(esp3$resp,esp3$esp,ex05.anova$Df[4],ex05.anova$Mean[4],alpha=0.05) ###---------------------------------------------------------------------------------------------------------------------#### ############################################### Esquema fatorial 2x2x2 ##################################################### ### Experimento no esquema fatorial com três fatores ### Adubação de cafeeiro ### O exemplo a seguir é de um delineamento em blocos completos casualizados com 6 repetições no esquema fatorial 2x2x2, no qual foi estudado o efeito da adubação mineral no cafeeiro, com Nitrogênio (N), Fósforo (P) e Potássio(K), sendo cada contendo 2 níveis: 0 (ausência) e 1 (presença). Cada parcela era constituída por 12 covas no espaçamento 3,5 x 2,5, com uma área de 105 m2, sendo que os dados obtidos foram transformados para Kg/ha.extraido "Experimentação Agricola" Kronha & Banzatto pg 112. rm(list=ls()) ex06<-read.table("exemplo06.txt",header=T) names(ex06) summary(ex06) is.factor(ex06$n) is.factor(ex06$p) is.factor(ex06$k) is.factor(ex06$bloc) ### definindo os fatores ex06$resp<-as.numeric(ex06$resp) ex06$n<-as.factor(ex06$n) ex06$p<-as.factor(ex06$p) ex06$k<-as.factor(ex06$k) ex06$bloc<-as.factor(ex06$bloc) summary(ex06) attach(ex06) ### Explorando os dados ex06.mn <- tapply(resp, n, mean) ex06.mn ex06.mp <- tapply(resp, p, mean) ex06.mp ex06.mk <- tapply(resp, k, mean) ex06.mk ex06.mNPK <- aggregate(resp,list(N=n,P=p,K=k), mean) ex06.mNPK par(mfrow=c(2,2)) interaction.plot(n, p, resp) title("Interação N vs P") interaction.plot(n, k, resp) title("Interação N vs K") interaction.plot(p, k, resp) title("Interação P vs K") par(mfrow=c(2,2)) boxplot(resp~n,main="Nitrogenio") boxplot(resp~k,main="Potassio") boxplot(resp~p,main="Fosforo") ### ANOVA ### Para o experimento fatorial conduzido no delineamento em blocos completos casualizados tem o seguinte comando: ex06.av <- aov(resp ~ bloc + n + p + k + n*p + n*k + p*k + n*p*k) ### ou ex06.av <- aov(resp ~ bloc + n*p*k) summary(ex06.av) ### O uso da função model.tables permite a obtenção de estatísticas para cada um dos fatores do modelo e a combinação entre eles. ex06.mt = model.tables(ex06.av, ty="means") ex06.mt ### Verificando os pressuposto par(mfrow=c(2,2)) plot(ex06.av) ### Normalidade hist(ex06.av$res) shapiro.test(ex06.av$res) ### Homocedacicidade par(mfrow=c(2,2)) boxplot(ex06.av$res~n,main="Nitrogenio") boxplot(ex06.av$res~k,main="Potassio") boxplot(ex06.av$res~p,main="Fosforo") boxplot(ex06.av$res~n*p*k,main="Interação n*p*k") ### Teste de Bartlett bartlett.test(ex06.av$res, n) bartlett.test(ex06.av$res, p) bartlett.test(ex06.av$res, k) ### Os pressupostos do modelo de análise de variância foram todos atendidos ### Portanto pode ser observado que o efeito da interação NK foi significativo. Portanto efetua-se o desdobramento dos graus de liberdade da interação. ### Análise da Interação par(mfrow=c(2,2)) interaction.plot(n, p, resp) title("Interação N vs P") interaction.plot(k, n, resp) title("Interação N vs K") interaction.plot(p, k, resp) title("Interação P vs K") ### Desdobramento da interação NK: interaction.plot(k, n, resp,main="Interação N vs K") ##-------------------------------------------------------------------------------------------## ### Inicialmente, deve-se utilizar o modelo com o efeito de N aninhado em K. ex06.avKN <- aov(resp~bloc+k/n+n*p*k-n) effects(ex06.avKN) ### Observe que apenas nos interessa interpretar os efeitos do desdobramento. summary(ex06.avKN,split=list("k:n"=list(N.dentro.K0=1, N.dentro.K1=2))) ### Inicialmente, deve-se utilizar o modelo com o efeito de N aninhado em K. ex06.avNK <- aov(resp~bloc+n/k+n*p*k-k) effects(ex06.avNK) ### Observe que apenas nos interessa interpretar os efeitos do desdobramento. summary(ex06.avNK,split=list("n:k"=list(K.dentro.N0=1, K.dentro.N1=2))) ##-------------------------------------------------------------------------------------------## ### Desdobramento da interação NP: par(mfrow=c(2,1)) interaction.plot(n, p, resp,main="Interação N vs P") interaction.plot(p, n, resp,main="Interação P vs N") ##-------------------------------------------------------------------------------------------## ### Inicialmente, deve-se utilizar o modelo com o efeito de N aninhado em P. ex06.avPN <- aov(resp~bloc+p/n+n*p*k-n) ### Observe que apenas nos interessa interpretar os efeitos do desdobramento. summary(ex06.avPN,split=list("p:n"=list(N.dentro.P0=1, N.dentro.P1=2))) ### Inicialmente, deve-se utilizar o modelo com o efeito de N aninhado em P. ex06.avNP <- aov(resp~bloc+n/p+n*p*k-p) ### Observe que apenas nos interessa interpretar os efeitos do desdobramento. summary(ex06.avNP,split=list("n:p"=list(P.dentro.N0=1, P.dentro.N1=2))) ##-------------------------------------------------------------------------------------------## ### Desdobramento da interação PK: par(mfrow=c(2,1)) interaction.plot(k, p, resp,main="Interação K vs P") interaction.plot(p,k, resp,main="Interação P vs K") ##-------------------------------------------------------------------------------------------## ### Inicialmente, deve-se utilizar o modelo com o efeito de K aninhado em P. ex06.avPK <- aov(resp~bloc+p/k+n*p*k-k) ### Observe que apenas nos interessa interpretar os efeitos do desdobramento. summary(ex06.avPK,split=list("p:k"=list(K.dentro.P0=1, K.dentro.P1=2))) ### Inicialmente, deve-se utilizar o modelo com o efeito de K aninhado em P. ex06.avKP <- aov(resp~bloc+k/p+n*p*k-p) ### Observe que apenas nos interessa interpretar os efeitos do desdobramento. summary(ex06.avKP,split=list("k:p"=list(P.dentro.K0=1, P.dentro.K1=2))) ##-------------------------------------------------------------------------------------------## ### Teste de Tukey para comparações múltiplas ex06.tk = TukeyHSD(ex06.avn, "k:n",ord=T) ex06.tk ### Tukey ex06.anova<-anova(ex06.av) ex06.anova #Execuraremos os teste de Tukey somente para a interação significativa, mais o procedimento e similar para as outras interações duplas require(agricolae) ##-----------------------------------------------------------------------## ## Tukey de K:N inteNK<-tapply(resp,list(N=n,K=k),mean) inteNK N0<-subset(ex06,ex06$n=="0") HSD.test(N0$resp,N0$k,ex06.anova$Df[9],ex06.anova$Mean[9],alpha=0.05) N1<-subset(ex06,ex06$n=="1") HSD.test(N1$resp,N1$k,ex06.anova$Df[9],ex06.anova$Mean[9],alpha=0.05) K0<-subset(ex06,ex06$k=="0") HSD.test(K0$resp,K0$n,ex06.anova$Df[9],ex06.anova$Mean[9],alpha=0.05) K1<-subset(ex06,ex06$k=="1") HSD.test(K1$resp,K1$n,ex06.anova$Df[9],ex06.anova$Mean[9],alpha=0.05) ##-------------------------------------------------------------------------------------------## # Em nosso exeomplo nossa interação tripla não foi significativa, porem, caso a tripla seja significativa, uma das alternativas e proceder conforme os comandos abixo onde estudamos as interação entre N e P dentro do K, a ordem dos fatores poderiam ser alterados: ## Tukey de N:P:K inteNPK<-tapply(resp,list(N=n,P=p,K=k),mean) inteNPK ##-----------------------------------------------------------------------## ##Dentro de K0 ##Compararando as medias de P dentro de N0KO K0N0<-subset(ex06,ex06$n=="0" & ex06$k=="0") HSD.test(K0N0$resp,K0N0$p,ex06.anova$Df[9],ex06.anova$Mean[9],alpha=0.05) ##Compararando as medias de P dentro de N1K0 K0N1<-subset(ex06,ex06$n=="1" & ex06$k=="0") HSD.test(K0N1$resp,K0N1$p,ex06.anova$Df[9],ex06.anova$Mean[9],alpha=0.05) ##Compararando as medias de N dentro de P0K0 K0P0<-subset(ex06,ex06$p=="0" & ex06$k=="0") HSD.test(K0P0$resp,K0P0$n,ex06.anova$Df[9],ex06.anova$Mean[9],alpha=0.05) ##Compararando as medias de N dentro de P1K0 K0P1<-subset(ex06,ex06$p=="1" & ex06$k=="0") HSD.test(K0P1$resp,K0P1$n,ex06.anova$Df[9],ex06.anova$Mean[9],alpha=0.05) ##-----------------------------------------------------------------------## ##Dentro de K1 ##Compararando as medias de P dentro de N0K1 K0N0<-subset(ex06,ex06$n=="0" & ex06$k=="1") HSD.test(K0N0$resp,K0N0$p,ex06.anova$Df[9],ex06.anova$Mean[9],alpha=0.05) ##Compararando as medias de P dentro de N1K1 K0N1<-subset(ex06,ex06$n=="1" & ex06$k=="1") HSD.test(K0N1$resp,K0N1$p,ex06.anova$Df[9],ex06.anova$Mean[9],alpha=0.05) ##Compararando as medias de N dentro de P0K1 K0P0<-subset(ex06,ex06$p=="0" & ex06$k=="1") HSD.test(K0P0$resp,K0P0$n,ex06.anova$Df[9],ex06.anova$Mean[9],alpha=0.05) ##Compararando as medias de N dentro de P1K1 K0P1<-subset(ex06,ex06$p=="1" & ex06$k=="1") HSD.test(K0P1$resp,K0P1$n,ex06.anova$Df[9],ex06.anova$Mean[9],alpha=0.05) ###---------------------------------------------------------------------------------------------------------------------#### ############################################### Parcela sub-dividida ####################################################### ### Experimentos em Parcelas Subdivididas ### 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. rm(list=ls()) ex07<- read.table("exemplo07.txt", header=T) ex07 summary(ex07) ex07$a <- as.factor(ex07$a) ex07$b <- as.factor(ex07$b) ex07$bloco <- as.factor(ex07$bloco) summary(ex07) attach(ex07) ###Análise de variância ex07.av <- aov(resp ~ bloco + a*b + Error(bloco:a)) summary(ex07.av) ### Visualização gráfica do comportamento de variedades dentro de tratamentos de sementes e vice-versa. par(mfrow=c(1,2)) interaction.plot(a, b, resp, xlab="Variedades") interaction.plot(b, a, resp,xlab="Trat. de Sementes") ### Em experimento de parcela subdividida, 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: ex07.avb <- aov(resp ~ bloco + a*b + bloco*a) summary(ex07.avb) ### As médias de cada fator podem ser obtidas com o comando model.tables ex07.m <- model.tables(ex07.av, ty="means") ex07.m ### Análise de resíduos par(mfrow=c(2,2)) plot(ex07.avb) ### Normalidade shapiro.test(ex07.avb$res) ### Homocedacicidade par(mfrow=c(2,2)) boxplot(ex07.avb$res~b) boxplot(ex07.avb$res~a) boxplot(ex07.avb$res~a:b) # O bartlet pode ser executado mais sempre lembramos de suas restrições bartlett.test(ex07.avb$res, b) bartlett.test(ex07.avb$res, a) ##-----------------------------------------------------------------------## #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 "a" ex07.anovaRA<-anova(ex07.av$'bloco:a') ex07.anovaRA #extraindo os QM e DF de "a" QMEa<-ex07.anovaRA$Mean[2] QMEa Dfea<-ex07.anovaRA$Df[2] Dfea HSD.test(resp,a,Dfea,QMEa,alpha=0.05) ##-----------------------------------------------------------------------## ###Fator secundario "b" ex07.anovaRB<-anova(ex07.av$'Within') ex07.anovaRB #extraindo os QM e DF de "a" QMEb<-ex07.anovaRB$Mean[2] QMEb Dfeb<-ex07.anovaRB$Df[2] Dfeb HSD.test(resp,b,Dfeb,QMEb,alpha=0.05) #Desdobrando e testanto a interação "a:b" ##Primeiro vamos estudar o fatores secundarios dentro dos principais interaction.plot(a, b, resp, xlab="Variedades") ex07.avISP<-aov(resp ~ bloco +a/b + b*a -b + Error(bloco:a)) summary(ex07.avISP) effects(ex07.avISP$'Within') summary(ex07.avISP,split=list("a:b"=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(ex07$resp,list(Aveia=ex07$a,Tratamento=ex07$b),mean) ## b dentro de A1 A1<-subset(ex07,ex07$a=="1") HSD.test(A1$resp,A1$b,DFRb,QMRb) ## b dentro de A2 A2<-subset(ex07,ex07$a=="2") HSD.test(A2$resp,A2$b,DFRb,QMRb) ## b dentro de A3 A3<-subset(ex07,ex07$a=="3") HSD.test(A3$resp,A3$b,DFRb,QMRb) ## b dentro de A4 A4<-subset(ex07,ex07$a=="4") HSD.test(A4$resp,A4$b,DFRb,QMRb) ##-----------------------------------------------------------------------## ##Agora vamos estudar a interação dos efeitos principais dentro dos primarios interaction.plot(b, a, resp, xlab="Tratamentos de Sementes") ex07.avIPS<-aov(resp ~ bloco +b/a + b*a -a) summary(ex07.avIPS) effects(ex07.avIPS) summary(ex07.avIPS,split=list("b:a"=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 esta 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(ex07.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) ) ##Com o novo erro e Gl podemos agora aplicar o teste F adequadamente e realizar também as comparações multiplas interaction.plot(b, a, resp, xlab="Tratamentos de Sementes") tapply(ex07$resp,list(Aveia=ex07$a,Tratamento=ex07$b),mean) ## a dentro de b1 B1<-subset(ex07,ex07$b=="1") HSD.test(B1$resp,B1$a,Glx,QMex) ## a dentro de B2 B2<-subset(ex07,ex07$b=="2") HSD.test(B2$resp,B2$a,Glx,QMex) ## a dentro de B3 B3<-subset(ex07,ex07$b=="3") HSD.test(B3$resp,B3$a,Glx,QMex) ## a dentro de B4 B4<-subset(ex07,ex07$b=="4") HSD.test(B4$resp,B4$a,Glx,QMex) ###---------------------------------------------------------------------------------------------------------------------#### #EXERCICIO: No arquivo "exemplo07E.txt", estão os dados retirados do livro "Experimentação Agrícola" Kronha & Banzatto, este experimento #foi conduzido em bloco ao acaso em parcela subdivididas, o tratamento principal foi adubos: #A1-Salitre do Chile #A2-Sulfato de amônio #A3-Uréia #A4-Calnitro #Como tratamento secundario foi adotado 3 Doses: #D1- 10 kg de N/ha #D2- 40 kg de N/ha #D3- 70 kg de N/ha #Estes tratamentos foram avaliados na cultura do milho, onde como variavel resposta foi mensurado a produtividade em kg/ha. #Obtenha a ANOVA,desdobramento da ANOVA, testes de comparação multiplas e conclusões sobre o experimento. ###---------------------------------------------------------------------------------------------------------------------#### ###---------------------------------------------------------------------------------------------------------------------#### ################################################ Delineamento em faixa ##################################################### ### O exemplo a seguir e retirados do livro "Experimentação Agrícola" Kronha & Banzatto pg 161, o experimento avaliou t/acre de # beteraba submetida a 0,80,160 e 320 libras/acre de N e 5 épocas de plantio conduzidos em um experimento em faixas. rm(list=ls()) ex08<- read.table("exemplo08.txt", header=T) ex08 ### Explorando os dados names(ex08) summary(ex08) ex08$nitro <- as.factor(ex08$nitro) ex08$epoca <- as.factor(ex08$epoca) ex08$bloco<- as.factor(ex08$bloco) summary(ex08) attach(ex08) hist(resp) ###Interação par(mfrow=c(1,2)) interaction.plot(epoca, nitro, resp) interaction.plot(nitro, epoca, resp) ### Análise de variância ex08.av<- aov(resp ~ bloco + nitro*epoca + Error(bloco:nitro+bloco:epoca)) summary(ex08.av) #Esse modelo armazena os resíduos (Note que este resido e a soma dos três residuos do modelo no objeto ex08.av ex08.av1= aov(resp ~ bloco + nitro*epoca,data=ex08) summary(ex08.av1) ### Análise de resíduos (geral) par(mfrow=c(2,2)) plot(ex08.av1) ### Normalidade shapiro.test(ex08.av1$res) ### Homocedacicidade bartlett.test(ex08.av1$res,ex08$nitro) bartlett.test(ex08.av1$res,ex08$epoca) par(mfrow=c(2,2)) plot(ex08$nitro,ex08.av1$res) plot(ex08$epoca,ex08.av1$res) plot(ex08$nitro:ex08$epoca,ex08.av1$res) ##-----------------------------------------------------------------------## ### Comparação Multipla ### Note que nenhum fator em estudo e significativo, portanto poderiamos parar a analise, porem apenas de forma didatica vamos realizar as comparações ### Tukey summary(ex08.av) names(ex08.av) ### Vamos testar todos os fatores, porem pela anova deveriamos testar apenas a interação ex08.anovaRA<-anova(ex08.av$'bloco:nitro') ex08.anovaRA ### Testando fator nitro HSD.test(resp,nitro,ex08.anovaRA$Df[2],ex08.anovaRA$Mean[2],alpha=0.05) ##-----------------------------------------------------------------------## ### Testando fator epoca ex08.anovaRB<-anova(ex08.av$'bloco:epoca') ex08.anovaRB ### Fator b HSD.test(resp,epoca,ex08.anovaRB$Df[1],ex08.anovaRB$Mean[1],alpha=0.05) ##-----------------------------------------------------------------------## ### Testando interação Epoca:Plantio tapply(ex08$resp,list(epoca=ex08$epoca,nitro=ex08$nitro),mean) ###Comparação de nitrogênio em cada epoca summary(ex08.av) QMa<-19.79 Dfa<-9 QMb<-2.3 Dfb<-12 QMc<-3.979 Dfc<-36 ##Tratamentos a<-4 b<-5 #Erro medio a e c QMeac<-(QMa+(b-1)*QMc)/b QMeac ##Graus de liberdade a e c Glac<-((QMa+(b-1)*QMc)^2) / ( ((QMa^2)/Dfa) + ((((b-1)*QMc)^2)/Dfc) ) Glac ##-----------------------------------------------------------------------## ## Nitrogênio dentro da E1 E1<-subset(ex08,ex08$epoca=="E1") HSD.test(E1$resp,E1$nitro,Glac,QMeac) ## Nitrogênio dentro da E2 E2<-subset(ex08,ex08$epoca=="E2") HSD.test(E2$resp,E2$nitro,Glac,QMeac) ## Nitrogênio dentro da E3 E3<-subset(ex08,ex08$epoca=="E3") HSD.test(E3$resp,E3$nitro,Glac,QMeac) ## Nitrogênio dentro da E4 E4<-subset(ex08,ex08$epoca=="E4") HSD.test(E4$resp,E4$nitro,Glac,QMeac) ## Nitrogênio dentro da E5 E5<-subset(ex08,ex08$epoca=="E5") HSD.test(E5$resp,E5$nitro,Glac,QMeac) ##-----------------------------------------------------------------------## #Erro medio b e c QMebc<-(QMb+(a-1)*QMc)/a QMebc ##Graus de liberdade a e c Glbc<-((QMb+(a-1)*QMc)^2) / ( ((QMb^2)/Dfb) + ((((a-1)*QMc)^2)/Dfc) ) Glbc ## Época dentro dos niveis de N0 N0<-subset(ex08,ex08$nitro=="N0") HSD.test(N0$resp,N0$epoca,Glbc,QMebc) ## Época dentro dos niveis de N1 N1<-subset(ex08,ex08$nitro=="N1") HSD.test(N1$resp,N1$epoca,Glbc,QMebc) ## Época dentro dos niveis de N2 N2<-subset(ex08,ex08$nitro=="N2") HSD.test(N2$resp,N2$epoca,Glbc,QMebc) ## Época dentro dos niveis de N3 N3<-subset(ex08,ex08$nitro=="N3") HSD.test(N3$resp,N3$epoca,Glbc,QMebc) ###---------------------------------------------------------------------------------------------------------------------#### #EXERCICIO:No arquivo "exemplo08E.txt" existem dados oriundos da produção de milho em kg/parcela em diferentes espaçamentos #e diferentes densidades, sendo o experimento conduzido em faixas. obtenha a ANOVA e as comparações caso for necessario. ###---------------------------------------------------------------------------------------------------------------------#### ################################################ Regressão linear ########################################################## ### Dados retirados da apostila do Prof° Daniel Furtado (UFLA) ### disponivel em: http://www.dex.ufla.br/~danielff/RRC0.pdf ### O dados referem-se sobre o crescimento de plantas (y) versus tempo de exposição a luz solar (x), ### em unidades apropiadas ### digitamos os dados mais poderiam ter sido carregados de um arquivo com fizemos em experimentos conjuntos, ### Usando a função "read.table()". rm(list=ls()) x<-c(0.1,0.2,0.3,0.5,0.8,1,1.5,2) y<-c(0.88,0.9,0.99,1.12,1.4,1.62,2.2,3.10) ### Explorando Visualmente os dados plot(x,y,xlab="Tempo Exposição Solar em hora",ylab="Crescimento",main="Exposição solar x Crescimento") ### 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 todos os ajustes lineares ajuste1<-lm(y~x) summary(ajuste1) names(ajuste1) ### Para verificar se o modelo é significativo utilizamos a função "anova", que neste caso testa se o modelo proposto ### explica ou não a variação dos dados anova(ajuste1) par(mfrow=c(2,2)) plot(ajuste1) ### Analise da normalidade dos resíduos shapiro.test(ajuste1$resi) ### Gerando Graficos plot(x,y,xlab="Tempo Exposição Solar em hora",ylab="Crescimento",main="Exposição solar x Crescimento") abline(ajuste1) ### Ajuste2 ajuste de uma função polinomial y=a+b*x+c*x^2, lembrando que isso também é uma FUNÇÃO LINEAR!!!!!!!!!! ### 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) ### NOrmalidade dos residuos shapiro.test(ajuste2$resi) ###Gerando Graficos xx<-seq(min(x),max(x),by=0.01) yy<-predict(ajuste2,list(x=xx)) plot(x,y,xlab="Tempo Exposição Solar em hora",ylab="Crescimento",main="Exposição solar x Crescimento") lines(xx,yy) ### limpando a memoria rm(list=ls()) ### Regressão linear Multivarida ### Dados relativos a araucaria angustifolia, mensurado o volume m³/acre Y), área Basal dm²(X1), ### area basal relativa % outra especie (X1), altura pés (X3) ### Fonte: Apostila de Daniel Furtado Ferreira UFLA Y<-c(65,78,82,86,87,90,93,96,104,113) X1<-c(41,71,90,80,93,90,87,95,100,101) X2<-c(79,48,80,81,61,70,96,84,78,96) X3<-c(35,53,64,59,66,64,62,67,70,71) arvore<-data.frame(Y,X1,X2,X3) plot(arvore) ### Verificando a estrutura de correlação entre as variaveis cor(arvore) ### Na sequencia ajustaremos um modelo para cada variavel de forma independente ### ajuste X1 ajuste.x1<-lm(Y~X1) summary(ajuste.x1) anova(ajuste.x1) par(mfrow=c(2,2)) plot(ajuste.x1) shapiro.test(ajuste.x1$res) ### ajuste X2 ajuste.x2<-lm(Y~X2) summary(ajuste.x2) anova(ajuste.x2) par(mfrow=c(2,2)) plot(ajuste.x2) shapiro.test(ajuste.x2$res) ### ajuste X3 ajuste.x3<-lm(Y~X3) summary(ajuste.x3) anova(ajuste.x2) par(mfrow=c(2,2)) plot(ajuste.x3) shapiro.test(ajuste.x3$res) ### Ajuste multivariado ### Para ajuste multivariados simplismente adicionamos as variaveis "X1 + X2 ... Xn-1 + Xn" no modelos "lm(Y~X1+X2...+Xn-1+Xn)" ajuste.m<-lm(Y~X1+X2+X3) summary(ajuste.m) anova(ajuste.m) ### Relembrando a estrutura de correlação cor(arvore) ### Analisando os residuos do ajuste par(mfrow=c(2,2)) plot(ajuste.m) shapiro.test(ajuste.m$res) ### Após analisar decidimos retirar a variavel X3(altura), do modelo, pelo sua forte correlação e pensando em motivos praticos no campo ajuste.x1x2<-lm(Y~X1+X2) summary(ajuste.x1x2) anova(ajuste.x1x2) ### Analisando os residuos par(mfrow=c(2,2)) plot(ajuste.x1x2) shapiro.test(ajuste.x1x2$res) ### plotando em 3d, para isso utilizaremos o pacote "scatterplot3d", caso não tenha ele instalado: #install.packages("scatterplot3d",dep=T) ### Requerendo o pacote para utilizar suas funções require(scatterplot3d) scatterplot3d(X1,X2,Y) ### Gerando series para plot a função ajustada xx1 <- seq(min(X1),max(X1),length=50) xx2 <- seq(min(X2),max(X2),length=50) ### Coeficientes para as funções names(ajuste.x1x2) ajuste.x1x2$coef ### com "function()" podemos excrever funções dentro do R f <- function(xx1,xx2) { r <-ajuste.x1x2$coef[1]+ajuste.x1x2$coef[2]*xx1+ajuste.x1x2$coef[3]*xx2 return(r) } yy <- outer(xx1,xx2,f) yy ### teste escrever "f(41,48)" e verifique o resultado a<-f(41,48) a ### Plotando os dados em pespectiva persp(xx1,xx2,yy,theta = 30, phi = 30,expand =0.7,col = "blue", shade = 0.75, ticktype = "detailed", xlab = "X1", ylab = "X2", zlab = "Y") ### Visualizando de outras formas image(xx1,xx2,yy) filled.contour(xx1,xx2,yy,color.palette = colorRampPalette(c("red","white","blue"))) contour(xx1,xx2,yy,nlevels=15) require(scatterplot3d) ### Bem legal a visualização s3d <- scatterplot3d(X1,X2,Y,type="h",highlight.3d=TRUE,angle=55, scale.y =0.7, pch=16) s3d$plane3d(ajuste.x1x2, col=4) ### OBS: Veja o help das funções para mais detalhes.... tem coisas legais ####---------------------------------------------------------------------------------------------------------------------#### ################################################ Regressão não linear ###################################################### ### Fonte: Apostila Prof° Paulo Justiniano UFPR ### a Equação de Van Genutchen, importante para caracterizar os aspestos do comportamento hidrico do solo ### pot significa potencial matrico do solo ### u significa umidade do solo g/g rm(list=ls()) pot<-c(10,19,30,45,63,64,75,89,105,138,490,3000,4100,5000,26300) u<-c(0.3071,0.2931,0.2828,0.2753,0.2681,0.2628,0.2522,0.2404,0.2272,0.2120,0.1655,0.1468,0.1205,0.1013,0.0730) plot(pot,u) plot(log10(pot),u, xlab=expression(log[10](psi[m])), ylab=expression(theta)) ### Utilizaremos a função "nls" para Ajuste não lineares help(nls) ### Realizando os ajustes ajuste<-nls(u~ur+(us-ur)/((1+(alpha*pot)^n)^(1-1/n)),start=list(us=0.2236,ur=0.0611,alpha=0.056,n=1.5351),trace=T) summary(ajuste) ### A função não calcula o R² de forma automatica, precisamos calcularmos de uma forma manual ### R^2 ### Soma de quadrados residual SQE <- summary(ajuste)$sigma^2*summary(ajuste)$df[2] SQE ### Soma de quadrado total corrigida SQT <- var(u)*(length(u)-1) SQT R2 <- 1 - SQE/SQT R2 ### Verificando os residuos resi<-u-predict(ajuste) resi hist(resi) plot(resi) shapiro.test(resi) ###PLotando xx<-seq(min(pot),max(pot),l=100) plot(log10(pot),u, xlab=expression(log[10](psi[m])), ylab=expression(theta), main="Equação de Van Genutchen") lines(log10(xx),predict(ajuste,list(pot=xx))) rm(list=ls()) ############################################################################################################# ### Uma pequena brincadeira... xx1 <- seq(-2*pi,2*pi,length=80) xx2 <- seq(-2*pi,2*pi,length=80) for(t in seq(1,50,by=0.3)) { f <- function(xx1,xx2) { r <-sin(((t-xx1^2)+(t-xx2^2)^0.8)) r } yy <- outer(xx1,xx2,f) persp(xx1,xx2,yy,theta = 30, phi = 35,expand =0.15,col =2,ticktype = "detailed", xlab = "x", ylab = "y", zlab = "z") } ### Ababou, fim.....