Nothing
#' @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))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.