# R/anova.elliptical.R In gwer: Geographically Weighted Elliptical Regression

#### Documented in anova.elliptical

#' @title Analysis of Deviance for Elliptical Regression Models
#' @method anova elliptical
#' @description Compute an analysis of deviance table for the fitted elliptical regression model.
#' @param object an object with the result of the fitted elliptical regression model.
#' @param dispersion the dispersion parameter for the fitting family. By default it is obtained from the object.
#' @param test a character string indicating the hypothesis test considered. Only chi-square test is implemented at moment.
#' @param ... arguments to be used to form the default control argument if it is not supplied directly.
#' @return Return an object of class \dQuote{anova}. This object contain the analysis of deviance table.
#' @references Cysneiros, F. J. A., Paula, G. A., and Galea, M. (2007). Heteroscedastic
#' symmetrical linear models. Statistics & probability letters, 77(11), 1084-1090.
#' \doi{10.1016/j.spl.2007.01.012}
#' @keywords Elliptical regression models
#' @keywords Hypothesis tests
#' @keywords ANOVA
#' @examples
#' data(luzdat)
#' y <- luzdat\$y
#' x1 <- luzdat\$x1 ; x1 <- factor(x1) ; x1 <- C(x1,treatment)
#' x2 <- luzdat\$x2
#' x3 <- (luzdat\$x2)^2
#' luz <- data.frame(y,x1,x2,x3)
#' elliptical.fitt <- elliptical(y ~ x1+x2+x3, family = Student(df=5),
#' data=luz)
#' anova(elliptical.fitt, test = "Chisq")
#' @export

anova.elliptical <- function(object, ..., dispersion = NULL, test = c("Chisq"))
{
test <- match.arg(test)
margs <- function(...) nargs()
if (margs(...))
return(anova.ellipticallist(list(object, ...), test = test))
Terms <- object\$terms
term.labels <- attr(Terms, "term.labels")
nt <- length(term.labels)
m <- model.frame(object)
family <- family(object)[[1]]
y <- model.extract(m, "response")
loglik <- double(nt + 1)
df.res <- loglik
if (nt) {
loglik[nt + 1] <- -2 * object\$loglik
df.res[nt + 1] <- object\$df.residual
fit <- object
for (iterm in seq(from = nt, to = 1, by = -1)) {
ww <- fit\$weights
argslist <- list(object = fit, formula = eval(parse(text = paste("~ . -",
term.labels[iterm]))))
fit <- do.call("update", argslist)
loglik[iterm] <- -2 * fit\$loglik
df.res[iterm] <- fit\$df.residual
}
dev <- c(NA, -diff(loglik))
df <- c(NA, -diff(df.res))
}
else {
loglik[1] <- -2 * object\$loglik
df.res[1] <- dim(y)[1] - attr(Terms, "intercept")
dev <- df <- as.numeric(NA)
}
heading <- c("Analysis of Deviance Table\n", paste("Error distribution:",
family), paste("Response:", as.character(formula(object))[2],
"\n", sep = ""), "Terms added sequentially (first to last)")
if (is.null(dispersion)) {
dispersion <- 1
df.dispersion <- if (dispersion == 1)
Inf
else object\$df.residual
}
else df.scale <- Inf
aod <- data.frame(Df = df, Deviance = dev, `Resid. Df` = df.res,
`-2*LL` = loglik, row.names = c("NULL", term.labels),
check.names = F)
class(aod) <- c("anova", "data.frame")
if (is.null(test))
return(aod)
else aod <- stat.anova(aod, test, scale = 1, df.scale = df.dispersion,
n = nrow(y))
}

anova.ellipticallist <- function (object, ...,test = c("Chisq"))
{
diff.term <- function(term.labels, i) {
t1 <- term.labels[[1]]
t2 <- term.labels[[2]]
m1 <- match(t1, t2, F)
m2 <- match(t2, t1, F)
if (all(m1)) {
if (all(m2))
return("=")
else return(paste(c("", t2[-m1]), collapse = "+"))
}
else {
if (all(m2))
return(paste(c("", t1[-m2]), collapse = "-"))
else return(paste(i - 1, i, sep = " vs. "))
}
}
test <- match.arg(test)
rt <- length(object)
if (rt == 1) {
object <- object[[1]]
UseMethod("anova")
}
forms <- sapply(object, function(x) as.character(formula(x)))
subs <- as.logical(match(forms[2, ], forms[2, 1], F))
if (!all(subs))
warning("Some fit objects deleted because response differs from the first model")
if (sum(subs) == 1)
stop("The first model has a different response from the rest")
forms <- forms[, subs]
object <- object[subs]
dfres <- sapply(object, "[[", "df.resid")
m2loglik <- -2 * sapply(object, "[[", "loglik")
tl <- lapply(object, labels)
rt <- length(m2loglik)
effects <- character(rt)
for (i in 2:rt) effects[i] <- diff.term(tl[c(i - 1, i)],
i)
dm2loglik <- -diff(m2loglik)
ddf <- -diff(dfres)
family <- family(object[[1]])[[1]]
fixed <- object[[1]]\$fixed
heading <- c("Analysis of Deviance Table\n", paste("Error distribution:",
family), paste("Response:", forms[2, 1], "\n"))
nmodels <- length(object)
#  varirest <- lapply(object, function(x) paste(deparse(x\$call\$restrict),
#                                               collapse = "\n"))
variables <- lapply(object, function(x) paste(deparse(formula(x)),
collapse = "\n"))
#  topnote <- paste("Model ", format(1:nmodels), ": ", variables,
#                   "     restrict =", varirest, sep = "", collapse = "\n")
aod <- data.frame(Terms = forms[3, ], `Resid. Df` = dfres,
`-2*LL` = m2loglik, Test = effects, Df = c(NA, abs(ddf)),
Deviance = c(NA, abs(dm2loglik)), check.names = F)
if (!is.null(test)) {
n <- length(object[[1]]\$residuals)
o <- order(dfres)
aod <- stat.anova(aod, test, 1, dfres[o[1]], n)
}
else aod