| # File src/library/stats/R/logLik.R |
| # Part of the R package, https://www.R-project.org |
| # |
| # Copyright (C) 2001-2014 The R Core Team |
| # |
| # This program is free software; you can redistribute it and/or modify |
| # it under the terms of the GNU General Public License as published by |
| # the Free Software Foundation; either version 2 of the License, or |
| # (at your option) any later version. |
| # |
| # This program is distributed in the hope that it will be useful, |
| # but WITHOUT ANY WARRANTY; without even the implied warranty of |
| # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| # GNU General Public License for more details. |
| # |
| # A copy of the GNU General Public License is available at |
| # https://www.R-project.org/Licenses/ |
| |
| ### originally from package nlne. |
| |
| |
| logLik <- function(object, ...) UseMethod("logLik") |
| |
| print.logLik <- function(x, digits = getOption("digits"), ...) |
| { |
| cat("'log Lik.' ", paste(format(c(x), digits = digits), collapse = ", "), |
| " (df=", format(attr(x,"df")), ")\n", sep = "") |
| invisible(x) |
| } |
| |
| str.logLik <- function(object, digits = max(2L, getOption("digits") - 3L), |
| vec.len = getOption("str")$vec.len, ...) |
| { |
| cl <- oldClass(object) |
| len <- length(co <- c(object)) |
| cutl <- len > vec.len |
| cat("Class", if (length(cl) > 1L) "es", |
| " '", paste(cl, collapse = "', '"), "' : ", |
| paste |
| (format(co[seq_len(min(len,vec.len))], digits = digits), |
| collapse = ", "), if(cutl) ", ...", |
| " (df=", format(attr(object,"df")), ")\n", sep = "") |
| } |
| |
| ## rather silly (but potentially used in pkg nlme): |
| as.data.frame.logLik <- function (x, ...) |
| as.data.frame(c(x), ...) |
| |
| ## >> logLik.nls() in nls.R |
| |
| ## from package:nlme |
| |
| ## log-likelihood for glm objects |
| logLik.glm <- function(object, ...) |
| { |
| if(!missing(...)) warning("extra arguments discarded") |
| fam <- family(object)$family |
| p <- object$rank |
| ## allow for estimated dispersion |
| if(fam %in% c("gaussian", "Gamma", "inverse.gaussian")) p <- p + 1 |
| val <- p - object$aic / 2 |
| ## Note: zero prior weights have NA working residuals. |
| attr(val, "nobs") <- sum(!is.na(object$residuals)) |
| attr(val, "df") <- p |
| class(val) <- "logLik" |
| val |
| } |
| |
| ## log-likelihood for lm objects |
| logLik.lm <- function(object, REML = FALSE, ...) |
| { |
| if(inherits(object, "mlm")) |
| stop("'logLik.lm' does not support multiple responses") |
| res <- object$residuals # not resid(object) because of NA methods |
| p <- object$rank |
| N <- length(res) |
| if(is.null(w <- object$weights)) { |
| w <- rep.int(1, N) |
| } else { |
| ## this is OK as both resids and weights are for the cases used |
| excl <- w == 0 # eliminating zero weights |
| if (any(excl)) { |
| res <- res[!excl] |
| N <- length(res) |
| w <- w[!excl] |
| } |
| } |
| N0 <- N |
| if(REML) N <- N - p |
| val <- .5* (sum(log(w)) - N * (log(2 * pi) + 1 - log(N) + |
| log(sum(w*res^2)))) |
| if(REML) val <- val - sum(log(abs(diag(object$qr$qr)[1L:p]))) |
| attr(val, "nall") <- N0 # NB, still omits zero weights |
| attr(val, "nobs") <- N |
| attr(val, "df") <- p + 1 |
| class(val) <- "logLik" |
| val |
| } |
| |
| ## Since AIC, BIC call logLik on everything |
| logLik.logLik <- function(object, ...) object |
| |
| nobs <- function(object, ...) UseMethod("nobs") |
| |
| ## also used for mlm fits |
| nobs.lm <- function(object, ...) |
| if(!is.null(w <- object$weights)) sum(w != 0) else NROW(object$residuals) |
| |
| nobs.glm <- function(object, ...) |
| if(!is.null(w <- object$prior.weights)) sum(w != 0) else length(object$residuals) |
| |
| nobs.logLik <- function(object, ...) { |
| res <- attr(object, "nobs") |
| if (is.null(res)) stop("no \"nobs\" attribute is available") |
| res |
| } |
| |
| nobs.nls <- function(object, ...) |
| if (is.null(w <- object$weights)) length(object$m$resid()) else sum(w != 0) |
| |
| ## it is probably too unsafe to use residuals generally, not least |
| ## because of e.g. weighted fits. |
| nobs.default <- function(object, use.fallback = FALSE, ...) |
| { |
| ## MASS::loglm and MASS::polr fits have an 'nobs' component |
| if((is.L <- is.list(object)) && !is.null(n <- object[["nobs"]])) n |
| ## cov.wt() unfortunately uses 'n.obs': |
| else if(is.L && !is.null(n <- object[["n.obs"]])) n |
| else if(use.fallback) { |
| if(!is.null(w <- object[["weights"]])) sum(w != 0) |
| else if("residuals" %in% names(object)) |
| NROW(object$residuals) # and not residuals(object) |
| ## perhaps sum(!is.na(object$residuals)) ? |
| else { |
| warning("no 'nobs' method is available") |
| 0L # which is what object$residuals used to give. |
| } |
| } else stop("no 'nobs' method is available") # or maybe NA_integer_ |
| } |