################################################################################ # Script utilizado para análise dos dados de solo da # # fazenda Tupã, Echaporan/SP # # # # Fonseca B.H.F. & Ribeiro Jr. P.J. # # 2008 # ################################################################################ ## Requerindo pacotes e carregando dados require(geoR) require(MASS) load("ph.RData") load("sat.RData") ## Análise Descritiva do pH summary(ph) sd(ph$data) boxplot(ph$data) points(ph) lines(ph$reg1) plot(ph) boxcox(ph,lambda=seq(-5,5)) ## Análise Descritiva da saturação por bases summary(sat) sd(sat$data) boxplot(sat$data) points(sat) lines(sat$reg1) plot(sat) boxcox(sat,lambda=seq(-2,6)) ## Modelagem Univariada ## Estimação por Máxima Verossimilhança do pH emv_ph1_1<-likfit(ph,ini.cov.pars=c(0.41,560),nugget=0.08, fix.kappa=TRUE,kappa=0.5) emv_ph1_2<-likfit(ph,ini.cov.pars=c(0.42,510.58),nugget=0.08, fix.kappa=TRUE,kappa=1) emv_ph1_3<-likfit(ph,ini.cov.pars=c(0.43,376.22),nugget=0.08, fix.kappa=TRUE,kappa=1.5) emv_ph1_4<-likfit(ph,ini.cov.pars=c(0.48,349.34),nugget=0.1, fix.kappa=TRUE,kappa=2) emv_ph1_5<-likfit(ph,ini.cov.pars=c(0.47,268.73),nugget=0.08, fix.kappa=TRUE,kappa=2.5) emv_ph2_1<-likfit(ph,trend=~AREA,ini.cov.pars=c(0.17,350),nugget=0.08, fix.kappa=TRUE,kappa=0.5) emv_ph2_2<-likfit(ph,trend=~AREA,ini.cov.pars=c(0.17,450),nugget=0.08, fix.kappa=TRUE,kappa=1) emv_ph2_3<-likfit(ph,trend=~AREA,ini.cov.pars=c(0.18,350),nugget=0.09, fix.kappa=TRUE,kappa=1.5) emv_ph2_4<-likfit(ph,trend=~AREA,ini.cov.pars=c(0.19,349.34),nugget=0.07, fix.kappa=TRUE,kappa=2) emv_ph2_5<-likfit(ph,trend=~AREA,ini.cov.pars=c(0.21,429.96),nugget=0.1, fix.kappa=TRUE,kappa=2.5) emv_ph3_1<-likfit(ph,trend=~coords[,1],ini.cov.pars=c(0.17,350),nugget=0.08, fix.kappa=TRUE,kappa=0.5) emv_ph3_2<-likfit(ph,trend=~coords[,1],ini.cov.pars=c(0.17,450),nugget=0.08, fix.kappa=TRUE,kappa=1) emv_ph3_3<-likfit(ph,trend=~coords[,1],ini.cov.pars=c(0.18,350),nugget=0.09, fix.kappa=TRUE,kappa=1.5) emv_ph3_4<-likfit(ph,trend=~coords[,1],ini.cov.pars=c(0.19,349.34),nugget=0.07, fix.kappa=TRUE,kappa=2) emv_ph3_5<-likfit(ph,trend=~coords[,1],ini.cov.pars=c(0.21,429.96),nugget=0.1, fix.kappa=TRUE,kappa=2.5) emv_ph4_1<-likfit(ph,trend=~coords[,1]+AREA,ini.cov.pars=c(0.17,350),nugget=0.08, fix.kappa=TRUE,kappa=0.5) emv_ph4_2<-likfit(ph,trend=~coords[,1]+AREA,ini.cov.pars=c(0.17,250),nugget=0.08, fix.kappa=TRUE,kappa=1) emv_ph4_3<-likfit(ph,trend=~coords[,1]+AREA,ini.cov.pars=c(0.18,150),nugget=0.09, fix.kappa=TRUE,kappa=1.5) emv_ph4_4<-likfit(ph,trend=~coords[,1]+AREA,ini.cov.pars=c(0.19,100),nugget=0.07, fix.kappa=TRUE,kappa=2) emv_ph4_5<-likfit(ph,trend=~coords[,1]+AREA,ini.cov.pars=c(0.21,100),nugget=0.1, fix.kappa=TRUE,kappa=2.5) ## Estimação por Máxima Verossimilhança do bases emv_sat1_1<-likfit(sat,ini.cov.pars=c(326.29,700),nugget=35, fix.kappa=TRUE,kappa=0.5) emv_sat1_2<-likfit(sat,ini.cov.pars=c(326.29,625.58),nugget=42.1, fix.kappa=TRUE,kappa=1) emv_sat1_3<-likfit(sat,ini.cov.pars=c(336.81,516.78),nugget=52.62, fix.kappa=TRUE,kappa=1.5) emv_sat1_4<-likfit(sat,ini.cov.pars=c(368.39,489.58),nugget=63.15, fix.kappa=TRUE,kappa=2) emv_sat1_5<-likfit(sat,ini.cov.pars=c(350,400),nugget=63.15, fix.kappa=TRUE,kappa=2.5) emv_sat2_1<-likfit(sat,trend=~AREA,ini.cov.pars=c(65,300),nugget=45, fix.kappa=TRUE,kappa=0.5) emv_sat2_2<-likfit(sat,trend=~AREA,ini.cov.pars=c(68.12,308.79),nugget=49.96, fix.kappa=TRUE,kappa=1) emv_sat2_3<-likfit(sat,trend=~AREA,ini.cov.pars=c(72.66,353.59),nugget=54.50, fix.kappa=TRUE,kappa=1.5) emv_sat2_4<-likfit(sat,trend=~AREA,ini.cov.pars=c(63.58,190.39),nugget=45.41, fix.kappa=TRUE,kappa=2) emv_sat2_5<-likfit(sat,trend=~AREA,ini.cov.pars=c(65,190.39),nugget=45.41, fix.kappa=TRUE,kappa=2.5) emv_sat3_1<-likfit(sat,trend=~coords[,1],ini.cov.pars=c(65,300),nugget=45, fix.kappa=TRUE,kappa=0.5) emv_sat3_2<-likfit(sat,trend=~coords[,1],ini.cov.pars=c(68.12,308.79),nugget=49.96, fix.kappa=TRUE,kappa=1) emv_sat3_3<-likfit(sat,trend=~coords[,1],ini.cov.pars=c(72.66,353.59),nugget=54.50, fix.kappa=TRUE,kappa=1.5) emv_sat3_4<-likfit(sat,trend=~coords[,1],ini.cov.pars=c(63.58,190.39),nugget=45.41, fix.kappa=TRUE,kappa=2) emv_sat3_5<-likfit(sat,trend=~coords[,1],ini.cov.pars=c(65,190.39),nugget=45.41, fix.kappa=TRUE,kappa=2.5) emv_sat4_1<-likfit(sat,trend=~coords[,1]+AREA,ini.cov.pars=c(50,250),nugget=45, fix.kappa=TRUE,kappa=0.5) emv_sat4_2<-likfit(sat,trend=~coords[,1]+AREA,ini.cov.pars=c(58,280),nugget=49.96, fix.kappa=TRUE,kappa=1) emv_sat4_3<-likfit(sat,trend=~coords[,1]+AREA,ini.cov.pars=c(60,250),nugget=54.50, fix.kappa=TRUE,kappa=1.5) emv_sat4_4<-likfit(sat,trend=~coords[,1]+AREA,ini.cov.pars=c(55,200),nugget=45.41, fix.kappa=TRUE,kappa=2) emv_sat4_5<-likfit(sat,trend=~coords[,1]+AREA,ini.cov.pars=c(56,200),nugget=45.41, fix.kappa=TRUE,kappa=2.5) ## Seleção de modelo do pH c(emv_ph1_1$AIC,emv_ph1_1$nospatial$AIC.ns) c(emv_ph1_2$AIC,emv_ph1_2$nospatial$AIC.ns) c(emv_ph1_3$AIC,emv_ph1_3$nospatial$AIC.ns) c(emv_ph1_4$AIC,emv_ph1_4$nospatial$AIC.ns) c(emv_ph1_5$AIC,emv_ph1_5$nospatial$AIC.ns) c(emv_ph2_1$AIC,emv_ph2_1$nospatial$AIC.ns) c(emv_ph2_2$AIC,emv_ph2_2$nospatial$AIC.ns) c(emv_ph2_3$AIC,emv_ph2_3$nospatial$AIC.ns) c(emv_ph2_4$AIC,emv_ph2_4$nospatial$AIC.ns) c(emv_ph2_5$AIC,emv_ph2_5$nospatial$AIC.ns) c(emv_ph3_1$AIC,emv_ph3_1$nospatial$AIC.ns) c(emv_ph3_2$AIC,emv_ph3_2$nospatial$AIC.ns) c(emv_ph3_3$AIC,emv_ph3_3$nospatial$AIC.ns) c(emv_ph3_4$AIC,emv_ph3_4$nospatial$AIC.ns) c(emv_ph3_5$AIC,emv_ph3_5$nospatial$AIC.ns) c(emv_ph4_1$AIC,emv_ph4_1$nospatial$AIC.ns) c(emv_ph4_2$AIC,emv_ph4_2$nospatial$AIC.ns) c(emv_ph4_3$AIC,emv_ph4_3$nospatial$AIC.ns) c(emv_ph4_4$AIC,emv_ph4_4$nospatial$AIC.ns) c(emv_ph4_5$AIC,emv_ph4_5$nospatial$AIC.ns) ## Seleção de modelo da Saturação c(emv_sat1_1$AIC,emv_sat1_1$nospatial$AIC.ns) c(emv_sat1_2$AIC,emv_sat1_2$nospatial$AIC.ns) c(emv_sat1_3$AIC,emv_sat1_3$nospatial$AIC.ns) c(emv_sat1_4$AIC,emv_sat1_4$nospatial$AIC.ns) c(emv_sat1_5$AIC,emv_sat1_5$nospatial$AIC.ns) c(emv_sat2_1$AIC,emv_sat2_1$nospatial$AIC.ns) c(emv_sat2_2$AIC,emv_sat2_2$nospatial$AIC.ns) c(emv_sat2_3$AIC,emv_sat2_3$nospatial$AIC.ns) c(emv_sat2_4$AIC,emv_sat2_4$nospatial$AIC.ns) c(emv_sat2_5$AIC,emv_sat2_5$nospatial$AIC.ns) c(emv_sat3_1$AIC,emv_sat3_1$nospatial$AIC.ns) c(emv_sat3_2$AIC,emv_sat3_2$nospatial$AIC.ns) c(emv_sat3_3$AIC,emv_sat3_3$nospatial$AIC.ns) c(emv_sat3_4$AIC,emv_sat3_4$nospatial$AIC.ns) c(emv_sat3_5$AIC,emv_sat3_5$nospatial$AIC.ns) c(emv_sat4_1$AIC,emv_sat4_1$nospatial$AIC.ns) c(emv_sat4_2$AIC,emv_sat4_2$nospatial$AIC.ns) c(emv_sat4_3$AIC,emv_sat4_3$nospatial$AIC.ns) c(emv_sat4_4$AIC,emv_sat4_4$nospatial$AIC.ns) c(emv_sat4_5$AIC,emv_sat4_5$nospatial$AIC.ns) ######### Modelos finais ############## #emv_sat3_5 & emv_ph3_5 ## Krigagem # Função para predição pred <- function(geo1,gr_pred,X_pred,parametros,kap,nug){ coords <- geo1$coords auxiliar <- as.numeric(geo1$covariate$AREA) auxiliar[auxiliar==1]<-0 auxiliar[auxiliar==2]<-1 n1 <- nrow(coords) npr <- nrow(gr_pred) nt <- npr + n1 k1 <- c(1:n1) k3 <- c(1:npr) k4 <- c((npr+1):nt) dados <- geo1$data if(ncol(X_pred) == 2){ if(X_pred[1,2] < 2){ X_dados <- matrix(c(rep(1,n1),auxiliar),n1,2) } if(X_pred[1,2] > 2){ X_dados <- matrix(c(rep(1,n1),coords[,1]),n1,2) } } if(ncol(X_pred)==1){ X_dados <- matrix(c(rep(1,n1)),n1,1) } X <- rbind(X_pred,X_dados) gr_dados <- unname(coords) gr <- rbind(gr_pred,gr_dados) dists <- matrix(0,nt,n1) for(i in 1:nt){ dists[i,] <- spDistsN1(gr_dados,gr[i,]) } betas <- parametros$beta mu <- as.vector(X%*%betas) s2 <- parametros$sigmasq phi1 <- parametros$phi kappa1 <- kap R <- (((2^(kappa1-1))*gamma(kappa1))^(-1)) * ((dists/phi1)^(kappa1)) * besselK(dists/phi1,nu=kappa1) R[dists < 0.005] <- 1 Sig <- s2*R[k4,] diag(Sig) <- diag(Sig) + nug Sig_c <- s2*R[k3,] m_z_dadoy <- mu[k3] - (Sig_c%*%ginv(Sig)%*%(mu[k4]-dados)) pondera_var <- Sig_c%*%tcrossprod(ginv(Sig),Sig_c) saida <- list() saida$pred <- m_z_dadoy saida$var_pred <- s2 - diag(pondera_var) return(saida) } ## Preparando constantes minimo<-min(ph$borders[,1]) minimo maximo<-max(ph$borders[,1]) minimoy<-min(ph$borders[,2]) maximoy<-max(ph$borders[,2]) gr<-expand.grid(seq(minimo,maximo,l=100),seq(minimoy,maximoy,l=100)) gr_pred <- matrix(locations.inside(locations=gr,borders=ph$borders,as.is=F),ncol=2) gr_pred_past <- matrix(locations.inside(locations=gr,borders=ph$reg_past,as.is=F),ncol=2) gr_pred_soja <- matrix(locations.inside(locations=gr,borders=ph$reg_soja,as.is=F),ncol=2) colar <- function(obj){ paste(obj,collapse = "") } b<-apply(gr,1,colar) d<-apply(gr_pred,1,colar) e<-apply(gr_pred_past,1,colar) f<-apply(gr_pred_soja,1,colar) controle1 <- is.element(b,d) controle2 <- is.element(d,e) controle3 <- is.element(d,f) npr <- nrow(gr_pred) X_pred <- matrix(0,nrow=npr,ncol=2) X_pred[,1] <- 1 X_pred[,2] <- gr_pred[,1] ## Krigagem pH kr_ph <- pred(geo1=ph,gr_pred=gr_pred,X_pred=X_pred,parametros=emv_ph3_5,kap=2.5,nug=emv_ph3_5$nugget) ## Krigagem Saturação kr_sat <- pred(geo1=sat,gr_pred=gr_pred,X_pred=X_pred,parametros=emv_sat3_5,kap=2.5,nug=emv_sat3_5$nugget) ## Análise descritiva das krigagens summary(kr_ph$pred) summary(kr_ph$var_pred) summary(kr_sat$pred) summary(kr_sat$var_pred) ## Análise de resíduos par(mfrow=c(1,2),pty="s") qqnorm(emv_sat3_5$model.components$residuals) qqline(emv_sat3_5$model.components$residuals) plot(density(emv_sat3_5$model.components$residuals)) shapiro.test(emv_sat3_5$model.components$residuals) par(mfrow=c(1,2),pty="s") qqnorm(emv_ph3_5$model.components$residuals) qqline(emv_ph3_5$model.components$residuals) plot(density(emv_ph3_5$model.components$residuals)) shapiro.test(emv_ph3_5$model.components$residuals) ## Modelagem bivariada ## Funções para estimação por máxima verossimilhança dos parâmetros da log(veros.) concentrada # BGCCM est_bgccm<-function(p,Ct,rui=F,kap){ if ( any(p < 1e-20)) return(.Machine$double.xmax^.5) if ( any(p > 1000)) return(.Machine$double.xmax^.5) k1 <- c(1:Ct$n1) k2 <- c((Ct$n1+1):Ct$n) R0 <- (((2^(kap[1]-1))*gamma(kap[1]))^(-1)) * ((Ct$dists/p[4])^(kap[1])) * besselK(Ct$dists/p[4],nu=kap[1]) R0[ Ct$dists < 0.005 ] <- 1 R1 <- (((2^(kap[2]-1))*gamma(kap[2]))^(-1)) * ((Ct$dists/p[5])^(kap[2])) * besselK(Ct$dists/p[5],nu=kap[2]) R1[ Ct$dists < 0.005 ] <- 1 R2 <- (((2^(kap[3]-1))*gamma(kap[3]))^(-1)) * ((Ct$dists/p[6])^(kap[3])) * besselK(Ct$dists/p[6],nu=kap[3]) R2[ Ct$dists < 0.005 ] <- 1 V1 <- R0[k1,k1] + p[2]^2 * R1[k1,k1] V2 <- p[1]^2 * R0[k2,k2] + p[3]^2 * R2[k2,k2] V3 <- p[1] * R0 [k1,k2] V <- cbind(rbind(V1,t(V3)),rbind(V3,V2)) if(rui == TRUE) {diag(V) <- diag(V) + c(rep((p[7]^2),Ct$n1),rep((p[8]^2),Ct$n2))} Vinv <- ginv(V) tVinvX <- crossprod(Vinv,Ct$X) Q_hat <- drop(sum(crossprod(Ct$dados,Vinv)*Ct$dados) - crossprod(Ct$dados,tVinvX)%*%ginv(crossprod(tVinvX,Ct$X))%*%crossprod(tVinvX,Ct$dados)) lQ_hat <- log(Q_hat) ldv <- log(det(V)) if(is.infinite(lQ_hat) | is.nan(lQ_hat) | is.infinite(ldv) | is.nan(ldv)) return(.Machine$double.xmax^.5) ll <- drop(-0.5*((Ct$n1+Ct$n2)*(log(2*pi)+lQ_hat-log(Ct$n1+Ct$n2)+1)+ldv)) if(is.infinite(ll)) return(.Machine$double.xmax^.5) if(is.nan(ll)) return(.Machine$double.xmax^.5) return(-ll) } # BCRM est_bcrm<-function(p,Ct,rui=F,kap,repet){ if ( any(p < 1e-20)) return(.Machine$double.xmax^.5) if ( any(p > 1000)) return(.Machine$double.xmax^.5) k1 <- c(1:Ct$n1) k2 <- c((Ct$n1+1):Ct$n) R1 <- (((2^(kap[1]-1))*gamma(kap[1]))^(-1)) * ((Ct$dists/p[3])^(kap[1])) * besselK(Ct$dists/p[3],nu=kap[1]) R1[ Ct$dists < 0.005 ] <- 1 R2 <- (((2^(kap[2]-1))*gamma(kap[2]))^(-1)) * ((Ct$dists/p[4])^(kap[2])) * besselK(Ct$dists/p[4],nu=kap[2]) R2[ Ct$dists < 0.005 ] <- 1 V1 <- R1[k1,k1] V2 <- p[1]^2 * R1[k2,k2] + p[2]^2 * R2[k2,k2] V3 <- p[1] * R1 [k1,k2] V <- cbind(rbind(V1,t(V3)),rbind(V3,V2)) if(rui == TRUE) {diag(V) <- diag(V) + c(rep((p[7]^2),Ct$n1),rep((p[8]^2),Ct$n2))} Vinv <- ginv(V) tVinvX <- crossprod(Vinv,Ct$X) Q_hat <- drop(sum(crossprod(Ct$dados,Vinv)*Ct$dados) - crossprod(Ct$dados,tVinvX)%*%ginv(crossprod(tVinvX,Ct$X))%*%crossprod(tVinvX,Ct$dados)) lQ_hat <- log(Q_hat) ldv <- log(det(V)) if(is.infinite(lQ_hat) | is.nan(lQ_hat) | is.infinite(ldv) | is.nan(ldv)) return(.Machine$double.xmax^.5) ll <- drop(-0.5*((Ct$n1+Ct$n2)*(log(2*pi)+lQ_hat-log(Ct$n1+Ct$n2)+1)+ldv)) if(is.infinite(ll)) return(.Machine$double.xmax^.5) if(is.nan(ll)) return(.Machine$double.xmax^.5) return(-ll) } ## Preparando entradas das funções y1<-sat y2<-ph ctes <- list() ctes$n1 <- length(y1$data) ctes$n2 <- length(y2$data) ctes$n <- ctes$n1 + ctes$n2 ctes$dados <- c(y1$data,y2$data) ctes$dists <- unname(as.matrix(dist(rbind(y1$coords,y2$coords),diag=TRUE,upper=TRUE))) auxiliar <- as.numeric(ph$covariate$AREA) auxiliar[auxiliar==1]<-0 auxiliar[auxiliar==2]<-1 # Para construção da matriz de delineamento e definição dos kappas # diversas combinações podem ser levadas em consideração, # Devido ao tempo computacional exigido, segue apenas um exemplo, # Esse exemplo é facilmente adaptável # Exemplo para o BGCCM # Entrando com kappa0, kappa1, kappa2 kappas <- matrix(0,36,3) kappas[1,] <- c(0.5,2.5,2.5) # Entrando com a matriz de delineamento # Nesse caso, média constante para as duas respostas ctes$X<-matrix(c(rep(1,67),rep(0,67),rep(0,67),rep(1,67)),134,2) X <- ctes$X # Valores iniciais para o métodos numérico par.ini <- 0 par.ini[1] <- .5/30 par.ini[2] <- 8/30 par.ini[3] <- .3/30 par.ini[4] <- 25 par.ini[5] <- 15 par.ini[6] <- 10 # Maximizando funçao de verossimilhança concentrada emv_bgccm <- optim(p=par.ini,fn=est_bgccm,Ct=ctes,kap=kappas[1,],control=list(maxit=10000),hessian=T) ## Calculando parâmetros originais kappa0 <- kappas[1,1] kappa1 <- kappas[1,2] kappa2 <- kappas[1,3] dists <- ctes$dists dados <- ctes$dados n1 <- ctes$n1 n2 <- ctes$n2 n <- n1 + n2 k1 <- c(1:n1) k2 <- c((n1+1):n) R0 <- (((2^(kappa0-1))*gamma(kappa0))^(-1)) * ((dists/emv_bgccm$par[4])^(kappa0)) * besselK(dists/emv_bgccm$par[4],nu=kappa0) R0[dists<0.005] <- 1 R1 <- (((2^(kappa1-1))*gamma(kappa1))^(-1)) * ((dists[k1,k1]/emv_bgccm$par[5])^(kappa1)) * besselK(dists[k1,k1]/emv_bgccm$par[5],nu=kappa1) R1[dists[k1,k1]<0.005] <- 1 R2 <- (((2^(kappa2-1))*gamma(kappa2))^(-1)) * ((dists[k2,k2]/emv_bgccm$par[6])^(kappa2)) * besselK(dists[k2,k2]/emv_bgccm$par[6],nu=kappa2) R2[dists[k2,k2]<0.005] <- 1 V1 <- R0[k1,k1] + emv_bgccm$par[2]^2 * R1 V2 <- emv_bgccm$par[1]^2 * R0[k2,k2] + emv_bgccm$par[3]^2 * R2 V3 <- emv_bgccm$par[1] * R0[k1,k2] V <- cbind(rbind(V1,t(V3)), rbind(V3,V2)) mu <- drop(ginv(crossprod(X,ginv(V))%*%X)%*% crossprod(X,ginv(V))%*%dados) sigmasq <- drop((1/n) * crossprod(dados-X%*%mu,ginv(V)%*%dados-X%*%mu)) sigma01 <- sqrt(sigmasq) sigma02 <- emv_bgccm$par[1] *sigma01 sigma1 <- emv_bgccm$par[2] * sigma01 sigma2 <- emv_bgccm$par[3] * sigma01 sigmas <- c(sigma01,sigma02,sigma1,sigma2) para_bgccm <- round(c(mu,sigmas,emv_bgccm$par[4],emv_bgccm$par[5],emv_bgccm$par[6],emv_bgccm$value,emv_bgccm$convergence),5) para_bgccm # Exemplo do emv BCRM # Valores iniciais par.ini <- 0 par.ini[1] <- .5/35 par.ini[2] <- .3/35 par.ini[3] <- 30 par.ini[4] <- 20 # Entrando com kappa1 e kappa2 kap<- c(0.5,0.5) # Pode-se mudar a matriz X tb, neste caso está considerada ainda médias constantes # Maximização da função de verossimilhança concentrada emv_bcrm <- optim(p=par.ini,fn=est_bcrm,Ct=ctes,kap=kap,control=list(maxit=10000),hessian=T) # Calculando parâmetros originais kappa1 <- kappas[1,1] kappa2 <- kappas[1,2] dists <- ctes$dists dados <- ctes$dados n1 <- ctes$n1 n2 <- ctes$n2 n <- n1 + n2 k1 <- c(1:n1) k2 <- c((n1+1):n) R1 <- (((2^(kappa1-1))*gamma(kappa1))^(-1)) * ((dists/emv_bcrm$par[3])^(kappa1)) * besselK(dists/emv_bcrm$par[3],nu=kappa1) R1[dists<0.005] <- 1 R2 <- (((2^(kappa2-1))*gamma(kappa2))^(-1)) * ((dists[k2,k2]/emv_bcrm$par[4])^(kappa2)) * besselK(dists[k2,k2]/emv_bcrm$par[4],nu=kappa2) R2[dists[k2,k2]<0.005] <- 1 V1 <- R1[k1,k1] V2 <- emv_bcrm$par[1]^2 * R1[k2,k2] + emv_bcrm$par[2]^2 * R2 V3 <- emv_bcrm$par[1] * R1[k1,k2] V <- cbind(rbind(V1,t(V3)), rbind(V3,V2)) mu <- drop(ginv(crossprod(X,ginv(V))%*%X)%*% crossprod(X,ginv(V))%*%dados) sigmasq <- drop((1/n) * crossprod(dados-X%*%mu,ginv(V)%*%dados-X%*%mu)) sigma11 <- sqrt(sigmasq) sigma12 <- emv_bcrm$par[1] *sigma11 sigma22 <- emv_bcrm$par[2] * sigma11 sigmas <- c(sigma11,sigma12,sigma22) para_bcrm <- round(c(mu,sigmas,emv_bcrm$par[3],emv_bcrm$par[4],emv_bcrm$value,emv_bcrm$convergence),5) para_bcrm ### Modelos finais ## Médias com tendência induzida pela área de manejo para as duas variáveis ## beta = (beta01,beta01,beta11,beta12) # BGCCM resultados_bgccm <- list() resultados_bgccm$betas <- c(beta0_sat=47.8336,beta1_sat=9.4583,beta0_ph=4.7216,beta1_ph=0.4048) resultados_bgccm$sigmas <- c(sigma01=6.8499,sigma1=3.4358,sigma02=0.3109,sigma2=1e-5) resultados_bgccm$phis <- c(phi0=26.7319,phi1=49.0824,phi2=86.9527) resultados_bgccm$kappas <- c(kappa0=1.5,kappa1=0.5,kappa2=0.5) resultados_bgccm$loglik <- 207.67317 # BCRM resultados_bcrm <- list() resultados_bcrm$betas <- c(beta0_sat=47.9261,beta1_sat=9.1249,beta0_ph=4.7177,beta1_ph=0.3961) resultados_bcrm$sigmas <- c(sigma11=7.8187,sigma12=0.2818,sigma22=0.1393) resultados_bcrm$phis <- c(phi1=53.3170,phi2=19.3394) resultados_bcrm$kappas <- c(kappa1=0.5,kappa2=2.5) resultados_bcrm$loglik <- 207.4801 ### Krigagem bivariada # Função para krigagem no BGCCM considerando os modelos finais predBGCCM <- function(geo1,geo2,gr_pred,X_pred,parametros,kap){ coords1 <- geo1$coords coords2 <- geo2$coords auxiliar1 <- as.numeric(geo1$covariate$AREA) auxiliar1[auxiliar1==1]<-0 auxiliar1[auxiliar1==2]<-1 auxiliar2 <- as.numeric(geo2$covariate$AREA) auxiliar2[auxiliar2==1]<-0 auxiliar2[auxiliar2==2]<-1 n1 <- nrow(coords1) n2 <- nrow(coords2) n <- n1+n2 npr <- nrow(gr_pred) nt <- npr + n k1 <- c(1:n1) k2 <- c(1:n2) k3 <- c(1:npr) k4 <- c((npr+1):(npr+n1)) k5 <- c((npr+1):(2*npr)) k6 <- c((npr+n1+1):nt) k7 <- c((n1+1):n) k8 <- c((npr+1):nt) k9 <- c(((2*npr)+1):((2*npr)+n)) dados <- c(geo1$data,geo2$data) X_dados <- matrix(c(rep(1,n1),rep(0,n2),rep(0,n1),rep(1,n2),auxiliar1,rep(0,n1),rep(0,n2),auxiliar2),n,4) X <- rbind(X_pred,X_dados) gr_dados <- unname(rbind(coords1,coords2)) gr <- rbind(gr_pred,gr_dados) dists <- matrix(0,nt,n) for(i in 1:nt){ dists[i,] <- spDistsN1(gr_dados,gr[i,]) } betas <- parametros[1:4] mu <- as.vector(X%*%betas) s01 <- parametros[5] s1 <- parametros[7] s02 <- parametros[6] s2 <- parametros[8] phi0 <- parametros[9] phi1 <- parametros[10] phi2 <- parametros[11] kappa0 <- kap[1] kappa1 <- kap[2] kappa2 <- kap[3] R0 <- (((2^(kappa0-1))*gamma(kappa0))^(-1)) * ((dists/phi0)^(kappa0)) * besselK(dists/phi0,nu=kappa0) R0[dists < 0.005] <- 1 R1 <- (((2^(kappa1-1))*gamma(kappa1))^(-1)) * ((dists[k4,k1]/phi1)^(kappa1)) * besselK(dists[k4,k1]/phi1,nu=kappa1) distaux <- dists[k4,k1] R1[distaux < 0.005] <- 1 R2 <- (((2^(kappa2-1))*gamma(kappa2))^(-1)) * ((dists[k6,k7]/phi2)^(kappa2)) * besselK(dists[k6,k7]/phi2,nu=kappa2) distaux <- dists[k6,k7] R2[distaux < 0.005] <- 1 Sig1 <- (s01^2)*R0[k4,k1] + (s1^2)*R1 Sig2 <- (s02^2)*R0[k6,k7] + (s2^2)*R2 Sig12 <- (s01*s02)*R0[k4,k7] Sig <- cbind(rbind(Sig1,t(Sig12)),rbind(Sig12,Sig2)) R1c <- (((2^(kappa1-1))*gamma(kappa1))^(-1)) * ((dists[k3,k1]/phi1)^(kappa1)) * besselK(dists[k3,k1]/phi1,nu=kappa1) distaux <- dists[k3,k1] R1c[distaux < 0.005] <- 1 Sig_c1 <- (s01^2)*R0[k3,k1] + (s1^2)*R1c Sig_c12 <- (s01*s02) * R0[k3,k7] Sig_c <- cbind(Sig_c1,Sig_c12) m_z_dadoy <- mu[k3] - (Sig_c%*%ginv(Sig)%*%(mu[k9]-dados)) pondera_var <- Sig_c%*%tcrossprod(ginv(Sig),Sig_c) saida <- list() saida$pred1 <- m_z_dadoy saida$var_pred1 <- ((s01^2) + (s1^2)) - diag(pondera_var) R2c <- (((2^(kappa2-1))*gamma(kappa2))^(-1)) * ((dists[k3,k7]/phi2)^(kappa2)) * besselK(dists[k3,k7]/phi2,nu=kappa2) distaux <- dists[k3,k7] R2c[distaux < 0.005] <- 1 Sig_c2 <- (s02^2)*R0[k3,k7] + (s2^2)*R2c Sig_c21 <- (s01*s02) * R0[k3,k1] Sig_c <- cbind(Sig_c21,Sig_c2) m_z_dadoy <- mu[k5] - (Sig_c%*%ginv(Sig)%*%(mu[k9]-dados)) pondera_var <- Sig_c%*%tcrossprod(ginv(Sig),Sig_c) saida$pred2 <- m_z_dadoy saida$var_pred2 <- ((s02^2) + (s2^2)) - diag(pondera_var) return(saida) } # Função para krigagem no BCRM predBCRM <- function(geo1,geo2,gr_pred,X_pred,parametros,kap){ coords1 <- geo1$coords coords2 <- geo2$coords auxiliar1 <- as.numeric(geo1$covariate$AREA) auxiliar1[auxiliar1==1]<-0 auxiliar1[auxiliar1==2]<-1 auxiliar2 <- as.numeric(geo2$covariate$AREA) auxiliar2[auxiliar2==1]<-0 auxiliar2[auxiliar2==2]<-1 n1 <- nrow(coords1) n2 <- nrow(coords2) n <- n1+n2 npr <- nrow(gr_pred) nt <- npr + n k1 <- c(1:n1) k2 <- c(1:n2) k3 <- c(1:npr) k4 <- c((npr+1):(npr+n1)) k5 <- c((npr+1):(2*npr)) k6 <- c((npr+n1+1):nt) k7 <- c((n1+1):n) k8 <- c((npr+1):nt) k9 <- c(((2*npr)+1):((2*npr)+n)) dados <- c(geo1$data,geo2$data) X_dados <- matrix(c(rep(1,n1),rep(0,n2),rep(0,n1),rep(1,n2),auxiliar1,rep(0,n1),rep(0,n2),auxiliar2),n,4) X <- rbind(X_pred,X_dados) gr_dados <- unname(rbind(coords1,coords2)) gr <- rbind(gr_pred,gr_dados) dists <- matrix(0,nt,n) for(i in 1:nt){ dists[i,] <- spDistsN1(gr_dados,gr[i,]) } betas <- parametros[1:4] mu <- as.vector(X%*%betas) s11 <- parametros[5] s21 <- parametros[6] s22 <- parametros[7] phi1 <- parametros[8] phi2 <- parametros[9] kappa1 <- kap[1] kappa2 <- kap[2] R1 <- (((2^(kappa1-1))*gamma(kappa1))^(-1)) * ((dists/phi1)^(kappa1)) * besselK(dists/phi1,nu=kappa1) R1[dists < 0.005] <- 1 R2 <- (((2^(kappa2-1))*gamma(kappa2))^(-1)) * ((dists[k6,k7]/phi2)^(kappa2)) * besselK(dists[k6,k7]/phi2,nu=kappa2) distaux <- dists[k6,k7] R2[distaux < 0.005] <- 1 Sig1 <- (s11^2)*R1[k4,k1] Sig2 <- (s21^2)*R1[k6,k7] + (s22^2)*R2 Sig12 <- (s11*s21)*R1[k4,k7] Sig <- cbind(rbind(Sig1,t(Sig12)),rbind(Sig12,Sig2)) Sig_c1 <- (s11^2)*R1[k3,k1] Sig_c12 <- (s11*s21) * R1[k3,k7] Sig_c <- cbind(Sig_c1,Sig_c12) m_z_dadoy <- mu[k3] - (Sig_c%*%ginv(Sig)%*%(mu[k9]-dados)) pondera_var <- Sig_c%*%tcrossprod(ginv(Sig),Sig_c) saida <- list() saida$pred1 <- m_z_dadoy saida$var_pred1 <- (s11^2) - diag(pondera_var) R2c <- (((2^(kappa2-1))*gamma(kappa2))^(-1)) * ((dists[k3,k7]/phi2)^(kappa2)) * besselK(dists[k3,k7]/phi2,nu=kappa2) distaux <- dists[k3,k7] R2c[distaux < 0.005] <- 1 Sig_c2 <- (s21^2)*R1[k3,k7] + (s22^2)*R2c Sig_c21 <- (s11*s21) * R1[k3,k1] Sig_c <- cbind(Sig_c21,Sig_c2) m_z_dadoy <- mu[k5] - (Sig_c%*%ginv(Sig)%*%(mu[k9]-dados)) pondera_var <- Sig_c%*%tcrossprod(ginv(Sig),Sig_c) saida$pred2 <- m_z_dadoy saida$var_pred2 <- ((s21^2) + (s22^2)) - diag(pondera_var) return(saida) } # Preparando entradas minimo<-min(ph$borders[,1]) minimo maximo<-max(ph$borders[,1]) minimoy<-min(ph$borders[,2]) maximoy<-max(ph$borders[,2]) gr<-expand.grid(seq(minimo,maximo,l=100),seq(minimoy,maximoy,l=100)) gr_pred <- matrix(locations.inside(locations=gr,borders=ph$borders,as.is=F),ncol=2) gr_pred_past <- matrix(locations.inside(locations=gr,borders=ph$reg_past,as.is=F),ncol=2) gr_pred_soja <- matrix(locations.inside(locations=gr,borders=ph$reg_soja,as.is=F),ncol=2) colar <- function(obj){ paste(obj,collapse = "") } b<-apply(gr,1,colar) d<-apply(gr_pred,1,colar) e<-apply(gr_pred_past,1,colar) f<-apply(gr_pred_soja,1,colar) controle1 <- is.element(b,d) controle2 <- is.element(d,e) controle3 <- is.element(d,f) npr <- nrow(gr_pred) k1 <- c(1:npr) k2 <- c((npr+1):(2*npr)) X_pred <- matrix(0,nrow=2*npr,ncol=4) X_pred[k1,1] <- 1 X_pred[k2,2] <- 1 auxiliar <- rep(0,npr) auxiliar[controle3 ==T] <- 1 X_pred[k1,3] <- auxiliar X_pred[k2,4] <- auxiliar ## A matriz de delineamento está considerando a área de manejo ## Krigagem do BGCCM #Parâmetros do modelo final com o BGCCM para_aux1 <- c(47.8336,4.7216,9.4583,0.4048,6.8499,0.3109,3.4358,1e-5,26.7319,49.0824,86.9527) kappas <- c(1.5,.5,.5) krig_bgccm <- predBGCCM(geo1=sat,geo2=ph,gr_pred=gr_pred,X_pred=X_pred,parametros=para_aux1,kap=kappas) ## Krigagem do BCRM #Parâmetros do modelo final com o BCRM para_aux2 <- c(47.9261,4.7177,9.1249,0.3961,7.8187,0.2818,0.1393,53.3170,19.3394) kappas <- c(0.5,2.5) krig_bcrm <- predBCRM(geo1=sat,geo2=ph,gr_pred=gr_pred,X_pred=X_pred,parametros=para_aux2,kap=kappas) ## Descrição dos resultados do BGCCM ## Saturação summary(krig_bgccm$pred1) summary(krig_bgccm$var_pred1) ## pH summary(krig_bgccm$pred2) summary(krig_bgccm$var_pred2) ## Descrição dos resultados do BCRM ## Saturação summary(krig_bcrm$pred1) summary(krig_bcrm$var_pred1) ## pH summary(krig_bcrm$pred2) summary(krig_bcrm$var_pred2) ### Mostrando o funcionamento dos modelos bivariados com um menor número de observações ### da saturação por bases # Retirando alguns dados da saturação por bases sat_aux <- sat sat_aux$data <- sat$data[-c(1,4,18,21,23,37,40,43,56,67,7,10,13,25,28,32,44,49,62,64)] sat_aux$coords <- sat$coords[-c(1,4,18,21,23,37,40,43,56,67,7,10,13,25,28,32,44,49,62,64),] sat_aux$covariate <- NULL sat_aux$covariate <- data.frame(AREA=sat$covariate[-c(1,4,18,21,23,37,40,43,56,67,7,10,13,25,28,32,44,49,62,64),]) # Gráficos comparando a saturação com todos os pontos e com a retirada de algumas locs. points(sat) lines(sat$reg1) x11() points(sat_aux) lines(sat_aux$reg1) # Estimação de máxima verossimilhança da saturação com menos localizações emv_sat_aux <-likfit(sat_aux,trend=~coords[,1],ini.cov.pars=c(65,190.39),nugget=45.41, fix.kappa=TRUE,kappa=2.5) # krigagem da saturação com o modelo univariado minimo<-min(ph$borders[,1]) minimo maximo<-max(ph$borders[,1]) minimoy<-min(ph$borders[,2]) maximoy<-max(ph$borders[,2]) gr<-expand.grid(seq(minimo,maximo,l=100),seq(minimoy,maximoy,l=100)) gr_pred <- matrix(locations.inside(locations=gr,borders=ph$borders,as.is=F),ncol=2) gr_pred_past <- matrix(locations.inside(locations=gr,borders=ph$reg_past,as.is=F),ncol=2) gr_pred_soja <- matrix(locations.inside(locations=gr,borders=ph$reg_soja,as.is=F),ncol=2) colar <- function(obj){ paste(obj,collapse = "") } b<-apply(gr,1,colar) d<-apply(gr_pred,1,colar) e<-apply(gr_pred_past,1,colar) f<-apply(gr_pred_soja,1,colar) controle1 <- is.element(b,d) controle2 <- is.element(d,e) controle3 <- is.element(d,f) npr <- nrow(gr_pred) X_pred <- matrix(0,nrow=npr,ncol=2) X_pred[,1] <- 1 X_pred[,2] <- gr_pred[,1] kr_sat_aux <- pred(geo1=sat_aux,gr_pred=gr_pred,X_pred=X_pred,parametros=emv_sat_aux, kap=2.5,nug=emv_sat_aux$nugget) # Comparando krigagens univariadas utilizando todos as locs e com 20 a menos summary(kr_sat$pred) summary(kr_sat_aux$pred) summary(kr_sat$var_pred) summary(kr_sat_aux$var_pred) # krigagem nos pontos retirados utilizando o modelo univariado kr_pts_ret<-pred(geo1=sat_aux,gr_pred=sat$coords[c(1,4,18,21,23,37,40,43,56,67,7,10,13,25,28,32,44,49,62,64),], X_pred=X_pred,parametros=emv_sat_aux,kap=2.5,nug=emv_sat_aux$nugget) # Comparando valores observados e preditos cbind(sat$data[c(1,4,18,21,23,37,40,43,56,67,7,10,13,25,28,32,44,49,62,64)],kr_pts_ret$pred) ## Ajustando modelos bivariados utilizando a saturação retirando alguns pontos ## Preparando entradas para estimação y1<-sat_aux y2<-ph ctes <- list() ctes$n1 <- length(y1$data) ctes$n2 <- length(y2$data) ctes$n <- ctes$n1 + ctes$n2 ctes$dados <- c(y1$data,y2$data) ctes$dists <- unname(as.matrix(dist(rbind(y1$coords,y2$coords),diag=TRUE,upper=TRUE))) auxiliar1 <- as.numeric(ph$covariate$AREA) auxiliar1[auxiliar1==1]<-0 auxiliar1[auxiliar1==2]<-1 auxiliar2 <- as.numeric(sat_aux$covariate$AREA) auxiliar2[auxiliar2==1]<-0 auxiliar2[auxiliar2==2]<-1 ctes$X<-matrix(c(rep(1,ctes$n1),rep(0,ctes$n2),rep(0,ctes$n1),rep(1,ctes$n2), auxiliar2,rep(0,ctes$n2),rep(0,ctes$n1),auxiliar1),(ctes$n1+ctes$n2),4) X <- ctes$X ## Estimação do BGCCM par.ini <- 0 par.ini[1] <- .5/30 par.ini[2] <- 8/30 par.ini[3] <- .3/30 par.ini[4] <- 25 par.ini[5] <- 15 par.ini[6] <- 10 kappas <- matrix(0,36,3) kappas[1,] <- c(1.5,0.5,0.5) emv_bgccm <- optim(p=par.ini,fn=est_bgccm,Ct=ctes,kap=kappas[1,],control=list(maxit=10000),hessian=F) # Calculando parâmetros originais kappa0 <- kappas[1,1] kappa1 <- kappas[1,2] kappa2 <- kappas[1,3] dists <- ctes$dists dados <- ctes$dados n1 <- ctes$n1 n2 <- ctes$n2 n <- n1 + n2 k1 <- c(1:n1) k2 <- c((n1+1):n) R0 <- (((2^(kappa0-1))*gamma(kappa0))^(-1)) * ((dists/emv_bgccm$par[4])^(kappa0)) * besselK(dists/emv_bgccm$par[4],nu=kappa0) R0[dists<0.005] <- 1 R1 <- (((2^(kappa1-1))*gamma(kappa1))^(-1)) * ((dists[k1,k1]/emv_bgccm$par[5])^(kappa1)) * besselK(dists[k1,k1]/emv_bgccm$par[5],nu=kappa1) R1[dists[k1,k1]<0.005] <- 1 R2 <- (((2^(kappa2-1))*gamma(kappa2))^(-1)) * ((dists[k2,k2]/emv_bgccm$par[6])^(kappa2)) * besselK(dists[k2,k2]/emv_bgccm$par[6],nu=kappa2) R2[dists[k2,k2]<0.005] <- 1 V1 <- R0[k1,k1] + emv_bgccm$par[2]^2 * R1 V2 <- emv_bgccm$par[1]^2 * R0[k2,k2] + emv_bgccm$par[3]^2 * R2 V3 <- emv_bgccm$par[1] * R0[k1,k2] V <- cbind(rbind(V1,t(V3)), rbind(V3,V2)) mu <- drop(ginv(crossprod(X,ginv(V))%*%X)%*% crossprod(X,ginv(V))%*%dados) sigmasq <- drop((1/n) * crossprod(dados-X%*%mu,ginv(V)%*%dados-X%*%mu)) sigma01 <- sqrt(sigmasq) sigma02 <- emv_bgccm$par[1] *sigma01 sigma1 <- emv_bgccm$par[2] * sigma01 sigma2 <- emv_bgccm$par[3] * sigma01 sigmas <- c(sigma01,sigma02,sigma1,sigma2) para_bgccm <- round(c(mu,sigmas,emv_bgccm$par[4],emv_bgccm$par[5],emv_bgccm$par[6],emv_bgccm$value,emv_bgccm$convergence),5) para_bgccm # Estimação do BCRM par.ini <- 0 par.ini[1] <- .5/35 par.ini[2] <- .3/35 par.ini[3] <- 30 par.ini[4] <- 20 kappas <- matrix(0,25,2) kappas[1,] <- c(0.5,2.5) emv_bcrm <- optim(p=par.ini,fn=est_bcrm,Ct=ctes,kap=kappas[1,],control=list(maxit=10000),hessian=F) # Calculando parâmetros originais kappa1 <- kappas[1,1] kappa2 <- kappas[1,2] dists <- ctes$dists dados <- ctes$dados n1 <- ctes$n1 n2 <- ctes$n2 n <- n1 + n2 k1 <- c(1:n1) k2 <- c((n1+1):n) R1 <- (((2^(kappa1-1))*gamma(kappa1))^(-1)) * ((dists/emv_bcrm$par[3])^(kappa1)) * besselK(dists/emv_bcrm$par[3],nu=kappa1) R1[dists<0.005] <- 1 R2 <- (((2^(kappa2-1))*gamma(kappa2))^(-1)) * ((dists[k2,k2]/emv_bcrm$par[4])^(kappa2)) * besselK(dists[k2,k2]/emv_bcrm$par[4],nu=kappa2) R2[dists[k2,k2]<0.005] <- 1 V1 <- R1[k1,k1] V2 <- emv_bcrm$par[1]^2 * R1[k2,k2] + emv_bcrm$par[2]^2 * R2 V3 <- emv_bcrm$par[1] * R1[k1,k2] V <- cbind(rbind(V1,t(V3)), rbind(V3,V2)) mu <- drop(ginv(crossprod(X,ginv(V))%*%X)%*% crossprod(X,ginv(V))%*%dados) sigmasq <- drop((1/n) * crossprod(dados-X%*%mu,ginv(V)%*%dados-X%*%mu)) sigma11 <- sqrt(sigmasq) sigma12 <- emv_bcrm$par[1] *sigma11 sigma22 <- emv_bcrm$par[2] * sigma11 sigmas <- c(sigma11,sigma12,sigma22) para_bcrm <- round(c(mu,sigmas,emv_bcrm$par[3],emv_bcrm$par[4],emv_bcrm$value,emv_bcrm$convergence),5) para_bcrm ### Krigagem ## Preparando entradas minimo<-min(ph$borders[,1]) minimo maximo<-max(ph$borders[,1]) minimoy<-min(ph$borders[,2]) maximoy<-max(ph$borders[,2]) gr<-expand.grid(seq(minimo,maximo,l=100),seq(minimoy,maximoy,l=100)) gr_pred <- matrix(locations.inside(locations=gr,borders=ph$borders,as.is=F),ncol=2) gr_pred_past <- matrix(locations.inside(locations=gr,borders=ph$reg_past,as.is=F),ncol=2) gr_pred_soja <- matrix(locations.inside(locations=gr,borders=ph$reg_soja,as.is=F),ncol=2) colar <- function(obj){ paste(obj,collapse = "") } b<-apply(gr,1,colar) d<-apply(gr_pred,1,colar) e<-apply(gr_pred_past,1,colar) f<-apply(gr_pred_soja,1,colar) controle1 <- is.element(b,d) controle2 <- is.element(d,e) controle3 <- is.element(d,f) npr <- nrow(gr_pred) k1 <- c(1:npr) k2 <- c((npr+1):(2*npr)) X_pred <- matrix(0,nrow=2*npr,ncol=4) X_pred[k1,1] <- 1 X_pred[k2,2] <- 1 auxiliar <- rep(0,npr) auxiliar[controle3 ==T] <- 1 X_pred[k1,3] <- auxiliar X_pred[k2,4] <- auxiliar ## Krigagem no BGCCM kappas <- c(1.5,.5,.5) krig_teste1 <- predBGCCM(geo1=sat_aux,geo2=ph,gr_pred=gr_pred,X_pred=X_pred,parametros=para_bgccm,kap=kappas) krig_bgccm_aux <- krig_teste1 ## Krigagem no BCRM kappas <- c(0.5,2.5) krig_teste2 <- predBCRM(geo1=sat_aux,geo2=ph,gr_pred=gr_pred,X_pred=X_pred,parametros=para_bcrm,kap=kappas) krig_bcrm_aux <- krig_teste2 ## Comparando krigagens bivariadas com saturação completa e com alguns pts retirados ## BGCCM # Saturação summary(krig_bgccm$pred1) summary(krig_bgccm_aux$pred1) summary(krig_bgccm$var_pred1) summary(krig_bgccm_aux$var_pred1) #pH summary(krig_bgccm$pred2) summary(krig_bgccm_aux$pred2) summary(krig_bgccm$var_pred2) summary(krig_bgccm_aux$var_pred2) ## BGCCM # Saturação summary(krig_bcrm$pred1) summary(krig_bcrm_aux$pred1) summary(krig_bcrm$var_pred1) summary(krig_bcrm_aux$var_pred1) #pH summary(krig_bcrm$pred2) summary(krig_bcrm_aux$pred2) summary(krig_bcrm$var_pred2) summary(krig_bcrm_aux$var_pred2) ### Krigando a saturação nos pontos retirados sob os modelos bivariados gr_aux <- sat$coords[c(1,4,18,21,23,37,40,43,56,67,7,10,13,25,28,32,44,49,62,64),] npr <- nrow(gr_aux) gr_pred_past <- matrix(locations.inside(locations=gr_aux,borders=ph$reg_past,as.is=F),ncol=2) gr_pred_soja <- matrix(locations.inside(locations=gr_aux,borders=ph$reg_soja,as.is=F),ncol=2) colar <- function(obj){ paste(obj,collapse = "") } d<-apply(gr_aux,1,colar) e<-apply(gr_pred_past,1,colar) f<-apply(gr_pred_soja,1,colar) controle2 <- is.element(d,e) controle3 <- is.element(d,f) k1 <- c(1:npr) k2 <- c((npr+1):(2*npr)) X_pred <- matrix(0,nrow=2*npr,ncol=4) X_pred[k1,1] <- 1 X_pred[k2,2] <- 1 auxiliar <- rep(0,npr) auxiliar[controle3 ==T] <- 1 X_pred[k1,3] <- auxiliar X_pred[k2,4] <- auxiliar kappas <- c(1.5,.5,.5) kr_pts_ret_bgccm <- predBGCCM(geo1=sat_aux,geo2=ph,gr_pred=gr_aux,X_pred=X_pred,parametros=para_bgccm,kap=kappas) kappas <- c(0.5,2.5) kr_pts_ret_bcrm <- predBCRM(geo1=sat_aux,geo2=ph,gr_pred=gr_aux,X_pred=X_pred,parametros=para_bcrm,kap=kappas) ## Comparando os resultados de krigagem univariado, bivariados com valor observado cbind(sat$data[c(1,4,18,21,23,37,40,43,56,67,7,10,13,25,28,32,44,49,62,64)],kr_pts_ret$pred, kr_pts_ret_bgccm$pred1,kr_pts_ret_bcrm$pred1) dif1 <- sat$data[c(1,4,18,21,23,37,40,43,56,67,7,10,13,25,28,32,44,49,62,64)]-kr_pts_ret$pred dif2 <- sat$data[c(1,4,18,21,23,37,40,43,56,67,7,10,13,25,28,32,44,49,62,64)]-kr_pts_ret_bgccm$pred1 dif3 <- sat$data[c(1,4,18,21,23,37,40,43,56,67,7,10,13,25,28,32,44,49,62,64)]-kr_pts_ret_bcrm$pred1 ## Imagens dos mapas de krigagem sob todas abordagens minimo<-min(ph$borders[,1]) minimo maximo<-max(ph$borders[,1]) minimoy<-min(ph$borders[,2]) maximoy<-max(ph$borders[,2]) gr<-expand.grid(seq(minimo,maximo,l=100),seq(minimoy,maximoy,l=100)) gr_pred <- matrix(locations.inside(locations=gr,borders=ph$borders,as.is=F),ncol=2) gr_pred_past <- matrix(locations.inside(locations=gr,borders=ph$reg_past,as.is=F),ncol=2) gr_pred_soja <- matrix(locations.inside(locations=gr,borders=ph$reg_soja,as.is=F),ncol=2) colar <- function(obj){ paste(obj,collapse = "") } b<-apply(gr,1,colar) d<-apply(gr_pred,1,colar) e<-apply(gr_pred_past,1,colar) f<-apply(gr_pred_soja,1,colar) controle1 <- is.element(b,d) controle2 <- is.element(d,e) controle3 <- is.element(d,f) npr <- nrow(gr_pred) X_pred <- matrix(0,nrow=npr,ncol=2) X_pred[,1] <- 1 X_pred[,2] <- gr_pred[,1] # Sob o modelo univariado kr_ph <- pred(geo1=ph,gr_pred=gr_pred,X_pred=X_pred,parametros=emv_ph3_5,kap=2.5,nug=emv_ph3_5$nugget) kr_sat <- pred(geo1=sat,gr_pred=gr_pred,X_pred=X_pred,parametros=emv_sat3_5,kap=2.5,nug=emv_sat3_5$nugget) aux <- krige.conv(ph,loc=gr,krige=krige.control(obj=emv_sat3_5)) aux$predict <- kr_ph$pred image(aux,xlim=c(minimo-20,maximo+250),ylim=c(minimoy-50,maximoy+50),xlab="Oeste - Leste",ylab="Sul - Norte",col=gray(100:1/100)) legend.krige(x.leg=c(maximo+100,maximo+150), y.leg=c(minimoy+10,maximoy-10),aux$predict, vert=TRUE,col=gray(100:1/100)) x11() aux$predict <- kr_sat$pred image(aux,xlim=c(minimo-20,maximo+250),ylim=c(minimoy-50,maximoy+50),xlab="Oeste - Leste",ylab="Sul - Norte",col=gray(100:1/100)) legend.krige(x.leg=c(maximo+100,maximo+150), y.leg=c(minimoy+10,maximoy-10),aux$predict, vert=TRUE,col=gray(100:1/100)) # Sob os modelos bivariados # BGCCM # Saturação aux$predict <- krig_bgccm$pred1 image(aux,xlim=c(minimo-20,maximo+250),ylim=c(minimoy-50,maximoy+50),xlab="Oeste - Leste",ylab="Sul - Norte",col=gray(100:1/100)) legend.krige(x.leg=c(maximo+100,maximo+150), y.leg=c(minimoy+10,maximoy-10),aux$predict, vert=TRUE,col=gray(100:1/100)) lines(sat$reg1) # ph aux$predict <- krig_bgccm$pred2 image(aux,xlim=c(minimo-20,maximo+250),ylim=c(minimoy-50,maximoy+50),xlab="Oeste - Leste",ylab="Sul - Norte",col=gray(100:1/100)) legend.krige(x.leg=c(maximo+100,maximo+150), y.leg=c(minimoy+10,maximoy-10),aux$predict, vert=TRUE,col=gray(100:1/100)) lines(ph$reg1) # BCRM # Saturação aux$predict <- krig_bcrm$pred1 image(aux,xlim=c(minimo-20,maximo+250),ylim=c(minimoy-50,maximoy+50),xlab="Oeste - Leste",ylab="Sul - Norte",col=gray(100:1/100)) legend.krige(x.leg=c(maximo+100,maximo+150), y.leg=c(minimoy+10,maximoy-10),aux$predict, vert=TRUE,col=gray(100:1/100)) lines(sat$reg1) # ph aux$predict <- krig_bcrm$pred2 image(aux,xlim=c(minimo-20,maximo+250),ylim=c(minimoy-50,maximoy+50),xlab="Oeste - Leste",ylab="Sul - Norte",col=gray(100:1/100)) legend.krige(x.leg=c(maximo+100,maximo+150), y.leg=c(minimoy+10,maximoy-10),aux$predict, vert=TRUE,col=gray(100:1/100)) lines(ph$reg1)