#### Script para fazer a comparação entre o modelo geoestatisco e o GAM dados <- read.table("dados.csv",header=TRUE,sep=";",dec=",") ## Mapa com os pontos de amostragem with(dados,plot(X~Y)) ## Requerendo os pacotes necessários require(geoR) require(MASS) require(mgcv) ## Ajuste do modelo geoestatístico dados.geo = as.geodata(dados) ## Analise descritiva plot(dados.geo) #boxcox(dados.geo, lam=seq(0,6,l=150)) ## Tem que fazer transformação modelo <- likfit(dados.geo, lambda=1, ini=c(0.8, 20)) modelo ## Fazendo a Krigagem gx <- seq(2,150,l=100) gy <- seq(2,120,l=100) gr <- expand.grid(gx,gy) kc <- krige.conv(dados.geo, loc=gr, krige=krige.control(obj=modelo)) pred.geo <- data.frame(kc$predict) par(mfrow=c(1,3)) image(gx,gy,matrix(pred.geo[1:100^2,],ncol=100,nrow=100), asp=0, col=gray(seq(1, 0, len=10)), main="Krigagem com 128 pontos", xlab="Coordenada X", ylab="Coordenada Y", xlim=c(0,155),ylim=c(-10,125)) points(dados$X,dados$Y,pch=19,cex=0.5,col="white") legend.krige(y.leg=c(-8,-4), x.leg=c(15,140),val=c(0,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.25),vertical=FALSE,col=gray(seq(1, 0, len=10))) ## Fazendo via Thin Plate Splines ajuste.gam <- gam(Produtividade ~ s(X, Y),data = dados) (ajuste.gam) ndados <- data.frame(X=gr$Var1, Y= gr$Var2) pred <- predict(ajuste.gam,ndados,type="response",se.fit=TRUE) pred <- data.frame(pred) image(gx,gy,matrix(pred$fit,ncol=100,nrow=100), asp=0, col=gray(seq(1, 0, len=10)), main="Thin Plate Splines com 128 pontos", xlab="Coordenada X", ylab="Coordenada Y", xlim=c(0,155),ylim=c(-10,125)) points(dados$X,dados$Y,pch=19,cex=0.5,col="white") legend.krige(y.leg=c(-8,-4), x.leg=c(15,140),val=c(0,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.25),vertical=FALSE,col=gray(seq(1, 0, len=10))) ## Fazendo via Tensor Product Splines ajuste.gam.te <- gam(Produtividade ~ te(X, Y),data = dados) summary(ajuste.gam.te) ndados <- data.frame(X=gr$Var1, Y= gr$Var2) pred.te <- predict(ajuste.gam.te,ndados,type="response",se.fit=TRUE) pred.te <- data.frame(pred.te) image(gx,gy,matrix(pred.te$fit,ncol=100,nrow=100), asp=0, col=gray(seq(1, 0, len=10)), main="Tensor Product Splines com 128 pontos", xlab="Coordenada X", ylab="Coordenada Y", xlim=c(0,155),ylim=c(-10,125)) points(dados$X,dados$Y,pch=19,cex=0.5,col="white") legend.krige(y.leg=c(-8,-4), x.leg=c(15,140),val=c(0,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.25),vertical=FALSE,col=gray(seq(1, 0, len=10))) ########### Obtendo os erros de predição por Thin Plate e Tensor Product########## dados = data.frame(TPS= data.frame(data.frame(fitted(ajuste.gam)), TEP = data.frame(fitted(ajuste.gam.te)),REAL = dados$Produtividade)) names(dados) <- c("TPS","TEP","REAL") ## Erro quadrático médio mean((dados$TPS - dados$REAL)^2) mean((dados$TEP - dados$REAL)^2) ## Erro absoluto médio mean(abs(dados$TPS - dados$REAL)) mean(abs(dados$TEP - dados$REAL)) ## Correlação entre preditos e observados cor(dados$TPS,dados$REAL) cor(dados$TEP,dados$REAL) ## Nível de cobertura TPS TPS = predict(ajuste.gam,se.fit=TRUE) TPS.erro = data.frame(minimo = TPS$fit - qnorm(0.975)*TPS$se.fit, maximo = TPS$fit + qnorm(0.975)*TPS$se.fit) TPS.erro$REAL <- dados$REAL TPS.erro$Cob <- ifelse(TPS.erro$REAL < TPS.erro$minimo | TPS.erro$REAL > TPS.erro$maximo,0,1) mean(TPS.erro$Cob) ## Nível de cobertura TEP TEP = predict(ajuste.gam.te,se.fit=TRUE) TEP.erro = data.frame(minimo = TEP$fit - qnorm(0.975)*TEP$se.fit, maximo = TEP$fit + qnorm(0.975)*TEP$se.fit) TEP.erro$REAL <- dados$REAL TEP.erro$Cob <- ifelse(TEP.erro$REAL < TEP.erro$minimo | TEP.erro$REAL > TEP.erro$maximo,0,1) mean(TEP.erro$Cob) ########################################################################################################################################### ### Ajustando com menos amostras ###################################################################################################################################################################################################################################################### comparativo <- function(dados,dim.gridi){ ## Montando o gridi de predicao gx <- seq(2,150,l=dim.gridi) gy <- seq(2,120,l=dim.gridi) gr <- expand.grid(gx,gy) ## Ajustando os modelos ## GEO dados.geo = as.geodata(dados) modelo <- likfit(dados.geo, lambda=1, ini=c(0.8, 20)) kc <- krige.conv(dados.geo, loc=gr, krige=krige.control(obj=modelo)) pred.geo <- data.frame(kc$predict) ## TPS ajuste.gam <- gam(Produtividade ~ s(X, Y),data = dados) ##TEP ajuste.gam.te <- gam(Produtividade ~ te(X, Y),data = dados) ## Plotando as superficie par(mfrow=c(1,3)) kc <- krige.conv(dados.geo, loc=gr, krige=krige.control(obj=modelo)) pred.geo <- data.frame(kc$predict) image(gx,gy,matrix(pred.geo[1:100^2,],ncol=100,nrow=100), asp=0, col=gray(seq(1, 0, len=10)), main="Krigagem com 128 pontos", xlab="Coordenada X", ylab="Coordenada Y", xlim=c(0,155),ylim=c(-10,125)) points(dados$X,dados$Y,pch=19,cex=0.5,col="white") legend.krige(y.leg=c(-8,-4), x.leg=c(15,140),val=c(0,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.25),vertical=FALSE,col=gray(seq(1, 0, len=10))) ndados <- data.frame(X=gr$Var1, Y= gr$Var2) pred <- predict(ajuste.gam,ndados,type="response",se.fit=TRUE) pred <- data.frame(pred) image(gx,gy,matrix(pred$fit,ncol=100,nrow=100), asp=0, col=gray(seq(1, 0, len=10)), main="Thin Plate Splines com 128 pontos", xlab="Coordenada X", ylab="Coordenada Y", xlim=c(0,155),ylim=c(-10,125)) points(dados$X,dados$Y,pch=19,cex=0.5,col="white") legend.krige(y.leg=c(-8,-4), x.leg=c(15,140),val=c(0,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.25),vertical=FALSE,col=gray(seq(1, 0, len=10))) ndados <- data.frame(X=gr$Var1, Y= gr$Var2) pred.te <- predict(ajuste.gam.te,ndados,type="response",se.fit=TRUE) pred.te <- data.frame(pred.te) image(gx,gy,matrix(pred.te$fit,ncol=100,nrow=100), asp=0, col=gray(seq(1, 0, len=10)), main="Tensor Product Splines com 128 pontos", xlab="Coordenada X", ylab="Coordenada Y", xlim=c(0,155),ylim=c(-10,125)) points(dados$X,dados$Y,pch=19,cex=0.5,col="white") legend.krige(y.leg=c(-8,-4), x.leg=c(15,140),val=c(0,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.25),vertical=FALSE,col=gray(seq(1, 0, len=10))) ## ## Erro quadrático médio eq.tps = mean((fitted(ajuste.gam) - dados$Produtividade)^2) eq.tep = mean((fitted(ajuste.gam.te) - dados$Produtividade)^2) ## Erro absoluto médio ea.tps = mean(abs(fitted(ajuste.gam) - dados$Produtividade)) ea.tep = mean(abs(fitted(ajuste.gam.te) - dados$Produtividade)) ## Correlação entre preditos e observados co.tps = cor(fitted(ajuste.gam),dados$Produtividade) co.tep = cor(fitted(ajuste.gam.te),dados$Produtividade) ## Nível de cobertura TPS TPS = predict(ajuste.gam,se.fit=TRUE) TPS.erro = data.frame(minimo = TPS$fit - qnorm(0.975)*TPS$se.fit, maximo = TPS$fit + qnorm(0.975)*TPS$se.fit) TPS.erro$REAL <- dados$Produtividade TPS.erro$Cob <- ifelse(TPS.erro$REAL < TPS.erro$minimo | TPS.erro$REAL > TPS.erro$maximo,0,1) cob.tps = mean(TPS.erro$Cob) ## Nível de cobertura TEP TEP = predict(ajuste.gam.te,se.fit=TRUE) TEP.erro = data.frame(minimo = TEP$fit - qnorm(0.975)*TEP$se.fit, maximo = TEP$fit + qnorm(0.975)*TEP$se.fit) TEP.erro$REAL <- dados$Produtividade TEP.erro$Cob <- ifelse(TEP.erro$REAL < TEP.erro$minimo | TEP.erro$REAL > TEP.erro$maximo,0,1) cob.tep = mean(TEP.erro$Cob) resultado = data.frame(eq.tps,eq.tep,ea.tps,ea.tep,co.tps,co.tep,cob.tps,cob.tep) return(resultado) } comparativo(dados,100)