1# Copyright 2005-2012 by Roger Bivand 2spautolm <- function(formula, data = list(), listw, weights, 3 na.action, family="SAR", method="eigen", verbose=NULL, trs=NULL, 4 interval=NULL, zero.policy=NULL, tol.solve=.Machine$double.eps, llprof=NULL, 5 control=list()) { 6# .Deprecated("spatialreg::spautolm", msg="Function spautolm moved to the spatialreg package") 7# if (!requireNamespace("spatialreg", quietly=TRUE)) 8# stop("install the spatialreg package") 9# if (requireNamespace("spatialreg", quietly=TRUE)) { 10# if (!missing(weights)) stop("run spatialreg::spautolm directly") 11# return(spatialreg::spautolm(formula=formula, data=data, listw=listw, na.action=na.action, family = family, method=method, verbose=verbose, trs=trs, interval=interval, zero.policy=zero.policy, tol.solve=tol.solve, llprof=llprof, control=control)) 12# } 13 warning("install the spatialreg package") 14# if (FALSE) { 15 timings <- list() 16 .ptime_start <- proc.time() 17 con <- list(tol.opt=.Machine$double.eps^(2/3), 18 fdHess=NULL, optimHess=FALSE, optimHessMethod="optimHess", 19 Imult=2, cheb_q=5, MC_p=16, MC_m=30, super=NULL, spamPivot="MMD", 20 in_coef=0.1, type="MC", 21 correct=TRUE, trunc=TRUE, SE_method="LU", nrho=200, 22 interpn=2000, small_asy=TRUE, small=1500, SElndet=NULL, 23 LU_order=FALSE, pre_eig=NULL) 24 nmsC <- names(con) 25 con[(namc <- names(control))] <- control 26 if (length(noNms <- namc[!namc %in% nmsC])) 27 warning("unknown names in control: ", paste(noNms, collapse = ", ")) 28 if (!inherits(listw, "listw")) 29 stop("No neighbourhood list") 30 if (is.null(verbose)) verbose <- get("verbose", envir = .spdepOptions) 31 stopifnot(is.logical(verbose)) 32 if (is.null(zero.policy)) 33 zero.policy <- get("zeroPolicy", envir = .spdepOptions) 34 stopifnot(is.logical(zero.policy)) 35 36 if (family == "SMA" && method != "eigen") stop("SMA only for eigen method") 37 if (method == "spam" || method == "spam_update") stop("spam not supported as method") 38 mf <- match.call(expand.dots = FALSE) 39 m <- match(c("formula", "data", "weights", "na.action"), names(mf), 0) 40 mf <- mf[c(1, m)] 41 mf$drop.unused.levels <- TRUE 42 mf[[1]] <- as.name("model.frame") 43 mf <- eval(mf, parent.frame()) 44 mt <- attr(mf, "terms") 45 46# mt <- terms(formula, data = data) 47# mf <- lm(formula, data, , weights, na.action=na.action, 48# method="model.frame") 49 na.act <- attr(mf, "na.action") 50 if (!is.null(na.act)) { 51 subset <- !(1:length(listw$neighbours) %in% na.act) 52 listw <- subset(listw, subset, zero.policy=zero.policy) 53 } 54 55 Y <- model.extract(mf, "response") 56 if (any(is.na(Y))) stop("NAs in dependent variable") 57 X <- model.matrix(mt, mf) 58 if (any(is.na(X))) stop("NAs in independent variable") 59 n <- nrow(X) 60 if (n != length(listw$neighbours)) 61 stop("Input data and neighbourhood list have different dimensions") 62 weights <- as.vector(model.extract(mf, "weights")) 63# set up default weights 64 if (!is.null(weights) && !is.numeric(weights)) 65 stop("'weights' must be a numeric vector") 66 if (is.null(weights)) weights <- rep(as.numeric(1), n) 67 if (any(is.na(weights))) stop("NAs in weights") 68 if (any(weights < 0)) stop("negative weights") 69 lm.base <- lm(Y ~ X - 1, weights=weights) 70 aliased <- is.na(coefficients(lm.base)) 71 cn <- names(aliased) 72 names(aliased) <- substr(cn, 2, nchar(cn)) 73 if (any(aliased)) { 74 nacoef <- which(aliased) 75# bug x for X Bjarke Christensen 090924 76 X <- X[,-nacoef] 77 } 78 can.sim <- FALSE 79 if (listw$style %in% c("W", "S")) 80 can.sim <- can.be.simmed(listw) 81 82 sum_lw <- sum(log(weights)) 83# env <- new.env(parent=globalenv()) 84 env <- new.env() 85 assign("Y", Y, envir=env) 86 assign("X", X, envir=env) 87 assign("n", n, envir=env) 88 assign("weights", weights, envir=env) 89 assign("can.sim", can.sim, envir=env) 90 assign("family", family, envir=env) 91 assign("method", method, envir=env) 92 assign("verbose", verbose, envir=env) 93 assign("listw", listw, envir=env) 94 assign("sum_lw", sum_lw, envir=env) 95 W <- as(listw, "CsparseMatrix") 96 if (family == "CAR") if (!isTRUE(all.equal(W, t(W)))) 97 warning("Non-symmetric spatial weights in CAR model") 98 assign("W", W, envir=env) 99 I <- as_dsCMatrix_I(n) 100 assign("I", I, envir=env) 101 Sweights <- as(as(Diagonal(x=weights), "symmetricMatrix"), 102 "CsparseMatrix") 103 assign("Sweights", Sweights, envir=env) 104 timings[["set_up"]] <- proc.time() - .ptime_start 105 .ptime_start <- proc.time() 106 107 if (verbose) cat(paste("\nJacobian calculated using ")) 108 109 interval <- jacobianSetup(method, env, con, pre_eig=con$pre_eig, trs=trs, 110 interval=interval) 111 assign("interval", interval, envir=env) 112 113# fix SMA bounds 114 if (family == "SMA") interval <- -rev(interval) 115 116 nm <- paste(method, "set_up", sep="_") 117 timings[[nm]] <- proc.time() - .ptime_start 118 .ptime_start <- proc.time() 119 120 if (!is.null(llprof)) { 121 if (length(llprof) == 1L) 122 llprof <- seq(interval[1], interval[2], length.out=llprof) 123 ll_prof <- numeric(length(llprof)) 124 for (i in seq(along=llprof)) 125 ll_prof[i] <- .opt.fit(llprof[i], env=env, tol.solve=tol.solve) 126 nm <- paste(method, "profile", sep="_") 127 timings[[nm]] <- proc.time() - .ptime_start 128 .ptime_start <- proc.time() 129 } 130 131 opt <- optimize(.opt.fit, interval=interval, maximum=TRUE, 132 tol = con$tol.opt, env=env, tol.solve=tol.solve) 133 lambda <- opt$maximum 134 if (isTRUE(all.equal(lambda, interval[1])) || 135 isTRUE(all.equal(lambda, interval[2]))) 136 warning("lambda on interval bound - results should not be used") 137 names(lambda) <- "lambda" 138 LL <- opt$objective 139 nm <- paste(method, "opt", sep="_") 140 timings[[nm]] <- proc.time() - .ptime_start 141 .ptime_start <- proc.time() 142 143# get GLS coefficients 144 fit <- .SPAR.fit(lambda=lambda, env, out=TRUE, tol.solve=tol.solve) 145# create residuals and fitted values (Cressie 1993, p. 564) 146 fit$signal_trend <- drop(X %*% fit$coefficients) 147 fit$signal_stochastic <- drop(lambda * W %*% (Y - fit$signal_trend)) 148 fit$fitted.values <- fit$signal_trend + fit$signal_stochastic 149 fit$residuals <- drop(Y - fit$fitted.values) 150 151# get null LL 152 LL0 <- .opt.fit(lambda=0, env, tol.solve=tol.solve) 153# NK null 154 LLNullLlm <- logLik(lm(Y ~ 1, weights=weights)) 155 nm <- paste(method, "output", sep="_") 156 timings[[nm]] <- proc.time() - .ptime_start 157 .ptime_start <- proc.time() 158# if (method != "eigen") { 159# if (con$small >= n && con$small_asy) do_asy <- TRUE 160# else do_asy <- FALSE 161# } else do_asy <- TRUE 162 do_asy <- FALSE 163 if (is.null(con$fdHess)) { 164 con$fdHess <- !do_asy #&& method != "eigen" 165 fdHess <- NULL 166 } 167 stopifnot(is.logical(con$fdHess)) 168 lambda.se <- NULL 169 170 if (con$fdHess) { 171 coefs <- c(lambda, fit$coefficients) 172 fdHess <- getVcovmat(coefs, env, tol.solve=tol.solve, 173 optim=con$optimHess, optimM=con$optimHessMethod) 174 lambda.se <- sqrt(fdHess[1, 1]) 175 } 176 177 timings[["fdHess"]] <- proc.time() - .ptime_start 178 rm(env) 179 GC <- gc() 180 res <- list(fit=fit, lambda=lambda, LL=LL, LL0=LL0, call=match.call(), 181 parameters=(ncol(X)+2), aliased=aliased, method=method, family=family, 182 zero.policy=zero.policy, weights=weights, interval=interval, trs=trs, 183 timings=do.call("rbind", timings)[, c(1, 3)], LLNullLlm=LLNullLlm, 184 fdHess=fdHess, lambda.se=lambda.se, X=X, Y=Y) 185 if (!is.null(na.act)) 186 res$na.action <- na.act 187 if (is.null(llprof)) res$llprof <- llprof 188 else { 189 res$llprof <- list(lambda=llprof, ll=ll_prof) 190 } 191 if (zero.policy) { 192 zero.regs <- attr(listw$neighbours, 193 "region.id")[which(card(listw$neighbours) == 0)] 194 if (length(zero.regs) > 0L) 195 attr(res, "zero.regs") <- zero.regs 196 } 197 198 class(res) <- "spautolm" 199 res 200} 201#} 202 203#if (FALSE) { 204.opt.fit <- function(lambda, env, tol.solve=.Machine$double.eps) { 205# fitting function called from optimize() 206 SSE <- .SPAR.fit(lambda=lambda, env=env, out=FALSE, tol.solve=tol.solve) 207 n <- get("n", envir=env) 208 s2 <- SSE/n 209 ldet <- do_ldet(lambda, env) 210 det <- ifelse(get("family", envir=env) == "CAR", 0.5*ldet, ldet) 211 ret <- (det + (1/2)*get("sum_lw", envir=env) - ((n/2)*log(2*pi)) - 212 (n/2)*log(s2) - (1/(2*(s2)))*SSE) 213 if (get("verbose", envir=env)) cat("lambda:", lambda, "function:", ret, "Jacobian", ldet, "SSE", SSE, "\n") 214 ret 215} 216 217 218.SPAR.fit <- function(lambda, env, out=FALSE, tol.solve=.Machine$double.eps) { 219 dmmf <- eval(parse(text=get("family", envir=env))) 220 if (get("family", envir=env) == "SMA") IlW <- dmmf((get("I", envir=env) + 221 lambda * get("W", envir=env)), get("Sweights", envir=env)) 222 else IlW <- dmmf((get("I", envir=env) - lambda * get("W", envir=env)), 223 get("Sweights", envir=env)) 224 X <- get("X", envir=env) 225 Y <- get("Y", envir=env) 226 imat <- base::solve(crossprod(X, as.matrix(IlW %*% X)), tol=tol.solve) 227 coef <- crossprod(imat, crossprod(X, as.matrix(IlW %*% Y))) 228 fitted <- X %*% coef 229 residuals <- Y - fitted 230 SSE <- c(crossprod(residuals, as.matrix(IlW %*% residuals))) 231 if (!out) return(SSE) 232 233 n <- get("n", envir=env) 234 s2 <- SSE/n 235# var <- s2 * diag(imat) 236 coef <- c(coef) 237 names(coef) <- colnames(X) 238 res <- list(coefficients=coef, SSE=c(SSE), s2=c(s2), imat=imat, 239 N=length(residuals)) 240 res 241} 242 243# Simultaneous autoregressive 244SAR <- function(IlW, weights) { 245 t(IlW) %*% weights %*% IlW 246} 247 248# Conditional autoregressive 249CAR <- function(IlW, weights) { 250 IlW %*% weights 251} 252 253# Spatial moving average 254SMA <- function(IlW, weights) { 255 IlW <- solve(IlW) 256 t(IlW) %*% weights %*% IlW 257} 258#} 259 260print.spautolm <- function(x, ...) { 261# warning("Method print moved to the spatialreg package") 262# if (!requireNamespace("spatialreg", quietly=TRUE)) 263# stop("install the spatialreg package") 264# if (requireNamespace("spatialreg", quietly=TRUE)) { 265# return(print(x=x, ...)) 266# } 267 warning("install the spatialreg package") 268# if (FALSE) { 269 if (isTRUE(all.equal(x$lambda, x$interval[1])) || 270 isTRUE(all.equal(x$lambda, x$interval[2]))) 271 warning("lambda on interval bound - results should not be used") 272 cat("\nCall:\n") 273 print(x$call) 274 cat("\nCoefficients:\n") 275 print(coef(x)) 276 cat("\nLog likelihood:", logLik(x), "\n") 277 invisible(x) 278 279} 280#} 281 282residuals.spautolm <- function(object, ...) { 283# warning("Method residuals moved to the spatialreg package") 284# if (!requireNamespace("spatialreg", quietly=TRUE)) 285# stop("install the spatialreg package") 286# if (requireNamespace("spatialreg", quietly=TRUE)) { 287# return(residuals(object=object, ...)) 288# } 289 warning("install the spatialreg package") 290# if (FALSE) { 291 if (is.null(object$na.action)) 292 object$fit$residuals 293 else napredict(object$na.action, object$fit$residuals) 294} 295#} 296 297fitted.spautolm <- function(object, ...) { 298# warning("Method fitted moved to the spatialreg package") 299# if (!requireNamespace("spatialreg", quietly=TRUE)) 300# stop("install the spatialreg package") 301# if (requireNamespace("spatialreg", quietly=TRUE)) { 302# return(fitted(object=object, ...)) 303# } 304 warning("install the spatialreg package") 305# if (FALSE) { 306 if (is.null(object$na.action)) 307 object$fit$fitted.values 308 else napredict(object$na.action, object$fit$fitted.values) 309} 310#} 311 312deviance.spautolm <- function(object, ...) { 313# warning("Method deviance moved to the spatialreg package") 314# if (!requireNamespace("spatialreg", quietly=TRUE)) 315# stop("install the spatialreg package") 316# if (requireNamespace("spatialreg", quietly=TRUE)) { 317# return(deviance(object=object, ...)) 318# } 319 warning("install the spatialreg package") 320# if (FALSE) { 321 object$SSE 322} 323#} 324 325coef.spautolm <- function(object, ...) { 326# warning("Method coef moved to the spatialreg package") 327# if (!requireNamespace("spatialreg", quietly=TRUE)) 328# stop("install the spatialreg package") 329# if (requireNamespace("spatialreg", quietly=TRUE)) { 330# return(coef(object=object, ...)) 331# } 332 warning("install the spatialreg package") 333# if (FALSE) { 334 c(object$fit$coefficients, object$lambda) 335} 336#} 337 338 339logLik.spautolm <- function(object, ...) { 340# warning("Method logLik moved to the spatialreg package") 341# if (!requireNamespace("spatialreg", quietly=TRUE)) 342# stop("install the spatialreg package") 343# if (requireNamespace("spatialreg", quietly=TRUE)) { 344# return(logLik(object=object, ...)) 345# } 346 warning("install the spatialreg package") 347# if (FALSE) { 348 LL <- c(object$LL) 349 class(LL) <- "logLik" 350 N <- object$fit$N 351 attr(LL, "nall") <- N 352 attr(LL, "nobs") <- N 353 attr(LL, "df") <- object$parameters 354 LL 355} 356#} 357 358LR1.spautolm <- function(object) { 359# .Deprecated("spatialreg::LR1.Spautolm", msg="Method LR1.spautolm moved to the spatialreg package") 360# if (!requireNamespace("spatialreg", quietly=TRUE)) 361# stop("install the spatialreg package") 362# if (requireNamespace("spatialreg", quietly=TRUE)) { 363# return(spatialreg::LR1.Spautolm(object=object)) 364# } 365 warning("install the spatialreg package") 366# if (FALSE) { 367 if (!inherits(object, "spautolm")) stop("Not a spautolm object") 368 LLx <- logLik(object) 369 LLy <- object$LL0 370 statistic <- 2*(LLx - LLy) 371 attr(statistic, "names") <- "Likelihood ratio" 372 parameter <- 1 373 attr(parameter, "names") <- "df" 374 p.value <- 1 - pchisq(abs(statistic), parameter) 375 estimate <- c(LLx, LLy) 376 attr(estimate, "names") <- c(paste("Log likelihood of spatial regression fit"), paste("Log likelihood of OLS fit", 377 deparse(substitute(y)))) 378 method <- "Likelihood Ratio diagnostics for spatial dependence" 379 res <- list(statistic=statistic, parameter=parameter, 380 p.value=p.value, estimate=estimate, method=method) 381 class(res) <- "htest" 382 res 383} 384#} 385 386summary.spautolm <- function(object, correlation = FALSE, adj.se=FALSE, 387 Nagelkerke=FALSE, ...) { 388# warning("Method summary moved to the spatialreg package") 389# if (!requireNamespace("spatialreg", quietly=TRUE)) 390# stop("install the spatialreg package") 391# if (requireNamespace("spatialreg", quietly=TRUE)) { 392# return(summary(object=object, correlation = correlation, adj.se=adj.se, 393# Nagelkerke=Nagelkerke, ...)) 394# } 395 warning("install the spatialreg package") 396# if (FALSE) { 397 N <- object$fit$N 398 adj <- ifelse (adj.se, N/(N-length(object$fit$coefficients)), 1) 399 object$fit$s2 <- object$fit$s2*adj 400 object$resvar <- object$fit$s2*object$fit$imat 401 rownames(object$resvar) <- colnames(object$resvar) <- 402 names(object$fit$coefficients) 403 object$adj.se <- adj.se 404 405 object$rest.se <- sqrt(diag(object$resvar)) 406 object$Coef <- cbind(object$fit$coefficients, object$rest.se, 407 object$fit$coefficients/object$rest.se, 408 2*(1-pnorm(abs(object$fit$coefficients/object$rest.se)))) 409 colnames(object$Coef) <- c("Estimate", "Std. Error", 410 ifelse(adj.se, "t value", "z value"), "Pr(>|z|)") 411 if (Nagelkerke) { 412 nk <- NK.sarlm(object) 413 if (!is.null(nk)) object$NK <- nk 414 } 415 if (correlation) { 416 object$correlation <- diag((diag(object$resvar)) 417 ^(-1/2)) %*% object$resvar %*% 418 diag((diag(object$resvar))^(-1/2)) 419 dimnames(object$correlation) <- dimnames(object$resvar) 420 } 421 object$LR1 <- LR1.spautolm(object) 422 rownames(object$Coef) <- names(object$fit$coefficients) 423 structure(object, class=c("summary.spautolm", class(object))) 424} 425#} 426 427print.summary.spautolm <- function(x, digits = max(5, .Options$digits - 3), 428 signif.stars = FALSE, ...) 429{ 430# warning("Method print moved to the spatialreg package") 431# if (!requireNamespace("spatialreg", quietly=TRUE)) 432# stop("install the spatialreg package") 433# if (requireNamespace("spatialreg", quietly=TRUE)) { 434# return(print(x=x, digits = digits, 435# signif.stars = signif.stars, ...)) 436# } 437 warning("install the spatialreg package") 438# if (FALSE) { 439 cat("\nCall: ", deparse(x$call), sep = "", fill=TRUE) 440 if (isTRUE(all.equal(x$lambda, x$interval[1])) || 441 isTRUE(all.equal(x$lambda, x$interval[2]))) 442 warning("lambda on interval bound - results should not be used") 443 cat("\nResiduals:\n") 444 resid <- residuals(x) 445 nam <- c("Min", "1Q", "Median", "3Q", "Max") 446 rq <- if (length(dim(resid)) == 2L) 447 structure(apply(t(resid), 1, quantile), dimnames = list(nam, 448 dimnames(resid)[[2]])) 449 else structure(quantile(resid), names = nam) 450 print(rq, digits = digits, ...) 451 if (x$zero.policy) { 452 zero.regs <- attr(x, "zero.regs") 453 if (!is.null(zero.regs)) 454 cat("\nRegions with no neighbours included:\n", 455 zero.regs, "\n") 456 } 457 cat("\nCoefficients:", x$coeftitle, "\n") 458 coefs <- x$Coef 459 if (!is.null(aliased <- x$aliased) && any(x$aliased)){ 460 cat(" (", table(aliased)["TRUE"], 461 " not defined because of singularities)\n", sep = "") 462 cn <- names(aliased) 463 coefs <- matrix(NA, length(aliased), 4, dimnames = list(cn, 464 colnames(x$Coef))) 465 coefs[!aliased, ] <- x$Coef 466 } 467 printCoefmat(coefs, signif.stars=signif.stars, digits=digits, 468 na.print="NA") 469 res <- x$LR1 470 cat("\nLambda:", format(signif(x$lambda, digits)), 471 "LR test value:", format(signif(res$statistic, digits)), 472 "p-value:", format.pval(res$p.value, digits), 473 "\n") 474 if (!is.null(x$lambda.se)) 475 cat("Numerical Hessian standard error of lambda:", 476 format(signif(x$lambda.se, digits)), "\n") 477 cat("\nLog likelihood:", logLik(x), "\n") 478 if (x$adj.se) cat("Residual variance (sigma squared): ") 479 else cat("ML residual variance (sigma squared): ") 480 cat(format(signif(x$fit$s2, digits)), ", (sigma: ", 481 format(signif(sqrt(x$fit$s2), digits)), ")\n", sep="") 482 cat("Number of observations:", x$fit$N, "\n") 483 cat("Number of parameters estimated:", x$parameters, "\n") 484 cat("AIC: ", format(signif(AIC(x), digits)), "\n", sep="") 485 if (!is.null(x$NK)) cat("Nagelkerke pseudo-R-squared:", 486 format(signif(x$NK, digits)), "\n") 487 correl <- x$correlation 488 if (!is.null(correl)) { 489 p <- NCOL(correl) 490 if (p > 1) { 491 cat("\nCorrelation of Coefficients:\n") 492 correl <- format(round(correl, 2), nsmall = 2, 493 digits = digits) 494 correl[!lower.tri(correl)] <- "" 495 print(correl[-1, -p, drop = FALSE], quote = FALSE) 496 } 497 } 498 cat("\n") 499 invisible(x) 500} 501#} 502 503#if (FALSE) { 504getVcovmat <- function(coefs, env, tol.solve=.Machine$double.eps, optim=FALSE, 505 optimM="optimHess") { 506 if (optim) { 507 if (optimM == "nlm") { 508 options(warn=-1) 509 opt <- nlm(f=f_spautolm_hess_nlm, p=coefs, env=env, hessian=TRUE) 510 options(warn=0) 511 mat <- opt$hessian 512# opt <- optimHess(par=coefs, fn=f_laglm_hess, env=env) 513# mat <- opt 514 } else if (optimM == "optimHess") { 515 mat <- optimHess(par=coefs, fn=f_spautolm_hess, env=env) 516 } else { 517 opt <- optim(par=coefs, fn=f_spautolm_hess, env=env, method=optimM, 518 hessian=TRUE) 519 mat <- opt$hessian 520 } 521# opt <- optimHess(par=coefs, fn=f_spautolm_hess, env=env) 522# mat <- opt 523 } else { 524 fd <- fdHess(coefs, f_spautolm_hess, env) 525 mat <- fd$Hessian 526 } 527 res <- solve(-(mat), tol.solve=tol.solve) 528 res 529} 530 531f_spautolm_hess_nlm <- function(coefs, env) { 532 ret <- f_spautolm_hess(coefs, env) 533 -ret 534} 535 536f_spautolm_hess <- function(coefs, env) { 537 lambda <- coefs[1] 538 int <- get("interval", envir=env) 539 if (lambda <= int[1] || lambda >= int[2]) return(-Inf) 540 beta <- coefs[-1] 541 X <- get("X", envir=env) 542 Y <- get("Y", envir=env) 543 fitted <- X %*% beta 544 residuals <- Y - fitted 545 dmmf <- eval(parse(text=get("family", envir=env))) 546 if (get("family", envir=env) == "SMA") IlW <- dmmf((get("I", envir=env) + 547 lambda * get("W", envir=env)), get("Sweights", envir=env)) 548 else IlW <- dmmf((get("I", envir=env) - lambda * get("W", envir=env)), 549 get("Sweights", envir=env)) 550 SSE <- c(crossprod(residuals, as.matrix(IlW %*% residuals))) 551 n <- get("n", envir=env) 552 s2 <- SSE/n 553 ldet <- do_ldet(lambda, env) 554 det <- ifelse(get("family", envir=env) == "CAR", 0.5*ldet, ldet) 555 ret <- (det + (1/2)*get("sum_lw", envir=env) - ((n/2)*log(2*pi)) - 556 (n/2)*log(s2) - (1/(2*(s2)))*SSE) 557 if (get("verbose", envir=env)) 558 cat("lambda:", lambda, "function:", ret, "Jacobian", ldet, "SSE", 559 SSE, "\n") 560 if (!is.finite(ret)) return(-Inf) 561 ret 562} 563 564#} 565 566 567 568