Nothing
#' @name methods_spsur
#' @title Methods for class spsur
#' @description The \code{anova()} function provides tables of fitted
#' spsur models including information criteria (AIC and BIC),
#' log-likelihood and degrees of freedom of each fitted model. The
#' argument \code{lrtest} allows to perform LR tests between nested models.
#' The \code{plot()} function allows the user to plot both beta and spatial
#' coefficients for all equations of the spsur model. The argument
#' \code{viewplot} is used to choose between interactive or non-interactive
#' plots. The \code{print()} function is used to print short tables including the values of beta and
#' spatial coefficients as well as p-values of significance test for each
#' coefficient. This can be used as an alternative to
#' \code{\link{summary.spsur}} when a brief output is needed.
#' The rest of methods works in the usual way.
#'
#' @param object a \code{spsur} object created by \code{\link{spsurml}},
#' \code{\link{spsur3sls}} or \code{\link{spsurtime}}.
#' @param x similar to \code{object} argument for \code{print()}
#' and \code{plot} functions.
#' @param digits number of digits to show in printed tables.
#' Default: max(3L, getOption("digits") - 3L).
#' @param lrtest logical value to compute likelihood ratio
#' test for nested models in `anova` method. Default = \code{TRUE}
#' @param ci confidence level for the intervals in `plot` method.
#' Default \code{ci = 0.95}
#' @param viewplot logical value to show interactively the plots.
#' Default = \code{TRUE}
#' @param ... further arguments passed to or from other methods.
#' @examples
#' rm(list = ls()) # Clean memory
#' data(spc)
#' Tformula <- WAGE83 | WAGE81 ~ UN83 + NMR83 + SMSA | UN80 + NMR80 + SMSA
#' spcsur.sim <-spsurml(formula = Tformula, data = spc, type = "sim")
#' ## Print Table
#' print(spcsur.sim)
#' \donttest{
#' spcsur.slm <-spsurml(formula = Tformula, data = spc, type = "slm",
#' listw = Wspc)
#' # ANOVA table and LR test for nested models:
#' anova(spcsur.sim, spcsur.slm)
#' ## Plot spatial and beta coefficients
#' # Interactive plot
#' plot(spcsur.slm)
#' # Non-interactive plot
#' if (require(gridExtra)) {
#' pl <- plot(spcsur.slm, viewplot = FALSE)
#' grid.arrange(pl$lplbetas[[1]], pl$lplbetas[[2]],
#' pl$pldeltas, nrow = 3)
#' }
#'}
#' @author
#' \tabular{ll}{
#' Fernando Lopez \tab \email{fernando.lopez@@upct.es} \cr
#' Roman Minguez \tab \email{roman.minguez@@uclm.es} \cr
#' Jesus Mur \tab \email{jmur@@unizar.es} \cr
#' }
NULL
#' @name anova
#' @rdname methods_spsur
#' @export
anova.spsur <- function(object, ..., lrtest = TRUE) {
object <- list(object, ...)
ancall <- sys.call()
nmodels <- length(object)
vtypes <- character(length = nmodels)
vLL <- vdf <- vector(mode = "numeric", length=nmodels)
vAIC <- vBIC <- vector(mode = "numeric", length=nmodels)
if (lrtest) {
vlrtest <- vpval <- vector(mode = "numeric", length=nmodels)
vlrtest[1] <- vpval[1] <- NA
}
for (i in 1:nmodels) {
obj_i <- object[[i]]
if (!inherits(obj_i, "spsur"))
stop("object not a fitted spsur model")
if (is.null(obj_i$LL))
stop("object not fitted by maximum likelihood")
LL <- logLik(obj_i)
vAIC[i] <- stats::AIC(LL)
vBIC[i] <- stats::BIC(LL)
vLL[i] <- LL
vdf[i] <- attr(LL, "df")
vtypes[i] <- paste("model ",i,": ",obj_i$type, sep="")
if (lrtest) {
if(i>1) {
vlrtest[i] <- 2*(vLL[i] - vLL[i-1])
vpval[i] <- pchisq(vlrtest[i],
df = vdf[i] - vdf[i-1],
lower.tail = FALSE)
}
}
}
if (lrtest) {
res <- data.frame(logLik = vLL,
df = vdf,
AIC = vAIC,
BIC = vBIC,
LRtest = vlrtest,
`p-val` = vpval,
row.names = vtypes)
} else {
res <- data.frame(`logLik` = vLL,
df = vdf,
AIC = vAIC,
BIC = vBIC,
row.names = vtypes)
}
class(res) <- c("anova", "data.frame")
return(res)
}
#' @name coef
#' @rdname methods_spsur
#' @export
coef.spsur <- function(object, ...)
{
ret <- NULL
if (!(object$type == "sim" || object$type == "slx"))
ret <- c(ret, object$deltas)
ret <- c(ret, object$coefficients)
ret
}
#' @name fitted
#' @rdname methods_spsur
#' @export
fitted.spsur <- function(object, ...)
{
if (is.null(object$na.action))
fits <- object$fitted.values
else fits <- napredict(object$na.action, object$fitted.values)
G <- object$G
if (!(length(fits) %% G == 0))
stop("Length of fitted values incompatible with number of equations")
N <- length(fits) %/% G
res <- vector("list", G)
indx <- 1
for (i in 1:G) {
res[[i]] <- fits[indx:(indx + N - 1)]
indx <- indx + N
}
names(res) <- paste0("fitted_values_eq",1:G)
res
}
#' @name logLik
#' @rdname methods_spsur
#' @export
logLik.spsur <- function(object, ...)
{
LL <- c(object$LL)
class(LL) <- "logLik"
N <- length(residuals(object))
attr(LL, "nall") <- N
attr(LL, "nobs") <- N
attr(LL, "df") <- object$parameters
LL
}
# #' @name print
# #' @rdname methods_spsur
# #' @export
# print.spsur <- function(x, ...)
# {
# if (!(x$type == "sim" || x$type == "slx"
# || x$type == "sarar")) {
# G <- x$G
# for (i in 1:G) {
# delta <- x$deltas[i]
# if (isTRUE(all.equal(unname(delta), x$interval[1]))
# || isTRUE(all.equal(unname(delta), x$interval[2]))) {
# cat("Warning in equation ",i,"\n")
# warning("spatial parameter on interval bound -
# results should not be used")
# }
# }
# }
# cat("\nCall:\n")
# print(x$call)
# cat("Type:", x$type, "\n")
# cat("\nCoefficients:\n")
# print(coef(x))
# cat("\nLog likelihood:", logLik(x), "\n")
# invisible(x)
# }
# #' @name linearHypothesis
# #' @rdname methods_spsur
# #' @importFrom car linearHypothesis linearHypothesis.default
# #' @method linearHypothesis spsur
# #' @examples
# #' ## Linear Hypothesis (Wald test)
# #' linearHypothesis(spcsur.slm, c("rho_1 = rho_2"))
# #' ## Linear Hypothesis in matrix form (first the spatial coefficients)
# #' coef(spcsur.slm)
# #' R1 <- matrix(c(1,-1,0,0,0,0,0,0,0,0), nrow=1)
# #' b1 <- c(0)
# #' linearHypothesis(spcsur.slm, hypothesis.matrix = R1, rhs = b1)
# #' @export
# linearHypothesis.spsur <- function(object, ...) {
# car::linearHypothesis.default(object, ..., test = "Chisq")
# }
# #' @name lrtest
# #' @rdname methods_spsur
# #' @importFrom lmtest lrtest lrtest.default
# #' @method lrtest spsur
# #' @export
# lrtest.spsur <- function(object, ...) {
# typespsur <- function(x) {
# if (!is.null(x$type)) {
# return(x$type)
# } else stop ("object must be of spsur class")
# }
# lmtest::lrtest.default(object, ...,
# name = typespsur)
# }
#' @name residuals
#' @rdname methods_spsur
#' @export
residuals.spsur <- function(object, ...)
{
if (is.null(object$na.action))
resids <- object$residuals
else resids <- napredict(object$na.action, object$residuals)
G <- object$G
if (!(length(resids) %% G == 0))
stop("Length of residuals incompatible with number of equations")
N <- length(resids) %/% G
res <- vector("list", G)
indx <- 1
for (i in 1:G) {
res[[i]] <- resids[indx:(indx + N - 1)]
indx <- indx + N
}
names(res) <- paste0("residuals_eq",1:G)
res
}
#' @name vcov
#' @rdname methods_spsur
#' @export
vcov.spsur <- function(object, ...)
{
res <- object$resvar
# For compatibility with linearHypothesis method
idxsigma <- grepl("sigma", rownames(res))
res <- res[!idxsigma, !idxsigma]
res
}
#' @name print
#' @rdname methods_spsur
#' @export
print.spsur <- function(x, digits = max(3L, getOption("digits") - 3L),
...) {
if (!inherits(x, "spsur")) stop("Argument must be a spsur object")
summx <- summary(x)
fulltable <- summx$coef_table
G <- length(fulltable)
fullrownames <- NULL
for (i in 1:G) {
rownames_i <- rownames(fulltable[[i]])
rownames_i <- gsub(paste("_",i,sep=""), "",
rownames_i)
newrownames <- which (!(rownames_i %in% fullrownames))
rownames_i <- rownames_i[newrownames]
fullrownames <- c(fullrownames, rownames_i)
}
mtable <- matrix(NA, nrow = length(fullrownames),
ncol = 2*G)
rownames(mtable) <- fullrownames
colnames(mtable) <- paste(rep(c("coeff_","pval_"), G),
rep(1:G, each = 2), sep="")
idxcol <- 1
for (i in 1:G) {
table_i <- fulltable[[i]][,c("Estimate", "Pr(>|t|)")]
table_i[, 2] <- round(table_i[, 2],
digits = min(digits, 3))
rownames(table_i) <- gsub(paste("_",i,sep=""), "",
rownames(table_i))
mtable[rownames(table_i), c(idxcol,idxcol+1)] <- table_i
idxcol <- idxcol + 2
}
if (any(grepl("rho",rownames(mtable)))) {
idxrowrho <- which(grepl("rho",
rownames(mtable)))
rowrho <- matrix(mtable[idxrowrho, ], nrow=1)
rownames(rowrho) <- c("rho")
mtable <- mtable[-idxrowrho, ]
mtable <- rbind(mtable, rowrho)
}
if (any(grepl("lambda",rownames(mtable)))) {
idxrowlambda <- which(grepl("lambda",
rownames(mtable)))
rowlambda <- matrix(mtable[idxrowlambda, ], nrow=1)
rownames(rowlambda) <- c("lambda")
mtable <- mtable[-idxrowlambda, ]
mtable <- rbind(mtable, rowlambda)
}
print(round(mtable, digits = digits))
invisible(x)
}
#' @name plot
#' @rdname methods_spsur
#' @export
plot.spsur <- function(x, ci = 0.95, viewplot = TRUE, ...) {
if (!inherits(x, "spsur")) stop("Argument must be a spsur object")
critval <- qnorm((1-ci)/2, lower.tail = FALSE)
G <- x$G
p <- x$p
eq <- NULL
for (i in 1:G) {
eq <- c(eq, rep(i,p[i]))
}
eq <- as.factor(eq)
dfx <- data.frame(eq = eq,
nbetas = names(x$coefficients),
betas = x$coefficients,
sebetas = x$rest.se,
lwbetas = x$coefficients - critval*x$rest.se,
upbetas = x$coefficients + critval*x$rest.se,
row.names = NULL)
lplbetas <- vector("list", G)
for (i in 1:G) {
plbetasi <- ggplot2::ggplot(dfx[dfx$eq==i, ],
ggplot2::aes(.data$nbetas,
.data$betas,
colour = i)) +
ggplot2::geom_pointrange(ggplot2::aes(ymin = .data$lwbetas,
ymax = .data$upbetas),
show.legend = FALSE) +
ggplot2::labs(title = paste("Equation ",i),
x = "", y = "betas")
lplbetas[[i]] <- plbetasi
if (viewplot) {
gridExtra::grid.arrange(lplbetas[[i]], newpage = TRUE)
readline(prompt="Press [enter] to continue")
}
}
if(!is.null(x$deltas)) {
eq <- rep(0, length(x$deltas))
for (i in 1:G) {
ni <- i*as.integer(grepl(paste("_", i, sep=""), names(x$deltas)))
eq <- eq + ni
}
eq <- as.factor(eq)
dfxsp <- data.frame(eq = eq,
ndeltas = names(x$deltas),
deltas = x$deltas,
sedeltas = x$deltas.se,
lwdeltas = x$deltas - critval*x$deltas.se,
updeltas = x$deltas + critval*x$deltas.se,
row.names = NULL)
pldeltas <- ggplot2::ggplot(dfxsp, ggplot2::aes(.data$ndeltas,
.data$deltas,
colour = eq)) +
ggplot2::geom_pointrange(ggplot2::aes(ymin = .data$lwdeltas,
ymax = .data$updeltas),
show.legend = FALSE) +
ggplot2::labs(title = "", x = "",
y = "Spatial Coefficients")
if (viewplot) gridExtra::grid.arrange(pldeltas, nrow = 1)
} else pldeltas <- NULL
plots <- list(lplbetas = lplbetas,
pldeltas = pldeltas)
}
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.