## Carregando pacotes require(INLA) require(mgcv) ## Lendo a base de dados dados <- read.table("dados.txt",header=TRUE) summary(dados) ## Ordenando a base teste = dados[order(dados$DATA),] ## Cortando por grupos reser <- split(teste,teste$UHE) reserva <- c(3,6,7,11,12,13,14,15) total <- c() for(i in 1:length(reserva)){ total[i] = table(reser[[reserva[i]]]$LOCAL)[1]} total <- as.numeric(total) for(i in 1:length(total)){ reser[[reserva[i]]]$tempo <- rep(1:total[i],each = 3)} ## Pegando apenas alguns reservatorios dados <- rbind(reser[[3]],reser[[6]],reser[[7]],reser[[11]],reser[[12]],reser[[13]],reser[[14]],reser[[15]]) dados$LOCAL <- factor(dados$LOCAL,levels=c("MONT","RESER","JUSA")) ## Transformando a resposta vari <- (dados$IQA/100)/(1-dados$IQA/100) dados$vari <- log(vari) dados$UHE <- as.character(dados$UHE) dados$UHE <- as.factor(dados$UHE) ## Análise exploratória source("descritiva.R") par(mfrow=c(2,3),mar=c(5,2.8,1.7,1),mgp=c(1.8,0.7,0),oma=c(0,0,0,0)) descritiva(dados,tit.grafico="IQA") ############################################################################################################################ # Modelo 1 INLA formu1 <- vari ~ 1 mod1 <- inla(formu1,family="gaussian",control.compute=list(dic=TRUE,cpo=TRUE,mlik=TRUE),data=dados) summary(mod1) plot(mod1) # Modelo 1 pelo GAM mod1.gam <- gam(formu1,data=dados) summary(mod1.gam) AIC(mod1.gam) data.frame(INLA = mod1$summary.fixed[,1], GAM = mod1.gam$coefficients) # Modelo 2 pelo INLA formu2 <- vari ~ UHE + LOCAL mod2 <- inla(formu2,family="gaussian",control.compute=list(dic=TRUE,cpo=TRUE,mlik=TRUE),data=dados) summary(mod2) plot(mod2) # Modelo 2 pelo gam mod2.gam <- gam(vari ~ UHE + LOCAL, data=dados) summary(mod2.gam) plot(mod2.gam) data.frame(INLA = mod2$summary.fixed[,1], GAM = mod2.gam$coefficients) AIC(mod2.gam) # Modelo 3 formu3 <- vari ~ f(UHE,model="iid",param=c(1,0.01)) + LOCAL + f(tempo,model="rw1",param=c(1,0.01)) mod3 <- inla(formu3,family="gaussian",control.compute=list(dic=TRUE,cpo=TRUE,mlik=TRUE),data=dados) summary(mod3) ## DIC np MV par(mfrow=c(1,3),mar=c(5,2.8,1.7,1),mgp=c(1.8,0.7,0),oma=c(0,0,0,0)) plot(mod3$marginals.fixed$"(Intercept)",type="l",main="Intercepto") plot(mod3$marginals.fixed$"LOCALRESER",type="l",main="Reservatorio") plot(mod3$marginals.fixed$"LOCALJUSA",type="l",main="Jusante") formu3 <- vari ~ f(UHE,model="iid",param=c(0.01,0.01)) + LOCAL + f(tempo,model="rw1",param=c(0.01,0.01)) mod3.priori <- inla(formu3,family="gaussian",control.compute=list(dic=TRUE,cpo=TRUE,mlik=TRUE),data=dados) plot(mod3.priori) summary(mod3.priori) par(mfrow=c(1,3),mar=c(5,2.8,1.7,1),mgp=c(1.8,0.7,0),oma=c(0,0,0,0)) plot(mod3$marginals.hyper$"Precision for the Gaussian observations",type="l",main="Precisao - Gaussiana") lines(mod3.priori$marginals.hyper$"Precision for the Gaussian observations",lty=2,col="red") legend("topright", c("a=1,b=0.01","a=0.01, b=0.01"), col=c("black","red"),lty=c(1,2),cex=0.7) plot(mod3.priori$marginals.hyper$"Precision for UHE",type="l",pch=2,col="red",main="Precisao - UHE") lines(mod3$marginals.hyper$"Precision for UHE") legend("topright", c("a=1,b=0.01","a=0.01, b=0.01"), col=c("black","red"),lty=c(1,2),cex=0.7) plot(mod3.priori$marginals.hyper$"Precision for tempo",type="l",lty=2,col="red",main="Precisao - Tempo") lines(mod3$marginals.hyper$"Precision for tempo") legend("topright", c("a=1,b=0.01","a=0.01, b=0.01"), col=c("black","red"),lty=c(1,2),cex=0.7) # Modelo 3 pelo gam mod3.gam <- gam(vari ~ UHE + LOCAL + s(tempo),data=dados) summary(mod3.gam) par(mfrow=c(2,2)) plot(mod3.gam) AIC(mod3.gam) ## Comparando os resultados data.frame(INLA = mod3$summary.fixed[,1], GAM = mod3.gam$coefficients[c(1,9,10)]) par(mfrow=c(1,3),mar=c(5,2.8,1.7,1),mgp=c(1.8,0.7,0),oma=c(0,0,0,0)) plot(mod3$marginals.fixed$"(Intercept)",type="l",main="Intercepto") plot(mod3$marginals.fixed$"LOCALRESER",type="l",main="Reservatorio") plot(mod3$marginals.fixed$"LOCALJUSA",type="l",main="Jusante") ## Sobrepondo os efeitos temporais ## Pegando o efeito temporal do GAM tempo.gam <- predict(mod3.gam,type="terms",se.fit=TRUE) tempo.gam <- data.frame(int.min = unique(tempo.gam$fit[,3] - qnorm(0.975)*tempo.gam$se.fit[,3]), media = unique(tempo.gam$fit[,3]), int.max = unique(tempo.gam$fit[,3] + qnorm(0.975)*tempo.gam$se.fit[,3])) ## Pegando o efeito temporal do INLA tempo.inla <- data.frame(int.min = mod3$summary.random$tempo$"0.025quant", media = mod3$summary.random$tempo$mean, int.max = mod3$summary.random$tempo$"0.975quant") ## Fazendo o gráfico plot(tempo.inla$media,type="n",xlab="",ylim=c(-0.6, 0.6),xaxt="n",ylab="Efeito temporal") dta <- sort(unique(reser[[3]]$DATA)) labels <- dta axis(1, at = 1:23, labels = labels, las = 2,cex.axis=0.8) matlines(tempo.inla, lty = c(1, 2, 1), lwd=c(1,2,1),col=1) matlines(tempo.gam, lty = c(1, 2, 1),lwd=c(1,2,1), col=2) legend("topleft", c("INLA","GAM"), col=c(1,2),lty=1,cex=0.7) ### Modelo de interação LOCAL:TEMPO mat.LOCAL <- diag(3) source("mat.temporal.R") mat.tempo <- mat.temporal(ncol=23,nrow=23) mat.inte <- kronecker(mat.LOCAL,mat.tempo) mat.zero <- matrix(0,ncol=69,nrow=69) mat.1 <- cbind(mat.inte,mat.zero,mat.zero,mat.zero,mat.zero,mat.zero,mat.zero,mat.zero) mat.2 <- cbind(mat.zero,mat.inte,mat.zero,mat.zero,mat.zero,mat.zero,mat.zero,mat.zero) mat.3 <- cbind(mat.zero,mat.zero,mat.inte,mat.zero,mat.zero,mat.zero,mat.zero,mat.zero) mat.4 <- cbind(mat.zero,mat.zero,mat.zero,mat.inte,mat.zero,mat.zero,mat.zero,mat.zero) mat.5 <- cbind(mat.zero,mat.zero,mat.zero,mat.zero,mat.inte,mat.zero,mat.zero,mat.zero) mat.6 <- cbind(mat.zero,mat.zero,mat.zero,mat.zero,mat.zero,mat.inte,mat.zero,mat.zero) mat.7 <- cbind(mat.zero,mat.zero,mat.zero,mat.zero,mat.zero,mat.zero,mat.inte,mat.zero) mat.8 <- cbind(mat.zero,mat.zero,mat.zero,mat.zero,mat.zero,mat.zero,mat.zero,mat.inte) mat.intera <- rbind(mat.1,mat.2,mat.3,mat.4,mat.5,mat.6,mat.7,mat.8) ## Matriz de restricoes source("restricoes.R") A <- restri(n.dados=69,n.espaco=3,n.tempo=23) A.zero <- matrix(0,nrow=26,ncol=69) A.1 <- cbind(A,A.zero,A.zero,A.zero,A.zero,A.zero,A.zero,A.zero) A.2 <- cbind(A.zero,A,A.zero,A.zero,A.zero,A.zero,A.zero,A.zero) A.3 <- cbind(A.zero,A.zero,A,A.zero,A.zero,A.zero,A.zero,A.zero) A.4 <- cbind(A.zero,A.zero,A.zero,A,A.zero,A.zero,A.zero,A.zero) A.5 <- cbind(A.zero,A.zero,A.zero,A.zero,A,A.zero,A.zero,A.zero) A.6 <- cbind(A.zero,A.zero,A.zero,A.zero,A.zero,A,A.zero,A.zero) A.7 <- cbind(A.zero,A.zero,A.zero,A.zero,A.zero,A.zero,A,A.zero) A.8 <- cbind(A.zero,A.zero,A.zero,A.zero,A.zero,A.zero,A.zero,A) mat.restricao <- rbind(A.1,A.2,A.3,A.4,A.5,A.6,A.7,A.8) head(mat.restricao) e <- rep(0,nrow(mat.restricao)) source("mat.INLA.R") mat.intera = mat.INLA(mat.intera) dados$int <- rep(1:552) forma6 = vari ~ f(UHE,model="iid")+ LOCAL + f(tempo,model="rw1") + f(int,model="generic0", Cmatrix=list(i=mat.intera$linha,j=mat.intera$coluna, Cij=mat.intera$corpo),extraconstr=list(A=mat.restricao,e=e),param=c(1,0.01),constr=TRUE) mod6 <- inla(forma6,family="gaussian",control.inla=list(strategy="GAUSSIAN"), control.compute=list(dic=TRUE,cpo=TRUE,mlik=TRUE),data=dados) summary(mod6) ## DIC NP MV plot(mod6) ### UM GAM quase equivalente seria CAT.1 <- as.numeric(dados$LOCAL== "MONT") CAT.2 <- as.numeric(dados$LOCAL == "RESER") CAT.3 <- as.numeric(dados$LOCAL == "JUSA") mod.gam.int <- gam(vari ~ UHE + LOCAL + s(tempo,by=CAT.1) + s(tempo,by=CAT.2) + s(tempo,by=CAT.3),data=dados) summary(mod.gam.int) plot(mod.gam.int) AIC(mod.gam.int) ##################################################################################################################################### Modelos com interação LOCAL e TEMPO mat.UHE <- diag(8) mat.LOCAL <- diag(3) source("mat.temporal.R") mat.tempo <- mat.temporal(ncol=23,nrow=23) mat1 <- kronecker(mat.UHE,mat.LOCAL) mat.inte <- kronecker(mat1,mat.tempo) source("mat.INLA.R") mat.intera = mat.INLA(mat.inte) dados$int <- 1:dim(dados)[1] ## Matriz de restrições source("restricoes.R") A <- restri(n.dados=552,n.espaco=24,n.tempo=23) e <- rep(0,47) ## Formula do modelo forma6 = vari ~ f(UHE,model="iid")+ LOCAL + f(tempo,model="rw1") + f(int,model="generic0", Cmatrix=list(i=mat.intera$linha,j=mat.intera$coluna, Cij=mat.intera$corpo),extraconstr=list(A=A,e=e),param=c(1,0.01),constr=TRUE) mod7 <- inla(forma6,family="gaussian",control.inla=list(strategy="GAUSSIAN"), control.compute=list(dic=TRUE,cpo=TRUE,mlik=TRUE),data=dados) summary(mod7) ## DIC NP MV plot(mod7) ### Tabela comparando os modelos via INLA pelo DIC cbind(mod1$dic,mod2$dic,mod3$dic,mod6$dic,mod7$dic) cbind(mod1$mlik,mod2$mlik,mod3$mlik,mod6$mlik,mod7$mlik) ## Histogram do PIT para os quatro modelos par(mfrow=c(3,2)) hist(mod1$pit) hist(mod2$pit) hist(mod3$pit) hist(mod6$pit) hist(mod7$pit) ### Tabela comparando os modelos via GAM pelo AIC data.frame(mod1=AIC(mod1.gam),mod2=AIC(mod2.gam),mod3=AIC(mod3.gam),mod4=AIC(mod.gam.int)) par(mfrow=c(2,2)) gam.check(mod1.gam) gam.check(mod2.gam) gam.check(mod3.gam) gam.check(mod.gam.int) ## Escolhe o modelo com efeitos principais formu3 <- vari ~ f(UHE,model="iid",param=c(1,0.01)) + LOCAL + f(tempo,model="rw1",param=c(1,0.01)) mod3.gauss <- inla(formu3,family="gaussian",control.compute=list(dic=TRUE,cpo=TRUE,mlik=TRUE), control.inla=list(strategy="gaussian"),data=dados) summary(mod3.gauss) ## DIC np MV mod3.lasim <- inla(formu3,family="gaussian",control.compute=list(dic=TRUE,cpo=TRUE,mlik=TRUE), control.inla=list(strategy="simplified.laplace"),data=dados) summary(mod3.lasim) ## DIC -769.02 NP 25.23 MV 357.30 mod3.laplace <- inla(formu3,family="gaussian",control.compute=list(dic=TRUE,cpo=TRUE,mlik=TRUE), control.inla=list(strategy="laplace"),data=dados) summary(mod3.laplace) ## DIC -769.03 NP 25.22 MV 357.30 plot(mod3.laplace) ## Comparando as posteriori pelo diferentes métodos ## Não da diferença entre os metodos de integração par(mfrow=c(3,2),mar=c(2.8,2.8,1.7,1),mgp=c(1.8,0.7,0),oma=c(0,0,0,0)) plot(mod3.gauss$marginals.fixed$"(Intercept)",type="l",main="Intercepto") lines(mod3.lasim$marginals.fixed$"(Intercept)",type="l",pch=2,col="red") lines(mod3.laplace$marginals.fixed$"(Intercept)",type="l",pch=2,col="blue") plot(mod3.gauss$marginals.fixed$"LOCALRESER",type="l",main="Montante - Reservatorio") lines(mod3.lasim$marginals.fixed$"LOCALRESER",type="l",pch=2,col="red") lines(mod3.laplace$marginals.fixed$"LOCALRESER",type="l",pch=2,col="blue") plot(mod3.gauss$marginals.fixed$"LOCALJUSA",type="l",main="Montante - Jusante") lines(mod3.lasim$marginals.fixed$"LOCALJUSA",type="l",pch=2,col="red") lines(mod3.laplace$marginals.fixed$"LOCALJUSA",type="l",pch=2,col="blue") plot(mod3.gauss$marginals.hyper$"Precision for the Gaussian observations",type="l",main="Precisao Gaussiana") lines(mod3.lasim$marginals.hyper$"Precision for the Gaussian observations",col="red") lines(mod3.laplace$marginals.hyper$"Precision for the Gaussian observations",col="blue") plot(mod3.gauss$marginals.hyper$"Precision for UHE",type="l",main="Precisao UHE") lines(mod3.lasim$marginals.hyper$"Precision for UHE",col="red") lines(mod3.laplace$marginals.hyper$"Precision for UHE",col="blue") plot(mod3.gauss$marginals.hyper$"Precision for tempo",type="l",main="Precisao Tempo") lines(mod3.lasim$marginals.hyper$"Precision for tempo",col="red") lines(mod3.laplace$marginals.hyper$"Precision for tempo",col="blue") ######### Comparando os efeitos temporais tempo.gam <- predict(mod3.gam,type="terms",se.fit=TRUE) tempo.gam <- data.frame(int.min = unique(tempo.gam$fit[,3] - qnorm(0.975)*tempo.gam$se.fit[,3]), media = unique(tempo.gam$fit[,3]), int.max = unique(tempo.gam$fit[,3] + qnorm(0.975)*tempo.gam$se.fit[,3])) ## Pegando o efeito temporal do INLA tempo.inla <- data.frame(int.min = mod3.laplace$summary.random$tempo$"0.025quant", media = mod3.laplace$summary.random$tempo$mean, int.max = mod3.laplace$summary.random$tempo$"0.975quant") ## Fazendo o gráfico plot(tempo.inla$media,type="n",xlab="",ylim=c(-0.17, 0.1),xaxt="n",ylab="Efeito temporal") dta <- sort(unique(reser[[3]]$DATA)) labels <- dta axis(1, at = 1:23, labels = labels, las = 2,cex.axis=0.8) matlines(tempo.inla, lty = c(1, 2, 1), lwd=c(1,2,1),col=1) matlines(tempo.gam, lty = c(1, 2, 1),lwd=c(1,2,1), col=2) legend("topleft", c("INLA","GAM"), col=c(1,2),lty=1,cex=0.7) ## Estimativas e erros padrões saida.gam <- summary(mod3.gam) saida = data.frame(INLA = mod3.laplace$summary.fixed[,1:2], GAM= cbind(saida.gam[1]$p.coeff[c(1,9,10)],saida.gam[2]$se[c(1,9,10)])) saida ### FIM DO CODIGO #####