> lagsarlm
function (formula, data = list(), listw, na.action, Durbin, type,
method = "eigen", quiet = NULL, zero.policy = NULL,
interval = NULL, tol.solve = .Machine$double.eps, trs = NULL,
control = list())
{
timings <- list()
.ptime_start <- proc.time()
con <- list(tol.opt = .Machine$double.eps^0.5, fdHess = NULL,
optimHess = FALSE, optimHessMethod = "optimHess",
compiled_sse = FALSE, Imult = 2, cheb_q = 5, MC_p = 16L,
MC_m = 30L, super = NULL, spamPivot = "MMD", in_coef = 0.1,
type = "MC", correct = TRUE, trunc = TRUE, SE_method = "LU",
nrho = 200, interpn = 2000, small_asy = TRUE, small = 1500,
SElndet = NULL, LU_order = FALSE, pre_eig = NULL, OrdVsign = 1)
nmsC <- names(con)
con[(namc <- names(control))] <- control
if (length(noNms <- namc[!namc %in% nmsC]))
warning("unknown names in control: ", paste(noNms,
collapse = ", "))
if (is.null(quiet))
quiet <- !get("verbose", envir = .spatialregOptions)
stopifnot(is.logical(quiet))
if (is.null(zero.policy))
zero.policy <- get.ZeroPolicyOption()
stopifnot(is.logical(zero.policy))
if (!inherits(formula, "formula"))
formula <- as.formula(formula)
mt <- terms(formula, data = data)
mf <- lm(formula, data, na.action = na.action, method = "model.frame")
na.act <- attr(mf, "na.action")
if (!inherits(listw, "listw"))
stop("No neighbourhood list")
can.sim <- FALSE
if (listw$style %in% c("W", "S"))
can.sim <- can.be.simmed(listw)
if (!is.null(na.act)) {
subset <- !(1:length(listw$neighbours) %in% na.act)
listw <- subset(listw, subset, zero.policy = zero.policy)
if (!is.null(con$pre_eig)) {
warning("NAs found, precomputed eigenvalues ignored")
con$pre_eig <- NULL
}
}
if (missing(type))
type <- "lag"
if (type == "Durbin")
type <- "mixed"
if (missing(Durbin))
Durbin <- ifelse(type == "lag", FALSE, TRUE)
if (listw$style != "W" && is.formula(Durbin)) {
Durbin <- TRUE
warning("formula Durbin requires row-standardised weights; set TRUE")
}
if (is.logical(Durbin) && isTRUE(Durbin))
type <- "mixed"
if (is.formula(Durbin))
type <- "mixed"
if (is.logical(Durbin) && !isTRUE(Durbin))
type <- "lag"
switch(type, lag = if (!quiet) cat("\nSpatial lag model\n"),
mixed = if (!quiet) cat("\nSpatial mixed autoregressive model\n"),
stop("\nUnknown model type\n"))
y <- model.extract(mf, "response")
x <- model.matrix(mt, mf)
if (NROW(x) != length(listw$neighbours))
stop("Input data and weights have different dimensions")
n <- NROW(x)
m <- NCOL(x)
stopifnot(is.logical(con$small_asy))
if (method != "eigen") {
if (con$small >= n && con$small_asy)
do_asy <- TRUE
else do_asy <- FALSE
}
else do_asy <- TRUE
if (is.null(con$fdHess)) {
con$fdHess <- method != "eigen" && !do_asy
fdHess <- NULL
}
stopifnot(is.logical(con$fdHess))
stopifnot(is.numeric(con$OrdVsign))
stopifnot(length(con$OrdVsign) == 1)
stopifnot(abs(con$OrdVsign) == 1)
xcolnames <- colnames(x)
K <- ifelse(xcolnames[1] == "(Intercept)", 2, 1)
wy <- lag.listw(listw, y, zero.policy = zero.policy)
if (any(is.na(wy)))
stop("NAs in lagged dependent variable")
dvars <- c(NCOL(x), 0L)
if (is.formula(Durbin) || isTRUE(Durbin)) {
prefix <- "lag"
if (isTRUE(Durbin)) {
WX <- create_WX(x, listw, zero.policy = zero.policy,
prefix = prefix)
}
else {
data1 <- data
if (!is.null(na.act) && (inherits(na.act, "omit") ||
inherits(na.act, "exclude"))) {
data1 <- data1[-c(na.act), ]
}
dmf <- lm(Durbin, data1, na.action = na.fail, method = "model.frame")
fx <- try(model.matrix(Durbin, dmf), silent = TRUE)
if (inherits(fx, "try-error"))
stop("Durbin variable mis-match")
WX <- create_WX(fx, listw, zero.policy = zero.policy,
prefix = prefix)
inds <- match(substring(colnames(WX), 5, nchar(colnames(WX))),
colnames(x))
if (anyNA(inds))
stop("WX variables not in X: ", paste(substring(colnames(WX),
5, nchar(colnames(WX)))[is.na(inds)], collapse = " "))
icept <- grep("(Intercept)", colnames(x))
iicept <- length(icept) > 0L
if (iicept) {
xn <- colnames(x)[-1]
}
else {
xn <- colnames(x)
}
wxn <- substring(colnames(WX), nchar(prefix) + 2,
nchar(colnames(WX)))
zero_fill <- NULL
if (length((which(!(xn %in% wxn)))) > 0L)
zero_fill <- length(xn) + (which(!(xn %in% wxn)))
}
dvars <- c(NCOL(x), NCOL(WX))
if (is.formula(Durbin)) {
attr(dvars, "f") <- Durbin
attr(dvars, "inds") <- inds
attr(dvars, "zero_fill") <- zero_fill
}
x <- cbind(x, WX)
m <- NCOL(x)
rm(WX)
}
lm.base <- lm(y ~ x - 1)
aliased <- is.na(coefficients(lm.base))
cn <- names(aliased)
names(aliased) <- substr(cn, 2, nchar(cn))
if (any(aliased)) {
if (is.formula(Durbin)) {
stop("Aliased variables found: ", paste(names(aliased)[aliased],
collapse = " "))
}
else {
warning("Aliased variables found: ", paste(names(aliased)[aliased],
collapse = " "))
nacoef <- which(aliased)
x <- x[, -nacoef]
}
}
LL_null_lm <- logLik(lm(y ~ 1))
m <- NCOL(x)
similar <- FALSE
lm.null <- lm(y ~ x - 1)
logLik_lm.model <- logLik(lm.null)
AIC_lm.model <- AIC(lm.null)
lm.w <- lm.fit(x, wy)
e.null <- lm.null$residuals
e.w <- lm.w$residuals
e.a <- t(e.null) %*% e.null
e.b <- t(e.w) %*% e.null
e.c <- t(e.w) %*% e.w
env <- new.env()
assign("y", y, envir = env)
assign("wy", wy, envir = env)
assign("x", x, envir = env)
assign("n", n, envir = env)
assign("m", m, envir = env)
assign("K", K, envir = env)
assign("e.a", e.a, envir = env)
assign("e.b", e.b, envir = env)
assign("e.c", e.c, envir = env)
assign("family", "SAR", envir = env)
assign("verbose", !quiet, envir = env)
assign("compiled_sse", con$compiled_sse, envir = env)
assign("can.sim", can.sim, envir = env)
assign("listw", listw, envir = env)
assign("similar", FALSE, envir = env)
assign("f_calls", 0L, envir = env)
assign("hf_calls", 0L, envir = env)
timings[["set_up"]] <- proc.time() - .ptime_start
.ptime_start <- proc.time()
if (!quiet)
cat("Jacobian calculated using ")
interval <- jacobianSetup(method, env, con, pre_eig = con$pre_eig,
trs = trs, interval = interval)
assign("interval", interval, envir = env)
nm <- paste(method, "set_up", sep = "_")
timings[[nm]] <- proc.time() - .ptime_start
.ptime_start <- proc.time()
opt <- optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE,
tol = con$tol.opt, env = env)
rho <- opt$maximum
if (isTRUE(all.equal(rho, interval[1])) || isTRUE(all.equal(rho,
interval[2])))
warning("rho on interval bound - results should not be used")
names(rho) <- "rho"
LL <- opt$objective
optres <- opt
nm <- paste(method, "opt", sep = "_")
timings[[nm]] <- proc.time() - .pti