R/aovGeo.R
In lrcastro/spANOVA: Spatial Analysis of Field Trials Experiments using Geostatistics and Spatial Autoregressive Model
Defines functions
aovGeo
aovGeo.spVariofitCRD
print.GEOcrd
summary.GEOcrd
anova.GEOcrd
aovGeo.spVariofitRCBD
print.GEOrcbd
summary.GEOrcbd
anova.GEOrcbd
Documented in
aovGeo
aovGeo.spVariofitRCBD
#' @name aovGeo
#'
#' @title Analysis of variance using a geostatistical model to handle spatial dependence
#'
#' @description Fit an analysis of variance model using geostatistics for modeling the
#' spatial dependence among the observations.
#'
#' @usage aovGeo(model, cutoff = 0.5, tol = 1e-3)
#'
#' @param model an object of class spVariofit.
#' @param cutoff a value in (0,1] interval representing the percentage of the maximum
#' distance adopted to estimate the variogram in the algorithm suggested by
#' Pontes & Oliveira (2004). See 'Details'.
#' @param tol the desired accuracy.
#'
#'
#' @details
#' Three assumptions are made about the error in the analysis of variance (ANOVA):
#'
#' 1. The errors come from a normal distribution.
#'
#' 2. The errors have the same variance.
#'
#' 3. The errors are uncorrelated.
#'
#' However, in many experiments, especially field trials, there is a type of correlation
#' generated by the sample locations known as spatial autocorrelation, and this condition
#' violates the independence assumption.
#'
#' In that way, we need to regard this spatial autocorrelation and include it in the final
#' model. It could be done adopting a geostatistical model to characterize the spatial
#' variability among the observations directly in the covariance matrix.
#'
#' The geostatistical modeling is based on the residuals of the standard model
#' (where the errors are assumed to be independent, uncorrelated and having a normal
#' distribution with mean 0 and constant variance). The basic idea is using them to
#' estimate the residuals of the spatially autocorrelated model in order to fit a
#' theoretical geostatistic model to build the covariance matrix.
#' As Pointed by Pontes & Oliveira (2004), this task can be done using the
#' following algorithm
#'
#' 1 - Extract the residuals from the standard model
#'
#' 2 - Fit a variogram based on residuals obtained in step 1.
#'
#' 3 - Fit a theoretical model to describe the spatial dependence observed in the variogram.
#'
#' 4 - On basis in the theoretical model fitted in step 3 and its parameter estimates,
#' create the covariance matrix.
#'
#' 5 - Estimate the residuals using the covariance matrix obtained in step 4 and use
#' them to create a variogram.
#'
#' 6 - Fit a theoretical model to the residual variogram obtained in step 5 and use
#' its parameters estimates to build a new covariance matrix.
#'
#' 7 - Repeat 5 to 6 until convergence.
#'
#' Step 1 is implemented in spVariog. Steps 2 and 3 are implemented in spVariofit.
#' aovGeo implements steps 4 to 7 and needs a cutoff argument to define the maximum
#' distance in the computation of the residual variogram described in step 6
#'
#' In presence of spatial trend, the model is modified
#' in order to include the effect of the spatial coordinates.
#'
#' @return
#' \code{aovGeo} returns an object of \code{\link[base]{class}} "GEOanova".
#' The functions \code{\link[base]{summary}} and anova are used to obtain and print a summary
#' and analysis of variance table of the results.
#' An object of class "GEOanova" is a list containing the following components:
#'
#' \item{DF}{degrees of freedom of coordinates (in presence of spatial trend), treatments, block (if RCBD), residual and total.}
#' \item{SS}{sum of squares of coordinates (in presence of spatial trend), treatments, block (if RCBD), residual and total.}
#' \item{MS}{mean of squares of coordinates (in presence of spatial trend), treatments, block (if RCBD), residual and total.}
#' \item{Fc}{F statistic calculated for coordinates (in presence of spatial trend), treatment and block (if RCBD).}
#' \item{p.value}{p-value associated to F statistic for coordinates (in presence of spatial trend), treatment and block (if RCBD).}
#' \item{residuals}{residuals extracted from the model}
#' \item{params}{variogram parameter estimates from Pontes & Oliveira (2004) algorithm.}
#' \item{type}{type of trend in the data. It can be "trend" or "cte".}
#' \item{model}{geostatistical model used to describe the spatial dependence.}
#' \item{data}{data set used in the analysis.}
#' \item{des.mat}{design matrix.}
#' \item{beta}{parameter estimates of the linear model taking into account the spatial dependence.}
#' \item{n}{number of observations.}
#' \item{vcov}{covariance matrix built taking into account the spatial dependence.}
#' \item{design}{character string indicating the name of the experimental design.}
#'
#' @references
#' Nogueira, C. H., de Liima, R. R., & de Oliveira, M. S. (2013).
#' Aprimoramento da Análise de Variância: A Influência da Proximidade Espacial.
#' Rev. Bras. Biom, 31(3), 408-422.
#'
#' Nogueira, C.H., et al. (2015). Modelagem espacial na análise de um plantio
#' experimental de candeia. Rev. Bras. Biom., São Paulo, v.33, n.1, p.14-29.
#'
#' Pontes, J. M., & de Oliveira, M. S. (2004).
#' An alternative proposal to the analysis of field experiments using geostatistics.
#' Ciência e Agrotecnologia, 28(1), 135-141.
#'
#' Gotway, C. A., & Cressie, N. A. (1990). A spatial analysis of variance applied
#' to soil water infiltration. Water Resources Research, 26(11), 2695-703.
#'
#' @examples
#' data("crd_simulated")
#'
#' #Geodata object
#' geodados <- as.geodata(crd_simulated, coords.col = 1:2, data.col = 3,
#' covar.col = 4)
#' h_max <- summary(geodados)[[3]][[2]]
#' dist <- 0.6*h_max
#'
#' # Computing the variogram
#' variograma <- spVariog(geodata = geodados,
#' trend = "cte", max.dist = dist, design = "crd",
#' scale = FALSE)
#'
#' plot(variograma, ylab = "Semivariance", xlab = "Distance")
#'
#' # Gaussian Model
#' ols <- spVariofit(variograma, cov.model = "gaussian", weights = "equal",
#' max.dist = dist)
#'
#' lines(ols, col = 1)
#'
#' # Compute the model and get the analysis of variance table
#' mod <- aovGeo(ols, cutoff = 0.6)
#' anova(mod)
#'
#' @importFrom MASS ginv
#' @importFrom Matrix rankMatrix
#' @importFrom graphics lines plot
#' @importFrom stats AIC anova aov as.formula dist model.matrix pchisq qtukey pf
#' @import geoR
#' @export
aovGeo <- function(model, cutoff = 0.5, tol = 1e-3){
if(cutoff <= 0 | cutoff > 1) stop("'cutoff' should be in (0,1] interval")
if(is.numeric(cutoff) == FALSE) stop("'cutoff' must be numeric")
if(length(cutoff) > 1) stop("'cutoff' must of length 1")
if(cutoff < 0 | cutoff > 1) stop("'cutoff' must be in (0,1) interval")
UseMethod("aovGeo", model)
}
#' @export
#' @method aovGeo spVariofitCRD
aovGeo.spVariofitCRD <- function(model, cutoff = 0.5, tol = 1e-3){
# variaveis
dados <- model$data.geo
resp <- dados$data
treat <- dados$covariate[ ,1]
coords <- dados$coords
# Atributos do modelo geoestatistico
covMod <- model$mod$cov.model
nugget <- model$mod$nugget
psill <- model$mod$cov.pars[1]
phi <- model$mod$cov.pars[2]
trend <- model$trend
#matriz de covariancia espacial
sigma <- varcov.spatial(coords = coords, cov.model = covMod, nugget = nugget,
kappa = 0.5, cov.pars = c(psill, phi))$varcov
#Matriz de covariancias do modelo - V
V <- sigma*(1/as.numeric(nugget+psill))
i.V <- chol2inv(chol(V))
#Construindo a matriz de delineamento X
trt <- length(unique(treat)) #numero de tratamentos
r <- tapply(resp, treat, length) # O nº de obs por categoria
Y <- resp # variavel resposta
n <- sum(r)
# Obtendo o theta com dependencia espacial
# 1 col da Matriz de delineamento
if(trend == "cte"){
X1 <- as.matrix(model$des.mat[ ,1], ncol = 1)
X2 <- model$des.mat[ ,-1]
X <- cbind(X1, X2)
theta_dep <- ginv(t(X)%*%i.V%*%X)%*%t(X)%*%i.V%*%Y
}else{
X1<-matrix(rep(1,n))
X2 <- model$des.mat[ ,1:trt]
X <- cbind(X1, X2, coords)
theta_dep <- ginv(t(X)%*%i.V%*%X)%*%t(X)%*%i.V%*%Y
}
# Erros considerando a dependencia espacial
erro <- Y-X%*%theta_dep
# algoritmo iterativo pontes (2002)
e1 = e2 = e3 = 1
int <- 0
while((abs(e1) > tol) & (abs(e2) > tol) & (abs(e3) > tol)){
#Converter os dados para a classe geodata
geodados <- as.geodata(data.frame(coords, erro), coords.col = 1:2, data.col = 3)
#Definir a h maxima
h_max <- summary(geodados)[[3]][[2]]
dist <- cutoff*h_max
#Obter o semiovariograma empirico
vario.res <- variog(geodados, max.dist = dist, trend = "cte", messages = FALSE)
#plot(vario.res)
#Ajustando o modelo teorico via MQO
ols <- variofit(vario.res, cov.model = covMod, max.dist = dist,
messages = FALSE, ini.cov.pars = c(psill, phi), nugget = nugget)
#lines(ols)
#Extrair os parametros do objeto ols
phi1 <- summary(ols)[[10]][[3]] #Alcance
psill1 <- summary(ols)[[10]][[2]] #contribuição
nugget1 <- summary(ols)[[10]][[1]] #Efeito pepita
#Obter a matriz de covariancias
#matriz de covariancia espacial
sigma<-varcov.spatial(coords = coords, cov.model = covMod, nugget = nugget1,
kappa = 0.5, cov.pars = c(psill1, phi1))$varcov
#Matriz de covariancias do modelo V
V <- sigma*(1/as.numeric(nugget1+psill1))
i.V <- solve(V)
#Theta considerando a dependencia espacial
theta_dep <- ginv(t(X)%*%i.V%*%X)%*%t(X)%*%i.V%*%Y
#Erros considerando a dependencia espacial
erro <- Y-X%*%theta_dep
e1 <- phi-phi1
e2 <- psill-psill1
e3 <- nugget-nugget1
phi <- phi1
psill <- psill1
nugget <- nugget1
int <- int + 1
message(paste("Iteration number: ", int), "\n")
}
#pb <- txtProgressBar(style = 3)
# Modelo sem tendencia
if(trend == "cte"){
#Forma como apresentado na tese
P <- i.V - i.V%*%X%*%ginv(t(X)%*%i.V%*%X)%*%t(X)%*%i.V
#setTxtProgressBar(pb, 1/10)
P1 <- i.V - i.V%*%X1%*%solve(t(X1)%*%i.V%*%X1)%*%t(X1)%*%i.V
#setTxtProgressBar(pb, 2/10)
gltrat <- rankMatrix(P1 - P)[1]
#setTxtProgressBar(pb, 3/10)
glres <- rankMatrix(P)[1]
#setTxtProgressBar(pb, 4/10)
gltot <- rankMatrix(P1)[1]
#setTxtProgressBar(pb, 5/10)
sqtrat <- t(Y)%*%(P1 - P)%*%Y
Sys.sleep(0.1)
#setTxtProgressBar(pb, 6/10)
sqres <- t(Y)%*%(P)%*%Y
#setTxtProgressBar(pb, 7/10)
sqtot <- t(Y)%*%(P1)%*%Y
#setTxtProgressBar(pb, 8/10)
QMTrat <- sqtrat/gltrat
QMRes <- sqres/glres
F.trat <- QMTrat/QMRes
pvalTrat <- round(pf(F.trat, gltrat, glres, lower.tail = F), 4)
#setTxtProgressBar(pb, 9/10)
} else {
# Modelo com tendencia
Mc <- cbind(X1, X2, coords)
M2 <- cbind(X1, X2)
M3 <- cbind(X1, coords)
Pc <- i.V%*%Mc%*%ginv(t(Mc)%*%i.V%*%Mc)%*%t(Mc)%*%i.V
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 1/10)
P2 <- i.V%*%M2%*%ginv(t(M2)%*%i.V%*%M2)%*%t(M2)%*%i.V
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 2/10)
P3 <- i.V%*%M3%*%solve(t(M3)%*%i.V%*%M3)%*%t(M3)%*%i.V
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 3/10)
sqcoord <- t(Y)%*%(Pc - P2)%*%Y
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 4/10)
sqtrat <- t(Y)%*%(Pc - P3)%*%Y
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 5/10)
sqres <- t(Y)%*%(i.V - Pc)%*%Y
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 6/10)
sqtot <- sqcoord + sqtrat + sqres
####Graus de liberdade###
glres <- rankMatrix(i.V - Pc)[1]
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 7/10)
gltrat <- rankMatrix(Pc - P3)[1]
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 8/10)
glcoord <- rankMatrix(Pc - P2)[1]
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 9/10)
gltot <- glres + gltrat + glcoord
### Construindo a analise de variancia ###
QMTrat <- sqtrat/gltrat
QMRes <- sqres/glres
QMCoord <- sqcoord/glcoord
F.trat <- QMTrat/QMRes
F.coord <- QMCoord/QMRes
pvalTrat <- pf(F.trat, gltrat, glres, lower.tail = FALSE)
pvalCoord <- pf(F.coord, glcoord, glres, lower.tail = FALSE)
}
#Estimacao das componentes do modelo
trend1 <- X%*%theta_dep
erro <- Y-trend1
sig0 <- varcov.spatial(coords, cov.model = covMod,
nugget = 0, cov.pars = c(psill, phi))$varcov
i.sig <- chol2inv(chol(sigma))
compEsp <- sig0%*%i.sig%*%erro
resid <- as.numeric(erro - compEsp)
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 10/10)
# Lista que armazenara a saida
if(trend != "cte"){
output <- list(DF = c(glcoord, gltrat, glres, gltot),
SS = c(sqcoord, sqtrat, sqres, sqtot),
MS = c(QMCoord, QMTrat, QMRes),
Fc = c(F.coord, F.trat),
p.value = c(pvalCoord, pvalTrat),
residuals = resid, params = c(sigsq = psill, phi = phi, tausq = nugget),
type = "trend", model = covMod, data = dados, des.mat = model$des.mat,
beta = theta_dep, n = n, vcov = sigma, design = "CRD")
} else {
output <- list(DF = c(gltrat, glres, gltot),
SS = c(sqtrat, sqres, sqtot),
MS = c(QMTrat, QMRes),
Fc = c(F.trat),
p.value = c(pvalTrat),
residuals = resid, params = c(sigsq = psill, phi = phi, tausq = nugget),
type = "cte", model = covMod, data = dados, des.mat = model$des.mat,
beta = theta_dep, n = n, vcov = sigma, design = "crd")
}
class(output) <- c("GEOanova", "GEOcrd")
return(output)
}
#' @export
#' @method print GEOcrd
print.GEOcrd <- function(x, ...){
if(x$type == "cte"){
cat("Terms:","\n")
trm <- data.frame(treat = c(as.character(round(x$SS[1],3)),as.character(x$DF[1])),
Residuals = c(as.character(round(x$SS[2],3)),as.character(x$DF[2])))
rownames(trm) <- c("Sum of Squares","Deg. of Freedom")
print(trm)
rse <- sqrt(x$MS[2])
cat("\n")
cat("Residual standard error:", rse, "\n")
cat("Covariance Model:", x$model, "\n")
} else {
cat("Terms:","\n")
trm <- data.frame(coordinates = c(as.character(round(x$SS[1],3)),as.character(x$DF[1])),
treat = c(as.character(round(x$SS[2],3)),as.character(x$DF[2])),
Residuals = c(as.character(round(x$SS[3],3)),as.character(x$DF[3]))
)
rownames(trm) <- c("Sum of Squares","Deg. of Freedom")
print(trm)
rse <- sqrt(x$MS[3])
cat("\n")
cat("Residual standard error:",rse, "\n")
cat("Covariance Model:", x$model, "\n")
}
}
#' @export
#' @method summary GEOcrd
summary.GEOcrd <- function(object, ...){
k <- nrow(object$beta)
if(object$type == "cte"){
cat("Terms:","\n")
trm <- data.frame(treat = c(as.character(round(object$SS[1],3)),as.character(object$DF[1])),
Residuals = c(as.character(round(object$SS[2],3)),as.character(object$DF[2])))
rownames(trm) <- c("Sum of Squares","Deg. of Freedom")
print(trm)
rse <- sqrt(object$MS[2])
cat("\n")
cat("Residual standard error:",rse, "\n")
r2 <- 1 - (object$SS[2]/object$SS[3])
r2adj <- 1 - ((object$n - 1)/(object$n - k))*(1-r2)
cat("Adjusted R-squared:", r2adj, "\n")
cat("Covariance Model:", object$model, "\n")
cat("Spatial Parameters:", "\n")
print(object$params)
} else {
cat("Terms:","\n")
trm <- data.frame(coordinates = c(as.character(round(object$SS[1],3)),as.character(object$DF[1])),
treat = c(as.character(round(object$SS[2],3)),as.character(object$DF[2])),
Residuals = c(as.character(round(object$SS[3],3)),as.character(object$DF[3]))
)
rownames(trm) <- c("Sum of Squares","Deg. of Freedom")
print(trm)
rse <- sqrt(object$MS[3])
cat("\n")
cat("Residual standard error:",rse, "\n")
r2 <- 1 - (object$SS[3]/object$SS[4])
r2adj <- 1 - ((object$n - 1)/(object$n - k))*(1-r2)
cat("Adjusted R-squared:", r2adj, "\n")
cat("Covariance Model:", object$model, "\n")
cat("Spatial Parameters:", "\n")
print(object$params)
}
}
#' @export
#' @method anova GEOcrd
anova.GEOcrd <- function(object, compare = FALSE, ...) {
if(is.logical(compare) == FALSE){
warning("'compare' must be logical. Assuming compare == FALSE")
compare = FALSE
}
if(object$type == "trend"){
cat("Analysis of Variance With Spatially Correlated Errors","\n")
cat("\n")
star1 <- stars.pval(object$p.value[1])
star2 <- stars.pval(object$p.value[2])
anova.p1 <- data.frame("DF" = round(object$DF[1:4],0),
"SS" = round(object$SS[1:4],4),
"MS" = c(round(object$MS[1:3],4), ""),
"Fc" = c(round(object$Fc[1:2],4), "", ""),
"Pv" = c(format.pval(object$p.value[1]), format.pval(object$p.value[2]), "", ""),
"St" = c(star1, star2, "", "")
)
colnames(anova.p1) <- c("Df", "Sum Sq", "Mean Sq", "F value", "Pr(>F)", "")
rownames(anova.p1) <- c("Coordinates", "Treatment", "Residuals", "Total")
print(anova.p1)
cat("---", "\n")
cat("Signif. codes: ", attr(star1, "legend"), "\n")
if(compare == TRUE){
cat("\n", "\n")
cat("---------------------------------------------------------------","\n")
cat("Standard Analysis of Variance", "\n")
cat("---------------------------------------------------------------")
cat("\n")
Treatment <- factor(object$data$covariate[ ,1])
Resp <- object$data$data
print(anova(aov(Resp ~ Treatment)))
}
} else {
cat("Analysis of Variance With Spatially Correlated Errors","\n")
cat("\n")
star <- stars.pval(object$p.value)
anova.p1 <- data.frame("DF" = round(object$DF[1:3],0),
"SS" = round(object$SS[1:3],4),
"MS" = c(round(object$MS[1:2],4), ""),
"Fc" = c(round(object$Fc,4), "", ""),
"Pv" = c(format.pval(object$p.value), "", ""),
"St" = c(star, "", "")
)
colnames(anova.p1) <- c("Df", "Sum Sq", "Mean Sq", "F value" ,"Pr(>F)", "")
rownames(anova.p1) <- c("Treatment", "Residuals", "Total")
print(anova.p1)
cat("---", "\n")
cat("Signif. codes: ", attr(star, "legend"),"\n")
if(compare == TRUE){
cat("\n", "\n")
cat("---------------------------------------------------------------","\n")
cat("Standard Analysis of Variance", "\n")
cat("---------------------------------------------------------------")
cat("\n")
Treatment <- factor(object$data$covariate[ ,1])
Resp <- object$data$data
print(anova(aov(Resp ~ Treatment)))
}
}
return(invisible(anova.p1))
}
#' @export
#' @method aovGeo spVariofitRCBD
#' @rdname aovGeo
aovGeo.spVariofitRCBD <- function(model, cutoff = 0.5, tol = 1e-3){
# variaveis
dados <- model$data.geo
resp <- dados$data
coords <- dados$coords
block <- dados$covariate[,2]
treat <- dados$covariate[,1]
# Atributos do modelo geoestatistico
covMod <- model$mod$cov.model
nugget <- model$mod$nugget
psill <- model$mod$cov.pars[1]
phi <- model$mod$cov.pars[2]
trend <- "cte"
#matriz de covariancia espacial
sigma <- varcov.spatial(coords = coords, cov.model = covMod, nugget = nugget,
kappa = 0.5, cov.pars = c(psill, phi))$varcov
trt <- length(unique(treat)) #numero de tratamentos
blk <- length(unique(block)) #numero de blocos
r <- tapply(resp, treat, length) #
n <- sum(r)
#Matriz de covariancias do modelo - V
V <- sigma*(1/as.numeric(nugget+psill))
i.V <- chol2inv(chol(V))
#Construindo a matriz de delineamento X
Y <- resp # variavel resposta
X1 <- as.matrix(model$des.mat[ ,1], ncol = 1)
X2 <- model$des.mat[ ,-c(1,(trt+2):(blk+trt+1))]
X <- model$des.mat
theta_dep <- ginv(t(X)%*%i.V%*%X)%*%t(X)%*%i.V%*%Y
#Erros considerando a dependência espacial
erro <- Y-X%*%theta_dep
# algoritmo iterativo pontes (2002)
e1 = e2 = e3 = 1
int <- 0
while((abs(e1) > tol) & (abs(e2) > tol) & (abs(e3) > tol)){
#Converter os dados para a classe geodata
geodados <- as.geodata(data.frame(coords, erro), coords.col = 1:2, data.col = 3)
#Definir a h maxima
h_max <- summary(geodados)[[3]][[2]]
dist <- cutoff*h_max
#Obter o semiovariograma empirico
vario.res <- variog(geodados, max.dist = dist, trend = "cte", messages = FALSE)
#plot(vario.res)
#Ajustando o modelo teorico via MQO
ols <- variofit(vario.res, cov.model = covMod, max.dist = dist,
messages = FALSE, ini.cov.pars = c(psill, phi), nugget = nugget)
#lines(ols)
#Extrair os parametros do objeto ols
phi1 <- summary(ols)[[10]][[3]] #Alcance
psill1 <- summary(ols)[[10]][[2]] #contribuição
nugget1 <- summary(ols)[[10]][[1]] #Efeito pepita
#Obter a matriz de covariancias
#matriz de covariancia espacial
sigma<-varcov.spatial(coords = coords, cov.model = covMod, nugget = nugget1,
kappa = 0.5, cov.pars = c(psill1, phi1))$varcov
#Matriz de covariancias do modelo V
V <- sigma*(1/as.numeric(nugget1+psill1))
i.V <- solve(V)
#Theta considerando a dependencia espacial
theta_dep <- ginv(t(X)%*%i.V%*%X)%*%t(X)%*%i.V%*%Y
theta_dep <- ginv(t(X)%*%i.V%*%X)%*%t(X)%*%i.V%*%Y
#Erros considerando a dependencia espacial
erro <- Y-X%*%theta_dep
e1 <- phi-phi1
e2 <- psill-psill1
e3 <- nugget-nugget1
phi <- phi1
psill <- psill1
nugget <- nugget1
int <- int + 1
message(paste("Iteration number: ", int), "\n")
}
#pb <- txtProgressBar(style = 3)
#Obtendo a soma de quadrados
P1 <- i.V%*%X1%*%solve(t(X1)%*%i.V%*%X1)%*%t(X1)%*%i.V
P <- i.V%*%X%*%ginv(t(X)%*%i.V%*%X)%*%t(X)%*%i.V
#setTxtProgressBar(pb, 1/10)
R <- i.V-P1
P2 <- R%*%X2%*%ginv(t(X2)%*%R%*%X2)%*%t(X2)%*%R
#setTxtProgressBar(pb, 2/10)
SQT <- t(Y)%*%(i.V-P1)%*%Y
SQRes <- t(Y)%*%(i.V-P)%*%Y
#setTxtProgressBar(pb, 3/10)
SQTrat <- t(Y)%*%(P2)%*%Y
SQBloco <- t(Y)%*%(P-P1-P2)%*%Y
#setTxtProgressBar(pb, 4/10)
#Graus de liberdade
glt <- rankMatrix(i.V-P1)[1] #Total
#setTxtProgressBar(pb, 5/10)
gle <- rankMatrix(i.V-P)[1] #Erro
#setTxtProgressBar(pb, 6/10)
gltra <- rankMatrix(P2)[1] #tratamento
#setTxtProgressBar(pb, 7/10)
glblo <- rankMatrix(P-P1-P2)[1] #Bloco
#setTxtProgressBar(pb, 8/10)
#Quadrados médios
QMTrat <- SQTrat/gltra
QMRes <- SQRes/gle
QMBlo <- SQBloco/glblo
#setTxtProgressBar(pb, 9/10)
#Estatística F
F0 <- QMTrat/QMRes
F0b <- QMBlo/QMRes
pvalTrat <- round(pf(F0, gltra, gle, lower.tail = F), 4)
pvalBlo <- round(pf(F0b, glblo, gle, lower.tail = F), 4)
#Estimação das componentes do modelo
trend <- X%*%theta_dep
erro <- Y-trend
sig0 <- varcov.spatial(coords, cov.model = covMod,
nugget = 0, cov.pars = c(psill,phi))$varcov
i.sig <- chol2inv(chol(sigma))
compEsp <- sig0%*%i.sig%*%erro
resid <- as.numeric(erro-compEsp)
#setTxtProgressBar(pb, 10/10)
output <- list(DF = c(gltra, glblo, gle, glt),
SS = c(SQTrat, SQBloco, SQRes, SQT),
MS = c(QMTrat, QMBlo, QMRes),
Fc = c(F0, F0b),
p.value = c(pvalTrat, pvalBlo),
residuals = resid, params = c(sigsq = psill, phi = phi, tausq = nugget),
type = "cte", model = covMod, data = dados, des.mat = model$des.mat,
beta = theta_dep, n = n, vcov = sigma, design = "rcbd")
class(output) <- c("GEOanova", "GEOrcbd")
return(output)
}
#' @export
#' @method print GEOrcbd
print.GEOrcbd <- function(x, ...){
cat("Terms:","\n")
trm <- data.frame(treat = c(as.character(round(x$SS[1], 3)), as.character(x$DF[1])),
block = c(as.character(round(x$SS[2], 3)), as.character(x$DF[2])),
Residuals = c(as.character(round(x$SS[3], 3)), as.character(x$DF[3])))
rownames(trm) <- c("Sum of Squares","Deg. of Freedom")
print(trm)
rse <- sqrt(x$MS[3])
cat("\n")
cat("Residual standard error:", rse, "\n")
cat("Covariance Model:", x$model, "\n")
}
#' @export
#' @method summary GEOrcbd
summary.GEOrcbd<- function(object, ...){
cat("Terms:","\n")
trm <- data.frame(treat = c(as.character(round(object$SS[1], 3)), as.character(object$DF[1])),
block = c(as.character(round(object$SS[2], 3)), as.character(object$DF[2])),
Residuals = c(as.character(round(object$SS[3], 3)), as.character(object$DF[3])))
rownames(trm) <- c("Sum of Squares","Deg. of Freedom")
print(trm)
rse <- sqrt(object$MS[3])
cat("\n")
cat("Residual standard error:", rse, "\n")
k <- nrow(object$beta)
r2 <- 1 - (object$SS[3]/object$SS[4])
r2adj <- 1 - ((object$n - 1)/(object$n - k))*(1-r2)
cat("Adjusted R-squared:", r2adj, "\n")
cat("Covariance Model:", object$model, "\n")
cat("Spatial Parameters:", "\n")
print(object$params)
}
#' @export
#' @method anova GEOrcbd
anova.GEOrcbd <- function(object, compare = FALSE, ...) {
if(is.logical(compare) == FALSE){
warning("'compare' must be logical. Assuming compare == FALSE")
compare = FALSE
}
cat("Analysis of Variance With Spatially Correlated Errors","\n")
cat("\n")
star1 <- stars.pval(object$p.value[1])
star2 <- stars.pval(object$p.value[2])
anova.p1 <- data.frame("DF" = round(object$DF[1:4],0),
"SS" = round(object$SS[1:4],4),
"MS" = c(round(object$MS[1:3],4), ""),
"Fc" = c(round(object$Fc[1:2],4), "", ""),
"Pv" = c(format.pval(object$p.value[1]), format.pval(object$p.value[2]), "", ""),
"St" = c(star1, star2, "", "")
)
colnames(anova.p1) <- c("Df", "Sum Sq", "Mean Sq", "F value" ,"Pr(>F)", "")
rownames(anova.p1) <- c("Treatment","Block","Residuals","Total")
print(anova.p1)
cat("---","\n")
cat("Signif. codes: ",attr(star1, "legend"), "\n")
if(compare == TRUE){
cat("\n", "\n")
cat("---------------------------------------------------------------","\n")
cat("Standard Analysis of Variance", "\n")
cat("---------------------------------------------------------------")
cat("\n")
Treatment <- as.factor(object$data$covariate[ ,1])
Block <- as.factor(object$data$covariate[ ,2])
Resp <- object$data$data
print(anova(aov(Resp ~ Treatment + Block)))
}
return(invisible(anova.p1))
}
lrcastro/spANOVA documentation built on Aug. 21, 2019, 11:45 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions aovGeo aovGeo.spVariofitCRD print.GEOcrd summary.GEOcrd anova.GEOcrd aovGeo.spVariofitRCBD print.GEOrcbd summary.GEOrcbd anova.GEOrcbd
Documented in aovGeo aovGeo.spVariofitRCBD
#' @name aovGeo
#'
#' @title Analysis of variance using a geostatistical model to handle spatial dependence
#'
#' @description Fit an analysis of variance model using geostatistics for modeling the
#' spatial dependence among the observations.
#'
#' @usage aovGeo(model, cutoff = 0.5, tol = 1e-3)
#'
#' @param model an object of class spVariofit.
#' @param cutoff a value in (0,1] interval representing the percentage of the maximum
#' distance adopted to estimate the variogram in the algorithm suggested by
#' Pontes & Oliveira (2004). See 'Details'.
#' @param tol the desired accuracy.
#'
#'
#' @details
#' Three assumptions are made about the error in the analysis of variance (ANOVA):
#'
#' 1. The errors come from a normal distribution.
#'
#' 2. The errors have the same variance.
#'
#' 3. The errors are uncorrelated.
#'
#' However, in many experiments, especially field trials, there is a type of correlation
#' generated by the sample locations known as spatial autocorrelation, and this condition
#' violates the independence assumption.
#'
#' In that way, we need to regard this spatial autocorrelation and include it in the final
#' model. It could be done adopting a geostatistical model to characterize the spatial
#' variability among the observations directly in the covariance matrix.
#'
#' The geostatistical modeling is based on the residuals of the standard model
#' (where the errors are assumed to be independent, uncorrelated and having a normal
#' distribution with mean 0 and constant variance). The basic idea is using them to
#' estimate the residuals of the spatially autocorrelated model in order to fit a
#' theoretical geostatistic model to build the covariance matrix.
#' As Pointed by Pontes & Oliveira (2004), this task can be done using the
#' following algorithm
#'
#' 1 - Extract the residuals from the standard model
#'
#' 2 - Fit a variogram based on residuals obtained in step 1.
#'
#' 3 - Fit a theoretical model to describe the spatial dependence observed in the variogram.
#'
#' 4 - On basis in the theoretical model fitted in step 3 and its parameter estimates,
#' create the covariance matrix.
#'
#' 5 - Estimate the residuals using the covariance matrix obtained in step 4 and use
#' them to create a variogram.
#'
#' 6 - Fit a theoretical model to the residual variogram obtained in step 5 and use
#' its parameters estimates to build a new covariance matrix.
#'
#' 7 - Repeat 5 to 6 until convergence.
#'
#' Step 1 is implemented in spVariog. Steps 2 and 3 are implemented in spVariofit.
#' aovGeo implements steps 4 to 7 and needs a cutoff argument to define the maximum
#' distance in the computation of the residual variogram described in step 6
#'
#' In presence of spatial trend, the model is modified
#' in order to include the effect of the spatial coordinates.
#'
#' @return
#' \code{aovGeo} returns an object of \code{\link[base]{class}} "GEOanova".
#' The functions \code{\link[base]{summary}} and anova are used to obtain and print a summary
#' and analysis of variance table of the results.
#' An object of class "GEOanova" is a list containing the following components:
#'
#' \item{DF}{degrees of freedom of coordinates (in presence of spatial trend), treatments, block (if RCBD), residual and total.}
#' \item{SS}{sum of squares of coordinates (in presence of spatial trend), treatments, block (if RCBD), residual and total.}
#' \item{MS}{mean of squares of coordinates (in presence of spatial trend), treatments, block (if RCBD), residual and total.}
#' \item{Fc}{F statistic calculated for coordinates (in presence of spatial trend), treatment and block (if RCBD).}
#' \item{p.value}{p-value associated to F statistic for coordinates (in presence of spatial trend), treatment and block (if RCBD).}
#' \item{residuals}{residuals extracted from the model}
#' \item{params}{variogram parameter estimates from Pontes & Oliveira (2004) algorithm.}
#' \item{type}{type of trend in the data. It can be "trend" or "cte".}
#' \item{model}{geostatistical model used to describe the spatial dependence.}
#' \item{data}{data set used in the analysis.}
#' \item{des.mat}{design matrix.}
#' \item{beta}{parameter estimates of the linear model taking into account the spatial dependence.}
#' \item{n}{number of observations.}
#' \item{vcov}{covariance matrix built taking into account the spatial dependence.}
#' \item{design}{character string indicating the name of the experimental design.}
#'
#' @references
#' Nogueira, C. H., de Liima, R. R., & de Oliveira, M. S. (2013).
#' Aprimoramento da Análise de Variância: A Influência da Proximidade Espacial.
#' Rev. Bras. Biom, 31(3), 408-422.
#'
#' Nogueira, C.H., et al. (2015). Modelagem espacial na análise de um plantio
#' experimental de candeia. Rev. Bras. Biom., São Paulo, v.33, n.1, p.14-29.
#'
#' Pontes, J. M., & de Oliveira, M. S. (2004).
#' An alternative proposal to the analysis of field experiments using geostatistics.
#' Ciência e Agrotecnologia, 28(1), 135-141.
#'
#' Gotway, C. A., & Cressie, N. A. (1990). A spatial analysis of variance applied
#' to soil water infiltration. Water Resources Research, 26(11), 2695-703.
#'
#' @examples
#' data("crd_simulated")
#'
#' #Geodata object
#' geodados <- as.geodata(crd_simulated, coords.col = 1:2, data.col = 3,
#' covar.col = 4)
#' h_max <- summary(geodados)[[3]][[2]]
#' dist <- 0.6*h_max
#'
#' # Computing the variogram
#' variograma <- spVariog(geodata = geodados,
#' trend = "cte", max.dist = dist, design = "crd",
#' scale = FALSE)
#'
#' plot(variograma, ylab = "Semivariance", xlab = "Distance")
#'
#' # Gaussian Model
#' ols <- spVariofit(variograma, cov.model = "gaussian", weights = "equal",
#' max.dist = dist)
#'
#' lines(ols, col = 1)
#'
#' # Compute the model and get the analysis of variance table
#' mod <- aovGeo(ols, cutoff = 0.6)
#' anova(mod)
#'
#' @importFrom MASS ginv
#' @importFrom Matrix rankMatrix
#' @importFrom graphics lines plot
#' @importFrom stats AIC anova aov as.formula dist model.matrix pchisq qtukey pf
#' @import geoR
#' @export
aovGeo <- function(model, cutoff = 0.5, tol = 1e-3){
if(cutoff <= 0 | cutoff > 1) stop("'cutoff' should be in (0,1] interval")
if(is.numeric(cutoff) == FALSE) stop("'cutoff' must be numeric")
if(length(cutoff) > 1) stop("'cutoff' must of length 1")
if(cutoff < 0 | cutoff > 1) stop("'cutoff' must be in (0,1) interval")
UseMethod("aovGeo", model)
}
#' @export
#' @method aovGeo spVariofitCRD
aovGeo.spVariofitCRD <- function(model, cutoff = 0.5, tol = 1e-3){
# variaveis
dados <- model$data.geo
resp <- dados$data
treat <- dados$covariate[ ,1]
coords <- dados$coords
# Atributos do modelo geoestatistico
covMod <- model$mod$cov.model
nugget <- model$mod$nugget
psill <- model$mod$cov.pars[1]
phi <- model$mod$cov.pars[2]
trend <- model$trend
#matriz de covariancia espacial
sigma <- varcov.spatial(coords = coords, cov.model = covMod, nugget = nugget,
kappa = 0.5, cov.pars = c(psill, phi))$varcov
#Matriz de covariancias do modelo - V
V <- sigma*(1/as.numeric(nugget+psill))
i.V <- chol2inv(chol(V))
#Construindo a matriz de delineamento X
trt <- length(unique(treat)) #numero de tratamentos
r <- tapply(resp, treat, length) # O nº de obs por categoria
Y <- resp # variavel resposta
n <- sum(r)
# Obtendo o theta com dependencia espacial
# 1 col da Matriz de delineamento
if(trend == "cte"){
X1 <- as.matrix(model$des.mat[ ,1], ncol = 1)
X2 <- model$des.mat[ ,-1]
X <- cbind(X1, X2)
theta_dep <- ginv(t(X)%*%i.V%*%X)%*%t(X)%*%i.V%*%Y
}else{
X1<-matrix(rep(1,n))
X2 <- model$des.mat[ ,1:trt]
X <- cbind(X1, X2, coords)
theta_dep <- ginv(t(X)%*%i.V%*%X)%*%t(X)%*%i.V%*%Y
}
# Erros considerando a dependencia espacial
erro <- Y-X%*%theta_dep
# algoritmo iterativo pontes (2002)
e1 = e2 = e3 = 1
int <- 0
while((abs(e1) > tol) & (abs(e2) > tol) & (abs(e3) > tol)){
#Converter os dados para a classe geodata
geodados <- as.geodata(data.frame(coords, erro), coords.col = 1:2, data.col = 3)
#Definir a h maxima
h_max <- summary(geodados)[[3]][[2]]
dist <- cutoff*h_max
#Obter o semiovariograma empirico
vario.res <- variog(geodados, max.dist = dist, trend = "cte", messages = FALSE)
#plot(vario.res)
#Ajustando o modelo teorico via MQO
ols <- variofit(vario.res, cov.model = covMod, max.dist = dist,
messages = FALSE, ini.cov.pars = c(psill, phi), nugget = nugget)
#lines(ols)
#Extrair os parametros do objeto ols
phi1 <- summary(ols)[[10]][[3]] #Alcance
psill1 <- summary(ols)[[10]][[2]] #contribuição
nugget1 <- summary(ols)[[10]][[1]] #Efeito pepita
#Obter a matriz de covariancias
#matriz de covariancia espacial
sigma<-varcov.spatial(coords = coords, cov.model = covMod, nugget = nugget1,
kappa = 0.5, cov.pars = c(psill1, phi1))$varcov
#Matriz de covariancias do modelo V
V <- sigma*(1/as.numeric(nugget1+psill1))
i.V <- solve(V)
#Theta considerando a dependencia espacial
theta_dep <- ginv(t(X)%*%i.V%*%X)%*%t(X)%*%i.V%*%Y
#Erros considerando a dependencia espacial
erro <- Y-X%*%theta_dep
e1 <- phi-phi1
e2 <- psill-psill1
e3 <- nugget-nugget1
phi <- phi1
psill <- psill1
nugget <- nugget1
int <- int + 1
message(paste("Iteration number: ", int), "\n")
}
#pb <- txtProgressBar(style = 3)
# Modelo sem tendencia
if(trend == "cte"){
#Forma como apresentado na tese
P <- i.V - i.V%*%X%*%ginv(t(X)%*%i.V%*%X)%*%t(X)%*%i.V
#setTxtProgressBar(pb, 1/10)
P1 <- i.V - i.V%*%X1%*%solve(t(X1)%*%i.V%*%X1)%*%t(X1)%*%i.V
#setTxtProgressBar(pb, 2/10)
gltrat <- rankMatrix(P1 - P)[1]
#setTxtProgressBar(pb, 3/10)
glres <- rankMatrix(P)[1]
#setTxtProgressBar(pb, 4/10)
gltot <- rankMatrix(P1)[1]
#setTxtProgressBar(pb, 5/10)
sqtrat <- t(Y)%*%(P1 - P)%*%Y
Sys.sleep(0.1)
#setTxtProgressBar(pb, 6/10)
sqres <- t(Y)%*%(P)%*%Y
#setTxtProgressBar(pb, 7/10)
sqtot <- t(Y)%*%(P1)%*%Y
#setTxtProgressBar(pb, 8/10)
QMTrat <- sqtrat/gltrat
QMRes <- sqres/glres
F.trat <- QMTrat/QMRes
pvalTrat <- round(pf(F.trat, gltrat, glres, lower.tail = F), 4)
#setTxtProgressBar(pb, 9/10)
} else {
# Modelo com tendencia
Mc <- cbind(X1, X2, coords)
M2 <- cbind(X1, X2)
M3 <- cbind(X1, coords)
Pc <- i.V%*%Mc%*%ginv(t(Mc)%*%i.V%*%Mc)%*%t(Mc)%*%i.V
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 1/10)
P2 <- i.V%*%M2%*%ginv(t(M2)%*%i.V%*%M2)%*%t(M2)%*%i.V
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 2/10)
P3 <- i.V%*%M3%*%solve(t(M3)%*%i.V%*%M3)%*%t(M3)%*%i.V
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 3/10)
sqcoord <- t(Y)%*%(Pc - P2)%*%Y
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 4/10)
sqtrat <- t(Y)%*%(Pc - P3)%*%Y
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 5/10)
sqres <- t(Y)%*%(i.V - Pc)%*%Y
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 6/10)
sqtot <- sqcoord + sqtrat + sqres
####Graus de liberdade###
glres <- rankMatrix(i.V - Pc)[1]
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 7/10)
gltrat <- rankMatrix(Pc - P3)[1]
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 8/10)
glcoord <- rankMatrix(Pc - P2)[1]
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 9/10)
gltot <- glres + gltrat + glcoord
### Construindo a analise de variancia ###
QMTrat <- sqtrat/gltrat
QMRes <- sqres/glres
QMCoord <- sqcoord/glcoord
F.trat <- QMTrat/QMRes
F.coord <- QMCoord/QMRes
pvalTrat <- pf(F.trat, gltrat, glres, lower.tail = FALSE)
pvalCoord <- pf(F.coord, glcoord, glres, lower.tail = FALSE)
}
#Estimacao das componentes do modelo
trend1 <- X%*%theta_dep
erro <- Y-trend1
sig0 <- varcov.spatial(coords, cov.model = covMod,
nugget = 0, cov.pars = c(psill, phi))$varcov
i.sig <- chol2inv(chol(sigma))
compEsp <- sig0%*%i.sig%*%erro
resid <- as.numeric(erro - compEsp)
#Sys.sleep(0.1)
#setTxtProgressBar(pb, 10/10)
# Lista que armazenara a saida
if(trend != "cte"){
output <- list(DF = c(glcoord, gltrat, glres, gltot),
SS = c(sqcoord, sqtrat, sqres, sqtot),
MS = c(QMCoord, QMTrat, QMRes),
Fc = c(F.coord, F.trat),
p.value = c(pvalCoord, pvalTrat),
residuals = resid, params = c(sigsq = psill, phi = phi, tausq = nugget),
type = "trend", model = covMod, data = dados, des.mat = model$des.mat,
beta = theta_dep, n = n, vcov = sigma, design = "CRD")
} else {
output <- list(DF = c(gltrat, glres, gltot),
SS = c(sqtrat, sqres, sqtot),
MS = c(QMTrat, QMRes),
Fc = c(F.trat),
p.value = c(pvalTrat),
residuals = resid, params = c(sigsq = psill, phi = phi, tausq = nugget),
type = "cte", model = covMod, data = dados, des.mat = model$des.mat,
beta = theta_dep, n = n, vcov = sigma, design = "crd")
}
class(output) <- c("GEOanova", "GEOcrd")
return(output)
}
#' @export
#' @method print GEOcrd
print.GEOcrd <- function(x, ...){
if(x$type == "cte"){
cat("Terms:","\n")
trm <- data.frame(treat = c(as.character(round(x$SS[1],3)),as.character(x$DF[1])),
Residuals = c(as.character(round(x$SS[2],3)),as.character(x$DF[2])))
rownames(trm) <- c("Sum of Squares","Deg. of Freedom")
print(trm)
rse <- sqrt(x$MS[2])
cat("\n")
cat("Residual standard error:", rse, "\n")
cat("Covariance Model:", x$model, "\n")
} else {
cat("Terms:","\n")
trm <- data.frame(coordinates = c(as.character(round(x$SS[1],3)),as.character(x$DF[1])),
treat = c(as.character(round(x$SS[2],3)),as.character(x$DF[2])),
Residuals = c(as.character(round(x$SS[3],3)),as.character(x$DF[3]))
)
rownames(trm) <- c("Sum of Squares","Deg. of Freedom")
print(trm)
rse <- sqrt(x$MS[3])
cat("\n")
cat("Residual standard error:",rse, "\n")
cat("Covariance Model:", x$model, "\n")
}
}
#' @export
#' @method summary GEOcrd
summary.GEOcrd <- function(object, ...){
k <- nrow(object$beta)
if(object$type == "cte"){
cat("Terms:","\n")
trm <- data.frame(treat = c(as.character(round(object$SS[1],3)),as.character(object$DF[1])),
Residuals = c(as.character(round(object$SS[2],3)),as.character(object$DF[2])))
rownames(trm) <- c("Sum of Squares","Deg. of Freedom")
print(trm)
rse <- sqrt(object$MS[2])
cat("\n")
cat("Residual standard error:",rse, "\n")
r2 <- 1 - (object$SS[2]/object$SS[3])
r2adj <- 1 - ((object$n - 1)/(object$n - k))*(1-r2)
cat("Adjusted R-squared:", r2adj, "\n")
cat("Covariance Model:", object$model, "\n")
cat("Spatial Parameters:", "\n")
print(object$params)
} else {
cat("Terms:","\n")
trm <- data.frame(coordinates = c(as.character(round(object$SS[1],3)),as.character(object$DF[1])),
treat = c(as.character(round(object$SS[2],3)),as.character(object$DF[2])),
Residuals = c(as.character(round(object$SS[3],3)),as.character(object$DF[3]))
)
rownames(trm) <- c("Sum of Squares","Deg. of Freedom")
print(trm)
rse <- sqrt(object$MS[3])
cat("\n")
cat("Residual standard error:",rse, "\n")
r2 <- 1 - (object$SS[3]/object$SS[4])
r2adj <- 1 - ((object$n - 1)/(object$n - k))*(1-r2)
cat("Adjusted R-squared:", r2adj, "\n")
cat("Covariance Model:", object$model, "\n")
cat("Spatial Parameters:", "\n")
print(object$params)
}
}
#' @export
#' @method anova GEOcrd
anova.GEOcrd <- function(object, compare = FALSE, ...) {
if(is.logical(compare) == FALSE){
warning("'compare' must be logical. Assuming compare == FALSE")
compare = FALSE
}
if(object$type == "trend"){
cat("Analysis of Variance With Spatially Correlated Errors","\n")
cat("\n")
star1 <- stars.pval(object$p.value[1])
star2 <- stars.pval(object$p.value[2])
anova.p1 <- data.frame("DF" = round(object$DF[1:4],0),
"SS" = round(object$SS[1:4],4),
"MS" = c(round(object$MS[1:3],4), ""),
"Fc" = c(round(object$Fc[1:2],4), "", ""),
"Pv" = c(format.pval(object$p.value[1]), format.pval(object$p.value[2]), "", ""),
"St" = c(star1, star2, "", "")
)
colnames(anova.p1) <- c("Df", "Sum Sq", "Mean Sq", "F value", "Pr(>F)", "")
rownames(anova.p1) <- c("Coordinates", "Treatment", "Residuals", "Total")
print(anova.p1)
cat("---", "\n")
cat("Signif. codes: ", attr(star1, "legend"), "\n")
if(compare == TRUE){
cat("\n", "\n")
cat("---------------------------------------------------------------","\n")
cat("Standard Analysis of Variance", "\n")
cat("---------------------------------------------------------------")
cat("\n")
Treatment <- factor(object$data$covariate[ ,1])
Resp <- object$data$data
print(anova(aov(Resp ~ Treatment)))
}
} else {
cat("Analysis of Variance With Spatially Correlated Errors","\n")
cat("\n")
star <- stars.pval(object$p.value)
anova.p1 <- data.frame("DF" = round(object$DF[1:3],0),
"SS" = round(object$SS[1:3],4),
"MS" = c(round(object$MS[1:2],4), ""),
"Fc" = c(round(object$Fc,4), "", ""),
"Pv" = c(format.pval(object$p.value), "", ""),
"St" = c(star, "", "")
)
colnames(anova.p1) <- c("Df", "Sum Sq", "Mean Sq", "F value" ,"Pr(>F)", "")
rownames(anova.p1) <- c("Treatment", "Residuals", "Total")
print(anova.p1)
cat("---", "\n")
cat("Signif. codes: ", attr(star, "legend"),"\n")
if(compare == TRUE){
cat("\n", "\n")
cat("---------------------------------------------------------------","\n")
cat("Standard Analysis of Variance", "\n")
cat("---------------------------------------------------------------")
cat("\n")
Treatment <- factor(object$data$covariate[ ,1])
Resp <- object$data$data
print(anova(aov(Resp ~ Treatment)))
}
}
return(invisible(anova.p1))
}
#' @export
#' @method aovGeo spVariofitRCBD
#' @rdname aovGeo
aovGeo.spVariofitRCBD <- function(model, cutoff = 0.5, tol = 1e-3){
# variaveis
dados <- model$data.geo
resp <- dados$data
coords <- dados$coords
block <- dados$covariate[,2]
treat <- dados$covariate[,1]
# Atributos do modelo geoestatistico
covMod <- model$mod$cov.model
nugget <- model$mod$nugget
psill <- model$mod$cov.pars[1]
phi <- model$mod$cov.pars[2]
trend <- "cte"
#matriz de covariancia espacial
sigma <- varcov.spatial(coords = coords, cov.model = covMod, nugget = nugget,
kappa = 0.5, cov.pars = c(psill, phi))$varcov
trt <- length(unique(treat)) #numero de tratamentos
blk <- length(unique(block)) #numero de blocos
r <- tapply(resp, treat, length) #
n <- sum(r)
#Matriz de covariancias do modelo - V
V <- sigma*(1/as.numeric(nugget+psill))
i.V <- chol2inv(chol(V))
#Construindo a matriz de delineamento X
Y <- resp # variavel resposta
X1 <- as.matrix(model$des.mat[ ,1], ncol = 1)
X2 <- model$des.mat[ ,-c(1,(trt+2):(blk+trt+1))]
X <- model$des.mat
theta_dep <- ginv(t(X)%*%i.V%*%X)%*%t(X)%*%i.V%*%Y
#Erros considerando a dependência espacial
erro <- Y-X%*%theta_dep
# algoritmo iterativo pontes (2002)
e1 = e2 = e3 = 1
int <- 0
while((abs(e1) > tol) & (abs(e2) > tol) & (abs(e3) > tol)){
#Converter os dados para a classe geodata
geodados <- as.geodata(data.frame(coords, erro), coords.col = 1:2, data.col = 3)
#Definir a h maxima
h_max <- summary(geodados)[[3]][[2]]
dist <- cutoff*h_max
#Obter o semiovariograma empirico
vario.res <- variog(geodados, max.dist = dist, trend = "cte", messages = FALSE)
#plot(vario.res)
#Ajustando o modelo teorico via MQO
ols <- variofit(vario.res, cov.model = covMod, max.dist = dist,
messages = FALSE, ini.cov.pars = c(psill, phi), nugget = nugget)
#lines(ols)
#Extrair os parametros do objeto ols
phi1 <- summary(ols)[[10]][[3]] #Alcance
psill1 <- summary(ols)[[10]][[2]] #contribuição
nugget1 <- summary(ols)[[10]][[1]] #Efeito pepita
#Obter a matriz de covariancias
#matriz de covariancia espacial
sigma<-varcov.spatial(coords = coords, cov.model = covMod, nugget = nugget1,
kappa = 0.5, cov.pars = c(psill1, phi1))$varcov
#Matriz de covariancias do modelo V
V <- sigma*(1/as.numeric(nugget1+psill1))
i.V <- solve(V)
#Theta considerando a dependencia espacial
theta_dep <- ginv(t(X)%*%i.V%*%X)%*%t(X)%*%i.V%*%Y
theta_dep <- ginv(t(X)%*%i.V%*%X)%*%t(X)%*%i.V%*%Y
#Erros considerando a dependencia espacial
erro <- Y-X%*%theta_dep
e1 <- phi-phi1
e2 <- psill-psill1
e3 <- nugget-nugget1
phi <- phi1
psill <- psill1
nugget <- nugget1
int <- int + 1
message(paste("Iteration number: ", int), "\n")
}
#pb <- txtProgressBar(style = 3)
#Obtendo a soma de quadrados
P1 <- i.V%*%X1%*%solve(t(X1)%*%i.V%*%X1)%*%t(X1)%*%i.V
P <- i.V%*%X%*%ginv(t(X)%*%i.V%*%X)%*%t(X)%*%i.V
#setTxtProgressBar(pb, 1/10)
R <- i.V-P1
P2 <- R%*%X2%*%ginv(t(X2)%*%R%*%X2)%*%t(X2)%*%R
#setTxtProgressBar(pb, 2/10)
SQT <- t(Y)%*%(i.V-P1)%*%Y
SQRes <- t(Y)%*%(i.V-P)%*%Y
#setTxtProgressBar(pb, 3/10)
SQTrat <- t(Y)%*%(P2)%*%Y
SQBloco <- t(Y)%*%(P-P1-P2)%*%Y
#setTxtProgressBar(pb, 4/10)
#Graus de liberdade
glt <- rankMatrix(i.V-P1)[1] #Total
#setTxtProgressBar(pb, 5/10)
gle <- rankMatrix(i.V-P)[1] #Erro
#setTxtProgressBar(pb, 6/10)
gltra <- rankMatrix(P2)[1] #tratamento
#setTxtProgressBar(pb, 7/10)
glblo <- rankMatrix(P-P1-P2)[1] #Bloco
#setTxtProgressBar(pb, 8/10)
#Quadrados médios
QMTrat <- SQTrat/gltra
QMRes <- SQRes/gle
QMBlo <- SQBloco/glblo
#setTxtProgressBar(pb, 9/10)
#Estatística F
F0 <- QMTrat/QMRes
F0b <- QMBlo/QMRes
pvalTrat <- round(pf(F0, gltra, gle, lower.tail = F), 4)
pvalBlo <- round(pf(F0b, glblo, gle, lower.tail = F), 4)
#Estimação das componentes do modelo
trend <- X%*%theta_dep
erro <- Y-trend
sig0 <- varcov.spatial(coords, cov.model = covMod,
nugget = 0, cov.pars = c(psill,phi))$varcov
i.sig <- chol2inv(chol(sigma))
compEsp <- sig0%*%i.sig%*%erro
resid <- as.numeric(erro-compEsp)
#setTxtProgressBar(pb, 10/10)
output <- list(DF = c(gltra, glblo, gle, glt),
SS = c(SQTrat, SQBloco, SQRes, SQT),
MS = c(QMTrat, QMBlo, QMRes),
Fc = c(F0, F0b),
p.value = c(pvalTrat, pvalBlo),
residuals = resid, params = c(sigsq = psill, phi = phi, tausq = nugget),
type = "cte", model = covMod, data = dados, des.mat = model$des.mat,
beta = theta_dep, n = n, vcov = sigma, design = "rcbd")
class(output) <- c("GEOanova", "GEOrcbd")
return(output)
}
#' @export
#' @method print GEOrcbd
print.GEOrcbd <- function(x, ...){
cat("Terms:","\n")
trm <- data.frame(treat = c(as.character(round(x$SS[1], 3)), as.character(x$DF[1])),
block = c(as.character(round(x$SS[2], 3)), as.character(x$DF[2])),
Residuals = c(as.character(round(x$SS[3], 3)), as.character(x$DF[3])))
rownames(trm) <- c("Sum of Squares","Deg. of Freedom")
print(trm)
rse <- sqrt(x$MS[3])
cat("\n")
cat("Residual standard error:", rse, "\n")
cat("Covariance Model:", x$model, "\n")
}
#' @export
#' @method summary GEOrcbd
summary.GEOrcbd<- function(object, ...){
cat("Terms:","\n")
trm <- data.frame(treat = c(as.character(round(object$SS[1], 3)), as.character(object$DF[1])),
block = c(as.character(round(object$SS[2], 3)), as.character(object$DF[2])),
Residuals = c(as.character(round(object$SS[3], 3)), as.character(object$DF[3])))
rownames(trm) <- c("Sum of Squares","Deg. of Freedom")
print(trm)
rse <- sqrt(object$MS[3])
cat("\n")
cat("Residual standard error:", rse, "\n")
k <- nrow(object$beta)
r2 <- 1 - (object$SS[3]/object$SS[4])
r2adj <- 1 - ((object$n - 1)/(object$n - k))*(1-r2)
cat("Adjusted R-squared:", r2adj, "\n")
cat("Covariance Model:", object$model, "\n")
cat("Spatial Parameters:", "\n")
print(object$params)
}
#' @export
#' @method anova GEOrcbd
anova.GEOrcbd <- function(object, compare = FALSE, ...) {
if(is.logical(compare) == FALSE){
warning("'compare' must be logical. Assuming compare == FALSE")
compare = FALSE
}
cat("Analysis of Variance With Spatially Correlated Errors","\n")
cat("\n")
star1 <- stars.pval(object$p.value[1])
star2 <- stars.pval(object$p.value[2])
anova.p1 <- data.frame("DF" = round(object$DF[1:4],0),
"SS" = round(object$SS[1:4],4),
"MS" = c(round(object$MS[1:3],4), ""),
"Fc" = c(round(object$Fc[1:2],4), "", ""),
"Pv" = c(format.pval(object$p.value[1]), format.pval(object$p.value[2]), "", ""),
"St" = c(star1, star2, "", "")
)
colnames(anova.p1) <- c("Df", "Sum Sq", "Mean Sq", "F value" ,"Pr(>F)", "")
rownames(anova.p1) <- c("Treatment","Block","Residuals","Total")
print(anova.p1)
cat("---","\n")
cat("Signif. codes: ",attr(star1, "legend"), "\n")
if(compare == TRUE){
cat("\n", "\n")
cat("---------------------------------------------------------------","\n")
cat("Standard Analysis of Variance", "\n")
cat("---------------------------------------------------------------")
cat("\n")
Treatment <- as.factor(object$data$covariate[ ,1])
Block <- as.factor(object$data$covariate[ ,2])
Resp <- object$data$data
print(anova(aov(Resp ~ Treatment + Block)))
}
return(invisible(anova.p1))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.