Nothing
RsqGLM <- function(model = NULL, obs = NULL, pred = NULL, use = "pairwise.complete.obs",
plot = TRUE, plot.type = "lollipop", ...) {
# version 1.9 (5 Jan 2023)
model.provided <- ifelse(is.null(model), FALSE, TRUE)
if (model.provided) {
if (!("glm" %in% class(model))) stop ("'model' must be of class 'glm'.")
if (!is.null(pred)) message("Argument 'pred' ignored in favour of 'model'.")
if (!is.null(obs)) message("Argument 'obs' ignored in favour of 'model'.")
obs <- model$y
pred <- model$fitted.values
} else { # if model not provided
if (is.null(obs) | is.null(pred)) stop ("You must provide either 'obs' and 'pred', or a 'model' object of class 'glm'.")
if (length(obs) != length(pred)) stop ("'obs' and 'pred' must have the same number of values (and in the same order).")
pred[pred == 0] <- 2e-16 # avoid NaN in log below
pred[pred == 1] <- 1 - 2e-16 # avoid NaN in log below
# new (15 Sep 2015):
dat <- data.frame(obs, pred)
n.in <- nrow(dat)
dat <- na.omit(dat)
n.out <- nrow(dat)
if (n.out < n.in) warning (n.in - n.out, " observations removed due to missing data; ", n.out, " observations actually evaluated.")
obs <- dat$obs
pred <- dat$pred
if (any(!(obs %in% c(0, 1)) | pred < 0 | pred > 1)) stop ("Sorry, 'obs' and 'pred' options only implemented for binomial GLMs (binary observations with values 0 or 1, and probabilistic predictions with values between 0 and 1). For GLMs of other types, you need to provide a 'model' object.")
logit <- log(pred / (1 - pred))
model <- glm(obs ~ logit, family = "binomial")
}
null.mod <- glm(obs ~ 1, family = family(model))
loglike.M <- as.numeric(logLik(model))
loglike.0 <- as.numeric(logLik(null.mod))
#like.M <- exp(loglike.M)
#like.0 <- exp(loglike.0)
#like.M <- like.M + 2.2e-16
#like.0 <- like.0 + 2.2e-16 # zeros were causing problems, but this too
N <- length(obs)
# from http://stackoverflow.com/questions/3337044/model-fit-statistics-for-a-logistic-regression/15949904#15949904:
CoxSnell <- 1 - exp((model$deviance - model$null.deviance) / 2 * N)
Nagelkerke <- CoxSnell / (1 - exp((- model$null.deviance) / N))
McFadden <- 1 - ((model$deviance / -2)/(model$null.deviance / -2))
# based on Nagelkerke 1991:
CoxSnell <- 1 - exp(-(2 / N) * (loglike.M - loglike.0))
Nagelkerke <- CoxSnell / (1 - exp((2 * N ^ (-1)) * loglike.0))
# based on Allison 2014:
#CoxSnell <- 1 - (loglike.0 / loglike.M) ^ (2 / N)
#Nagelkerke <- CoxSnell / (1 - loglike.0 ^ (2 / N))
McFadden <- 1 - (loglike.M / loglike.0)
if (family(model)$family == "binomial" & family(model)$link == "logit")
#Tjur <- abs(diff(unique(ave(pred, obs))))
#Tjur <- diff(tapply(pred, obs, mean))
Tjur <- mean(pred[obs == 1]) - mean(pred[obs == 0])
else {
Tjur <- NA
message("NOTE: Tjur R-squared applies only to binomial GLMs")
}
sqPearson <- cor(obs, pred, use = use) ^ 2
if (plot) {
plot.type <- match.arg(plot.type, c("lollipop", "barplot"))
m <- data.frame(CoxSnell = CoxSnell, Nagelkerke = Nagelkerke, McFadden = McFadden, Tjur = Tjur, sqPearson = sqPearson)
if (plot.type == "barplot") barplot(as.matrix(m), las = 2, ...)
if (plot.type == "lollipop") lollipop(unlist(as.vector(m)), las = 2, ymin = 0, ylab = "", ...)
else message("Invalid 'plot.type'. Plot not produced.")
}
return(list(CoxSnell = CoxSnell, Nagelkerke = Nagelkerke, McFadden = McFadden, Tjur = Tjur, sqPearson = sqPearson))
}
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.