###################################################### #####Sistematico - Arranjo Circular - trat fator require(geoR) a <- read.csv ("Coords.csv") #mu<-0.43 n <- 360 nrep <- n/10 sigma <- 0.013 phi <- 1 nugget <- 0.001 nsim <- 1000 #D<-as.matrix(read.csv("contraste.csv")) #D #verdc<-D%*%rbind(0.16,0.22,0.24,0.33,0.39,0.45,0.53,0.59,0.68,0.69) #verdef<-rbind(0.16,0.22,0.24,0.33,0.39,0.45,0.53,0.59,0.68,0.69) set.seed(1974) #### Funções Gerais #### Gerando os dados geoestatiticos tratamentos <- rep(1:10,each=nrep) efeitos <- function(tratamentos){ trat <- tratamentos vetor<- numeric(length(tratamentos)) vetor[trat== 1] <- 0.16 vetor[trat== 2] <- 0.22 vetor[trat== 3] <- 0.24 vetor[trat== 4] <- 0.33 vetor[trat== 5] <- 0.39 vetor[trat== 6] <- 0.45 vetor[trat== 7] <- 0.53 vetor[trat== 8] <- 0.59 vetor[trat== 9] <- 0.68 vetor[trat==10] <- 0.69 grid<-matrix(c(a$x,a$y),360,2) S <- grf(n,grid=grid,cov.model="matern",kappa=1.5,cov.pars=c(sigma,phi),message=FALSE) Z <- rnorm(n,0,nugget) data <- vetor + S$data + Z resposta <- data.frame(cbind(S$coords,data,trat)) resposta$trat <- factor(resposta$trat) return(resposta) } ## Erros quadráticos - Anova eq.an <- function(X){ v <- ((sum((an$residuals)^2))/(n-10)) var <- v*(solve(crossprod(X))) e.qua <- sqrt(diag(var)) return(e.qua) } ## Erros quadráticos - MQG eq.var <- function(X,inv){ var <- solve(crossprod(X,inv)%*%X) e.qua1 <- sqrt(diag(var)) return(e.qua1) } ##Função MQG - efeitos eftratmqg <- function(X,inv,y){ efmqg <- solve(crossprod(X,inv)%*%X)%*%crossprod(X,inv)%*%y return(efmqg) } #### Objetos que guardam os resultados gerados em cada simulação ## ANOVA efeitos efanova <- matrix (0,10,nsim) dimnames(efanova) <- list(paste("t",1:10,sep=""),NULL) ## ANova Intervalo de COnfianca eq.anova <- matrix(0,10,nsim) dimnames(eq.anova) <- list(paste("t",1:10,sep=""),NULL) ic1.anova <- NULL ic2.anova <- NULL ic3.anova <- NULL ic4.anova <- NULL ic5.anova <- NULL ic6.anova <- NULL ic7.anova <- NULL ic8.anova <- NULL ic9.anova <- NULL ic10.anova <- NULL #MV - efeitos efmv <- matrix (0,10,nsim) dimnames(efmv) <- list(paste("t",1:10,sep=""),NULL) ## MV Intervalo de Confianca eq.mv <- matrix(0,10,nsim) dimnames(eq.mv) <- list(paste("t",1:10,sep=""),NULL) ic1.mv <- NULL ic2.mv <- NULL ic3.mv <- NULL ic4.mv <- NULL ic5.mv <- NULL ic6.mv <- NULL ic7.mv <- NULL ic8.mv <- NULL ic9.mv <- NULL ic10.mv <- NULL #MQG - efeitos efmqg <- matrix (0,10,nsim) dimnames(efmqg) <- list(paste("t",1:10,sep=""),NULL) # MQG Intervalo de COnfianca eq.variog <- matrix(0,10,nsim) dimnames(eq.variog) <- list(paste("t",1:10,sep=""),NULL) ic1.var <- NULL ic2.var <- NULL ic3.var <- NULL ic4.var <- NULL ic5.var <- NULL ic6.var <- NULL ic7.var <- NULL ic8.var <- NULL ic9.var <- NULL ic10.var <- NULL #### cta<-0 repeat{ cta<-cta + 1 cat(paste("simulaçao ", cta, "de ", nsim, "\n")) dados<- efeitos(tratamentos) dados.geo <- as.geodata(dados,covar.col=4) y <- (dados$data) trat <- factor(dados$trat) an <- lm(data~(trat-1),data=dados) res <- residuals(an) m <- cbind(dados.geo$coords,res) resigeo <- as.geodata(m,coords.col=1:2,data.col=3) X <- model.matrix(an) ### Anova - efeitos efanova1 <- an$coefficients efanova[,cta] <- efanova1 # Erro quadrático - ANOVA erro.q1 <- eq.an(X) eq.anova[,cta] <- erro.q1 # Intervalo de Confiança - Anova ic1.ano <- efanova1[1]+ qt(c(0.025,0.975),(n-10))*erro.q1[1] ic1.anova <- rbind(ic1.anova,ic1.ano) ic2.ano <- efanova1[2]+ qt(c(0.025,0.975),(n-10))*erro.q1[2] ic2.anova <- rbind(ic2.anova,ic2.ano) ic3.ano <- efanova1[3]+ qt(c(0.025,0.975),(n-10))*erro.q1[3] ic3.anova <- rbind(ic3.anova,ic3.ano) ic4.ano <- efanova1[4]+ qt(c(0.025,0.975),(n-10))*erro.q1[4] ic4.anova <- rbind(ic4.anova,ic4.ano) ic5.ano <- efanova1[5]+ qt(c(0.025,0.975),(n-10))*erro.q1[5] ic5.anova <- rbind(ic5.anova,ic5.ano) ic6.ano <- efanova1[6]+ qt(c(0.025,0.975),(n-10))*erro.q1[6] ic6.anova <- rbind(ic6.anova,ic6.ano) ic7.ano <- efanova1[7]+ qt(c(0.025,0.975),(n-10))*erro.q1[7] ic7.anova <- rbind(ic7.anova,ic7.ano) ic8.ano <- efanova1[8]+ qt(c(0.025,0.975),(n-10))*erro.q1[8] ic8.anova <- rbind(ic8.anova,ic8.ano) ic9.ano <- efanova1[9]+ qt(c(0.025,0.975),(n-10))*erro.q1[9] ic9.anova <- rbind(ic9.anova,ic9.ano) ic10.ano <- efanova1[10]+ qt(c(0.025,0.975),(n-10))*erro.q1[10] ic10.anova <- rbind(ic10.anova,ic10.ano) ##### Construção do Variograma dist.eu <- max(dist(dados.geo$coords)) varemp <- variog(resigeo,coords=resigeo$coords,resigeo$data,max.dist=(dist.eu)*0.75,message=FALSE) varteo <- variofit(varemp,ini.cov.pars=c(0.01,10),nugget=0.01,cov.model="matern", kappa=1.5,limits=pars.limits(phi=c(0,dist.eu)),messages=FALSE) #### MV mv<-likfit(dados.geo,ini.cov.pars=cbind(c(0.005, 0.010, 0.015),c(1,3,6)),nugget=c(0.0005, 0.001, 0.005),trend=~(trat-1),cov.model="matern", kappa=1.5,messages=FALSE) ## Efeitos dos tratamentos para MV eftrat<-mv$beta efmv[,cta] <- eftrat # Erros quadráticos MV erro.q3 <- sqrt(diag(mv$beta.var)) eq.mv[,cta] <- erro.q3 # Intervalo de Confianca - MV ic1mv <- eftrat[1] + qnorm(c(0.025,0.975))*erro.q3[1] ic1.mv <- rbind(ic1.mv,ic1mv) ic2mv <- eftrat[2] + qnorm(c(0.025,0.975))*erro.q3[2] ic2.mv <- rbind(ic2.mv,ic2mv) ic3mv <- eftrat[3] + qnorm(c(0.025,0.975))*erro.q3[3] ic3.mv <- rbind(ic3.mv,ic3mv) ic4mv <- eftrat[4] + qnorm(c(0.025,0.975))*erro.q3[4] ic4.mv <- rbind(ic4.mv,ic4mv) ic5mv <- eftrat[5] + qnorm(c(0.025,0.975))*erro.q3[5] ic5.mv <- rbind(ic5.mv,ic5mv) ic6mv <- eftrat[6] + qnorm(c(0.025,0.975))*erro.q3[6] ic6.mv <- rbind(ic6.mv,ic6mv) ic7mv <- eftrat[7] + qnorm(c(0.025,0.975))*erro.q3[7] ic7.mv <- rbind(ic7.mv,ic7mv) ic8mv <- eftrat[8] + qnorm(c(0.025,0.975))*erro.q3[8] ic8.mv <- rbind(ic8.mv,ic8mv) ic9mv <- eftrat[9] + qnorm(c(0.025,0.975))*erro.q3[9] ic9.mv <- rbind(ic9.mv,ic9mv) ic10mv <- eftrat[10] + qnorm(c(0.025,0.975))*erro.q3[10] ic10.mv <- rbind(ic10.mv,ic10mv) #### MQG ## Construção da inversa da matriz de covariância inv<- varcov.spatial(m[,1:2],cov.pars=varteo$cov.pars, inv=TRUE, nug=varteo$nugget) ## Efeitos de tratamentos para MQG efmqg1<- eftratmqg(X,inv$inverse,y) efmqg[,cta] <- efmqg1 # Erros quadráticos do MQG erro.q2 <- eq.var(X, inv$inverse) eq.variog[,cta] <- erro.q2 # Intervalo de Confiança para MQG ic1.variog <- efmqg1[1] + qnorm(c(0.025,0.975))*erro.q2[1] ic1.var <- rbind(ic1.var,ic1.variog) ic2.variog <- efmqg1[2] + qnorm(c(0.025,0.975))*erro.q2[2] ic2.var <- rbind(ic2.var,ic2.variog) ic3.variog <- efmqg1[3] + qnorm(c(0.025,0.975))*erro.q2[3] ic3.var <- rbind(ic3.var,ic3.variog) ic4.variog <- efmqg1[4] + qnorm(c(0.025,0.975))*erro.q2[4] ic4.var <- rbind(ic4.var,ic4.variog) ic5.variog <- efmqg1[5] + qnorm(c(0.025,0.975))*erro.q2[5] ic5.var <- rbind(ic5.var,ic5.variog) ic6.variog <- efmqg1[6] + qnorm(c(0.025,0.975))*erro.q2[6] ic6.var <- rbind(ic6.var,ic6.variog) ic7.variog <- efmqg1[7] + qnorm(c(0.025,0.975))*erro.q2[7] ic7.var <- rbind(ic7.var,ic7.variog) ic8.variog <- efmqg1[8] + qnorm(c(0.025,0.975))*erro.q2[8] ic8.var <- rbind(ic8.var,ic8.variog) ic9.variog <- efmqg1[9] + qnorm(c(0.025,0.975))*erro.q2[9] ic9.var <- rbind(ic9.var,ic9.variog) ic10.variog <- efmqg1[10] + qnorm(c(0.025,0.975))*erro.q2[10] ic10.var <- rbind(ic10.var,ic10.variog) if(cta==nsim) break } ## Médias medanova <- apply (efanova, 1, mean) medmv <- apply (efmv, 1, mean) medmqg <- apply (efmqg, 1, mean) Manova <- matrix(medanova,10,1) Mmv <- matrix(medmv,10,1) Mmqg <- matrix(medmqg,10,1) Sb1 <- data.frame (Manova, Mmv, Mmqg) Sb1 write.table(Sb1, file = "sb1.csv", sep = ",", col.names = NA) ## Intervalo de Cobertura int.cob <- function(ic,ef){ vetor <- numeric(nsim) vetor[ic[,1]< ef & ef