# limpar todas as variáveis rm(list = ls(all.names = TRUE)) # Pacotes necessários library(xtable) library(lattice) library(plyr) library(plotrix) set.seed(4143) # file.save file.save <- "C:\\Users\\USER\\cnaber\\Disciplinas\\Pos\\1_semestre_2021\\Regressao MI 406\\aulas\\Transformacao e outros metodos de estimacao\\" # # Arquivos necessários source("C:\\Users\\USER\\cnaber\\Disciplinas\\Pos\\1_semestre_2021\\Regressao MI 406\\programas\\diag2_norm.r") source("C:\\Users\\USER\\cnaber\\Disciplinas\\Pos\\1_semestre_2021\\Regressao MI 406\\programas\\envel_norm.r") source("C:\\Users\\USER\\cnaber\\Disciplinas\\Pos\\1_semestre_2021\\Regressao MI 406\\programas\\Funcoes Auxiliares Pos Reg MI406 1S 2021.r") # Exemplo 1 # Lendo os dados braga<-scan("C:\\Users\\USER\\cnaber\\Disciplinas\\Pos\\1_semestre_2021\\Regressao MI 406\\dados\\Braga1998.txt",what=list(Etiologia="",CARGA=0,VO2=0)) # Construindo as variáveis etio <-braga$Etiologia # etiologia carga <-braga$CARGA # carga vo2<-braga$VO2 # consumo # gerando uma variável qualitativa para ser usada na função lm etiofac<-factor(etio,levels=c("C","CH","ID","IS")) # Modelo com a variável explicativa não centrada fit.model<-result<-lm(vo2~-1+etiofac+carga:etiofac) res.fit.model<-res.aju.mod.lm(fit.model,0.95,3)$result xtable(res.fit.model,digits=4) # Análise de variância anores<-anova(result) xtable(anores) # Diagnóstico pdf(paste(file.save,sep="","Exemplo1Diag.pdf")) diag2norm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo1Envel.pdf")) par(mfrow=c(1,1)) envelnorm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo1AI.pdf")) par(mfrow=c(1,2)) anainflu_norm(fit.model,0.10) dev.off() # Modelo com a variável explicativa não centrada fit.model<-result<-lm(log(vo2)~-1+etiofac+carga:etiofac) res.fit.model<-res.aju.mod.lm(fit.model,0.95,3)$result xtable(res.fit.model,digits=4) # Análise de variância anores<-anova(result) xtable(anores) # Diagnóstico pdf(paste(file.save,sep="","Exemplo1TLNDiag.pdf")) diag2norm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo1TLNEnvel.pdf")) par(mfrow=c(1,1)) envelnorm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo1TLNAI.pdf")) par(mfrow=c(1,2)) anainflu_norm(fit.model,0.10) dev.off() # Exemplo 2 (fósforo) # definindo as variáveis de interesse # produção de milho prodmi <- c(2.38,6.77,3.50,5.94,6.15,8.78,8.99,9.10,9.07,8.73,6.92,8.48,9.55,8.95,10.24,8.66,9.14,10.17,9.75,9.50) # fósforo fosforoq <- c(matrix(0,4,1),matrix(25,4,1),matrix(50,4,1),matrix(75,4,1),matrix(100,4,1)) # Gráfico de dispersão plot(fosforoq,prodmi,cex=1.2,pch=19,cex.lab=1.4,cex.axis=1.4,cex.main=1.4,main="Produção de milho (kg/parcela) em função da quantidade de fósforo (kg/ha)",xlab="quantidade de fósforo(kg/ha)",ylab="produção de milho (kg/parcela)") # Medidas descritivas por nível de fósforo # Variável fósforo como um fator fosfac <- c(matrix("0 kg/ha",4,1),matrix("25 kg/ha",4,1),matrix("50 kg/ha",4,1),matrix("75 kg/ha",4,1),matrix("100 kg/ha",4,1)) datafosf<-data.frame(prodmi,fosfac) cfosf <-ddply(datafosf,.(fosfac),summarise,media=mean(prodmi),dp=sqrt(var(prodmi)),cv=100*((sqrt(var(prodmi))/mean(prodmi))),vari=var(prodmi),minimo=min(prodmi),maximo=max(prodmi)) cfosf <- rbind(cfosf[1,],cfosf[3:5,],cfosf[2,]) xtable(cfosf,digits=3) # Modelo linear de regressão fit.model<-lm(prodmi~fosforoq)# ajuste do modelo res.fit.model<-res.aju.mod.lm(fit.model,0.95,3)$result xtable(res.fit.model,digits=4) # Diagnóstico pdf(paste(file.save,sep="","Exemplo3MLDiag.pdf")) diag2norm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo3MLEnvel.pdf")) par(mfrow=c(1,1)) envelnorm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo3MLAI.pdf")) par(mfrow=c(1,2)) anainflu_norm(fit.model,0.10) dev.off() # Modelo linear de regressão fit.model<-lm(log(prodmi)~fosforoq)# ajuste do modelo res.fit.model<-res.aju.mod.lm(fit.model,0.95,3)$result xtable(res.fit.model,digits=4) # Diagnóstico pdf(paste(file.save,sep="","Exemplo3TLNMLDiag.pdf")) diag2norm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo3TLNMLEnvel.pdf")) par(mfrow=c(1,1)) envelnorm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo3TLNMLAI.pdf")) par(mfrow=c(1,2)) anainflu_norm(fit.model,0.10) dev.off() # Modelo quadrático de regressão fosforoq2 <- fosforoq^2 # ajuste do modelo # fit.model<-lm(prodmi~fosforoq+fosforoq2)# ajuste do modelo res.fit.model<-res.aju.mod.lm(fit.model,0.95,3)$result xtable(res.fit.model,digits=4) m.Xq<-model.matrix(fit.model) # matriz de planejamento do modelo quadrático sigma2 <- ((summary(fit.model))$sigma)^2 # Diagnóstico pdf(paste(file.save,sep="","Exemplo3MQDiag.pdf")) diag2norm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo3MQEnvel.pdf")) par(mfrow=c(1,1)) envelnorm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo3MQAI.pdf")) par(mfrow=c(1,2)) anainflu_norm(fit.model,0.10) dev.off() fit.model<-lm(log(prodmi)~fosforoq+fosforoq2)# ajuste do modelo res.fit.model<-res.aju.mod.lm(fit.model,0.95,3)$result xtable(res.fit.model,digits=4) m.Xq<-model.matrix(fit.model) # matriz de planejamento do modelo quadrático sigma2 <- ((summary(fit.model))$sigma)^2 # Diagnóstico pdf(paste(file.save,sep="","Exemplo3TLNMQDiag.pdf")) diag2norm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo3TLNMQEnvel.pdf")) par(mfrow=c(1,1)) envelnorm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo3TLNMQAI.pdf")) par(mfrow=c(1,2)) anainflu_norm(fit.model,0.10) dev.off() # Exemplo 5 (absorbancia) solvente <- c(matrix("E50",5,1),matrix("EAW",5,1),matrix("MAW",5,1),matrix("E70",5,1),matrix("M1M",5,1)) mabsor <- c(0.5553,0.5623,0.5585,0.5096,0.5110,0.5436,0.5660,0.5860,0.5731,0.5656,0.4748,0.4321,0.4309,0.5010,0.4094,0.6286,0.6143,0.5826,0.6079,0.6060,0.1651,0.1840,0.2144,0.2249,0.1954) solvfac <- C(as.factor(solvente)) # é possível escolher outro solvente de referência, por exemplo factor(solvente,levels=c("EAW","E50","E70","MAW","M1M")) fit.model<- lm(mabsor~solvfac) res.fit.model<-res.aju.mod.lm(fit.model,0.95,3)$result xtable(res.fit.model,digits=4) m.X<-model.matrix(fit.model) anovatab<-anova(fit.model) # Diagnóstico pdf(paste(file.save,sep="","Exemplo5Diag.pdf")) diag2norm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo5Envel.pdf")) par(mfrow=c(1,1)) envelnorm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo5EnvelAI.pdf")) par(mfrow=c(1,2)) anainflu_norm(fit.model,0.10) dev.off() fit.model<- lm(log(mabsor/(1-mabsor))~solvfac) res.fit.model<-res.aju.mod.lm(fit.model,0.95,3)$result xtable(res.fit.model,digits=4) # Diagnóstico pdf(paste(file.save,sep="","Exemplo5TLOGDiag.pdf")) diag2norm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo5TLOGEnvel.pdf")) par(mfrow=c(1,1)) envelnorm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo5TLOGAI.pdf")) par(mfrow=c(1,2)) anainflu_norm(fit.model,0.10) dev.off() # Simulação para comparar a estimação por MV e por MQ no modelo lognormal # Estimação por MV para beta0,beta1 para o modelo log-normal R <- 10000 #(número de réplicas) vlnorm <- function(parm,y,x) { beta0<-parm[1] beta1<-parm[2] logl<-sum(dlnorm(y,meanlog=beta0*x^(beta1),sdlog=1,log=TRUE)) return(-logl) } beta0<-0.5 beta1<-1.05 n<-100 x <- runif(n,1,10) y<-rlnorm(n,meanlog=beta0*x^(beta1),sdlog=1) # amq <- lm(log(y)~log(x)) amv <- optim(as.numeric(c(exp(coef(amq)[1]),coef(amq)[2])),vlnorm,y=y,x=x,hessian=T) mestMQ <- as.numeric(c(exp(coef(amq)[1]),coef(amq)[2])) mestMV <- c(amv$par) mvarMQ <- as.numeric(diag(vcov(amq))) mvarMV <- diag(solve(amv$hessian)) for (i in 2:R) { y<-rlnorm(n,meanlog=beta0*x^(beta1),sdlog=1) # amq <- lm(log(y)~log(x)) amv <- optim(as.numeric(c(exp(coef(amq)[1]),coef(amq)[2])),vlnorm,y=y,x=x,hessian=T) mestMQ <- rbind(mestMQ,as.numeric(c(exp(coef(amq)[1]),coef(amq)[2]))) mestMV <- rbind(mestMV,c(amv$par)) mvarMQ <- rbind(mvarMQ, as.numeric(diag(vcov(amq)))) mvarMV <- rbind(mvarMV, diag(solve(amv$hessian))) } mestMQ <-matrix(mestMQ,R,2) mestMV <-matrix(mestMV,R,2) # Histogramas pdf(paste(file.save,sep="","ResultSim.pdf")) par(mfrow=c(2,2)) a1<-hist(mestMQ[,1],plot=FALSE) a1$density <- a1$counts/sum(a1$counts)*100 plot(a1,freq=F,main=expression(paste("Estimador de MQ ",beta,0)),ylab="percentual",cex=1.4,cex.axis=1.4,cex.lab=1.4,xlab="estimativa") abline(v=beta0,lwd=3,col=2) a1<-hist(mestMQ[,2],plot=FALSE) a1$density <- a1$counts/sum(a1$counts)*100 plot(a1,freq=F,main=expression(paste("Estimador de MQ ",beta,1)),ylab="percentual",cex=1.4,cex.axis=1.4,cex.lab=1.4,xlab="estimativa",xlim=c(min(c(beta1,mestMQ[,2])),max(c(beta1,mestMQ[,2])))) abline(v=beta1,lwd=3,col=2) a1<-hist(mestMV[,1],plot=FALSE) a1$density <- a1$counts/sum(a1$counts)*100 plot(a1,freq=F,main=expression(paste("Estimador de MV ",beta,0)),ylab="percentual",cex=1.4,cex.axis=1.4,cex.lab=1.4,xlab="estimativa") abline(v=beta0,lwd=3,col=2) a1<-hist(mestMV[,2],plot=FALSE) a1$density <- a1$counts/sum(a1$counts)*100 plot(a1,freq=F,main=expression(paste("Estimador de MV ",beta,1)),ylab="percentual",cex=1.4,cex.axis=1.4,cex.lab=1.4,xlab="estimativa") abline(v=beta1,lwd=3,col=2) dev.off() # mMQ <- apply(mestMQ,2,mean) mMV <- apply(mestMV,2,mean) vMQ <- apply(mestMQ,2,var) vMV <- apply(mestMV,2,var) bMQ <- mMQ - c(beta0,beta1) bMV <- mMV - c(beta0,beta1) reMQ <- sqrt(vMQ+bMQ^2) reMV <- sqrt(vMV+bMV^2) # mmvarMQ <- apply(mvarMQ,2,mean) mmvarMV <- apply(mvarMV,2,mean) # #round(mmvarMQ,digits=3) #round(mmvarMV,digits=3) xtable(cbind(mmvarMQ,mmvarMV),digits=3) estat <- rbind(cbind(mMQ,bMQ,vMQ,reMQ),cbind(mMV,bMV,vMV,reMV)) round(estat,3) xtable(estat,digits=3) #amv <- optim(c(0.5,1.2),vlnorm,y=y,x=x,method="BFGS",hessian=T) # Lendo os dados sinand<-read.table("C:\\Users\\USER\\cnaber\\Disciplinas\\Pos\\1_semestre_2021\\Regressao MI 406\\dados\\Singer&Andrade1997.txt") # Construindo as variáveis de interesse hae <- sinand[,3] # IPB pré-escovação: escova Hugger hde <- sinand[,4] # IPB pós-escovação: escova Hugger cae <- sinand[,5] # IPB pré-escovação: escova tradicional cde <- sinand[,6] # IPB pós-escovação: escova tradicional # Dispersão par(mfrow=c(1,2)) plot(hae,hde,pch=19,cex=1.4,cex.lab=1.4,cex.axis=1.4,cex.main=1.4,xlab="IPB antes da escovação",ylab="IPB depois da escovação",main="Escova - Hugger",ylim=c(0,2),xlim=c(0,4)) plot(cae,cde,pch=19,cex=1.4,cex.lab=1.4,cex.axis=1.4,cex.main=1.4,xlab="IPB antes da escovação",ylab="IPB depois da escovação",main="Escova - Covencional",xlim=c(0,4)) # Montando a base de dados ipbpre <- c(cae,hae) # IPB pré-escovação ipbpos <- c(cde,hde) # IPB pós-escovação tesc <- c(rep("Convencional",26),rep("Hugger",26))# tipo de escova tescf <- factor(tesc,levels=c("Convencional","Hugger")) ipbposaux <- ifelse(ipbpos==0,0.001,ipbpos) # Ajustando o modelo (sem transformação) fit.model<-result <- lm(ipbposaux~-1+ tescf+ipbpre:tescf) res.fit.model<-res.aju.mod.lm(fit.model,0.95,3)$result xtable(res.fit.model,digits=4) m.X<-(model.matrix(result)) # Diagnóstico pdf(paste(file.save,sep="","Exemplo2Diag.pdf")) diag2norm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo2Envel.pdf")) par(mfrow=c(1,1)) envelnorm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo2AI.pdf")) par(mfrow=c(1,2)) anainflu_norm(fit.model,0.10) dev.off() # Ajustando o modelo (com transformação) fit.model<-result <- lm(log(ipbposaux)~-1+ tescf+log(ipbpre):tescf) res.fit.model<-res.aju.mod.lm(fit.model,0.95,3)$result xtable(res.fit.model,digits=4) m.X<-(model.matrix(result)) # Diagnóstico pdf(paste(file.save,sep="","Exemplo2TLDiag.pdf")) diag2norm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo2TLEnvel.pdf")) par(mfrow=c(1,1)) envelnorm(fit.model) dev.off() pdf(paste(file.save,sep="","Exemplo2TLAI.pdf")) par(mfrow=c(1,2)) anainflu_norm(fit.model,0.10) dev.off()