########################################################################################################################### # Curso de Estatistica Experimental com Software R -- XXXII Ciclo de Atuzlização em Ciências Agrárias # Éder David Borges da Silva / Renato Gonçalves de Oliveira # 08/05/2010 ########################################################################################################################### ###---------------------------------------------------------------------------------------------------------------------#### ####################################### Manipulação de Banco de Dados ###################################################### ###---------- RODBC -----------### ### OBS: RODBC não funciona em LINUX ### Manipulando dados em formato mdb (ACESS) require(RODBC) con <- odbcConnectAccess("CHUVAS.mdb") con tabelas <- sqlTables(con)$"TABLE_NAME" tabelas ### extraindo uma tabela COMPLETAMENTE chuvas <- sqlFetch(con, "Chuvas") head(chuvas) class(chuvas) odbcCloseAll() ### Manipulando dados em formato xls (Excel) con <- odbcConnectExcel("tripla.xls") con plan <- sqlTables(con)$"TABLE_NAME" plan dados <- sqlFetch(con,"tripla") head(dados) summary(dados) str(dados) odbcClose(con) ###---------- GDATA -----------### #colocar exemplo require(gdata) dados <- read.xls("tripla.xls",sheet="tripla") ### Além destes formatos de BD, o R tem suporte a BD geograficos e bancos mais complexos com SQL, ORACLE, etc.. ###---------------------------------------------------------------------------------------------------------------------#### ##################################################### 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) search() ex02.m <-tapply(prod, trat, mean) ex02.m ex02.v <-tapply(prod, trat, var) ex02.v plot(ex02) 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) ex02.avSB<-aov(prod~trat) summary(ex02.avSB) #Obtendo o CV ex02.anova<-anova(ex02.av) CV<-((ex02.anova$Mean[3]^0.5)/(model.tables(ex02.av,"means")$tables$Gran))*100 CV #Verificando os residuos graficamente par(mfrow=c(2,2)) plot(ex02.av) ###Homocedaciade, testando através do teste de Bartlett bartlett.test(ex02.av$res~trat) ###Normalidade Verificada através do teste de Shapiro-Wilk 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) # Ou tuk<-HSD.test(prod,trat,42,71476.3,alpha=0.05) #Gravando um resultado sink("ANalisebloco.txt") tuk sink() 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) detach(ex02) ###---------------------------------------------------------------------------------------------------------------------#### ############################################## 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) # Verificar residuos graficamente par(mfrow=c(2,2)) plot(ex03.av) ##Homocedasticidade bartlett.test(ex03$resp, ex03$trat) ##Normalidade 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 ## Comparando dados transformados com originais par(mfrow=c(1,2)) hist(ex03$resp) hist(ex03$respt) attach(ex03) # Realizando a ANOVA com os dados transformados ex03.avt <- aov(respt ~ trat) anova(ex03.avt) ### Analisando os residuos graficamente plot(ex03.avt) ### Comparando os dois residuos graficamente par(mfrow=c(4,2)) plot(ex03.avt) plot(ex03.av) ##Homocedasticidade plot.default(ex03$trat,ex03.avt$res) bartlett.test(ex03$respt, ex03$trat) ##Normalidade qqnorm(ex03.avt$res,ylab="Residuos", main="Grafico Normal de Probabilidade dos Resíduos") qqline(ex03.avt$res) 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 detach(ex03) ###---------------------------------------------------------------------------------------------------------------------#### ############################################### Esquema fatorial 3x2x2 ##################################################### ### Experimento retirados do livro do banzato e Kronha 'Experimentação Agricola' pg120. e diferentes doses de ### calcario (Ca) com diferentes doses de fosforo com três cultivares de trigo. dados<-read.csv("tripla.csv",header=T,sep=";") names(dados) summary(dados) attach(dados) ### Explorando a Interação par(mfrow=c(2,2)) interaction.plot(cultivares,calcio, resp,main='Calcio dentro das Cultivares') interaction.plot(cultivares,fosforo, resp,main='Fosforo dentro das Cultivares') interaction.plot(calcio,fosforo, resp,main='Fosforo dentro de Calcio') ### Modelo de ANOVA dados.av <- aov(resp~cultivares*calcio*fosforo) summary(dados.av) ###Convertendo para tabela da ANOVA dados.anova <- anova(dados.av) ### Verificando os pressuposto graficamente par(mfrow=c(2,2)) plot(dados.av) ### Normalidade hist(dados.av$res) shapiro.test(dados.av$res) ######-------------------------------------------------------------------------- ## Desdobrando as Duplas (Fosforo dentro de Cultivares) interaction.plot(cultivares,fosforo, resp,main='Fosforo dentro das Cultivares') dados.avCuFo <- aov(resp~cultivares/fosforo+cultivares*calcio*fosforo-fosforo) summary(dados.avCuFo) effects(dados.avCuFo) ### Observe que apenas nos interessa interpretar os efeitos do desdobramento. summary(dados.avCuFo,split=list("cultivares:fosforo"=list( P.dentro.Cul1=1, P.dentro.Cul2=2, P.dentro.Cul3=3))) ######-------------------------------------------------------------------------- ## Desdobrando as Duplas (Calcio dentro de Cultivares) interaction.plot(cultivares,calcio, resp,main='Calcio dentro das Cultivares') dados.avCuCa <- aov(resp~cultivares/calcio+cultivares*calcio*fosforo-calcio) summary(dados.avCuCa) effects(dados.avCuCa) ### Observe que apenas nos interessa interpretar os efeitos do desdobramento. summary(dados.avCuCa,split=list("cultivares:calcio"=list( Ca.dentro.Cul1=1, Ca.dentro.Cul2=2, Ca.dentro.Cul3=3))) ######--------------------------------------------------------------------------## ### Desdrobramento da interação tripla ### Inicialmente, cultivares dentro de combinações calcio e fosforo dados.avCF <- aov(resp~calcio*fosforo/cultivares) summary(dados.avCF) effects(dados.avCF)##Verifica posição dento das interações triplas para informar no desdobramento summary(dados.avCF,split=list("calcio:fosforo:cultivares"=list( Cultivares.d.C0P0=c(1,5), Cultivares.d.C0P1=c(3,7), Cultivares.d.C1P0=c(2,6), Cultivares.d.C1P1=c(4,8)))) inteCF<-tapply(resp,list(cultivares,calcio,fosforo),mean) inteCF # P0 P1 # Ca0 Ca1 Ca0 Ca1 # C1 1211.00 704.5 530.0 427.5 # C2 1491.75 983.5 586.0 405.0 # C3 1664.00 777.5 635.5 489.5 ##-----------------------------------------------------------------------## require(agricolae) ##Cultivares dentro de Ca0P0 Ca0P0<-subset(dados,dados$calcio=="Ca0" & dados$fosforo=="P0") Ca0P0 HSD.test(Ca0P0$resp,Ca0P0$cultivares,dados.anova$Df[8],dados.anova$Mean[8],alpha=0.05) ##Cultivares dentro de Ca0P1 Ca0P1<-subset(dados,dados$calcio=="Ca0" & dados$fosforo=="P1") Ca0P1 HSD.test(Ca0P1$resp,Ca0P1$cultivares,dados.anova$Df[8],dados.anova$Mean[8],alpha=0.05) ##Cultivares dentro de Ca1P0 Ca1P0<-subset(dados,dados$calcio=="Ca1" & dados$fosforo=="P0") Ca1P0 HSD.test(Ca1P0$resp,Ca1P0$cultivares,dados.anova$Df[8],dados.anova$Mean[8],alpha=0.05) ##Cultivares dentro de Ca1P1 Ca1P1<-subset(dados,dados$calcio=="Ca1" & dados$fosforo=="P1") Ca1P1 HSD.test(Ca1P1$resp,Ca1P1$cultivares,dados.anova$Df[8],dados.anova$Mean[8],alpha=0.05) ######--------------------------------------------------------------------------######################## ### depois, Calcio dentro de combinações cultivares e fosforo dados.avCulP <- aov(resp~cultivares*fosforo/calcio) summary(dados.avCulP) effects(dados.avCulP) summary(dados.avCulP,split=list("cultivares:fosforo:calcio"=list( Calcio.d.C1P0=1, Calcio.d.C1P1=4, Calcio.d.C2P0=2, Calcio.d.C2P1=5, Calcio.d.C3P0=3, Calcio.d.C3P1=6))) inteCulP<-tapply(resp,list(calcio,cultivares,fosforo),mean) inteCulP ### Para tukey seguir os procedimentos adotados para o outro desdobramento, mas note que pelos contrastes já se tira as conclusões ######--------------------------------------------------------------------------######################## ### depois, Calcio dentro de combinações cultivares e fosforo dados.avCC <- aov(resp~cultivares*calcio/fosforo) summary(dados.avCC) effects(dados.avCC) summary(dados.avCC,split=list("cultivares:calcio:fosforo"=list( Fosforo.d.C1Ca0=1, Fosforo.d.C1Ca1=4, Fosforo.d.C2Ca0=2, Fosforo.d.C2Ca1=5, Fosforo.d.C3Ca0=3, Fosforo.d.C3Ca1=6))) inteCulP<-tapply(resp,list(fosforo,calcio,cultivares),mean) inteCulP ### Para tukey seguir os procedimentos adotados para o outro desdobramento, mas note que pelos contrastes já se tira as conclusões #### ------------------------------------------------------------------------ ### Comparação via ScottKnott require(ScottKnott) head(dados) ### Fator Cultivares Cul.SK <- SK(x=dados[,1:3], y=dados[,4], model='resp~cultivares*calcio*fosforo',which='cultivares') summary(Cul.SK) plot(Cul.SK) ### Fator Calcio Ca.SK <- SK(x=dados[,1:3], y=dados[,4], model='resp~cultivares*calcio*fosforo',which='calcio') summary(Ca.SK) plot(Ca.SK) ### Fator fosforo P.SK <- SK(x=dados[,1:3], y=dados[,4], model='resp~cultivares*calcio*fosforo',which='fosforo') summary(P.SK) plot(P.SK) ### Duplas # # Ca0 Ca1 # P0 1455.5833 821.8333 # P1 583.8333 440.6667 ## Calcio:Fosforo CaP.SK1 <- SK.nest(x=dados[,1:3], y=dados[,4], model='resp~cultivares*calcio*fosforo',which='calcio:fosforo',fl2=1) summary(CaP.SK1) plot(CaP.SK1) CaP.SK2 <- SK.nest(x=dados[,1:3], y=dados[,4], model='resp~cultivares*calcio*fosforo',which='calcio:fosforo',fl2=2) summary(CaP.SK2) plot(CaP.SK2) ### Ainda podemos estudar as outras duplas # cultivares:calcio # Ca0 Ca1 # C1 870.500 566.00 # C2 1038.875 694.25 # C3 1149.750 633.50 # cultivares:fosforo # P0 P1 # C1 957.750 478.75 # C2 1237.625 495.50 # C3 1220.750 562.50 # Mais isso fica como exercicio # Triplas # C1 C2 C3 # Ca0 Ca1 Ca0 Ca1 Ca0 Ca1 # P0 1211 704.5 1491.75 983.5 1664.0 777.5 # P1 530 427.5 586.00 405.0 635.5 489.5 CaPCul.SK1 <- SK.nest(x=dados[,1:3], y=dados[,4], model='resp~calcio:fosforo:cultivares',which='calcio:fosforo:cultivares',fl2=1,fl3=1) summary(CaPCul.SK1) plot(CaPCul.SK1) ## Proceder de maneira analoga para as outras comparações detach(dados) ###---------------------------------------------------------------------------------------------------------------------#### ############################################### 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$aveia <- as.factor(ex07$aveia) ex07$trat <- as.factor(ex07$trat) ex07$bloco <- as.factor(ex07$bloco) summary(ex07) attach(ex07) ###Análise de variância ex07.av <- aov(resp ~ bloco + aveia*trat + Error(bloco:aveia)) 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(aveia, trat, resp, xlab="Variedades") interaction.plot(trat, aveia, 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 + aveia*trat + bloco*aveia) 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~trat) boxplot(ex07.avb$res~aveia) boxplot(ex07.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(ex07.avb$res, trat) bartlett.test(ex07.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 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,aveia,Dfea,QMEa,alpha=0.05) ##-----------------------------------------------------------------------## ###Fator secundario trat 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,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") ex07.avISP<-aov(resp ~ bloco +aveia/trat + aveia*trat -trat + Error(bloco:aveia)) summary(ex07.avISP) effects(ex07.avISP$'Within') summary(ex07.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(ex07$resp,list(Aveia=ex07$aveia,Tratamento=ex07$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(ex07,ex07$aveia=="1") A1 HSD.test(A1$resp,A1$trat,DFRb,QMRb) ## b dentro de A2 A2<-subset(ex07,ex07$aveia=="2") HSD.test(A2$resp,A2$trat,DFRb,QMRb) ## b dentro de A3 A3<-subset(ex07,ex07$aveia=="3") HSD.test(A3$resp,A3$trat,DFRb,QMRb) ## b dentro de A4 A4<-subset(ex07,ex07$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 ex07.avIPS<-aov(resp ~ bloco + trat/aveia + trat*aveia -aveia) summary(ex07.avIPS) effects(ex07.avIPS) summary(ex07.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(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) ) 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(ex07$resp,list(Aveia=ex07$aveia,Tratamento=ex07$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(ex07,ex07$trat=="1") B1 HSD.test(B1$resp,B1$aveia,Glx,QMex) ## aveia dentro de trat2 B2<-subset(ex07,ex07$trat=="2") HSD.test(B2$resp,B2$aveia,Glx,QMex) ## aveia dentro de trat3 B3<-subset(ex07,ex07$trat=="3") HSD.test(B3$resp,B3$aveia,Glx,QMex) ## aveia dentro de trat4 B4<-subset(ex07,ex07$trat=="4") HSD.test(B4$resp,B4$aveia,Glx,QMex) detach(ex07) ###---------------------------------------------------------------------------------------------------------------------#### ################################################ 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) abline(ajuste1) text(1.5,1,'y=0.828+0.40x²+0.36x³') ### limpando a memoria rm(list=ls()) ###-------------------------------------------------------------------------------### ############################## Regressão linear Multipla ############################ ### 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) arvore plot(arvore) ### Verificando a estrutura de correlação entre as variaveis cor(arvore) cor.test(X1,X2) require(agricolae) correlation(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) ###-------------------------------------------------------------------------------### ### 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(Y~X1+X2+X3) summary(ajuste.m) anova(ajuste.m) ### 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) 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 yy <- outer(xx1,xx2,function(xx1,xx2){predict(ajuste.x1x2,list(X1=xx1,X2=xx2))}) ### Função para máximo fmax<-function(x,y,z){ MC<-apply(z,2,max) j<-which.max(MC) i<-which.max(z[,j]) xm<-x[i] ym<-y[j] return(print(c(xm,ym,z[i,j]))) } ### Função para minimo fmin<-function(x,y,z){ MC<-apply(z,2,min) j<-which.min(MC) i<-which.min(z[,j]) xm<-x[i] ym<-y[j] return(print(c(xm,ym,z[i,j]))) } ###Encontrando os pontos Máximo e Minimo Zmax<-fmax(xx1,xx2,yy) Zmin<-fmin(xx1,xx2,yy) #contour(xx1,xx2,yy,nlevels=15) #abline(v=Zmax[1],h=Zmax[2],col='red') ### 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) ### Seleção de modelos library(MASS) ##?stepAIC modelo1 <- lm(Y~1) result1 <- stepAIC(modelo1,direction="forward",scope=list(upper = ~X1+X2+X3+I(X3^2), lower = ~1),trace=1) summary(result1) ####-----------------------------------------------------------------------------#### ############################## 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 <- nls(y~a*b^x,data=data,start=list(a=0.5,b=2),alg='plinear',trace=TRUE) summary(model) ## Algoritimo de port model <- nls(y~a*b^x,data=data,start=list(a=0.5,b=2),alg='port',trace=TRUE) summary(model) ### Outros argumentos podem ser modificados na "nls" ### ?nls ### ?nls.control control=nls.control(maxiter = 50,### Numero de intetraçõ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)) #parA <- c(0.5,0.3579661 , 0.6471207 , 1.9011656 , 2.340576 , 2.4143688 , 2.5178626 , 2.5248613 , 2.5249049 ) #parB <- c(2, 1.6940760 , 1.0308940 , 0.3765157 , 0.757061 , 0.6336439 , 0.5980226 , 0.5968229 , 0.5968138 ) #lines(parA,parB) 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 atribuição 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) ###-------------------------------------------------------------------------------### ############################################################################################################# ### 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=1)) { 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.....