# limpar todas as variáveis rm(list = ls(all.names = TRUE)) # Carregando os pacotes necessários #library(nlme) library(car) #library(ggplot2) library(xtable) #library(plotrix) #library(plyr) library(tidyverse) library(investr) library(ggpubr) # Arquivos com alguns preditores não lineares source("C:\\Users\\USER\\cnaber\\Disciplinas\\Pos\\1_semestre_2021\\Regressao MI 406\\programas\\Modelos nao lineares Pos Reg MI406 1S 2021.r") # Funções auxiliares source("C:\\Users\\USER\\cnaber\\Disciplinas\\Pos\\1_semestre_2021\\Regressao MI 406\\programas\\Funcoes Auxiliares Pos Reg MI406 1S 2021.r") source("C:\\Users\\USER\\cnaber\\Disciplinas\\Pos\\1_semestre_2021\\Regressao MI 406\\programas\\diag_nls.r") source("C:\\Users\\USER\\cnaber\\Disciplinas\\Pos\\1_semestre_2021\\Regressao MI 406\\programas\\envel_nls.r") # Onde salvar file.save<-"C:\\Users\\USER\\cnaber\\Disciplinas\\Pos\\1_semestre_2021\\Regressao MI 406\\aulas\\Modelos nao lineares\\" x=c(11.5,13.0,14.3,15.6,16.0,17.3,19.3,21.1,21.5,22.6,22.6,24.0,24.0,24.6,25.2,25.5,26.3,27.9,28.3,28.4,28.6,30.9,31.9, 34.5,40.1,40.1,43.0,44.1,46.5,47.3,48.7,52.9,56.9,59.9,60.2,60.3,60.5,62.1,62.8,66.5,67.0,67.1,67.9,68.8,75.4,100.5) # resposta y=c(3280,5046,1563,4707,977,2834,2266,2208,1040,700,1583,482,804,1093,1125,884,1300,852,580,1066,1114,386,745,736,750, 316,456,552,355,242,190,127,185,255,195,283,212,327,373,125,187,135,245,137,200,190) pdf(paste(file.save,sep="","DispEx13.pdf")) plot(x,y,xlab="trabalho por ciclo",ylab="número de ciclos até ocorrer a falha",cex=1.2,cex.axis=1.2,cex.lab=1.2,pch=19) dev.off() fadiga <- data.frame(cbind(x,y)) colnames(fadiga) <- c("trabalho","ciclos") # Comparação entre os dados e curvas "simuladas" pdf(paste(file.save,sep="","DispOCEEx13.pdf")) par(mfrow=c(2,2)) # Modelo 1 plot(x,y,xlab="trabalho por ciclo",ylab="número de ciclos até ocorrer a falha",cex=1.2,cex.axis=1.2,cex.lab=1.2,pch=19) lines(x,(exp1(x,0.01,400,30)),type="p",pch=17,col=2,cex=1.2) title("Modelo 1",cex=1.2) legend(40,3000,col=c(1,2),c("observado","predito"),pch=c(19,17),cex=1.2,bty="n") # M2 plot(x,y,xlab="trabalho por ciclo",ylab="número de ciclos até ocorrer a falha",cex=1.2,cex.axis=1.2,cex.lab=1.2,pch=19) lines(x,(frac(x,10,-40000,0.01)),type="p",pch=17,col=2,cex=1.2) title("Modelo 2",cex=1.2) legend(40,3000,col=c(1,2),c("observado","predito"),pch=c(19,17),cex=1.2,bty="n") # M3 plot(x,y,xlab="trabalho por ciclo",ylab="número de ciclos até ocorrer a falha",cex=1.2,cex.axis=1.2,cex.lab=1.2,pch=19) lines(x,(potencia(x,10000,-0.8)),type="p",pch=17,col=2,cex=1.2) title("Modelo 3",cex=1.2) legend(40,3000,col=c(1,2),c("observado","predito"),pch=c(19,17),cex=1.2,bty="n") legend(40,3000,col=c(1,2),c("observado","predito"),pch=c(19,17),cex=1.2,bty="n") dev.off() # M4 #plot(x,y,xlab="trabalho por ciclo",ylab="número de ciclos até ocorrer a falha",cex=1.2,cex.axis=1.2,cex.lab=1.2,pch=19) #lines(SSlogis(x,20000,10,-10),type="p") # Modelo 5 #plot(x,y,xlab="trabalho por ciclo",ylab="número de ciclos até ocorrer a falha",cex=1.2,cex.axis=1.2,cex.lab=1.2,pch=19) #lines(x,(exp2(x,10,100,10,100)),type="p",pch=17,col=2,cex=1.2) # Modelo M1 result1<-nls(ciclos~exp1(trabalho,alpha,beta,gamma),start=c(alpha=0.01,beta=5000,gamma=30),data=fadiga) pdf(paste(file.save,sep="","DiagmEx13M1.pdf")) diagnls(result1) dev.off() # pdf(paste(file.save,sep="","EnvelmEx13M1.pdf")) par(mfrow=c(1,1)) envelnls(result1) dev.off() # Modelo M2 result2<-nls(ciclos~frac(trabalho,alpha,beta,gamma),start=c(alpha=1,beta=-10000,gamma=1),data=fadiga) pdf(paste(file.save,sep="","DiagmEx13M2.pdf")) diagnls(result2) dev.off() # pdf(paste(file.save,sep="","EnvelmEx13M2.pdf")) par(mfrow=c(1,1)) envelnls(result2) dev.off() # Modelo M3 result3<-nls(ciclos~potencia(trabalho,alpha,beta),data=fadiga,start=c(alpha=10000,beta=-0.3)) pdf(paste(file.save,sep="","DiagmEx13M3.pdf")) diagnls(result3) dev.off() # pdf(paste(file.save,sep="","EnvelmEx13M3.pdf")) par(mfrow=c(1,1)) envelnls(result3) dev.off() # Valores preditos pred1 <- predict(result1,type="response",se.fit=TRUE) pred2 <- predict(result2,type="response") pred3 <- predict(result3,type="response") # pdf(paste(file.save,sep="","VOPEx13M3.pdf")) par(mfrow=c(1,1)) plot(x,y,xlab="trabalho por ciclo",ylab="número de ciclos até ocorrer a falha",cex=1.2,cex.axis=1.2,cex.lab=1.2,pch=19) lines(x,pred1,col=2,type="b",pch=19,cex=1.2) lines(x,pred2,col=3,type="b",pch=19,cex=1.2) lines(x,pred3,col=4,type="b",pch=19,cex=1.2) legend(60,4000,col=c(2,3,4),legend=c("modelo 1","modelo 2","modelo 3"),lty=c(1,1,1),bty="n",lwd=c(2,2,2)) dev.off() # Resultados summary(result1) summary(result2) summary(result3) #summary(result4) estat<-rbind(cbind(AIC(result1),BIC(result1)),cbind(AIC(result2),BIC(result2)),cbind(AIC(result3),BIC(result3))) xtable(estat) # Estatísticas de comparação de modelos result_estat<-data.frame(calc.estat.mod.comp.MRNNLH(result1,3), calc.estat.mod.comp.MRNNLH(result2,3), calc.estat.mod.comp.MRNNLH(result3,2)) # colnames(result_estat)<-c("modelo 1","modelo 2","modelo 3") xtable(result_estat) # Valores preditos com intervalos de confiancça interval1 <- as_tibble(predFit(result1, interval = "confidence", level= 0.95)) interva1aux <- data.frame(interval1,fadiga) p1 <- ggplot(fadiga) + geom_point(aes(x=trabalho, y=ciclos),size=2, colour="black") + xlab("trabalho por ciclo") + ylab("ciclos") p1aux <- p1+ geom_line(data=interva1aux, aes(x = trabalho, y = fit ))+ geom_ribbon(data=interva1aux, aes(x=trabalho, ymin=lwr, ymax=upr), alpha=0.5, inherit.aes=F, fill="blue")+ theme_classic()+ggtitle("Modelo 1") # interval2 <- as_tibble(predFit(result2, interval = "confidence", level= 0.95)) interva2aux <- data.frame(interval2,fadiga) p2 <- ggplot(fadiga) + geom_point(aes(x=trabalho, y=ciclos),size=2, colour="black") + xlab("trabalho por ciclo") + ylab("ciclos") p2aux <-p2 + geom_line(data=interva2aux, aes(x = trabalho, y = fit ))+ geom_ribbon(data=interva2aux, aes(x=trabalho, ymin=lwr, ymax=upr), alpha=0.5, inherit.aes=F, fill="blue")+ theme_classic()+ggtitle("Modelo 2") # interval3 <- as_tibble(predFit(result3, interval = "confidence", level= 0.95)) interva3aux <- data.frame(interval3,fadiga) p3 <- ggplot(fadiga) + geom_point(aes(x=trabalho, y=ciclos),size=2, colour="black") + xlab("trabalho por ciclo") + ylab("ciclos") p3aux <-p3+ geom_line(data=interva3aux, aes(x = trabalho, y = fit ))+ geom_ribbon(data=interva3aux, aes(x=trabalho, ymin=lwr, ymax=upr), alpha=0.5, inherit.aes=F, fill="blue")+ theme_classic()+ggtitle("Modelo 3") pdf(paste(file.save,sep="","VOPICEUSx13M3.pdf")) ggarrange(p1aux,p2aux,p3aux) dev.off() # Em um gráfico só pdf(paste(file.save,sep="","VOPICEx13M3.pdf")) p1 + geom_line(data=interva1aux, aes(x = trabalho, y = fit))+ geom_ribbon(data=interva1aux, aes(x=trabalho, ymin=lwr, ymax=upr,fill="M1"), alpha=0.5, inherit.aes=F) + geom_line(data=interva2aux, aes(x = trabalho, y = fit ))+ geom_ribbon(data=interva2aux, aes(x=trabalho, ymin=lwr, ymax=upr,fill="M2"), alpha=0.5, inherit.aes=F) + geom_line(data=interva3aux, aes(x = trabalho, y = fit ))+ geom_ribbon(data=interva3aux, aes(x=trabalho, ymin=lwr, ymax=upr,fill="M3"), alpha=0.5, inherit.aes=F)+ scale_color_manual(name = " ", breaks = c("M1", "M2", "M3"), values = c("M1" = "blue", "M2" = "green", "M3" = "red"))+ #scale_colour_manual(values = c("y1" = "darkred", "y2" = "red" )) theme_classic() dev.off() # Modelo 3 summary(result3) ajust3<-res.aju.mod.nls(result3,0.95,digitos=4)$result xtable(ajust3,digits=3)